Integral of controllability gramian

I am having trouble finding a way to integrate the controllability gramian in Matlab. My system is unstable so I can't use the built in function. My system also has eigenvalues on the imaginary axis so I can't use a function I found online for unstable systems. So i'm trying to integrate the controllability gramian for some finite time interval on MATLAB, but it seems impossible. Here's what the integral looks like.

Where A is a 6x6 matrix and B is a 6x1 matrix. The matrix exponential in the equation is what's causing me the most trouble.

Any ideas?

Answer 1

[Star Strider]

Try this with your matrix and vector:

A = rand(6);                                        % Create Data
B = rand(6, 1);                                     % Create Data
f = @(tau) expm(A*tau)*B*B'*expm(A'*tau);           % Integrand
W = @(t) integral(f, 0, t, 'ArrayValued',1);        % Controllability Gramian
Wt = W(1)

The integral function was introduced in R2012a. Before that, I believe the appropriate function is quadv.

[Meena Farag]

This worked! Thank you so much. One questions though, I'm using quadv since I have MATLAB 2011. When I input the last two arguments

'ArrayValued',1);

I get another three lines under my gramian... It looks like this:

W =

      0.0071   -0.0083    0.0572    0.0081   -0.0120   -0.0278
     -0.0083    0.0097   -0.0642   -0.0094    0.0140    0.0344
      0.0572   -0.0642    0.5600    0.0679   -0.0990   -0.1558
      0.0081   -0.0094    0.0679    0.0093   -0.0138   -0.0300
     -0.0120    0.0140   -0.0990   -0.0138    0.0205    0.0462
     -0.0278    0.0344   -0.1558   -0.0300    0.0462    0.1647

         9     1.0000000000    2.71580000e-001     0.0001603555
        11     1.2715800000    4.56840000e-001     0.0007011706
        13     1.7284200000    2.71580000e-001     0.0009244803

What are those last three lines?

[Star Strider]

My pleasure!

The quadv function does not support the 'ArrayValued',1 name-value pair. Instead, it would interpret those as the tolerance and trace, with unpredictable results. It is best not to use them.

Instead, just use:

W = @(t) quadv(f, 0, t);                            % Controllability Gramian
Wt = W(1)

[Meena Farag]

Awesome. This simple problem has been stressing me out for a while, lol. Thanks again!

[Star Strider]

As always, my pleasure!

I was surprised that your question hasn’t been asked before.

Answer 2

[Sheng Cheng]

I have a different way for computing the controllability gramian matrix based on an early paper, 'Computing integrals involving the matrix exponential', by Charles Van Loan. The paper can be found here: https://www.cs.cornell.edu/cv/ResearchPDF/computing.integrals.involving.Matrix.Exp.pdf

I will skip the mathematical rigorous proof in the paper. In fact, all the result you need to compute the gramian matrix is written in the left column on the first page. Especially, equation (1.2) is the form we are looking for. (Please read the paper for the extremely simple structure of this integral (actually the controllability gramian is indeed an integral involving matrix exponential)).

The code is just in two lines

    A = rand(6);                                        % Create Data
    B = rand(6, 1);                                     % Create Data
    temp = expm([-A B*B';zeros(6,6) A']);               % Coming from the first equation below (1.4)
    Wc = temp(7:12,7:12)'*temp(1:6,7:12);               % Coming from the second equation below (1.4)

Here, you don't need to define a function and an integral like the one suggested by Star Strider. All you need is expm and then some very simple matrix operation.

[Star Strider]

Interesting. Noted.

Thank you.

[Sheng Cheng]

You are welcome!

[Star Strider]

In all my control courses, I never encountered that. The paper you attached is in my collection.

+1 Vote!

[Meena Farag]

It's been five years since I originally posted this question while in grad school. I can't believe i'm only coming across this alternative method now. What a great paper. Thanks for sharing!

[Martín]

Just read this and I had to comment this is a beatuful way to solve the problem.

Answer 3

[Ahmed Rashid]

Why don't you check the rank of the controllability matrix?

C = rank([B AB A^2B ... A^(n-1)*B])

If C has full rank, then the system is controllable.

[Meena Farag]

The rank test is generally a weak test for controllability. The gramian provides a degree of controllability. My system is in fact completely controllable, but I need to know how controllable it is, and the controllability gramian will tell me that.

[Ahmed Rashid]

Calculating the integral is equivalent to calculating the Lyapunov equation.

A*W+W*A^T+B*B^T=0

So you can find it by solving the above system of equations in Matlab.

[Meena Farag]

It is only equivalent for stable systems, my system is unstable. I found that I have no choice but to compute the integral for some finite time.

[Bryan Jevon]

May you help me with source or literature that provide information about how weak rank test by using controllability and observability matrix?

Answer 4

[Roger Stafford]

In matlab there is a very important difference between e.^((A.’)*τ) and e^((A.’)*τ) (without the dot.) The first of these is an element-wise exponentiation and the second a matrix exponentiation. If you use exp((A.’)*τ), it will produce the element-wise version. I rather suspect you want the element-wise version.

[Meena Farag]

Yes, I am using the element-wise version (and expm, not exp). I am still not able to integrate. I still get the same error:

Error using  .* 
Matrix dimensions must agree.

Error in @(t)(expm(A.*t)*B1)*(expm(A.*t)*B1)'

Error in quad (line 76)
y = f(x, varargin{:});

Error in DOC (line 63)
quad(fn(:,1),0,1)

[Star Strider]

Actually, I looked this up in my control reference (and Wikipedia Controllability Gramian). It’s matrix exponentiation, expm.

Answer 5

[Rajani Metri]

How to calculate Minimum control u*(t) required to state transfer from x1(t) to x2(t) and from it the states x1*(t) and x2(*)? also how to Plot them?

Thank you

[Star Strider]

Post this as a new Question.

No one will respond to it here.