Solving an integral when the parameters is a matrix - Controllability Gramian

2 次查看（过去 30 天）
mpz
mpz 2022-11-28
评论： Paul 2022-11-29
Ran in:
Hi,
I need solving the controllability gramian below given B and psi
Psi is calculated as below as phi_s. B matrix is give below. [to tf]=[0 10]
clear all;clc;close all
syms t;
A = [0 1 0 0;0 0 -4.91 0;0 0 0 1;0 0 73.55 0];
B = [0;2;0;-10];
[M,J] = jordan(A);
phi_s=M * expm(J*t) * inv(M);
figure(3)
fplot([phi_s(1,1) phi_s(1,2) phi_s(1,3) phi_s(1,4) ...
phi_s(2,1) phi_s(2,2) phi_s(2,3) phi_s(2,4) ...
phi_s(3,1) phi_s(3,2) phi_s(3,3) phi_s(3,4) ...
phi_s(4,1) phi_s(4,2) phi_s(4,3) phi_s(4,4)],[0 2])

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-11-28
Ran in:
You mean
syms t t0 tf
A = [0 1 0 0;0 0 -4.91 0;0 0 0 1;0 0 73.55 0];
B = [0;2;0;-10];
[M,J] = jordan(A);
phi_s=M * expm(J*t) * inv(M);
Mat = (phi_s*B) * (phi_s*B).';
G_c = int(Mat,t0,tf)
G_c = 
size(G_c)
ans = 1×2
4 4
?
  3 个评论
显示 1更早的评论
Torsten
Torsten 2022-11-28
It's done automatically if you leave away the ";" after the line of the generated symbolic expression.
Paul
Paul 2022-11-29
Ran in:
Another option is shown in this Answer by @Sheng Chen
Torsten's solution by direct, symbolic integration
syms t t0 tf real
assumeAlso(tf >= t0);
A = sym([0 1 0 0;0 0 -4.91 0;0 0 0 1;0 0 73.55 0]);
B = sym([0;2;0;-10]);
[M,J] = jordan(A);
phi_s=M * expm(J*t) * inv(M);
Mat = (phi_s*B) * (phi_s*B).';
G_c(t0,tf) = int(Mat,t0,tf)
G_c(t0, tf) = 
Sheng's solution using symbolic expm
S(t0,tf) = G_c1(t0,tf,A,B);
Show they are equivalent
E = simplify(G_c - S)
E(t0, tf) = 
The latter can, in principle, also be used for direct numerical evaluation, (assuming no overflows, etc.), compare to vpa of the symbolic result
vpa(S(1,2.5))
ans = 
format long e
G_c1(1,2.5,double(A),double(B))
ans = 4×4
1.0e+00 * 3.705635009574892e+14 3.178000391865468e+15 -5.550904953993540e+15 -4.760528063354170e+16 3.178000391865468e+15 2.725494136524754e+16 -4.760527702652360e+16 -4.082690284319699e+17 -5.550904953993540e+15 -4.760527702652360e+16 8.315051463151344e+16 7.131095950415761e+17 -4.760528063354170e+16 -4.082690284319699e+17 7.131095950415761e+17 6.115720351147812e+18
function out = G_c1(t0,tf,A,B)
Qs = expm(A*t0)*B;
Q = Qs*Qs.';
temp = (tf-t0)*[-A Q; 0*A A.'];
temp = expm(temp);
if isa(temp,'sym')
out = simplify(expand(simplify(expand(temp(end/2+1:end,end/2+1:end).')) * simplify(expand(temp(1:end/2,end/2+1:end)))));
else
out = temp(end/2+1:end,end/2+1:end).' * temp(1:end/2,end/2+1:end);
end
end

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by