A text version of the notebook
A Stability Analysis Method for Non-linear Systems
Abstract: Liapunov functions allow the determination of stability for non-linear systems without the need to find exact solutions. The Variable Gradient Method provides a systematic approach to determining a suitable Liapunov function. This package automates the steps involved in generating a Liapunov Function using the Variable Gradient Method. The notebook provides an overview of the theory behind non-linear stability analysis and a discussion of how to interact with the package, as well as some examples.
Theoretical Overview
The concept of stability has different definitions depending on whether the discussion is primarily linear time-invariant systems or non-linear systems.
Linear time-invariant systems are modeled as differential equations with constant coefficients. A determination of stability for these types of systems can be found through the use of Laplace Transform techniques.[3] If the poles of the system lie in the left-hand plane, then the time domain solution will fall to zero as time approaches infinity. If, on the other hand, the poles lie in the right hand plane, then the time domain solution will increase without bound as time approaches infinity. Stability, in the linear sense, is thought of in terms of the exact solution to the defining differential equation.
For non-linear differential equations this definition would also apply, if we could find exact solutions to the differential equations. Since in only very rare cases is this possible, the need to search for an alternate definition of stability is required. One approach is to look at the basic definitions in elementary physics of what constitutes stable and unstable equilibrium points. This is developed by utilizing energy concepts.
A quick example: We have a marble and a bowl. We place the bowl upside-down with the marble place at the center. The marble is in equilibrium, however, if it is moved - even slightly- it will fall off of the bowl. This is unstable equilibrium. Now flip the bowl over to it’s normal configuration and place the marble at the center. The marble is again in equilibrium. If the marble is moved slightly, it returns to the center position. This is stable equilibrium. Notice that we could increase the marbles speed sufficiently so that it escapes from the bowl.
For a simple one dimensional mechanical system, the total energy is the simple sum of the kinetic energy and the potential energy. The kinetic energy is the energy of motion, computed as the square of the speed, while potential energy is the energy due to a position in a force field. The total energy of the object must always be greater than or equal to the potential energy. (Note: the implications if it were not are that either the objects speed is an imaginary number or the mass is negative, neither one of which is naturally occurring.) For a stable system, the total energy of the object must be less than the maximum potential energy of the "bowl’s rim." This can be illustrated by taking some arbitrary potential energy function graphed against its position, as shown in figure 1. The lowest total energy possible occurs at the point x0, which translates to the potential energy U(x0). The total energy on the graph is represented as horizontal lines. Therefore, if the total energy of the object is less than E1, then it could not escape the confined energy well centered around x0. This object would then be thought of as being in stable equilibrium. If, on the other hand, the object had a total energy greater than E1, then it could escape the energy well around x0.
This concept of a one-dimensional energy well, can be expanded to two dimensions, as in the case of the bowl and the marble. It can also be extended to n-dimensional space, where the energy wells are surfaces.
If the total energy of the system is a function of time, then an additional requirement must be placed on it in order to insure that it remains in stable equilibrium. The constraint is that E(t) must be a decreasing function of time, or dE(t)/dt < 0.
System Models
Mathematical models of systems normally result in nth order differential equations. In the linear time-invariant case these equations all have constant coefficients. These constant coefficient differential equations can be solved by utilizing either classical or transform methods. For non-linear nth order differential equations there are no systematic solution methods. However, by converting the problem from an nth order differential equation to n 1st order differential equations, we can view possible solutions in a state space. For example, if we are interested in a second order system, then we could plot the system solutions in terms of the objects position vs. its speed. Take an ideal pendulum for example. The solution to its equation yields a family of circles in this space. Higher order systems are more difficult to visualize, but are conceptually similar.
The conversion process proceeds as follows [3]:
1) Let a system be described by the differential equation:
d^n y/dt ^n + a_1 d^(n-1) y/dt^(n-1) + ......+ a_n-1 dy/dt + a_n y = u(t)
where the an’s can be functions of both time and y or its derivatives.
2) Make the substitutions:
x[1]_dot = x[2] (= dx[1]/dt)
x[2]_dot = x[3]
.......
x[n-1]_dot = x[n]
x[n]_dot = -a_n x[1] -.....- a_1 x[n] + u(t)
We now have transformed the problem to n 1st order differential equations, where the n variables are known as state variables. These transformed equations can be thought of as existing in an n dimensional space and the equations written in matrix notation as:
(dX/dt =) X_dot = A X
Non-Linear Stability Criteria
If we can find an energy function (here after referred to as the V function) in the n dimensional state space that describes the system of interest, then by using the concepts of energy and equilibrium we can determine in what region of the state space is the system stable. This concept was put forth by A. M. Liapunov near the end of the nineteenth century in Russia.[1]
If we move an equilibrium point to the origin of the state space by linear transformations, then the following stability criteria apply:
1) If V(X) is positive for all values of the state variables within a region except at the origin, where V(0) = 0, then this scalar function is called positive definite in that region. [3]
2) If the time derivative of V(X) is negative for all values of the state variables within a region except the origin where V(0) = 0, then the scalar function is called negative definite in that region. [3]
If both criteria are true for a V function, then the system is considered to be asymptotically stable within that region. If 1 and 2 are true in the entire region of the state space, then the system is globally asymptotically stable. Simply put, this criteria is a restatement of the energy requirements of stability for simple mechanical systems. Criteria 2 indicates that as time increases the system must not increase its energy, since it could then drift out of the stable region. It should be noted that if V(X) and/or dV(X)/dt are not functions of all the state variables, then they are considered to be semi-definite.
So the problem of determining stability for non-linear systems is reduced to determining the energy function V, better known as a Liapunov Function. Finding a suitable V function was for the most part a matter of guessing until 1962. Then D.G. Schultz put forth the systematic method of using the gradient of V to determine a suitable Liapunov Function for a system.[1]
Variable Gradient Method
In order to meet the stability criteria set forth by Liapunov, a V function must be found and then it’s time derivative. Since the state variables are implicit functions of time, then: (partial is the partial derivative)
dV/dt = partial V/partial x[1] x[1]_dot + .......+ partial V/partial x[n] x[n]_dot
dV/dt = Grad[V] X_dot
where Grad[ ] is the Gradient Operator
V can also be found from the gradient of V by a path integration through state space. This path integration will be independent of the path if :
partial Grad[V][[i]]/partial x[j] = partial Grad[V] [[j]]/partial x[i]
which for the three dimensional case reduces to the vector identity:
Curl[Grad[V] ] = 0
The gradient of V has the form:
Grad[V] = alpha X
(where alpha is an n x n matrix of unknown coefficients)
Methodology
To generate the V function, using the Variable Gradient Method the following steps are performed. [1,4]
1) Form dV/dt as:
dV/dt = Grad[V] X_dot = alpha X X_dot = alpha X A X
2) Utilizing the vector identity and the constraint that dV/dt must be negative definite, obtain values for the alpha[i,j]’s. (At minimum, force dV/dt to be negative semi-definite).
This is a step which involves making quasi-arbitrary choices for the alpha matrix depending on the two constraints.
3) Form Gradient of V from the alpha matrix
4) Perform a path integral over all variables in the state space to determine the Liapounov function V.
5) Apply the appropriate theorems to determine stability.
While the method is straightforward, the process is long and arduous. This Mathematica implementation performs steps 1, allows the user to interact with the constraint equations until step 2 is completed, then performs steps 3 and 4.
Bibliography
1.Gibson, J.E. and Schultz, D.G. "The Variable Gradient Method of Generating Liapunov Functions with applications to Automatic Control Systems." Doctorial Dissertation, Purdue University, April 1962.
2. Halliday, D. and Resnick, R. "Fundamentals of Physics" Second Edition. Section 8.5 "Mechanical Energy and Potential Energy Curve." New York: Wiley & Sons 1981.
3. Ogata, K. "Modern Control Engineering" Section 6-5 and Chapter 15. Englewood Cliffs: Prentice-Hall, 1970.
4. Schultz, D.G. "The Generation of Liapunov Functions" Advances in Control Systems Vol.2, pages 1 - 63. New York: Academic Press, 1965.
Implementation
Functions
VariableGradient[A] determines the time derivative of the
Liapounov Function for an nth order system described by
the matrix A. Any functions in A must be specified with
the variables x[1],x[2],... The output is a numbered list
of constraints that should be forced to zero by direct
substitution. To aid in the process use the
FindConstraint function.
FindConstraint[number,variable] solves the numbered
constraint equation for the desired variable.
Curl determines the additional curl contraints on Grad[V] so
that the V function can be found by any path integration.
It computes D[GradV[i],x[j]] - D[GradV[j],x[i]] for all
possible combinations of i and j. It's output consists of
these constraints and the form of Grad[V].
FindV[ Path List ] integrates the Grad[V] function along the
desired straight-line path. If the path list is {1,2,3},
it integrates along x[1], with x[2] = x[3] = 0; then along
x[2], with x[1] = x[1] and x[3] = 0, etc. The output is
the resulting Liapounov Function V.
ResetVG eliminates all remaining variables, so that the
process can be restarted.
VariableGradient process flow:
a) Define Matrix A from system.
b) Invoke VariableGradient[ ] function.
1) the output of this function is:
i) list of constraints (cross terms of dV/dt) equation.
ii) explicit form of dV/dt
c) Set alpha[i,j]’s to eliminate constraint equations.
1) FindConstraint functions aids in determining alpha[i,j] values.
2) an update of the VariableGradient output can be displayed by invoking VDOT
d) Repeat step c until dV/dt is negative semi-definite or negative definite.
e) Invoke Curl to list Curl Constraints and Gradient of V.
f) Set alpha[i,j]’s to eliminate constraint equations.
1) an update of the Curl output can be displayed by invoking DELV.
g) Repeat step f until all DELV constraints are 0.
h) Invoke FindV[ ] function. This performs the path integration, the argument is a path list. (i.e.: If the path is x[2], x[1], x[3]; then the path list would be {2,1,3})
i) Invoke ResetVG to clear all variables for reuse.
Example
A system represented by a 2 x 2 Matrix
A system is represented by the Matrix Equation A1, defined as:
INPUT
A1 = {{-1,2 x[1]^2},{0,-1}}
INPUT
VariableGradient[A1]
OUTPUT
Constraints:
1 -alpha[1, 2] - alpha[2, 1]
V dot:
(-alpha[1, 2] - alpha[2, 1]) x[1] x[2] +
3 2
2 alpha[1, 1] x[1] x[2] - alpha[2, 2] x[2] +
2 2
x[1] (-alpha[1, 1] + 2 alpha[1, 2] x[2] )
INPUT
alpha[1,2]=alpha[2,1]=0;
INPUT
VDOT
OUTPUT
Constraints:
1 0
V dot:
2 3
-(alpha[1, 1] x[1] ) + 2 alpha[1, 1] x[1] x[2] -
2
alpha[2, 2] x[2]
alpha[1,1] = 1;
INPUT
Curl
OUTPUT
Constraints:
0
Grad V:
x[1]
alpha[2, 2] x[2]
INPUT
FindV[{1,2}]
OUTPUT
2 2
x[1] alpha[2, 2] x[2]
{----- + -----------------}
2 2
INPUT
ResetVG