53 views (last 30 days)

Show older comments

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?

Star Strider
on 1 Jun 2016

Edited: Star Strider
on 1 Jun 2016

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.

Star Strider
on 2 Jun 2016

As always, my pleasure!

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

Sheng Cheng
on 16 Feb 2017

Edited: Sheng Cheng
on 20 Feb 2017

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.

Ahmed Rashid
on 31 May 2016

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.

Bryan Jevon
on 30 Oct 2018

Roger Stafford
on 1 Jun 2016

Star Strider
on 1 Jun 2016

Rajani Metri
on 1 Dec 2018

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

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

Start Hunting!