Global Stability

The notion of a fixed point being stable i.s.L. is inherently a local notion of stability (defined with $\epsilon$ and $\delta$ balls around the origin), but the notions of asymptotic and exponential stability can be applied globally. The Lyapunov theorems work for this case, too, with only minor modification.

The new condition, on the behavior as $||\bx|| \rightarrow \infty$ is known as "radially unbounded", and is required to make sure that trajectories cannot diverge to infinity even as $V$ decreases; it is only required for global stability analysis. This radially unbounded condition is sufficient but is not necessary; other conditions can potentially be invoked to guarantee that the system avoids this pathology. More essential is that the resulting trajectory $\bx(t)$ is bounded for all $t$.

LaSalle's Invariance Principle

Perhaps you noticed the disconnect between the statement above and the argument that we made for the stability of the pendulum. In the pendulum example, using the mechanical energy resulted in a Lyapunov function time derivative that was only negative semi-definite, but we eventually argued that the fixed points were asymptotically stable. That took a little extra work, involving an argument about the fact that the fixed points were the only place that the system could stay with $\dot{E}=0$; every other state with $\dot{E}=0$ was only transient. We can formalize this idea for the more general Lyapunov function statements--it is known as LaSalle's Theorem.

To be clear, an invariant set, ${\cal G}$, of the dynamical system is a set for which $\bx(0)\in{\cal G} \Rightarrow \forall t>0, \bx(t) \in {\cal G}$. In other words, once you enter the set you never leave. The "largest invariant set" with $\dot{V}(\bx) = 0$ need not be connected; in fact for the pendulum example each fixed point is an invariant set with $\dot{V} = 0$, so the largest invariant set is the union of all the fixed points of the system. There are also variants of LaSalle's Theorem which work over a region. For example, the total energy of the pendulum does not actually satisfy the radially-unbounded condition (the energy is finite as $\theta \rightarrow \infty$), but observing that $\bx(t)$ is bounded for all initial conditions is sufficient to invoke a variant of LaSalle's theorem Khalil01.

Finding a Lyapunov function which $\dot{V} \prec 0$ is more difficult than finding one that has $\dot{V} \preceq 0$. LaSalle's theorem gives us the ability to make a statement about asymptotic stability even in this case. In the pendulum example, every state with $\dot\theta=0$ had $\dot{E}=0$, but only the fixed points are in the largest invariant set.

Relationship to the Hamilton-Jacobi-Bellman equations

At this point, you might be wondering if there is any relationship between Lyapunov functions and the cost-to-go functions that we discussed in the context of dynamic programming. After all, the cost-to-go functions also captured a great deal about the long-term dynamics of the system in a scalar function. We can see the connection if we re-examine the HJB equation \[ 0 = \min_\bu \left[ \ell(\bx,\bu) + \pd{J^*}{\bx}f(\bx,\bu). \right] \]Let's imagine that we can solve for the optimizing $\bu^*(\bx)$, then we are left with $ 0 = \ell(\bx,\bu^*) + \pd{J^*}{\bx}f(\bx,\bu^*) $ or simply \[ \dot{J}^*(\bx) = -\ell(\bx,\bu^*) \qquad \text{vs} \qquad \dot{V}(\bx) \preceq 0. \] In other words, in optimal control we must find a cost-to-go function which matches this gradient for every $\bx$; that's very difficult and involves solving a potentially high-dimensional partial differential equation. By contrast, Lyapunov analysis is asking for much less - any function which is going downhill (at any rate) for all states. This can be much easier, for theoretical work, but also for our numerical algorithms. Also note that if we do manage to find the optimal cost-to-go, $J^*(\bx)$, then it can also serve as a Lyapunov function so long as $\ell(\bx,\bu^*(\bx)) \succeq 0$.

Lyapunov functions for estimating regions of attraction

There is another very important connection between Lyapunov functions and the concept of an invariant set: any sublevel set of a Lyapunov function is also an invariant set. This gives us the ability to use sublevel sets of a Lyapunov function as approximations of the region of attraction for nonlinear systems. For a locally attractive fixed point, $\bx^*$, the region of attraction to $\bx^*$ is the largest set ${\cal D} \subseteq {\cal X}$ for which $\bx(0) \in {\cal D} \Rightarrow \lim_{t\rightarrow \infty} \bx(t) = \bx^*.$ Regions of attraction are always open, connected, and invariant sets (with boundaries ).

Region of attraction for a one-dimensional system

Consider the first-order, one-dimensional system $\dot{x} = -x + x^3.$ We can quickly understand this system using our tools for graphical analysis.

In the vicinity of the origin, $\dot{x}$ looks like $-x$, and as we move away it looks increasingly like $x^3$. There is a stable fixed point at the origin and unstable fixed points at $\pm 1$. In fact, we can deduce visually that the region of attraction to the stable fixed point at the origin is $\bx \in (-1,1)$. Let's see if we can demonstrate this with a Lyapunov argument.

First, let us linearize the dynamics about the origin and use the Lyapunov equation for linear systems to find a candidate $V(\bx)$. Since the linearization is $\dot{x} = -x$, if we take ${\bf Q}=1$, then we find ${\bf P}=\frac{1}{2}$ is the positive definite solution to the algebraic Lyapunov equation (\ref{eq:algebraic_lyapunov}). Proceeding with $$V(\bx) = \frac{1}{2} x^2,$$ we have $$\dot{V} = x (-x + x^3) = -x^2 + x^4.$$ This function is zero at the origin, negative for $|x| < 1$, and positive for $|x| > 1$. Therefore we can conclude that the sublevel set $V < \frac{1}{2}$ is invariant and the set $x \in (-1,1)$ is inside the region of attraction of the nonlinear system. In fact, this estimate is tight.

Robustness analysis using "common Lyapunov functions"

While we will defer most discussions on robustness analysis until later in the notes, the idea of a common Lyapunov function, which we introduced briefly for linear systems in the example above, can be readily extended to nonlinear systems and region of attraction analysis. Imagine that you have a model of a dynamical system but that you are uncertain about some of the parameters. For example, you have a model of a quadrotor, and are fairly confident about the mass and lengths (both of which are easy to measure), but are not confident about the moment of inertia. One approach to robustness analysis is to define a bounded uncertainty, which could take the form of $$\dot{\bx} = f_\alpha(\bx), \quad \alpha_{min} \le \alpha \le \alpha_{max},$$ with $\alpha$ a vector of uncertain parameters in your model. Richer specifications of the uncertainty bounds are also possible, but this will serve for our examples.

In standard Lyapunov analysis, we are searching for a function that goes downhill for all $\bx$ to make statements about the long-term dynamics of the system. In common Lyapunov analysis, we can make many similar statements about the long-term dynamics of an uncertain system if we can find a single Lyapunov function that goes downhill for all possible values of $\alpha$. If we can find such a function, then we can use it to make statements with all of the variations we've discussed (local, regional, or global; in the sense of Lyapunov, asymptotic, or exponential).

A one-dimensional system with gain uncertainty

Let's consider the same one-dimensional example used above, but add an uncertain parameter into the dynamics. In particular, consider the system: $$\dot{x} = -x + \alpha x^3, \quad \frac{3}{4} < \alpha < \frac{3}{2}.$$ Plotting the graph of the one-dimensional dynamics for a few values of $\alpha$, we can see that the fixed point at the origin is still stable, but the robust region of attraction to this fixed point (shaded in blue below) is smaller than the region of attraction for the system with $\alpha=1$.

Taking the same Lyapunov candidate as above, $V = \frac{1}{2} x^2$, we have $$\dot{V} = -x^2 + \alpha x^4.$$ This function is zero at the origin, and negative for all $\alpha$ whenever $x^2 > \alpha x^4$, or $$|x| < \frac{1}{\sqrt{\alpha_{max}}} = \sqrt{\frac{2}{3}}.$$ Therefore, we can conclude that $|x| < \sqrt{\frac{2}{3}}$ is inside the robust region of attraction of the uncertain system.

Not all forms of uncertainty are as simple to deal with as the gain uncertainty in that example. For many forms of uncertainty, we might not even know the location of the fixed points of the uncertain system. In this case, we can often still use common Lyapunov functions to give some guarantees about the system, such as guarantees of robust set invariance. For instance, if you have uncertain parameters on a quadrotor model, you might be ok with the quadrotor stabilizing to a pitch of $0.01$ radians, but you would like to guarantee that it definitely does not flip over and crash into the ground.

A one-dimensional system with additive uncertainty

Now consider the system: $$\dot{x} = -x + x^3 + \alpha, \quad -\frac{1}{4} < \alpha < \frac{1}{4}.$$ Plotting the graph of the one-dimensional dynamics for a few values of $\alpha$, this time we can see that the fixed point is not necessarily at the origin; the location of the fixed point moves depending on the value of $\alpha$. But we should be able to guarantee that the uncertain system will stay near the origin if it starts near the origin, using an invariant set argument.

Taking the same Lyapunov candidate as above, $V = \frac{1}{2} x^2$, we have $$\dot{V} = -x^2 + x^4 + \alpha x.$$ This Lyapunov function allows us to easily verify, for instance, that $V \le \frac{1}{3}$ is a robust invariant set, because whenever $V = \frac{1}{3}$, we have $$\forall \alpha \in [\alpha_{min},\alpha_{max}],\quad \dot{V}(x,\alpha) < 0.$$ Therefore $V$ can never start at less than one-third and cross over to greater than one-third. To see this, see that $$ V=\frac{1}{3} \quad \Rightarrow \quad x = \pm \sqrt{\frac{2}{3}} \quad \Rightarrow \quad \dot{V} = -\frac{2}{9} \pm \alpha \sqrt{\frac{2}{3}} < 0, \forall \alpha \in \left[-\frac{1}{4},\frac{1}{4} \right]. $$ Note that not all sublevel sets of this invariant set are invariant. For instance $V < \frac{1}{32}$ does not satisfy this condition, and by visual inspection we can see that it is in fact not robustly invariant.

Barrier functions

Lyapunov functions are used to certify stability or to establish invariance of a region. But the same conditions can be used to certify that the state of a dynamical system will not visit some region of state space. This can be useful for e.g. verifying that your autonomous car will not crash or, perhaps, that your walking robot will not fall down. When used this way, we call the associated functions, $\mathcal{B}(\bx)$, "barrier functions".

Natural extensions exist for certifying these conditions over a region, considering piecewise-polynomial dynamics, worse-case robustness, etc.

Linear systems

The Lyapunov conditions we've explored above take a particularly nice form for linear systems.

Note that the radially-unbounded condition is satisfied by ${\bf P} \succ 0$, and that the derivative condition is equivalent to the matrix condition $${\bf PA} + {\bf A}^T {\bf P} \prec 0.$$

For stable linear systems the existence of a quadratic Lyapunov function is actually a necessary (as well as sufficient) condition. Furthermore, if $\bA$ is stable, then a Lyapunov function can always be found by finding the positive-definite solution to the matrix Lyapunov equation \begin{equation}{\bf PA} + {\bf A}^T{\bf P} = - {\bf Q},\label{eq:algebraic_lyapunov} \end{equation} for any ${\bf Q}={\bf Q}^T\succ 0$.

This is a very powerful result - for nonlinear systems it will be potentially difficult to find a Lyapunov function, but for linear systems it is straightforward. In fact, this result is often used to propose candidate Lyapunov functions for nonlinear systems, e.g., by linearizing the equations and solving a local linear Lyapunov function which should be valid in the vicinity of a fixed point.

Lyapunov analysis for linear systems has an extremely important connection to convex optimization. In particular, we could have also formulated the Lyapunov conditions for linear systems above using semi-definite programming (SDP). Semidefinite programming is a subset of convex optimization -- an extremely important class of problems for which we can produce efficient algorithms that are guaranteed find the global optima solution (up to a numerical tolerance and barring any numerical difficulties).

If you don't know much about convex optimization or want a quick refresher, please take a few minutes to read the optimization preliminaries in the appendix. The main requirement for this section is to appreciate that it is possible to formulate efficient optimization problems where the constraints include specifying that one or more matrices are positive semi-definite (PSD). These matrices must be formed from a linear combination of the decision variables. For a trivial example, the optimization $$\min_a a,\quad \subjto \begin{bmatrix} a & 0 \\ 0 & 1 \end{bmatrix} \succeq 0,$$ returns $a = 0$ (up to numerical tolerances).

The value in this is immediate for linear systems. For example, we can formulate the search for a Lyapunov function for the linear system $\dot\bx = {\bf A} \bx$ by using the parameters ${\bf p}$ to populate a symmetric matrix ${\bf P}$ and then write the SDP: \begin{equation} \find_{\bf p} \quad \subjto \quad {\bf P} \succeq 0, \quad {\bf PA} + {\bf A}^T {\bf P} \preceq 0.\label{eq:lyap} \end{equation} Note that you would probably never use that particular formulation, since there specialized algorithms for solving the simple Lyapunov equation which are more efficient and more numerically stable. But the SDP formulation does provide something new -- it will give us our first example of robustness analysis. In particular, we can now easily formulate the search for a "common Lyapunov function" for uncertain linear systems.

Common Lyapunov analysis for linear systems

Suppose you have a system governed by the equations $\dot\bx = {\bf A}\bx$, where the matrix ${\bf A}$ is unknown, but its uncertain elements can be bounded. There are a number of ways to write down this uncertainty set; let us choose to write this by describing ${\bf A}$ as the convex combination of a number of known matrices, $${\bf A} = \sum_{i} \beta_i {\bf A}_i, \quad \sum_i \beta_i = 1, \quad \forall i, \beta_i > 0.$$ This is just one way to specify the uncertainty; geometrically it is describing a polygon of uncertain parameters (in the space of elements of ${\bf A}$ with each ${\bf A}_i$ as one of the vertices in the polygon.

Now we can formulate the search for a common Lyapunov function using \[ \find_{\bf p} \quad \subjto \quad {\bf P} \succeq 0, \quad \forall_i, {\bf PA}_i + {\bf A}_i^T {\bf P} \preceq 0.\] The solver will then return a matrix ${\bf P}$ which satisfies all of the constraints, or return saying "problem is infeasible". It can easily be verified that if ${\bf P}$ satisfies the Lyapunov condition at all of the vertices, then it satisfies the condition for every ${\bf A}$ in the set: \[ {\bf P}(\sum_i \beta_i {\bf A}_i) + (\sum_i \beta_i {\bf A}_i)^T {\bf P} = \sum_i \beta_i ({\bf P A}_i + {\bf A}_i^T {\bf P}) \preceq 0, \] since $\forall i$, $\beta_i > 0$. Note that, unlike the simple Lyapunov equation for a known linear system, this condition being satisfied is a sufficient but not a necessary condition -- it is possible that the set of uncertain matrices ${\bf A}$ is robustly stable, but that this stability cannot be demonstrated with a common quadratic Lyapunov function.

You can try this example for yourself in .

As always, make sure that you open up the code and take a look.

Note that the argument above does not actually make use of the $\sum_i \beta_i = 1;$ that assumption can be safely removed. So really, our common Lyapunov function, when found, has certified the convex cone with the $A_i$ matrices as the extreme rays. This is not surprising, since if we have any linear system/Lyapunov quadratic form pair $A, P$, then $P$ also confirms that $\beta A$ is stable for any $\beta > 0.$ Multiplying a system matrix by a scalar will not change the sign of its eigenvalues.

There are many small variants of this result that are potentially of interest. For instance, a very similar set of conditions can certify "mean-square stability" for linear systems with multiplicative noise (see e.g. Boyd94, § 9.1.1).

This example is very important because it establishes a connection between Lyapunov functions and (convex) optimization. But so far we've only demonstrated this connection for linear systems where the PSD matrices provide a magical recipe for establishing the positivity of the (quadratic) functions for all $\bx$. Is there any hope of extending this type of analysis to more general nonlinear systems? Surprisingly, it turns out that there is.

Global analysis for polynomial systems

Sums of squares (SOS) optimization provides a natural generalization of the SDP formulation to nonlinear systems. In SDP, we were allowed to write constraints of the form, ${\bf P} \succeq 0$. In sums of squares optimization, we generalize this to writing constraints of the form "$p_\balpha(\bx)$ is SOS", where $p(\bx)$ is a polynomial with coefficients $\balpha$. Let's examine this in the specific context of Lyapunov functions.

In linear systems, we searched over $V(\bx) = \bx^T {\bf P} \bx, {\bf P} \succeq 0.$ In the nonlinear generalization, we search over $$V(\bx) = {\bf m}^T(\bx) {\bf P} {\bf m}(\bx),\quad {\bf P} \succeq 0,$$ with ${\bf m}(\bx)$ a set of nonlinear basis functions. This is a linear (in the parameters ${\bf P}$) function approximator which is guaranteed to be strictly positive by virtue of being a sum of squares. While the positivity of this sums of squares is valid for any choice of ${\bf m}(\bx)$, the recipe is particularly powerful for polynomials, because for any polynomial, $p(\bx)$, we can confirm its positivity by creating a surrogate sums-of-squares representation, where we have strong recipes for automatically choosing good bases ${\bf m}(\bx)$ given a polynomial $p(\bx)$, and we can constrain that $p(\bx) = {\bf m}^T(\bx) {\bf P} {\bf m}(\bx)$ with linear constraints relating the coefficients of $p(\bx)$ with the coefficient in ${\bf P}$.

This suggests that it may be possible to generalize the optimization approach using SDP to search for Lyapunov functions for linear systems to searching for Lyapunov functions for nonlinear systems: $\dot\bx = \bf{f}(\bx),$ where $\bf{f}$ is a vector-valued polynomial function. Importantly, if $V(\bx)$ is polynomial, and $\bf{f}(\bx)$ is polynomial, then $\dot{V}(\bx) = \pd{V}{\bx}{\bf f}(\bx)$ is also polynomial.

If we parameterize a fixed-degree Lyapunov candidate as a polynomial with unknown coefficients, e.g., \[ V_\alpha(\bx) = \alpha_0 + \alpha_1 x_1 + \alpha_2 x_2 + \alpha_3 x_1x_2 + \alpha_4 x_1^2 + ..., \] then we can formulate the optimization: \begin{align*} \find_\alpha, \quad \subjto \quad& V_\alpha(\bx) \sos, \quad V_\alpha(0) = 0 \\ & -\dot{V}_\alpha(\bx) = -\pd{V_\alpha}{\bx} f(\bx) \sos. \end{align*} Because this is a convex optimization, the solver will return a solution if one exists, and this solution would certify stability i.s.L. of the origin. Note that I don't actually need to include the constraint $\dot{V}_\alpha({\bf 0}) = 0$; if $V_\alpha$ is a smooth function, is zero at zero, and is positive for all $\bx$, then this already implies that $\pd{V_\alpha}{\bx} = 0.$

If we wish to demonstrate (global) asymptotic stability, then we can write the constraint $\dot{V}_\alpha(\bx) \prec 0$ as $$-\dot{V}_\alpha(\bx) - \epsilon \bx^T\bx \text{ is SOS,}$$ where $\epsilon$ is a small positive constant that need only be larger than the numerical tolerances of the solver. The conditions for exponential stability fit easily, too.

It is important to remember that there are a handful of gaps which make the existence of this solution a sufficient condition (to guarantee that $V(\bx)$ is a Lyapunov function for $\bf{f}(\bx)$) instead of a necessary one. First, there is no guarantee that a stable polynomial system can be verified using a polynomial Lyapunov function (of any degree), and in fact there are known counter-examples Ahmadi11a; here we are only searching over a fixed-degree polynomial. Second, even if a polynomial Lyapunov function does exist, there is a gap between the SOS polynomials and the positive polynomialsBlekherman13.

Despite these caveats, I have found this formulation to be surprisingly effective in practice. Intuitively, I think that this is because there is relatively a lot of flexibility in the Lyapunov conditions -- if you can find one function which is a Lyapunov function for the system, then there are also many "nearby" functions which will satisfy the same constraints. This is the fundamental benefit of changing the equality in the HJB into an inequality in Lyapunov analysis.

Region of attraction estimation for polynomial systems

Now we have arrived at the tool that I believe can be a work-horse for many serious robotics applications. Most of our robots are not actually globally stable. That's not because they are robots -- if you push me hard enough, I will fall down, too! Asking for global stability is too much, but we can instead certify that our closed-loop systems behave well over a sufficiently large region of state space (and potentially for a range of possible parameters or disturbances). One of the simplest versions of this is estimating the region of attraction to a fixed point.

Sums-of-squares optimization effectively gives us an oracle which we can ask: is this polynomial positive for all $\bx$? To use this for regional analysis, we have to figure out how to modify our questions to the oracle so that the oracle will say "yes" or "no" when we ask if a function is positive over a certain region which is a subset of $\Re^n$. That trick is called the S-procedure. It is closely related to the Lagrange multipliers from constrained optimization, and has deep connections to "Positivstellensatz" from algebraic geometry.

Convex outer approximations

All of our region of attraction approximations so far have been "inner approximations" -- the certified region is guaranteed to be contained within the true region of attraction for the system. This is what we want if we are to give some guarantee of stability. As we increase the degree of our polynomial and multipliers, we expect that the estimated regions of attraction will converge on the true regions.

It turns out that if we instead consider outer approximations, which converge on the true ROA from the other side, we can write formulations that enable searching for the Lyapunov function directly as a convex optimization. As we will see below, this approach also allows one to write convex formulations for controller design. These methods have been explored in a very nice line of work by Henrion and Korda (e.g. Henrion14), using the method of moments from Lasserre10 also called "occupation measures". Their treatment emphasizes infinite-dimensional linear programs and heirarchies of LMI relaxations; these are just dual formulations to the SOS optimizations that we've been writing here. I find the SOS form more clear, so will stick to it here. Some of the occupation measure papers include the dual formulation of the infinite-dimension linear program which will look quite similar; my coauthors and I typically call out the sums-of-squares version in our papers on the topic (e.g. Posa17).

To find an outer approximation, instead of solving for a Lyapunov function that certifies convergence to the origin, we instead search for a Lyapunov-like "barrier certificate", $\mathcal{B}(\bx).$ Like a Lyapunov function, we'd like $\dot{\mathcal{B}}(\bx) \leq 0$; this time we'll ask for this to be true everywhere (or at least in some set that is sufficiently larger than the ROA of interest). Then we will set $\mathcal{B}(0) = 1 > 0$. If we can find such a function, then certainly any state for which $\mathcal{B}(\bx) < 0$ is outside of the region of attraction of the fixed point -- since $\mathcal{B}$ cannot increase it can never reach the origin which has $\mathcal{B} > 0.$ The superlevel set $\{ \bx | \mathcal{B}(\bx) \ge 0\}$ is an outer approximation of the true region of attraction: $$-\dot{\mathcal{B}}(\bx) \text{ is SOS} \quad \text{and} \quad \mathcal{B}(0) = 1.$$

In order to find the smallest such outer approximation (given a fixed degree polynomial), we choose an objective that tries to "push-down" on the value of $\mathcal{B}(\bx)$. We typically accomplish this by introducing another polynomial function $W(\bx)$ with the requirements that $$W(\bx) \ge 0 \quad \text{and} \quad W(\bx) \ge \mathcal{B}(\bx) +1,$$ implemented as SOS constraints. Then we minimize the integral $\int_\bx W(\bx)d\bx$, which can be readily computed over a compact set like a ball in $\Re^n$ Folland01 and is linear in the coefficients. More sophisticated alternatives also existHenrion09. Majumdar16a writes the objective as minimizing the volume of the covering ellipse, which is convexBoyd04a.

Outer approximation for the cubic polynomial

Let's revisit my favorite scalar dynamical system, $\dot{x} = -x + x^3,$ which has a region of attraction for the fixed point at zero of $\{x \| |x|< 1 \}.$ This time, we'll estimate the ROA using the outer approximation. We can accomplish this with the following program: \begin{align*} \min_{\mathcal{B}(x), W(x)} \quad & \int_{-2}^2 W(x)dx, \\ \subjto \quad & -\dot{\mathcal{B}}(x) & \text{is SOS}, \\ & W(x) &\text{ is SOS}, \\ & W(x) - \mathcal{B}(x) - 1.0 & \text{ is SOS}, \\ & \mathcal{B}(0)= 0.1.\end{align*}

To make the problem a little numerically better, you'll see in the code that I've asked for $\dot{\mathcal{B}}(x)$ to be strictly negative definite, and I've chosen to only include even-degree monomials in $\mathcal{B}(x)$ and $W(x)$. Plotting the solution reveals:

As you can see, the superlevel set, $\mathcal{B}(x) \ge 0$ is a tight outer-approximation of the true region of attraction. In the limit of increasing the degrees of the polynomials to infinity, we would expect that $W(\bx)$ would converge to the indicator function that is one inside the region of attraction, and zero outside (we are quite far from that here, but nevertheless have a tight approximation from $\mathcal{B}(\bx) \ge 0$).

Go ahead and try adjusting some of these decisions. For example, if you set $B(0)=1$, then you can increase the degree of $B(x)$ to get tighter and tighter approximations.

x = Variable("x")
sys = SymbolicVectorSystem(state=[x], dynamics=[-x+x**3])
context = sys.CreateDefaultContext()
V = RegionOfAttraction(sys, context)
          

Robustness analysis using the S-procedure

Piecewise-polynomial systems

Rigid-body dynamics are (rational) polynomial

We've been talking a lot in this chapter about numerical methods for polynomial systems. But even our simple pendulum has a $\sin\theta$ in the dynamics. Have I been wasting your time? Must we just resort to polynomial approximations of the non-polynomial equations? It turns out that our polynomial tools can perform exact analysis of the manipulation equation for almost all^††the most notable exception to this is if your robot has helical/screw joints (see Wampler11).of our robots. We just have to do a little more work to reveal that structure.

Let us first observe that rigid-body kinematics are polynomial (except the helical joint). This is fundamental -- the very nature of a "rigid body" assumption is that Euclidean distance is preserved between points on the body; if $\bp_1$ and $\bp_2$ are two points on a body, then the kinematics enforce that $|\bp_1 - \bp_2|_2^2$ is constant -- these are polynomial constraints. Of course, we commonly write the kinematics in a minimal coordinates using $\sin\theta$ and $\cos\theta$. But because of rigid body assumption, these terms only appear in the simplest forms, and we can simply make new variables $s_i = \sin\theta_i, c_i = \cos\theta_i$, and add the constraint that $s_i^2 + c_i^2 = 1.$ For a more thorough discussion see, for instance, Wampler11 and Sommese05. Since the potential energy of a multi-body system is simply an accumulation of weight times the vertical position for all of the points on the body, the potential energy is polynomial.

If configurations (positions) of our robots can be described by polynomials, then velocities can as well: forward kinematics $\bp_i = f(\bq)$ implies that $\dot\bp_i = \frac{\partial f}{\partial \bq}\dot\bq,$ which is polynomial in $s, c, \dot\theta$. Since the kinetic energy of our robot is given by the accumulation of the kinetic energy of all the mass, $T = \sum_i \frac{1}{2} m_i v_i^Tv_i,$ the kinetic energy is polynomial, too (even when we write it with inertial matrices and angular velocities).

Finally, the equations of motion can be obtained by taking derivatives of the Lagrangian (kinetic minus potential). These derivatives are still polynomial!

Global stability of the simple pendulum via SOS

We opened this chapter using our intuition about energy to discuss stability on the simple pendulum. Now we'll replace that intuition with convex optimization (because it will also work for more difficult systems where our intuition fails).

Let's change coordinates from $[\theta,\dot\theta]^T$ to $\bx = [s,c,\dot\theta]^t$, where $s \equiv \sin\theta$ and $c \equiv \cos\theta$. Then we can write the pendulum dynamics as $$\dot\bx = \begin{bmatrix} c \dot\theta \\ -s \dot\theta \\ -\frac{1}{m l^2} \left( b \dot\theta + mgls \right) \end{bmatrix}.$$

Now let's parameterize a Lyapunov candidate $V(s,c,\dot\theta)$ as the polynomial with unknown coefficients which contains all monomials up to degree 2: $$V = \alpha_0 + \alpha_1 s + \alpha_2 c + ... \alpha_{9} s^2 + \alpha_{10} sc + \alpha_{11} s\dot\theta.$$ Now we'll formulate the feasibility problem: \[ \find_{\bf \alpha} \quad \subjto \quad V \sos, \quad -\dot{V} \sos.\] In fact, this is asking too much -- really $\dot{V}$ only needs to be negative when $s^2+c^2=1$. We can accomplish this with the S-procedure, and instead write \[ \find_{{\bf \alpha},\lambda} \quad \subjto \quad V \sos, \quad -\dot{V} -\lambda(\bx)(s^2+c^2-1) \sos.\] (Recall that $\lambda(\bx)$ is another polynomial with free coefficients which the optimization can use to make terms arbitrarily more positive when $s^2+c^2 \neq 1$.) Finally, for style points, in the code example in we ask for exponential stability:

As always, make sure that you open up the code and take a look. The result is a Lyapunov function that looks familiar (visualized as a contour plot here):

Aha! Not only does the optimization finds us coefficients for the Lyapunov function which satisfy our Lyapunov conditions, but the result looks a lot like mechanical energy. In fact, the result is a little better than energy... there are some small extra terms added which prove exponential stability without having to invoke LaSalle's Theorem.

The one-degree-of-freedom pendulum did allow us to gloss over one important detail: while the manipulator equations $\bM(\bq) \ddot{\bq} + \bC(\bq, \dot\bq)\dot{\bq} = ...$ are polynomial, in order to solve for $\ddot{\bq}$ we actually have to multiply both sides by $\bM^{-1}$. This, unfortunately, is not a polynomial operation, so in fact the final dynamics of the multibody systems are rational polynomial. Not only that, but evaluating $\bM^{-1}$ symbolically is not advised -- the equations get very complicated very fast. But we can actually write the Lyapunov conditions using the dynamics in implicit form, e.g. by writing $V(\bq,\dot\bq,\ddot\bq)$ and asking it to satisfy the Lyapunov conditions everywhere that $\bM(\bq)\ddot\bq + ... = ... + {\bf B}\bu$ is satisfied, using the S-procedure.

Alternatives for obtaining polynomial equations

In practice, you can also Taylor approximate any smooth nonlinear function using polynomials. This can be an effective strategy in practice, because you can limit the degrees of the polynomial, and because in many cases it is possible to provide conservative bounds on the errors due to the approximation.

One final technique that is worth knowing about is a change of coordinates, often referred to as the stereographic projection, that provides a coordinate system in which we can replace $\sin$ and $\cos$ with polynomials:

By projecting onto the line, and using similar triangles, we find that $p = \frac{\sin\theta}{1 + \cos\theta}.$ Solving for $\sin\theta$ and $\cos\theta$ reveals $$\cos\theta = \frac{1-p^2}{1+p^2}, \quad \sin\theta = \frac{2p}{1+p^2}, \quad \text{and} \quad \frac{\partial p}{\partial \theta} = \frac{1+p^2}{2},$$ where $\frac{\partial p}{\partial \theta}$ can be used in the chain rule to derive the dynamics $\dot{p}$. Although the equations are rational, they share the denominator $1+p^2$ and can be treated efficiently in mass-matrix form. Compared to the simple substitution of $s=\sin\theta$ and $c=\cos\theta$, this is a minimal representation (scalar to scalar, no $s^2+c^2=1$ required); unfortunately it does have a singularity at $\theta=\pi$, so likely cannot be used for global analysis.

Verifying dynamics in implicit form

Typically we write our differential equations in the form $\dot\bx = {\bf f}(\bx, \bu).$ But let us consider for a moment the case where the dynamics are given in the form $${\bf g}(\bx, \bu, \dot\bx ) = 0.$$ This form is strictly more general because I can always write ${\bf g}(\bx, \bu, \dot\bx) = f(\bx \bu) - \dot\bx$, but importantly here I can also write the bottom rows of ${\bf g}$ as $\bM(\bq)\ddot\bq + \bC(\bq,\dot\bq)\dot\bq - \btau_g - \bB \bu$. This form can also represent differential algebraic equations (DAEs) which are more general than ODEs; $\bg$ could even include algebraic constraints like $s_i^2 + c_i^2 - 1$. Most importantly, for manipulators, ${\bf g}$ can be polynomial, even if ${\bf f}$ would have been rational. provides access to continuous-time dynamics in implicit form via the CalcImplicitTimeDerivativesResidual method.

We typically write a linear system in implicit form as $${\bf E}\dot\bx = \bA\bx + \bB\bu.$$ Systems in this form are particularly interesting when ${\bf E}$ is singular, and are known as descriptor systems, semistate systems, generalized state-space systems or just singular systems (c.f. Mehrmann91). Just like our standard approach to linearization, we can potentially obtain the matrices ${\bf E}, \bA, \bB$ from a first-order Taylor approximation of the nonlinear equations in ${\bf g}(\bx,\dot\bx,\bu).$ When it comes to Lyapunov analysis, linear systems are often analyzed using the so-called "generalized Lyapunov equation".

This can be verified by taking $$V(\bx) = \bx^T {\bf E}^T {\bf P}{\bf E} \bx,$$ as the Lyapunov function. A corresponding approach can be used (e.g. via sums-of-squares) for implicit nonlinear systems like our manipulator equations where ${\bf g}(\bx,\dot\bx)=0$ can be written in the form ${\bf E}(\bx) \dot\bx = {\bf f}(\bx).$

Alternatively, we can check the Lyapunov conditions, $\dot{V}(\bx) \le 0$, directly on a system in its implicit form, $\bg(\bx,\dot\bx)=0$. Simply define a new function $Q(\bx, \bz) = \frac{\partial V(\bx)}{\partial \bx} \bz.$ If we can show $Q(\bx, \bz) \le 0, \forall \bx,\bz \in \{ \bx, \bz | \bg(\bx,\bz) = 0 \}$ using SOS, then we have verified that $\dot{V}(\bx) \le 0$, albeit at the non-trivial cost of adding indeterminates $\bz$ and an additional S-procedure.

To combine the trigonometric substitutions with the implicit form for a simple second-order system (where the number of positions is equal to the number of velocities), we use ${\bf z}$ to represent the implicit $\ddot{\bq}$ and convert \begin{gather*} V(\bq, \dot\bq) \Rightarrow V({\bf s}, {\bf c}, \dot\bq) \\ {\bf g}(\bq, \dot\bq, {\bf z}) \Rightarrow {\bf g}({\bf s}, {\bf c}, \dot\bq, {\bf z}), \end{gather*} and write the Lyapunov conditions using $\dot{V} = \sum_i \pd{V}{s_i} c_i \dot{q}_i - \sum_i \pd{V}{c_i}s_i \dot{q}_i + \pd{V}{\dot\bq} \bz.$ The most common case for the number of positions is not equal to the number of velocities is when we deal with 3D rotations represented by unit quaternions; this case can be handled using the implicit form, too (we use the implicit form to relate $\dot{q}$ to the velocities and verify the Lyapunov conditions only on the manifold defined by the quaternion unit norm constraint).

Time-varying dynamics and Lyapunov functions

Before focusing on the finite time, let us first realize that the basic (infinite-time) Lyapunov analysis can also be applied to time-varying systems: $\dot{\bx} = f(t, \bx).$ We can analyze the stability (local, regional, or global) of this system with very little change. If we find \begin{gather*}V(\bx) \succ 0, \\ \forall t,\forall \bx \ne 0, \quad \dot{V}(t, \bx) = \pd{V}{\bx}f(t,\bx) < 0, \quad \dot{V}(t, 0) = 0, \end{gather*} then all of our previous statement still hold. In our SOS formulations, $t$ is simply one more indeterminate.

Similarly, even for a time-invariant system, it is also possible to define a time-varying Lyapunov function, $V(t, \bx)$ and establish local, regional, or global stability using the almost the same conditions: \begin{gather*}\forall t, \forall \bx \ne 0, \quad V(t, \bx) > 0, \quad V(t, 0) = 0, \\ \forall t, \forall \bx \ne 0, \quad \dot{V}(t, \bx) = \pd{V}{\bx}f(\bx) + \pd{V}{t} < 0, \quad \dot{V}(t, 0) = 0. \end{gather*}

These two ideas each stand on their own, but they very often go together, as time-varying dynamics are perhaps the best motivator for studying a time-varying Lyapunov function. Afterall, we do know that a stable time-invariant system must have a time-invariant Lyapunov function that demonstrates this stability (from the "converse Lyapunov function" theorems). But we do not know apriori how to represent this function; as an example, remember we know that there are stable polynomial systems that cannot be verified with a polynomial Lyapunov function. For global stability analysis, the time-varying Lyapunov analysis does not add modeling power: since the conditions must be satisfied for all $t$, we could have just set $t$ to a constant and used the time-invariant function. But there may be cases where a time-varying analysis could afford a different analysis for the regional or local stability cases. We will see a good example of this when we study limit cycles.

Finite-time reachability

Finite-time reachability analysis is an important concept for control design and verification, where we seek to understand the behavior of a system over only a finite time interval, $[t_1, t_2]$. It is almost always a region-based analysis, and attempts to make a statement of the form: $$\bx(t_1) \in \mathcal{X}_1 \Rightarrow \bx(t_2) \in \mathcal{X}_2,$$ where $\mathcal{X}_1, \mathcal{X}_2 \subset \Re^n$ are regions of state space. More generally, we might like to understand the time-varying reachable set $\forall t\in [t_1, t_2], \mathcal{X}(t).$

Reachability analysis can be done forward in time: we choose $\mathcal{X}_1$ and try to find the smallest region $\mathcal{X}_2$ for which we can make the statement hold. $\mathcal{X}(t)$ would be called the forward-reachable set (FRS), and can be very useful for certifying e.g. a motion plan. For instance, you might like to prove that your UAV does not crash into a tree in the next 5 seconds. In this case $\mathcal{X}_1$ might be take to be a point representing the current state of the vehicle, or a region representing an outer-approximation of the current state if the state in uncertain. In this context, we would call $\mathcal{X}_2$ a forward-reachable set. In this use case, we would typically choose any approximations in our Lyapunov analysis to certify that an estimate of the reachable region is also an outer-approximation: $\mathcal{X}_2 \subseteq \hat{\mathcal{X}}_2.$

Reachability analysis can also be done backward in time: we choose $\mathcal{X}_2$ and try to maximize the region $\mathcal{X}_1$ for which the statement can be shown to hold. Now $\mathcal{X}(t)$ is called the backward-reachable set (BRS), and for robustness we typically try to certify that our estimates are an inner-approximation, $\hat{\mathcal{X}}_1 \subseteq \mathcal{X}_1.$ The region-of-attraction analysis we studied above can be viewed as a special case of this, with $\mathcal{X}_2$ taken to be the fixed-point, $t_2 = 0$ and $t_1 = -\infty$. But finite-time BRS also have an important role to play, for instance when we are composing multiple controllers in order to achieve a more complicated task, which we will study soon.

Reachability via Lyapunov functions

Lyapunov functions can be used to certify finite-time reachability, even for continuous-time systems. The basic recipe is to certify the Lyapunov conditions over a (potentially time-varying) invariant set. Once again, we typically represent this as a (time-varying) level set of the Lyapunov function containing the origin, $V(\bx) \le \rho(t),$ where $\rho(t)$ is now a positive scalar function of time. Since we already have time as a decision variable, we can easily accommodate time-varying dynamics and Lyapunov functions, as well, so: \begin{align*} \forall t\in[t_1, t_2], \quad \forall \bx, \qquad& V(t,\bx) > 0, \quad& V(t, 0) = 0 \\ \forall t \in [t_1, t_2], \quad \forall \bx \in \{\bx | V(t,\bx) = \rho(t)\}, \qquad& \dot{V}(t, x) \le \dot\rho(t).\end{align*} Note that finite-time reachability is about proving invariance of the set, not stability, and the $\dot{V}$ condition need only be certified at the boundary of the level set. If $V$ is decreasing at the boundary, then trajectories can never leave. One can certainly ask for more -- we may want to show that the system is converging towards $V(t,\bx) = 0$, perhaps even at some rate -- but only invariance is required to certify reachability.

Again, for polynomial systems and dynamics, sums-of-squares optimization is a powerful tool for certifying these properties numerically, and optimizing the volume of the estimated regions. Having taken $t$ as an indeterminate, we can use the S-procedure again to certify the conditions $\forall t \in [t_1, t_2].$

Like in the case for region of attraction, we have many formulations. We can certify an existing Lyapunov candidate, $V(t,\bx)$, and just try to maximize/minimize $\rho(t)$. Or we can search for the parameters of $V(t,\bx)$, too. Again, we can initialize that search using the time-varying version of the Lyapunov equation, or the solutions to a time-varying LQR Riccati equation.

In practice, we often certify the Lyapunov conditions over $\bx$ at only a finite set of samples $t_i \in [t_1, t_2]$. I don't actually have anything against sampling in one dimension; there are no issues with scaling to higher dimensions, and one can make practical rigorous statement about bounding the sampling error. And in these systems, adding $t$ into all of the equation otherwise can dramatically increase the degree of the polynomials required for the SOS certificates. All of this was written up nicely in Tobenkin10b, and robust variants of it were developed in Majumdar16a.

Control design via alternations

For linear optimal control, we found controllers of the form, $\bu = -\bK \bx.$ The natural generalization of this to polynomial analysis will be to look for controllers of the form $\bu = \bK(\bx)$, where $\bK(\bx)$ is a polynomial vector. (For mechanical systems like the pendulum above, we might make $\bK(x)$ polynomial in $s$ and $c$.)

State feedback for linear systems

There is a clever reparameterization of the original problem that can make the co-design of Lyapunov functions and stabilizing controllers jointly convex in some cases.

Let's re-examine the Lyapunov conditions for linear systems from Eq. \ref{eq:lyap}, but now add in a linear state feedback $\bu = \bK\bx$, resulting in the closed-loop dynamics $\dot{\bx} = (\bA + \bB\bK)\bx.$ One can show that the set of all stabilizing $\bK$ can be characterized by the following two matrix inequalities Boyd94: \[{\bf P} = {\bf P}^T \succ 0, \qquad {\bf P}(\bA + \bB\bK) + (\bA + \bB\bK)^T{\bf P} \prec 0. \] Unfortunately, these inequalities are bilinear in the decision variables, ${\bf P}$ and $\bK$, and therefore non-convex. It turns out that we can turn these into linear matrix inequalities (LMIs) through a simple change of coordinates, $\bQ = {\bf P}^{-1}, {\bf Y} = \bK{\bf P}^{-1}$, resulting in \[ \bQ = \bQ^T \succ 0, \qquad \bA\bQ + \bQ\bA^T + \bB{\bf Y} + {\bf Y}^T \bB^T \prec 0. \] Furthermore, given matrices $\bA$ and $\bB$, there exists a matrix $\bK$ such that $(\bA + \bB\bK)$ is stable if and only if there exist matrices $\bQ$ and ${\bf Y}$ which satisfy this (strict) linear matrix inequality.

For global stability, there is a nice extension of this using sums of squares Prajna04b to the case of $$\dot{\bx} = {\bf f}_1(\bx){\bf m}(\bx) + {\bf f}_2(\bx)\bu,$$ where ${\bf m}(\bx)$ is a vector of monomials with ${\bf m}(\bx)=0$ iff $\bx=0$, and the elements of ${\bf f}_1(\bx)$ and ${\bf f}_2$ are polynomial. The search is over polynomial feedback controllers of the form $\bu = \bK(\bx){\bf m}(\bx).$

Control-Lyapunov Functions

Another core idea connecting Lyapunov analysis with control design is the "control-Lyapunov function". Given a system $\dot{\bx} = f(\bx, \bu)$, a control-Lyapunov function $V$ is a positive-definite function for which $$\forall \bx \neq 0, \exists \bu \quad \dot{V}(\bx,\bu) = \pd{V}{\bx} f(\bx, \bu) < 0 \quad \text{and} \quad \exists \bu \quad \dot{V}(0, \bu) = 0.$$ In words, for all $\bx$, there exists a control input that would allow the $V$ to decrease. Once again, we adorn this basic property with extra conditions (e.g. radially unbounded, or satisfied over a control- invariant set) in order to construct global or regional stability statements. What is important to understand is that we can design control-Lyapunov functions without explicitly parametrizing the controller; often the control action is not even determined until execution time by finding a particular $\bu$ that goes downhill.

Our sums-of-squares toolkit is well-suited for addressing questions with the $\forall$ quantifier over indeterminates. Working with the $\exists$ quantifier is much more difficult; we don't have an S-procedure-like solution for it. Interestingly, there is one approach that we've discussed above that effectively turns the $\exists$ into a $\forall$ -- that is the outer approximation approach to region of attraction analysis.

In the outer-approximation, we produce a barrier certificate to find the set of states where the controller cannot go (for any $\bu$). Our barrier certificate now has the form $$\forall \bx, \forall \bu, \quad \dot{\mathcal{B}}(\bx, \bu) \le 0.$$ Certainly, $\mathcal{B}(\bx) \ge 0$ is now an outer-approximation of the true "backward-reachable set" (BRS) of the fixed point Henrion14+Korda13. Again Posa17 has some nice examples of this written directly in sums-of-squares form.

Now here is where it gets a little frustrating. Certainly, sublevel sets of $\mathcal{B}(\bx)$ are control-invariant (via the proper choice of $\bu$). But we do not (cannot) expect that the entire estimated region (the 0-superlevel set) can be rendered invariant. The estimated region is an outer approximation of the backward reachable set. it's more than that... \dot{B}=0 happens somewhere in the interior, and it's conservative due to the forall u. Majumdar13b gave a recipe for extracting a polynomial control law from the BRS; an inner approximation of this control law can be obtained via the original SOS region of attraction tools. This is unfortunatley suboptimal/conservative. Although we would like to certify the control-Lyapunov conditions directly in an inner approximation, the $\exists$ quantifier remains as an obstacle.

Approximate dynamic programming with SOS

We have already established the most important connection between the HJB conditions and the Lyapunov conditions: $$\dot{J}^*(\bx) = -\ell(\bx,\bu^*) \qquad \text{vs} \qquad \dot{V}(\bx) \preceq 0.$$ The HJB involves solving a complex PDE; by changing this to an inequality, we relax the problem and make it amenable to convex optimization.

Sums-of-Squares Dynamic Programming

The linear programming approach works natively in discrete (or sampled) time/state/action settings. Sums-of-squares provides the a beautiful extension to continuous time/state/actions.

In the simple case, let us assume that $f(\bx,\bu)$ and $\ell(\bx,\bu)$ are polynomials (this assumption is discussed at length in the Lyapunov chapter). We can search for a polynomial $\hat{J}_\alpha(\bx)$ which satisfies the HJB inequality $$\forall \bx, \forall \bu, \quad 0 \le \ell(\bx,\bu) + \pd{\hat{J}_\alpha}{\bx} f(\bx,\bu),$$ using the following convex optimization: \begin{align*} \max_\alpha & \int_{\bf X} \hat{J}_\alpha(\bx) d\bx \\ \subjto \quad & \ell(\bx,\bu) + \pd{\hat{J}_\alpha}{\bx} f(\bx,\bu) \text{ is SOS }(\forall \bx, \forall \bu) \\ & \hat{J}_\alpha({\bf 0}) = 0. \end{align*} Since $\hat{J}_\alpha(\bx)$ is polynomial, the integral objective (over a bounded domain of interest) can be computed exactly, and is still linear in the coefficients $\alpha$. As in the linear programming approach, here the SOS constraint ensures that $\hat{J}$ is a lower bound on the cost to go, and the objective "pushes up" on this lower bound. As we increase the degree of the polynomial representing $\hat{J}$, we expect better and better approximations of the true cost-to-go.

Cubic polynomial optimal control

Consider a system with (scalar) dynamics: $$\dot{x} = f(x,u) = x - 4x^3 + u, u \in [-1, 1]$$ and the running cost $$\ell(\bx,\bu) = x^2 + u^2.$$ I've provided code that optimizes the cost-to-go over the interval $\bx \in [-1, 1].$

Pendulum swing-up

We can use this approach to design an approximation to the optimal controller for the (torque-limited) swing-up of the pendulum!

Cart-pole swing-up

This final example is adapted from Yang23. Even the lower-bound controller can effectively swing-up the cart-pole system!

I've only provided simple implementations above; their effectiveness is limited by the poor numerics of using polynomials with exceedingly high degree ($x^{20}$ is a very small number around $x=0$ and a very large number around $x=2\pi$). We should be able to improve the numerics dramatically by choosing more intelligent polynomial bases. I personally think this method has a lot of offer and it would be very interesting to see how far it could go with mature implementations.

Sampling Quotient-Ring Sum-of-Squares

Shen20

"Satisfiability modulo theories" (SMT)

Satisiability modulo theories (SMT). dReal is available in .

Mixed-integer programming (MIP) formulations

Continuation methods