How to manually simulate linear system without lsim

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)

Fig. 2: output from manual simulation (actual - code below)

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!

Answer 1

[Paul]

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.

[Valeriy]

I still can't get it working. Perhaps it's because I don't know what "discretization" means - should I identify the system again using a Discrete State Space model, so the discretization is built-in?

The output doesn't need to match lsim exactly - just look "similar" to the original measurement data in the CSV file. There is a column called "measurement", and I am expecting this output to somewhat resemble the data in that column, except for the noise:

I am getting the following from the code below. Notably, my output should be between 100 and 400, roughly, but you can see below it's between -1 x 10^163 and 3 x 10^163.

Here's my updated code - I simplified it with the assumption that the time step is constant (0.02).

% 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 (3x1 vector)
x0 = [ ...
  -0.0458;
  0.0099;
  -0.0139;
];

% Read input file
csv = readmatrix('https://raw.githubusercontent.com/01binary/systemid/main/input.csv');
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)
  u = input(i);
  [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

  % x = Ax + Bu + Ke
  dx = ...
    A * x + ... % Add contribution of state
    B * u + ... % Add contribution of input
    e * K;      % Add contribution of disturbance
end

[Valeriy]

OK, I re-identified the same system using Discrete State-Space Model and got different A, B, C, D, K constants which are much smaller. Am I correct to assume that discretization is baked into these matrices now, so I can use the equations directly as long as my time step remains 0.02, the time step the system was identified with?

I still can't get this working, don't understand why. Tried both with and without multiplying by time step.

% Same code but ABCDK constants are for a discrete state-space system

% A weights (3x3 matrix)
A = [ ...
  0.9988,   0.05193,  -0.02261;
  0.02222,   -0.01976,  0.7353;
  0.0009856,  -0.2093, -0.5957;
];

% B weights (3x1 vector)
B = [ ...
  -2.663e-06;
  5.727e-05;
  -0.0001872;
];

% C weights (1x3 vector)
C = [ -5317, 24.87, 105.9 ];

% D weight (scalar)
D = 0;

% K weights (3x1 vector)
K = [ ...
  -0.0001655;
  -0.001508;
  6.209e-06;
];

% Initial state (3x1 vector)
x0 = [ ...
  -0.0458;
  0.0099;
  -0.0139;
];

% Read input file
csv = readmatrix('https://raw.githubusercontent.com/01binary/systemid/main/input.csv');
time = csv(:,1);
measurement = csv(:,2);
input = csv(:,3);

% Simulate
x = x0;
output = zeros(length(input), 1);

for i = 1:length(input)
  u = input(i);
  [y, dx] = systemModel(A, B, C, D, K, x, u, 0);
  x = x + dx;
  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 = Ax + Bu + Ke
  dx = ...
    A * x + ... % Add contribution of state
    B * u + ... % Add contribution of input
    e * K;      % Add contribution of disturbance
end

[Paul]

Please post the complete code, including reading the data from the .csv file and generating the identified model.

[Torsten]

Use

% Simulate
x = x0;
output = zeros(length(input), 1);
time = [0;time];
for i = 1:length(input)
  u = input(i);
  [y, dx] = systemModel(A, B, C, D, K, x, u, 0);
  x = x + dx*(time(i+1)-time(i));
  output(i) = y;
end
time = time(2:end);
% Plot against original measurements
plot(time, measurement, time, output);

instead of

% Simulate
x = x0;
output = zeros(length(input), 1);
for i = 1:length(input)
  u = input(i);
  [y, dx] = systemModel(A, B, C, D, K, x, u, 0);
  x = x + dx;
  output(i) = y;
end
% Plot against original measurements
plot(time, measurement, time, output);

[Valeriy]

Thanks for the reply! I tried this, but I am still getting output that looks nothing like the system:

Actual result: Min value is around zero and max is around 12 x 10^5, with some kind of exponential growth

Expected result: something that has 2-3 stair-step peaks, between 100 and 400.

[Torsten]

The curve you get for y looks smooth to me. So I guess there must be something wrong with your input matrices/vectors - the numerical computation with Euler method seems stable.

Note that your matrix A has 1 as eigenvalue. This will make the solution explode as t -> Inf.

[Valeriy]

@Paul - thanks for taking your time to look at this question.

The files can be found here:

https://github.com/01binary/systemid

The input file:

https://github.com/01binary/systemid/blob/main/input.csv

The script file (updated with my system re-identified as "discrete" instead of continuous):

https://github.com/01binary/systemid/blob/main/simulate.m

I recorded a video of how I identified this system by using System Identification Toolbox, and thus filled in the ABCDK constants in the script file:

https://youtu.be/A47aR6UAewo

[Torsten]

Could you include the lsim curve for y derived from your system identification ?

[Paul]

Just copy/paste the code into a comment like you did previously, but include the steps that read the data from the file and generate the identified system.

[Valeriy]

@Torsten I attempted to simulate my (re-identified) now-discrete system with lsim, but it kept giving me error messages about time step not matching. I ended up re-doing all the code and found that one of the matrices still had values for the continuous system while the rest had new values for the discrete system.

After I fixed that, I also noticed that the equations for the discrete system do not accumulate the input. The input is directly the value from the previous step.

With those two adjustments, it started working!

I will accept Paul's answer because it sent me in the direction of re-identifying the system as discrete.

[Valeriy]

Final working code for the accepted solution:

• According to accepted solution by Paul, I had to re-identify the system as a discrete-time state-space model. Originally I had a continuous-time state-space model.
• Re-identifying the system as discrete with time step 0.02 has made it possible to use simple equations to update system state and calculate predictions (see code below).

% A weights (3x3 matrix)
A = [ ...
  0.9988,   0.05193,  -0.02261;
  0.02222,   -0.01976,  0.7353;
  0.0009856,  -0.2093, -0.5957;
];

% B weights (3x1 vector)
B = [ ...
  -0.00000266;
  0.0000572747;
  -0.0001872152;
];

% C weights (1x3 vector)
C = [ ...
  -5316.903919, ...
  24.867656, ...
  105.92416 ...
];

% D weight (scalar)
D = 0;

% K weights (3x1 vector)
K = [ ...
  -0.0001655;
  -0.001508;
  6.209e-06;
];

% Initial state (3x1 vector)
x0 = [ ...
  -0.0458;
  0.0099;
  -0.0139;
];

% Read input
csv = readmatrix('https://raw.githubusercontent.com/01binary/systemid/main/input.csv');
time = csv(:,1);
measurement = csv(:,2);
input = csv(:,3);

% Simulate
x = x0;
output = zeros(length(input), 1);

for i = 1:length(input)
  u = input(i);
  [y, x] = systemModel(A, B, C, D, K, x, u, 0);
  output(i) = y;
end

% Plot
plot(time, measurement, time, output);

function [y, x] = 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

  % x = Ax + Bu + Ke
  x = ...
    A * x + ... % Add contribution of state
    B * u + ... % Add contribution of input
    e * K;      % Add contribution of disturbance
end

[Valeriy]

@Paul Re: "Just copy/paste the code into a comment like you did previously, but include the steps that read the data from the file and generate the identified system. "

I didn't do this in a script. I identified the system by using Matlab UI, see the video I recorded here:

https://youtu.be/A47aR6UAewo

I assume you saw I accepted your answer. The complete working solution is in my previous comment on your accepted answer.

[Torsten]

% x = Ax + Bu + Ke
  x = ...
    A * x + ... % Add contribution of state
    B * u + ... % Add contribution of input
    e * K;      % Add contribution of disturbance
  

So the time stepsize of 0.02 is already included in A and B by the System Identification Toolbox ?

[Paul]

Read in the data and plot it

T = readtable('input.csv')

T = 423x3 table
      time      reading    PWM 
    ________    _______    ____

    0.019995      243         0
    0.040017      243         0
    0.060022      251         0
    0.079997      244         0
     0.10001      228      -111
     0.12001      228      -111
     0.13999      241      -111
     0.16002      232      -111
     0.18002      220      -111
     0.20001      213      -111
     0.22001      214      -159
     0.24009      233      -159
     0.26008      234      -159
     0.28002      234      -159
     0.30002      234      -159
     0.32002      258      -207

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

sysd =
  Discrete-time identified state-space model:
    x(t+Ts) = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)
 
  A = 
              x1         x2         x3
   x1     0.9988    0.05193   -0.02261
   x2    0.02222   -0.01976     0.7353
   x3  0.0009856    -0.2093    -0.5957
 
  B = 
               u1
   x1  -2.663e-06
   x2   5.727e-05
   x3  -0.0001872
 
  C = 
          x1     x2     x3
   y1  -5317  24.87  105.9
 
  D = 
       u1
   y1   0
 
  K = 
               y1
   x1  -0.0001655
   x2    0.001508
   x3   6.209e-06
 
Sample time: 0.02 seconds

Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: none
   Disturbance component: estimate
   Number of free coefficients: 18
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                          
Estimated using N4SID on time domain data "data".
Fit to estimation data: 74.87% (prediction focus)
FPE: 230.3, MSE: 217.6                           
 

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))

ans = 3x1

0.9999 0.4078 0.4078

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

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)

ans = 
  -3.1047e-03 + 0.0000e+00i
  -4.4850e+01 + 1.2140e+02i
  -4.4850e+01 - 1.2140e+02i

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))

ans = 3x1

9.9994e-01 2.4301e+00 2.4301e+00

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

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

[Valeriy]

Wow, what a cool analysis, I'm blown away. Thanks for providing this extra information for anyone interested in this thread.

Answer 2

[Sam Chak]

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

[Valeriy]

Thanks for posting this! I need to do this completely manually, without solvers, using only matrix and vector arithmetic, by looping through the input vector and outputting results into an output vector.

I updated my original question to indicate this, and I also updated the code to more closely matched your implementation. Unfortunately I still can't get it working, as I cannot use ODE3 or Runge-Kutta 4.

[Sam Chak]

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.

[Valeriy]

@Sam Chak my situation is actually the opposite of a homework assignment. I need to get a Kalman filter working to track the position of a robot arm for my drumming robot, but I never studied Calculus or Differential Equations and don't have any engineering background.

I decided to write an article about the Kalman filter as I'm learning about it, because presenting or teaching about a topic forces me to understand it at a great level of depth than simply getting it working and moving on.

If math & engineering are something you feel excited about, I would appreciate feedback. I am writing the article on a branch in the repository of my blog:

https://github.com/01binary/zeroweb/blob/kalman-filter/src/articles/kalman-filter.md

(I haven't incorporated new insights I learned from this thread yet, still testing)