Solving an integral when the parameters is a matrix - Controllability Gramian
2 views (last 30 days)
Show older comments
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])
0 Comments
Accepted Answer
Torsten
on 28 Nov 2022
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)
size(G_c)
?
3 Comments
Torsten
on 28 Nov 2022
It's done automatically if you leave away the ";" after the line of the generated symbolic expression.
Paul
on 29 Nov 2022
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)
Sheng's solution using symbolic expm
S(t0,tf) = G_c1(t0,tf,A,B);
Show they are equivalent
E = simplify(G_c - S)
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))
format long e
G_c1(1,2.5,double(A),double(B))
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
More Answers (0)
See Also
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!