Matlab Coder Error: The output of the ODE function must be a column vector

I tried to generate C code from the following m-file:

function [t1,y1,t2,y2,t3,y3] = FMD(m, d, c)
    %definition of the transfer function of a spring mass damper system
    num = 1;
    denum = [m d c];
    %tranformation to a state space system
    [A,B,C,D] = tf2ss(num, denum);
    
    %definition of three input signals (sine, exponential and step)
    ts = 0:0.01:8;
    u1 = sin(10*ts);
    u2 = exp(ts);
    u3 = zeros(1,length(ts));
    u3(ts>=1 & ts<max(ts)) = 1;
    
    %time span for the numerical solver
    tspan = [2 5];
    %initial condition of the state space system
    iniCon = [0;0];
    
    %numerical solving of the ODEs for the three different input signals
    [t1, x1] = ode45(@(t,x) sys(t,x,u1), tspan, iniCon);
    [t2, x2] = ode45(@(t,x) sys(t,x,u2), tspan, iniCon);
    [t3, x3] = ode45(@(t,x) sys(t,x,u3), tspan, iniCon);
  
    
    %preallocation of output variables
    y1 = zeros(length(x1),1);
    y2 = zeros(length(x2),1);
    y3 = zeros(length(x3),1);
    
    %calculation of the output values for every time step (element-wise
    %output function of the state space system)
    for i=1:1:length(x1)
        y1(i,:) = C*transpose(x1(i,:));
    end
    
    for i=1:1:length(x2)
        y2(i,:)=C*transpose(x2(i,:));
    end
    
    for i=1:1:length(x3)
        y3(i,:)=C*transpose(x3(i,:));
    end

    %nested function for the definition of the state space system
    function  dx = sys(t, x, u)
       %interpolate the data set (ts,u3) at time t
       u3I = interp1(ts, u, t);
       dx = A*x + B*u3I;
    end

end

When I try to run the code in Matlab by calling

clear;
[t1,y1,t2,y2,t3,y3]=FMD(2,5,10);
plot(t1,y1);

everything works as expected. However when I want to generate C code out the function FMD, the first error that I get is the following one:

.

As the program works fine in Matlab without having trouble with the resulting vector, I do not get the problem and I would be very glad if someone could help me to understand and solve the problem.

Answer 1

[Mike Hosea]

I have to admit the problem here is pretty hard to diagnose. The problem actually traces back to the definitions of A and B:

[A,B,C,D] = tf2ss(num, denum);

After this call, codegen thinks A is :2-by-:2 and B is :2-by-:1. This is because we have

>> [A,B,C,D] = tf2ss(1,[1,1,1]); size(A), size(B)
ans =
     2     2
ans =
     2     1
>> [A,B,C,D] = tf2ss(1,[0,1,1]); size(A), size(B)
ans =
     1     1
ans =
     1     1
>> [A,B,C,D] = tf2ss(1,[0,0,1]); size(A), size(B)
ans =
     0     0
ans =
     0     0

But the assumption built into the rest of the program is clearly just the first case, so I assume that we know for certain that m is nonzero, A is 2-by-2, and B is 2-by-1. Codegen isn't so sure. All it knows is that A is :2-by-:2 and B is :2-by-:1. It also knows that inside sys x is 2-by-1 and t is scalar, hence u3I is scalar. I guess the problem theoretically occurs if A were 1-by-2 and B were 2-by-0, then A*x is scalar, and B*u3I is 2-by-0. Consequently, A*x + B*u3I is 2-by-0, which is not a column vector. That can't happen, but codegen doesn't know this.

There are multiple ways of fixing this. One way is to add the coder.noImplicitExpansionInFunction declaration to sys:

    function  dx = sys(t, x, u)
    coder.noImplicitExpansionInFunction;
    %interpolate the data set (ts,u3) at time t
    u3I = interp1(ts, u, t);
    dx = A*x + B*u3I;
    end

I think that works because by turning off implicit expansion, you also eliminate support for run-time scalar expansion, so the codegen will require and enforce that A*x and B*u3I have the same number of rows and columns (because they are added).

If that seems too mysterious, you can fix it where B is used:

    function  dx = sys(t, x, u)
       %interpolate the data set (ts,u3) at time t
       u3I = interp1(ts, u, t);
       dx = A*x + B(:,1)*u3I;
    end

My recommendation, however is to fix it where B is defined:

    [Atmp,Btmp,C,D] = tf2ss(num, denum);
    A = Atmp(1:2,1:2);
    B = Btmp(1:2,1);

because now A and B are fixed in size. That should generate better code.

[Mike Hosea]

BTW, if you don't mind having the same t values for each signal, you might find some extra performance from the generated code by combining the integrations as follows. This generalizes to any number of input signals.

function [t,y] = FMD(m, d, c)
%definition of the transfer function of a spring mass damper system
num = 1;
denum = [m d c];

%tranformation to a state space system
[Atmp,Btmp,C,D] = tf2ss(num, denum);
A = Atmp(1:2,1:2);
B = Btmp(1:2,1);

%definition of three input signals (sine, exponential and step)
ts = 0:0.01:8;
u1 = sin(10*ts);
u2 = exp(ts);
u3 = zeros(1,length(ts));
u3(ts>=1 & ts<max(ts)) = 1;

U = [u1(:),u2(:),u3(:)]; % U has one column for each signal.

%time span for the numerical solver
tspan = [2 5];
%initial condition of the state space system
iniCon = [0;0];
% Stack the initial conditions.
iniConStack = repmat(iniCon,size(U,2),1);

[t, x] = ode45(@(t,x) sys(t,x,U), tspan, iniConStack);
% Now x has numel(t) rows and 2*size(U,2) columns.
% Separate the adjacent pairs of columns into pages of 3D array, so that
% x(:,1:2) --> x(:,:,1),
% x(:,3:4) --> x(:,:,2), and so forth.
x = reshape(x,numel(t),2,size(U,2));
% Now we have the columns corresponding to the k-th signal as x(:,:,k).
% Allocate the return value.
y = zeros(numel(t),size(U,2));
% Populate the return value.
for k = 1:size(U,2)
    y(:,k) = x(:,:,k)*C.';
end

    function  dx = sys(t, x, U)
        ui = interp1(ts, U, t); % ui is 1-by-size(U,2).
        xr = reshape(x,2,size(U,2)); % Unstack x.
        dx = A*xr + B.*ui; % Compute the dx values for each column of xr.
        dx = dx(:); % Stack the results back up as column vector.
    end

end

[SKRE]

Thank you very much for the quick and detailed answer, including the tip on how to increase performance.