Solution of linear-time varying system in matlab

I tried the given below to get the solution of the linear-time varying equation using for-loop, but the solution is not right what I am getting using ode45. I don’t know where I am getting wrong while implementing this. Then I want to estimate the parameter g(t) using the gradient method.

r0 = 0.05;
L = 0.1;
d = 0.005;
w0 = 1.5;
Ts=10;
t=[0:Ts:800];
x0=[0 0 1]';
y0=1;
x_value=[];

for k=1:(length(t)) % Number of Iterations
    x_value=[x_value x0];
    
    g(k) =  (2*r0*L*sinh((d*(k))/2))./(d*cosh((d*(k))/2)+L*sinh((d*(k))/2));
    A0 = [ -0.5*g(k)   -w0    0 ; 
        w0     -g(k)*0.5   0 ; 
        0        0     -g(k)];
    B0 =[0; 0; -g(k)];
    Q=integral(@(u) expm(A0*((k+1)*Ts)-u)*B0,(k*Ts),((k+1)*Ts), 'ArrayValued', true);
    
    x0=expm(A0*Ts)*x0+Q;
end
plot(t,x_value(1,:),'r-','linewidth',1);

[Sam Chak]

Can you show the ode45 results?

[ShooQ]

Thanks, please see the attached result

[ShooQ]

This is system equation

[Torsten]

What are your equations written in continuous form that you used for the ODE45 solution ?

[ShooQ]

This is my ode45 code and system equations are attached

function [Tout, Xout] = run_odesq

% Some Parameters

clear all

close all

clc

r0 = 0.05;

L = 0.1;

d = 0.005;

w0 = 1.5;

Ts=0.1;

t=[0:Ts:800];

%Time Varying Elastance

g = @(t) (2*r0*L*sinh((d*t)/2))./(d*cosh((d*t)/2)+L*sinh((d*t)/2));

plot(w0*t,g(t),'LineWidth', 2)

title('\bf \gamma(t)')

xlabel('\omega_0 t')

ylabel('\gamma(t)')

grid on;

% Set up simulation

Ts = 0.1;

dx = @(t,x)ode_sys(Ts,x,g(t));

dy = @(t,x)ode_sys(Ts,x,g(t));

t_end = 800;

t = 0:Ts:t_end;

x0=[0 0 1];

y0=0;

[Tout, Xout] = ode45(dx,t,x0);

figure

subplot(311)

h=plot(Tout, Xout(:,1),'r')

set (h, 'LineWidth', 2);

xlabel(' t')

ylabel('x1')

grid on

subplot(312)

h=plot(Tout, Xout(:,2),'r')

xlabel('t')

ylabel('x2')

set (h, 'LineWidth', 2);

grid on

%

subplot(313)

h=plot(Tout, Xout(:,3),'b')

set (h, 'LineWidth', 2);

xlabel('t')

ylabel('x3')

grid on

end

% Define system equations

function [dx,dy] = ode_sys(Ts,x,g)

w0 = 1.5;

A = [ -0.5*g -w0 0 ;

w0 -g*0.5 0 ;

0 0 -g];

B = [0; 0; -g];

C=[0 0 1];

dx = A*x + B;

dy=C*x;

end

[Torsten]

I didn't check whether your numerical method using matrix exponentials is correct to give the solution of the ODE system

dx/dt = A*x + b, x(0) = [0 0 1]

with

r0 = 0.05;

L = 0.1;

d = 0.005;

w0 = 1.5;

A = [ -0.5*g(t) -w0 0 ;

w0 -g(t)*0.5 0 ;

0 0 -g(t)];

B = [0; 0; -g(t)];

g(t) = (2*r0*L*sinh((d*t)/2))./(d*cosh((d*t)/2)+L*sinh((d*t)/2))

I doubt it since the results are too different.

[T,X] = run_odesq();

h = 
  Line with properties:

              Color: [1 0 0]
          LineStyle: '-'
          LineWidth: 0.5000
             Marker: 'none'
         MarkerSize: 6
    MarkerFaceColor: 'none'
              XData: [0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2 2.1000 2.2000 2.3000 2.4000 2.5000 2.6000 2.7000 2.8000 2.9000 3 3.1000 3.2000 3.3000 3.4000 … ]
              YData: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 … ]

  Show all properties

h = 
  Line with properties:

              Color: [1 0 0]
          LineStyle: '-'
          LineWidth: 0.5000
             Marker: 'none'
         MarkerSize: 6
    MarkerFaceColor: 'none'
              XData: [0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2 2.1000 2.2000 2.3000 2.4000 2.5000 2.6000 2.7000 2.8000 2.9000 3 3.1000 3.2000 3.3000 3.4000 … ]
              YData: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 … ]

  Show all properties

h = 
  Line with properties:

              Color: [0 0 1]
          LineStyle: '-'
          LineWidth: 0.5000
             Marker: 'none'
         MarkerSize: 6
    MarkerFaceColor: 'none'
              XData: [0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2 2.1000 2.2000 2.3000 2.4000 2.5000 2.6000 2.7000 2.8000 2.9000 3 3.1000 3.2000 3.3000 3.4000 … ]
              YData: [1 1.0000 0.9998 0.9996 0.9992 0.9988 0.9982 0.9976 0.9969 0.9961 0.9952 0.9942 0.9931 0.9919 0.9907 0.9893 0.9879 0.9864 0.9848 0.9831 0.9813 0.9795 0.9776 0.9756 0.9735 0.9713 0.9691 0.9668 0.9644 0.9620 0.9595 0.9569 0.9542 … ]

  Show all properties

function [Tout, Xout] = run_odesq
% Some Parameters
clear all
close all
clc
r0 = 0.05;
L = 0.1;
d = 0.005;
w0 = 1.5;
Ts=0.1;
t=[0:Ts:800];
%Time Varying Elastance
g = @(t)  (2*r0*L*sinh((d*t)/2))./(d*cosh((d*t)/2)+L*sinh((d*t)/2));
plot(w0*t,g(t),'LineWidth', 2)
title('\bf \gamma(t)')
xlabel('\omega_0 t')
ylabel('\gamma(t)')
grid on;
% Set up simulation
%Ts = 0.1;
fun =@(t,x)fun_ode(t,x,g);
%dx = @(t,x)ode_sys(Ts,x,g(t));
%dy = @(t,x)ode_sys(Ts,x,g(t));
t_end = 800;
t = 0:Ts:t_end;
x0=[0 0 1];
%y0=0;
[Tout, Xout] = ode45(fun,t,x0);
figure
subplot(311)
h=plot(Tout, Xout(:,1),'r')
set (h, 'LineWidth', 2); 
xlabel(' t')
ylabel('x1')
grid on
subplot(312)
h=plot(Tout, Xout(:,2),'r')
xlabel('t')
ylabel('x2')
set (h, 'LineWidth', 2); 
grid on
% 
subplot(313)
h=plot(Tout, Xout(:,3),'b')
set (h, 'LineWidth', 2); 
xlabel('t')
ylabel('x3')
grid on
end
% Define system equations
function [dx] = fun_ode(t,x,g)
w0 = 1.5;
A = [ -0.5*g(t)   -w0     0 ; 
    w0     -g(t)*0.5   0 ; 
    0        0     -g(t)];
B = [0; 0; -g(t)];
%C=[0 0 1];
dx = A*x + B;
%dy=C*x;
end

[ShooQ]

I am also confused, I think the result of numerical method using matrix exponentials is not correct but with ODE45 is giving the right result. I think the I am missing something or the code is not right. But with ode45 the result are same what i want to.

[Sam Chak]

Yup, something is wrong. Most likely the integral part that you are trying to solve the continuous-time state matrix in the discrete-time manner.

% parameters
r0 = 0.05;
L  = 0.1;
d  = 0.005;
w0 = 1.5;
Ts = 0.1;
t  = 0:Ts:800;
x0 = [0 0 1]';

% time-varying gamma function
gamma = @(t)  (2*r0*L*sinh(d*t/2))./(d*cosh(d*t/2) + L*sinh(d*t/2));
plot(t, gamma(t),'LineWidth', 2), grid on,
xlabel('t')
ylabel('\gamma(t)')

x_value = [];
for k = 1:(length(t))   % Number of Iterations
    x_value = [x_value x0];
    
    g(k) = (2*r0*L*sinh(d*(k)/2))./(d*cosh(d*(k)/2) + L*sinh(d*(k)/2));
    A0   = [-0.5*g(k) -w0 0; 
             w0 -g(k)*0.5 0; 
              0    0   -g(k)];
    B0   = [0; 0; -g(k)];
    Q    = integral(@(u) expm(A0*((k + 1)*Ts) - u)*B0, (k*Ts), ((k + 1)*Ts), 'ArrayValued', true);
    
    x0   = expm(A0*Ts)*x0 + Q;
end

subplot(311)
plot(t, x_value(1, :), 'r', 'LineWidth', 2)
xlabel('t')
ylabel('x1')
grid on

subplot(312)
plot(t, x_value(2, :), 'r', 'LineWidth', 2)
xlabel('t')
ylabel('x2')
grid on

subplot(313)
plot(t, x_value(3, :), 'b', 'LineWidth', 2)
xlabel('t')
ylabel('x3')
grid on

[Torsten]

In your matrix exponential version you define

g(k) =  (2*r0*L*sinh((d*(k))/2))./(d*cosh((d*(k))/2)+L*sinh((d*(k))/2));

Shouldn't this be

g(k) =  (2*r0*L*sinh((d*(t(k)))/2))./(d*cosh((d*(t(k)))/2)+L*sinh((d*(t(k)))/2));

?

[ShooQ]

What I have shown in equation they are descritzed and I want to slove the system in descrete-time manner using numerical exponential matrix method.

[ShooQ]

g(k) = (2*r0*L*sinh((d*(t(k)))/2))./(d*cosh((d*(t(k)))/2)+L*sinh((d*(t(k)))/2)); Even though it doesn't work with this.

[Torsten]

What I have shown in equation they are descritzed and I want to slove the system in descrete-time manner using numerical exponential matrix method.

I know.

But how is X(k) (your discrete time method with matrix exponential) related to X(t(k)) (solution of ODE45 at time t(k)) ?

They cannot be the same since k is not t(k).

[Sam Chak]

Hi @ShooQ

@Torsten is right. The discrete-time system is not the same as the continuous-time system. You need to discretizes the continuous-time system using zero-order hold on the inputs and with a sample rate.

[ShooQ]

@Torsten@Sam Chak Thank so much for prompt reply and suggestions. Yes you both are right but I am confused how I can do in my code (discrete time method with matrix exponential)? Beacuse I am not so good in coding.

[Torsten]

It's not a coding problem, but a problem of how your discrete formulation of the problem (sys.png) is to be transformed into a continuous formulation (for ODE45).

I can't help you in this respect, sorry.

[ShooQ]

@Torsten Thankyou for your precious time. Sir, I don't want the result of ODE45 and I only want the result of descrete system (sys.png). But the response of discrete system is not right what I want to. Do you think the code (discrete time method with matrix exponential) is right to simulate (sys.png)? I have doubt in this. Thanks

[Torsten]

I don't have doubts, but I'm sure that your definition of g must be wrong.

Your g and thus your A0 and b0 must depend on t and not on an arbitrary loop index k. So my guess would be that you have to replace

g(k) =  (2*r0*L*sinh((d*(k))/2))./(d*cosh((d*(k))/2)+L*sinh((d*(k))/2));

by

g(k) =  (2*r0*L*sinh((d*(t(k)))/2))./(d*cosh((d*(t(k)))/2)+L*sinh((d*(t(k)))/2));

as I already suggested.

But I have no experience with discrete time evolution systems - so I can't tell with certainty.

[ShooQ]

This is what I want, I want to to update A0 and B0 like A0(k) or A0(t(k)) and same as B0(k) or B0(t(k)) so that it update iteratively. I tried but it doesn't work, e.g.,

r0 = 0.05;

L = 0.1;

d = 0.005;

w0 = 1.5;

Ts=10;

t=[0:Ts:800];

x0=[0 0 1]';

y0=1;

A0=zeros(3,3);

B0=zeros(3,1);

x_value=[];

A_value=[];

B_value=[];

for k=1:(length(t)) % Number of Iterations

x_value=[x_value x0];

A_value=[A_value A0];

B_value=[B_value B0];

for k=1:(length(t)) % Number of Iterations

x_value=[x_value x0];

g(k) = (2*r0*L*sinh((d*(k))/2))./(d*cosh((d*(k))/2)+L*sinh((d*(k))/2));

A0(k) = [ -0.5*g(k) -w0 0 ;

w0 -g(k)*0.5 0 ;

0 0 -g(k)];

B0(k) =[0; 0; -g(k)];

Q(k)=integral(@(u) expm(A0(k)*((k+1)*Ts)-u)*B0(k),(k*Ts),((k+1)*Ts), 'ArrayValued', true);

x0(k+1)=expm(A0*Ts)*x0(k)+Q(k);

end

plot(t,x_value(1),'r-','linewidth',1);

[Torsten]

If these are the correct update formulae, the original code you posted was correct. You can't write A0(k)=..., B0(k) = ...; x0(k+1) = ... since the objects on the right-hand side (the ...) aren't scalar values, but vectors or matrices.

This code works, but must be wrong for the reason I already mentionned (g only depends on the loop index, but not on t):

r0 = 0.05;
L = 0.1;
d = 0.005;
w0 = 1.5;
Ts = 10;
t = 0:Ts:800;
x0=[0 0 1]';
A0=zeros(3,3);
B0=zeros(3,1);
x_value=[];
for k=1:length(t) % Number of Iterations
    x_value=[x_value x0];
    g =  (2*r0*L*sinh((d*(k))/2))./(d*cosh((d*(k))/2)+L*sinh((d*(k))/2));
    A0 = [ -0.5*g   -w0    0 ; 
           w0     -g*0.5   0 ; 
           0        0     -g];
    B0 =[0; 0; -g];
    Q = integral(@(u) expm(A0*((k+1)*Ts)-u)*B0,(k*Ts),((k+1)*Ts), 'ArrayValued', true);    
    x0 = expm(A0*Ts)*x0+Q;
end
plot(t,x_value(3,:),'r-','linewidth',1);