What’s in this article?
The Sum of Squares (SOS) programming is an interesting optimization problem. It was introduced by Pablo A. Parrilo in his Ph.D. thesis. There are many magical applications of the SOS programming, including global optimization and proving a dynamic system is stable.
The SOS programming can be formidable from literature. It’s actually easy understand if we build the knowledge from the beginning.
I want to use simple language and examples to explain the SOS programming and implement some interesting applications of the SOS programming.
This article has,
- A simple explanation of the Sum of Squares (SOS) Programming.
- Using the SOS to solve a polynomial global optimization problem.
- Using the SOS to find a Lyapunov function for a control system.
- Using the SOS to compute Funnels (region of attractions) for a controller.
- SOS can be solved by Semidefinite Programming. I will discuss Semidefinite Programming.
Note: SOS is not Least Squares.
1. What is the Sum of Squares programming?
The SOS programming asks a simple question,
Can a polynomial function f be factored as Sum of Squares of polynomial basis b?
Examples of Sum of Squares.
\(b1 = (1, x, x^2, x^3)\). sum of squares of b1 is \(||b1||^2 = 1 + x^2 + x^4 + x^6\).
\(b2 = (x, 3y, 4xy)\). sum of squares of b2 is \(||b2||^2 = x^2 + 9y^2 + 16x^2y^2\).
Not all polynomial can be expressed as sum of squares. \(f(x) = -x^2\) can’t.
If a polynomial f can be refactored into sum of squares of b, then the polynomial is always positive. i.e. checking f is SOS is a sufficient condition for f > 0.
Consider a example for the SOS program.
Given a polynomial function,
\(f(x) = x^2 + 8 x^4\)
And polynomial basis,
\(b(x) =\begin{bmatrix} x \\ x^2 \end{bmatrix}\)
Define a quadratic function of the basis,
\( b^T Q b\)
Where \(Q\) is the coefficients for the SOS polynomial. The \(b^T Q b\) is a polynomial function. For example,
\(b^T Q b = \begin{bmatrix} x \\ x^2 \end{bmatrix}^T
\begin{bmatrix} q1 & q2 \\ q3 & q4 \end{bmatrix}
\begin{bmatrix} x \\ x^2 \end{bmatrix}
\\ = q1 x^2 + (q2 + q3) x^3 + q4 x^4\)
Because \(q2, q3\) are the same coefficients for term \(xy\), we let \(q2 = q3\). As a result \(Q\) is a symmetric matrix.
Go back to the SOS problem:
Can \(f(x)\) be factored as Sum of Squares (SOS) of polynomial basis \(b(x)\)?
Translate into math,
\(f(x) = b^T Q b \stackrel{?}{\in} SOS \)
I abused math notation a little bit. The symbol means,
- We select a basis \(b(x)\), and use \(b^T Q b\) to represent \(f(x)\). Q is the optimization variable.
- \(f(x) = b^T Q b\) are constraints for Q. Because \(b^T Q b\) and \(f(x)\) are polynomials, the constrain is affine for Q. It will be clear later in a example.
- Because \(b^T Q b = f(x)\), checking if \(b^T Q b\) is SOS is a sufficient condition (since we pre-select a basis) for \(f(x)\) is SOS.
Apply this equation to the example,
\(f(x) = b^T Q b\)
Where \(f(x) = x^2 + 8 x^4\), \(b(x) = \begin{bmatrix} x \\ x^2 \end{bmatrix}\), \(Q = \begin{bmatrix} q1 & q2 \\ q2 & q4 \end{bmatrix}\).
\(x^2 + 8 x^4 = q1 x^2 + (q2 + q2) x^3 + q4 x^4 \stackrel{?}{\in} SOS\)
The problem become: Can we find a symmetric matrix \(Q\) such that \(f(x) = b^T Q b \in SOS\)?
For this problem, I can easily see that \(q1 = 1, q2 = q3 = 0, q4 = 8\) makes \(f(x) = b^T Q b = ||(x^1, \sqrt{8} x^2)||^2 \in SOS\).
But Eye-balling a solution is not salable.
A better solution is to formulate “Can we find a symmetric matrix \(Q\) such that \(f(x) = b^T Q b \in SOS\)” as a Semidefinite program.
2. Semidefinite Programming
Semidefinite programming (SDP) is a convex optimization problem. It’s defined as
\(\text{minimize}_{X \in S^n} \; \text{trace}(CX)\)
\(\text{subject to:}
\\ \; \; \; \; \; X \succeq 0
\\ \; \; \; \; \; \text{trace}(A_i X) = b_i ,\; \forall i\)
I know, it’s confusing. Let’s me explain.
\(X \in S^n\) means \(X\) is a symmetry matrix.
\(X \succeq 0 \) means the matrix variable \(X\) needs to be positive-semidefinite.
Trace is the sum of diagonal elements in a square matrix. It is the linear operator for matrices.
For example, If we want to make the (0,0) element of \(\begin{bmatrix} x & y \\ z & k \end{bmatrix}\) equal to b.
We can do,
\(\text{trace}(\begin{bmatrix} 1 & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} x & y \\ z & k \end{bmatrix}) = b\)
It implies “I want a element of X to be b” can be expressed in the form of \(\text{trace}(A_i X) = b_i \), which is the SDP constraints.
SDP programs are be solved by Newton’s Method. In the end of this article, I will give a brief overview of how to solve a SDP by Primal-Dual Interior Point Methods.
There are free SDP solvers. In the implementation section, I will use the CVXPY solver for SDP programs.
3. Convert SOS to SDP
Go back to the example SOS problem, we want to check,
\(f(x) = b^T Q b \stackrel{?}{\in} SOS \)
\(x^2 + 8 x^4 = q1 x^2 + (q2 + q2) x^3 + q4 x^4 \stackrel{?}{\in} SOS\)
i.e. Can we find a symmetric matrix Q such that \(b^T Q b \in SOS\) with the constraint \(x^2 + 8 x^4 = q1 x^2 + (q2 + q2) x^3 + q4 x^4\)?
A sufficient condition for \(b^T Q b\) is sum of squares is \(Q \succeq 0\). i.e. Q is a positive-semidefinite matrix.
Assuming there is a \(Q \succeq 0\), then \(Q\) can be factored as \(L^T L \). (Assuming Q is not singular, it is a well-know property of PSD matrices. Check this link).
\(b^T Q b = b^T L^T L b = ||L b||^2 \in SOS\)
b is a vector of polynomial basis. \(L b\) is a linear transformation of b. If have found a \(Q \succeq 0\), then \(b^T Q b\) is sum of squares.
This idea can be formulated as a SDP.
We want to find a \(Q \succeq 0\) subject to the coefficient of the polynomial function \(b^T Q b\) is the same as the coefficients of the polynomial function \(f(x)\).
In math,
\(\text{minimize}_{Q \in S^n} \; 1\)
\(\text{subject to:}
\\ \; \; \; \; \; Q \succeq 0
\\ \; \; \; \; \; b^T Q b = f(x) \)
We can convert it to the standard SDP form by matching the coefficient for \(b^T Q b = f(x)\)
The objective function minimize 1 mean it’s a feasibility SDP problem. e.g. We want to found a \(Q\) that satisfies all constraints.
By solving this SDP, we know if \(Q \succeq 0\) exist. If we have a \(Q \succeq 0\). we find a symmetric matrix \(Q\) such that \(f(x) = b^T Q b \in SOS\)
Go back to the example,
\(x^2 + 8 x^4 = q1 x^2 + 2 q2 x^3 + q4 x^4 \stackrel{?}{\in} SOS\)
By matching coefficients \((q1, q2, q4)\) of the \(b^T Q b\) with the coefficient of \( f(x) = x^2 + 8 x^4\),
\(\text{minimize}_{Q \in S^n} \; 1\)
\(\text{subject to:}
\\ \; \; \; \; \; Q \succeq 0
\\ \; \; \; \; \; q1 = 1
\\ \; \; \; \; \; q2 = 0
\\ \; \; \; \; \; q4 = 8 \)
It’s a SDP problem, we can use a SDP solver to check if there exist a Q that satisfy the problem.
Now, we know what SOS program does and how to solve a SOS program by SDP. Let’s discuss some magical applications of SOS programming.
4. SOS application: Polynomial Function Global Optimization
Optimize a non-convex function is a NP hard-problem. However for a non-convex polynomial function, we can do the global optimization efficiently.
4.1 Problem Setup
I will explain what happened in this slide from Pablo A. Parrilo, who is the godfather for SOS programming.
Given a non-convex polynomial function,
\(F(x,y) = 4x^2 – \frac{21}{10}x^4 + \frac{1}{3}x^6 + xy – 4y^2 + 4y^4\)
We want to optimization \(F(x,y)\) globally.
Obviously gradient based methods find a local minimum. The best we can do to find global minimum is to try many initial guesses and select the best local minimum. The procedure is NP hard.
Alternatively, we can formulate this problem as a SOS program.
We can define \(\gamma\) to be the lower bound of \(F(x,y)\). Then we check if \(F(x, y) – \gamma > 0\) by checking a sufficient condition \(F(x, y) – \gamma \in SOS\) with a predefined polynomial basis \(b(x)\).
One way to minimum the lower bound \(\gamma\) is doing a binary search for highest \(\gamma^*\) which makes \(F(x, y) – \gamma \in SOS\).
But we can easily prove that,
\(\text{minimize} \; \gamma\)
\(\text{subject to:}
\\F(x,y) – \gamma \in \text{SOS from b(x)} \)
Is a SOS program.
So, we can solve for \(\gamma^*\) directly in a SOS program.
Why we can optimize a non-convex function globally? My intuition is we convert the optimization in \((x, y)\) space into an optimization in the space of polynomial coefficients \(Q\). In the \(Q\) space, the optimization problem is convex. (with the caveat that we solved a sufficient condition of the original problem.)
4.2 Non-convex Optimization Example and Implementation
Let’s solve the problem given in the slide.
Given \(F(x,y) = 4x^2 – \frac{21}{10}x^4 + \frac{1}{3}x^6 + xy – 4y^2 + 4y^4\), the goal is,
\(\text{minimize}_{x,y} F(x, y)\)
As we discuss above, we are going to solve for the sufficient condition of the lower-bound of \(F(x, y)\).
\(\text{minimize}_{Q \in S^n} \; \gamma\)
\(\text{subject to:}
\\F(x,y) – \gamma = b(x,y)^T Q b(x, y) \in \text{SOS} \)
I choose the polynomial basis as \(b(x, y) = (1, x, x^2, x^3, y, y^2, y^3, xy, x^2y, xy^2)\)
I converted the above SOS program into a SDP program. And I used the CVXPY to solve the SDP program.
Symbolic Coefficient Matching
In the SOS program, we let \(b^T Q b = F(x, y) – \gamma\).
To convert \(b^T Q b = F(x, y) – \gamma\) to the form of SDP constraints, e.g. \(Q(i, j) = b_i \;, \forall i, j\) and \(Q \succeq 0\). We need to match the coefficient of polynomial \(b^T Q b\) with the coefficient of polynomial \(F(x, y) – \gamma\). (Go back to the simple example if you get confused.)
Doing it by hand is way too hard. I did the coefficient matching by symbolic computation.
def coefficient_symbolic_match(self):
x, y, gamma = sympy.symbols('x y gamma')
# f(x, y) = 4 x^2 - 21/10* x^4 + 1/3 x^6 + xy - 4y^2 + 4y^4
f_monomials = [x**2, x**4, x**6, x*y, y**2, y**4]
f_coeffs = [4., -21/10., 1/3., 1., -4., 4.]
# b^T Q b
w = sympy.Matrix([1, x, x**2, x**3, y, y**2, y**3, x*y, x*y*y, x*x*y])
Q = sympy.MatrixSymbol('Q', 10, 10)
V_dot_SOS = (w.T @ Q @ w).as_explicit()
V_dot_SOS_poly = sympy.Poly(V_dot_SOS[0], x, y)
print('V_dot_SOS_poly:', V_dot_SOS_poly)
constraint_list_poly = []
for f_monomial, f_coeff in zip(f_monomials, f_coeffs):
Q_coeff = V_dot_SOS_poly.coeff_monomial(f_monomial)
constrain = '{}=={}'.format(Q_coeff, f_coeff)
print('constrain:', constrain)
constraint_list_poly.append(constrain)
MAX_ORDER = 10
constraint_list_zero = []
for x_order in range(0, MAX_ORDER + 1):
for y_order in range(0, MAX_ORDER + 1):
# skip symmetry. not sure how to do it.
# having duplicate constraints seem ok :)
# skip constant, gamma will do it
if y_order == 0 and x_order == 0:
continue
monomial = x**x_order * y ** y_order
# skip non-zero coef
if monomial in f_monomials:
continue
coeff = V_dot_SOS_poly.coeff_monomial(monomial)
if not coeff is sympy.S.Zero:
constrain = '{} == 0'.format(coeff)
print('constrain:', constrain, 'for coef:',
x**x_order * y ** y_order)
constraint_list_zero.append(constrain)
print('constraint_poly:', ','.join(constraint_list_poly))
print('constraint_zero:', ','.join(constraint_list_zero))
return constraint_list_poly, constraint_list_zeroThe code matches polynomial coefficients of \(b^T Q b\) with the coefficients of \(F(x, y) – \gamma\).
For instance, we got \(b^T Q b = … + (Q(1,2) + Q(0,0)) x^2 y\) and \(F(x,y) = 1024 x^2 y\). Then we got a constrain \(Q(1,2) + Q(0, 0) == 1024\).
In the end, I outputted the constrains. I copied them to the CVXPY solver manually. I can do a compiling script to do it automatically. But I’d prefer to keep thing simple.
Solving the SDP Program
def solve_sos_as_sdp(self):
num_var_w = 10
Q = cp.Variable((num_var_w, num_var_w), symmetric=True)
gamma = cp.Variable()
# sufficient condition
Epsilon = 0
constraints = [Q >> Epsilon * np.identity(num_var_w)]
constraints += [Q[0, 0] == -gamma]
constraints += [Q[0, 2] + Q[1, 1] + Q[2, 0] == 4.0, Q[1, 3] + Q[2, 2] + Q[3, 1] == -2.1, Q[3, 3] == 0.3333333333333333,
Q[0, 7] + Q[1, 4] + Q[4, 1] + Q[7, 0] == 1.0, Q[0, 5] + Q[4, 4] + Q[5, 0] == -4.0, Q[4, 6] + Q[5, 5] + Q[6, 4] == 4.0]
constraints += [Q[0, 4] + Q[4, 0] == 0, Q[0, 6] + Q[4, 5] + Q[5, 4] + Q[6, 0] == 0, Q[5, 6] + Q[6, 5] == 0, Q[6, 6] == 0, Q[0, 1] + Q[1, 0] == 0, Q[0, 8] + Q[1, 5] + Q[4, 7] + Q[5, 1] + Q[7, 4] + Q[8, 0] == 0, Q[1, 6] + Q[4, 8] + Q[5, 7] + Q[6, 1] + Q[7, 5] + Q[8, 4] == 0, Q[5, 8] + Q[6, 7] + Q[7, 6] + Q[8, 5] == 0, Q[6, 8] + Q[8, 6] == 0, Q[0, 9] + Q[1, 7] + Q[2, 4] + Q[4, 2] + Q[7, 1] + Q[9, 0] == 0, Q[1, 8] + Q[2, 5] + Q[4, 9] + Q[5, 2] + Q[7, 7] + Q[8,
1] + Q[9, 4] == 0, Q[2, 6] + Q[5, 9] + Q[6, 2] + Q[7, 8] + Q[8, 7] + Q[9, 5] == 0, Q[6, 9] + Q[8, 8] + Q[9, 6] == 0, Q[0, 3] + Q[1, 2] + Q[2, 1] + Q[3, 0] == 0, Q[1, 9] + Q[2, 7] + Q[3, 4] + Q[4, 3] + Q[7, 2] + Q[9, 1] == 0, Q[2, 8] + Q[3, 5] + Q[5, 3] + Q[7, 9] + Q[8, 2] + Q[9, 7] == 0, Q[3, 6] + Q[6, 3] + Q[8, 9] + Q[9, 8] == 0, Q[2, 9] + Q[3, 7] + Q[7, 3] + Q[9, 2] == 0, Q[3, 8] + Q[8, 3] + Q[9, 9] == 0, Q[2, 3] + Q[3, 2] == 0, Q[3, 9] + Q[9, 3] == 0]
prob = cp.Problem(cp.Minimize(-gamma),
constraints)
prob.solve(verbose=False)
# Print result.
print("status:", prob.status)
print("The optimal value is", prob.value)
print("The low bound is", gamma.value)With the symbolic computation, I just need to copy the symbolic constraints to the SDP solver and solve the problem.
The output
----------------------------------------------------------------------------
SCS v2.1.1 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
----------------------------------------------------------------------------
Lin-sys: sparse-direct, nnz in A = 111
eps = 1.00e-04, alpha = 1.50, max_iters = 5000, normalize = 1, scale = 1.00
acceleration_lookback = 10, rho_x = 1.00e-03
Variables n = 56, constraints m = 83
Cones: primal zero / dual free vars: 28
sd vars: 55, sd blks: 1
Setup time: 2.55e-04s
----------------------------------------------------------------------------
Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
0| 3.64e+19 3.17e+19 1.00e+00 -3.29e+20 1.51e+20 2.23e+20 1.64e-04
80| 7.04e-05 5.44e-05 7.87e-06 1.03e+00 1.03e+00 4.66e-16 6.04e-03
----------------------------------------------------------------------------
Status: Solved
Timing: Solve time: 6.06e-03s
Lin-sys: nnz in L factor: 250, avg solve time: 8.65e-07s
Cones: avg projection time: 2.31e-05s
Acceleration: avg step time: 4.65e-05s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 1.1876e-09, dist(y, K*) = 2.3932e-09, s'y/|s||y| = -1.2132e-12
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 7.0445e-05
dual res: |A'y + c|_2 / (1 + |c|_2) = 5.4439e-05
rel gap: |c'x + b'y| / (1 + |c'x| + |b'y|) = 7.8653e-06
----------------------------------------------------------------------------
c'x = 1.0316, -b'y = 1.0317
============================================================================
status: optimal
The optimal value is 1.0316461331878024
The low bound is -1.0316461331878024The lower-bound \(\gamma \) of \(F(x,y)\) is computed as -1.0316. Which is the same as the low-bound -1.0316 from the slide, and (luckly) it’s also the global minimum.
5. Prove the Stability of a Nonlinear Control System by Finding a Lyapunov Function
Given a dynamical system (a different equation), if we can find a Lyapunov function for it, the system is stable (converge to a fixed state).
I know, it’s confusing. Let me discuss some backgrounds needed.
5.1 Stability
Given a differential equations, \(\dot{x} = f(x)\). We want to know whether the differential equation is bounded (stable).
For example, consider a discrete case \(\Delta x_i = 2 x_i\). We define the propagation as \(x_{i+1} = \Delta x_i + x_i\)
Consider propagating the system in time.
\(x_0 = 1. \\ x_1 = \Delta x_0 + x_0 = 2 x_0 + x_0 = 3 x_0 \\ x_2 = \Delta x_1 + x_1 = 2 x_1 + x_1 = 9 x_0 \\ …. \\ x_n = 3^n x_0\)
As \(n \Rightarrow \inf\), \(x_n \Rightarrow \inf\). The system is not stable.
Another example, \(\Delta y_i = -0.1 y_i\). The propagation is defined as \(y_{i+1} = \Delta y_i + y_i\).
We can propagate this system,
\(y_0 = 1. \\ y_1 = y_0 – 0.1 y_0 = 0.9 y_0 \\ y_2 = y_1 – 0.1 y_1 = 0.9 y_1 = 0.81 y_0 \\ …. \\ y_n = 0.9^n y_0\)
As \(n \Rightarrow \inf\), \(y_n \Rightarrow 0\). The system is stable.
The same analysis also applied for complex dynamical control system.
For example, given a complex system (e.g. rocket landing)
\(\dot x = f(x, u)\)
Where, \(x\) is the state. \(\dot x\) is the time derivative of the state. \(u\) is the control.
For control engineers, they want to design a control law \(u = \{u_0, u_1, u_2, u_3, … \}\) such that, when time goes, the state \(x_n = x_{target}\).
There are many classical methods to design the control law. And we always want to make sure that the dynamical system under the control law is stable.
How to prove the stability of a control system?
The existence of a Lyapunov function proves it.
5.2 Lyapunov Function
Given a differential equation \(\dot{x} = f(x)\),
If we can find a non-negative scalar function V(x) such at
\(V(x) > 0 \; , \forall x\).
And the V(x) decrease along the dynamic system \(\dot{x} = f(x)\).
\(\frac{\partial V(x)}{\partial t} = \frac{\partial V(x)}{\partial x} \frac{\partial x}{\partial t} = \frac{\partial V(x)}{\partial x} \dot{x} < 0\)
Then the differential equation \(\dot{x} = f(x)\) is stable.
In English, the existence of \(V(x)\) certified that \(x\) converging to a lower \(V(f(x))\) along the time coordination. Because \(V(x) > 0\) on the domain x, we can argue that \(x_{n->\inf} \) converge to a local minimum of \(V(x)\). It means \(f(x)\) is stable.
\(V(x)\) is the Lyapunov function for the dynamical system \(f(x)\).
In sum, the existence of Lyapunov Function prove the stability of a dynamic system.
5.3 Finding a Lyapunov Function for a Polynomial Nonlinear Dynamic System
For linear system, a Lyapunov function can be found by solving the linear Lyapunov equation.
For a Polynomial dynamical system, we can find a Lyapunov function by SOS programming.
I will explain what happened in this slide from Pablo A. Parrilo.
Given a Jet engine dynamical system,
\(\dot{x} = -y + \frac{3}{2}x^2 – \frac{1}{2} x^3 \\ \dot{y} = 3x – y\)
We want to prove the stability of the system by finding a Lyapunov function \(V(x)\).
We want,
\(\text{minimize}_{Q \in S^n} \; 1\)
\(\text{subject to:}
\\ V(x, y, Q) \leq 0
\\ – \nabla V(x, y, Q) f(x, y) \leq 0\)
Since the dynamical system is polynomial, we can do it by SOS programming, which is a sufficient condition for the original problem.
\(\text{minimize}_{Q \in S^n} \; 1\)
\(\text{subject to:}
\\ V(x, y, Q) \in -SOS \leq 0
\\ – \nabla V(x, y, Q) f(x, y) \in -SOS \leq 0\)
Concretely, we need to
- Choose a polynomial basis \(b1(x, y)\).
- The optimization variables are the symmetric matrix \(Q\).
- Express \(V(x, y)\) by \(b1(x, y)^T Q b1(x,y)\). It can be modeled as affine constraints for the optimization variable Q. \(V(x, y, Q) = b1^T Q b1\).
- Enforcing this first Lyapunov condition \(V(x, y, Q) = b1^T Q b1 > 0\) by enforcing a sufficient condition \(b1^T Q b1 \in SOS\). The SOS can be directly convert to SDP \(Q > 0\).
- Enforce the second Lyapunov condition. \(– \nabla V(x, y, Q) f(x, y) > 0\) by enforcing a sufficient condition \(– \nabla V(x, y, Q) f(x, y) \in SOS\).
- However, because the optimization variable is \(Q\) we need to convert the constraints to the form of an affine function of \(Q\). Firstly, we need to choose a polynomial basis b2(x,y) for \(– \nabla V(x, y, Q) f(x, y) \in SOS\). And we do the math to convert \(– \nabla V(x, y, Q) f(x, y) \in SOS\) to \(b2^T \text{affine-function(Q)} b2 \in SOS\) by matching coefficients. Then, we can convert the SOS to SDP constraints : \(\text{affine-function(Q)} > 0\).
Now, we have,
\(\text{minimize}_{Q \in S^n} \; 1\)
\(\text{subject to:}
\\ V(x, y, Q) = b1(x,y)^T Q b1(x, y) \in \text{SOS}
\\ – \nabla V(x, y, Q) f(x, y) = b2(x, y)^T \text{affine-function}(Q) b2(x, y) \in \text{SOS}\)
Because all constraints are affine in Q, we can convert this SOS program into SDP.
\(\text{minimize}_{Q \in S^n} \; 1\)
\(\text{subject to:}
\\ Q \succeq 0
\\ \text{affine-function}(Q) \succeq 0\)
The previous derivation can be confusing without a example.
Example 1. A analytical example to find a Lyapunov function.
Example 2. Algebra doesn’t scale. We should convert the problem to SDP and let compute do it for us.
5.4 Implementation
I used CVXPY to do the SDP programming, and the sympy to match coefficients for polynomials.
# Import packages.
import cvxpy as cp
import numpy as np
import sympy
from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import
import matplotlib.pyplot as plt
from matplotlib import cm
class LyapunovSOS:
def __init__(self):
# Warning: you need to copy constaints!
x, y = sympy.symbols('x y')
self.f = sympy.Matrix([
[-y - 3/2*x**2 - 1/2*x**3],
[3*x - y],
])
self.b1 = sympy.Matrix(
[[x, y, x**2, x*y, y**2, x**3, x**2*y, x*y**2, y**3]]).transpose()
def polynomial_arrangement(self):
x, y = sympy.symbols('x y')
Q = sympy.MatrixSymbol('Q', 9, 9)
dzdx = sympy.diff(self.b1, x)
dzdy = sympy.diff(self.b1, y)
B = sympy.Matrix([
[dzdx, dzdy]
])
V = (self.b1.T @ Q @ self.b1).as_explicit()
V_poly = sympy.Poly(V[0], x, y)
V_dot = (- 2 * self.b1.T @ Q @ B @ self.f).as_explicit()
V_dot_poly = sympy.Poly(V_dot[0], x, y)
b2 = sympy.Matrix(
[[x, y, x**2, x*y, y**2, x**3, x**2*y, x*y**2, y**3]]).transpose()
G = sympy.MatrixSymbol('G', 9, 9)
SOS_of_V_dot = (b2.T @ G @ b2).as_explicit()
SOS_of_V_dot_poly = sympy.Poly(SOS_of_V_dot[0], x, y)
# Very tricky!
b2_sqr = [b*b for b in b2]
for b_sqr in b2_sqr:
coeff = V_dot_poly.coeff_monomial(b_sqr)
if coeff is sympy.S.Zero:
print('WARNING: digonal coeff is zero for',
b_sqr, ' not PSD for sure!')
constraint_list = []
for max_order in range(10):
for x_order in range(0, max_order + 1):
y_order = max_order - x_order
monomial = x ** x_order * y ** y_order
SOS_coeff = SOS_of_V_dot_poly.coeff_monomial(
monomial)
if SOS_coeff is sympy.S.Zero:
continue
V_dot_coeff = V_dot_poly.coeff_monomial(monomial)
if V_dot_coeff is not sympy.S.Zero:
constrain = '{}=={}'.format(V_dot_coeff, SOS_coeff)
print('constrain:', constrain, " of ", monomial)
constraint_list.append(constrain)
else:
constrain = '{}==0.'.format(SOS_coeff)
print('constrain:', constrain, " of ", monomial)
constraint_list.append(constrain)
print('Constraints (copy this!):', ','.join(constraint_list))
def solve_sos_as_sdp(self):
num_var_q = 9
Q = cp.Variable((num_var_q, num_var_q), symmetric=True)
num_var_w = 9
G = cp.Variable((num_var_w, num_var_w), symmetric=True)
# sufficient condition
Epsilon = 1e-8
constraints = [Q >> Epsilon * np.identity(num_var_q)]
constraints += [G >> Epsilon * np.identity(num_var_w)]
constraints += [2*Q[1, 0] + 2*Q[1, 1] == G[1, 1], 2*Q[0, 0] + 2*Q[0, 1] - 6*Q[1, 1] == G[0, 1] + G[1, 0], -6*Q[0, 1] == G[0, 0], 2*Q[1, 3] + 4*Q[1, 4] + 2*Q[4, 0] + 2*Q[4, 1] == G[1, 4] + G[4, 1], 2*Q[0, 3] + 4*Q[0, 4] + 4*Q[1, 2] + 2*Q[1, 3] - 12*Q[1, 4] + 2*Q[3, 0] + 2*Q[3, 1] - 6*Q[4, 1] == G[0, 4] + G[1, 3] + G[3, 1] + G[4, 0], 4*Q[0, 2] + 2*Q[0, 3] - 12*Q[0, 4] + 3.0*Q[1, 0] - 6*Q[1, 3] + 2*Q[2, 0] + 2*Q[2, 1] - 6*Q[3, 1] == G[0, 3] + G[1, 2] + G[2, 1] + G[3, 0], 3.0*Q[0, 0] - 6*Q[0, 3] - 6*Q[2, 1] == G[0, 2] + G[2, 0], 2*Q[1, 7] + 6*Q[1, 8] + 2*Q[4, 3] + 4*Q[4, 4] + 2*Q[8, 0] + 2*Q[8, 1] == G[1, 8] + G[4, 4] + G[8, 1], 2*Q[0, 7] + 6*Q[0, 8] + 4*Q[1, 6] + 4*Q[1, 7] - 18*Q[1, 8] + 2*Q[3, 3] + 4*Q[3, 4] + 4*Q[4, 2] + 2*Q[4, 3] - 12*Q[4, 4] + 2*Q[7, 0] + 2*Q[7, 1] - 6*Q[8, 1] == G[0, 8] + G[1, 7] + G[3, 4] + G[4, 3] + G[7, 1] + G[8, 0], 4*Q[0, 6] + 4*Q[0, 7] - 18*Q[0, 8] + 3.0*Q[1, 3] + 6*Q[1, 5] + 2*Q[1, 6] - 12*Q[1, 7] + 2*Q[2, 3] + 4*Q[2, 4] + 4*Q[3, 2] + 2*Q[3, 3] - 12*Q[3, 4] + 3.0*Q[4, 0] - 6*Q[4, 3] + 2*Q[6, 0] + 2*Q[6, 1] - 6*Q[7, 1] == G[0, 7] + G[1, 6] + G[2, 4] + G[3, 3] + G[4, 2] + G[6, 1] + G[7, 0], 3.0*Q[0, 3] + 6*Q[0, 5] + 2*Q[0, 6] - 12*Q[0, 7] + 1.0*Q[1, 0] + 6.0*Q[1, 2] - 6*Q[1, 6] + 4*Q[2, 2] + 2*Q[2, 3] - 12*Q[2, 4] + 3.0*Q[3, 0] - 6*Q[3, 3] + 2*Q[5, 0] + 2*Q[5, 1] - 6*Q[6, 1] == G[0, 6] + G[1, 5] + G[2, 3] + G[3, 2] + G[5, 1] + G[6, 0], 1.0*Q[0, 0] + 6.0*Q[0, 2] - 6*Q[0, 6] + 3.0*Q[2, 0] - 6*Q[2, 3] - 6*Q[5, 1] == G[0, 5] + G[2, 2] + G[5, 0], 2*Q[4, 7] + 6*Q[4, 8] + 2*Q[8, 3] + 4*Q[8, 4] == G[4, 8] + G[8, 4], 2*Q[3, 7] + 6*Q[3, 8] + 4*Q[4, 6] + 4*Q[4, 7] - 18*Q[4, 8] + 2*Q[7, 3] + 4*Q[7, 4] + 4*Q[8, 2] + 2*Q[8, 3] - 12*Q[8, 4] == G[3, 8] + G[4, 7] + G[7, 4] + G[8, 3], 3.0*Q[1, 7] + 2*Q[2, 7] + 6*Q[2, 8] + 4*Q[3, 6] + 4*Q[3, 7] - 18*Q[3, 8] + 3.0*Q[4, 3] + 6*Q[4, 5] + 2*Q[4, 6] - 12*Q[4, 7] + 2*Q[6, 3] + 4*Q[6, 4] +
4*Q[7, 2] + 2*Q[7, 3] - 12*Q[7, 4] + 3.0*Q[8, 0] - 6*Q[8, 3] == G[2, 8] + G[3, 7] + G[4, 6] + G[6, 4] + G[7, 3] + G[8, 2], 3.0*Q[0, 7] + 1.0*Q[1, 3] + 6.0*Q[1, 6] + 4*Q[2, 6] + 4*Q[2, 7] - 18*Q[2, 8] + 3.0*Q[3, 3] + 6*Q[3, 5] + 2*Q[3, 6] - 12*Q[3, 7] + 1.0*Q[4, 0] + 6.0*Q[4, 2] - 6*Q[4, 6] + 2*Q[5, 3] + 4*Q[5, 4] + 4*Q[6, 2] + 2*Q[6, 3] - 12*Q[6, 4] + 3.0*Q[7, 0] - 6*Q[7, 3] == G[2, 7] + G[3, 6] + G[4, 5] + G[5, 4] + G[6, 3] + G[7, 2], 1.0*Q[0, 3] + 6.0*Q[0, 6] + 2.0*Q[1, 2] + 9.0*Q[1, 5] + 3.0*Q[2, 3] + 6*Q[2, 5] + 2*Q[2, 6] - 12*Q[2, 7] + 1.0*Q[3, 0] + 6.0*Q[3, 2] - 6*Q[3, 6] + 4*Q[5, 2] + 2*Q[5, 3] - 12*Q[5, 4] + 3.0*Q[6, 0] - 6*Q[6, 3] == G[2, 6] + G[3, 5] + G[5, 3] + G[6, 2], 2.0*Q[0, 2] + 9.0*Q[0, 5] + 1.0*Q[2, 0] + 6.0*Q[2, 2] - 6*Q[2, 6] + 3.0*Q[5, 0] - 6*Q[5, 3] == G[2, 5] + G[5, 2], 2*Q[8, 7] + 6*Q[8, 8] == G[8, 8], 2*Q[7, 7] + 6*Q[7, 8] + 4*Q[8, 6] + 4*Q[8, 7] - 18*Q[8, 8] == G[7, 8] + G[8, 7], 3.0*Q[4, 7] + 2*Q[6, 7] + 6*Q[6, 8] + 4*Q[7, 6] + 4*Q[7, 7] - 18*Q[7, 8] + 3.0*Q[8, 3] + 6*Q[8, 5] + 2*Q[8, 6] - 12*Q[8, 7] == G[6, 8] + G[7, 7] + G[8, 6], 1.0*Q[1, 7] + 3.0*Q[3, 7] + 1.0*Q[4, 3] + 6.0*Q[4, 6] + 2*Q[5, 7] + 6*Q[5, 8] + 4*Q[6, 6] + 4*Q[6, 7] - 18*Q[6, 8] + 3.0*Q[7, 3] + 6*Q[7, 5] + 2*Q[7, 6] - 12*Q[7, 7] + 1.0*Q[8, 0] + 6.0*Q[8, 2] - 6*Q[8, 6] == G[5, 8] + G[6, 7] + G[7, 6] + G[8, 5], 1.0*Q[0, 7] + 2.0*Q[1, 6] + 3.0*Q[2, 7] + 1.0*Q[3, 3] + 6.0*Q[3, 6] + 2.0*Q[4, 2] + 9.0*Q[4, 5] + 4*Q[5, 6] + 4*Q[5, 7] - 18*Q[5, 8] + 3.0*Q[6, 3] + 6*Q[6, 5] + 2*Q[6, 6] - 12*Q[6, 7] + 1.0*Q[7, 0] + 6.0*Q[7, 2] - 6*Q[7, 6] == G[5, 7] + G[6, 6] + G[7, 5], 2.0*Q[0, 6] + 3.0*Q[1, 5] + 1.0*Q[2, 3] + 6.0*Q[2, 6] + 2.0*Q[3, 2] + 9.0*Q[3, 5] + 3.0*Q[5, 3] + 6*Q[5, 5] + 2*Q[5, 6] - 12*Q[5, 7] + 1.0*Q[6, 0] + 6.0*Q[6, 2] - 6*Q[6, 6] == G[5, 6] + G[6, 5], 3.0*Q[0, 5] + 2.0*Q[2, 2] + 9.0*Q[2, 5] + 1.0*Q[5, 0] + 6.0*Q[5, 2] - 6*Q[5, 6] == G[5, 5]]
# It is even possible to pick a particularly elegant solution, given by a quadratic Lyapunov function.
# This can be achieved by minimizing the sum of diagonal elements corresponding to
# the nonquadratic terms, subject to the LMI constraints.
C = np.identity(num_var_q)
C[0:2, 0:2] = 0.
prob = cp.Problem(cp.Minimize(cp.trace(C @ Q)),
constraints)
prob.solve(verbose=False)
# Print result.
print("status:", prob.status)
print("The optimal value is", prob.value)
print("A solution Q is")
print(Q.value)
V_optimal = self.b1.T @ Q.value @ self.b1
x, y = sympy.symbols('x y')
V_func = sympy.utilities.lambdify([x, y], V_optimal)
self.plot(V_func)
def plot(self, V_func):
...
def main():
lyapunov = LyapunovSOS()
lyapunov.polynomial_arrangement()
lyapunov.solve_sos_as_sdp()
if __name__ == "__main__":
main()
Result
Pablo A. Parrilo’s result,
6. The Generalized S-procedure for SOS programming
We talked about how to prove the global stability by finding a global Lyapunov function. In reality, most control systems are not globally stable (otherwise, all control engineers will lost their job).
We are interested in proving the stability of a dynamical system in a region close to the current state.
In the context of the SOS programming, it’s done by the generalized S-procedure.
6.1 Polynomial Lagrangian
Given \(p(x)\) and \(g(x)\), we want to show \(p(x) > 0\) on the region defined by \(\{x | g(x) > 0 \}\).
To prove \(p(x) > 0 \text{ on } \{x | g(x) > 0 \}\), we can do a Lagrangian multiplier.
\(L(x) = p(x) – l(x)g(x)\)
where, \(l(x)\) is a polynomial Lagrangian multiplier.
Assuming we can show,
- \(L(x) > 0 \; , \forall x\)
- \(l(x) > 0 \; , \forall x\)
Then, we can argue on \(\{x | g(x) > 0 \}\), \(-l(x)g(x) < 0\). Then,
\(L(x) = p(x) – l(x)g(x) > 0\)
\(L(x) = p(x) > l(x)g(x)\)
Since \(l(x)g(x) > 0\) on \(\{x | g(x) > 0 \}\),
\(L(x) = p(x) > 0, \; \text{ on } \{x | g(x) > 0 \}\)
We showed if we can prove \(L(x) > 0 \; ,\forall x\) and \(l(x) > 0 \; ,\forall x\), then \(p(x) > 0\) on the region defined by \(\{x | g(x) > 0 \}\).
However, showing \(L(x) > 0 \; ,\forall x\) and \(l(x) > 0 \; ,\forall x\) is hard. But a sufficient condition is,
- \(L(x) \in SOS \rightarrow L(x) > 0 \; \forall x\)
- \(l(x) \in SOS \rightarrow l(x) > 0 \; \forall x\)
Check SOS is computationally feasible. If we can prove the above statements, we can also show \(p(x) > 0\) on the region defined by \(\{x | g(x) > 0 \}\).
6.2 Example
With loss of generality, let me assume,
\(p(x) = (x-3)^2 – 1\)
\(g(x) = 1 – x^2\)
We want to show \(p(x) > 0\) on the region defined by \(\{x | g(x) > 0 \}\).
Select the polynomial Lagrange multiplier to be \(l(x) = 1^T q 1 = q\). (For how to choice the degrees of the polynomial, please read Pablo’s slides)
To show \(p(x) > 0\) on the region defined by \(\{x | g(x) > 0 \}\), we need convert the problem into 2 SOS programs.
- \(L(x) \in SOS \rightarrow L(x) > 0 \; \forall x\)
- \(l(x) \in SOS \rightarrow l(x) > 0 \; \forall x\)
To make things clear, I will convert the 2 SOS into a SDP.
Now, we just need to type the SDP program into CVXPY.
6.3 Implementation
# Import packages.
import cvxpy as cp
import numpy as np
# verify p(x) >= 0 on g(x) > 0
# p(x) = (x-3)^2 - 1
# g(x) = 1 - x^2
#
# L = p(x) - l(x)g(x)
# l(x) := 1^T Q 1, sum of squares of 1
def solve_simple_constraint_sos_1():
num_var = 1
Q = cp.Variable((num_var, num_var), symmetric=True)
slack = cp.Variable((2, 2), symmetric=True)
# Q.value = np.identity(num_var)
# sufficient condition
Epsilon = 1e-4 * np.identity(num_var)
constraints = [Q >> Epsilon]
constraints += [slack >> Epsilon]
constraints += [-Q + 8 == slack[0,0]]
constraints += [-3 == slack[1,0]]
constraints += [1+Q == slack[1,1]]
prob = cp.Problem(cp.Minimize(1),
constraints)
prob.solve(verbose = False)
# Print result.
print("status:", prob.status)
print("The optimal value is", prob.value)
print("A solution Q is")
print(Q.value)
...
def main():
solve_simple_constraint_sos_1()
...
if __name__ == "__main__":
main()yimu@yimu-mate:~/Desktop/blog$ python3 /blog/yimu-blog/random/lyapunov/sos_s_procedure_simple.py # problem 1 status: optimal The optimal value is 1.0 A solution Q is [[2.31821375]] # problem 2 status: infeasible The optimal value is inf A solution Q is None
In the Python program, the SDP found a feasible solution. To do a comparison, I also did a SDP for a infeasible problem.
6.4 Finding a local Lyapunov function in Action
For a real controller, people want to prove it is stable in a region. In this paper, authors used the SOS programming and Lyapunov analysis to show the stability of the controller of an aircraft in the falling-leaf mode.
Firstly, we need a couple of definitions,
Region of attraction,
Sub-level set parameterized by \(\gamma\),
The author wanted to find the largest sub-level set \(\gamma^*\) which make the controller stable \(\nabla V(x) f(x) = \dot V(x) < 0\). To make things simple for analysis, the author found the largest ellipsoid within the sub-level set.
For (7), the author didn’t mention how they solve the optimization problem. My guess is, the author wanted to find the largest \(\gamma\), such that \(\nabla V(x) f(x) < 0\) on \(V(x) < \gamma\). Apply the generalized S-procedure, we have a SOS program for (7),
\(\text{maximize } \gamma \)
\(\text{subject to:} \\ \nabla V(x) f(x) + l(x) (V(x) – \gamma) < 0 \in -SOS \\ l(x) > 0 \in SOS\)
For (8), given \(\gamma^*\), the author formulated the problem of “find the largest ellipsoid within the sub-level set” as a SOS.
7. Paper unlock, The Funnels
The LQR-Trees is a planning algorithm similar to the RRT.
In the RRT planning algorithm, we randomly grow a tree until a target state is reached. The search space can be huge.
In the LQR-Trees, we grow the tree by doing LQR and compute a region of attraction. The idea is within the region of attraction, we can reuse the existing LQR control law. So, we don’t need to sample states inside existing region of attractions. As a result, the search space is smaller.
The paper mainly discussed how to apply the SOS programming and use the generalized S-procedure to compute the region of attractions, which is defined as Funnels.
7.1 Definition of Funnels
The definition of funnels is similar to the definition of the Lyapunov Function.
The \(B(t)\) is the funnel. It is a region of attraction of the dynamical system \(F(x,t)\) at time t. e.g. for \(x_0 \in B(0)\), Let \(x_n = F \circ F \circ F … F(x_0)\) (forward integrate n times), then \(x_n \in B(n)\).
The region-of-attraction \(B(t)\) is parameterized as a sub-level set.
Using this parameterization, authors used a sufficient condition for the region-of-attraction,
Here, \(\dot V(x, t) \leq \dot \rho(t)\) states, on the boundary of the region, the dynamic grows slower than the region. As a result, a state which starts in the region stays inside the region because when a state moves to the boundary, the dynamic pushes it back.
Note: the condition is not the the Lyapunov condition, which prove the stability of a system. The dynamic can be unstable in the region. e.g. oscillating within the region.
For the first condition \(V(x,t) \leq 0\), it is a direct result of the LQR regulator. i.e. When a LQR problem is solvable, it’s cost function is positive-definite.
7.2 Use SOS to compute funnels
The interesting part is using the SOS programming to solve \(\dot V(x, t) \leq \dot \rho(t)\).
Here the authors want:
- satisfy \(q(x,t) = \dot V(x, t) – \dot \rho(t) \leq 0\).
- inside the time region \(g_1(t) = t_{k+1} – t > 0\) and \(g_2(t) = t – t_k > 0\).
- on the boundary of \(B(t)\), i.e. \(e(x,t) = J(x,t) – \rho(t) = 0\).
The authors applied the generalized S-procedure, the problem become,
\(L(x,t) = q(x,t) +h_1(x,t)e(x,t) + h_2(x,t)g_1(t) + h_3(x,t)g_2(t) < 0\)
\(h_2(x,t) > 0\)
\(h_3(x,t) > 0\)
Note \(h_1(x,t)\) is a polynomial multiplier for equality constraints.
The author checked the above conditions by checking a sufficient SOS program.
\(L(x,t) = q(x,t) +h_1(x,t)e(x,t) + h_2(x,t)g_1(t) + h_3(x,t)g_2(t) \in -SOS\)
\(h_2(x,t) \in SOS\)
\(h_3(x,t) \in SOS\)
Note, the \(J(x,t) = V(x,t)\) is the output of LQR. It’s a discrete function of time. But the SOS need continuous polynomial inputs.
The author did an interpolation to convert \(J(x,t)\) to a continuous polynomial function \(J^*(x,t)\) and they were able to do the SOS program.
The “size” of the Funnels \(\rho_i\) is parameterized by a piecewise-linear function.
The final objective is to find the maximum the size of a Funnel subject to the Funnel constraints.
By solving the above SOS program, the authors found a Funnel for the dynamic system with a LQR controller applied.
Go back to the tree-grow algorithm, since there is a LQR controller which works in a Funnel, we don’t need to sample state from existing Funnels.
8. What is the S-procedure?
We discussed the generalized S-procedure. A natural question is what is the S-procedure?
Lieven Vandenberghe and Stephen Boyd had the answer (again).
The problem is to find a Lyapunov function for a linear dynamic system with bounded control.
The authors used the definition of Lyapunov function \(\dot V(x) < 0\) and do some algebra for the bounded control constrains \(|u(t)| < |y(t)|\). They got a condition for the stability of the system,
So the problem became, for all \((x_i,u_i)\) satisfy \(|u(t)| < |y(t)|\), can we show \(\dot V(x) < 0\)?
It’s hard.
But checking a sufficient condition for the problem, which is similar to the Lagrange method, is feasible.
It’s the S-procedure. It was invented in the Soviet literature in 1944.
9. How to Solve Semidefinite Programs (SDP)?
In previous discussion, we converted a SOS program into a SDP, and used the CVXPY solver to solve the SDP.
How do we solve the SDP ? This paper gave a good overview.
A SDP is defined as,
The dual of SDP is,
Note that \(A \in S^n, X \in S^n, C \in S^n\). i.e. There are symmetric matrices.
In the SDP, the center path is,
Recall from the primal-dual interior points methods, the central path is the optimal condition of the barrier terms. It is basically a numerically friendly Complementary Slackness Condition.
The (relaxed) KKT condition (without inequality constraints) for SDP is given by,
Please compare the KKT conditions for SDP with KKT conditions for a scalar function.
The Primal Dual Interior Points Methods linearized the KKT conditions and solve for a update step. Then, it does a backtracking line search to ensure inequality constraints. Read this article for a simple explanation of the Primal Dual Interior Points Methods.
For the SDP, the (2.1) looks like a system of equations, can we use Newton’s method to solve it?
Yes, we can. With some considerations for the symmetric matrix (read the paper for details), the SDP KKT conditions can be linearized as,
Then, we can solve for the update step for primal and dual variables. And do a linear search to ensure PSD for matrix variables.
10. Final Thoughts
The Sum of Squares programming is a very interesting topic. It has deep connection to classical polynomial theory. If you want to know more, I recommend watching Pablo Parrilo’s lecture.
Personally, I feel like the SOS programming works very well for simple problems, but it may not able to scale to larger problems. For example, selecting basis for the SOS program can be tricky; The number of optimization variables grow quadratically with the size of basis.
Well, the SOS programming may not to very practical but it is cute, sexy and fun 🙂
8 thoughts on “Sum of Squares Programming, Global Non-Convex Optimization, Lyapunov Function and Semidefinite Programming”