Neural Lyapunov Redesign

Studying the stability region of autonomous systems and designing controllers (policies) to drive a non-autonomous system towards a target behavior are of fundamental importance in many disciplines, such as aviations (Liao and Wang, 2004), autonomous driving (Wen-Xing and Li-Dong, 2018), and robotics (Pierson and Gashler, 2017). An indisputable goal for a controller is to stabilize the system. Unlike the global nature of linear systems, stability is a local property in nonlinear systems. Knowledge of the stability region is essential in many applications, such as stability of power systems (Xin et al., 2007), design of associative memory in artificial neural networks (Hopfield, 1994), robotics (Westervelt et al., 2003), and biology (Baer et al., 2006).

Controllers that enhance the stability region of a system, also known as Region of Attraction (RoA), are highly desired as they make more clever use of the inherent nonlinear structure of the system. For example, an autonomous driving system will remain safe under more diverse and potentially harsh conditions (Imani Masouleh and Limebeer, 2018).

A great number of methods for designing a controller for nonlinear systems and determining their RoA (Isidori, 2014; Khalil and Grizzle, 2002) have been proposed in the literature mainly based on Lyapunov’s theory of stability (Liapounoff, 1907). However, the design of a stabilizing controller has been ad-hoc for every class of nonlinear systems. Our work is inspired by a classic method called Lyapunov Redesign in control theory that uses a given Lyapunov function to design a controller such that the closed-loop system becomes stable when assessed by that Lyapunov function. Here, we relax the necessity to know the Lyapunov function by learning it together with the controller.

Machine learning and control theory have contributed to each other in various ways. Function approximators such as kernel methods (Pillonetto et al., 2014) and neural networks (Chu et al., 1990) have been successfully used for system identification. Moreover, ideas from information theory and statistics have been adopted for model selection in system identification (Solowjow et al., 2018; Mehrjou et al., 2018). Moreover, recent advances in implicit generative models have proven effective in filtering applications where the states of the system are inferred from noisy observations (Mehrjou and Schölkopf, 2018a). The development of automatic differentiation packages for machine learning has opened up a new approach to stability analysis and estimation of functions of interest in control theory such as Hamiltonian Greydanus et al. (2019) or Lyapunov function Mehrjou and Schölkopf (2019).

The analytical tools from control theory have helped to analyze the dynamics of machine learning algorithms. For example, the algorithms such as generative adversarial networks (Goodfellow et al., 2014) and $n$-player cooperative games have been analyzed as dynamical systems to obtain some insight into the dynamics of learning (Mehrjou and Schölkopf, 2018b; Mehrjou, 2018; Mescheder et al., 2017; Mehrjou et al., 2019) and the emergence of the playing strategies (Balduzzi et al., 2018). In the present work, analytical tools from control theory are employed to study the strength of the learning signal during the course of training.

The popularity of multi-player learning algorithms in machine learning after Generative Adversarial Networks (Goodfellow et al., 2014) have been inspired in areas of control theory where such multi-player setting can be defined. Mehrjou et al. (2020) have proposed a two-player setting where one player learns the Lyapunov function and the region of attraction while the other player learns the governing equations of the system. These players inform each other and their collaboration gives a faster convergence compared with learning each of them separately. The present work can be seen as a two-player algorithm where one player learns the Lyapunov function and the stability region while the other player improves the controller. The information exchange between the two players is the key factor of the proposed algorithm.

As the Lyapunov function is one of the most important objects in control theory, machine learning tools have been employed to learn it for nonlinear systems (Richards et al., 2018). A given Lyapunov function is used in (Berkenkamp et al., 2016; Berkenkamp et al., 2017) to derive a controller with a safety guarantee. Chang et al. (2019) solve an optimization problem at every stage to find the states that violate the Lyapunov condition. Our work proposes improvements in various directions. We improve the Lyapunov learning algorithm of Richards et al. (2018) by a theoretically motivated method and show that hyperbolicity is a needed feature for nonlinear systems whose controllers are learned iteratively. Moreover, Berkenkamp et al. (2016); Berkenkamp et al. (2017) assume the Lyapunov function is given while our work co-learns the Lyapunov function and the controller. Chang et al. (2019)’s algorithm needs to be solved until the end that can be too costly as a global optimization problem has to be solved multiple times. However, our work is a growing method that improves the controller while the system is in action.

We bring together various tools from control theory and machine learning to build an iterative algorithm that redesigns both controller and the Lyapunov function to enlarge the stability region. Our contributions are as follows: 1) Improving the learning of the Lyapunov function, 2) Interlacing the Lyapunov learning with the policy update to enlarge the RoA iteratively, and 3) Providing theoretical results for the tractable class of systems and analyzing the learning signal. Section 2 goes over the preliminary materials and Section 3 describe the problem. The proposed algorithm and its theoretical discussions are presented in Section 4. Finally, empirical evaluation comes in Section 5, followed by related work and conclusion in Section 6. To stay within the space limit, some proofs, theoretical discussions, and extended experimental results are provided in the Appendix.

Here we present the definitions, notations, and theoretical results that are used in later sections.

Both $W^{s}$ and $W^{u}$ are proved to be invariant sets (Palis and De Melo, 2012).

We consider a discrete-time TIA system as in 1, where the control signal is produced by a feedback-controller $\pi(\cdot;\psi):\mathcal{X}\to\mathcal{U}$ parameterized by $\psi$. Therefore, the closed-loop dynamics denoted by $f_{\pi}$ is a functional of the controller and is consequently parameterized by $\psi$ as $\mathbf{x}_{k+1}=f(\mathbf{x}_{k},\pi(\mathbf{x}_{k};\psi))=f_{\pi}(\mathbf{x}_{k};\psi)$. Without loss of generality, we assume that the equilibrium point of interest is located at the origin $\bar{\mathbf{x}}=\mathbf{0}$. The policy $\pi$ induces a RoA around the equilibrium point denoted by $\mathcal{R}_{\pi}$.²² 2 We drop $\bar{\mathbf{x}}$ from the superscript of $\mathcal{R}_{\pi}^{\bar{\mathbf{x}}}$, since we always assume $\bar{\mathbf{x}}=\mathbf{0}$, unless otherwise stated.

Each control task can be broken down into two subtasks: 1) Controller synthesis and 2) Closed-loop response evaluation. The controller is designed to optimize some measure of performance. In this work, the performance measure is the size of the stability region.

We endow the state space $\mathcal{X}$ with a measure $\mu$ to obtain the measure space $(\mathcal{X},\mathcal{B}(\mathcal{X}),\mu)$ with Borel sigma-algebra $\mathcal{B}(\mathcal{X})$. To prevent pathological cases, we assume $\mathcal{X}$ to be compact with $\mu(\mathcal{X})<\mu_{\infty}<\infty$. Let $\Pi=\{\pi:\mathcal{X}\to\mathcal{U}:\pi\in C^{1}\;\text{and bounded}\}$ be the set of all functions from which the policy is chosen. The goal is to find a member $\pi^{*}$ of the equivalence class of optimal policies $\Pi^{*}\subseteq\Pi$, where $\Pi^{*}$ is defined as $\Pi^{*}=\{\pi^{*}\in\Pi:\mu(\mathcal{R}_{\pi^{*}})=\max_{\pi\in\Pi}\mu(\mathcal{R}_{\pi})\}$.

The main challenge in this optimization problem is the fact that there is no analytical or straightforward way to infer how changing $\pi$ changes $\mu(\mathcal{R}_{\pi})$. If there exists a differentiable map from $\pi$ to $\mu(\mathcal{R}_{\pi})$, one could locally increase the stability region by perturbing the policy in the direction of $\partial\mu(\mathcal{R}_{\pi})/\partial\pi$. However, except for extremely simple systems, such a map cannot be derived analytically. In this work, we construct a bridge between these two objects by an auxiliary function, which is an evolving Lyapunov function that is learned alongside the policy.

Our goal is to construct a sequence of policies $(\pi_{1},\pi_{2},\ldots)$ that gives a sequence of RoAs $(\mathcal{R}_{\pi_{1}},\mathcal{R}_{\pi_{2}},\ldots)$, such that $\mathcal{R}_{\pi_{n}}\xrightarrow{\mu}\bar{\mathcal{R}}$ as $n\to\infty$, where $\bar{\mathcal{R}}\subseteq\mathcal{D}$ is the largest achievable RoA that is constrained by the physical limitations of the system (see Appendix A for the characterization of $\bar{\mathcal{R}}$ using the concept of control Lyapunov function). To achieve this goal, we need to address two challenges: 1) Approximating $\mathcal{R}_{\pi_{n}}$ for a fixed $\pi_{n}$ and 2) Using $(\mathcal{R}_{\pi_{n}},f_{\pi_{n}})$ to find $\pi_{n+1}$. Next section, explains our proposed method to address these two challenges.

For a function $g$, let’s define its sublevel set with level value $a$ as ${\mathcal{S}}_{a}(g):=\{\mathbf{x}\in\mathbb{R}^{d}:g(\mathbf{x})<a\}$. The index $n\in\mathbb{N}\cup\{0\}$ refers to a phase of the algorithm. At phase $n$, let $\mathcal{R}_{\pi_{n}}$ be the RoA of the closed-loop system 1 that is induced by the state-feedback policy $\pi_{n}$. It can be shown that there exists an optimal Lyapunov function $V_{\pi_{n}}$ with a level value $c_{n}$, such that $\mathcal{R}_{\pi_{n}}=S_{c_{n}}(V_{\pi_{n}})$ (Vannelli and Vidyasagar, 1985). Therefore, the information of $\mathcal{R}_{\pi_{n}}$ is encoded in $(V_{\pi_{n}},c_{n})$.

Starting with a conservative controller (e.g., a quadratic controller for the linearized system), our method inductively constructs a sequence of policies $(\pi_{1},\pi_{2},\ldots)$ that eventually converges to a policy with maximal achievable RoA. See Figure 1(Left) for an overall schematic of the algorithm.

Each step of this inductive process is called a phase of the algorithm. Each phase consists of two sub-phases: 1) Learning the Lyapunov function and the RoA corresponding to the policy of that phase 2) Updating the policy to enlarge the RoA. The RoA estimation sub-phase finds a Lyapunov function $V_{\pi_{n}}$ and level value $c_{n}$, such that ${\mathcal{S}}_{c_{n}}(V_{\pi_{n}})=\mathcal{R}_{\pi_{n}}$. Then, the policy update sub-phase learns a new policy $\pi_{n+1}$, such that ${\mathcal{S}}_{c_{n}}(V_{\pi_{n}})\subseteq\mathcal{R}_{\pi_{n+1}}$. These conditions need to be satisfied for every $n$ that consequently limits the class of treatable systems. Before delving into the algorithmic implementation, Section 4.1 provides necesary theoretical insights into this favorable class of systems along with other theoretical considertations for the proposed multi-phase growing algorithm.

We consider an inverted pendulum system defined as $\dot{\theta}=\omega$ and $\dot{\omega}=\frac{g}{l}\sin(\theta)+\frac{u}{I}-\mu_{f}\frac{\omega}{I}$ where the state vector $\mathbf{x}=(\theta,\omega)$ consists of the angle and angular velocity. Moreover $(g=0.81,l=0.5,I={\rm mass}\times l^{2}=0.25,\mu_{f}=0)$ are the acceleration of gravity, length, inertia, and friction coefficient. The scalar $u$ is the input force. The open-loop system (with $u=0$) has equilibrium points at $(\theta,\omega)=(k\pi,0)$ with $k\in\mathbb{Z}$. We focus on the equilibrium point $(0,0)$ in the frictionless setting where the system shows oscillatory behavior and consequently is not asymptotically stable (see Section G.2 for system and modeling details). First, an LQR controller $K$ is designed for the linearized system around the origin (see the vector field and the initial RoA in Figure 5 in the Appendix). The control signal passes through a loose saturation function as $u=\pi_{0}(\mathbf{x};\psi)={\rm SAT}_{\psi}([\theta,\omega]^{\mkern-1.5mu\mathsf{T}}K)$. The function ${\rm SAT}$ is parameterized by $\psi=(a,b,m_{a},m_{b})$ as illustrated in Figure 4. In the first experiment, the slopes $m_{a}=0$ and $m_{b}=0$ are kept fixed where $a$ and $b$ are trainable parameters of the policy, i.e., $\psi=(a,b)$. The Lyapunov function $V(\cdot;\theta)$ is realized by a $3-$layer neural network parameterized by $\theta$. Each layer has $64$ neurons with a special architecture (see Section G.1) inspired by (Richards et al., 2018) followed by tanh activation function that imposes the positive definiteness of the entire network as is required by the Lyapunov conditions in 2. The chosen hyper-parameters for the SGD training can be found in Section G.3.

The initial controller gives a small RoA since it is designed for the local linear approximation of the system. As the initial policy is LQR designed for the locally linearized model, $V(\cdot;\theta)$ is pre-trained by the quadratic function $0.1\theta^{2}+0.1\omega^{2}$.

After pre-training, sub-phases of the algorithm of Section 4.2 are run alternately to capture the RoA and improve the policy. The green step-like plot in Figure 2(a) and Figure 2(b) shows the true size of the RoA. Each jump in the green plots shows one iteration of the policy update algorithm resulting in an increased RoA. The fluctuating red plot in Figure 2(a) shows the size of the estimated RoA without our improvement over the RoA estimation algorithm of Richards et al. (2018) while Figure 2(b) shows the outcome of the presence of our proposed additional term in the loss function 5. It shows that the added term results in a less fluctuating estimate of the RoA, and when combined with the policy update sub-phase, gives a faster convergence to a larger RoA ($35.68\%\text{ vs }27.50\%$ fraction of the domain volume after $7$ policy updates).

As stated in Section 4.2, $\bar{c}=1$ is not necessarily equal to $c_{n}$. After learning $V_{n}$ with $\bar{c}$ in 5, the algorithm searches for a value of $c_{n}$ such that the Lyapunov decrease condition is met for all states within the sublevel set ${\mathcal{S}}_{c_{n}}(V_{n})$. It can be seen in Figure 2(c) that these values converge to $\bar{c}=1$ that can be perceived as a sign of the stable training of the algorithm.

The trace of the parameters of the policy is shown in Figure 2(d). As these policy’s trainable parameters represent the upper and lower limits of the loose threshold function, the policy learning algorithm updates them in the directions that decrease their suppressing effect. This is what we also expect from the physics of the system. Graphical visualization of the policy update and RoA estimation sub-phase is shown in Figure 1(Right). Each row shows one phase of the algorithm. The policies are updated along the rows from top to bottom. Within one row, the policy is fixed and RoA is estimated from left to right (See Figure 7 in the appendix for a larger visualization).

In the second experiment, the threshold limits $a=-0.2$ and $b=0.2$ are kept fixed while the slopes $\psi=(m_{a},m_{b})$ are trainable parameters. The rest of the training setting remains the same as the previous experiment. Figures 3(a) and 3(b) shows that the policy update sub-phase enlarges the RoA (green plot) of the system while the RoA estimation sub-phase (red plot) manages to follow the new RoA after each policy update. Our improved RoA estimation algorithm results in a more monotonic convergence of the estimated RoA that ultimately learns a controller that induces a larger RoA ($24.86\%\text{ vs }22.11\%$). Similar to Figure 2(c), convergence of $c_{n}$ values to $\bar{c}=1$ can be seen in Figure 3(c). The trace of the upper and lower slopes are shown in Figure 3(d). Almost equal learned values for upper and lower slopes are expected due to the structural symmetry of the saturation function (Figure 4) that appears in the closed-loop system.

We have proposed a two-player collaborative and iterative algorithm that iterates over two sub-phases that learn the Lyapunov function and use it to learn a controller to enlarge the RoA of the system.

The existing approaches that are close to the purpose of this paper are those that simultaneously synthesize a controller and maximize the stability region. Our work does not put any limit on the class of considered systems except the generic condition of hyperbolicity. For restricted classes such as polynomial systems, Sums-of-Squares (SOS) method leads to a bilinear optimization that is solved by some form of alternation (Jarvis-Wloszek et al., 2003; Majumdar et al., 2013). As a dual to Lyapunov-based methods, Majumdar et al. (2014) uses the notion of occupation measure to optimize a feedback controller for a polynomial system, but it has scalability issues due to its reliance on SDP optimization toolbox. Our work is different from this class of methods as our method is not limited to polynomial systems. Moreover, our method uses automatic differentiation that is naturally combined with neural networks to enjoy their superior scalability combined with SGD optimization.

Among data-driven approaches, Berkenkamp et al. (2017) use statistical models of the system to learn a controller with the assumption that the Lyapunov function is given. Our method is different as we learn the Lyapunov function and the controller together in an alternating fashion. Neural Lyapunov Control by Chang et al. (2019) is closer to our work while their approach is different from our method in multiple ways: They need to solve a costly global optimization problem for a component called falsifier to find the states on which the Lyapunov conditions are violated. We rather use the geometric properties of the level sets of the Lyapunov function to find the potentially unstable states that are chosen from a ring around the current stable sublevel set. In addition, Unlike Chang et al. (2019)’s method that needs to be done until the end to be usable on the system, the growing nature of our work allows the system be in action while the controller keeps improving.

In the context of reinforcement learning, our method can be seen as an actor-critic approach (Grondman et al., 2012; Bhasin et al., 2013; Lillicrap et al., 2015) where the actor tries to stabilize the system while the critic estimates the size of the RoA induced by the controller. The actor is the policy network and the critic is the network that implements the Lyapunov function.

In control theory, as the title of our work suggests, the proposed algorithm can be seen as an automatic version of the celebrated Lyapunov redesign method (Khalil and Grizzle, 2002; Hwang et al., 2013) where the Lyapunov function and the closed-loop system are re-designed together iteratively with the purpose of enlarging the stability region.

In this work, we assume the model of the system is given. This condition can be relaxed as suggested by Remark 3. Investigating this relaxation can be considered for the future. For example, learning a local model that adapts at each iteration could be one solution that adds another component to the algorithm and turns it into a three-player collaborative game.