State space modeling of LTI

I have a LTI system is given in the attached picture. For which I have A=[0 1;1 2] B=[1;1], C=[3,4], D=0.1 and L=[0 1]. I need to find z(t) over an interval [0,20] and initial condition is assumed as x0=[1;-1] with an input w=0.5sin(2pi)t and v=t, how can I make a code for this??

Answer 1

[Ameer Hamza]

Such LTI system can be solved using lsim

C=[3,4]; 
D=0.1;
L=[0 1];
sys = ss(A,B,C,D);

t = linspace(0,20,1000)';
x0 = [1; -1];
[t,x] = ode45(@odeFun, t, x0);

v = t;
y = C*x' + D*v';

z = L*x';
plot(t,z)

function dxdt = odeFun(t, x)
    A=[0 1;1 2];
    B=[1;1];
    w = 0.5*sin(2*pi*t);
    
    dxdt = A*x+B*w;
end

However, the system is unstable and the output diverges to infinity.

[Muhammad Atif]

thanks sir, but how can I use initial condition here and what's about disturbance 'v'??

[Ameer Hamza]

Ok. I just noticed that your LTI system is not of a standard form. So regular tools for LTI systems will not work here. I gave a solution by solving the differential equations. Please check the updated answer. It also considers disturbance v and the initial conditions.

[Muhammad Atif]

ok, let me try this, thanks for your help

[Ameer Hamza]

Glad to be of help.

[Muhammad Atif]

Hi Sir, I modified the code that you created above, In this I added the input w(t) and noise signal v(t). The plot of (z=L*v' ) should start when input starts to take effects. I attached two images for references. This plot I obtained from this code, while I wanna start it from t=5. How can I change this?

%%function to get dxdt
function dxdt = odeFun(t, x)
ms = 973;
ks = 42720;
cs = 3000;
ku = 101115;
mu = 114;
G = 512*10^(-6)
q = 0.1;
v1 = 12.5;
A = [0 0 1 -1 ; 0 0 0 1 ; -ks/ms 0 -cs/ms  cs/ms; ks/mu -ku/mu cs/mu -cs/mu];
B = [0 0 -2*pi*q*(G*v1)^(1/2) 0]';
a = 0.2 ;
l = 5 ;
%%input w(t)=1 for 5<=t<=10 and w(t)=-1 for 15<=t<=20
w_int1 = t>=5 & t<=10; 
w_int2 = t>=15 & t<=20;
w = zeros(1,numel(t));%add this
rng(0);
w(w_int1)=1;
w(w_int2)=-1;
dxdt = A*x+B*w';
end

C = [0 0 0 1];
D = 0.1;
L = [0 0 1 0];
t = linspace(0,30,1000)';
x0 = [0.1; 0; -0.1; -0.2];
[t,x] = ode45(@odeFun, t, x0);
T = linspace(0,20,1000)';
%%noise v(t) distributed over an interval [0,20]
v_int = t>=0 & t<=20; 
v = zeros(1,numel(t)); 
rng(0);
v(v_int)= 0.2*rand(1,length(t(v_int)));
y = C*x' + D*v;
z = L*x';
plot(t,z)

[Ameer Hamza]

Start from initial condition of zero to see the effect. Alsi it is better to use if-block to calculate w based on t

C = [0 0 0 1];
D = 0.1;
L = [0 0 1 0];
t = linspace(0,30,1000)';
x0 = [0; 0; -0; -0];
[t,x] = ode45(@odeFun, t, x0);
T = linspace(0,20,1000)';
%%noise v(t) distributed over an interval [0,20]
% v_int = t>=0 & t<=20; 
% v = zeros(1,numel(t)); 
% rng(0);
% v(v_int)= 0.2*rand(1,length(t(v_int)));
% y = C*x' + D*v;
% z = L*x';
plot(t,x)
    
%%function to get dxdt
function dxdt = odeFun(t, x)
ms = 973;
ks = 42720;
cs = 3000;
ku = 101115;
mu = 114;
G = 512*10^(-6);
q = 0.1;
v1 = 12.5;
A = [0 0 1 -1 ; 0 0 0 1 ; -ks/ms 0 -cs/ms  cs/ms; ks/mu -ku/mu cs/mu -cs/mu];
B = [0 0 -2*pi*q*(G*v1)^(1/2) 0]';
a = 0.2 ;
l = 5 ;
%%input w(t)=1 for 5<=t<=10 and w(t)=-1 for 15<=t<=20
w = 0;
if t>=5 && t<=10
    w = 1;
elseif t>=15 && t<=20
    w = -1;
end
dxdt = A*x+B*w';
end

[Muhammad Atif]

Hi sir, As in the perviously question we found z(t) and y(t); Now in this question I need to use the y(t) as the in put of this state-space equation for the given Af = [-2.8682 0.5240 0.1093 0.4576; 0.6676 -2.3337 -0.1285 -0.5333; 4.2978 1.3131 -13.1636 3.3939; -30.1511 55.6790 48.9362 -218.6780];

Bf = [0.0167; -0.0400; -0.1676; -3.6279];

Cf = [1.6042 1.1359 -0.7995 -0.1854];

and xf=[0;0;0;0] is the intial state for this system

and the state space equation is given in the picture

[Ameer Hamza]

You can use lsim for this system

Af = [-2.8682 0.5240 0.1093 0.4576; 0.6676 -2.3337 -0.1285 -0.5333; 4.2978 1.3131 -13.1636 3.3939; -30.1511 55.6790 48.9362 -218.6780];
Bf = [0.0167; -0.0400; -0.1676; -3.6279];
Cf = [1.6042 1.1359 -0.7995 -0.1854];
Df = 0;

sys = ss(Af,Bf,Cf,Df);
lsim(sys,y_hat,t); % t should be a column vector and y_hat should have same number of rows as t and 2 columns

[Muhammad Atif]

Af = [-2.8682 0.5240 0.1093 0.4576; 0.6676 -2.3337 -0.1285 -0.5333; 4.2978 1.3131 -13.1636 3.3939; -30.1511 55.6790 48.9362 -218.6780];
Bf = [0; -0; -0; 0];
Cf = [1.6042 1.1359 -0.7995 -0.1854];
Df = 0;
sys = ss(Af,Bf,Cf,Df);
y_hat = zeros(size(t));
lsim(sys,y_hat,t)

why this system gives zero output when I choose y_hat=0??

[Ameer Hamza]

Because the input to the system is zero so the states cannot change. You need to give some non-zero input to change the state.

[Muhammad Atif]

Sorry to bother you again. Now I have no idea how to use the the state of xf0 for the other system. It's a sort of switching condition in which two subsystems. When 't' reaches the defined interval the system one will run, otherwise the 2nd system will run. But how can I use the state xf0 of the first system for system 2 when system 1 run for the last time??

clc
clear all
Af1 = [-2.8682 0.5240 0.1093 0.4576; 0.6676 -2.3337 -0.1285 -0.5333; 4.2978 1.3131 -13.1636 3.3939; -30.1511 55.6790 48.9362 -218.6780];
Bf1 = [0.0167; -0.0400; -0.1676; -3.6279];
Bf2 = [0;0;0;0];
Af2 = [-2.8682 0.5240 0.1093 0.4576; 0.6676 -2.3337 -0.1285 -0.5333; 4.2978 1.3131 -13.1636 3.3939; -30.1511 55.6790 48.9362 -218.6780];
Cf1 = [1.6042 1.1359 -0.7995 -0.1854];
Cf2 = [0.9791 0.7264 -0.6448 -0.0632];
Df1=0;
Df2=0;
C = [0 0 0 1];
D = 0.1;
L = [0 0 1 0];
t = linspace(0,30,1000);
x0 = [0; 0; -0.00; -0.00];
[t,x] = ode45(@odeFun12, t, x0);
%%noise v(t) distributed over an interval [0,20]
v = 0.05*rand(size(t));
v = v.*(t<20);
y = C*x' + D*v';
z = L*x';
%plot(t,y)
%grid on
%%system 1
if t<=1 | (t>=2 & t<=7) | (t>=12 & t<=17) | (t>=22 & t<=27)
sys = ss(Af1,Bf1,Cf1,Df1);
  xf0 = [0; 0; 0.00; 0.00];
  lsim(sys,y,t,xf0);
%%system 2 I need to use xf0 of the 1st system when program switches to system 2 
else
    sys1 = ss(Af2,Bf2,Cf2,Df2);
    y1=zeros(size(t));
    lsim(sys1,y1,t,xf0);
end

[Ameer Hamza]

What is the definition of the function odeFun12?

[Muhammad Atif]

function dxdt = odeFun12(t, x)
ms = 973;
ks = 42720;
cs = 3000;
ku = 101115;
mu = 114;
G = 512*10^(-6);
q = 0.1;
v1 = 12.5;
A = [0 0 1 -1; 0 0 0 1; -ks/ms 0 -cs/ms  cs/ms; ks/mu -ku/mu cs/mu -cs/mu];
B = [0 -2*pi*q*(G*v1)^(1/2) 0 0]';
a = 0.2 ;
l = 25;
%%input w = a*pi*v/l.*sin(2*pi*v/l.*t) for 5<=t<=l/v and zero otherwise, while 'v' is uniformely distributed noise between[0,0.1] for an interval 0 to 20

w = zeros(size(t));
mask = t>=5 & t<=5+l/v1;
w = a*pi*v1/l.*sin(2*pi*v1*t/l).*mask;
dxdt = A*x+B*w';
end

[Muhammad Atif]

It's the same as before

[Ameer Hamza]

I guess something like this will work

Cf1 = [1.6042 1.1359 -0.7995 -0.1854];
Cf2 = [0.9791 0.7264 -0.6448 -0.0632];
Df1=0;
Df2=0;

t = linspace(0,30,1000);
x0 = [0; 0; -0.00; -0.00];
[t,x] = ode45(@odeFun12, t, x0);
%%noise v(t) distributed over an interval [0,20]
mask = t<=1 | (t>=2 & t<=7) | (t>=12 & t<=17) | (t>=22 & t<=27);
y(mask) = Cf1*x(mask,:).';
y(~mask) = Cf2*x(~mask,:).';


function dxdt = odeFun12(t, x)

Af1 = [-2.8682 0.5240 0.1093 0.4576; 0.6676 -2.3337 -0.1285 -0.5333; 4.2978 1.3131 -13.1636 3.3939; -30.1511 55.6790 48.9362 -218.6780];
Bf1 = [0.0167; -0.0400; -0.1676; -3.6279];
Bf2 = [0;0;0;0];
Af2 = [-2.8682 0.5240 0.1093 0.4576; 0.6676 -2.3337 -0.1285 -0.5333; 4.2978 1.3131 -13.1636 3.3939; -30.1511 55.6790 48.9362 -218.6780];

ms = 973;
ks = 42720;
cs = 3000;
ku = 101115;
mu = 114;
G = 512*10^(-6);
q = 0.1;
v1 = 12.5;
a = 0.2 ;
l = 25;
%%input w = a*pi*v/l.*sin(2*pi*v/l.*t) for 5<=t<=l/v and zero otherwise, while 'v' is uniformely distributed noise between[0,0.1] for an interval 0 to 20

w = a*pi*v1/l.*sin(2*pi*v1*t/l);

if  t<=1 || (t>=2 && t<=7) || (t>=12 && t<=17) || (t>=22 && t<=27)
    dxdt = Af1*x+Bf1*w';
else
    dxdt = Af2*x+Bf2*w';
end
end

[Muhammad Atif]

I think I didn't explain well and you couldn't get me. y(t) is already calculated from the ode fuction. In the last question I need to use y(t) as an input for the next two systems. In which y(t) for the 2nd system is zero, while the first system will use the y(t) calculated before. lsim(sys,y,t,xf0),as you see 'y' is used as an input for this system. y1=zeros(size(t)),lsim(sys1,y1,t,xf0); in this it is zero, but this system will use the last state of first system when it will run

[Ameer Hamza]

I think something like this will work. You can extend it to the full time duration.

clc
clear all
Af1 = [-2.8682 0.5240 0.1093 0.4576; 0.6676 -2.3337 -0.1285 -0.5333; 4.2978 1.3131 -13.1636 3.3939; -30.1511 55.6790 48.9362 -218.6780];
Bf1 = [0.0167; -0.0400; -0.1676; -3.6279];
Bf2 = [0;0;0;0];
Af2 = [-2.8682 0.5240 0.1093 0.4576; 0.6676 -2.3337 -0.1285 -0.5333; 4.2978 1.3131 -13.1636 3.3939; -30.1511 55.6790 48.9362 -218.6780];
Cf1 = [1.6042 1.1359 -0.7995 -0.1854];
Cf2 = [0.9791 0.7264 -0.6448 -0.0632];
Df1=0;
Df2=0;
C = [0 0 0 1];
D = 0.1;
L = [0 0 1 0];
t = linspace(0,30,1000);
x0 = [0; 0; -0.00; -0.00];
[t,x] = ode45(@odeFun12, t, x0);
%%noise v(t) distributed over an interval [0,20]
v = 0.05*rand(size(t));
v = v.*(t<20);
y = C*x' + D*v';
z = L*x';
%plot(t,y)
%grid on
%%system 1
% if t<=1 | (t>=2 & t<=7) | (t>=12 & t<=17) | (t>=22 & t<=27)
sys1 = ss(Af1,Bf1,Cf1,Df1);
sys2 = ss(Af2,Bf2,Cf2,Df2);
xf0 = [0; 0; 0.00; 0.00];
t1 = t(t<=1);
y1 = y(t<=1);
[~,Y1,X1] = lsim(sys1,y1,t1,xf0);
t2 = t(1<t & t<=2);
y2 = 0*y(1<t & t<=2);
[~,Y2,X2] = lsim(sys2,y2,t2,X1(end,:));
t3 = t(2<t & t<=7);
y3 = y(2<t & t<=7);
[~,Y3,X3] = lsim(sys1,y3,t3,X2(end,:));
t4 = t(7<t & t<=12);
y4 = 0*y(7<t & t<=12);
[~,Y4,X4] = lsim(sys2,y4,t4,X3(end,:));

Y = [Y1;Y2;Y3;Y4];
X = [X1;X2;X3;X4];

function dxdt = odeFun12(t, x)
ms = 973;
ks = 42720;
cs = 3000;
ku = 101115;
mu = 114;
G = 512*10^(-6);
q = 0.1;
v1 = 12.5;
A = [0 0 1 -1; 0 0 0 1; -ks/ms 0 -cs/ms  cs/ms; ks/mu -ku/mu cs/mu -cs/mu];
B = [0 -2*pi*q*(G*v1)^(1/2) 0 0]';
a = 0.2 ;
l = 25;
%%input w = a*pi*v/l.*sin(2*pi*v/l.*t) for 5<=t<=l/v and zero otherwise, while 'v' is uniformely distributed noise between[0,0.1] for an interval 0 to 20

w = zeros(size(t));
mask = t>=5 & t<=5+l/v1;
w = a*pi*v1/l.*sin(2*pi*v1*t/l).*mask;
dxdt = A*x+B*w';
end

[Muhammad Atif]

That's it, finally I got the results. Thank you very much for your help

[Ameer Hamza]

I am glad to be of help.