Clear Filters
Clear Filters

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

3 views (last 30 days)
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.

Accepted Answer

Mike Hosea
Mike Hosea on 17 Oct 2022
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.
  2 Comments
Mike Hosea
Mike Hosea on 17 Oct 2022
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
SKRE on 18 Oct 2022
Thank you very much for the quick and detailed answer, including the tip on how to increase performance.

Sign in to comment.

More Answers (0)

Categories

Find more on MATLAB Coder in Help Center and File Exchange

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!