How to manually simulate linear system without lsim

Question

Valeriy 2024-4-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2112301-how-to-manually-simulate-linear-system-without-lsim

评论： Valeriy 2024-4-28

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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Paul 2024-4-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2112301-how-to-manually-simulate-linear-system-without-lsim#answer_1448931

在 MATLAB Online 中打开

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 个评论
显示 13更早的评论隐藏 13更早的评论

Valeriy 2024-4-27

在 MATLAB Online 中打开

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 2024-4-27

编辑：Valeriy 2024-4-27

在 MATLAB Online 中打开

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

Valeriy 2024-4-28

编辑：Valeriy 2024-4-28

在 MATLAB Online 中打开

input.csv

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

Paul 2024-4-28

在 MATLAB Online 中打开

input.csv

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 2024-4-28

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 2024-4-27

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2112301-how-to-manually-simulate-linear-system-without-lsim#answer_1448801

在 MATLAB Online 中打开

input.csv

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 个评论
显示 1更早的评论隐藏 1更早的评论

Sam Chak 2024-4-28

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 2024-4-28

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

请先登录，再进行评论。