To check the stability of an autonomous linear time invariant system, LMI can be used to solve the Lyapunov equation.
\[ \min tr(X) \\ \text{subject to} \quad X=X^T>0 \\ A^T X+XA+Q<0 \]using SCS, JuMP, LinearAlgebra, ControlSystems
A=[
0.0 1.0 0.0
0.0 0.0 1.0
-1.0 -2.0 -3.0
]
Q=Matrix(1.0I,3,3)
n1,=size(A);
model = Model(SCS.Optimizer)
@variable(model,X[1:n1,1:n1],PSD)
@objective(model, Min, tr(X))
@constraint(model, -(A'*X+X*A+Q) in PSDCone())
optimize!(model)
@show P=value.(X)
@show eigvals(A'*P+P*A)
@show lyap(A',Q)------------------------------------------------------------------
SCS v3.2.4 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 6, constraints m: 12
cones: s: psd vars: 12, ssize: 2
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 100000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 10
lin-sys: sparse-direct-amd-qdldl
nnz(A): 21, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.83e+01 1.01e+00 3.82e+01 -1.85e+01 1.00e-01 1.03e-04
175| 3.86e-04 5.08e-06 8.11e-06 7.50e+00 1.00e-01 4.04e-04
------------------------------------------------------------------
status: solved
timings: total: 4.05e-04s = setup: 6.51e-05s + solve: 3.40e-04s
lin-sys: 3.46e-05s, cones: 2.63e-04s, accel: 9.32e-06s
------------------------------------------------------------------
objective = 7.499986
------------------------------------------------------------------
P = value.(X) = [2.300508454499796 2.100191834840075 0.5001406354987649; 2.100191834840075 4.59952292379182 1.2999994479883699; 0.5001406354987649 1.2999994479883699 0.5999588292756685]
eigvals(A' * P + P * A) = [-1.0003734956415766, -0.9999457528420768, -0.9993302244644763]
lyap(A', Q) = [2.2999999999999985 2.0999999999999988 0.4999999999999999; 2.0999999999999983 4.599999999999994 1.299999999999999; 0.5 1.2999999999999994 0.6000000000000001]
To find the stabilizing control,
\[\begin{align*} \max \quad &tr(X^{-1}) \\ \text{subject to} \quad &X^{-1}=\left( X^{-1} \right)^T>0 \\ &\begin{bmatrix} X^{-1}A^T+L^T B^T + AX^{-1} +BL & X^{-1}\\ X^{-1} & -Q^{-1} \end{bmatrix} <0 \end{align*} \]using JuMP, LinearAlgebra, ControlSystems
# using SCS
# using ProxSDP
using Clarabel
A=[
0.0 1.0 0.0
0.0 0.0 1.0
1.0 -2.0 -3.0
]
B=[
0.0
0.0
1.0
]
n1,=size(A);
Q=Matrix(1.0I,n1,n1)
model = Model(Clarabel.Optimizer)
@variable(model,Xi[1:n1,1:n1],PSD)
@variable(model,L[1,1:n1])
# @objective(model, Max, 0)
@objective(model, Max, tr(Xi))
# @constraint(model, Xi.==Xi')
# @constraint(model, Xi in PSDCone())
@constraint(model, -[Xi*A'+L'*B'+A*Xi+B*L Xi;Xi -inv(Q)] in PSDCone())
# @constraint(model, -(Xi*A'+L'*B'+A*Xi+B*L+Q) in PSDCone())
# @constraint(model, -[Xi*A'+L'*B'+A*Xi+B*L Xi*Q;Xi*Q -I] in PSDCone())
# JuMP.optimize!(model)
optimize!(model)
Pi=value.(Xi)
LL=Array(value.(L))
K=LL*inv(Pi)-------------------------------------------------------------
Clarabel.jl v0.7.1 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 9
constraints = 27
nnz(P) = 0
nnz(A) = 33
cones (total) = 2
: PSDTriangle = 2, numel = (6,21)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-05, max_scale = 1.0e+05
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 0.0000e+00 -4.3264e+00 4.33e+00 5.29e-01 2.77e+00 1.00e+00 1.69e+00 ------
1 -1.2847e+00 -2.0704e+00 6.12e-01 6.99e-02 3.81e-01 2.78e-01 3.68e-01 8.01e-01
2 -6.1478e+00 -7.3642e+00 1.98e-01 2.35e-02 1.67e-01 1.58e+00 2.07e-01 7.19e-01
3 -4.4963e+00 -4.6823e+00 4.14e-02 7.05e-03 5.14e-02 3.78e-01 6.22e-02 8.96e-01
4 -7.3643e+00 -7.1921e+00 2.40e-02 1.24e-03 1.13e-02 3.87e-01 1.16e-02 8.42e-01
5 -1.7585e+01 -1.5980e+01 1.00e-01 4.65e-04 6.76e-03 2.00e+00 4.21e-03 8.96e-01
6 -1.5895e+01 -1.5303e+01 3.87e-02 2.21e-04 3.35e-03 7.70e-01 2.02e-03 6.87e-01
7 -2.5345e+01 -2.4554e+01 3.22e-02 3.06e-05 5.28e-04 8.48e-01 2.83e-04 8.69e-01
8 -5.8900e+01 -5.4477e+01 8.12e-02 9.06e-06 1.73e-04 4.51e+00 8.64e-05 9.41e-01
9 -6.7226e+01 -6.5241e+01 3.04e-02 1.54e-06 2.96e-05 2.01e+00 1.49e-05 9.26e-01
10 -1.5145e+02 -1.4327e+02 5.71e-02 2.47e-07 4.78e-06 8.20e+00 2.44e-06 9.87e-01
11 -2.2087e+02 -2.1214e+02 4.12e-02 4.69e-08 9.04e-07 8.74e+00 4.65e-07 8.34e-01
12 -5.0810e+02 -4.7048e+02 8.00e-02 1.15e-08 2.22e-07 3.76e+01 1.15e-07 9.72e-01
13 -5.2979e+02 -5.1121e+02 3.63e-02 4.09e-09 7.86e-08 1.86e+01 4.09e-08 7.48e-01
14 -7.8134e+02 -7.5546e+02 3.43e-02 9.47e-10 1.82e-08 2.59e+01 9.50e-09 8.30e-01
15 -1.9001e+03 -1.7814e+03 6.66e-02 1.62e-10 3.11e-09 1.19e+02 1.63e-09 9.90e-01
16 -2.3204e+03 -2.2428e+03 3.46e-02 3.99e-11 7.65e-10 7.75e+01 4.02e-10 7.93e-01
17 -4.3760e+03 -4.1447e+03 5.58e-02 9.63e-12 1.85e-10 2.31e+02 9.70e-11 9.53e-01
18 -6.5122e+03 -6.2461e+03 4.26e-02 1.87e-12 3.58e-11 2.66e+02 1.88e-11 8.35e-01
19 -1.3745e+04 -1.2815e+04 7.26e-02 4.63e-13 8.91e-12 9.30e+02 4.68e-12 9.54e-01
20 -1.4292e+04 -1.3696e+04 4.36e-02 2.30e-13 4.43e-12 5.97e+02 2.33e-12 6.05e-01
21 -2.0360e+04 -1.9646e+04 3.63e-02 5.40e-14 1.05e-12 7.14e+02 5.53e-13 8.15e-01
22 -4.7166e+04 -4.4148e+04 6.84e-02 1.52e-13 3.56e-13 3.02e+03 1.09e-13 9.85e-01
23 -5.4313e+04 -5.2417e+04 3.62e-02 4.84e-14 1.13e-13 1.90e+03 3.48e-14 7.42e-01
24 -8.0705e+04 -7.7713e+04 3.85e-02 6.93e-14 3.99e-14 2.99e+03 9.28e-15 8.34e-01
25 -1.6546e+05 -1.5671e+05 5.58e-02 6.74e-14 7.40e-14 8.75e+03 1.47e-15 9.54e-01
26 -2.4752e+05 -2.3538e+05 5.16e-02 3.07e-14 1.26e-14 1.21e+04 3.72e-16 8.23e-01
27 -4.8859e+05 -4.5607e+05 7.13e-02 2.60e-14 7.59e-15 3.25e+04 7.96e-17 9.53e-01
28 -5.1190e+05 -4.8646e+05 5.23e-02 2.19e-14 1.24e-14 2.54e+04 4.98e-17 4.50e-01
29 -7.1854e+05 -6.8854e+05 4.36e-02 1.26e-14 7.63e-15 3.00e+04 1.34e-17 8.01e-01
30 -1.4800e+06 -1.3833e+06 7.00e-02 7.51e-15 1.12e-13 9.68e+04 3.28e-18 9.57e-01
31 -1.5259e+06 -1.4622e+06 4.35e-02 7.29e-15 6.50e-14 6.37e+04 1.80e-18 5.70e-01
32 -2.2086e+06 -2.1334e+06 3.52e-02 4.21e-15 1.44e-14 7.52e+04 3.96e-19 8.27e-01
33 -4.7801e+06 -4.4703e+06 6.93e-02 3.04e-15 3.56e-15 3.10e+05 1.05e-19 9.62e-01
34 -5.3184e+06 -5.1445e+06 3.38e-02 1.56e-15 2.05e-15 1.74e+05 3.52e-20 7.65e-01
35 -7.7125e+06 -7.4995e+06 2.84e-02 1.22e-15 3.90e-16 2.13e+05 7.33e-21 8.57e-01
36 -1.4132e+07 -1.3416e+07 5.33e-02 7.65e-15 3.00e-16 7.15e+05 2.17e-21 9.90e-01
37 -2.1317e+07 -2.0862e+07 2.18e-02 7.17e-16 1.48e-15 4.55e+05 2.63e-22 9.40e-01
38 -2.2966e+07 -2.2933e+07 1.44e-03 1.03e-15 2.71e-15 3.30e+04 1.61e-23 9.90e-01
39 -2.8528e+07 -2.8525e+07 9.53e-05 1.26e-15 3.79e-15 2.72e+03 9.85e-25 9.84e-01
40 -2.8636e+07 -2.8635e+07 5.52e-06 6.04e-14 3.87e-15 1.58e+02 5.43e-26 9.60e-01
41 -2.8198e+07 -2.8198e+07 1.19e-07 6.70e-14 3.82e-15 3.37e+00 1.04e-27 9.90e-01
42 -2.8042e+07 -2.8042e+07 3.89e-09 1.05e-14 3.80e-15 1.09e-01 3.43e-29 9.67e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 1.47s