27 views (last 30 days)

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

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;

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.

David Goodmanson
on 30 Nov 2020

David Goodmanson
on 30 Nov 2020

Cristina Magnetti Gisolo
on 30 Nov 2020

David Goodmanson
on 30 Nov 2020

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

Start Hunting!
## 0 Comments

Sign in to comment.