現在この質問をフォロー中です
- フォローしているコンテンツ フィードに更新が表示されます。
- コミュニケーション基本設定に応じて電子メールを受け取ることができます。
Solution of linear-time varying system in matlab
15 ビュー (過去 30 日間)
古いコメントを表示
ShooQ
2022 年 6 月 27 日
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 件のコメント
Torsten
2022 年 6 月 27 日
What are your equations written in continuous form that you used for the ODE45 solution ?
ShooQ
2022 年 6 月 27 日
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
2022 年 6 月 27 日
編集済み: Torsten
2022 年 6 月 27 日
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
2022 年 6 月 27 日
編集済み: ShooQ
2022 年 6 月 27 日
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
2022 年 6 月 27 日
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
2022 年 6 月 27 日
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
2022 年 6 月 27 日
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
2022 年 6 月 27 日
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
2022 年 6 月 27 日
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
2022 年 6 月 27 日
ShooQ
2022 年 6 月 27 日
ShooQ
2022 年 6 月 27 日
@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
2022 年 6 月 27 日
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
2022 年 6 月 27 日
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
2022 年 6 月 27 日
編集済み: Torsten
2022 年 6 月 27 日
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);
回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!エラーが発生しました
ページに変更が加えられたため、アクションを完了できません。ページを再度読み込み、更新された状態を確認してください。
Web サイトの選択
Web サイトを選択すると、翻訳されたコンテンツにアクセスし、地域のイベントやサービスを確認できます。現在の位置情報に基づき、次のサイトの選択を推奨します:
また、以下のリストから Web サイトを選択することもできます。
最適なサイトパフォーマンスの取得方法
中国のサイト (中国語または英語) を選択することで、最適なサイトパフォーマンスが得られます。その他の国の MathWorks のサイトは、お客様の地域からのアクセスが最適化されていません。
南北アメリカ
- América Latina (Español)
- Canada (English)
- United States (English)
ヨーロッパ
- 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)
アジア太平洋地域
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)