You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Solution of linear-time varying system in matlab
7 views (last 30 days)
Show older comments
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);
19 Comments
Torsten
on 27 Jun 2022
What are your equations written in continuous form that you used for the ODE45 solution ?
ShooQ
on 27 Jun 2022
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
on 27 Jun 2022
Edited: Torsten
on 27 Jun 2022
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
on 27 Jun 2022
Edited: ShooQ
on 27 Jun 2022
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
on 27 Jun 2022
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
on 27 Jun 2022
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
on 27 Jun 2022
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
on 27 Jun 2022
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
on 27 Jun 2022
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
on 27 Jun 2022
ShooQ
on 27 Jun 2022
ShooQ
on 27 Jun 2022
@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
on 27 Jun 2022
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
on 27 Jun 2022
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
on 27 Jun 2022
Edited: Torsten
on 27 Jun 2022
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);
Answers (0)
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)