# Integral of controllability gramian

62 views (last 30 days)
Meena Farag on 31 May 2016
Commented: Meena Farag on 12 Oct 2021 at 19:36
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!

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.
Meena Farag on 12 Oct 2021 at 19:36
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!

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
May you help me with source or literature that provide information about how weak rank test by using controllability and observability matrix?

Roger Stafford on 1 Jun 2016
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.
##### 2 CommentsShowHide 1 older comment
Star Strider on 1 Jun 2016
Actually, I looked this up in my control reference (and Wikipedia Controllability Gramian). It’s matrix exponentiation, expm.

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
Star Strider on 1 Dec 2018
Post this as a new Question.
No one will respond to it here.