MATLAB Answers

Choose Lyapunov function of a linear system

27 views (last 30 days)
Cristina Magnetti Gisolo
Cristina Magnetti Gisolo on 27 Nov 2020
Commented: David Goodmanson on 30 Nov 2020
Hello everyone,
I would like to perform the Lyapunov stability of the following linear system. It is the linearization of a quite complex nonlinear system around the equilibrium in the origin. The Jacobian matrix (A below) has a null eigenvalue, so I can't say anything about the stability of the equilibrium point.
The equilibrium point is .
In matrix form it is like: , where and
A = 0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
-c c 0 -b 0 0
0 -c c 0 -b 0
c 0 -c 0 0 -b
C and b are positive numbers.
I would like to perform the Lyapunov stability, but I'm not sure about the kind of Lyapunov function to select. I tried to employ quadratic functions where all the state variables are squared and summed together, everyone with its own parameter, but it doesn't seem a good idea. I tried also to have other quadratic functions combining the state variables in a more complex way, but it didn't work.
Does anyone have any idea? Just about the kind of function to select. I keep on trying...
Thank you very much in advance


Sign in to comment.

Answers (2)

David Goodmanson
David Goodmanson on 29 Nov 2020
Edited: David Goodmanson on 29 Nov 2020
Hi Cristina,
The first part of this is probably telling you what you already know, but denoting
C = c*[1 -1 0; 0 1 -1; -1 0 1]
you have the block matrix equation
d/dt [z] = A * [z] = [ 0 I] * [z]
[y] [y] [-C -bI] [y]
where A is made up of a 3x3 blocks as you have shown. Let V and Lambda be the eigenvectors and eigenvalues of A where the columns of V are eigenvectors, and Lambda is diagonal. Suppose the Initial starting point is [z0;y0]. [z0;y0] is a linear combination of eigenvectors
[z0;y0] = V*g0
where the 6x1 vector of coefficients vector g0 is determined by
g0 = inv(V)*[z0;y0]
< this is for illustration, in actual code you would use the better method g0 = V \ [z0;y0] >
For the time domain solution each individual eigenvector is multiplied by exp(lambda*t) using its corresponding lambda. The time domain solution is
[z;y] = V*exp(Lambda*t)*g0 = V*exp(Lambda*t)*inv(V)*[z0;y0]
Any eigenvector with eigenvalue = 0 will not change with time.
Let col1 = [1 1 1]' and col0 = [0 0 0]'. Note that
C*col1 = col0 = 0
which is the key. It follows that
d/dt [col1] = [ 0 I] * [col1] = [col0] = 0
[col0] [-C -bI] [col0] [col0]
so [col1; col0] = [1 1 1 0 0 0]' is the (unnormalized) eigenvector with eigenvalue 0 and is constant in time.
There is also another (unnormalized) eigenvector, [col1; -b*col1] = [1 1 1 -b -b -b]' with eigenvalue -b. It gains a factor of exp(-b*t) and takes a straight-line drop down to the origin.
That leaves a 4-d space perpendicular to those eigenvectors, to determine what is going on. Take a look at
b = 4;
c = 1;
C = c*[ 1 -1 0
0 1 -1
-1 0 1];
A = [zeros(3,3) eye(3,3)
-C -b*eye(3,3)];
lambda = eig(A)
lambda =
0.0000 + 0.0000i
-0.3960 + 0.2700i
-0.3960 - 0.2700i
-3.6040 + 0.2700i
-3.6040 - 0.2700i
-4.0000 + 0.0000i
So far so good. The eigenvalues 0 and -4 are there, and each of the eigenvalue has a negative real part. In the cases with complex conjugate eigenvalues the coordinates that span that 2d subspace spiral down to the origin. All together five eigenvectors go to zero, leaving just the part of [z;x] that involves [1 1 1 0 0 0]'. So the situation is basically stable (but see below).
However, if you try
b = 4;
c = 40;
lambda =
-4.2190 + 7.8054i
-4.2190 - 7.8054i
0.0000 + 0.0000i
0.2190 + 7.8054i
0.2190 - 7.8054i
-4.0000 + 0.0000i
Now two eigenvectors have positive real part and their coordinates spiral out, so that case is unstable.
You get basically stable solutions for c small enough, and it turns out that c can actually be fairly large. The condition turns out to be
0 < c < 2*b^2
Supposing that c is in that interval, then the trajectory will ultimately end up at rest somewhere on the line [z;y] = m*[1 1 1 0 0 0]' where m can be any value. I suppose you could call it a line of fixed points. But just where you end up on that line depends on the initial conditions [z0;y0] , and on b. Interestingly, N does not depend on the value of c (as long as you are in the allowed interval). Let's say the end point for the trajectory is N*[1 1 1 0 0 0]'. If the initial velocites y0 are large, then N is large and can increase without bound if the velocities do. In that sense the system is not exactly stable although it always settles out at to one particular point.
For c in the right range, it should be possible to come up with a Lyapunov function (taking into account that the bottom of the well occurs at N*[1 1 1 0 0 0]'), although that might not be an easy task.


Show 3 older comments
David Goodmanson
David Goodmanson on 30 Nov 2020
Hi Cristina, would it be enough to show that in the linear model, five of the eigenvalues have negative real part, which shows that the solution converges on a single point? Or do you want to find a Lyapunov function in the hope that it can be shown to still work in the full nonlinear problem?
Cristina Magnetti Gisolo
Cristina Magnetti Gisolo on 30 Nov 2020
I was looking to prove the stability of the linear system and, since one of the eigenvalues is null, I chose to study the Lyapunov stability.
Anyway, I don't get any solution, I don't know if it is not able to find it or if there's no solution. The result gives empty variables. I'll give up...
Thank you very much for your help and your time.
David Goodmanson
David Goodmanson on 30 Nov 2020
Hi Christina, before you give up, if you are still interested, I want to mention that it is possible to prove algebraically that with 0 < c < 2*b^2, five of the eigenvalues have negative real part. I don't know how to do that with symbolics, but with the archaic method of pencil and paper, yes.

Sign in to comment.

Cristina Magnetti Gisolo
Cristina Magnetti Gisolo on 30 Nov 2020
Thank you very much. Actually, I have b and c fixed, so I see that 5 eigenvalues have negative real part.

  1 Comment

David Goodmanson
David Goodmanson on 30 Nov 2020
good sense to proceed on. Although the algebraic proof works, it is a bit of a slog and I would say not very enlightening.

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