## Archive for June, 2011

### Sum of squares decomposition

June 30, 2011

Consider the following quartic form in $\mathbb{R}[x, y]$

$F(x, y) = 2 x^4 + 2 x^3 y - x^2 y^2 + 5 y^4$

and suppose we would like to decide if $F$ is globally nonnegative. In other words, we want to know whether the proposition

$F(x, y) \geq 0, \qquad{} \forall x, y \in \mathbb{R}$

is true. A sufficient condition for $F(x,y)$ to be globally nonnegative is the existence of a sum of squares decomposition [1, 2]

$F(x, y) = \displaystyle \sum_{i} f_i^2 (x, y)$

where $f_i \in \mathbb{R}[x, y]$. If such a decomposition does indeed exist, how do we compute it? How do we find the $f_i$ polynomials?

Note that $F(x, y)$ is  a quartic form, i.e., a polynomial whose monomials all have total degree 4. Using Parrilo’s method [1, 2], let us try to write $F(x, y)$ as a quadratic form in all the monomials of total degree equal to 2, as follows

$F(x,y) = \left[\begin{array}{c} x^2\\ y^2\\ x y\end{array}\right]^T \left[\begin{array}{ccc} q_{11} & q_{12} & q_{13}\\ q_{12} & q_{22} & q_{23}\\ q_{13} & q_{23} & q_{33}\end{array}\right] \left[\begin{array}{c} x^2\\ y^2\\ x y\end{array}\right]$

where the real matrix

$Q := \left[\begin{array}{ccc} q_{11} & q_{12} & q_{13}\\ q_{12} & q_{22} & q_{23}\\ q_{13} & q_{23} & q_{33}\end{array}\right]$

is symmetric. Hence, we obtain

$F(x,y) = q_{11} x^4 + 2 q_{13} x^3 y + (2 q_{12} + q_{33}) x^2 y^2 + 2 q_{23} x y^3 + q_{22} y^4$

and, matching coefficients, we get the equality constraints

$q_{11} = 2$,   $q_{13} = 1$,   $2 q_{12} + q_{33} = -1$,   $q_{23} = 0$,   $q_{22} = 5$.

Let us make $q_{12} = - \lambda$, where $\lambda \in \mathbb{R}$, so that $q_{33} = 2 \lambda - 1$. Hence, we obtain a family of symmetric $Q$ matrices parameterized by $\lambda$

$Q := \left[\begin{array}{ccc} 2 & -\lambda & 1\\ -\lambda & 5 & 0\\ 1 & 0 & 2 \lambda - 1\end{array}\right]$

For each choice of $\lambda$, we get a matrix $Q$. We then conclude that the representation $F(x, y) = z^T Q z$, where $z := (x^2, y^2, x y)$, is not unique. Since there are infinitely many admissible $Q$ matrices, if we are able to find a positive semidefinite one, then $F(x, y) = z^T Q z \geq 0$ for all $z$, and global nonnegativity of the given polynomial $F$ is guaranteed!

For instance, let us try $\lambda = 3$. We then get

$Q = \left[\begin{array}{ccc} 2 & -3 & 1\\ -3 & 5 & 0\\ 1 & 0 & 5\end{array}\right]$

which is positive semidefinite. Computing a Cholesky factorization of $Q$, i.e., finding a matrix $L$ such that $Q = L^T L$, we get

$L = \displaystyle \frac{1}{\sqrt{2}}\left[\begin{array}{ccc} 2 & -3 & 1\\ 0 & 1 & 3\end{array}\right]$

and, thus, $F(x, y) = z^T Q z = z^T L^T L z = \| L z \|_2^2$, which finally yields a sum of squares (SOS) decomposition

$F(x, y) = \displaystyle \frac{1}{2} \left(2 x^2 -3 y^2 + x y\right)^2 + \frac{1}{2} \left(y^2 + 3 x y\right)^2$.

To make sure that this SOS decomposition is correct, we can use a little Python / SymPy script:


from sympy import *

# symbolic variables
x = Symbol('x')
y = Symbol('y')

# define polynomial
F = 2*(x**4) + 2*(x**3)*y - (x**2)*(y**2) + 5*(y**4)

# define SOS polynomials
f1 = (sqrt(2)/2) * (2*(x**2) - 3*(y**2) + x*y)
f2 = (sqrt(2)/2) * (y**2 + 3*x*y)

# construct residual polynomial
R = F - (f1**2 + f2**2)

print "Residual polynomial: %s" % R.expand()



The residual polynomial is $R (x, y) = 0$. The SOS decomposition is correct! We can then conclude that the given polynomial $F$ is globally nonnegative.

__________

Remark: this post is based on example 4.1 of Parrilo’s doctoral dissertation [1] and also on example 3.5 of Parrilo’s paper [2].

__________

References

[1] Pablo Parrilo, Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization, Dissertation (Ph.D.), California Institute of Technology, 2000.

[2] Pablo Parrilo, Semidefinite programming relaxations for semialgebraic problems, Mathematical Programming, 2003.

### Safety verification via barrier certificates

June 19, 2011

Suppose that we are given a dynamical system of the form $\dot{x} (t) = f ( x (t) )$, where $x : [0, \infty) \to \mathbb{R}^n$ is the state trajectory, and $f : \mathbb{R}^n \to \mathbb{R}^n$ is a known vector field [1]. Suppose further that we are given an initial set $\mathcal{X}_0 \subset \mathbb{R}^n$, and an unsafe set $\mathcal{X}_u \subset \mathbb{R}^n$.

If there exists a state trajectory $x : [0, \infty) \to \mathbb{R}^n$ with initial condition $x(0) \in \mathcal{X}_0$ such that $x(T) \in \mathcal{X}_u$ for some $T \geq 0$, we say that the system is unsafe. An example of an unsafe system is depicted below

If there exists no state trajectory $x : [0, \infty) \to \mathbb{R}^n$ with initial condition $x(0) \in \mathcal{X}_0$ such that $x(T) \in \mathcal{X}_u$ for some $T \geq 0$, we say that the system is safe. Note that if $\mathcal{X}_0 \cap \mathcal{X}_u \neq \emptyset$, we can conclude immediately that the system is unsafe, as there exists an initial condition $x(0) \in \mathcal{X}_u$.

The safety verification problem can then be stated as follows:

Given a vector field $f : \mathbb{R}^n \to \mathbb{R}^n$, an initial set $\mathcal{X}_0 \subset \mathbb{R}^n$ and an unsafe set $\mathcal{X}_u \subset \mathbb{R}^n$, decide whether the dynamical system $\dot{x} = f (x)$ is safe.

To be more precise, safety is not a property of the dynamical system $\dot{x} = f (x)$ per se, but rather of the ordered triple $(f, \mathcal{X}_0, \mathcal{X}_u)$.

How do we go about solving the safety verification problem? The most obvious approach would be to “propagate” the initial set $\mathcal{X}_0$ forwards in time, i.e., to compute (forward) reachable sets. However, computing exact reachable sets is seldom possible, which forces us to compute approximations of reachable sets.

When performing safety verification, we are asking a binary question: is $(f, \mathcal{X}_0, \mathcal{X}_u)$ safe? Hence, computation of the reachable sets seems to provide more information than what we actually do need. This is a personal opinion, not a fact.

__________

Barrier certificates

Probably inspired on Lyapunov stability theory [1], Stephen Prajna introduced a few years ago a new method for safety verification of continuous-time dynamical systems that relies on what he termed barrier certificates [2]. A barrier certificate is a function of state $B : \mathbb{R}^n \to \mathbb{R}$ that satisfies the following inequalities

$\begin{array}{rl} B(x) \leq 0 \qquad{} \forall x \in \mathcal{X}_0 \\ B(x) > 0 \qquad{} \forall x \in \mathcal{X}_u\\ \displaystyle\frac{\partial B}{\partial x} (x) f (x) \leq 0 \qquad{} \forall x \in \mathbb{R}^n\end{array}$

Let us now introduce a function of time $b (t) := B (x (t))$ and fix the initial condition $x(0)$. Taking the derivative of $b$ with respect to time, we obtain

$\dot{b} (t) = \displaystyle\frac{\partial B}{\partial x} (x(t)) \dot{x} (t) = \displaystyle\frac{\partial B}{\partial x} (x(t)) f (x(t)) \leq 0$

which follows from the non-positivity of the Lie derivative of $B$ along the flow of $f$. From $\dot{b}(t) \leq 0$, we conclude that $b(t) \leq b(0)$ for all $t \geq 0$, i.e., we have that $B(x(t)) \leq B(x(0))$ for all $t \geq 0$. Hence, the following sublevel set

$\Omega_0 := \{ x \in \mathbb{R}^n \mid B(x) \leq B(x(0))\}$

is positively invariant [1]. In other words, a solution starting in $\Omega_0$ will remain in $\Omega_0$ for all $t \geq 0$. Since $B(x(0)) \leq 0$ for all $x(0) \in \mathcal{X}_0$, we conclude that $B(x(t)) \leq 0$ for all $t \geq 0$. However, since $B(x) > 0$ for all $x \in \mathcal{X}_u$, we can finally conclude that the safety of $(f, \{x(0)\}, \mathcal{X}_u)$ is guaranteed. Alternatively, one could argue that, since, by construction, $\Omega_0 \cap \mathcal{X}_u = \emptyset$, then safety is certified.

and once we know that $\Omega_0$ is a positively invariant set (from $B(x(t)) \leq B(x(0))$), we conclude that the state trajectory must live in $\Omega_0$, which is painted in dark gray below

Please do recall that $B(x(0)) \leq 0$, and also note that the unsafe set is contained in the superlevel set $\{ x \in \mathbb{R}^n \mid B(x) > 0\}$.

So far we assumed that the initial condition was fixed. Henceforth, we shall make the initial condition free to take any values in the initial set $\mathcal{X}_0$. Let us define

$b_0 := \displaystyle\sup_{x_0 \in \mathcal{X}_0} B(x_0)$

so that we get an upper bound on $B(x(0))$, as follows

$B(x(t)) \leq B(x(0)) \leq b_0 \leq 0$

where we used the inequality $B(x) \leq 0$ for all $x \in \mathcal{X}_0$ to conclude that $b_0 \leq 0$. This allows us to conclude that the sublevel set

$\tilde{\Omega}_0 := \{ x \in \mathbb{R}^n \mid B(x) \leq b_0\}$

is positively invariant. Since $\tilde{\Omega}_0 \cap \mathcal{X}_u = \emptyset$, safety of $(f, \mathcal{X}_0, \mathcal{X}_u)$ is guaranteed. Note that the sublevel set $\tilde{\Omega}_0$ is the union of all possible $\Omega_0$ sublevel sets, i.e.,

$\tilde{\Omega}_0 = \displaystyle\bigcup_{x_0 \in \mathcal{X}_0} \{ x \in \mathbb{R}^n \mid B(x) \leq B(x_0)\}$.

This new sublevel set is depicted below, painted in dark gray

Note that $\mathcal{X}_0 \subset \tilde{\Omega}_0$. Essentially, the zero level set of the barrier certificate $B$ separates an unsafe set, $\mathcal{X}_u$, from all system trajectories starting in the initial set $\mathcal{X}_0$, which are contained in $\tilde{\Omega}_0$, thus providing a proof of safety without requiring the computation of reachable sets. Personally, I like to think of this approach as a sophisticated generalization of the well-known separating hyperplane theorem.

At this point, you may be thinking that this framework is too good to be true. Instead of computing the flow or the reachable sets, all we need to do is find a function $B : \mathbb{R}^n \to \mathbb{R}$ that satisfies three inequalities, and then (voilà!) infinite-time safety of $(f, \mathcal{X}_0, \mathcal{X}_u)$ is certified! As it usually happens, the devil is in the details…

If finding Lyapunov functions [1] for low-dimensional dynamical systems is often “hard”, why should we expect the hunt for barrier certificates to be an “easy” endeavor? Quel panache! Are we merely transferring the complexity of the problem from the computation of reachable sets to the satisfaction of three inequalities? Will we fail to find barrier certificates due to excessive conservativeness? Should we abandon all hope? Please be patient, dear reader, for we will address these issues and alleviate your doubts in future posts ;-)

Next installment will be devoted to an example. до скорой встречи!

_________

References

[1] Hassan K. Khalil, Nonlinear Systems, 3rd edition, Prentice Hall.

[2] Stephen Prajna, Optimization-based methods for nonlinear and hybrid systems verification, Dissertation (Ph.D.), California Institute of Technology, 2005.