How to manually simulate linear system without lsim
21 views (last 30 days)
Show older comments
I have a 3rd order continuous system of the following form identified with System Identification Toolbox:
y = Cx + Du + e
dx/dt = Ax + Bu + Ke
I can simulate it by using lsim as follows:
simulatedOutput = lsim(systemModel, inputVector, timeVector, initialStateVector)
Now I need to convert this model into code that uses only loops and simple matrix/vector arithmetic. I cannot use lsim, ODE, Runge-Kutta, or any other Matlab-specific solver or functions (I used only interp1 which is easy to convert into code).
I wrote a simple script that will read the same input file used to identify the model, but the output is completely wrong, and I cannot figure out why.
Fig. 1: output from lsim (expected)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1680636/image.png)
Fig. 2: output from manual simulation (actual - code below)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1680641/image.png)
The following is my attempt to simulate the identified system by using a Matlab script with only basic vector/matrix arithmetic and a loop:
% A weights (3x3 matrix)
A = [ ...
0.0356, -3.4131, -14.9525;
-1.0591, 85.8128, 334.3098;
1.5729, -95.1123, -175.5517;
];
% B weights (3x1 vector)
B = [ ...
-0.0019;
0.0403;
-0.0225;
];
% C weights (1x3 vector)
C = [ -5316.9, 24.87, 105.92 ];
% D weight (scalar)
D = 0;
% K weights (3x1 vector)
K = [ ...
-0.0025;
-0.0582;
0.0984;
];
% Initial state
x0 = [ ...
-0.0458;
0.0099;
-0.0139;
];
% Read input file
csv = readmatrix('https://raw.githubusercontent.com/01binary/systemid/main/input.csv');
% Select column vectors
time = csv(:,1);
measurement = csv(:,2);
input = csv(:,3);
% Simulate
x = x0;
timeStep = 0.02;
output = zeros(length(input), 1);
for i = 1:length(input)
t = time(1) + (i - 1) * timeStep;
u = interp1(time, input, t);
[y, dx] = systemModel(A, B, C, D, K, x, u, 0);
x = x + dx * timeStep;
output(i) = y;
end
% Plot against original measurements
plot(time, measurement, time, output);
function [y, dx] = systemModel(A, B, C, D, K, x, u, e)
% y = Cx + Du + e
y = ...
C * x + ... % Add contribution of state
D * u + ... % Add contribution of input
e; % Add disturbance
% dx/dt = Ax + Bu + Ke
dx = ...
A * x + ... % Add contribution of state
B * u + ... % Add contribution of input
e * K; % Add contribution of disturbance
end
I am trying to keep everything really simple and minimal... could someone spot what I am missing?
Thank you!
0 Comments
Accepted Answer
Paul
on 27 Apr 2024
Hi Valeriy,
As indicated, the system in question is a differential equation model of a continuous-time sytem.
This line of code
% x = systemState(x, u, 0, A, B, K) * timeStep;
should be
% x = x + systemState(x, u, 0, A, B, K) * timeStep;
because systemState resturns xdot.
But that correction won't match lsim. As the doc page states: "For continuous-time systems, lsim first discretizes the system using c2d, and then propagates the resulting discrete-time state-space equations. Unless you specify otherwise with the method input argument, lsim uses the first-order-hold discretization method when the input signal is smooth, and zero-order hold when the input signal is discontinuous, such as for pulses or square waves."
So if you want to exactly match lsim, you'll have to implement the c2d transformation and then propagate the resulting discrete-time equation.
15 Comments
Paul
on 28 Apr 2024
Read in the data and plot it
T = readtable('input.csv')
figure
plot(T.time,T.PWM,T.time,T.reading),grid
legend('PWM (input)','reading (output)')
Create iddata object. We'll assume the samples are space at 0.02, though I think there is a sample or two towards the end of the data that has a corrupted sample time (or maybe a data point is missing)
data = iddata(T.reading,T.PWM,0.02);
Identify a third order, discrete-time model with the sampling period the same as the sampling period of the data.
sysd = n4sid(data,3);
sysd
To answer @Torsten's question, we see that the form of the model already has the 0.02 "built-in," i.e., the model equations don't include an explicit term of Ts.
The model is stable, even though we did not set the EnforceStability option to true (the default is false)
abs(eig(sysd.A))
Identify a continuous-time model
sysc = n4sid(data,3,'Ts',0);
This model is also stable, with two poles that have a very fast time constant.
format short e
eig(sysc.A)
If we Euler discretize this model, we see that the poles of the discretized system are well outside the unit circle, so we'd expect an unstable result, which is exactly what was seen from using lsim with inputs spaced by 0.02 sec (lsim doesn't use Euler integration, but the idea is the same).
abs(eig(eye(3) + sysc.A*0.02))
Propagation equations for the discrete-time model
x = sysd.X0;
for ii = 1:numel(T.time)
yd(ii) = sysd.C*x;
x = sysd.A*x + sysd.B*T.PWM(ii);
end
Propagation equations for the continuous-time model
% x = sysc.X0;
% for ii = 1:numel(T.time)
% yc(ii) = sysc.C*x;
% x = x + 0.02*(sysc.A*x + sysc.B*T.PWM(ii));
% end
figure
plot(T.time-T.time(1),T.reading, ...
T.time-T.time(1),lsim(sysd,T.PWM,(0:numel(T.time)-1)*0.02,sysd.X0), ...
T.time-T.time(1),yd),grid %...
legend('reading','lsim','yd')
% T.time-T.time(1),lsim(sysc,T.PWM,(0:numel(T.time)-1)*0.02,sysc.X0), ...
% T.time-T.time(1),yc),grid
More Answers (1)
Sam Chak
on 27 Apr 2024
Hi @Valeriy
Other than lsim, you have two options for solving the differential equations: you can use the built-in ODE solver or a custom numerical solver like Runge-Kutta 4. Below, you can find a comparison of results between the lsim command and the ode45 solver.
%% data extraction
csv = readmatrix('input.csv');
time = csv(:,1);
measurement = csv(:,2);
input = csv(:,3);
%% data processing (lsim requires evenly-spaced input signal)
ta = time(1);
tb = time(end);
tspan = linspace(ta, tb, (tb-ta)/0.01 + 1);
u = interp1(time, input, tspan);
%% state-space system
A = [0.0356, -3.4131, -14.9525;
-1.0591, 85.8128, 334.3098;
1.5729, -95.1123, -175.5517];
B = [-0.0019;
0.0403;
-0.0225];
C = [-5316.9, 24.87, 105.92];
D = 0;
% x'= A·x + B·u
% y = C·x + D·u
sys = ss(A, B, C, D);
%% initial values
x0 = [-0.0458;
0.0099;
-0.0139];
Approach #1: Run simulation using 'lsim()' command
%% response of state-space system via lsim
lsim(sys, u, tspan, x0), grid on
Approach #2: Run simulation using 'ode45()' solver
%% call ode45 solver
options = odeset('Reltol', 1e-9, 'Abstol', 1e-6);
[t, x] = ode45(@(t, x) odesys(t, x, A, B, C, D, time, input), tspan, x0, options);
[~, y] = odesys(t', x', A, B, C, D, time, input);
plot(t, y), grid on, xlim([0 round(tb)])
xlabel('t'), ylabel('y(t)'), title('ode45 Simulation Results')
%% Identified State-Space System
function [dx, y] = odesys(t, x, A, B, C, D, time, input)
u = interp1(time, input, t);
y = C*x + D*u(:,1);
dx = A*x + B*u(:,1);
end
3 Comments
Sam Chak
on 28 Apr 2024
Hi @Valeriy
I'm glad to hear that you have resolved the problem. Initially, I misunderstood your intention and thought you wanted to replicate the response of the continuous-time system generated by lsim. That's why I suggested using a custom-designed numerical algorithm, which doesn't necessarily have to be RK4.
Previously, I was unaware that this is a homework exercise aimed at testing a specific time integration method for solving the initial value problem of the identified state-space system. By the way, you mentioned that the Runge-Kutta method is not allowed, but it's worth noting that the standard Euler method is generally considered as a 1st-order Runge-Kutta method.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!