Lyapunov

EXAMPLE: Energy and stability of pendulum system

Figure 5 Pendulum system with parameters

$J\ddot{\theta }\left(t\right)=-\mu \dot{\theta }\left(t\right)-Mgl\sin \left(\theta \left(t\right)\right) \left(2\right)$

Let $\mathit{x}\left(t\right)={\left[\begin{array}{cc}{x}_1\left(t\right)& x_2\left(t\right)\end{array}\right]}^T={\left[\begin{array}{cc}\theta \left(t\right)& \dot{\theta }\left(t\right)\end{array}\right]}^T$ , denote the state variables. Then equation (2) is rewritten as

$\left[\begin{array}{c}{\dot{x}}_1\left(t\right)\\ {\dot{x}}_2\left(t\right)\end{array}\right]=\left[\begin{array}{c}{x}_2\left(t\right)\\ -\frac{Mgl}{J}\sin \left(x_1\left(t\right)\right)-\frac{\mu }{J}{x}_2\left(t\right)\end{array}\right] \left(3\right)$

Equilibrium point is obtained as

$\left[\begin{array}{c}{\dot{x}}_1\left(t\right)\\ {\dot{x}}_2\left(t\right)\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]\Rightarrow \left\{\begin{array}{l}0=x_{2e}\\ 0=-\frac{Mgl}{J}\sin \left(x_{1e}\right)-\frac{\mu }{J}{x}_{2e}\end{array}\Rightarrow \left\{\begin{array}{l}{x}_{2e}=0\\ \sin \left(x_{1e}\right)=0\end{array}\right.\right. \left(4\right)$

Therefore, ${\mathit{x}}_e=0$ is one of the Equilibrium point . Besides, the mechanical energy of pendulum system equation (3) is

$\varphi \left(x\right(t\left)\right)=\frac{1}{2}J{x}_2^2\left(t\right)+Mgl(1-\cos \left(x_1\left(t\right)\right)) \left(5\right)$

where the first term at right hand is called kinetic energy , the second one is called potential energy as is shown in Figure 6

Figure 6 Mechanical energy of pendulum

When $\mathit{x}\left(t\right)=0$ , that the position of pendulum is six o'clock, and velocity of it is 0, then $\varphi \left(\mathit{x}\right(t\left)\right)=0$ . For any other state $\mathit{x}\left(t\right)$ there is $\varphi \left(\mathit{x}\right(t\left)\right)>0$ . The derivative of $\varphi \left(\mathit{x}\right(t\left)\right)$ with equation (3) as simultaneous equations in order to cancel derivatives of states, the following formula holds as

$\begin{array}{rcl}\dot{\varphi }\left(\mathit{x}\right(t\left)\right)& =& J{x}_2\left(t\right){\dot{x}}_2\left(t\right)+Mgl{\dot{x}}_1\left(t\right)\sin \left(x_1\left(t\right)\right)\\ & =& \left(J{\dot{x}}_2\right(t)+Mgl\sin \left(x_1\left(t\right)\right)x_2(t)=-\mu x_{{}_2}^2(t)\end{array} \left(6\right)$

Assume

Matlab code of machanical energy is:

    [X,Y] = meshgrid(-7:0.1:7,-10:0.1:10);
    Z = Y.^2 + 9.8*(1-cos(X));
    surf(X,Y,Z)
    xlabel('angle')
    ylabel('angular velocity')

Figure 7 Mechanical energy of pendulum (Matlab chart)

Matlab code of $\dot{\varphi }\left(\mathit{x}\right(t\left)\right)$ is:

    [X,Y] = meshgrid(-7:0.1:7,-10:0.1:10);
    Z = -Y.^2;
    surf(X,Y,Z)
    xlabel('angle')
    ylabel('angular velocity')

Figure 8 Derivative of pendulum mechanical energy (Matlab chart)

Note that $\varphi \left(\mathit{x}\right(t\left)\right)$ is a positive-definite function (PDF), ( $\varphi \left(\mathit{x}\right(t\left)\right)>0, \forall \mathit{x}\left(t\right)\ne 0\in {\mathrm{ℝ}}^2$ ), while $\dot{\varphi }\left(\mathit{x}\right(t\left)\right)\le 0$ . Unlike single-variable functions, whose time-derivatives indicate monotonicity, multi-variable function as shown in equation (5) , whose time-derivative is utilized to exploit the stability of an equilibrium point. ^[ref] Generally,

$\frac{dV}{dt}=\frac{\partial V}{\partial x_1}\frac{d{x}_1}{dt}+\frac{\partial V}{\partial x_2}\frac{d{x}_2}{dt}+\cdots +\frac{\partial V}{\partial x_n}\frac{d{x}_n}{dt} \left(7\right)$

which in a scalar (dot) product can be written as

$\frac{dV}{dt}=\left(grad V, \frac{d\mathit{X}}{dt}\right), grad V=\left(\frac{\partial V}{\partial x_1}, \frac{\partial V}{\partial x_2}, \dots , \frac{\partial V}{\partial x_n}\right) , \frac{d\mathit{X}}{dt}=\left(\frac{d{x}_1}{dt}, \frac{d{x}_2}{dt}, \dots , \frac{d{x}_n}{dt}\right) \left(8\right)$

Here, the first vector is the gradient of $V\left(\mathit{X}\right)$ , id est, it’s always directed toward the greatest increase in $V\left(\mathit{X}\right)$ . The second vector is velocity vector, it is always tagent to the phase trajectory .

(a)

(b)

Figure 9 (a) Lyapunov function. (b) Phase trajectory.

Take pendulum system equation (3) again for example, the gradient and phase trajectory of which is obtained by Matlab code as below ^[ref] :

    x=-pi:pi/20:pi;
    y=x';
    z = y.^2 + 9.8*(1-cos(x));
    [px,py] = gradient(z);
    figure
    contour(x,y,z)
    hold on
    quiver(x,y,px,py)
    hold off
    xlabel('angle')
    ylabel('angular velocity')

    x=-pi:pi/20:pi;
    y=x';
    z = y.^2 + 9.8*(1-cos(x));
    vx = y*ones(1,41);
    vy=-9.8*sin(x)-y*ones(1,41);
    figure
    contour(x,y,z)
    hold on
    quiver(x,y,vx,vy)
    hold off
    xlabel('angle')
    ylabel('angular velocity')

(a)

(b)

For phase trajectory and state trajectory of system (2) ^{[ref] [ref] [ref] [ref]}

    syms x(t)
    [V] = odeToVectorField(diff(x, 2) == -diff(x,1) - 9.8*sin(x));
    M = matlabFunction(V,'vars', {'t','Y'});
    sol = ode45(M,[0 20],[2 2]);
    z = linspace(0,20,1000);
    x1 = deval(sol,z,1);
    x2 = deval(sol,z,2);
    figure
    hold on;
    plot(x1,x2);
    xlabel('angle')
    ylabel('angular velocity')
    hold off;
    figure
    hold on;
    fplot(@(x)deval(sol,x,1), [0, 20])
    xlabel('time')
    ylabel('angle')
    hold off;
    figure
    hold on;
    fplot(@(x)deval(sol,x,2), [0, 20])
    xlabel('time')
    ylabel('angular velocity')
    hold off;

Figure 11 Phase portrait

Figure 12(a) angle trajectory

Figure 12(b) angular velocity trajectory

Test initial condition for [5,5],[3.14,5],[3.14,0],[3.14,-5],[pi,0].
Back to equation (6) , when $x_1\left(t\right)\ne 0$ , if $x_2\left(t\right)=0$ , there are chances that $\mathit{x}\left(t\right)=0$ is stable. Whatsoever,

$\left\{\begin{array}{l}{x}_2\left(T\right)={\dot{x}}_1\left(T\right)=0\\ {\dot{x}}_2\left(T\right)=-\frac{Mgl}{J}\sin \left(x_1\left(T\right)\right)\end{array}\right. \left(9\right)$

When $x_1\left(T\right)\ne 0$ , then $\sin \left(x_1\left(T\right)\right)\ne 0\Rightarrow {\dot{x}}_2\left(T\right)\ne 0$ . Such that, $\dot{\varphi }\left(\mathit{x}\right(t\left)\right)=-\mu x_2^2\left(t\right)\ne 0$ , when $t>T$ . Therefore, such "equilibriums" are not stable.
Figure (13) and figure (14) represent mechnical energy of pendulum system (2) and its time derivatives respectively. The matlab code is

    ...
    fplot(@(x) 0.5*deval(sol,x,2).^2+9.8-9.8*cos(deval(sol,x,1)), [0, 20])
    fplot(@(x) -deval(sol,x,2).^2, [0, 20])

Figure 13 $\theta \left(\mathit{x}\right(t\left)\right), \mathit{x}\left(0\right)=[2,5]$

Figure 14 $\dot{\theta }\left(\mathit{x}\right(t\left)\right), \mathit{x}\left(0\right)=[2,5]$