Choose Lyapunov function of a linear system

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

Answer 1

[David Goodmanson]

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;
then
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.

[Cristina Magnetti Gisolo]

Hi David,

thank you very much for your help, I'm going to study your analysis. My c and b parameters actually satisfy the condition you wrote above.

Now I'm trying with a Lyapunov function like this:

I saw that, in this way, the elements that come from the derivative of V(x) combine with one another, and I found a system of equalities and inequalities to impose in order to have negative.

Which tool/function shoud I use to solve a system of equations and inequalities?

Maybe symbolic toolbox and function solve?

[David Goodmanson]

Hi Christina, I modied the answer in light of what I found out today. It's interesting that you found a set of conditions for Vdot(x). Hopefully they will work.

Any V that you find will have to have its minimum at the the point N*[1 1 1 0 0 0]' now mentioned in the answer, and I don't think that the value of N is easy to calculate except as follows. In the expression

[z;y] = V*exp(Lambda*t)*g0  = V*exp(Lambda*t)*inv(V)*[z0;y0]

let t--> inf. Since exp(Lambda*t) is diagonal and five of the Lambdas have a negative real part, that means that the diagonal element corresponding to lambda = 0 will be 1, and every other element of that matrix is 0. Calling that matrix Q, then

[z;y] (t=inf) = V*Q*inv(V)*[z0;y0]

Easy to calculate numerically, not easy to find an algebraic solution.

Whether it's that calculation or one that you are doing on the Vdot calculation, it's possible, but I don't think likely, that symbolics will supply an answer. That's because you will probably will be getting into roots of polynomials of degree greater than 4, and there are no analytic solutions for those.

Rather than a Lyapunov function, maybe it can be shown that five of the eigenvalues have a negative real part, which would be the equivalent result, practically speaking.

[Cristina Magnetti Gisolo]

Hi, thank you very much, you're very kind.

Indeed I'm having some problems in the computation: I modified a bit the candidate function, which still has 12 unknowns:

I derived 9 equations that have degree equal to 2 and they allow to have strictly negative. I'm using Symbolic Math Toolbox and solve() to find a solution, but it seems to take too much time and I fear it won't be able to find it.

*Just got the result: unable to solve symbolically. I'll try with vpasolve

** Nope, I get that all the coefficients are zero and the message that solutions might be lost. Actually, if the coefficients are zero, the equations are solved, but it's of course a solution I don't like

[David Goodmanson]

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]

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]

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.

Answer 2

[Cristina Magnetti Gisolo]

Thank you very much. Actually, I have b and c fixed, so I see that 5 eigenvalues have negative real part.

[David Goodmanson]

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

Answer 3

[Sam Chak]

Hi @Cristina Magnetti Gisolo

Although my answer comes late, I hope it provides a closure on this matter. @David Goodmanson's analysis was helpful. The characteristic polynomial for

is given by

.

Since one λ can be factored out of every term in the polynomial, then you will get a zero eigenvalue, and it simply means that the matrix is singular. In other words, the system is marginally stable, or bounded when there are perturbations.

Now, if , , and in the matrix

then it can be shown that the system is exponentially stable.

d = 1.0001;
c = 1;
b = 1;
a = [0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1; -d*c c 0 -b 0 0; 0 -d*c c 0 -b 0; c 0 -d*c 0 0 -b]
rank(a)
eig(a)