source: AAE 666 Class Notes (Martin Corless) pg. 178/ pg. 265
Consider a disturbed linear system:
\begin{align*} \dot{x} = Ax + Bw\\ z = Cx + Dw \end{align*}where all eigenvalues of A have negative real part and w is a bounded input. Suppose there exiss a positive real scalar $\alpha$ such that:
\begin{align*} \begin{pmatrix} PA + A^TP + \alpha P & PB\\ B^TP & -\alpha \mu_1 I \end{pmatrix} < 0\\ \begin{pmatrix} C^TC - P & C^T D \\ D^TC & D^TD - \mu_2 I \end{pmatrix} < 0 \end{align*}Then $\lim \sup_{t \rightarrow \infty} ||z(t)||_{\infty} \leq \gamma ||w(t)||_{\infty}$ where
\begin{align*} \gamma = \sqrt{\mu_1 + \mu_2} \end{align*}Throughout this paper we will use the following python toolboxes.
In [136]:
import picos as pic
import cvxopt as cvx
import numpy as np
import scipy.linalg
import scipy.optimize
import control
from IPython.core.display import Image
In [14]:
def solve_bounded_disturbance():
def solve_lmi(A_data, B_data, C_data, D_data, alpha, verbose=False):
sdp = pic.Problem()
# variables
P = sdp.add_variable('P', (2,2), vtype='symmetric')
alpha = pic.new_param('alpha', alpha)
mu_1 = sdp.add_variable('mu_1')
mu_2 = sdp.add_variable('mu_2')
# parameters
A = pic.new_param('A', cvx.matrix(A_data))
B = pic.new_param('B', cvx.matrix(B_data))
C = pic.new_param('C', cvx.matrix(C_data))
D = pic.new_param('D', cvx.matrix(D_data))
I = pic.new_param('I', cvx.sparse(cvx.matrix(np.eye(2))))
sdp.add_constraint(
(P*A + A.T*P + alpha*P & P*B) //
(B.T*P & -alpha*mu_1*I) << 0)
sdp.add_constraint(
(C.T*C - P & C.T*D) //
(D.T*C & D.T*D - mu_2*I) << 0)
sdp.add_constraint(P >> 1e-20*I)
sdp.solve(verbose=verbose)
return sdp
A = np.array([[-1,0],[0,-2]])
B = np.eye(2)
C = np.eye(2)
D = np.zeros((2,2))
sdp = solve_lmi(A, B, C, D, 10.0)
# we use fmin to solve for minimum alpha by iteratively solving lmi
alpha_opt = scipy.optimize.fmin(lambda alpha: -alpha if (
alpha > 0 and solve_lmi(A, B, C, D, alpha).status == 'optimal')
else 0.5*abs(alpha), 0)[0]
sdp = solve_lmi(A, B, C, D, alpha_opt)
print sdp
print 'optimal alpha: ', alpha_opt
print 'gamma: ', np.sqrt(
sdp.variables['mu_1'].value[0,0] +
sdp.variables['mu_2'].value[0,0])
solve_bounded_disturbance()
source: AAE 666 Class Notes (Martin Corless)
Consider the system described by \begin{align*} \dot{x} = A(x)x \end{align*} where the n-vector $x$ is the state. The state dependent matrix $A(x)$ has the following structure: \begin{align*} A(x) = A_0 + \Psi(x)\Delta A \end{align*} where $A_0$ and $\Delta A$ are constant $n \times n$ matrices and $\Psi$ is a scalar value function $x$ which is bounded above and below, that is \begin{align*} a \leq \Psi(x) \leq b \end{align*} for some constants a and b. Examples of functions satisfying the above conditions are given by $\Psi(x) = g(c(x))$ where $c(x)$ is a scalar and $g(y)$ is given by $\sin y$, $\cos y$, $e^{−y}$ or $\text{sgm} (y)$. The signum function is useful for modeling switching systems.
Example: \begin{align*} A_0 = \begin{pmatrix} 0 & 1\\ -2 & -1 \end{pmatrix} \\ \Delta A = \begin{pmatrix} 0 & 0\\ 1 & 0 \end{pmatrix} \end{align*} and \begin{align*} \Psi(x) = \begin{cases} \gamma sin x_1/x_1 & x_1 \neq 0 \\ \gamma & x_1=0 \end{cases} \end{align*}
Since $|sin x_1| \leq |x_1|$, we have \begin{align*} -\gamma \leq \Psi(x) \leq \gamma \end{align*} hence $a = -\gamma$ and $b = \gamma$ \begin{align*} A_1 \equiv A_0 + a \Delta A, & A_2 \equiv A_0 + b \Delta A \end{align*}
Suppose there exists a positive-definite symmetric matrix which satisfies the following inequalities: \begin{align*} P A_1 + A_1^T P < 0 P A_2 + A_2^T P < 0 \end{align*} Then the system is globally exponentially stable (GES) about the origin with Lyapunov matrix $P$.
In [131]:
def solve_uncertain(verbose=False):
A = np.array([[0,1],[-2,-1]])
Delta_A = np.array([[0,0],[1,0]])
sdp = pic.Problem()
P = sdp.add_variable('P', (2,2), vtype='symmetric')
A_1 = pic.new_param('A_1', A + Delta_A)
A_2 = pic.new_param('A_2', A - Delta_A)
I = pic.new_param('I', cvx.sparse(cvx.matrix(np.eye(2))))
sdp.add_constraint(P*A_1 + A_1.T*P << 0)
sdp.add_constraint(P*A_2 + A_2.T*P << 0)
sdp.add_constraint(P >> 1e-10*I)
sol = sdp.solve(verbose=verbose)
return sdp
sdp = solve_uncertain()
print sdp
print 'P=\n', sdp.variables['P'].value
source: AAE 666 Class Notes (Martin Corless)
Consider the dynamics given by: \begin{align} \dot{x} &= A(t,x)x + B(t,x)u\\ A(t,x) &= A_0 + \Psi(t,x) \Delta A\\ B(t,x) &= B_0 + \Psi(t,x) \Delta B\\ \end{align} where \begin{align} a \leq \Psi(t,x) \leq b\\ \end{align}
The system can be globally stabilized with rate of convergence $\alpha$ by solving the following optimization problem:
Minimize $\beta$ subject to
\begin{align*} A_1 S + B_1 L + S A_1^T + L^T B_1^T + 2\alpha S \leq 0\\ A_2 S + B_2 L + S A_2^T + L^T B_2^T + 2\alpha S \leq 0\\ I \leq S\\ \begin{pmatrix} \beta I & L^T\\ L & S \end{pmatrix} \ge 0 \end{align*}where $u = LS^{-1}x$
In [133]:
def solve_uncertain_feedback():
def solve_lmi(A, B, Delta_A, Delta_B, alpha, verbose=False):
sdp = pic.Problem()
nx = B.shape[0]
nu = B.shape[1]
S = sdp.add_variable('S', (nx,nx), vtype='symmetric')
L = sdp.add_variable('L', (nu,nx))
beta = sdp.add_variable('beta')
A_1 = pic.new_param('A_1', A + Delta_A)
A_2 = pic.new_param('A_2', A - Delta_A)
B_1 = pic.new_param('B_1', B + Delta_B)
B_2 = pic.new_param('B_2', B - Delta_B)
I_nx = pic.new_param('Inx', cvx.sparse(cvx.matrix(np.eye(nx))))
I_nu = pic.new_param('Inu', cvx.sparse(cvx.matrix(np.eye(nu))))
alpha = pic.new_param('alpha', alpha)
sdp.add_constraint(
A_1*S+ B_1*L + S*A_1.T + L.T*B_1.T + 2*alpha*S << 0)
sdp.add_constraint(#
A_2*S+ B_2*L + S*A_2.T + L.T*B_2.T + 2*alpha*S << 0)
sdp.add_constraint(
(beta*I_nu & L)//
(L.T & S) >> 0)
sdp.add_constraint(S >> 1e-10*I_nx)
sdp.set_objective('min', beta)
K = np.NAN
try:
sdp.solve(verbose=verbose)
if (sdp.status == 'optimal'):
L = np.array(sdp.variables['L'].value)
S = np.array(sdp.variables['S'].value)
K = L.dot(np.linalg.inv(S))
except Exception as e:
pass
return sdp, K
A = np.array([[0,1],[0, 0]])
Delta_A = np.array([[0,0],[1,0]])
B = np.array([[0],[1]])
Delta_B = np.array([[0],[0]])
# we set alpha to set the convergence rate, and then
# compute the gains
for alpha in [0.1, 1, 10, 40, 100]:
sdp, K = solve_lmi(A, B, Delta_A, Delta_B, alpha=alpha)
print 'alpha:\t', alpha, '\tK:\t', K
solve_uncertain_feedback()
Consider the LTI continuous-time system $\dot{x}(t) = Ax(t) + Bu(t)$ where $x \in \mathbb{R}^n$ and $u \in \mathbb{R}^m$ are the system state and the contrl input respectively.
Consider aperiodic sampled-data based state-feedback control laws of the form: $u(t) = Kx(t_k)$, $t \in [t_k, t_{k+1}]$ (2) where $K$ is the controler gain. In order to refine the stability conditions, a fragmentation of the interval $[t_k, t_{k+1}]$ in $N > 0$ disjoint parts is performed and a functional is built for the obtained fragmented system.
The systems (1)-(2) is asymptotically stable for any time-varying sampling period in $[T^-, T^+]$ if there exist constant matrices: $P=P^T \ge 0$, $R_i = R_i^T > 0$, $i=0,\ldots, N-1$, and $U = U^T \in \mathbb{R}^{n \times n}$, $S=S^T \in \mathbb{R}^{nN \times nN}$, $Q \in \mathbb{R}^{nN \times n}$, and $Y \in \mathbb{R}^{n(N+1) \times nN}$ such that the LMIs: \begin{align*} \Psi_1 + T^-(\Psi_2 + \Psi_3) \leq 0\ \Psi_1 + T^+(\Psi_2 + \Psi_3) \leq 0\ \begin{bmatrix} \Psi_1 - T^-\Psi_3 & T^-Y\
with:
\begin{align*} M_i &= \begin{bmatrix} 0_{n \times in} & A & 0_{n \times (N-i-1)n} BK \end{bmatrix}\\ M &= \text{col}_{i=0}^N\{M_i\}\\ N_0 &= \begin{bmatrix} I_n & 0_{n \times N n} \end{bmatrix}\\ N_2 &= \begin{bmatrix} 0_{n \times N n} & I_n \end{bmatrix}\\ \Lambda_1 &= \begin{bmatrix} I_{(N+1)n} & 0_{(N+1)n \times n} \end{bmatrix}\\ \Lambda_2 &= \begin{bmatrix} 0_{(N+1)n\times n} & I_{(N+1)n} \end{bmatrix}\\ \Lambda_{12} &= \Lambda_1 - \Lambda_2\\ \Delta_N &= \text{diag}_{i=0,\ldots,N} \left\{ \sqrt{\frac{N-i}{N}}I_n\right\}\\ \end{align*}Then the system is asmpytotically stable.
This method had some typos in the paper, but I started to code it here:
In [4]:
def solve_async():
n = 2
N = 4
A = np.array([[0, 1], [0, -0.1]])
BK = np.array([[0, 0], [-0.375, -1.15]])
N_0 = np.concatenate((np.eye(n), np.zeros((n,(N)*n))), axis=1)
N_2 = np.concatenate((np.zeros((n,(N)*n)), np.eye(n)), axis=1)
Lambda_1 = np.concatenate((np.eye((N)*n), np.zeros(((N)*n,n))), axis=1)
Lambda_2= np.concatenate((np.zeros(((N)*n,n)), np.eye((N)*n)), axis=1)
Lambda_12 = Lambda_1 - Lambda_2
Delta_N = scipy.linalg.block_diag(*( np.sqrt((N-i)/float(N))*np.eye(n) for i in range(N+1) ))
M = np.zeros((n*N, n*N))
#for i in range(N):
# M[i*n:(i+1)*n, :] = np.concatenate((np.zeros((n,i*n)),A, np.zeros((n,(N-i-1)*n)), BK), axis=1)
M_0 = M[0:n,:]
P = np.ones((n,n))
S = np.ones((n*N,n*N))
Q = np.ones((n*N,n))
Y = np.ones((n*(N+1),n*N))
U = np.ones((n,n))
R_bar = np.ones((n*N,n*N))
R_sig = Delta_N.T.dot(Lambda_1.T.dot(R_bar).dot(Lambda_1) + Lambda_2.T.dot(R_bar).dot(Lambda_2)).dot(Delta_N)
#Psi_1 = 2*N_0.T.dot(P).dot(M_0) - Lambda_12.T.dot(S.dot(Lambda_12) + 2*Q.dot(N_2)) - 2*Y.dot(Lambda_12)
#Psi_2 = M.T.dot(R_sig).dot(M) + 2*(Lambda_12.dot(Delta_N.dot(Delta_N)).dot(M)).T.dot(S.dot(Lambda_12) + 2*Q.dot(N_2))
#Psi_3 = N_2.T.dot(U).dot(N_2)
Source: Destabilizing Strategies for uncertain linear systems
Quadratic Lyapunov functions are conservative. So they are sufficient for stability but not necessary. An interesting application of LMIs is the calculation of polyhedral Lyapunov functions. Polyhedral lyapunov functions can more accurately characterize the dynamics. They are of great theoretical interest because the existence of polyhedral Lyapunov function is necessary and sufficient for a switched or uncretain linear system.
Consider the dynamics
\begin{align*} \dot{x} = A(t)x \end{align*}with $x \in \mathbb{R}^n$, and uncertainty described by a polytope of matrices
\begin{align*} A(t) = \sum_{i=1}^K u_i(t)A_i \end{align*}where \begin{align*} u_i(t) \ge 0,& & \sum_{i=1}^K u_i(t) = 1 \end{align*}
The polyhedral Lyapunov function can be expressed as:
\begin{align*} F(m) = \left\{ x: \begin{bmatrix} W \\ -W \end{bmatrix} x \leq \mathbb{1}^-(M+m) \right\} \end{align*}where each row $w_i$ of $W$ corresponds to the nonempty and nondegenerate face of the polyhedron $F(m)$, and $\mathbb{1}^-(M+m)$ denotes a $2M$ dimensional column vector with all elements equal to 1 except the $M+m$ element which is equal to $-1$. The rows of $\mathbb{1}^-$ will be denoted $h_j$. Finding a destabilizig strategy can be formulated as LMI as detailed in Polanski's paper.
The code for this analysis is too large to include in this notebook, but the plots/ LP output have been included below for the case in the paper.
In [153]:
Image('13/quadratic.png')
Out[153]:
Below, a polyhedral Lyapunov function is shown for the same dynamics. Note that the polyhedral Lyapunov function more accurately captures the stability.
In [139]:
Image('13/polanski_switch.png')
Out[139]:
To constrain control input when the initial condition is known, an LMI can be added to a problem to a quadratic stability problem easily. Consider $u(t) = Kx(t)$. Pick $Q > 0$ and $Y$ which satisfy the quadratic stabilizability condtion $AQ + QA^T + B_u^T \leq 0$, and in addition $x(t)^T Q^{-1} x(0) \leq 1$.
\begin{align*} \max_{t \geq 0} ||u(t)|| &= \max_{t \geq 0} || Y Q^{-1} x(t) ||\\ &\leq max_{x \in \xi} || Y Q^{-1} x ||\\ &= \lambda_{max} (Q^{-1/2}Y^TYQ^{-1/2}) \end{align*}Therefore, the constraint $||u(t)|| \leq \mu$ is enforced at all times $t \geq 0$ if the following LMIs hold: \begin{align*} \begin{bmatrix} 1 & x(0)^T \\ x(0) & Q \end{bmatrix} \geq 0 & & \begin{bmatrix} Q & Y^T \\ Y & \mu^2 I \end{bmatrix} \geq 0 \end{align*}
Source Linear Matrix Inequalities in System and Control Theory (Boyd)
The condition number of a matrix is the ratio of its largest and smallest singular value:
\begin{align*} \kappa(M) \equiv \left(\frac{\lambda_{max}(M^T M)}{\lambda_{min}(M^T M)} \right)^{1/2} \end{align*}Which simplifies for square invertible matrices to: $\kappa(M) = ||M|| ||M^{-1}||$
We want to scale the matrix by $L$ and $R$ which are diagonal and nonsingular.
minimize:
$\kappa (LMR)$
subject to:
$L \in \mathbb{R}^{p\times p}$, diagonal and nonsingular $R \in \mathbb{R}^{q\times q}$, diagonal and nonsingular
This is equivalent to the existence of diagonal $P \in \mathbb{R}^{p\times p}$, $Q \in \mathbb{R}^{q\times q}$ with $P > 0$ and $Q >0$ and $Q \leq M^T P M \leq \gamma^2 Q$. (see Boyd)
Which can the be formulated as a standard LMI:
minimize: $\gamma^2$
subject to:
\begin{align*} P \in \mathbb{R}^{p\times p}, & \text{diagonal}, & P > 0\\ Q \in \mathbb{R}^{q\times q}, & \text{diagonal}, & Q > 0\\ Q \leq &M^T P M \leq \gamma^2 Q \end{align*}Source Linear Matrix Inequalities in System and Control Theory (Boyd)
Consider the system described by the delay-differential equation
\begin{align*} \dot{x}(t) = A x(t) + \sum_{i=1}^L A_ix(t - \tau_i) \end{align*}where $x(t) \in \mathbb{R}^n$, and $\tau_i > 0$.
If the functional \begin{align*} V(x,t) = x(t)^TPx(t) + \sum_{i=1}^L \int_0^{\tau_i} x(t-s)^TP_i x(t-s) ds \end{align*} where $P>0$, $P_1 > 0$, $\ldots$, $P_L > 0$, satisfies $dV(x,t)/dt < 0$ for every x satisfying the dynamics is stable, it can be verified that $dV(x,t)/dt = y(t)^TWy(t)$, where
\begin{align*} W = \begin{bmatrix} A^TP + PA + \sum_{i=1}^L P_i & P A_1 & \ldots & PA_L\\ A_1^TP & -P_1 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ A_L^T P & 0 & \ldots & -P_L\\ \end{bmatrix} \end{align*}\begin{align*} y(t) = \begin{bmatrix} x(t) \\ x(t - \tau_q) \\ \vdots \\ x(t - \tau_L) \end{bmatrix} \end{align*}Therefore we can prove stability of the system using Lyapunov functiionals of the form $V(x,t)$, by solving the LMIP $W<0$, $P>0$, $P_1>0$, $\ldots$, $P_L$ > 0
In [161]:
%%bash
curdir=$(pwd)
author="James Goppert"
date="4/14/2014"
title="Control Applications of LMIs"
notebook="Quiz 13.ipynb"
mkdir -p tmp; cd tmp
echo """
((*- extends 'article.tplx' -*))
((* block margins *))
\geometry{top=1cm,bottom=2cm,left=1cm,right=1cm}
((* endblock margins *))
((* block py_err *))
((* endblock py_err *))
((* block title *))
\\title{$title}
((* endblock title *))
((* block date *))
\date{$date}
((* endblock date *))
((* block author *))
\author{$author}
((* endblock author *))
""" > notebook.tplx
ipython nbconvert --quiet --to=latex \
--template="notebook.tplx" \
--post=PDF "$curdir/$notebook" #&> /dev/null
In [162]:
%%bash
evince ./tmp/Quiz\ 13.pdf &> /dev/null
In [ ]: