Lyapunov stability

Given a linear dynamic system x' = Ax, our goal is to prove stability by finding a symmetric matrix P satisfying

Define a matrix A and the symmetric matrix P.

A = [-1 2;-3 -4];
P = sdpvar(2,2);

Having P, we are ready to define the constraints.

F = set(P > 0) + set(A'*P+P*A < 0);

To avoid the zero solution or an unbounded solution, we constrain the trace of the matrix (Of course, this is not the only way. We could have used, e.g., the constraint P>I instead)

F = F + set(trace(P) == 1);

At this point, we are ready to solve our problem. But first, we display the collection of constraints to see what we have defined.

F
+++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|      Constraint|                      Type|
+++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Numeric value|     Matrix inequality 2x2|
|   #2|   Numeric value|     Matrix inequality 2x2|
|   #3|   Numeric value|   Equality constraint 1x1|
+++++++++++++++++++++++++++++++++++++++++++++++++++

We only need a feasible solution, so one argument is sufficient when we call solvesdp to solve the problem.

solvesdp(F);
P_feasible = double(P);

The resulting constraint satisfaction can be conveniently investigated with checkset.

checkset(F)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|      Constraint|                  Type|   Primal residual|   Dual residual|   Compl. slack|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Numeric value|                   LMI|           0.32029|     2.7929e-014|     4.983e-014|
|   #2|   Numeric value|                   LMI|            1.6706|     2.5781e-015|     1.936e-014|
|   #3|   Numeric value|   Equality constraint|      -2.0151e-014|    -1.0991e-014|    2.2148e-028|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Minimizing, e.g., the top-left element of P is done by specifying an objective function.

F = set(P > 0) + set(A'*P+P*A < 0);
solvesdp(F,P(1,1));