How do you write a system of linear inequalities represented by the graph?

If you're seeing this message, it means we're having trouble loading external resources on our website.

If you're behind a web filter, please make sure that the domains *.kastatic.org and *.kasandbox.org are unblocked.

The solution of a linear inequality in two variables like Ax + By > C is an ordered pair (x, y) that produces a true statement when the values of x and y are substituted into the inequality.


Example

Is (1, 2) a solution to the inequality

$$2x+3y>1$$

$$2\cdot 1+3\cdot 2\overset{?}{>}1$$

$$2+5\overset{?}{>}1$$

$$7>1$$

The graph of an inequality in two variables is the set of points that represents all solutions to the inequality. A linear inequality divides the coordinate plane into two halves by a boundary line where one half represents the solutions of the inequality. The boundary line is dashed for > and < and solid for ≤ and ≥. The half-plane that is a solution to the inequality is usually shaded.

To graph a linear inequality in two variables (say, x and y ), first get y alone on one side. Then consider the related equation obtained by changing the inequality sign to an equality sign. The graph of this equation is a line.

If the inequality is strict ( < or > ), graph a dashed line. If the inequality is not strict ( ≤ or ≥ ), graph a solid line.

Finally, pick one point that is not on either line ( (0,0) is usually the easiest) and decide whether these coordinates satisfy the inequality or not. If they do, shade the half-plane containing that point. If they don't, shade the other half-plane.

Graph each of the inequalities in the system in a similar way. The solution of the system of inequalities is the intersection region of the solutions of the three inequalities.

Optimal control provides a powerful framework for formulating control problems using the language of optimization. But solving optimal control problems for nonlinear systems is hard! In many cases, we don't really care about finding the optimal controller, but would be satisfied with any controller that is guaranteed to accomplish the specified task. In many cases, we still formulate these problems using computational tools from optimization, and in this chapter we'll learn about tools that can provide guaranteed control solutions for systems that are beyond the complexity for which we can find the optimal feedback.

There are many excellent books on Lyapunov analysis; for instanceSlotine90is an excellent and very readable reference andKhalil01can provide a rigorous treatment. In this chapter I will summarize (without proof) some of the key theorems from Lyapunov analysis, but then will also introduce a number of numerical algorithms... many of which are new enough that they have not yet appeared in any mainstream textbooks.

Let's start with our favorite simple example.

Recall that the equations of motion of the damped simple pendulum are given by \[ ml^2 \ddot{\theta} + mgl\sin\theta = -b\dot{\theta}, \] which I've written with the damping on the right-hand side to remind us that it is an external torque that we've modeled.

These equations represent a simple second-order differential equation; in chapter 2 we discussed at some length what was known about the solutions to this differential equation--in practice we do not have a closed-form solution for $\theta(t)$ as a function of the initial conditions. Since we couldn't provide a solution analytically, in chapter 2 we resorted to a graphical analysis, and confirmed the intuition that there are fixed points in the system (at $\theta = k\pi$ for every integer $k$) and that the fixed points at $\theta = 2\pi k$ are asymptotically stable with a large basin of attraction. The graphical analysis gave us this intuition, but can we actually prove this stability property? In a way that might also work for much more complicated systems?

One route forward was from looking at the total system energy (kinetic + potential), which we can write down: \[ E(\theta,\dot{\theta}) = \frac{1}{2} ml^2\dot{\theta}^2 - mgl \cos\theta. \] Recall that the contours of this energy function are the orbits of the undamped pendulum.

A natural route to proving the stability of the downward fixed points is by arguing that energy decreases for the damped pendulum (with $b>0$) and so the system will eventually come to rest at the minimum energy, $E = -mgl$, which happens at $\theta=2\pi k$. Let's make that argument slightly more precise.

Evaluating the time derivative of the energy reveals \[ \frac{d}{dt} E = - b\dot\theta^2 \le 0. \] This is sufficient to demonstrate that the energy will never increase, but it doesn't actually prove that the energy will converge to the minimum when $b>0$ because there are multiple states(not only the minimum) for which $\dot{E}=0$. To take the last step, we must observe that set of states with $\dot\theta=0$ is not an invariant set; that if the system is in, for instance $\theta=\frac{\pi}{4}, \dot\theta=0$ that it will not stay there, because $\ddot\theta \neq 0$. And once it leaves that state, energy will decrease once again. In fact, the fixed points are the only subset the set of states where $\dot{E}=0$ which do form an invariant set. Therefore we can conclude that as $t\rightarrow \infty$, the system will indeed come to rest at a fixed point (though it could be any fixed point with an energy less than or equal to the initial energy in the system, $E(0)$).

This is an important example. It demonstrated that we could use a relatively simple function -- the total system energy -- to describe something about the long-term dynamics of the pendulum even though the actual trajectories of the system are (analytically) very complex. It also demonstrated one of the subtleties of using an energy-like function that is non-increasing (instead of strictly decreasing) to prove asymptotic stability.

Lyapunov functions generalize this notion of an energy function to more general systems, which might not be stable in the sense of some mechanical energy. If I can find any positive function, call it $V(\bx)$, that gets smaller over time as the system evolves, then I can potentially use $V$ to make a statement about the long-term behavior of the system. $V$ is called a Lyapunov function.

Recall that we defined three separate notions for stability of a fixed-point of a nonlinear system: stability i.s.L., asymptotic stability, and exponential stability. We can use Lyapunov functions to demonstrate each of these, in turn.

Given a system $\dot{\bx} = f(\bx)$, with $f$ continuous, and for some region ${\cal D}$ around the origin (specifically an open subset of $\mathbf{R}^n$ containing the origin), if I can produce a scalar, continuously-differentiable function $V(\bx)$, such that \begin{gather*} V(\bx) > 0, \forall \bx \in {\cal D} \setminus \{0\} \quad V(0) = 0, \text{ and} \\ \dot{V}(\bx) = \pd{V}{\bx} f(\bx) \le 0, \forall \bx \in {\cal D} \setminus \{0\} \quad \dot{V}(0) = 0, \end{gather*} then the origin $(\bx = 0)$ is stable in the sense of Lyapunov (i.s.L.). [Note: the notation $A \setminus B$ represents the set $A$ with the elements of $B$ removed.]

If, additionally, we have $$\dot{V}(\bx) = \pd{V}{\bx} f(\bx) < 0, \forall \bx \in {\cal D} \setminus \{0\},$$ then the origin is (locally) asymptotically stable. And if we have $$\dot{V}(\bx) = \pd{V}{\bx} f(\bx) \le -\alpha V(x), \forall \bx \in {\cal D} \setminus \{0\},$$ for some $\alpha>0,$ then the origin is (locally) exponentially stable.

Note that for the sequel we will use the notation $V \succ 0$ to denote a positive-definite function, meaning that $V(0)=0$ and $V(\bx)>0,$ for all $\bx \ne 0$. We can use "positive definite in ${\cal D}$" if this applies $\forall x \in {\cal D}$. Similarly, we'll use $V \succeq 0$ for positive semi-definite, $V \prec 0$ for negative-definite functions, etc.

The intuition here is exactly the same as for the energy argument we made in the pendulum example: since $\dot{V}(x)$ is always zero or negative, the value of $V(x)$ will only get smaller (or stay the same) as time progresses. Inside the subset ${\cal D}$, for every $\epsilon$-ball, I can choose a $\delta$ such that $|x(0)|^2 < \delta \Rightarrow |x(t)|^2 < \epsilon, \forall t$ by choosing $\delta$ sufficiently small so that the sublevel set of $V(x)$ for the largest value that $V(x)$ takes in the $\delta$ ball is completely contained in the $\epsilon$ ball. Since the value of $V$ can only get smaller (or stay constant) in time, this gives stability i.s.L.. If $\dot{V}$ is strictly negative away from the origin, then it must eventually get to the origin (asymptotic stability). The exponential condition is implied by the fact that $\forall t>0, V(\bx(t)) \le V(\bx(0)) e^{-\alpha t}.$

Notice that the system analyzed above, $\dot{\bx}=f(\bx)$, did not have any control inputs. Therefore, Lyapunov analysis is used to study either the passive dynamics of a system or the dynamics of a closed-loop system (system + control in feedback). We will see generalizations of the Lyapunov functions to input-output systems later in the text.

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.

Given a system $\dot{\bx} = f(\bx)$, with $f$ continuous, if I can produce a scalar, continuously-differentiable function $V(\bx)$, such that \begin{gather*} V(\bx) \succ 0, \\ \dot{V}(\bx) = \pd{V}{\bx} f(\bx) \prec 0, \text{ and} \\ V(\bx) \rightarrow \infty \text{ whenever } ||\bx||\rightarrow \infty,\end{gather*} then the origin $(\bx = 0)$ is globally asymptotically stable (G.A.S.).

If additionally we have that $$\dot{V}(\bx) \preceq -\alpha V(\bx),$$ for some $\alpha>0$, then the origin is globally exponentially stable.

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.

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.

Given a system $\dot{\bx} = f(\bx)$ with $f$ continuous. If we can produce a scalar function $V(\bx)$ with continuous derivatives for which we have $$V(\bx) \succ 0,\quad \dot{V}(\bx) \preceq 0,$$ and $V(\bx)\rightarrow \infty$ as $||\bx||\rightarrow \infty$, then $\bx$ will converge to the largest invariant set where $\dot{V}(\bx) = 0$.

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.

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.

Recall the example of using partial-feedback linearization to generate a swing-up controller for the cart-pole system. We first examined the dynamics of the pole (pendulum) only, by writing it's energy: $$E(\bx) = \frac{1}{2}\dot\theta^2 - \cos\theta,$$ desired energy, $E^d = 1$, and the difference $\tilde{E}(\bx) = E(\bx) - E^d.$ We were able to show that our proposed controller produced $$\dot{\tilde{E}} = -k \dot\theta^2 \cos^2\theta \tilde{E},$$ where $k$ is a positive gain that we get to choose. And we said that was good!

Now we have the tools to understand that really we have a Lyapunov function $$V(\bx) = \frac{1}{2}\tilde{E}^2(\bx),$$ and what we have shown is that $\dot{V} \le 0$. By LaSalle, we can only argue that the closed-loop system will converge to the largest invariant set, which here is the entire homoclinic orbit: $\tilde{E}(\bx) = 0$. We have to switch to the LQR controller in order to stabilize the upright.

Lyapunov function: $V(\bx) = \frac{1}{2}\tilde{E}^2(\bx).$Time-derivative of the Lyapunov function: $\dot{V}(\bx).$

As you can see from the plots, $\dot{V}(\bx)$ ends up being a quite non-trivial function! We'll develop the computational tools for verifying the Lyapunov/LaSalle conditions for systems of this complexity in the upcoming sections.

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$.

Include instability results, as in Briat15 Theorem 2.2.5

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^*.$

Given a system $\dot{\bx} = f(\bx)$ with $f$ continuous, if we can find a scalar function $V(\bx) \succ 0$ and a bounded sublevel set $${\cal G}: \{ \bx | V(\bx) \le \rho \}$$ on which $$\forall \bx \in {\cal G}, \dot{V}(\bx) \preceq 0,$$ then ${\cal G}$ is an invariant set. By LaSalle, $\bx$ will converge to the largest invariant subset of ${\cal G}$ on which $\dot{V}=0$.

Furthermore, if $\dot{V}(\bx) \prec 0$ in ${\cal G}$, then the origin is locally asymptotically stable and the set ${\cal G}$ is inside the region of attraction of this fixed point. Alternatively, if $\dot{V}(\bx) \preceq 0$ in ${\cal G}$ and $\bx = 0$ is the only invariant subset of ${\cal G}$ where $\dot{V}=0$, then the origin is asymptotically stable and the set ${\cal G}$ is inside the region of attraction of this fixed point.

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.

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).

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.

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.

robust quadrotor example

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), which means that understanding the regions where a particular controller can be guaranteed to work can be of critical importance.

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.

Consider a scalar polynomial, $p(\bx)$, and a semi-algebraic set $g(\bx) \preceq 0$, where $g$ is a vector of polynomials. If I can find a polynomial "multiplier", $\lambda(\bx)$, such that \[ p(\bx) + \lambda^T(\bx) g(\bx) \sos, \quad \text{and} \quad \lambda(\bx) \sos, \] then this is sufficient to demonstrate that $$p(\bx)\ge 0 \quad \forall \bx \in \{ \bx | g(\bx) \le 0 \}.$$ To convince yourself, observe that when $g(\bx) \le 0$, it is only harder to be positive, but when $g(\bx) > 0$, it is possible for the combined function to be SOS even if $p(\bx)$ is negative. We will sometimes find it convenient to use the short-hand: $$g(\bx) \le 0 \Rightarrow p(\bx) \ge 0$$ to denote the implication certified by the S-procedure (e.g. "whenever $g(\bx) \le 0$, we have $p(\bx) \ge 0$").

We can also handle equality constraints with only a minor modification -- we no longer require the multiplier to be positive. If I can find a polynomial "multiplier", $\lambda(\bx)$, such that \[p(\bx) + \lambda^T(\bx) g(\bx) \sos \] then this is sufficient to demonstrate that $$p(\bx)\ge 0 \quad \forall \bx \in \{ \bx | g(\bx) = 0 \}.$$ Here the intuition is that $\lambda(\bx)$ can add arbitrary positive terms to help me be SOS, but those terms contribute nothing precisely when $g(\bx)=0$.

The S-procedure gives us the tool we need to evaluate positivity only over a region of state space, which is precisely what we need to certify the Lyapunov conditions for a region-of-attraction analysis. Let us start with a positive-definite polynomial Lyapunov candidate, $V(\bx) \succ 0$, then we can write the Lyapunov conditions: $$\dot{V}(\bx) \prec 0 \quad \forall \bx \in \{ \bx | V(\bx) \le \rho \},$$ using sums-of-squares and the S-procedure: $$-\dot{V}(\bx) + \lambda(\bx)(V(\bx) - \rho) \text{ is SOS,} \quad \text{and} \quad \lambda(\bx) \text{ is SOS},$$ where $\lambda(\bx)$ is a multiplier polynomial with free coefficients that are to be solved for in the optimization.

I think it's easiest to see the details in an example.

Let's return to our example from above: \[ \dot{x} = -x + x^3 \] and try to use SOS optimization to demonstrate that the region of attraction of the fixed point at the origin is $x \in (-1,1)$, using the Lyapunov candidate $V = x^2.$

First, define the multiplier polynomial, \[ \lambda(x) = c_0 + c_1 x + c_2 x^2. \] Then define the optimization \begin{align*} \find_{\bf c} \quad & \\ \subjto \quad& - \dot{V}(x) + \lambda(x) (V(x)-1) \sos \\ & \lambda(x) \sos \end{align*}

You can try this example for yourself in.

In this example, we only verified that the one-sublevel set of the pre-specified Lyapunov candidate is negative (certifying the ROA that we already understood). Even more useful is if you are able to search for the largest $\rho$ that can satisfy these conditions. Unfortunately, in this first formulation, optimizing $\rho$ directly would make the optimization problem non-convex because we would have terms like $\rho c_0$, $\rho c_1 x$, ... which are bilinear in the decision variables; we need the sums-of-squares constraints to be only linear in the decision variables.

Fortunately, because the problem is convex with $\rho$ fixed (and therefore can be solved reliably), and $\rho$ is just a scalar, we can perform a simple line search on $\rho$ to find the largest value for which the convex optimization returns a feasible solution. This will be our estimate for the region of attraction.

There are a number of variations to this basic formulation; I will describe a few of them below. There are also important ideas like rescaling and degree-matching that can have a dramatic effect on the numerics of the problem, and potentially make them much better for the solvers. But you do not need to master them all in order to use the tools effectively.

In, we have packaged most of the work in setting up and solving the sums-of-squares optimization for regions of attraction into a single method RegionOfAttraction(system, context, options). This makes it as simple as, for instance:

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

Remember that although we have tried to make it convenient to call these functions, they are not a black box. I highly recommend opening up the code implementing the RegionOfAttraction method and understanding how it works (even if you don't love C++, the MathematicalProgram modeling language in Drake makes these formulations pretty readable). There are lots of different options / formulations, and numerous numerical recipes to improve the numerics of the optimization problem.

In order to develop our intuition for these tools, it is helpful to understand how "tight" the inner approximation can be for each choice of the fixed degree of the Lyapunov function / Lagrange multipliers. I've included a few of my favorite examples of nonlinear systems which use some clever construction to have known regions of attraction.

We will study the Van der Pol oscillator in more detail when we study limit cycles. For now, it suffices to know that the dynamics are polynomial, and that the stable limit cycle of the system when simulated forward in time produces a known (via numerical integration) boundary for the region of attraction of the fixed point at the origin when we simulate the system backwards in time.

Figure here

I've provided an example here for solving this problem with the RegionOfAttraction in Drake, but we will also work you through writing the sums-of-squares problem yourself in the exercises.

Here is one important variation for finding the level set of a candidate region of attraction, which turns an inequality in the S-procedure into an equality. This formulation is jointly convex in $\lambda(\bx)$ and $\rho$, so one can optimize them in a single convex optimization. It appeared informally in an example inParrilo00, and was discussed a bit more inShen20.

Under the assumption that the Hessian of $\dot{V}(\bx)$ is negative-definite at the origin (which is easily checked), we can write \begin{align*} \max_{\rho, \lambda(\bx)} & \quad \rho \\ \subjto & \quad (\bx^T\bx)^d (V(\bx) - \rho) + \lambda(\bx) \dot{V}(\bx) \text{ is SOS},\end{align*} with $d$ a fixed positive integer. As you can see, $\rho$ no longer multiplies the coefficient of $\lambda(\bx)$. But why does this certify a region of attraction?

You can read the sums-of-squares constraint as certifying the implication that whenever $\dot{V}(x) = 0$, we have that either $V(\bx) \ge \rho$ OR $\bx = 0$. Multiplying by some multiple of $\bx^T\bx$ is just a clever way to handle the "or $\bx=0$" case, which is necessary since we expect $\dot{V}(0) = 0$. This implication is sufficient to prove that $\dot{V}(\bx) \le 0$ whenever $V(\bx) \le \rho$, since $V$ and $\dot{V}$ are smooth polynomials; we examine this claim in full detail in one of the exercises.

Using the S-procedure with equalities instead of inequalities also has the potential advantage of removing the SOS constraint on $\lambda(\bx).$ But perhaps the biggest advantage of this formulation is the possibility of dramatically simplifying the problem using the quotient ring of this algebraic variety, and in particular some recent results for exact certification using sampling varietiesShen20.

The machinery so far has used optimization to find the largest region of attraction that can be certified given a candidate Lyapunov function. This is not necessarily a bad assumption. For most stable fixed points, we can certify the local stability with a linear analysis, and this linear analysis gives us a candidate quadratic Lyapunov function that can be used for nonlinear analysis over a region. Practically speaking, when we solve an LQR problem, the cost-to-go from LQR is a good candidate Lyapunov function for the nonlinear system. If we are simply analyzing an existing system, then we can obtain this candidate by solving a Lyapunov equation (Eq \ref{eq:algebraic_lyapunov}).

But what if we believe the system is regionally stable, despite having an indefinite linearization? Or perhaps we can certify a larger volume of state space by using a Lyapunov candidate with degree greater than two. Can we use sums-of-squares optimization to find that in the region of attraction case, too?

To accomplish this, we will now make $V(x)$ a polynomial (of some fixed degree) with the coefficients as decision variables. First we will need to add constraints to ensure that $$V(0) = 0 \quad {and} \quad V(\bx) - \epsilon \bx^T\bx \text{ is SOS},$$ where $\epsilon$ is some small positive constant. This $\epsilon$ term simply ensures that $V$ is strictly positive definite. Now let's consider our basic formulation: $$-\dot{V}(\bx) + \lambda(\bx)(V(\bx) - 1) \text{ is SOS,} \quad \text{and} \quad \lambda(\bx) \text{ is SOS}.$$ Notice that I've replaced $\rho=1$; now that we are searching for $V$ the scale of the level-set can be set aribtrarily to 1. In fact, it's better to do so -- if we did not set $V(0)=0$ and $V(\bx) \le 1$ as the sublevel set, then the optimization problem would be underconstrained, and might cause problems for the solver.

Unfortunately, the derivative constraint is now nonconvex (bilinear) in the decision variables, since we are searching for both the coefficients of $\lambda$ and $V$, and they are multiplied together. Our equality-constrained formulation doesn't get us around this one, either (since the coefficients of $V$ also appear in $\dot{V}$). Here we have to resort to some weaker form of optimization. In practice, we have had good practical success using bilinear alternations: start with an initial $V$ (e.g. from LQR), and search for $\lambda$; then fix $\lambda$ and search for $V$ and repeat until convergence (see, for instanceTedrake10+Majumdar16a). A typical choice for the objective is to maximize some convex surrogate for the volume of the certified region of attraction in state space; typically we use the determinant of the quadratic form describing a contained ellipse.

Barring numerical issues, this algorithm is guaranteed to have recursive feasibility and tends to converge in just a few alternations, but there is no guarantee that it finds the optimal solution.

(Details coming soon...)

Although we are using convex optimization to find regions of attraction, this does not imply that the regions we are certifying need to be convex. To demonstrate that, let's make a system with a known, non-convex region of attraction. We'll do this by taking some interesting potential function $U(x) \in SOS$ and setting the dynamics to be $\dot{x} = (U(x)-1) \frac{\partial U}{\partial x}^T$, which has $U(x) <= 1$ as the region of attraction.

Slightly more general is to write $\dot{x} = (U(x)-1) {\bf R}(\theta) \frac{\partial U}{\partial x}^T$, where ${\bf R}(\theta) = \begin{bmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta\end{bmatrix}$ is the 2D rotation matrix, and $\theta<\pi$ is a constant parameter (not a decision variable nor indeterminate), which still has the same region of attraction.

Searching for Lyapunov functions of only degree 2 (quadratics) will, of course, limit our certified regions to convex regions. But if we increase the degree of the Lyapunov function to degree 4, then we can tighten our inner approximation, and indeed certify non-convex regions. While we can use the Lyapunov method for linear systems to initialize quadratic Lyapunov functions, the ability to search for the parameters of the Lyapunov functions really shines when the Lyapunov candidates are higher degree (such as the quartic we use here).

Figures here

So far we have only considered verifying regions of attraction for systems of the form $\dot{\bx} = f(\bx)$ (e.g. with no control input); the implication has been that one should design a feedback controller first and then analyze the stability of the closed-loop dynamics. But some control design techniques produce artifacts that can be useful for this verification. Often jointly designing a controller and certificate can be more tractable than trying to certify an arbitrary controller.

One nice example of this that we've already seen is the linear quadratic regulator. Even when it is applied to (the linearization of) nonlinear systems, it still produces a quadratic cost-to-go function $\bar\bx^T {\bf S} \bar\bx$ which can be immediately used as a candidate Lyapunov function for the nonlinear system.

To demonstrate this...

Numerical example of LQR for a polynomial systems

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 fromLasserre10also 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 use a very related set of conditions to 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) > 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) > 0.$$

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$Folland01and is linear in the coefficients. More sophisticated alternatives also existHenrion09.

Let's revisit my favorite simple 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) \ge 0.\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, for $\mathcal{B}(0) \ge 0.1$, 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$).

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 (seeWampler11).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,Wampler11andSommese05. 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!

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 inwe 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):

How do you write a system of linear inequalities represented by the graph?

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.

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".

Given a linear system in the form ${\bf E}\dot\bx = \bA\bx,$ the origin is asymptotically stable if there exists a matrix ${\bf P}$ such that $${\bf P} = {\bf P}^T \succ 0, \quad \bA^T {\bf P}{\bf E} + {\bf E}^T {\bf P} \bA \prec 0.$$ Furthermore, if such a ${\bf P}$ exists, then we can obtain one by solving the generalized Lyapunov equation: $$\bA^T {\bf P}{\bf E} + {\bf E}^T {\bf P} \bA = -{\bf Q},$$ for any ${\bf Q} \succ 0$Stykel02.

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).$

Flesh this out and add an example. $V(\bx) = \bx^T {\bf E}^T(\bx) {\bf P}{\bf E}(\bx) \bx$ seems a natural choice of Lyapunov candidates. V = v'M(q)'PM(q)v + q'Qq could bea natural form for manipulators (do I need q-v cross terms?); it looks a bit like kinetic energy...

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).

There are a few things that do break this clean polynomial view of the world. Rotary springs, for instance, if modeled as $\tau = k (\theta_0 - \theta)$ will mean that $\theta$ appears alongside $\sin\theta$ and $\cos\theta$, and the relationship between $\theta$ and $\sin\theta$ is sadly not polynomial. As we have presented it so far, even linear feedback from LQR actually looks like a linear spring.

Looking more closely, though, there is a natural alternative. Consider the mechanical energy of the simple pendulum , $E = \frac{1}{2}ml^2 \dot\theta^2 + mgl(1- \cos\theta),$ where we've taken the potential energy so that $E \ge 0$). Observe that using half-angle formulas we can write this as a quadratic form: $$E = \frac{1}{2} \begin{bmatrix} 2\sin(\frac{\theta}{2}) \\ \dot\theta \end{bmatrix}^T \begin{bmatrix} mgl & 0 \\ 0 & ml^2 \end{bmatrix} \begin{bmatrix} 2\sin(\frac{\theta}{2}) \\ \dot\theta \end{bmatrix}.$$ The term $2 \sin(\frac{\theta}{2})$ has slope one at $\theta = 0$. This isn't unique to the pendulum, $2\sin^2(\frac{\theta}{2})$ is a natural choice for a Lyapunov candidate in the trigonometric coordinates that is periodic in $2\pi.$ In the case of LQR applied to the linearization of a nonlinear mechanical system, one could replace $-k \theta$ with $-2 k \sin(\frac{\theta}{2})$; after all their linearization is the same.

example of design + certification of LQR using $u = -\bK \sin\theta$

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:

How do you write a system of linear inequalities represented by the graph?
The stereographic projection.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.

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

Details coming soon. For some recent work with a fairly extensive bibliography, seeDai21.

Control contraction metrics

other topics/ideas: verifying neural network control relation to adversarial examples. make the point that sampling doesn't scale (even with polynomials). "almost lyapunov functions" will happen when we talk about robust control. hopfield enegy function. scaling (e.g. quadratic in high dimensions with just one strictly negative eigenvalue); would be very hard to verify with samples.

For the system \begin{align*} \dot x_1 &=-\frac{6x_1}{(1+x_1^2)^2}+2x_2, \\ \dot x_2 &=-\frac{2(x_1+x_2)}{(1+x_1^2)^2}, \end{align*} you are given the positive definite function $V(\bx) =\frac{x_1^2}{1 + x_1^2}+ x_2^2$ (plotted here) and told that, for this system, $\dot V(\bx)$ is negative definite over the entire space. Is $V(\bx)$ a valid Lyapunov function to prove global asymptotic stability of the origin for the system described by the equations above? Motivate your answer.

You are given a dynamical system $\dot \bx = f(\bx)$, with $f$ continuous, which has a fixed point at the origin. Let $B_r$ be a ball of (finite) radius $r > 0$ centered at the origin: $B_r = \{ \bx : \| \bx \| \leq r \}$. Assume you found a continuously-differentiable scalar function $V(\bx)$ such that: $V(0) = 0$, $V(\bx) > 0$ for all $\bx \neq 0$ in $B_r$, and $\dot V(\bx) < 0$ for all $\bx \neq 0$ in $B_r$. Determine whether the following statements are true or false. Briefly justify your answer.

  1. $B_r$ is an invariant set for the given system, i.e.: if the initial state $\bx(0)$ lies in $B_r$, then $\bx(t)$ will belong to $B_r$ for all $t \geq 0$.
  2. $B_r$ is a subset of the ROA of the fixed point $\bx = 0$, i.e.: if $\bx(0)$ lies in $B_r$, then $\lim_{t \rightarrow \infty} \bx(t) = 0$.

If $V_1(\bx)$ and $V_2(\bx)$ are valid Lyapunov functions that prove global asymptotic stability of the origin, does $V_1(\bx)$ necessarily equal $V_2(\bx)$?

Consider the system given by \begin{align*} \dot x_1 &= x_2 - x_1^3, \\ \dot x_2 &= - x_1 - x_2^3. \end{align*} Show that the Lyapunov function $V(\bx) = x_1^2 + x_2^2$ proves global asymptotic stability of the origin for this system.

We can use Lyapunov analysis to analyze the behavior of optimization algorithms like gradient descent. Consider an objective function $\ell:\mathbf{R}^n\rightarrow \mathbf{R}$, and we want to minimize $\ell(x)$. The gradient flow is a continuous-time analog of gradient descent and is defined as $\dot{x} = -\frac{\partial \ell}{\partial x}$.

  1. Show that the objective function $\ell(x)-\ell(x^*)$ is a Lyapunov function of the gradient flow where $x^*$ is a unique minimizer.
  2. We can use Lyapunov to argue that an optimization problem will converge to a global optimum, even if it is non-convex. Suppose that the Lyapunov function $\ell$, has negative definite $\dot{\ell}$. Show that the objective function $\ell$ has a unique minimizer at the origin.
  3. Consider the objective function in the figure below. Could this be a valid Lyapunov function for showing global asymptotic stability?
  4. How do you write a system of linear inequalities represented by the graph?
    Lyapunov function candidate for asymptotic stability.
  5. Consider the objective function $\ell(x) = x_1^4-8x_1^3+18x_1^2+x_2^2$ with $x\in\mathbf{R}^2$. Find the largest $\rho$ such that $\{x\mid \ell(x)<\rho \}$ is a valid region of attraction for the origin.

In this exercise, we examine the idea of a control-Lyapunov function to drive a wheeled robot, implementing the controller proposed inAicardi95.

Similar to this previous example, we use a kinematic model of the robot. We represent with $z_1$ and $z_2$ its Cartesian position and with $z_3$ its orientation. The controls are the linear $u_1$ and angular $u_2$ velocities. The equations of motion read \begin{align*} \dot z_1 &= u_1 \cos z_3, \\ \dot z_2 &= u_1 \sin z_3, \\ \dot z_3 &= u_2.\end{align*} The goal is to design a feedback law $\pi(\bz)$ that drives the robot to the origin $\bz=0$ from any initial condition. As pointed out inAicardi95, this problem becomes dramatically easier if we analyze it in polar coordinates. As depicted below, we let $x_1$ be the radial and $x_2$ the angular coordinate of the robot, and we define $x_3 = x_2 - z_3$. Analyzing the figure, basic kinematic considerations lead to \begin{align*} \dot x_1 &= u_1 \cos x_3, \\ \dot x_2 &= - \frac{u_1 \sin x_3}{x_1}, \\ \dot x_3 &= - \frac{u_1 \sin x_3}{x_1} - u_2.\end{align*}

Wheeled robot with Cartesian $\bz$ and polar $\bx$ coordinate system.
  1. For the candidate Lyapunov function $V(\bx) = V_1(x_1) + V_2(x_2, x_3)$, with $V_1(x_1) = \frac{1}{2} x_1^2$ and $V_2(x_2, x_3) = \frac{1}{2}(x_2^2 + x_3^2)$, compute the time derivatives $\dot V_1 (\bx, u_1)$ and $\dot V_2(\bx, \bu)$.
  2. Show that the choice \begin{align*} u_1 &= \pi_1(\bx) = - x_1 \cos x_3, \\ u_2 &= \pi_2(\bx) = x_3 + \frac{(x_2 + x_3) \cos x_3 \sin x_3}{x_3}, \end{align*} makes $\dot V_1 (\bx, \pi_1(\bx)) \leq 0$ and $\dot V_2 (\bx, \pi(\bx)) \leq 0$ for all $\bx$. (Technically speaking, $\pi_2(\bx)$ is not defined for $x_3=0$. In this case, we let $\pi_2(\bx)$ assume its limiting value $x_2 + 2 x_3$, ensuring continuity of the feedback law.)
  3. Explain why Lyapunov's direct method does not allow us to establish asymptotic stability of the closed-loop system.
  4. Substitute the control law $\bu = \pi (\bx)$ in the equations of motion, and derive the closed-loop dynamics $\dot \bx = f(\bx, \pi(\bx))$. Use LaSalle's theorem to show (global) asymptotic stability of the closed-loop system.
  5. In we set up a simulation environment for you to try the controller we just derived. Type the control law from point (b) in the dedicated cell, and use the notebook plot to check your work.
  1. Are there positive definite functions that are not representable as sums of squares?
  2. If a fixed point of our dynamical system does not admit a SOS Lyapunov function, what can we conclude about its stability?

In this exercise you will use SOS optimization to approximate the ROA of the time-reversed Van der Pol oscillator (a variation of the classical Van der Pol oscillator which evolves backwards in time). In , you are asked to test the following SOS formulations.

  1. The one from the example above, augmented with a line search that maximizes the area of the ROA.
  2. A single-shot SOS program that can directly maximize the area of the ROA, without any line search.
  3. An improved version of the previous, where less SOS constraints are imposed in the optimization problem.

In this exercise we carefully investigate the procedure for designing stabilizing controller via alternations. Consider the problem of jointly designing a polynomial control function $u = K(x)$ and Lyapunov certificate $V(x)$ for the dynamic system $$\dot{x} = f_{1}(x) + f_{2}(x)K(x)$$ We can write the search for the maximal region of attraction as: \begin{align*} \max_{\rho, K(x), V(x), \lambda(x)}& \qquad \rho \\ \text{subject to } & \qquad V(x),~\lambda(x) \text{ are SOS} \\ & -\frac{\partial V}{\partial x} (f_{1}(x) + f_{2}(x)K(x)) + \lambda(x)( V(x) - \rho) \text{ is SOS} \end{align*} Assume we know some $\rho_{0},~ V_{0}(x),~ K_{0}(x)$, and $\lambda_{0}(x)$ which are feasible for the above program.

How do you write a system of linear inequalities by graphing?

SOLVE A SYSTEM OF LINEAR INEQUALITIES BY GRAPHING..
Graph the first inequality. Graph the boundary line. ... .
On the same grid, graph the second inequality. Graph the boundary line. ... .
The solution is the region where the shading overlaps..
Check by choosing a test point..