Reinforcement learning is promising to highly nonlinear control robotic systems with large state and action space [1]. Until recently, significant progress has been made by combining advances in deep learning with reinforcement learning. Impressive results are obtained in a series of highdimensional robotic control tasks where sophisticated and hard-to-engineer behaviors can be achieved [2], [3], [4], [5]. However, the performance of an RL agent is by large evaluated through trial-and-error and RL could hardly provide any guarantee for the reliability of the learned control policy. Given a control system, regardless of which controller design method is used, the first and most important property of a system needs to be guaranteed is stability, because an unstable control system is typically useless and potentially dangerous [6]. A stable system is guaranteed to converge to the equilibrium or reference signal and it could recover to these targets even in the presence of parametric uncertainties and disturbances [7]. Thus stability is closely related to the robustness, safety, and reliability of the robotic systems. The most useful and general approach for studying the stability of robotic systems is Lyapunov’s method [8], which
trol Science and Engineering, Harbin Institute of Technology, China. mhhan@hit.edu.cn,lixianzhang@hit.edu.cn.
is dominant in control engineering [9], [10]. In Lyapunov’s method, a scalar “energy-like” function called Lyapunov function L is constructed to analyze the stability of the system. The controller is designed to ensure that the difference of Lyapunov function along the state trajectory is negative definite for all time instants so that the state goes in the direction of decreasing the value of Lyapunov function and eventually converges to the equilibrium [11], [12]. In learning methods, as the dynamic model is unknown, the “energy decreasing” condition has to be verified by trying out all possible consecutive data pairs in the state space, i.e., to verify infinite inequalities . Obviously, the “infinity” requirement is impossible thus making the direct exploitation of Lyapunov’s method impossible.
In this paper, we propose a data-based stability theorem and an actor-critic reinforcement learning algorithm to jointly learn the controller/policy and a Lyapunov critic function both of which are parameterized by deep neural networks, with a focus on stabilization and tracking tasks in robotic systems. The contribution of our paper can be summarized as follows: 1) a novel data-based stability theorem where only one inequality on the expected value over the state space needs to be evaluated; 2) a sample approximation of the stability condition proposed above is exploited to derive an actor-critic algorithm to search for a controller with asymptotic stability guarantee (in the number of data points); 3) we show through experiments that the learned controller could stabilize the systems when interfered by uncertainties such as unseen disturbances and system parametric variations of a certain extent. In our experiment, we show that the stability guaranteed controller is more capable of handling uncertainties compared to those without such guarantees in nonlinear control problems including classic CartPole stabilization tasks, control of 3D legged robots and manipulator, and reference tracking tasks for synthetic biology gene regulatory networks.
A. Related Works
In model-free reinforcement learning, stability is rarely addressed due to the formidable challenge of analyzing and designing the closed-loop system dynamics in a model-free manner [13], and the associated stability theory in model-free RL remains as an open problem [13], [14].
Recently, Lyapunov analysis is used in model-free RL to solve control problems with safety constraints [15], [16]. In [15], the Lyapunov-based approach for solving constrained Markov decision processes is proposed with a novel way of constructing the Lyapunov function through linear programming. In [16], the above results were further generalized to continuous control tasks. Even though Lyapunov-based methods were adopted in these results, neither of them addressed the stability of the system. In this paper, the sufficient conditions for a dynamic system being stable are derived. Furthermore, it is shown that these conditions can be verified through sampling and ensured through model-free learning.
Other interesting results on the stability of learning-based control systems are reported in recent years. In [17], an initial result is proposed for the stability analysis of deterministic nonlinear systems with optimal controller for infinite-horizon discounted cost, based on the assumption that discount is sufficiently close to 1. However, in practice, it is rather diffi-cult to guarantee the optimality of the learned policy unless certain assumptions on the system dynamics are made [18]. Furthermore, the exploitation of multi-layer neural networks as function approximators [19] only adds to the impracticality of this requirement. In this paper, it is shown that it is sufficient to ensure stability by satisfying the Lyapunov criterion that is evaluated on samples and thus one is exempt from finding the optimal/suboptimal solutions. In [20], local stability of Lipschitz continuous dynamic systems is analyzed by validating the “energy decreasing” condition on discretized points in the subset of state space with the help of a learned model (Gaussian process). Nevertheless, the discretization technique may become infeasible as the dimension and space of interest increases, limiting its application to rather simple and low-dimensional systems. In this paper, the proposed method is applicable to the general class of stochastic dynamic systems modeled by MDP and does not need to learn a model for stability analysis and controller design.
In this paper, we focus on the stabilization and tracking tasks for systems modeled by Markov decision process (MDP). The state of a robot and its environment at time t is given by the state , where S denotes the state space. The robot then takes an action
according to a stochastic policy
, resulting in the next state
. The transition of the state is modeled by the transition probability
. In both stabilization and tracking tasks, there always is a cost function
to measure how good or bad a state-action pair is.
In stabilization tasks, the goal is to find a policy such that the norm of state
goes to zero eventually, where
denotes the Euclidean norm. In this case, cost function
. In tracking tasks, we divide the state s into two vectors,
is composed of elements of s that are aimed at tracking the reference signal r, while
contains the rest. The reference signal could be the desired velocity, path and even the picture of grasping an object in a certain pose. For tracking tasks,
From a control perspective, both stabilization and tracking tasks are related to the asymptotic stability of the closed-loop system (or error system) under
, i.e., starting from an initial point, the trajectories of state always converge to the origin or reference trajectory. Let
denote the cost function under the policy
, the definition of stability studied in this paper is given as follows.
Definition 1: The stochastic system is said to be stable in mean cost if holds for any initial condition
is arbitrarily large then the stochastic system is globally stable in mean cost. The above definition is equivalent to the mean square stability [21], [22] when the cost c is chosen to be the norm of the state; it is also equivalent to the partial stability [23], [24] when
. Thus the stabilization and tracking tasks can be collectively summarized as finding a policy
such that the closed-loop system is stable in mean cost according to Definition 1.
Before proceeding, some notations are to be defined. denotes the distribution of starting states. The closed-loop transition probability is denoted as
. We also introduce the closed-loop state distribution at a certain instant t as
, which could be defined iteratively:
In this section, we propose the main assumptions and a new theorem for stability analysis of stochastic systems. We assume that the Markov chain induced by policy is ergodic with a unique stationary distribution
as commonly exploited by many RL literature [25], [26], [27], [28].
In Definition 1, stability is defined in relation to the set of starting states, which is also called the region of attraction (ROA). If the MSS system starts within the ROA, its trajectory will be surely attracted to the equilibrium. To build a data-based stability guarantee, we need to ensure that the states in ROA are accessible for the stability analysis. Thus the following assumption is made to ensure that every state in ROA has a chance to be sampled as the starting state.
Our approach is to construct/find a Lyapunov function of which the difference along the state trajectory is negative definite, so that the state goes in the direction of decreasing the value of Lyapunov function and eventually converges to the origin. The Lyapunov’s method has long been used for stability analysis and controller design in control theory [29], but mostly exploited along with a known model so that the energy decreasing condition on the entire state space could be transformed into one inequality regarding model parameters [6], [30]. In the following, we show that without a dynamic model, this “infinity” problem could be solved through sampling and sufficient conditions for a stochastic system to be stable in mean cost are given.
Theorem 1: The stochastic system is stable in mean cost if there exists a function and positive constants
, such that
where
First, on the left-hand-side, for all
S according to (1). Since the probability density function
is (assumed to be) a bounded function on S for all t, thus there exists a constant M such that
Second, the sequence converges point-wise to the function
. According to the Lebesgue’s Dominated convergence theorem [31], if a sequence
converges point-wise to a function f and is dominated by some integrable function g in the sense that,
Then we get
Thus taking the relations above into consideration, (3) infers
Since is a finite value and L is positive definite, it follows that
Suppose that there exists a state positive constant d such that
, or
. Since
for all starting states in
(Assumption 1), it follows that
, which is contradictory with (5). Thus
,
. Thus the system is stable in mean cost by Definition 1.
(1) directs the choice and construction of Lyapunov function, of which the details are deferred to Section IV. (2) is called the energy decreasing condition and is the major criterion for determining stability.
Remark 1: This remark is on the connection to previous results concerning the stability of stochastic systems. It should be noted that the stability conditions of Markov chains have been reported in [21], [32], however, of which the validation requires verifying infinite inequalities on the state space if S is continuous. On the contrary, our approach solely validates one inequality (2) related to the sampling distribution , which further enables data-based stability analysis and policy learning.
In this section, we propose an actor-critic RL algorithm to learn stability guaranteed policies for the stochastic system. First, we introduce the Lyapunov critic function how it is constructed. Then based on the maximum entropy actor-critic framework, we use the Lyapunov critic function in the policy gradient formulation.
A. Lyapunov Critic Function
In our framework, the Lyapunov critic plays a role in both stability analysis and the learning of the actor. To enable the actor-critic learning, the Lyapunov critic is designed to be dependent on s and a and satisfies
with the Lyapunov function L(s), such that it can be exploited in judging the value of (2). In view of the requirement above,
should be a non-negative function of the state and action,
. In this paper, we construct Lyapunov critic with the following parameterization technique,
where is the output vector of a fully connected neural network with parameter
. This parameterization ensures the positive definiteness of
, which is necessary since L(s) is positive definite according to (1) and L(s) is the expectation of
over the distribution of actions.
Theoretically, some functions, such as the norm of state and value function, naturally satisfy the basic requirement of being a Lyapunov function (1). These functions are referred to as Lyapunov candidates. However, Lyapunov candidates are conceptual functions without any parameterization, thus their gradient with respect to the controller is intractable and are not directly applicable in an actor-critic learning process. In the proposed framework, the Lyapunov candidate acts as a supervision signal during the training of . During training,
is updated to minimize the following objective function,
where is the approximation target related to the chosen Lyapunov candidate,
and D is the set of collected transition pairs. In [15] and [20], the value function has been proved to be a valid Lyapunov candidate where the approximation target is
where is the target network parameterized by
as typically used in the actor-critic methods [33], [19].
same structure as
, but the parameter
is updated through exponentially moving average of weights of
controlled by a hyperparameter
In addition to value function, the sum of cost over a finite time horizon could also be employed as Lyapunov candidate, which is exploited in model predictive control literature [34], [10] for stability analysis. In this case,
Here, the time horizon N is a hyperparameter to be tuned, of which the influence will be demonstrated in the experiment in Section V.
The choice of the Lyapunov candidate plays an important role in learning a policy. Value function evaluates the infinite time horizon and thus offers a better performance in general but is rather difficult to approximate because of significant variance and bias [35]. On the other hand, the finite horizon sum of cost provides an explicit target for learning a Lyapunov function, thus inherently reduces the bias and enhances the learning process. However, as the model is unknown, predicting the future costs based on the current state and action inevitably introduces variance, which grows as the prediction horizon extends. In principle, for tasks with simple dynamics, the sum-of-cost choice enhances the convergence of learning and robustness of the trained policies, while for complicated systems the choice of value function generally produces better performance. In this paper, we use both value function and sum-of-cost as Lyapunov candidates. Later in Section V, we will show the influence of these different choices upon the performance and robustness of trained policies.
B. Lyapunov Actor-Critic Algorithm
In this subsection, we will focus on how to learn the controller in a novel actor-critic framework called Lyapunov Actor-Critic (LAC), such that the inequality (2) is satisfied. The policy learning problem is summarized as the following constrained optimization problem,
where the stochastic policy is parameterized by a deep neural network
that is dependent on s and a Gaussian noise
. The constraint (11) is the parameterized inequality (2), which can be evaluated through sampling. One may be curious why in the second term of (13), only one Lyapunov critic is explicitly dependent on the stochastic policy, while the other dependent on the samples of the action. First, note that this estimator is also unbiased estimation of (2), although the variance may be increased compared to replacing a with
. From a more practical perspective, having the second Lyapunov critic explicitly dependent on
will introduce a term in the policy gradient that updates
to increase the value of L(s), which is contradictory to our goal of stabilization. (12) is the minimum entropy constraint borrowed from the maximum entropy RL framework to improve the exploration in the action space [33], and
is the desired bound. Solving the above constrained optimization problem is equivalent to minimizing the following objective function,
The value of Lagrange multipliers and
are adjusted by gradient ascent to maximize the following objectives respectively while being clipped to be positive,
During training, the Lagrange multipliers are updated by
where is the learning rate. The pseudo-code of the proposed algorithm is shown in Algorithm 1.
Algorithm 1 Lyapunov-based Actor-Critic (LAC)
Require: Maximum episode length N and maximum update steps M repeat
until (11) is satisfied
Fig. 1. Cumulative control performance comparison. The Y-axis indicates the total cost during one episode and the X-axis indicates the total time steps in thousand. The shadowed region shows the 1-SD confidence interval over 10 random seeds. Across all trials of training, LAC converges to stabilizing solution with comparable or superior performance compared with SAC and SPPO.
In this section, we illustrate five simulated robotic control problems to demonstrate the general applicability of the proposed method. First of all, the classic control problem of CartPole balancing from control and RL literature [36] is illustrated. Then, we consider more complicated highdimensional continuous control problems of 3D robots, e.g., HalfCheetah and FetchReach, using MuJoCo physics engine [37], a multi-joint Swimmer robot [38], and a full quadruped (Minitaur) simulated by PyBullet platform [39]. Last, we extend our approach to control autonomous systems in the cell, i.e., molecular synthetic biological gene regulatory networks (GRN). Specifically, we consider the problem of reference tracking for two GRNs [40].
The proposed method is evaluated for the following aspects:
• Convergence: does the proposed training algorithm converge with random parameter initialization and does the stability condition (2) hold for the learned policies;
• Performance: can the goal of the task be achieved or the cumulative cost be minimized;
• Stability: if (2) hold, are the closed-loop systems stable indeed and generating stable state trajectories;
• Robustness: how do the trained policies perform when faced with uncertainties unseen during training, such as parametric variation and external disturbances;
• Generalization: can the trained policies generalize to follow reference signals that are different from the one
seen during training. We compare our approach with soft actor-critic (SAC) [33], one of the state-of-the-art actor-critic algorithms that outperform a series of RL methods such as DDPG [41], PPO [42] on the continuous control benchmarks. The variant of safe proximal policy optimization (SPPO) [16], a Lyapunov-based method, is also included in the comparison. The original SPPO is developed to deal with constrained MDP, where safety constraints exist. In our experiments, we modify it to apply the Lyapunov constraints on the MDP tasks and see whether it can achieve the same stability guarantee as LAC. In CartPole, we also compare with the linear quadratic regulator (LQR), a classical model-based optimal control method for stabilization. For both algorithms, the hyperparameters are tuned to reach their best performance. The outline of this section is as follows. In Section V-A, the performance, and convergence of LAC are demonstrated and compared with SAC. In Section V-B, a straight forward demonstration of stability is made by comparing with the baseline method. In Section V-C, the ability of generalization and robustness of the trained policies are evaluated and analyzed. Finally, in Section V-D, we show the influence of choosing different Lyapunov candidates upon the performance and robustness of trained policies. The hyperparameters of LAC and the detailed experiment setup are deferred to Appendix [43]. The code for reproduction can be found in our GitHub repository 1.
A. Performance
In each task, both LAC, SAC, and SPPO are trained 10 times with random initialization, average total cost, and its variance during training are demonstrated in Figure 1. In the examples (a)-(c) and (e), SAC and LAC perform comparably in terms of the total cost at convergence and the speed of convergence, while SPPO could converge in Cartpole, FetcheReach, and Swimmer. In GRN and CompGRN (see Figure 1(d) and Fig. S9(b) in the supplementary material), SAC is not always able to find a policy that is capable of completing the control objective, resulting in the bad average performance. In the Minitaur example (see Figure 1(f)), SAC and SPPO can only converge to suboptimal solutions. On the contrary, LAC performs stably regardless of the random initialization. As shown in Figure 1, LAC converges stably in all experiments.
B. Evaluation of Stability
In this part, further comparison between the stabilityassured method (LAC) and that without such guarantee (SAC) is made, by demonstrating the closed-loop system dynamic with the trained policies. A distinguishing feature of the stability assured policies is that it can force and sustain the state or tracking error to zero. This could be intuitively demonstrated by the state trajectories of the closed-loop system.
We evaluated the trained policies in the GRN and CompGRN and the results are shown in Figure 2. In our experiments, we found that the LAC agents stabilize the systems well. All the state trajectories converge to the reference signal eventually (see Figure 2 a and c). On the contrary, without stability guarantee, the state trajectories either diverge (see Figure 2 b), or continuously oscillate around the reference trajectory (see Figure 2 d).
Fig. 2. State trajectories over time under policies trained by LAC and SAC in the GRN and CompGRN. In each experiment, the policies are tested over 20 random initial states and all the resulting trajectories are displayed above. The X-axis indicates the time and Y-axis shows the concentration of the target protein— Protein 1.
C. Empirical Evaluation on Robustness and Generalization
It is well-known that over-parameterized policies are prone to become overfitted to a specific training environment. The ability of generalization is the key to the successful implementation of RL agents in an uncertain real-world environment. In this part, we first evaluate the robustness of policies in the presence of parametric uncertainties and process noise. Then, we test the robustness of controllers against external disturbances. Finally, we evaluate whether the policy is generalizable by setting different reference signals. To make a fair comparison, we removed the policies that did not converge in SAC and only evaluate the ones that perform well during training. During testing, we found that SPPO appears to be prone to variations in the environment, thus the evaluation results are contained in Appendix [43].
1) Robustness to dynamic uncertainties: In this part, during the inference, we vary the system parameters and introduce process noises in the model/simulator to evaluate the algorithm’s robustness. In CartPole, we vary the length of pole l. In GRN, we vary the promoter strength and dissociation rate
. Due to stochastic nature in gene expression, we also introduce uniformly distributed noise ranging from
(we indicate the noise level by
) to the dynamic of GRN. The state trajectories of the closed-loop system under LAC and SAC agents in the varied environment are demonstrated in Figure 3.
As shown in Figure 3 (a, c), the policies trained by LAC are very robust to the dynamic uncertainties and achieve high
Fig. 3. LAC and SAC agents in the presence of dynamic uncertainties. The solid line indicates the average trajectory and shadowed region for the 1-SD confidence interval. In (a) and (b), the pole length is varied during the inference. In (c) and (d), three parameters are selected to reflect the uncertainties in gene expression. The X-axis indicates the time and Y-axis shows the angle of the pole in (a,b) and concentration of target protein in (c,d), respectively. The dashed line indicates the reference signal. The line in orange indicates the dynamic in the original environment. For each curve, only the noted parameter is different from the original setting.
tracking precision in each case. On the other hand, though SAC performs well in the original environment (Figure 3 (b, d)), it fails in all of the varied environments.
2) Robustness to disturbances: An inherent property of a stable system is to recover from perturbations such as external forces and wind. To show this, we introduce periodic external disturbances with different magnitudes in each environment and observe the performance difference between policies trained by LAC and SAC. We also include LQR as the model-based baseline. In CartPole, the agent may fall over when interfered by an external force, ending the episode in advance. Thus in this task, we measure the robustness of controllers through the death-rate, i.e., the probability of falling over after being disturbed. For other tasks where the episodes are always of the same length and we measure the robustness of controller by the variation in the cumulative cost. Under each disturbance magnitude, the policies are tested for 100 trials and the performance is shown in Figure 4.
As shown in Figure 4, the controllers trained by LAC outperform SAC and LQR to a great extent in CartPole and GRN (lower death rate and cumulative cost). In HalfCheetah, SAC and LAC are both robust to small external disturbances while LAC is more reliable to larger ones. In FetchReach, SAC and LAC perform reliably across all of the external disturbances. The difference between SAC and LAC becomes obvious in GRN, Swimmer, and Minitaur, where the dynamics are more vulnerable to the external disturbances. In all of the experiments, SPPO agents could hardly sustain any external disturbances.3) Generalization over different tracking references: In this part, we introduce four different reference signals that are unseen during training in the GRN: sinusoids with periods of 150 (brown) and 400 (blue), and the constant reference of
Fig. 4. Performance of LAC, SAC, SPPO and LQR in the presence of persistent disturbances with different magnitudes. The X-axis indicates the magnitude of the applied disturbance. The Y-axis indicates the death rate in CartPole (a) and the cumulative cost in other examples (b)-(d). All of the trained policies are evaluated for 100 trials in each setting.
Fig. 5. State trajectories under policies trained by LAC and SAC when tracking different reference signals. The solid line indicates the average trajectory and shadowed region for the 1-SD confidence interval. The Xaxis indicates the time and Y-axis shows the concentration of protein to be controlled. Dashed lines in different colors are the different reference signals: sinusoid with a period of 150 (brown); sinusoid with a period of 200 (sky-blue); sinusoid with a period of 400 (blue); constant reference of 8 (red); constant reference of 16 (green).
8 (red) and 16 (green). We also show the original reference signal used for training (skyblue) as a benchmark. Reference signals are indicated in Figure 5 by the dashed line in respective colors. All of the trained policies are tested for 10 times with each reference signal. The average dynamics of the target protein are shown in Figure 5 with the solid line, while the variance of dynamic is indicated by the shadowed area.
As shown in Figure 5, the policies trained by LAC could generalize well to follow previously unseen reference signals with low variance (dynamics are very close to the dashed lines), regardless of whether they share the same mathematical form with the one used for training. On the other hand, though SAC tracks the original reference signal well after the training trials without convergence being removed (see the sky-blue lines), it is still unable to follow some of the reference signals (see the brown line) and possesses larger variance than LAC.
D. Influence of Different Lyapunov candidates
As an independent interest, we evaluate the influence of choosing different Lyapunov candidates in this part. First, we adopt candidates of different time horizon to train policies in the CartPole example and compare their performance in terms of cumulative cost and robustness. Here,
implies using value function as the Lyapunov candidate. Both of the Lyapunov critics are parameterized as (6). For evaluation of robustness, we apply an impulsive force
instant and observe the death-rate of trained policies. The results are demonstrated in Figure 6.
Fig. 6. Influence of different Lyapunov candidates. In (a), the Y-axis indicates cumulative cost during training and the X-axis indicates the total time steps in thousand. (b) shows the death-rate of policies in the presence of instant impulsive force F ranging from 80 to 150 Newton.
As shown in Figure 6, both choices of Lyapunov candidates converge fast and achieve comparable cumulative cost at convergence. However, in terms of robustness, the choice of N plays an important role. As observed in Figure 6 (b), the robustness of the controller decreases as the time horizon N increases. Besides, it is interesting that LQR is more robust than SAC when faced with instant impulsive disturbance.
In this paper, we proposed a data-based approach for analyzing the stability of discrete-time nonlinear stochastic systems modeled by Markov decision process, by using the classic Lyapunov’s method in control theory. By employing the stability condition as a critic, an actor-critic algorithm is proposed to learn a controller/policy to ensure the closed-loop stability in stabilization and tracking tasks. We evaluated the proposed method in various examples and show that our method achieves not only comparable or superior performance compared with the state-of-the-art RL algorithm but also outperforms impressively in terms of robustness to uncertainties such as model parameter variations and external disturbances. For future work, it might be interesting to extend this method to constrained Markov decision process in which state and action constraints are considered. Also, to quantify the robustness induced by the stability will be investigated.
[1] J. Kober, J. A. Bagnell, and J. Peters, “Reinforcement learning in robotics: A survey,” The International Journal of Robotics Research, vol. 32, no. 11, pp. 1238–1274, 2013.
[2] J. Peters and S. Schaal, “Reinforcement learning of motor skills with policy gradients,” Neural networks, vol. 21, no. 4, pp. 682–697, 2008.
[3] S. L¨ockel, J. Peters, and P. Van Vliet, “A probabilistic framework for imitating human race driver behavior,” IEEE Robotics and Automation Letters, 2020.
[4] Y. Zhu, R. Mottaghi, E. Kolve, J. J. Lim, A. Gupta, L. Fei-Fei, and A. Farhadi, “Target-driven visual navigation in indoor scenes using deep reinforcement learning,” in 2017 IEEE international conference on robotics and automation (ICRA). IEEE, 2017, pp. 3357–3364.
[5] S. Gu, E. Holly, T. Lillicrap, and S. Levine, “Deep reinforcement learning for robotic manipulation with asynchronous off-policy updates,” in 2017 IEEE international conference on robotics and automation (ICRA). IEEE, 2017, pp. 3389–3396.
[6] J.-J. E. Slotine, W. Li, et al., Applied nonlinear control. Prentice hall Englewood Cliffs, NJ, 1991, vol. 199, no. 1.
[7] M. Vidyasagar, Nonlinear systems analysis. Siam, 2002, vol. 42.
[8] A. M. Lyapunov, The general problem of the stability of motion (in Russian). PhD Dissertation, Univ. Kharkov, 1892.
[9] K. J. ˚Astr¨om and B. Wittenmark, Adaptive control. Courier Corporation, 1989.
[10] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. Scokaert, “Constrained model predictive control: Stability and optimality,” Automatica, vol. 36, no. 6, pp. 789–814, 2000.
[11] M. Corless and G. Leitmann, “Continuous state feedback guaranteeing uniform ultimate boundedness for uncertain dynamic systems,” IEEE Transactions on Automatic Control, vol. 26, no. 5, pp. 1139–1144, 1981.
[12] A. Thowsen, “Uniform ultimate boundedness of the solutions of uncertain dynamic delay systems with state-dependent and memoryless feedback control,” International Journal of control, vol. 37, no. 5, pp. 1135–1143, 1983.
[13] L. Bus¸oniu, T. de Bruin, D. Toli´c, J. Kober, and I. Palunko, “Reinforcement learning for control: Performance, stability, and deep approximators,” Annual Reviews in Control, 2018.
[14] D. Gorges, “Relations between model predictive control and reinforcement learning,” IFAC-PapersOnLine, vol. 50, no. 1, pp. 4920–4928, 2017.
[15] Y. Chow, O. Nachum, E. Duenez-Guzman, and M. Ghavamzadeh, “A lyapunov-based approach to safe reinforcement learning,” arXiv preprint arXiv:1805.07708, 2018.
[16] Y. Chow, O. Nachum, A. Faust, M. Ghavamzadeh, and E. DuenezGuzman, “Lyapunov-based safe policy optimization for continuous control,” arXiv preprint arXiv:1901.10031, 2019.
[17] R. Postoyan, L. Bus¸oniu, D. Neˇsi´c, and J. Daafouz, “Stability analysis of discrete-time infinite-horizon optimal control with discounted cost,” IEEE Transactions on Automatic Control, vol. 62, no. 6, pp. 2736–2749, 2017.
[18] J. J. Murray, C. J. Cox, and R. E. Saeks, “The adaptive dynamic programming theorem,” in Stability and control of dynamical systems with applications. Springer, 2003, pp. 379–394.
[19] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” arXiv preprint arXiv:1509.02971, 2015.
[20] F. Berkenkamp, M. Turchetta, A. Schoellig, and A. Krause, “Safe model-based reinforcement learning with stability guarantees,” in Advances in neural information processing systems, 2017, pp. 908–918.
[21] L. Shaikhet, “Necessary and sufficient conditions of asymptotic mean square stability for stochastic linear difference equations,” Applied Mathematics Letters, vol. 10, no. 3, pp. 111–115, 1997.
[22] P. Bolzern, P. Colaneri, and G. De Nicolao, “Markov jump linear systems with switching transition rates: mean square stability with dwell-time,” Automatica, vol. 46, no. 6, pp. 1081–1088, 2010.
[23] V. I. Vorotnikov, “Partial stability and control: The state-of-the-art and development prospects,” Automation and Remote Control, vol. 66, no. 4, pp. 511–561, 2005.
[24] W. M. Haddad et al., “Finite-time partial stability and stabilization, and optimal feedback control,” Journal of the Franklin Institute, vol. 352, no. 6, pp. 2329–2357, 2015.
[25] R. S. Sutton, H. R. Maei, and C. Szepesv´ari, “A convergent o(n) temporal-difference algorithm for off-policy learning with linear function approximation,” in Advances in neural information processing systems, 2009, pp. 1609–1616.
[26] N. Korda and P. La, “On td (0) with function approximation: Concentration bounds and a centered variant with exponential convergence,” in International Conference on Machine Learning, 2015, pp. 626–634.
[27] J. Bhandari, D. Russo, and R. Singal, “A finite-time analysis of temporal difference learning with linear function approximation,” arXiv preprint arXiv:1806.02450, 2018.
[28] S. Zou, T. Xu, and Y. Liang, “Finite-sample analysis for sarsa with linear function approximation,” in Advances in Neural Information Processing Systems, 2019, pp. 8665–8675.
[29] E. Boukas and Z. Liu, “Robust stability and h/sub/spl infin//control of discrete-time jump linear systems with time-delay: an lmi approach,” in Decision and Control, 2000. Proceedings of the 39th IEEE Conference on, vol. 2. IEEE, 2000, pp. 1527–1532.
[30] S. Sastry, Nonlinear systems: analysis, stability, and control. Springer Science & Business Media, 2013, vol. 10.
[31] H. L. Royden, Real analysis. Krishna Prakashan Media, 1968.
[32] S. P. Meyn and R. L. Tweedie, Markov chains and stochastic stability. Springer Science & Business Media, 2012.
[33] T. Haarnoja, A. Zhou, K. Hartikainen, G. Tucker, S. Ha, J. Tan, V. Kumar, H. Zhu, A. Gupta, P. Abbeel, et al., “Soft actor-critic algorithms and applications,” arXiv preprint arXiv:1812.05905, 2018.
[34] D. Q. Mayne and H. Michalska, “Receding horizon control of nonlinear systems,” IEEE Transactions on automatic control, vol. 35, no. 7, pp. 814–824, 1990.
[35] J. Schulman, P. Moritz, S. Levine, M. Jordan, and P. Abbeel, “Highdimensional continuous control using generalized advantage estimation,” arXiv preprint arXiv:1506.02438, 2015.
[36] A. G. Barto, R. S. Sutton, and C. W. Anderson, “Neuronlike adaptive elements that can solve difficult learning control problems,” IEEE transactions on systems, man, and cybernetics, no. 5, pp. 834–846, 1983.
[37] E. Todorov, T. Erez, and Y. Tassa, “Mujoco: A physics engine for model-based control,” in 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems. IEEE, 2012, pp. 5026–5033.
[38] Y. Tassa, Y. Doron, A. Muldal, T. Erez, Y. Li, D. d. L. Casas, D. Budden, A. Abdolmaleki, J. Merel, A. Lefrancq, et al., “Deepmind control suite,” arXiv preprint arXiv:1801.00690, 2018.
[39] E. Coumans and Y. Bai, “Pybullet, a python module for physics simulation for games, robotics and machine learning,” GitHub repository, 2016.
[40] M. B. Elowitz and S. Leibler, “A synthetic oscillatory network of transcriptional regulators,” Nature, vol. 403, no. 6767, p. 335, 2000.
[41] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” arXiv preprint arXiv:1509.02971, 2015.
[42] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” arXiv preprint arXiv:1707.06347, 2017.
[43] M. Han, L. Zhang, J. Wang, and W. Pan, “Actor-critic reinforcement learning for control with stability guarantee,” arXiv preprint arXiv:2004.14288, 2020.
[44] G. Brockman, V. Cheung, L. Pettersson, J. Schneider, J. Schulman, J. Tang, and W. Zaremba, “Openai gym,” arXiv preprint arXiv:1606.01540, 2016.
[45] N. Strelkowa and M. Barahona, “Switchable genetic oscillator operating in quasi-stable mode,” Journal of The Royal Society Interface, vol. 7, no. 48, pp. 1071–1082, 2010.
[46] A. Sootla, N. Strelkowa, D. Ernst, M. Barahona, and G.-B. Stan, “On periodic reference tracking using batch-mode reinforcement learning with application to gene regulatory network control,” in 52nd IEEE conference on decision and control. IEEE, 2013, pp. 4086–4091.
[47] Y. Bengio, N. Boulanger-Lewandowski, and R. Pascanu, “Advances in optimizing recurrent networks,” in 2013 IEEE International Conference on Acoustics, Speech and Signal Processing. IEEE, 2013, pp. 8624– 8628.
[48] I. Bello, B. Zoph, V. Vasudevan, and Q. V. Le, “Neural optimizer search with reinforcement learning,” in Proceedings of the 34th International Conference on Machine Learning-Volume 70. JMLR. org, 2017, pp. 459–468.
[49] Z. Wang, T. Schaul, M. Hessel, H. Van Hasselt, M. Lanctot, and N. De Freitas, “Dueling network architectures for deep reinforcement learning,” arXiv preprint arXiv:1511.06581, 2015.
We set up the experiment using OpenAi Gym [44], DeepMind Control Suite [38] and PyBullet physics simulation platforms [39]. A snapshot of the adopted environments in this paper can be found in Figure S7.
Fig. S7. Snapshot of environments using OpenAI Gym.
A. CartPole
In this experiment, the controller is expected to maintain the pole vertically at a target position x = 0. This is a modified version of CartPole in [44] with continuous action space. The action is the horizontal force applied upon the cart (and
represents the maximum of position and angle, respectively,
and
. The episode ends if
or
and the episodes end in advance. Cost function
. The episodes are of length 250. For robustness evaluation in Section V-C, we apply an impulsive disturbance force F on the cart every 20 seconds, of which the magnitude ranges from 80 to 150 and the direction is opposite to the direction of control input. In Section V-D, the impulsive disturbance has the same magnitude range and direction with that in Section V-C, but only applied once at instant t = 100.
B. HalfCheetah
HalfCheetah is a modified version of that in Gym’s robotics environment [44]. The task is to control a HalfCheetah (a 2-legged simulated robot) to run at the speed of 1 m/s. The reward is where v is the forward speed of the HalfCheetah. The control input is the torque applied on each joint, ranging from -1 to 1. The episodes are of length 200.
For robustness evaluation in Section V-C, we apply an impulsive disturbance torque on each joint every 20 seconds, of which the magnitude ranges from 0.2 to 2.0 and the direction is opposite to the direction of control input.
C. FetchReach
We modify the FetchReach in Gym’s robotics environment [44] to a cost version, where the controller is expected to control the manipulator’s end effector to reach a random goal position. The cost is designed as c = d, where d is the distance between goal and end-effector. The control input is the torque applied on each joint, ranging from -1 to 1. The episodes are of length 200.
For robustness evaluation in Section V-C, we apply an impulsive disturbance torque on each joint every 20 seconds, of which the magnitude ranges from 0.2 to 2.0 and the direction is opposite to the direction of control input.
D. Swimmer
Swimmer is a modified version based on the environment in the DeepMind control suite [38]. The task is to control a multi-joint snake robot to run at a speed of 1 m/s. The reward is is the forward speed of the Swimmer. The control input is the torque applied on each joint, ranging from -1 to 1. The episodes are of length 250.
For robustness evaluation in Section V-C, we apply an impulsive disturbance torque on each joint every 20 seconds, of which the magnitude ranges from 0.2 to 1.0 and the direction is opposite to the direction of control input.
E. Minitaur
This example is borrowed from the PyBullet environment [39]. In Minitaur, the agent controls the Ghost Robotics Minitaur quadruped to run forward at the speed of 1 m/s. The reward is is the forward speed of the robot. The control input is the desired pose of each actuator. The episodes are of length 500.
For robustness evaluation in Section V-C, we apply an impulsive disturbance on the input channel every 20 seconds, of which the magnitude ranges from 0.2 to 1.0 and the direction is opposite to the direction of control input.
Since the gene regulatory networks (GRN) considered here is in nano-scale whose physical property is different from the ones considered in Section S1. Also, GRN can exhibit interesting oscillatory behavior. We illustrate this example separately in this section.
A. Mathematical model of GRN
In this example, we consider a classical dynamical system in systems/synthetic biology, the repressilator, which we use to illustrate the reference tracking task at hand. The repressilator is a synthetic three-gene regulatory network where the dynamics of mRNAs and proteins follow an oscillatory behavior [40]. A discrete-time mathematical description of the repressilator, which includes both transcription and translation dynamics, is given by the following set of discrete-time equations:
Here, (resp.
) denote the concentrations of the mRNA transcripts (resp. proteins) of genes 1, 2, and 3, respectively.
are i.i.d. uniform noise ranging from
. During training,
and for evaluation
respectively in Section V-C.
denote the maximum promoter strength for their corresponding gene,
denote the mRNA degradation rates,
denote the protein degradation rates,
denote the protein production rates, and
are the dissociation constants. The set of equations in Eq.(15) corresponds to a topology where gene 1 is repressed by gene 2, gene 2 is repressed by gene 3, and gene 3 is repressed by gene 1. dt is the discretization time step.
In practice, only the protein concentrations are observed and given as readouts, for instance via fluorescent markers (e.g., green fluorescent protein, GFP or red fluorescent protein, mCherry). The control scheme will be implemented by light control signals which can induce the expression of genes through the activation of their photo-sensitive promoters. To simplify the system dynamics and as it is usually done for the repressilator model [40], we consider the corresponding parameters of the mRNA and protein dynamics for different genes to be equal. More background on mathematical modeling and control of synthetic biology gene regulatory networks can be referred to [45], [46]. In this example, the parameters are as follows:
In FigS8, a single snapshot of the state temporal evolution without is depicted. We uniformly initialized between 0 to 5, i.e.,
, which is the range we train the policy in Section V, persistent oscillatory behavior is also exhibiting similar to the snapshot in Fig S8.
Fig. S8. A snapshot of the natural oscillatory behavior of a repressilator system consisting of 3 genes. The oscillations have a period of approximately 150 arbitrary time units. The X-axis denotes time and Y-axis denotes the value/concentration of each state.
Fig. S9. Cumulative control performance comparison for synthetic GRN. The Y-axis indicates the total cost during one episode and the X-axis indicates the total time steps in thousand. The shadowed region shows the 1-SD confidence interval over 10 random seeds. Across all trials of training, LAC converges to stabilizing solution with comparable or superior performance compared with SAC and SPPO. Figure.S9 (a) is identical to Figure.1 (d).
TABLE S1 HYPERPARAMETERS OF LAC
For LAC, there are two networks: the policy network and the Lyapunov critic network. For the policy network, we use a fully-connected MLP with two hidden layers of 256 units, outputting the mean and standard deviations of a Gaussian distribution. As mentioned in section IV, it should be noted that the output of the Lyapunov critic network is a square term, which is always non-negative. More specifically, we use a fully-connected MLP with two hidden layers and one output layer with different units as in Table S1, outputting the feature vector . The Lyapunov critic function is obtained by
. All the hidden layers use Relu activation function and we adopt the same invertible squashing function technique as [33] to the output layer of the policy network.
Fig. S10. Value of Lagrange multiplier during the training of LAC policies. The Y-axis indicates the value of
and the X-axis indicates the total time steps in thousand. The shadowed region shows the 1-SD confidence interval over 10 random seeds. The value of
gradually drops and becomes zero at convergence, which implies the satisfaction of stability condition.
The convergence of LAC and validation of the stability guarantee can also be checked by observing the value of Lagrange multipliers. When (11) is satisfied, will continuously decrease until it becomes zero. Thus by checking the value of
, the satisfaction of stability condition during training and at convergence could be validated. In practice, the value of
is clipped under the maximum value 1 in case that
grows too large due to the violation of stability condition during the early training stage, resulting in the inappropriate step length for the policy update. Clipping is a useful technique to prevent instability of optimization, especially in gradient-based methods, see [42], [47], [48], [49]. In Figure S10, the value of
during training is demonstrated. Across all training trials in the experiments,
converges to zero eventually, which implies that the stability guarantee is valid. It is also shown that the clipping technique only makes effect in the early training stage in the FetchReach (see (c)) while in other experiments it is not triggered.
In this part, we evaluate the robustness and generalization ability of policies trained by SPPO in the same. First, the robustness of the policies is tested by perturbing the parameters and adding noise in the Cartpole and Repressilator environment, as described in Section V-C.1. Generalization of the policies is evaluated by setting reference signals that are unseen during training. State trajectories of the above experiments are demonstrated in Figure S11 and Figure S12, respectively. As demonstrated in the figures, the SPPO policies could hardly deal with previously unseen uncertainty or reference signals and failed in all of the Repressilator experiments.
The SPPO algorithm is originally developed for the control tasks with safety constraints, i.e. keeping the expectation of discounted cumulative safety cost below a certain threshold. Though Lyapunov’s method is exploited, the approach is not aimed at providing stability guarantee.
Fig. S11. State trajectories over time under policies trained by SPPO and tested in the presence of parametric uncertainties and process noise, for CartPole and Repressilator. The setting of the uncertainties is the same as in Section V-C.1.
Fig. S12. State trajectories under policies trained by SPPO when tracking different reference signals. The setting of the uncertainty is the same as in Section V-C.3.