Using matlab to solve Lyapunov function

Consider the following system

dx/dt=z
dz/dt=-x-z^3

Using the Lyapunov function as follows

V(x,z)=(x^2+z^2)

1. By applying LaSalle's invariance principle what can be said about the stability of the system?
2. What kind of mechanical system we can think of here?
3. If the output is y=z, construct an observer for the system.

[Steven Lord]

This sounds like a homework assignment. If that's correct, show us the code you've written so far to try to solve the problem and ask a specific question about where exactly you're having difficulty and you may receive some guidance.

Answer 1

[Sam Chak]

Despite this is a standard homework question, professors almost never use MATLAB to demo a rigorous Lyapunov or LaSalle stability proof for this kind of nonlinear system. However, this is how I would use MATLAB to assist me in the stability analysis.

To begin with, I'll just state that the only equilibrium point for this system is at the origin .

syms x z real

% Nonlinear system
f    = [z; -x - z^3]            % dx/dt = f(x)

f = 

% Lyapunov function candidate
V    = 1/2*x^2 + 1/2*z^2

V = 

% Derivative of V along f
Vdot = jacobian(V, [x, z])*f    % dV/dt = (dV/dx)*(dx/dt)

Vdot = 

Vdot = simplify(Vdot)

Vdot = 

% Test stability
assume(z ~= 0)              % for some unknown reasons, this assumption doesn't work
TF0 = isAlways(Vdot < 0)    % returns logical 1 (if true)

Warning: Unable to prove '-z^4 < 0'.

TF0 = logical
   0

assume(z < 0) 
TF1 = isAlways(Vdot < 0)    % returns logical 1 (if true)

TF1 = logical
   1

assume(z > 0) 
TF2 = isAlways(Vdot < 0)    % returns logical 1 (if true)

TF2 = logical
   1

is negative semi-definite because whenever even if . We look for the region where .

This set E is the entire x-axis (where ).

.

For a trajectory to stay inside E permanently, the derivative of with respect to time must also stay at zero ().

By substituting into the system dynamics:

,

For to equal 0 so the system remains trapped on the axis, x must also equal 0.

% Apply LaSalle's invariance principle
assume(z, 'clear')
sol = solve([f(1) == 0, subs(f(2), z, 0) == 0], [x, z])

sol = struct with fields:
    x: 0
    z: 0

The only trajectory that can stay identically inside E forever is the origin itself, .

By LaSalle's Invariance Principle, every trajectory approaches M as . Because our candidate function is also radially unbounded (meaning ) as ), this proof holds globally and implies that the equilibrium point at the origin is globally asymptotically stable as well.

However, despite the Lyapunov math says that the system converges to exactly zero as time approaches infinity, in practice, because of that weak cubic damping , the system converges so very slowly that the ode45 solver would have to run for an immense amount of time to get anywhere near absolute zero.

% Nonlinear system
function dx = system(t, x)
    dx      = zeros(2, 1);
    dx(1)   = x(2);
    dx(2)   = - x(1) - x(2)^3;
end

[t, x]  = ode45(@system, [0 3500], [1; 0]);
plot(t, x(:,1)), grid on
xlabel('Time')
ylabel('Amplitude')

Warning: Graphics acceleration hardware is unavailable. Graphics quality and performance might be diminished. See MATLAB System Requirements.

[Steven Lord]

The reason for your test stability comment:

syms z real
Vdot = -z^4;
% Test stability
assume(z ~= 0)              % for some unknown reasons, this assumption doesn't work
TF0 = isAlways(Vdot < 0)    % returns logical 1 (if true)

Warning: Unable to prove '-z^4 < 0'.

TF0 = logical
   0

is because your assume call overwrites the assumption that z is real.

syms x real
assumptions(x)

ans = 

assume(x ~= 0)
assumptions(x)

ans = 

To keep the assumption of reality and add another, use assumeAlso.

syms y real
assumptions(y)

ans = 

assumeAlso(y ~= 0)
assumptions(y)

ans = 

If you do this, now MATLAB can prove that -y^4 is always strictly less than 0.

V = -y^4;
isAlways(V < 0)

ans = logical
   1

[Sam Chak]

Hi @Steven Lord, thanks for showing how to set assumptions correctly.