# why are the output simulation results of the following LTI system different by using 'lsim' command and theorical method calculations?

2 views (last 30 days)
sara elmi on 17 Dec 2021
Commented: sara elmi on 20 Dec 2021
Dear all,
I would like to obtain the output response of an LTI state space model to an exponential input compatible with system invariant zero frequency and adjusted amplitude to get a zero oputput in all the time (input signal : u(t)=-2*exp(-3*t); initial state vector =x0=[1;-1;0];)
so I wrote the following code by two different ways :
1-by using 'lsim 'command
2- by using output response theoretical formulation.
unfortunatley ,the system output results weren't compatible with each other. what is the reason for this difference?is there any limitation on using "lsim "command for LTI systems?
clc
clear all
A=[1 4 0;0 -1 0;0 2 -3]
B=[0;-1;-1]
C=[-1 -1 0];
D=0;
sys1=ss(A,B,C,D);
G1=tf(sys1)
time=0:0.1:10;
for i=1:length(time)
u(i)=-2*exp(-3*time(i));
end
x0=[1;-1;0]; %% initial state vector
% % % % % % % % % % % % % % % %simulation 1
y=lsim(sys1,u',time,x0);
figure(1)
plot(time,y)
xlabel('time')
ylabel('output')
title('output response to exponetial inout(simulation 1)')
grid on
% % % % % % % % % % % % % % % %simulation 2
syms t to
expm(A*t)
y3=(C*expm(A*t)*(x0))+(C*expm(A*t)*(int(expm(-A*to)*B*-2*exp(-3*to),0,t)))
for i=1:length(time)
y3n(i)=subs(y3,t,time(i));
end
y3t=double(y3n)
figure(2)
plot(time,y3t,'r');
xlabel('time')
ylabel('output')
title('output response to exponetial inout(simulation 2)')
grid on

Paul on 17 Dec 2021
Edited: Paul on 17 Dec 2021
I suspect the problem is that lsim() is ineherently an approximation to the LTI system response. I think that lsim() generates a discrete time approximation, and so its accuracy will be a function of the time step. And of course there will always be numerical inaccuracy. I'm pretty sure, but not certain, that lsim computes the output due to the input and the initial conditions separately, and then adds them together at the end. For this example, it would be subtracting two very large numbers (because the system is unstable) in the hopes of getting zero as the net output.
Let's look at the code:
A=[1 4 0;0 -1 0;0 2 -3];
B=[0;-1;-1];
C=[-1 -1 0];
D=0;
sys1=ss(A,B,C,D);
G1=tf(sys1);
time = 0:0.1:10;
u = -2*exp(-3*time);
x0 = [1;-1;0]; %% initial state vector
% % % % % % % % % % % % % % % %simulation 1
y = lsim(sys1,u,time,x0);
figure
plot(time,y)
xlabel('time')
ylabel('output')
title('output response to exponential inout(simulation 1)')
grid on Now try it with a much smaller time step
time = 0:0.0001:10;
u = -2*exp(-3*time);
x0 = [1;-1;0]; %% initial state vector
% % % % % % % % % % % % % % % %simulation 1
y = lsim(sys1,u,time,x0);
figure
plot(time,y)
xlabel('time')
ylabel('output')
grid on Although the shape of the curve is the same, note the scale on the y-axis where the order of magnitude is much lower.
% % % % % % % % % % % % % % % %simulation 2
syms t to
expm(A*t);
y3(t) =(C*expm(A*t)*(x0))+(C*expm(A*t)*(int(expm(-A*to)*B*-2*exp(-3*to),0,t)));
y3(t)
ans = which can be simplified to
y3(t) = simplify(y3(t))
y3(t) =
0
which is the expected result. The plot with the un-simplified version of y3(t) was jsut suffering from numerical inaccuracy.
sara elmi on 20 Dec 2021
Dear Paul,
I really appreciate that you took the time to answer my question. It was comprehensive and exactly what I needed to know. Thank you so much.