MATLAB Answers

Integral of controllability gramian

53 views (last 30 days)
Meena Farag
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?

Accepted Answer

Star Strider
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.
  4 Comments
Star Strider
Star Strider on 2 Jun 2016
As always, my pleasure!
I was surprised that your question hasn’t been asked before.

Sign in to comment.

More Answers (4)

Sheng Cheng
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.
  4 Comments
Meena Farag
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!

Sign in to comment.


Ahmed Rashid
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.
  4 Comments
Bryan Jevon
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?

Sign in to comment.


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

Sign in to comment.


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

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by