Solve differnetial equation of state space

Hello everyone,

I want to solve a differential equation of state space dx/dt = A*x + B*u, so I have a simple model with a Ground movement.

My equation is y'' = -d/m * y' - k/m * y + d/m * u' + k/m * u. I have the y (Output) and I want to find u, thats the Input to my system.

To have a first order equation i wrote dy(1) = y(2) and than dy(2) = -(d/m)*y(1)-(k/m)*y + (d/m)*u(1) + (k/m)*u(2)

m = 50;
d = 1000;
k = 30000;
A = [0,1;-k/m,-d/m];
B = [k/m;d/m];
C = [1,0];
D = [0];

Thats the code I've tried but it doesent function.

function du = fun(u,y)
dy = zeros(2,1)
dy(1) = y(2)
du = zeros(1,1)
du(1) = u2
%dy(2) = -(d/m)*y(1)-(k/m)*y + (d/m)*u(1) + (k/m)*u(2)
du(1) = -(m/d)*dy(2)-y(1)-(k/d)*y + (k/d)*u(1)

 end
tspan = [0 10]
y0 = 1;
[u,y] = ode45(@(u,y) tspan, u0);

[Ameer Hamza]

Thats not how we usually solve differential equations. Normally, u is given, and the goal is to find 'y'. Can you give an example of how vector y is given?

[Bleron Kroni]

I have a simple state space model, Mass,Dämping and Spring, so I know the mass and dämping and stifness coefficients.

The equation that describes the system ist this one x'' = -d/m * x' - k/m * x + d/m * u' + k/m * u, where u ist the input(movement from ground) of the system and y = c*x, where c ist a matrix [1 0]. y ist given and its a vector 1001x1. I can substitute the y to the x so i get this equation when i substitute. y'' = -d/m * y' - k/m * y + d/m * u' + k/m * u. The target is to find u. Maybe the differential equations is not the right way. I thought I could solve it in this way, but i think I have to calculate the transfer function or I have to move to the freqequency domain.

The other idea ist to make a dekonvolution process, between the impulse response of the system and the output y.

[Ameer Hamza]

I misinterpreted the question initially. Try the code in my answer.

Answer 1

[Ameer Hamza]

try this code. It assumes that y is know and u is unknown. I set y(t) to sin(t) as an example.

time = linspace(0,10,1000);
y = sin(time);
dy = gradient(y,time);
ddy = gradient(dy,time);

u0 = 0;
[t, u] = ode45(@(t,u) odeFun(t,u,time,y,dy,ddy), time, u0);
plot(t,u);

function dudt = odeFun(t, u, time, y, dy, ddy)
    m = 50;
    d = 1000;
    k = 30000;
    
    y = interp1(time, y, t);
    dy = interp1(time, dy, t);
    ddy = interp1(time, ddy, t);
    
    dudt = m/d*(ddy + d/m*dy + k/m*y - k/m*u);
end

[Bleron Kroni]

Hi Ameer,

i have the state space model with Matrix A,B,C,D and as inuput I give a gaussian white noise r and I get an output Y (Y = lsim(sysd,r,t)), so the code that you wrote I want to give the output of the system Y and as output u is going to be r. The idea ist to invert the system in someway.

m = 50;

d = 1000;

k = 30000;

A = [0,1;-k/m,-d/m];

B = [k/m;d/m];

C = [1,0];

D = [0];

Ts = 0.01;

sysc = ss(A,B,C,D);

[b,a] = ss2tf(A,B,C,D);

sysd = c2d(sysc,Ts);

A1 = sysd.A;

B1 = sysd.B;

C1 = sysd.C;

D1 = sysd.D;

[nu,de] = ss2tf(A1,B1,C1,D1);

mySys_tf=tf(nu,de)

%H = impulse(sysd);

t = 0:0.01:10;

mean = 0; %mean

Ts=0.01; %Given sampling period

SNRdb = 75; %Given SNRdb

variance = 1/(Ts*(10^(SNRdb/10)));

%variance = 2;

standarddeviation = sqrt(variance);

r = mean + standarddeviation*randn(1001,1);

var(r);

Y = lsim(sysd,r,t);

time = 0:0.01:10;

time = transpose(time);

y = Y;

dy = gradient(y,time);

ddy = gradient(dy,time);

u0 = 0;

[t, u] = ode45(@(t,u) odeFun(t,u,time,y,dy,ddy), time, u0);

plot(t,u);

function dudt = odeFun(t, u, time, y, dy, ddy)

m = 50;

d = 1000;

k = 30000;

y = interp1(time, y, t);

dy = interp1(time, dy, t);

ddy = interp1(time, ddy, t);

dudt = m/d*(ddy + d/m*dy + k/m*y - k/m*u);

end