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

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])

Answer 1

[Torsten]

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

?

[mpz]

@Torsten thanks. How did you have matlab to spill out the the output in the format for alpha1, alpha2,...I'm trying to do that instead of using LATEX.

[Torsten]

It's done automatically if you leave away the ";" after the line of the generated symbolic expression.

[Paul]

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