| title | author | date |
|---|---|---|
Convex Optimization: A Practical Guide |
Behzad Samadi |
February 2, 2014 |
Behzad Samadi
http://www.mechatronics3d.com/
February 2, 2014
Optimization problems happen in many fields including engineering, economics and medicine. The objective of many engineering problems is to design "the best" product. "The best" needs a formal definition that is not subject to personal judgement as much as possible. In a mathematical optimization, the objective function defines "the best". The objective function can be the description of a cost function we want to minimize or it can be a quantified description of the product's fitness, which is required to be maximized. We live in a real world with lots of constraints. The energy is not free. The time we can spend to design the product is limited. The computation power we have is limited. The size, weight and price of the product is limited. Therefore, the optimization problems are usually defined as minimizing or maximizing an objective function considering a set of constraints. In this text, we focus on a certain class of optimization problems: convex optimization. The main importance of convex optimization problems is that there is no locally optimum point. If a given point is locally optimal then it is globally optimal. In addition, there exist effective numerical methods to solve convex optimization problems. In addition, it is possible to convert many nonconvex optimization problems to convex problems by changing the variables or introducing new variables. In this text, we will review what convex sets, functions and optimization problems are. Also, we show you numerical examples and applications of convex optimization in control systems. Also, there are many examples with the corresponding code in Python to help the reader understand how the problems are solved in practive. The Python code is based on cvxpy.
In the following, the definitions are taken from [1] unless otherwise stated. The reader is referred to the Convex Optimization book by Stephen Boyd and Lieven Vandenberghe for a detailed review of the theory of convex optimization and applications.
Definition. Convex combination: Given
Definition. Convex set: A set
Definition. Convex hull: The convex hull of a set
Definition. Affine combination:
Definition. Affine set: A set
Definition. Cone (nonnegative) combination: Cone combination of two points
Definition. Convex cone: A set
Definition. Hyperplane: A hyperplane is a set of the form
Definition. Halfspace: A halfspace is a set of the form
Definition. Polyhedron: A polyhedron is the intersection of finite number of
hyperplanes and halfspaces. A polyhedron can be written as:
\begin{equation}
\mathcal{P}={x| Ax \preceq b, Cx=d }
\end{equation}
where
Definition. Euclidean ball: A ball with center
Definition. Ellipsoid: An ellipsoid is defined as:
\begin{equation}
{x | (x-x_c)^\text{T}P^{-1}(x-x_c)\leq 1}
\end{equation}
where
Definition. Proper code: A cone is proper if it is closed (contains its boundary), solid (has nonempty interior) and pointed (contains no lines).
The nonnegative orthant of
Definition. Generalized inequality: A generalized inequality is defined by a
proper cone
In this context, we deal with the following inequalities:
- The inequality on real numbers is defined based on the proper cone of
nonnegative real numbers
$K=\mathbb{R}_+$ . - The componentwise inequality on real vectors in
$\mathbb{R}^n$ is defined based on the nonnegative orthant$K=\mathbb{R}^n_+$ . - The matrix inequality is defined based on the proper cone of positive
semidefinite matrices
$K=S^n_+$ .
Definition. Convex function: A function
A mathematical optimization is convex if the objective is a convex function and the feasible set is a convex set. The standard form of a convex optimization problem is: \begin{align} \text{minimize } & f_0(x) \nonumber\ \text{subject to } & f_i(x) \leq 0,\ i=1,\ldots,m\nonumber\ & h_i(x) = 0,\ i=1,\ldots,p \end{align}
Linear programming (LP) is one of the best known forms of convex optimization. A
LP problem can be written as:
\begin{align}\label{LP}
\text{minimize }&c^\text{T}x\nonumber\
\text{subject to }&a_i^\text{T}x\leq b_i,\ i=1,\ldots,m
\end{align}
where
If some of the entries of
In general, the feasible set of a linear programming is a polyhedron. The
objective function defines a family of parallel hyperplanes. The optimal value
for the objective function is the lowest value corresponding to a hyperplane
that has a non-empty intersection with the feasible set polyhedron. The
intersection can be a vertice or edge or any higher dimensional faces.
Therefore, the optimal value of the objective function is unique but the optimal
solution,
Example: Consider the following LP problem (LP1):
\begin{align} \text{maximize: } & x + y\nonumber\ \text{Subject to: } & x + y \geq -1 \ \text{} & \frac{x}{2}-y \geq -2\nonumber\ \text{} & 2x-y \leq -4\nonumber \end{align}
In order to solve this LP problem in Python, we need to import the required modules:
import numpy as np
from pylab import *
import matplotlib as mpl
import cvxopt as co
import cvxpy as cp
%pylab --no-import-all inline
Populating the interactive namespace from numpy and matplotlib
The next step is to define the optimization variables:
x = cp.Variable(1)
y = cp.Variable(1)
The constraints are then added:
constraints = [ x+y >= -1.,
0.5*x-y >= -2.,
2.*x-y <= 4.]
Then, the objective function and the optimization problem are defined as:
objective = cp.Maximize(x+y)
p = cp.Problem(objective, constraints)
The solution of the LP problem is computed with the following command:
result = p.solve()
print(round(result,5))
8.0
The optimal solution is now given by:
x_star = x.value
print(round(x_star,5))
4.0
y_star = y.value
print(round(y_star,5))
4.0
The feasible set of the LP problem (ref{LP1}) is shown in Figure ref{LPfeas}, which is drawn using the following commands:
xp = np.array([-3, 6])
# plot the constraints
plt.plot(xp, -xp-1, xp, 0.5*xp+2, xp, 2*xp-4)
# Draw the lines
plt.plot(xp, -xp+4, xp, -xp+8)
# Draw the feasible set (filled triangle)
path = mpl.path.Path([[4, 4], [1, -2], [-2, 1], [4, 4]])
patch = mpl.patches.PathPatch(path, facecolor='green')
# Add the triangle to the plot
plt.gca().add_patch(patch)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Feasible Set of a Linear Program')
plt.xlim(-3,6)
plt.ylim(-5,5)
plt.text(2, -4.5, "x+y=-1")
plt.text(3.5, -0.75, "x+y=4")
plt.text(4.5, 2.25, "x+y=8")
plt.show()
Now, to solve the following LP problem (LP2):
\begin{align} \text{minimize: } & x + y\nonumber\ \text{Subject to: } & x + y \geq -1 \ \text{} & \frac{x}{2}-y \leq -2\nonumber\ \text{} & 2x-y \leq -4\nonumber \end{align}
we change the objective function in the code:
objective = cp.Minimize(x+y)
p = cp.Problem(objective, constraints)
result = p.solve()
print(round(result,5))
-1.0
The optimal solution is now given by:
x_star = x.value
print(round(x_star,5))
0.49742
y_star = y.value
print(round(y_star,5))
-1.49742
In this case the optimzal value of the objective function is unique. However, it can be seen in Figure ref{LPfeas} that any point on the line connecting the two points (-2,1) and (1,-2) including the point (0.49742,-1.49742) can be the optimal solution. Therefore, the LP problem ref{LP2} has infinite optimal solutions. The code, however, returns just one of the optimal solutions.
Example: Finding the Chebyshev center of a polyhedron is an example of
optimization problems that can be solved using LP \cite{cvx}. However, the
original description of the problem is not in LP form. Consider the following
polyhedron:
\begin{equation}
\mathcal{P} = {x | a_i^Tx \leq b_i, i=1,...,m }
\end{equation}
The Chebyshev center of
r = cp.Variable(1)
xc = cp.Variable(2)
a1 = co.matrix([-1,-1], (2,1))
a2 = co.matrix([-0.5,1], (2,1))
a3 = co.matrix([2,-1], (2,1))
b1 = 1
b2 = 2
b3 = 4
constraints = [ a1.T*xc + np.linalg.norm(a1, 2)*r <= b1,
a2.T*xc + np.linalg.norm(a2, 2)*r <= b2,
a3.T*xc + np.linalg.norm(a3, 2)*r <= b3 ]
objective = cp.Maximize(r)
p = cp.Problem(objective, constraints)
result = p.solve()
The radius of the ball is:
print r.value
1.52896116777
and the Chebyshev center is located at:
print xc.value
[ 5.81e-01]
[ 5.81e-01]
The triangle and the largest circle that it can include are depicted in Figure ref{Cheb} using the following commands:
xp = np.linspace(-3, 5, 256)
theta = np.linspace(0,2*np.pi,100)
# plot the constraints
plt.plot( xp, -xp*a1[0]/a1[1] + b1/a1[1])
plt.plot( xp, -xp*a2[0]/a2[1] + b2/a2[1])
plt.plot( xp, -xp*a3[0]/a3[1] + b3/a3[1])
# plot the solution
plt.plot( xc.value[0] + r.value*cos(theta), xc.value[1] + r.value*sin(theta) )
plt.plot( xc.value[0], xc.value[1], 'x', markersize=10 )
plt.title('Chebyshev Center')
plt.xlabel('x1')
plt.ylabel('x2')
plt.axis([-3, 5, -3, 5])
plt.show()
A linear fractional program is defined as (LFP):
\begin{align}
\text{minimize }&f_0(x)\nonumber\
\text{subject to }&a_i^\text{T}x\leq b_i,\ i=1,\ldots,m\nonumber\
&h_i^\text{T}x = g_i,\ i=1,\ldots,p
\end{align}
with
Another approach is introduce auxiliary variables
A closely related class of optimization problems is the generalized linear
fraction problems formulated as (GLFP):
\begin{align}
\text{minimize }&f_0(x)\nonumber\
\text{subject to }&a_i^\text{T}x\leq b_i,\ i=1,\ldots,m\nonumber\
&h_i^\text{T}x = g_i,\ i=1,\ldots,p
\end{align}
where:
\begin{equation}
f_0(x)=\max_{i=1,\ldots,r} \frac{c_i^\text{T}x+d_i}{e_i^\text{T}x+f_i}
\end{equation}
and the domain of
A quadratic programming (QP) optimization problem is described as:
\begin{align}
\text{minimize }&\frac{1}{2}x^\text{T}Px+q^\text{T}x+r\nonumber\
\text{subject to }&Gx \preccurlyeq h \nonumber\
&Ax = b
\end{align}
If the objective function is quadratic and the constraints include quadratic
constraints, then we have a quadratically constrained quadratic program (QCQP):
\begin{align}
\text{minimize }&\frac{1}{2}x^\text{T}P_0x+q_0^\text{T}x+r_0\nonumber\
\text{subject to }&\frac{1}{2}x^\text{T}P_ix+q_i^\text{T}x+r_i\leq 0,
i=1\cdots,m \nonumber\
&Ax = b
\end{align}
where
Example: Consider the set of linear equations
A = np.array([[1,1],[2,1],[3,2]])
b = np.array([[2], [3], [4]])
xs = np.dot(np.linalg.pinv(A),b)
print(xs)
[[ 1. ]
[ 0.66666667]]
A similar result can be reached by solving the following QP problem:
A = co.matrix(A)
b = co.matrix(b)
x = cp.Variable(2)
I = np.identity(3)
objective = cp.Minimize( cp.quad_form(A*x-b, I) )
p = cp.Problem(objective)
The optimal value of the objective function is:
result = p.solve()
print(result)
0.333333326323
and the optimal solution is:
print(x.value)
[ 1.00e+00]
[ 6.67e-01]
Now, we can add linear constraints and find the optimal solution by solving the QP problem:
x = cp.Variable(2)
objective = cp.Minimize( b.T*b - 2*b.T*A*x + cp.quad_form(x, A.T*A) )
constraints = [ -0.9 <= x <= 0.9]
p = cp.Problem(objective, constraints)
The optimal cost function is equal to:
result = p.solve()
print(result)
0.33838157573
which is more than what it was without the linear constraints. The optimal solution is equal to:
print(x.value)
[ 9.00e-01]
[ 8.18e-01]
Example (Linear Program with a Stochastic Objective Function): Consider a
random vector
One way to formulate the problem so that it is practically solveable is to set
the objective function as:
\begin{equation}
{\bar c}^\text{T}x+\gamma x^\text{T}\Sigma x
\end{equation}
where
Now, the optimization can be solved with the following code:
Sigma = co.matrix([ 5, 1,
1, 4,], (2,2))
cb = co.matrix([1, 1], (2,1))
G = co.matrix([ -1, -0.5, 2,
-1, 1, -1], (3,2))
h = co.matrix([ 1,
2,
4],(3,1))
gamma = 0.5
x = cp.Variable(2)
objective = cp.Minimize( cb.T * x + gamma * cp.quad_form(x, Sigma) )
constraints = [ G*x <= h ]
p = cp.Problem(objective, constraints)
The optimal value of the objective function is:
result = p.solve()
print(result)
-0.184210521637
The optimal solution, in this case, is inside of the feasible set:
print(x.value)
[-1.58e-01]
[-2.11e-01]
In the following figure, the feasible set and the contours of the objective function are drawn.
xp = np.array([-3, 6])
# plot the constraints
plt.plot(xp, -xp-1, xp, 0.5*xp+2, xp, 2*xp-4)
# Draw the feasible set (filled triangle)
path = mpl.path.Path([[4, 4], [1, -2], [-2, 1], [4, 4]])
patch = mpl.patches.PathPatch(path, facecolor='green')
# Add the triangle to the plot
plt.gca().add_patch(patch)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Feasible Set of a Linear Program')
plt.xlim(-3,6)
plt.ylim(-5,5)
delta = 0.025
xc = np.arange(-3.0, 6.0, delta)
yc = np.arange(-5.0, 6.0, delta)
X, Y = np.meshgrid(xc, yc)
Z = cb[0]*X + cb[1]*Y + Sigma[0,0]*X*X + 2*Sigma[0,1]*X*Y + Sigma[1,1]*Y*Y
plt.contour(X, Y, Z)
X, Y = np.meshgrid(x, y)
plt.show()
Example (Distance between polyhedra): Consider the following two polyhedra:
\begin{equation}
\mathcal{P}_1={x| A_1x \preccurlyeq b_1},\ \mathcal{P}_2={x| A_2x
\preccurlyeq b_2}
\end{equation}
The distance between
# Draw the triangle
path = mpl.path.Path([[4, 4], [1, -2], [-2, 1], [4, 4]])
patch = mpl.patches.PathPatch(path, facecolor='green')
# Add the triangle to the plot
plt.gca().add_patch(patch)
# Draw the rectangle
path = mpl.path.Path([[-3, 4], [-2, 4], [-2, 2], [-3, 2], [-3, 4]])
patch = mpl.patches.PathPatch(path, facecolor='green')
# Add the rectangle to the plot
plt.gca().add_patch(patch)
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(-4,6)
plt.ylim(-3,5)
plt.show()
The distance between these two polygons is computed with the following QP optimization problem:
x1 = cp.Variable(2)
x2 = cp.Variable(2)
I = np.identity(2)
# Triangle
A1 = co.matrix([ -1, -0.5, 2,
-1, 1, -1], (3,2))
b1 = co.matrix([ 1,
2,
4],(3,1))
# Rectangle
A2 = co.matrix([ -1, 1, 0, 0,
0, 0, -1, 1], (4,2))
b2 = co.matrix([ 3,
-2,
-2,
4],(4,1))
objective = cp.Minimize( cp.quad_form(x1-x2, I) )
constraints = [ A1*x1<= b1, A2*x2<=b2]
p = cp.Problem(objective, constraints)
The distance between the two polygons is:
result=p.solve()
print(result)
0.799999994069
The correspondin point in the triangle is:
print(x1.value)
[-1.60e+00]
[ 1.20e+00]
and the corresponding point in the rectangle is:
print(x2.value)
[-2.00e+00]
[ 2.00e+00]
A second order cone program (SOCP) is defined as:
\begin{align}
\text{minimize }& f^\text{T}x\nonumber\
\text{subject to }&|A_ix+b_i|_2\leq c_i^\text{T}x+d_i,
i=1,\ldots,m\nonumber\
&Fx = g
\end{align}
Note that:
- If
$c_i=0$ for$i=1,\ldots,m$ , the SOCP is equivalent to a QP. - If
$A_i=0$ for$i=1,\ldots,m$ , the SOCP is equivalent to a LP.
Example (robust linear program): Consider the following LP problem:
\begin{align}
\text{minimize }&c^\text{T}x\nonumber\
\text{subject to }& a_i^\text{T}x \leq b_i,\ i=1,\ldots,m
\end{align}
where the parameters are assumed to be uncertain. For simplicity, let us assume
that
i=1,\ldots,m
\end{align}
Example (stochastic linear program): The same robust LP problem can be
addressed in a stochastic framework. In this framework,
Assuming that
Example: Consider the equation
x = cp.Variable(2)
t = cp.Variable(1)
t1 = cp.Variable(1)
t2 = cp.Variable(1)
gamma = 0.5
A = co.matrix([[1,2,3],[1,1,2]])
b = co.matrix([2, 3, 4],(3,1))
objective = cp.Minimize(t)
constraints = [cp.norm(A*x-b)+gamma*(t1+t2) <= t, -t1 <= x[0] <= t1, -t2 <= x[1] <= t2 ]
p = cp.Problem(objective, constraints)
The optimal value of the objective function is:
result=p.solve()
print(result)
1.36037960817
and the optimal solution is:
print(x.value)
[ 1.32e+00]
[ 1.40e-01]
Thanks to cvxpy, the same problem can be solved using a much shorter code:
p = cp.Problem(cp.Minimize(cp.norm(A*x-b)+gamma*cp.norm(x,1)), [])
result=p.solve()
print(x.value)
[ 1.32e+00]
[ 1.40e-01]
Generalized inequalities can be defined based on propoer cones. Till now we have
seen inequalities for real numbers and elementwise inequalities real vectors.
The former type of inequality is defined by the propoer cone of nonnegative real
numbers. The later type is defined by the proper cone of the nonnegative orthant
(
A semidefinite program (SDP) is defined as:
\begin{align}
\text{minimize} & c^\text{T}x\nonumber\
\text{subject to}& x_1F_1+x_2F_2+\cdots+x_nF_n+G \leq 0\nonumber\
& Ax=b
\end{align}
where $x =\left[\begin{matrix}x_1&x_2&\cdots&x_n\end{matrix}\right]^\text{T}$ is
a vector in
Note that the standard LP problem:
\begin{align}
\text{minimize}&c^\text{T}x\nonumber\
&Ax\preccurlyeq b
\end{align}
can be written as the following SDP problem:
\begin{align}
\text{minimize}&c^\text{T}x\nonumber\
&\text{diag}(Ax) \leq b
\end{align}
where
Also the standard SOCP problem: \begin{align} \text{minimize}&c^\text{T}x\nonumber\ &|A_ix+b_i|_2 \leq c_i^\text{T}x+d_i,\ i=1,\ldots,m \end{align} can be written as the following SDP problem: \begin{align} \text{minimize}&c^\text{T}x\nonumber\ & \left[\begin{matrix} (c_i^\text{T}x+d_i)I & A_ix+b_i\ A_ix+b_i & c_i^\text{T}x+d_i\end{matrix}\right] \geq 0,\ i=1,\ldots,m \end{align}
Example (eigenvalue minimization): Consider minimizing the maximum eigenvalue
of matrix
x1 = cp.Variable()
x2 = cp.Variable()
x3 = cp.Variable()
t = cp.Variable()
X = cp.SDPVar(4)
Y = cp.SDPVar(4)
A0 = co.matrix([[0,0,0,-2],[0,0,0,-3],[0,0,0,-4],[-2,-3,-4,0]])
A1 = co.matrix([[0,0,0,1],[0,0,0,2],[0,0,0,3],[1,2,3,0]])
A2 = co.matrix([[0,0,0,1],[0,0,0,1],[0,0,0,2],[1,1,2,0]])
A3 = co.matrix([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])
I = co.matrix(np.identity(4))
objective = cp.Minimize(t)
Ax = A0+A1*x1+A2*x2+A3*x3
constraints = [ Ax == X, t*I-Ax == Y ]
p = cp.Problem(objective, constraints)
The optimal largest eigenvalue of
result=p.solve()
print(result)
1.15470050608
The optimal solution for
print([x1.value,x2.value,x3.value])
[1.0000000000000004, 0.6666666666666659, 0.5773502530380881]
The value of
print(Ax.value)
[ 5.77e-01 0.00e+00 0.00e+00 -3.33e-01]
[ 0.00e+00 5.77e-01 0.00e+00 -3.33e-01]
[ 0.00e+00 0.00e+00 5.77e-01 3.33e-01]
[-3.33e-01 -3.33e-01 3.33e-01 5.77e-01]
In this case, eigenvalues of
EigenValues = np.linalg.eig(Ax.value)[0]
print(EigenValues)
[ -1.61515377e-08 5.77350253e-01 1.15470052e+00 5.77350253e-01]
Many problems in control theory can be formulated as convex optimization problems. It is beyond the scope of this text to cover all of them. However, in the following, you will see an introduction about applications of convex optimization in control theory.
Consider the following autonomous system:
\begin{equation}
\dot x = f(x),\ x(0)=x_0
\end{equation}
with
Definition [2]: An equilibrium point
-
stable (in the sense of Lyapunov) if for any given
$\epsilon>0$ , there exists$\delta(\epsilon)>0$ such that: \begin{equation} |x_0-x^\star|\leq \delta \Rightarrow \lim_{t\rightarrow \infty}|\phi(t,x_0)-x^\star|=0 \end{equation} -
attractive if there exists
$\delta>0$ such that: \begin{equation} |x_0-x^\star|\leq\delta \Rightarrow \lim_{t\rightarrow \infty}|\phi(t,x_0)-x^\star|=0 \end{equation} - asymptotically stable (in the sense of Lyapunov) if it is both stable and attractive.
Theorem [3]: If there exists a continuous function
Consider the linear system:
\begin{equation}
\dot x=Ax, x(0)=x_0
\end{equation}
where
Example: Consider a linear system with:
\begin{equation}
A=\left[\begin{matrix}0&1\-1&-2\end{matrix}\right]
\end{equation}
The eigenvalues of
A = np.matrix('0 1; -1 -2')
np.linalg.eig(A)[0]
array([-1., -1.])
Therefore, the system is stable. In the following, we now verify the stability of the system by solving the following LMIs:
A = cp.Parameter(2,2)
P = cp.SDPVar(2)
objective = cp.Minimize( 0 )
prob = cp.Problem(objective, [A.T*P+P*A == -cp.SDPVar(2,2), P == cp.SDPVar(2,2)])
A.value = co.matrix(np.matrix('0 1; -1 -2'))
prob.solve()
0.0
where SDPVar(n,n) is an auxiliary positive semi-definite matrix. Since the stability LMIs are a feasibility problem as opposed to an optimization problem, we have set the objective function to be a constant value. The optimal solution of the above LMIs is:
print(P.value)
[ 0.00e+00 0.00e+00]
[ 0.00e+00 0.00e+00]
However, this is the trivial answer of these LMIs: \begin{align} P\geq 0\nonumber\ A^\text{T}P+PA\leq 0 \end{align} In order to find a feasible answer for strict inequalities using non-strict inequalities, we can rewrite the inequalities as: \begin{align} P-\epsilon I\geq 0\nonumber\ A^\text{T}P+PA+\alpha P\leq 0 \end{align} Let us now solve the following LMIs to find a valid Lyapunov function for the linear system:
I = np.identity(2)
prob = cp.Problem(objective, [A.T*P+P*A+0.5*P == -cp.SDPVar(2,2), P - 0.1*I == cp.SDPVar(2,2)])
prob.solve()
0.0
print(P.value)
[ 1.13e+00 1.88e+00]
[ 4.93e-02 9.40e-01]
It is also possible to add an objective function to the LMIs, for example
objective = cp.Minimize( cp.lambda_max(P) )
prob = cp.Problem(objective, [A.T*P+P*A+0.5*P == -cp.SDPVar(2,2), P - 0.1*I == cp.SDPVar(2,2)])
prob.solve()
0.1777777617758615
The optimal solution is now equal to:
print(P.value)
[ 1.39e-01 3.89e-02]
[ 3.89e-02 1.39e-01]
We can verify the inequalities by computing the eigenvalues:
np.linalg.eig(P.value)[0]
array([ 0.17777777, 0.1 ])
np.linalg.eig( (A.T*P+P*A).value )[0]
array([-0.06318906, -0.4923496 ])
Consider the following linear system:
\begin{equation}
\dot x=A(\alpha)x, x(0)=x_0
\end{equation}
where
Example: Consider the uncertain linear system () with
A_2=\left[\begin{matrix}1&2\-2&-2\end{matrix}\right]
\end{equation}
The eigenvalues of
A1 = np.matrix('1 -2; 2 -2')
np.linalg.eig(A1)[0]
array([-0.5+1.32287566j, -0.5-1.32287566j])
A2 = np.matrix('1 2; -2 -2')
np.linalg.eig(A2)[0]
array([-0.5+1.32287566j, -0.5-1.32287566j])
However, an uncertain linear system with
np.linalg.eig(0.5*A1+0.5*A2)[0]
array([ 1., -2.])
That is a proof for the following LMIs to be infeasible:
A1 = cp.Parameter(2,2)
A2 = cp.Parameter(2,2)
P = cp.SDPVar(2)
I = np.identity(2)
objective = cp.Minimize( cp.lambda_max(P) )
prob = cp.Problem(objective, [P - 0.1*I == cp.SDPVar(2,2), A1.T*P+P*A1+0.5*P == -cp.SDPVar(2,2), A2.T*P+P*A2+0.5*P == -cp.SDPVar(2,2)])
A1.value = co.matrix(np.matrix('1 -2; 2 -2'))
A2.value = co.matrix(np.matrix('1 2; -2 -2'))
prob.solve()
'infeasible'
Consider the following linear system:
\begin{equation}
\dot x=Ax+Bu
\end{equation}
where
We know that the eigenvalues of a matrix and its transpose are the same.
Therefore, the closed loop system () is stable if and only if its dual system is
stable:
\begin{equation}
\dot x = (A+BK)x
\end{equation}
Now, let us write the stability inequalities for the dual system:
\begin{align}
{\color{red}Q}>0\nonumber\
(A+B{\color{red}K}){\color{red}Q}+{\color{red}Q}(A+B{\color{red}K})^\text{T} < 0
\end{align}
This is still a BMI. However, if we define
Example: Consider the following linear system:
\begin{equation}
\dot x = \left[\begin{matrix}1&0.1\0&-2\end{matrix}\right]x+\left[\begin{matrix
}0\1\end{matrix}\right]u
\end{equation}
The objective is to find
A = cp.Parameter(2,2)
B = cp.Parameter(2,1)
Q = cp.SDPVar(2)
Y = cp.Variable(1,2)
I = np.identity(2)
objective = cp.Minimize( cp.lambda_max(Q) )
prob = cp.Problem(objective, [Q - 0.1*I == cp.SDPVar(2,2), A*Q+Q*A.T+B*Y+Y.T*B.T+0.5*Q == -cp.SDPVar(2,2)])
A.value = co.matrix(np.matrix('1 0.1; 0 -2'))
B.value = co.matrix(np.matrix('0; 1'))
prob.solve()
62.699830858476346
The controller gain can now be computed as:
K = np.dot(Y.value,np.linalg.inv(Q.value))
print(K)
[[-54.46273825 -1.34901911]]
The closed loop system is stable:
Acl = (A.value+B.value*K)
np.linalg.eig(Acl)[0]
array([-1.17450955+0.84722018j, -1.17450955-0.84722018j])
Similar to stability for autonomous systems, there is a concept called
dissipativity for dynamic systems with input. Consider the following dynamical
system:
\begin{align}
\dot x=&f(x,w)\nonumber\
z=&g(x,w)
\end{align}
where
Definition: The system () is said to be dissipative with storage function
Now, consider the following linear system:
\begin{align}
\dot x=&Ax+Bw\nonumber\
z=&Cx+Dw
\end{align}
with
The following statements are equivalent:
- The system is strictly dissipative with the supply rate: \begin{equation} s(w,z)=\left[\begin{matrix} w\\z\end{matrix}\right]^\text{T} \left[\begin{matrix} Q&S\\S^\text{T}& R\end{matrix}\right] \left[\begin{matrix} w\\z \end{matrix}\right] \end{equation}
- For all
$\omega\in\mathbb{R}\cup{\infty}$ there holds: \begin{equation} \left[\begin{matrix} I\\T(j\omega)\end{matrix}\right]^\star \left[\begin{matrix} Q& S\\S^\text{T}&R\end{matrix}\right]\left[\begin{matrix} I\\ T(j\omega)\end{matrix}\right] \gt 0 \end{equation} - There exists
$P>0$ satisfying the LMI: \begin{equation} \left[\begin{matrix} A^\TR P+PA & PB\\B^\TR P & 0 \end{matrix}\right]-\left[\begin{matrix} 0&I\\C&D\end{matrix}\right]^\text{T} \left[\begin{matrix} Q& S\\S^\text{T}&R\end{matrix}\right]\left[\begin{matrix} 0&I\\ C&D\end{matrix}\right] \lt 0 \end{equation}
The following statements are equivalent:
- System () is strictly dissipative with the supply rate: \begin{equation} s(w,z)= w^\text{T} z+z^\text{T} w=2w^\text{T} z \end{equation}
- For all
$\omega\in\mathbb{R}$ with$\det(j\omega I-A)\neq 0$ , there holds: \begin{equation} T(j\omega)^\star+T(j\omega)>0 \end{equation} - There exists
$P>0$ satisfying the LMI: \begin{equation} \left[\begin{matrix} A^\text{T} P+PA & PB-C^\text{T}\\B^\text{T} P-C & D+D^\text{T} \end{matrix}\right] < 0 \end{equation} - If
$D=0$ , there exists$P>0$ satisfying: \begin{align} A^\text{T}P+PA<0\nonumber\ PB=C^\text{T} \end{align} - System () is RLC realizable, i.e. there exists an RLC network with transfer
function
$T(j\omega)$ . - For SISO systems: \begin{equation} |\angle G(j\omega)| < 90^\circ, \text{Re}(G(j\omega))>0 \end{equation}
- Gain margin of system () is infinite.
One of the properties of passive systems is that the feedback connection of two passive systems is always stable (the loop phase is less than 180 degrees).
The following statements are equivalent:
- The system is strictly dissipative with the supply rate: \begin{equation} s(w,z)=\gamma^2 w^\text{T} w-z^\text{T} z \end{equation}
- For all
$\omega\in\mathbb{R}$ : \begin{equation} |T(j\omega)|\infty=\sup{0<|w|_2<\infty}\frac{|z|_2}{|w|_2} \lt \gamma \end{equation} - There exists
$P>0$ satisfying the LMI: \begin{equation} \left[\begin{matrix} A^\text{T} P+PA+C^\text{T} C & PB+C^\text{T} D\\B^\text{T} P+D^\text{T} C & D^\text{T} D-\gamma^2 I \end{matrix}\right]< 0 \end{equation} - There exists
$P>0$ satisfying the LMI: \begin{equation} \left[\begin{matrix} A^\text{T} P+PA & PB & C^\text{T}\B^\text{T} P & -\gamma^2 I & D^\text{T}\C&D&-I \end{matrix}\right]< 0 \end{equation}
Consider the following linear system:
\begin{align}
\dot x=&Ax+B_ww+B_uu\nonumber\
z =&C_zx+D_ww+D_uu
\end{align}
where
Proof: It can easily be shown that the
If
- System () is strictly dissipative with the supply rate: \begin{equation} s(w,z)=w^\text{T}w \end{equation}
- For all
$\omega\in\mathbb{R}$ : \begin{equation} |T(j\omega)|{2,\infty}=\sup{0 < |w|2 < \infty}\frac{|z|\infty}{|w|_2} < \gamma \end{equation} - There exists
$P$ satisfying the LMIs: \begin{equation} \left[\begin{matrix} A^\text{T} P+PA&PB\\B^\text{T} P& -I\end{matrix}\right] \lt 0\nonumber\ \left[\begin{matrix} P&C^\text{T}\\C&\gamma^2 I\end{matrix}\right] \gt 0 \end{equation}
Consider the following linear system:
\begin{align}
\dot x=&Ax+B_ww+B_uu\nonumber\
z =&C_zx+D_uu
\end{align}
where
- Boyd, Stephen P. and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004.
- Scherer, Carsten and Weiland, Siep. Linear Matrix Inequalities in Control, Delft Center for Systems and Control, 2005.
- Boyd, Stephen P. Linear matrix inequalities in system and control theory. Vol. 15. Siam, 1994.
- Scherer, Carsten, Pascal Gahinet, and Mahmoud Chilali. "Multiobjective output-feedback control via LMI optimization." Automatic Control, IEEE Transactions on 42.7 (1997): 896-911.