Safe Reinforcement Learning Using Robust MPC

TL;DR

This paper proposes a method combining Reinforcement Learning and Model Predictive Control to achieve safe, stable, and near-optimal control in uncertain systems, addressing safety concerns often overlooked in RL.

Contribution

It introduces a novel integration of RL with MPC to enhance safety and stability while maintaining optimality in control tasks.

Findings

The combined RL-MPC controller ensures safety and stability in uncertain environments.

Simulation results demonstrate improved safety guarantees over standard RL methods.

The approach balances optimality and robustness in control performance.

Abstract

Reinforcement Learning (RL) has recently impressed the world with stunning results in various applications. While the potential of RL is now well-established, many critical aspects still need to be tackled, including safety and stability issues. These issues, while partially neglected by the RL community, are central to the control community which has been widely investigating them. Model Predictive Control (MPC) is one of the most successful control techniques because, among others, of its ability to provide such guarantees even for uncertain constrained systems. Since MPC is an optimization-based technique, optimality has also often been claimed. Unfortunately, the performance of MPC is highly dependent on the accuracy of the model used for predictions. In this paper, we propose to combine RL and MPC in order to exploit the advantages of both and, therefore, obtain a controller which is optimal and safe. We illustrate the results with a numerical example in simulations.

Equations

$$\mathbb{P}\left[{\mathrm{s}}_{+}|\mathrm{s},\mathrm{a}\right].$$

$${\pi }_{\star }:=\arg \underset{\pi }{\min }J(\pi )$$

$$V_{\star }({\mathrm{s}}_0)$$

$$Q_{\star }(\mathrm{s},\mathrm{a}):=\ell (\mathrm{s},\mathrm{a})+\gamma \mathbb{E}\left[V_{\star }({\mathrm{s}}_{+})|\mathrm{s},\mathrm{a}\right],$$

$${\theta }^{\star }:=\underset{\theta }{\min }$$

$$\psi ({\mathrm{s}}_{k+1},{\mathrm{s}}_k,{\mathrm{a}}_k,\theta )={\left(\ell ({\mathrm{s}}_k,{\mathrm{a}}_k)+\gamma V_{{\theta }_k}({\mathrm{s}}_{k+1})-Q_{\theta }({\mathrm{s}}_k,{\mathrm{a}}_k)\right)}^2$$

$${\mathbb{E}}_{{\tau }^{{\pi }_{\theta }}}\left[\psi ({\mathrm{s}}_{k+1},{\mathrm{s}}_k,{\mathrm{a}}_k,\theta )\right]=J({\pi }_{\theta })$$

$$\xi (\mathrm{s},\hat{\pi }(\mathrm{s}))\le 0,$$

$${\mathrm{S}}_{+}\left(\mathrm{s},\mathrm{a}\right)=\{{\mathrm{s}}_{+}|\mathbb{P}\left[{\mathrm{s}}_{+}|\mathrm{s},\mathrm{a}\right]>0\}.$$

$${\hat{\mathrm{S}}}_{+}(\mathrm{s},\mathrm{a},\theta ):=\{{\mathrm{s}}_{+}|{\mathrm{g}}_{\theta }({\mathrm{s}}_{+},\mathrm{s},\mathrm{a})\le 0\}.$$

$${\hat{\mathrm{S}}}_{+}(\mathrm{s},\mathrm{a},\theta )\supseteq {\mathrm{S}}_{+}(\mathrm{s},\mathrm{a}),$$

$$\chi \left(\theta \right)=\underset{\mathrm{s},\mathrm{a}}{\sup }\mathbb{E}\left[{\mathrm{s}}_{+}\notin {\hat{\mathrm{S}}}_{+}(\mathrm{s},\mathrm{a},\theta )|\mathrm{s},\mathrm{a}\right]\in [0,1],$$

$$\eta \left(\mathcal{D}\right)=\mathbb{P}\left[\chi \left(\theta \right)=0|\mathcal{D}\right],$$

$$\eta \left(\mathcal{D}\right)=1$$

$${\mathrm{s}}_{k+1}\in {\hat{\mathrm{S}}}_{+}({\mathrm{s}}_k,{\mathrm{a}}_k,\theta ),$$

$${\mathrm{S}}_{k+1}^{{\pi }_{\theta }}$$

$$\xi ({\mathrm{S}}_k^{{\pi }_{\theta }},{\pi }_{\theta }({\mathrm{S}}_k^{{\pi }_{\theta }}))\le 0,$$

$${\theta }^{\star }:=\underset{\theta }{\min }$$

$$s.t.$$

$$C\mathrm{s}+D\mathrm{a}+\bar{\mathrm{c}}\le 0.$$

$$Q_{\theta }(\mathrm{s},\mathrm{a}):=$$

$$\displaystyle\min_{\boldsymbol{\mathrm{z}}}\ \$$

$$+{\gamma }^N\left({\mathrm{x}}_N^{\top }P{\mathrm{x}}_N+{\mathrm{p}}^{\top }{\mathrm{x}}_N\right)$$

$$\displaystyle\mathrm{s.t.}\ \$$

$${\mathrm{x}}_{k+1}=A{\mathrm{x}}_k+B{\mathrm{u}}_k+\mathrm{b},$$

$$C{\mathrm{x}}_k+D{\mathrm{u}}_k+{\mathrm{c}}_k\le 0,$$

$$G{\mathrm{x}}_N+\mathrm{g}\le 0,$$

$$V_{\theta }(\mathrm{s}):=\underset{\mathrm{a}}{\min }{Q}_{\theta }(\mathrm{s},\mathrm{a}),$$

$$\theta =\left\{H,\mathrm{h},P,\mathrm{p},A,B,\mathrm{b},\bar{\mathrm{c}},K,{\theta }_{\mathrm{W}}\right\}.$$

$${\hat{\mathrm{S}}}_{+}(\mathrm{s},\mathrm{a},\theta )=A\mathrm{s}+B\mathrm{a}+\mathrm{b}+{\mathrm{W}}_{\theta }\supseteq {\mathrm{S}}_{+}(\mathrm{s},\mathrm{a}).$$

$${\mathrm{E}}_{k+1}$$

$${\mathrm{S}}_{k+1}^{{\pi }_{\theta }^{\mathrm{MPC}}}={\hat{\mathrm{S}}}_{+}({\mathrm{S}}_k^{{\pi }_{\theta ,k}^{\mathrm{MPC}}},{\mathrm{a}}_k-K{\mathrm{E}}_k,\theta ),$$

$$C_K:=C-DK,A_K:=A-BK.$$

$$C{\mathrm{x}}_k+D{\mathrm{u}}_k+C_K{\mathrm{e}}_k+\bar{\mathrm{c}}\le 0,$$

Full text

I Introduction

Reinforcement Learning (RL) is a technique for solving problems involving Markov Decision Processes (MDP) [1]. In RL, rather than modeled state transition probabilities, samples and observed costs (or rewards) are used. RL algorithms enabled computers beating Chess and Go masters [2], and robots learning to walk or fly without supervision [3, 4].

For each state $\boldsymbol{\mathrm{s}}$ the optimal action $\boldsymbol{\mathrm{a}}$ is computed as the the optimal feedback policy $\boldsymbol{\mathrm{\pi}}(\boldsymbol{\mathrm{s}})$ for the real system either directly (policy search methods) [5, 6] or indirectly (SARSA, $Q$-learning) [7]. In the latter, the optimal policy $\boldsymbol{\mathrm{\pi}}(\boldsymbol{\mathrm{s}})$ is indirectly obtained as the minimizer of the so-called action-value function $Q(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{a}})$ over the action or input $\boldsymbol{\mathrm{a}}$. In both cases, either $\boldsymbol{\mathrm{\pi}}$ or $Q$ are typically approximated by a function approximator: Deep Neural Network (DNN) are very commonly used for that purpose in recent applications.

While RL has demonstrated in practice a huge potential, properties that are typically expected from a controller, such as, e.g., some form of stability and safety, are hard to guarantee, especially when relying on a DNN as a function approximator. In general, any safety-enforcing approach requires either (a) to collect data from the real system thus incurring catastrophic events, (b) to use a model of the system to generate a sufficient amount of stochastic simulations, or (c) to model the uncertainty underlying the system, as we will detail later. Some approaches have been developed in order to guarantee some form of safety: see, e.g., the excellent survey in [8] and references therein. Most approaches, however, do not strictly guarantee that a given set of constraints is never violated, but rather that violations are rare events. Some approaches propose to project the action resulting from a DNN onto a safe set: in [9] each constraint is approximated as a ReLu or an additional cost and in [10] a QP is used. In both approaches the nominal prediction is used and uncertainty is neglected in the predictions, similarly to the approach of [11]. Projection approaches have been analyzed in the general case in [12].

The combination of learning and control techniques has been proposed in, e.g., [13, 14, 15, 16, 17]. The combination of RL and the linear quadratic regulator has been presented in [18, 19]. To the best of our knowledge, [11, 20, 21] are the first works proposing to use NMPC as a function approximator in RL. While strategies for providing some form of safety have been developed [8], to the best of the authors’ knowledge, none of these approaches is able to strictly satisfy some set of constraints at all time. Rather, constraint violation is strongly penalized in [11] and some of the approaches in [8]. The only contribution providing robust constraint satisfaction guarantees is [22], where a linear feedback policy is learned.

In this paper we propose an RL formulation based on MPC which addresses the issue of safety, which we did not rigorously enforce in [11, 20]. We summarize next two existing approaches and the scheme we propose, which combines them:

Robust MPC, Business-as-usual

a low-dimensional computationally tractable uncertainty set is first identified and then used to formulate a robust MPC problem, see (20).

RL, Business-as-usual

penalizing violations with a suitably high cost lets the optimization procedure yield a policy which tends to not violate the constraints. Safety is typically not strictly guaranteed and only few results provide weak guarantees.

Safe RL MPC

In the approach we propose, a robust MPC problem is formulated similarly to (a). Similarly to (b), RL updates the parametrization of the robust MPC scheme, and of the safety constraint to reduce conservatism while preserving safety.

Safe RL-MPC is based on the approach first advocated in [11, 20]. The scheme can be seen from two alternative points of view: (a) MPC is used as a function approximator within RL in order to provide safety and stability guarantees; and (b) RL is used in order to tune the MPC parameters, thus improving closed-loop performance in a data-driven fashion. Since safety is fundamental not only during exploitation, but also during exploration, we also address the issue of guaranteeing constraint satisfaction during this phase.

Another important contribution of this paper is the development of an efficient way to deal with the enormous amount of data typically collected by autonomous systems. In particular, the introduction of a nominal (potentially inaccurate) linear model allows one to significantly reduce the amount of stored data. Further efficiency is obtained by exploiting convexity and using a low-dimensional approximation of the uncertainty set.

We first present RL at a conceptual level and propose adaptations in order to make RL applicable to the safety-enforcing setup. The main issue to be tackled is related to the safety constraints, which need to be enforced when updating the function approximator parameter. We formulate the update by resorting to a constrained optimization problem, similarly to what has been proposed in [20]. The proposed safe RL can be directly applied to $Q$-learning, but actor-critic techniques require some adaptation when the input space is continuous and restricted by safety requirements as discussed in [23].

Contributions

This paper proposes an approach to combine MPC and RL so as to guarantee safety. In order to limit the complexity of the algorithm, we (a) rely on a linear system with bounded disturbances for predictions, and (b) propose an ad-hoc formulation of the robust constraint satisfaction which relies on the linear model to greatly reduce the amount of data to be stored. Finally, adaptations to the standard RL algorithms to guarantee that the parameter update does not jeopardize safety are introduced.

The paper is structured as follows. We introduce the problem of safe RL in general terms in Section II, while the rest of the paper specializes on the case of linear systems. We propose a tailored function approximator based on robust MPC in Section III and discuss the efficient use of data in Section IV. The necessary modifications to the standard RL algorithms are proposed in Section V and the whole framework is tested in simulations in Section VI. Conclusions and an outline for future research are given in Section VII.

Notation

$a$ is scalar, $\boldsymbol{\mathrm{a}}\in\mathbb{R}^{n_{\boldsymbol{\mathrm{a}}}}$ is a vector with components $\boldsymbol{\mathrm{a}}_{i}$, $A$ is a matrix with rows $A_{i}$ and $\boldsymbol{\mathrm{A}}$, $\mathcal{A}$ are sets. For any set, $|\cdot|$ defines its cardinality. For any function $\boldsymbol{\mathrm{f}}(\boldsymbol{\mathrm{x}})$, we define $\boldsymbol{\mathrm{f}}(\boldsymbol{\mathrm{X}}):=\{\,\boldsymbol{\mathrm{f}}(\boldsymbol{\mathrm{x}})\,|\,\boldsymbol{\mathrm{x}}\in\boldsymbol{\mathrm{X}}\,\}$ and denote $\boldsymbol{\mathrm{f}}(\boldsymbol{\mathrm{x}})\leq 0,\ \forall\,\boldsymbol{\mathrm{x}}\in\boldsymbol{\mathrm{X}}$ as $\boldsymbol{\mathrm{f}}(\boldsymbol{\mathrm{X}})\leq 0$. The only exceptions are $J$, $V$, $Q$ which are scalar functions, but denoted by capital letters in the literature.

II Background and Safe RL Formulation

In this paper, we consider real system dynamics described as a Markov Process (MP) with continuous state $\boldsymbol{\mathrm{s}}$ and action $\boldsymbol{\mathrm{a}}$, with state transitions $\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{a}}\to\boldsymbol{\mathrm{s}}_{+}$ having the underlying conditional probability density

[TABLE]

We furthermore consider a deterministic policy delivering the control input as $\boldsymbol{\mathrm{a}}=\boldsymbol{\mathrm{\pi}}\left(\boldsymbol{\mathrm{s}}\right),$ resulting in state distribution $\tau^{\boldsymbol{\mathrm{\pi}}}$. The RL problem then reads as

[TABLE]

where $\ell$ is called stage cost in optimal control and $-\ell$ instantaneous reward in RL. The scalar $\gamma$ is a discount factor, typically smaller than $1$ in RL, and $1$ in MPC. Note that in (2) we provide the definition for non-episodic settings, but the developments of the paper readily apply to episodic settings. The value function $V_{\star}(\boldsymbol{\mathrm{s}})$ is the optimal cost, obtained by applying the optimal policy $\boldsymbol{\mathrm{\pi}}_{\star}$, i.e.,

[TABLE]

and the Bellman Equation defines the action-value function

[TABLE]

where the expectation is taken over the state transition (1). The various forms of RL use parametric approximations $V_{\theta}$, $Q_{\theta}$, $\boldsymbol{\mathrm{\pi}}_{\theta}$ of $V_{\star}$, $Q_{\star}$, $\boldsymbol{\mathrm{\pi}}_{\star}$ in order to find the parameter $\boldsymbol{\mathrm{\theta}}^{\star}$ which (approximately) solves (2) either directly or by approximately solving (4) in a sampled-based fashion. For both approaches we summarize this as

[TABLE]

where function $\psi$ depends on the specific algorithm, e.g.,

[TABLE]

and $n=1$ in recursive $Q$ learning formulations; and

[TABLE]

in policy gradient approaches [1]. To optimize performance, the system ought to be controlled by using the best available policy $\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$, this is referred to as exploitation. However, in order for Problem (5) to be well-posed in general, it is necessary to also collect data by deviating from $\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$ and implementing a different policy $\boldsymbol{\mathrm{\pi}}^{\mathrm{e}}$, this is referred to as exploration.

Among others, one of the main difficulties related with RL is safety enforcement, e.g., collision avoidance. In this paper we define safety through a set of constraints

[TABLE]

where we note $\mathrm{supp}\left(\tau^{\boldsymbol{\mathrm{\hat{\pi}}}}\right)$ the support of the distribution of the MP (1) subject to policy $\hat{\boldsymbol{\mathrm{\pi}}}$. Ideally, condition (8) should be satisfied at all times with unitary probability, both during exploitation, i.e., $\boldsymbol{\mathrm{\hat{\pi}}}=\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$, and exploration, i.e., $\boldsymbol{\mathrm{\hat{\pi}}}=\boldsymbol{\mathrm{\pi}}^{\mathrm{e}}\neq\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$. Note that (8) can only hold if process (1) has bounded support.

Enforcing (8) poses two major challenges: (i) either the support of (1) or the support of $\tau^{\boldsymbol{\mathrm{\hat{\pi}}}}$ must be known or estimated; (ii) given knowledge on either supports, a policy satisfying (8) must be designed.

Problem (i) is fundamental, since one can never have the guarantee of being able to observe the full support of (1). Arguably, a reasonable approach can be to approximate the support based on the information extracted from the available samples, or on a prior, or on both. We assume that collected data is informative such that, in the limit for an infinite amount of data, the support is reconstructed exactly for policy $\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$. A theoretical justification of this is beyond the scope of this paper. In order to construct an approximation of the MP support, we first define the dispersion set confining the state transitions as:

[TABLE]

Since $\boldsymbol{\mathrm{S}}_{+}$ is not known a priori, we introduce a parametrized approximation $\boldsymbol{\mathrm{\hat{S}}}_{+}$ based on function $\boldsymbol{\mathrm{g}}_{\boldsymbol{\mathrm{\theta}}}$ given by:

[TABLE]

In order to enforce safety, $\boldsymbol{\mathrm{\hat{S}}}_{+}$ must be an outer approximation of set $\boldsymbol{\mathrm{S}}_{+}$, i.e., $\boldsymbol{\mathrm{\theta}}$ must be chosen such that:

[TABLE]

We label this condition Safe-Design Constraint (SDC), since it restricts the values that $\boldsymbol{\mathrm{\theta}}$ can take based on safety concerns. In order to discuss safety in mathematically simple terms, let us introduce the worst-case mass of (1) outside of $\boldsymbol{\mathrm{\hat{S}}}_{+}$, defined as:

[TABLE]

where the expected value is taken over (1). For a given $\boldsymbol{\mathrm{\theta}}$, (11) is ensured if $\chi\left(\boldsymbol{\mathrm{\theta}}\right)=0$. However, $\chi\left(\boldsymbol{\mathrm{\theta}}\right)$ is known only if we assume the real system dynamics (1) are known. Otherwise, in the Bayesian context, $\chi\left(\boldsymbol{\mathrm{\theta}}\right)$ ought to be treated as a random variable, reflecting our imperfect knowledge of it, conditioned on the current data $\mathcal{D}$ (i.e., knowledge) we have of the system. We therefore consider

[TABLE]

which provides a formal definition of the probability that $\boldsymbol{\mathrm{\hat{S}}}_{+}$ is capturing the real system dispersion provided the data $\mathcal{D}$. While computing $\eta$ is difficult in the general case, its definition allows us to discuss in rigorous terms the fundamental limitation provided next.

Fundamental Limitation 1

For any $\mathcal{D}$, if $\boldsymbol{\mathrm{\hat{S}}}_{+}\subset\mathbb{R}^{n_{\boldsymbol{\mathrm{s}}}}$, it is impossible to guarantee that

[TABLE]

without introducing any additional assumption on (1).

In other words, given a finite set of samples and, possibly, prior but not absolute knowledge of the system, one does not have enough information to construct a set containing all future samples, unless this set is $\mathbb{R}^{n_{\boldsymbol{\mathrm{s}}}}$. This is a fundamental limitation of robust constraint satisfaction and, therefore, it is independent of the proposed approach. In the following, we will therefore discuss $\eta$-safety, underlining that (14) cannot be achieved in practice.

Arguably, by introducing additional assumptions restricting the function space to which (1) can belong, one can envision overcoming this limitation. The investigation of this topic is beyond the scope of this paper and will be the subject of future research.

In order to formulate some form of SDC with limited information, we introduce the sample-based form of (11) as

[TABLE]

For problem (ii), the main challenge is to find values of $\boldsymbol{\mathrm{\theta}}$ such that the policy $\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$ strictly satisfies the safety constraints.

We define safety based on the set dispersion propagation under policy $\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$, which reads as

[TABLE]

Definition 1** ($\eta$-safe Policy)**

For a given data set $\mathcal{D}$ and a given set of initial conditions $\mathcal{S}_{0}$, a policy $\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$ is labeled as $\eta$-safe for initial states $\boldsymbol{\mathrm{s}}_{0}\in\mathcal{S}_{0}$ if it satisfies

[TABLE]

In general, there can exist initial states for which a safe policy cannot exist [24], and Definition 1 characterizes policies that preserve safety for a given initial state.

Providing safety guarantees is arguably an open problem when using DNNs as function approximators. However, this problem has been studied in control theory and one successful design technique is robust MPC [25, 26, 27]. Therefore, instead of building the function approximations based on the commonly used DNN approaches, we will use robust MPC, within an extended version of the RL-MPC scheme proposed in [11, 20]. Note that it has been proven in [11] that the optimal policy $\boldsymbol{\mathrm{\pi}}$, value and action-value functions $V_{\star}$ and $Q_{\star}$ can be recovered exactly by function approximations based on MPC, provided that their parametrization is rich enough, even in case the model used in MPC is different from (1).

In order to guarantee safety, robust MPC would ideally rely on the propagation of the dispersion set (16). Unfortunately, this poses severe computational challenges and an auxiliary (time-varying) policy $\boldsymbol{\mathrm{\pi}}^{\mathrm{MPC}}_{\boldsymbol{\mathrm{\theta}},k}$ is preferred in order to recover a computationally tractable formulation [26]. Policy $\boldsymbol{\mathrm{\pi}}^{\mathrm{MPC}}_{\boldsymbol{\mathrm{\theta}},k}$ is typically selected as an open-loop input profile corrected by an affine feedback, see Section III-A. We remark that MPC delivers $\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}=\boldsymbol{\mathrm{\pi}}^{\mathrm{MPC}}_{\boldsymbol{\mathrm{\theta}},0}$. By construction, robust MPC delivers a policy $\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$ which satisfies constraints $\boldsymbol{\mathrm{\xi}}$ at all future times, provided that the SDC (15) holds for all $k$. Even though the dispersion set propagation $\boldsymbol{\mathrm{S}}_{k}^{\boldsymbol{\mathrm{\pi}}^{\mathrm{MPC}}_{\boldsymbol{\mathrm{\theta}},k}}$ is computed based on $\boldsymbol{\mathrm{\pi}}^{\mathrm{MPC}}_{\boldsymbol{\mathrm{\theta}},k}$, the constraints are guaranteed to hold also for $\boldsymbol{\mathrm{S}}_{k}^{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}$.

Since the shape of $\boldsymbol{\mathrm{\hat{S}}}_{+}$ impacts the closed-loop performance, one can let safe RL adapt $\boldsymbol{\mathrm{\hat{S}}}_{+}$. However, this requires one to enforce (15) explicitly in RL. The safe RL problem is then formulated as

[TABLE]

By enforcing (18b), RL is explicitly made aware that some parameter updates are unsafe and, therefore, not feasible. Provided that the SDC holds, MPC delivers a safe policy by construction. Though in principle the SDC has to be enforced for each sample, in Section IV we propose an approach to largely reduce the amount of constraints.

Remark 1

The RL Problem (18) is typically solved using sensitivity-based methods, hence we need to differentiate the function approximator with respect to the parameter $\boldsymbol{\mathrm{\theta}}$. In our case, we need to differentiate the robust MPC problem. This will be detailed in Section III-B.

The proposed safe RL framework performs the following steps: (a) at every time instant MPC is solved and differentiated; the MPC input is applied to the system; state transitions are observed and data is collected to form the sample-based SDC (15); (b) the RL problem (18) is solved (possibly at a lower sampling rate than MPC) and parameter $\boldsymbol{\mathrm{\theta}}$ is updated whenever possible. A scheme is displayed in Figure 1, with reference to the section where each component is discussed.

In this section, we have established the safe RL framework; in the next sections, we will discuss (a) how to implement robust MPC, (b) how to differentiate it in order to be able to solve Problem (18), and (c) how to manage constraint (18b) in a data-efficient fashion. In this paper, we address (a)-(c) by relying on a linear model of (1) to enforce safety. On the one hand, this choice makes the robust constraint satisfaction tractable and not excessively demanding in terms of computations. On the other hand, any nonlinearity present in the system will be accounted for as a perturbation, therefore introducing some conservatism. We remark that nonlinear robust MPC formulations have been developed and can be deployed within the proposed algorithmic framework. The main drawbacks of these formulations are: (a) some form of conservatism cannot be avoided and (b) the computational burden typically becomes prohibitive.

III Robust MPC Based on Invariant Sets

Since guaranteeing robust constraint satisfaction in the general nonlinear case is extremely difficult [26], in the remainder of the paper we focus on the case of an affine model. We describe safety constraints (8) as the inner approximation

[TABLE]

We formulate a $Q$ function approximator based on the classic robust linear MPC [26, 25]:

[TABLE]

where $\boldsymbol{\mathrm{z}}:=(\boldsymbol{\mathrm{x}}_{0},\boldsymbol{\mathrm{u}}_{0},\ldots,\boldsymbol{\mathrm{u}}_{N-1},\boldsymbol{\mathrm{x}}_{N})$. Note that we use $\boldsymbol{\mathrm{x}}$, $\boldsymbol{\mathrm{u}}$ to distinguish the MPC predictions from the actual state and action realizations $\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{a}}$ which define the initial constraint (20i). The dynamic constraints (20j) assume a nominal model without any perturbation. The tube-based approach then treats the system stochasticity, model uncertainties and safety constraint (19) by performing a suitable tightening of the path constraints (20k), i.e., $\boldsymbol{\mathrm{c}}_{k}\geq\boldsymbol{\mathrm{\bar{c}}}$ is used. The terminal constraints (20l) are introduced to guarantee that the path constraints will never be violated at all future times $k>N$.

The value function and optimal policy are obtained as [11]

[TABLE]

In practice $V_{\boldsymbol{\mathrm{\theta}}}$ and $\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$ are computed jointly by solving (20) without enforcing the constraint $\boldsymbol{\mathrm{u}}_{0}=\boldsymbol{\mathrm{a}}$. Parameter $\boldsymbol{\mathrm{\theta}}$ to be adapted by RL may include any of the vector and matrices defining the MPC scheme (20), i.e., generally

[TABLE]

Parameters $H,P$ are typically assumed positive-definite to guarantee the solvability of (20). Parameters $K$ and $\boldsymbol{\mathrm{\theta}}_{\boldsymbol{\mathrm{W}}}$ do not appear explicitly in (20), and are used to compute the constraint tightening and the terminal set, which will be introduced next: we will first present the computation of the constraint tightening, and then discuss the computation of the sensitivities of $V_{\boldsymbol{\mathrm{\theta}}}$, $Q_{\boldsymbol{\mathrm{\theta}}}$, $\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$ with respect to $\boldsymbol{\mathrm{\theta}}$.

III-A Recursive Robust Constraint Satisfaction

In order to guarantee constraint satisfaction for the real system (1) using predictions given by the nominal model (20j), robust MPC explains the difference between predictions and actual state transitions by means of additive noise $\boldsymbol{\mathrm{w}}\in\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$, with $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ satisfying

[TABLE]

We remark that, by using the affine model (20j), conservatism is introduced as any nonlinearity present in (1) will be accounted for by $\boldsymbol{\mathrm{w}}$. Set $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ is parametrized by parameter $\boldsymbol{\mathrm{\theta}}_{\boldsymbol{\mathrm{W}}}$. Common choices in robust MPC are to parametrize $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ using ellipsoids or polytopes; in this paper we consider the latter.

In order to perform constraint tightening, i.e., the computation of $\boldsymbol{\mathrm{c}}_{k}$, we rely on the approach first proposed by [25]. We introduce the prediction error of the nominal model (20j):

[TABLE]

where set $\boldsymbol{\mathrm{E}}_{k}$ predicts an outer approximation of the dispersion set around the predicted trajectory, i.e., $\boldsymbol{\mathrm{S}}_{k}^{\boldsymbol{\mathrm{\pi}}^{\mathrm{MPC}}_{\boldsymbol{\mathrm{\theta}},k}}\subseteq\boldsymbol{\mathrm{x}}_{k}+\boldsymbol{\mathrm{E}}_{k}$, $k=0,\ldots,N$, where $\boldsymbol{\mathrm{\pi}}^{\mathrm{MPC}}_{\boldsymbol{\mathrm{\theta}},k}:=\boldsymbol{\mathrm{a}}_{k}-K\boldsymbol{\mathrm{e}}_{k}$, $\boldsymbol{\mathrm{e}}_{k}\in\boldsymbol{\mathrm{E}}_{k}$ and

[TABLE]

Feedback matrix $K$ is introduced in order to model the fact that any closed-loop strategy will compensate for perturbations on the nominal model. For ease of notation, we define

[TABLE]

Robust constraint satisfaction is then obtained provided that

[TABLE]

such that $\boldsymbol{\mathrm{c}}_{k}$ is obtained by adding the worst-case realization of $C_{K}\boldsymbol{\mathrm{E}}_{k}$ to $\boldsymbol{\mathrm{\bar{c}}}$. For each component $i$ of the path constraint at time $k$ we define

[TABLE]

We lump all components $\boldsymbol{\mathrm{d}}_{i,k}$ in vector $\boldsymbol{\mathrm{d}}_{k}$. Then, constraint satisfaction is obtained for all $\boldsymbol{\mathrm{w}}_{k}\in\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ if

[TABLE]

If $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ is a polytope, then (III-A) can be formulated as an LP; this implies that: (a) constraint tightening is relatively cheap to compute; (b) as detailed in Section IV, the SDC enforcement becomes easier to derive.

In order to guarantee that Problem (20) remains feasible at all times for all $\boldsymbol{\mathrm{w}}_{k}\in\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$, one needs to impose ad-hoc terminal conditions. More specifically, the terminal set $\mathcal{X}_{\mathrm{f}}:=\{\,\boldsymbol{\mathrm{x}}\,|\,G\boldsymbol{\mathrm{x}}+\boldsymbol{\mathrm{g}}\leq 0\,\}$ should be robustly invariant and output admissible, i.e., there exists a terminal control law $\boldsymbol{\mathrm{\kappa}}_{\mathrm{f}}$ such that $\boldsymbol{\mathrm{s}}_{k+1}\in\mathcal{X}_{\mathrm{f}}$ and $C\boldsymbol{\mathrm{s}}_{k}+D\boldsymbol{\mathrm{\kappa_{\mathrm{f}}}}(\boldsymbol{\mathrm{s}}_{k})+\boldsymbol{\mathrm{\bar{c}}}\leq\boldsymbol{\mathrm{0}}$ for every $\boldsymbol{\mathrm{s}}_{k}\in\mathcal{X}_{\mathrm{f}}$. We consider a linear control law $\boldsymbol{\mathrm{\kappa}}_{\mathrm{f}}(\boldsymbol{\mathrm{s}})=-K\boldsymbol{\mathrm{s}}$, coinciding with the linear feedback used to stabilize the prediction error $e$.

In order to compute set $\mathcal{X}_{\mathrm{f}}$, we define

[TABLE]

Note that, by (24)-(III-A), $\boldsymbol{\mathrm{x}}_{0}\in\mathcal{X}_{k}$ implies

[TABLE]

Set $\mathcal{X}_{\infty}$ is Maximal Robust Positive Invariant (MRPI) [28]:

[TABLE]

Additionally, whenever the system is stable and the origin is in the interior of the constraint set, the MRPI set is finitely determined [28, Theorem 6.3], i.e., $\exists\ k^{\prime}<\infty$ s.t. $\mathcal{X}_{k^{\prime}}\equiv\mathcal{X}_{k^{\prime}+i}$, for all $i\geq 1$. The stability requirement further motivates the introduction of feedback through matrix $K$.

As proven in [28], for $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ polyhedral, $G$, $\boldsymbol{\mathrm{\bar{g}}}$ are given by

[TABLE]

where we stress that $\boldsymbol{\mathrm{c}}_{k}=\boldsymbol{\mathrm{\bar{c}}}+\boldsymbol{\mathrm{d}}_{k}$, such that, by using $\boldsymbol{\mathrm{\kappa}}_{\mathrm{f}}=-K\boldsymbol{\mathrm{e}}$ one can spare a large amount of computations, $\boldsymbol{\mathrm{\bar{g}}}$ being already computed. The condition $G\boldsymbol{\mathrm{x}}_{N}+\boldsymbol{\mathrm{\bar{g}}}\leq 0$ would then guarantee robust constraint satisfaction for all future times if $\boldsymbol{\mathrm{s}}_{N}=\boldsymbol{\mathrm{x}}_{N}$. However, $\boldsymbol{\mathrm{s}}_{N}=\boldsymbol{\mathrm{x}}_{N}+\boldsymbol{\mathrm{e}}_{N}$, such that also the terminal constraint (20l) must be tightened. Analogously to the case of path constraints, we define $\boldsymbol{\mathrm{g}}=\boldsymbol{\mathrm{\bar{g}}}+\boldsymbol{\mathrm{h}}$ with

[TABLE]

Remark 2

Typically, constraints that can never become active are removed from (34), so as to reduce the dimension.

Safety is guaranteed by the following result on tube MPC.

Proposition 1** (Recursive Feasibility)**

Assume that $\mathcal{X}_{\mathrm{f}}:=\{\,\boldsymbol{\mathrm{x}}\,|\,G\boldsymbol{\mathrm{x}}+\boldsymbol{\mathrm{g}}\leq 0\,\}$ is RPI and Problem (20) is feasible at time $k=0$. Then Problem (20) is feasible for all $\boldsymbol{\mathrm{w}}_{k}\in\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ and all times $k\geq 0$. If moreover (23) holds, the real system (1) satisfies the safety constraint (19) at all times $k\geq 0$.

Proof:

The proof can be found in, e.g., [25]. ∎

This proof is valid in case $\boldsymbol{\mathrm{\theta}}$ is kept fixed. A discussion on how to enforce recursive feasibility upon updates of $\boldsymbol{\mathrm{\theta}}$ is proposed in Section V, Remark 9.

While the framework of robust linear MPC based on MRPI sets is well established, the computation of the parametric sensitivities of an MPC problem required to deploy most RL methods is not common. In particular, in the case of a tube-based formulation, also the constraint tightening procedure needs to be differentiated. This is also not common and deserves to be discussed in detail. We therefore devote the next subsection to the computation of the derivative of the value and action-value function with respect to parameter $\boldsymbol{\mathrm{\theta}}$.

III-B Differentiability

In order to be able to deploy RL algorithms to adapt parameter $\boldsymbol{\mathrm{\theta}}$, we need to be able to differentiate the MPC scheme and, therefore, the constraint definition with respect to $\boldsymbol{\mathrm{\theta}}$. In principle, $\boldsymbol{\mathrm{\theta}}$ could include $A,B,K$, but also all other parameters of the MPC formulation (22). In order to compute $\nabla_{\boldsymbol{\mathrm{\theta}}}\boldsymbol{\mathrm{c}}_{k}$, $\nabla_{\boldsymbol{\mathrm{\theta}}}C_{N}$, $\nabla_{\boldsymbol{\mathrm{\theta}}}\boldsymbol{\mathrm{c}}_{N}$ one can use results from parametric optimization to obtain the following lemmas [29].

Consider a parametric NLP with cost $\phi_{\boldsymbol{\mathrm{\theta}}}$, primal-dual variable $\boldsymbol{\mathrm{y}}$ and parameter $\boldsymbol{\mathrm{\theta}}$. We refer to [30] for the definition of Lagrangian $l_{\boldsymbol{\mathrm{\theta}}}^{0}(\boldsymbol{\mathrm{y}})$, KKT conditions, Linear Independence Constraint Qualification (LICQ), Second-Order Sufficient Conditions (SOSC) and Strict Complementarity (SC). For a fixed active set the KKT conditions reduce to the equality $\boldsymbol{\mathrm{r}}^{0}_{\boldsymbol{\mathrm{\theta}}}(\boldsymbol{\mathrm{y}})=0.$

Lemma 1

Consider a parametric optimization problem with optimal primal-dual solution $\boldsymbol{\mathrm{y}}^{\star}$. Assume that LICQ, SOSC and SC hold at $\boldsymbol{\mathrm{y}}^{\star}$. Then, the following holds

[TABLE]

Proof:

The result can be found in, e.g., [29]. ∎

Corollary 1

Assume that LICQ, SOSC and SC hold at the optimal solution of (20). Then, the value function $V_{\boldsymbol{\mathrm{\theta}}}$, action-value function $Q_{\boldsymbol{\mathrm{\theta}}}$ and optimal solution $\boldsymbol{\mathrm{y}}^{\star}$ (therefore also policy $\boldsymbol{\mathrm{\pi}}$) are differentiable with respect to parameter ${\boldsymbol{\mathrm{\theta}}}$, with:

[TABLE]

where $l_{\boldsymbol{\mathrm{\theta}}}$ is the Lagrangian of Problem (20), $\bar{l}_{\boldsymbol{\mathrm{\theta}}}$ is the Lagrangian when constraint $\boldsymbol{\mathrm{u}}_{0}=\boldsymbol{\mathrm{a}}$ is eliminated, and $\boldsymbol{\mathrm{r}}_{\boldsymbol{\mathrm{\theta}}}$ denotes the KKT conditions for the optimal active set.

Remark 3

When solving an LP, QP or NLP using a second-order method, e.g., active-set or interior-point, the most expensive operation is the factorization of the KKT matrix $\frac{\partial\boldsymbol{\mathrm{r}}_{\boldsymbol{\mathrm{\theta}}}}{\partial\boldsymbol{\mathrm{y}}}$. Once the matrix is factorized, the solution of the linear system is computationally inexpensive. Therefore, the sensitivities of the solution are in general much cheaper to evaluate than solving the problem itself. The sensitivity of the optimal value function is even simpler to compute, since it consists in the differentiation of the Lagrangian, see (36a).

In (36), $\boldsymbol{\mathrm{r}}_{\boldsymbol{\mathrm{\theta}}}$, $l_{\boldsymbol{\mathrm{\theta}}}$, and $\bar{l}_{\boldsymbol{\mathrm{\theta}}}$ depend on $\boldsymbol{\mathrm{c}}_{k}$, $\boldsymbol{\mathrm{g}}$ which, in turn, depend on $\boldsymbol{\mathrm{\theta}}$ as they are optimal values of parametric optimization problems (III-A) and (35). Consequently, one needs to evaluate $\nabla_{\boldsymbol{\mathrm{\theta}}}\boldsymbol{\mathrm{c}}_{k}$, $\nabla_{\boldsymbol{\mathrm{\theta}}}\boldsymbol{\mathrm{g}}$. In the following, we further detail the application of Lemma 1 to this case.

We consider only Problem (III-A), since the derivation for (35) is analogous. First, we state separability and, therefore, parallelizability of the computation of $\boldsymbol{\mathrm{d}}_{k}$ in the following Lemma.

Lemma 2

Each component of $\boldsymbol{\mathrm{d}}_{k}$ can be computed as

[TABLE]

Proof:

Each term in the sum $\sum_{j=0}^{k-1}A_{K}^{j}\boldsymbol{\mathrm{w}}_{j}$ depends only on variable $\boldsymbol{\mathrm{w}}_{j}$, and the problem is fully separable. ∎

Then, Lemma 1 can be applied to obtain

[TABLE]

where $l^{\boldsymbol{\mathrm{d}}}_{\boldsymbol{\mathrm{\theta}},i,j}$ is the Lagrangian of Problem (37). Provided that a second-order method is used for solving Problem (37), then the matrix factorization is available and can be reused to compute the sensitivities at a negligible cost. For any function $\boldsymbol{\mathrm{f}}(\boldsymbol{\mathrm{c}}_{k}(\boldsymbol{\mathrm{\theta}}))$, the chain rule yields

[TABLE]

Remark 4

As underlined above, Problems (37) can be solved in parallel not only for each prediction time $k$, but also for each component $i$. Moreover, if an active-set solver is used, the active set can be initialized and the matrix factorization reused, such that often there will be no need for recomputing the factorization and computations can be done in an extremely efficient manner. Finally, Problems (37) are very low dimensional and are therefore solved extremely quickly.

III-C Guaranteeing MPC Feasibility and LICQ

We detail next two issues that can be easily encountered when deploying RL based on MPC (20). Since these are common, a simple solution that has become a standard in MPC is readily available.

MPC Feasibility

Since the set of possible perturbations is not known a priori but rather approximated as $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ based on the collected samples, it cannot be excluded that some future sample $\boldsymbol{\mathrm{w}}_{k}$ will not be inside the set, i.e., $\boldsymbol{\mathrm{w}}_{k}\notin\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$. The set approximation must then be modified to include the new sample (see Section IV-B and (48)), but recursive feasibility is potentially lost, no action is computed and the controller stops working.

Sensitivity Computation

The sensitivity computation is valid only if LICQ holds. However, Problem (20) is not guaranteed to satisfy LICQ.

We propose to address both issues by using a common approach in MPC, i.e., a constraint relaxation for Constraints (20k) and (20l) with an exact penalty [31]: variables $\boldsymbol{\mathrm{\sigma}}_{k}$ are introduced, the constraints are modified as

[TABLE]

and the term $\sum_{k=1}^{N}\boldsymbol{\mathrm{\rho}}^{\top}\boldsymbol{\mathrm{\sigma}}_{k}$ is added to the cost.

Remark 5

Constraints only involving the controls do not need to be relaxed, while any constraint involving the states should not be imposed at $k=0$, since then LICQ can not be guaranteed even if the relaxation proposed above is deployed.

We formalize the result in the next proposition.

Proposition 2

Assume that Constraints (20k) and (20l) are relaxed as per (38) and the term $\sum_{k=1}^{N}\boldsymbol{\mathrm{\rho}}^{\top}\boldsymbol{\mathrm{\sigma}}_{k}$ is added to the cost. Then, if $\boldsymbol{\mathrm{\rho}}<\infty$ is large enough, the solution is unchanged whenever feasible and recursive feasibility and LICQ are guaranteed.

Proof:

The first result, i.e., solution equivalence whenever feasible and recursive feasibility is well-known [31] and [32, Theorem 14.3.1]. Regarding LICQ, we first note that in a formulation without Constraints (20k) and (20l) LICQ holds by construction: Constraints (20i) and (20j) can be eliminated by condensing [33], yielding a problem with $Nn_{\boldsymbol{\mathrm{u}}}$ unconstrained variables. The introduction of any linearly independent set of pure control constraints does then not jeopardize LICQ by construction. Assume now to introduce a linearly-dependent constraint of the form (20k) or (20l) with Jacobian $\boldsymbol{\mathrm{\nu}}$. By introducing slack variable $\sigma$, the new Jacobian becomes \left[\begin{array}[]{cc}\boldsymbol{\mathrm{\nu}}&1\end{array}\right], which is by construction linearly independent with \left[\begin{array}[]{cc}\boldsymbol{\mathrm{\nu}}&0\end{array}\right] and, consequently, with the other constraints in the problem. ∎

As discussed in Section II, safety holds with probability $\mathbb{P}\left[\,\eqref{eq:sdc}\,|\,\boldsymbol{\mathrm{\theta}}\,\right]=1$. This is an intrinsic issue of any safety-enforcing control scheme, and the proposed relaxation only solves the issue of avoiding infeasibility for the MPC scheme. This is further discussed in Section V-B.

IV Safe Design Constraint and Data Management

As explained in the previous section, safety is obtained if all possible state transitions are correctly captured in the MPC formulation, i.e., if $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ and, therefore, $\boldsymbol{\mathrm{g}}_{\boldsymbol{\mathrm{\theta}}}$ is correctly identified. In other words, based on the collected data the SDC must be enforced in order to guarantee that the uncertainty described by $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ is representative of the real system (1). In this section we propose a sample-based formulation of the SDC (11) to be used within RL formulations (5).

Many control systems are typically operated at high sampling rates and, consequently, data is collected at high rates. In order to be able to deal in real-time with the large amounts of data that cumulate, it is necessary to (a) retain only strictly relevant data, and (b) compress the available information using appropriately defined data structures. In the following, we first discuss the data structures involved in RL-MPC and then discuss how to make an efficient use of data in the context of the proposed MPC formulation.

IV-A Set Membership and SDC

Consider the (possibly very large) set of state transitions observed on the real system:

[TABLE]

The problem of enforcing the SDC (15) is related to the one of estimating the dispersion set $\boldsymbol{\mathrm{S}}_{+}$ from data, which has been studied in the context of set-membership system identification (SMSI) [34]. Essentially, the dispersion set must satisfy

[TABLE]

i.e., $\boldsymbol{\mathrm{\theta}}$ must satisfy the SDC (15), and therefore belongs to set

[TABLE]

We should underline here the difference between $\boldsymbol{\mathrm{\hat{S}}}_{+}$ defined in (10) and $\mathcal{S}_{\mathcal{D}}$. The former is an outer approximation of the set of all realizations $\boldsymbol{\mathrm{s}}_{+}$, given state and action $\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{a}}$, parametrized by $\boldsymbol{\mathrm{\theta}}$. The latter instead describes the set of parameters $\boldsymbol{\mathrm{\theta}}$ such that all state transitions from $\mathcal{D}$ are contained in $\boldsymbol{\mathrm{\hat{S}}}_{+}$.

In SMSI parameter $\boldsymbol{\mathrm{\theta}}$ is typically selected so as to obtain the smallest possible set $\boldsymbol{\mathrm{S}}_{+}$. In the absence of specific information on the control task to be executed this is arguably a very reasonable approach. However, given a specific control task to be executed, better performance might be obtained by selecting parameter $\boldsymbol{\mathrm{\theta}}$ to approximate some part of $\boldsymbol{\mathrm{S}}_{+}$ accurately even at the cost of increasing the volume of $\boldsymbol{\mathrm{\hat{S}}}_{+}$. We therefore provide the following definition.

Definition 2** (Set Membership Optimality)**

Equivalently to (2), we define the closed-loop cost $J(\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}})$ associated with the policy $\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$ learned by RL when relying on the set $\boldsymbol{\mathrm{\hat{S}}}_{+}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{a}},\boldsymbol{\mathrm{\theta}})$. Then, parameter $\boldsymbol{\mathrm{\theta}}$ is (locally) optimal for the control task iff $\nexists\ \boldsymbol{\mathrm{\bar{\theta}}}\ \in\mathcal{B}_{\epsilon}(\boldsymbol{\mathrm{\theta}})\ \mathrm{s.t.}\ J(\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\bar{\theta}}}})<J(\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}})$, where $\mathcal{B}_{\epsilon}(\cdot)$ denotes a ball of radius $\epsilon$ centered at $\boldsymbol{\mathrm{\theta}}$.

In other words, the optimal parameter minimizes the cost subject to the SDC (15). Based on this definition, the optimal set approximation is obtained by an SMSI which selects $\boldsymbol{\mathrm{\theta}}$ to maximize the closed-loop performance of the policy $\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$. In this paper, we aim at doing so by means of RL.

Every time the RL problem (18) is solved, one must ensure that $\boldsymbol{\mathrm{\theta}}\in\mathcal{S}_{\mathcal{D}}$, i.e., $|\mathcal{D}|$ constraints need to be included in the problem formulation. With large amounts of data, this could make the problem computationally intractable. In the following, we will analyze how to tackle this issue.

While the previous developments did not require any assumption on function $\boldsymbol{\mathrm{g}}_{\boldsymbol{\mathrm{\theta}}}$, in order to efficiently manage data we will assume hereafter that $\boldsymbol{\mathrm{g}}_{\boldsymbol{\mathrm{\theta}}}$ is affine in $\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{a}}$. Note that this choice is consistent with the choice of an affine model. The SDC then becomes

[TABLE]

for some $M,\boldsymbol{\mathrm{m}}$ possibly part of $\boldsymbol{\mathrm{\theta}}$. Since both $\mathcal{S}_{\mathcal{D}}$ and $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ are defined based on $\boldsymbol{\mathrm{g}}_{\boldsymbol{\mathrm{\theta}}}$, parameters $M,\boldsymbol{\mathrm{m}}$ for the two sets coincide. Therefore, the SDC directly defines the uncertainty set $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ used by MPC to compute safe policies. Note that $A$, $B$, $\boldsymbol{\mathrm{b}}$, $M$, $\boldsymbol{\mathrm{m}}$ can all be included in $\boldsymbol{\mathrm{\theta}}$ and, therefore, adapted by RL. Which parameter to adapt or keep constant is a design choice, see also Remark 7.

IV-B Model-Based Data Compression

In this section we discuss how to compress the available data without loss of information to significantly reduce the complexity of safe RL. It will become clear that the nominal model (20j) plays a very important role in this context, as it makes it possible to organize data such that it can be efficiently exploited. Unfortunately, this efficiency is lost if the model parameters are updated. Approaches to circumvent this issue can be devised and are the subject of ongoing research. We begin the analysis by providing the following definition.

Definition 3** (Optimal Data Compression)**

Given the selected parametrization and MPC formulation, an optimal data compression selects a dataset $\mathcal{\bar{D}}\subseteq\mathcal{D}$ such that:

[TABLE]

Hence an optimal data compression retains the minimum amount of data required to represent the set $\mathcal{S}_{\mathcal{D}}$. In the following, we assume that $A$, $B$, $\boldsymbol{\mathrm{b}}$ are fixed and exploit the model to achieve an optimal data compression.

The introduction of the nominal model (20j) allows us to restructure the data and only store the noise $\mathcal{W}:=\{\,\boldsymbol{\mathrm{w}}_{0},\ldots,\boldsymbol{\mathrm{w}}_{n}\,\}$, obtained from

[TABLE]

By using (43), we rewrite the SDC as

[TABLE]

with $\boldsymbol{\mathrm{\theta}}_{\boldsymbol{\mathrm{W}}}=(M,\boldsymbol{\mathrm{m}})$ a component of $\boldsymbol{\mathrm{\theta}}$. By exploiting the model we can therefore reduce the dimension of the space of the dataset from $2n_{\boldsymbol{\mathrm{s}}}+n_{\boldsymbol{\mathrm{a}}}$ to $n_{\boldsymbol{\mathrm{s}}}$, since $(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{a}},\boldsymbol{\mathrm{s}}_{+})\in\mathbb{R}^{n_{\boldsymbol{\mathrm{s}}}\times n_{\boldsymbol{\mathrm{a}}}\times n_{\boldsymbol{\mathrm{s}}}}$ and $\boldsymbol{\mathrm{w}}\in\mathbb{R}^{n_{\boldsymbol{\mathrm{s}}}}$. Note that this is only possible if one neglects the dependence of $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ (and therefore $\mathcal{W}$) on $\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{a}}$.

The dispersion set approximation related with (44) is then

[TABLE]

such that (III-A), (35) used in constraint tightening are LPs. Additionally, it is cheap to (a) evaluate if $\boldsymbol{\mathrm{s}}_{+}\in\boldsymbol{\mathrm{\hat{S}}}_{+}$ by using (43) and (b) verifying that $\boldsymbol{\mathrm{\theta}}\in\mathcal{S}_{\mathcal{W}}$, i.e., $M\boldsymbol{\mathrm{w}}_{k}\leq\boldsymbol{\mathrm{m}}$ holds, since both operations only require few matrix-vector operations. However, while (a) requires a single evaluation of the inequality in (45), (b) requires to evaluate the inequality for each data point in (44).

The use of the model and the convexity assumption make it possible to further compress data: any point in the interior of the convex hull of set $\mathcal{W}$ does not provide any additional information regarding (44), such that the convex hull

[TABLE]

carries all necessary information. Constructing the convex hull facet representation can be a rather expensive operation, which can in general not be done online. However, checking whether a sample $\boldsymbol{\mathrm{w}}_{k}$ lies inside the convex hull $\mathcal{\bar{W}}$ of a set of samples $\mathcal{W}$ can be done without building the facet representation of the convex hull. To this end, we define the LP

[TABLE]

and exploit the following result.

Proposition 3

Assume that $\boldsymbol{\mathrm{0}}\in\mathrm{Conv}(\mathcal{W})$, then

[TABLE]

with $\zeta$ from (46).

Proof:

The definition of convex hull implies $\boldsymbol{\mathrm{\hat{w}}}\in\mathrm{Conv}(\mathcal{W})\Rightarrow\exists\ \boldsymbol{\mathrm{z}}\geq 0,\|\boldsymbol{\mathrm{z}}\|_{1}=1\ \mathrm{s.t.}\ \boldsymbol{\mathrm{\hat{w}}}=\sum_{i=1}^{|\mathcal{\bar{W}}|}z_{i}\boldsymbol{\mathrm{w}}_{i}$, which proves $\zeta\leq 1\Leftarrow\boldsymbol{\mathrm{\hat{w}}}\in\mathrm{Conv}(\mathcal{W})$. This also covers the implication $\zeta=1\Rightarrow\boldsymbol{\mathrm{\hat{w}}}\in\mathrm{Conv}(\mathcal{W})$. The implication $\zeta<1\Rightarrow\boldsymbol{\mathrm{\hat{w}}}\in\mathrm{Conv}(\mathcal{W})$ is proven by noting that $\zeta\boldsymbol{\mathrm{w}}_{i}\in\mathrm{Conv}(\mathcal{W}),\ \forall\ \zeta\in[0,1]$. The case $\zeta<1$ is thus immediately reconducted to the case $\zeta=1$. ∎

Formulation (46) is very convenient, as it is always feasible, provided that $\mathcal{W}$ spans the full space $\mathbb{R}^{n_{\boldsymbol{\mathrm{w}}}}$, which is a minimum reasonable requirement in this context, since it guarantees that every vector in $\mathbb{R}^{n_{\boldsymbol{\mathrm{w}}}}$ can be obtained as a linear combination of $\boldsymbol{\mathrm{w}}_{i}$. Note also that the assumption $\boldsymbol{\mathrm{0}}\in\mathrm{Conv}(\mathcal{W})$ is a rather mild requirement on the accuracy of the model, as it always holds when standard system identification techniques are deployed to estimate the model parameters $A,B,\boldsymbol{\mathrm{b}}$.

We prove the efficiency of the convex hull approach in the following theorem.

Theorem 1** (Convex Hull Optimality)**

Given the dispersion set $\boldsymbol{\mathrm{\hat{S}}}_{+}$ defined in (45), the convex hull $\mathcal{\bar{W}}$ of the state transition noise $\mathcal{W}$ is an optimal data compression.

Proof:

The definition of $\mathcal{\bar{W}}$ implies that any point in $\mathcal{W}$ can be obtained as the convex combination of points in $\mathcal{\bar{W}}$. By definition we have that $\boldsymbol{\mathrm{g}}_{\boldsymbol{\mathrm{\theta}}}^{\boldsymbol{\mathrm{w}}}(\boldsymbol{\mathrm{w}}_{1})\leq 0$, $\boldsymbol{\mathrm{g}}_{\boldsymbol{\mathrm{\theta}}}^{\boldsymbol{\mathrm{w}}}(\boldsymbol{\mathrm{w}}_{2})\leq 0$, for all $\boldsymbol{\mathrm{w}}_{1},\boldsymbol{\mathrm{w}}_{2}\in\mathcal{W}$. For $\beta\in[0,1]$ we define $\boldsymbol{\mathrm{w}}_{\beta}:=\beta\boldsymbol{\mathrm{w}}_{1}+(1-\beta)\boldsymbol{\mathrm{w}}_{2}$ such that

[TABLE]

since convexity of $\mathcal{\bar{W}}$ is equivalent to convexity of $\boldsymbol{\mathrm{g}}_{\boldsymbol{\mathrm{\theta}}}^{\boldsymbol{\mathrm{w}}}$. This entails that $\boldsymbol{\mathrm{w}}_{\beta}\in\mathcal{W}$, such that any set $\boldsymbol{\mathrm{\hat{S}}}_{+}$ computed using $\mathcal{\bar{W}}$ satisfies

[TABLE]

i.e., $\forall\,\boldsymbol{\mathrm{w}}\in\mathcal{W}$. Finally, by removing any data point from $\mathcal{\bar{W}}$ one cannot guarantee that $\mathcal{S}_{\mathcal{\bar{W}}}=\mathcal{S}_{\mathcal{W}}$, such that $\mathcal{\bar{W}}$ solves (42). ∎

In order to construct the convex hull $\mathcal{\bar{W}}$, once a new sample $\boldsymbol{\mathrm{w_{k}}}$ is available, we check whether $\boldsymbol{\mathrm{w}}_{k}\in\mathcal{\bar{W}}$ by solving (46). In case $\zeta\leq 1$, $\mathcal{\bar{W}}$ does not need to be updated; otherwise, the new point is added to $\mathcal{W}$, such that $\mathcal{\bar{W}}$ is implicitly updated. Note that the LP (46) has linear complexity in the amount of vertices $|\mathcal{\bar{W}}|$.

Remark 6

As data are collected, $|\mathcal{D}|$ becomes indefinitely large. Even though $\mathcal{\bar{W}}$ is optimal, also $|\mathcal{\bar{W}}|$ can become indefinitely large, though arguably at a lower rate than $|\mathcal{D}|$. This is a fundamental issue of any sample-based SDC or set membership approach, unless additional assumptions are introduced. Simple strategies such as limiting the maximum amount of vertices/facets to be used for an approximation $\mathcal{\hat{W}}$ of the convex hull could be devised, at the price of introducing conservatism. Such investigations are beyond the scope of this paper and require future research.

Remark 7

One must take extra care if the model parameters $A,B$ are updated. Indeed, if one updates $A,B,\boldsymbol{\mathrm{b}}$ to $\tilde{A},\tilde{B},\boldsymbol{\mathrm{\tilde{b}}}$, then the new noise satisfies:

[TABLE]

such that the noise update

[TABLE]

is state-action dependent for all $\tilde{A}\neq A$, $\tilde{B}\neq B$. Therefore, any change in those parameters requires one to recompute the noise vector $\boldsymbol{\mathrm{w}}_{k}$ for all recorded state-action pairs. Updates in parameter $\boldsymbol{\mathrm{b}}$ instead can be performed without much complication and simply entail a shift of the noise, which is state-action independent, such that $\mathcal{\tilde{W}}=\mathcal{W}+\boldsymbol{\mathrm{b}}-\boldsymbol{\mathrm{\tilde{b}}}$.

IV-C Further Observations on the Sample-Based SDC

The convex hull $\mathcal{\bar{W}}$ is the smallest set encompassing all observed samples, i.e., it is the smallest set satisfying the SDC (15). Therefore, it is optimal both in terms of volume and in terms of cost, i.e., in the sense of Definition 2. Hence, one might be tempted to select $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}=\mathcal{\bar{W}}$ to reduce conservatism in MPC (20) as much as possible. However, $\mathcal{\bar{W}}$ is typically composed of a very large amount of facets, which renders the constraint tightening procedure very costly and results in a terminal constraint of high dimension. In practice, a set $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ of fixed and low complexity is preferred, hence the importance of enforcing set membership optimality as per Definition 2, i.e., through RL.

Thus far we have not discussed how the set $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ is represented. However, the choice of representation becomes important when dealing with large amounts of data. Moreover, the parametrization of $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ should be selected consistently with the algorithm used for solving the robust MPC Problem. Convex polytopes can be parametrized using the so-called (a) facet or (b) vertex representation. In case (a), the set is parametrized as $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}:=\{\,\boldsymbol{\mathrm{w}}\,|\,M\boldsymbol{\mathrm{w}}\leq\boldsymbol{\mathrm{m}}\,\}$ with parameter $\boldsymbol{\mathrm{\theta}}_{\boldsymbol{\mathrm{W}}}=(M,\boldsymbol{\mathrm{m}})$. In case (b), the set is parametrized as $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}:=\{\,\boldsymbol{\mathrm{w}}\,|\,\boldsymbol{\mathrm{w}}\in\mathrm{Conv}(\{\boldsymbol{\mathrm{v}}_{0},\ldots,\boldsymbol{\mathrm{v}}_{m}\})\,\}$ with parameter $\boldsymbol{\mathrm{\theta}}_{\boldsymbol{\mathrm{W}}}=(\boldsymbol{\mathrm{v}}_{0},\ldots,\boldsymbol{\mathrm{v}}_{m})$, i.e., the vertices the polytope.

Differently from $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$, for the convex hull $\mathcal{\bar{W}}$ we use the vertex representation. This allows a simpler and less computationally demanding construction and incremental update of $\mathcal{\bar{W}}$. This advantage, however, results in a more costly evaluation of $\boldsymbol{\mathrm{w}}_{k}\in\mathcal{\bar{W}}$ with respect to a facet representation. Finally, this makes it simple to enforce the SDC (44), which becomes

[TABLE]

The question on which representation is the most convenient for the convex hull $\mathcal{\bar{W}}$ is still open and will be further investigated in future research, which will also consider a combination of both the facet and vertex representation. We provide next some fundamental observations regarding $\mathcal{\bar{W}}$, its cardinality and its relationship with the nominal model.

Remark 8

As a consequence of Fundamental Limitation 1, in any sample-based context it is possible that a new sample $\boldsymbol{\mathrm{w}}_{k}$ falls out of the convex hull of previous samples, i.e., $\boldsymbol{\mathrm{w}}_{k}\notin\mathcal{\bar{W}}$, such that potentially $M_{i}\boldsymbol{\mathrm{w}}_{k}\geq\boldsymbol{\mathrm{m}}_{i}$ for some $i$. In this case, one needs to instantaneously adapt the SDC (15). A straightforward adaptation of $\mathcal{S}_{\mathcal{\bar{W}}}$ is obtained as

[TABLE]

This enlargement of the uncertainty set $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ entails an enlargement of the dispersion set, such that the constraints will be further tightened. This can jeopardize recursive feasibility of MPC (20) such that it is necessary to deploy the constraint relaxation proposed in Section III-C.

V Safe RL MPC Implementation

After having introduced all necessary components of the RL-MPC scheme, in this section we focus on the safe RL problem, discuss more in detail the RL problem and present some open research questions.

Most recursive RL approaches use only the current sample and update the parameter recursively as

[TABLE]

with $\boldsymbol{\mathrm{\theta}}^{*}=\boldsymbol{\mathrm{\theta}}_{k}+\nabla_{\boldsymbol{\mathrm{\theta}}}\psi(\boldsymbol{\mathrm{s}}_{k+1},\boldsymbol{\mathrm{s}}_{k},\boldsymbol{\mathrm{a}}_{k},\boldsymbol{\mathrm{\theta}})$ and $\psi$ defined, e.g., as per (6) or (7).

This means that the update is the solution (or a step of stochastic gradient descent) of an unconstrained optimization problem. Since we need to enforce the SDC, the safe RL problem (18) yielding $\boldsymbol{\mathrm{\theta}}^{*}$ is constrained. By relying on the developments of Section IV, we formulate (18) as

[TABLE]

Problem (50a) (without constraints) can be seen as a standard RL formulation [20]. Positive-definite constraints (50b) are introduced to make sure that MPC is properly formulated and easily and efficiently solvable. The SDC (15) is imposed in (50c). If the model parameters $A$, $B$ are updated, the SDC cannot be implemented as (15), but rather as $\boldsymbol{\mathrm{\theta}}\in\mathcal{S}_{\mathcal{D}}$, with $\mathcal{S}_{\mathcal{D}}$ given by (41).

Remark 9

Any update in parameter $\boldsymbol{\mathrm{\theta}}$ results in a modification of set $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$, such that the existence of a solution to the MPC problem could be jeopardized. Several options can be envisioned, including: (a) updating $\boldsymbol{\mathrm{\theta}}$ only when feasible, (b) reducing the step size until feasibility is recovered, (c) enforcing feasibility as an additional constraint in (50). Any of the strategies (a)-(c) guarantees the feasibility of the updated robust MPC scheme since, provided that MPC is correctly formulated, by Proposition 1 initial feasibility entails recursive feasibility. In turn, this entails safety of the parameter update. It is worth mentioning that we did not encounter feasibility issues in the simulations we performed. Nevertheless, future research will investigate this issue in depth.

Remark 10

Both $Q_{\boldsymbol{\mathrm{\theta}}}$ and $J(\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}})$ depend on the constraint tightening procedure, such that their first-order sensitivities can be discontinuous. Consequently, also the first-order sensitivities of $\psi$ are non-smooth. This is the case when SC does not hold either in the MPC problem (20) or in the constraint tightening problems (III-A), (35), i.e., when some constraint(s) are weakly active, such that infinitesimal perturbations could cause an active-set change. In principle, this could create problems to algorithms for continuous optimization. However, the set on which SC does not hold has zero measure, such that the probability that one sample falls onto one of these points is zero, and the RL solution is unaffected.

Since the main concern of this paper is safety, we present next how to guarantee safety also during exploration, i.e., when the action applied to the system is not given by (20). Afterwards, we will further discuss open research questions and possible ways of addressing them.

V-A Safe Exploration

The proposed formulation guarantees safety during exploitation, i.e., whenever the policy is given by (21). However, specific care needs to be taken during exploration, i.e., when the optimal policy is perturbed, since also in this phase constraint satisfaction must not be jeopardized. Exploration is typically performed by picking a random action among all feasible actions with a given probability distribution. The main difficulty in enforcing safety is related to the complexity of the set of safe actions

[TABLE]

Since this set is implicitly approximated within robust MPC, this issue can be tackled by using a modified version of the robust MPC Problem (20):

[TABLE]

with either $f\left(\boldsymbol{\mathrm{u}}_{0},\boldsymbol{\mathrm{q}}\right):=\rho\|\boldsymbol{\mathrm{u}}_{0}-\boldsymbol{\mathrm{q}}\|$, or $f\left(\boldsymbol{\mathrm{u}}_{0},\boldsymbol{\mathrm{q}}\right):=\boldsymbol{\mathrm{q}}^{\top}\boldsymbol{\mathrm{u}}_{0}$, and a randomly chosen $\boldsymbol{\mathrm{q}}$. Note that performing constrained exploration in policy gradient methods can produce a biased gradient estimation. Simple strategies to tackle this issue are proposed in [23].

Theorem 2

Consider a convergent RL scheme solving Problem (50) where (a) robust MPC (20) is used as function approximator; (b) available data is handled as detailed in Section IV; and (c) exploration is performed according to (51). Assume that all new data yields $\boldsymbol{\mathrm{w}}_{k}\in\mathcal{\bar{W}}$. Then, the RL scheme is safe in the sense of Definition 1, i.e., $\mathbb{P}\left[\exists\,k\,\mathrm{s.t.}\,\boldsymbol{\mathrm{\xi}}(\boldsymbol{\mathrm{s}}_{k},\boldsymbol{\mathrm{a}}_{k})>0\right]\leq 1-\eta(\mathcal{D})$. Additionally, the scheme is optimal in the sense of Definition 3.

Proof:

We observe that $\eta\left(\mathcal{D}\right)$ quantifies the probability that it is possible that a scenario is unaccounted for in the construction of $\boldsymbol{\mathrm{\hat{S}}}_{+}$. Moreover, we observe that if $\boldsymbol{\mathrm{\hat{S}}}_{+}$ does not account for all possible scenarios, then $\boldsymbol{\mathrm{w}}_{k}\notin\mathcal{\bar{W}}$ can happen for some $k$. As a result, one can verify that

[TABLE]

We note that, by Proposition 1, $\boldsymbol{\mathrm{w}}_{k}\in\mathcal{\bar{W}}$ implies that robust MPC (20) is recursively feasible. Moreover, exploration is performed using (51), i.e., (20) with a modified cost, such that recursive feasibility is also preserved by Proposition 1. This entails that a constraint violation can only occur if $\boldsymbol{\mathrm{w}}_{k}\notin\mathcal{\bar{W}}$. The converse, however, is not true in general, since there might exist $\boldsymbol{\mathrm{w}}_{k}\notin\mathcal{\bar{W}}$ which does not entail a constraint violation. Consequently, the probability that the constraints be violated at any time is bounded from above by $1-\eta(\mathcal{D})$:

[TABLE]

Data compression optimality (Definition 3) is a direct consequence of Theorem 1. ∎

Remark 11

We ought to stress that assuming that RL converges to a minimum of $\mathbb{E}[\psi]$ is standard. However, assuming convergence to a local minimum of $J$ can be a strong assumption, which is typically not met by $Q$-learning and SARSA, while actor-critic methods are typically expected to converge to a local minimum of $J$ for the given parametrization. Optimality in the sense of Definition 2 (set membership optimality) could then be claimed for convergent actor-critic methods, while it is reasonable to expect a certain degree of suboptimality with $Q$-learning.

We provide an overview of the computations performed by RL-MPC in Algorithm 1. Note that we introduced lines 8 and 10 since the RL problem and, consequently, the recomputation of the constraint tightening and terminal set might in principle be updated less often than the feedback sampling time. In that case, a batch RL approach could be used instead of a recursive one and the real-time requirements would only concern MPC. Additionally, we stress that in our Matlab implementation which was only partially made efficient, the constraint tightening procedure on line 11 required approximately twice the computation time of MPC on lines 2 and 4. The computation of the RL step on line 9 was approximately $5$ times faster.

V-B Discussion on the Proposed Approach

We discuss next a few open research questions and comment on possible extensions to the proposed framework.

V-B1 Safety

Our safety definition is based on the assumption that all future samples satisfy $\boldsymbol{\mathrm{w}}_{k}\in\mathcal{\bar{W}}$. In practice this assumption, though standard in robust MPC and SMSI, can be rather strong. However, this is a fundamental problem of safety: one can never guarantee a priori that a bounded set contains all possible future realizations of an unknown stochastic process. In Section III-C we have proposed a simple and practical approach to retain MPC feasibility even in the case $\boldsymbol{\mathrm{w}}_{k}\notin\mathcal{\bar{W}}$. However, with this approach safety is potentially lost every time a sample falls out of the convex hull $\mathcal{\bar{W}}$. Safety, however, is typically quickly recovered, as the RL parameter is instantaneously adjusted to account for the new sample. As also highlighted in Section II and Remark 8 this temporary loss of safety is a fundamental issue for any sample-based approach, unless additional assumptions are introduced. While one could envision the derivation of some measure of reliability of the identified set to be used to provide stronger guarantees, such investigation is beyond the scope of this paper and will be the subject of future research.

V-B2 Approximation Quality

It has been proven in [11, Theorem 1] that, provided that the MPC parametrization is rich enough, the correct value and action-value functions can be recovered exactly. In the context of this paper, however, the parametrization of the noise set is approximate by construction. Consider partitioning the parameter vector as $\boldsymbol{\mathrm{\theta}}=(\boldsymbol{\mathrm{\theta}}_{\mathrm{c}},\boldsymbol{\mathrm{\theta}}_{\boldsymbol{\mathrm{W}}})$, where $\boldsymbol{\mathrm{\theta}}_{\mathrm{c}}$ is the parameter vector directly defining the cost and constraint functions, similarly to the parametrization of [11], while $\boldsymbol{\mathrm{\theta}}_{\boldsymbol{\mathrm{W}}}$ is the parameter vector parametrizing function $\boldsymbol{\mathrm{g}}_{\boldsymbol{\mathrm{\theta}}}:=\boldsymbol{\mathrm{g}}_{\boldsymbol{\mathrm{\theta}}_{\boldsymbol{\mathrm{W}}}}$. Then, provided that the parametrization of the cost and constraint functions through $\boldsymbol{\mathrm{\theta}}_{\mathrm{c}}$ is rich enough, the value, action-value functions and policy can be recovered exactly on the feasible domain of the robust MPC (20). Therefore, as opposed to the result of [11], the equivalence only holds on the subset of the feasible domain of the RL problem which can be described by the chosen parametrization of $\boldsymbol{\mathrm{g}}_{\boldsymbol{\mathrm{\theta}}}$.

Since the hypothesis of a perfect parametrization is unrealistic in most relevant applications, $Q$ learning and SARSA will in general not deliver the best performance that can be attained with the selected parametrization, since these algorithms aim at fitting the action-value function rather than directly optimizing performance. It is therefore appealing to resort to policy gradient approaches, which seek the direct minimization of $J$ by manipluating $\boldsymbol{\mathrm{\theta}}$. While in principle the proposed approach is well suited for policy gradient methods (both stochastic and deterministic), the main difficulty that needs extra care to be handled is related to the exploration strategy to be deployed, which in general introduces a bias in the gradients, such that convergence to a local optimum is hindered. This topic is investigated in [23].

V-B3 Model Adaptation

In principle, one could choose to let RL update any of the parameters (22) of the MPC scheme (20), including the model parameters $A$, $B$, $\boldsymbol{\mathrm{b}}$. The model used by MPC plays an important role in (a) the definition of the value and action-value function and (b) the conservatism of the uncertainty set approximation $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$. As observed in Remark 7, letting RL adapt the model parameters results in a large increase of computations. The feedback matrix $K$, however, can be adapted by RL without major issues. Since it has an impact on the size of the uncertainty propagation and, through $\boldsymbol{\mathrm{c}}_{k}$ and $\boldsymbol{\mathrm{g}}$, on the cost, letting RL adapt $K$ is expected to further reduce the closed-loop cost by reducing the tightening of the specific constraint components $\boldsymbol{\mathrm{c}}_{k,i}$, $\boldsymbol{\mathrm{g}}_{j}$ corresponding to constraints, which are active in the execution of the the control task. Note that selecting a $K$ that is optimal for the given task is an open issue which can be addressed in a data-driven fashion using the approach we propose.

In [11, Theorem 1] and [20, Corollary 2] it is proven that RL can find $Q_{\boldsymbol{\mathrm{\theta}}}=Q$, $V_{\boldsymbol{\mathrm{\theta}}}=V$ and $\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}=\boldsymbol{\mathrm{\pi}}$ even without adjusting the model, provided that the parametrization is descriptive enough. However, typically the parametrization is low-dimensional and cannot be expected to fulfill this requirement. Therefore, the question on whether adapting the model provides an advantage and on how to best adapt it is still open.

One possibility to improve the approximation quality could be to carry a fixed model parametrization $A_{0}$, $B_{0}$, $\boldsymbol{\mathrm{b}}_{0}$ to be used for safety and one or several (possibly nonlinear) models $\boldsymbol{\mathrm{f}}_{i,\boldsymbol{\mathrm{\theta}}}(\boldsymbol{\mathrm{x}}_{i},\boldsymbol{\mathrm{u}}_{i})$, each associated with its cost to construct a more accurate prediction of the future cost distribution in a scenario-tree fashion. The investigation of alternative formulations will be the subject of future research.

V-B4 Stability

The results on robust MPC are stronger than the one provided in Proposition 1: asymptotic stability is proven, under the assumptions of having $\gamma=1$ and a positive-definite stage cost, usually called tracking cost, as opposed to economic cost. In case of a tracking cost, the presence of a discount factor in the MPC formulation complicates the stability analysis. While one could foresee formulating (20) with $\gamma=1$, even though $\gamma<1$ in the RL problem, it is not immediately clear to which extent this inconsistency will impact on the ability to learn the correct policy. While providing approximation guarantees for $Q$-learning is arguably nontrivial, policy gradient methods will at least provide the best policy for the selected parametrization. In this context that would be be best policy with stability guarantees. By exploiting the results of [35] it should be possible to prove that stabilizing linear control laws can be recovered exactly for linear systems. However, a thorough investigation of this topic is the subject of ongoing research.

The distinction between tracking and economic cost is relevant in relation to the stability properties of the closed-loop system: a thorough discussion on economic MPC can be found in [36, 37, 38]. The linear-quadratic case, which is particularly relevant for our framework, has been analyzed in, e.g., [35, 39, 40, 41, 42]. We only recall here that all the conclusions drawn in [11] apply directly to our setup, i.e., the economic case can be reduced to the tracking case by additionally learning an initial cost.

VI Simulation Results

We test our approach with two examples in simulations. First we consider a linear system to easily interpret the results. Then we consider a realistic example from the chemical industry: a nonlinear evaporation process.

VI-A Linear System

We consider a linear system having $2$ states, such that we can easily visualize the behavior of RL-MPC. Consider a simple linear system with dynamics and stage cost

[TABLE]

where $\boldsymbol{\mathrm{s}}=(p,v)$. We formulate a problem with prediction horizon $N=20$ and introduce the state and control constraints $-\boldsymbol{\mathrm{1}}\leq\boldsymbol{\mathrm{s}}\leq\boldsymbol{\mathrm{1}}$, $-10\leq a\leq 10$. The real noise set is selected as a regular octagon. In order to illustrate the ability of RL to adapt the approximation of the noise set, we select to parametrize $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ as a polytope with 4 facets. We compute the terminal cost matrix $P$ and feedback gain $K$ using the LQR formulation resulting from the nominal model and stage cost. The terminal feedback $K$ is used for constraint tightening, as per Section III.

We simulate the RL-MPC scheme over $200$ time steps in a scenario without exploration in which the reference is

[TABLE]

Since the setpoint reference is moving, for simplicity we impose the terminal constraint centered around the origin. While this does not affect recursive feasibility, practical stability is harder to prove in this case. A thorough analysis of this aspect goes beyond the scope of this paper, and we simply recall that such a terminal constraint induces a leaving arc in the MPC optimal control problem. This situation has been analyzed in [41, 40, 43] in the context of economic MPC, concluding that under suitable assumptions verified by the system considered here, practical stability is obtained.

We update $\boldsymbol{\mathrm{\theta}}=(M,\boldsymbol{\mathrm{m}})$ according to (49)-(50) with $\alpha=0.1$, $\psi$ given by (6) ($Q$ learning). Problem (50) is solved to full convergence at each step such that: (a) globalization techniques guarantee that $\boldsymbol{\mathrm{\theta}}^{*}$ is a better fit than $\boldsymbol{\mathrm{\theta}}_{k}$ for the sample at hand; (b) positive-definiteness of the cost yields a well-posed MPC formulation; and (c) the choice of parameter $\alpha$ is directly related to the horizon of a moving-average approximation of the expected value.

The simulation results are displayed in Figures 2 and 3 in the form of snapshots comparing a representation of the solver at two different time instants. In particular, the constraint, RPI and terminal sets are represented together with the trajectories and constraint tightening. Associated with these quantities, we also display the uncertainty set with the drawn samples, their convex hull and the approximation that RL learned. Initially, RL adapts the rather conservative approximation, therefore enlarging the terminal set by reducing the required constraint tightening. When the setpoint moves, there is initially no gain in modifying the set adaptation, until at $k=34$ the tightened constraint $p\leq 1$ becomes active. At this point, RL starts adapting the set approximation to better capture the shape of the uncertainty set in its top-right part, which is opposite to the bottom-left corner which was better approximated before. Consequently, the terminal set is shifted towards the reference and the $p\leq 1$ is tightened less, with an infinite-horizon difference of approximately $0.019$. Since the tightening steady-state is quickly reached, this difference is visible in Figure 3.

Note that $\boldsymbol{\mathrm{m}}$ does in principle not need to be adjusted, as any adjustment in $\boldsymbol{\mathrm{m}}$ can be equivalently obtained by suitably rescaling $M$. We performed the same simulations with $\theta=M$ and obtained qualitatively equivalent results, the detailed presentation of which we therefore omit. Finally, note that the convex hull of $\mathcal{W}$ is a polytope with $28$ facets, as opposed to the used approximation $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$, which only has $4$ facets.

We ran an additional simulation in which we let $\boldsymbol{\mathrm{\theta}}=(M,K)$, such that RL also adapts the feedback matrix $K$ used for constraint tightening. We remark that the problem of designing the terminal feedback and corresponding set $\mathcal{X}_{\infty}$ is nontrivial and many approaches have been developed. However, to the best of the authors’ knowledge, none of these approaches explicitly accounts for the specific control task to be executed and rather aim at minimizing constraint tightening or maximizing the size of the RPI set. The simulation results are similar to the previous case. However, RL acts on $K$ to enlarge the RPI and terminal set while reducing the required constraint tightening, therefore obtaining an increase in closed-loop performance.

VI-B Evaporation Process

Consider the evaporation process modeled in [44, 45] and used in [46, 39] to demonstrate the potential of economic MPC in the nominal case. The model has states $\boldsymbol{\mathrm{x}}=(X_{2},\,P_{2})$ (concentration and pressure); controls $\boldsymbol{\mathrm{u}}=(P_{100},\,F_{200})$ (pressure and flow); and dynamics given by [46]

[TABLE]

where $T_{2}=aP_{2}+bX_{2}+c$, $T_{3}=dP_{2}+e$, $\lambda F_{4}=Q_{100}-F_{1}C_{\mathrm{p}}(T_{2}-T_{1})$, $T_{100}=fP_{100}+g$, $Q_{100}=UA_{1}(T_{100}-T_{2})$, $UA_{1}=h(F_{1}+F_{3})$, $Q_{100}=UA_{1}(T_{100}-T_{2})$, $UA_{1}=h(F_{1}+F_{3})$, $Q_{200}=\frac{UA_{2}(T_{3}-T_{200})}{1+UA_{2}/(2C_{\mathrm{p}}F_{200})}$, $F_{100}=\frac{Q_{100}}{\lambda_{\mathrm{s}}}$, $\lambda F_{5}=Q_{200}$, $F_{2}=F_{1}-F_{4}$. The model parameters are given in [39].

Concentration $X_{1}$, flow $F_{2}$, and temperatures $T_{1},T_{200}$ are stochastic. In this example, we model them as a uniform distribution with $X_{1}=5\pm 0.5$, $F_{1}=10\pm 0.5$, $T_{1}=40\pm 4$, $T_{200}=25\pm 5$. Additionally, the controller must satisfy bounds $(25,40)\leq(X_{2},P_{2})\leq(100,80)$ on the states and $100\leq(P_{100},F_{200})\leq 400$ on the controls. The stage cost is given by

[TABLE]

which entails that the nominal model is optimally operated at the steady state $\boldsymbol{\mathrm{x}}_{\mathrm{s}}=(25,49.74)$, $\boldsymbol{\mathrm{u}}_{\mathrm{s}}=(191.71,215.89)$, with stage cost $\ell(\boldsymbol{\mathrm{x}}_{\mathrm{s}},\boldsymbol{\mathrm{u}}_{\mathrm{s}})=\ell_{\mathrm{s}}$.

The system is linearized at $\boldsymbol{\mathrm{x}}_{\mathrm{s}},\boldsymbol{\mathrm{u}}_{\mathrm{s}}$ to obtain a linear nominal model. The terminal set is centered at $\boldsymbol{\mathrm{x}}_{\mathrm{c}}=(29,53.57)$, $\boldsymbol{\mathrm{u}}_{\mathrm{c}}=(223.76,221.61)$, which is a steady state for the linearized dynamics. Since the safety constraints are already linear, we introduce discount factor $\gamma=0.99$ and formulate a robust linear MPC problem of the form (20). To that end, we use the linearization at $x_{\mathrm{s}}$, $u_{\mathrm{s}}$, the quadratic cost obtained by applying the tuning procedure proposed in [39]. With the given stage cost and linear model, we obtain the LQR feedback $K$ which we use to stabilize the model error in computing the constraint tightening (III-A) and the terminal set (35). We then use RL to adjust parameter $\boldsymbol{\mathrm{\theta}}=(\boldsymbol{\mathrm{h}},\boldsymbol{\mathrm{p}},M,K)$, i.e., the cost gradient, the uncertainty set and the feedback matrix $K$. Since we formulate the robust MPC problem (20) using a positive-definite cost, while the RL cost $\ell$ is economic, we also learn a quadratic initial cost, as briefly discussed in Section V-B and fully justified in [11].

We let the $Q$-learning algorithm run for $10^{4}$ samples with $\alpha=10^{-2}$. We note that the parameters are no longer adjusted towards the end of the learning procedure, indicating convergence of the algorithm. The RPI and terminal set, as well as the uncertainty set approximation $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ are displayed in Figure 4 at the beginning and at the end of the learning process. One can note that $\boldsymbol{\mathrm{W}}_{\boldsymbol{\mathrm{\theta}}}$ becomes larger in order to better approximate the top left part of the uncertainty set. In combination with the adjustment of the feedback matrix $K$, this allows to shift the terminal set towards the reference, as shown in the top plot.

VII Conclusions and Future Work

In this paper we have presented an RL algorithm which is guaranteed to be safe in the sense of strictly satisfying a set of prescribed constraints given the available data. We have discussed both an innovative function approximation based on robust MPC and an efficient management of data which makes it possible to deal with very large amounts of data in real-time. While linear MPC can be successfully applied also to nonlinear systems, as demonstrated with the second example, in case of strong nonlinearities linear MPC might fail at providing satisfactory performance and even safety.

The proposed framework paves the road for several extensions: (a) one can easily foresee the use of robust MPC based on scenario-trees (which is also suitable for nonlinear models); (b) the use of stochastic or deterministic policy gradient is expected to further improve performance, as discussed in Section II; (c) a formulation using the computational geometry approach for robustness and scenario trees for a refined cost approximation could be envisioned. These research directions are the subject of ongoing research.

References and Notes

• [1]

R. S. Sutton and A. G. Barto, Introduction to Reinforcement Learning, 1st ed.   Cambridge, MA, USA: MIT Press, 1998.

• [2]

D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. van den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot, S. Dieleman, D. Grewe, J. Nham, N. Kalchbrenner, I. Sutskever, T. Lillicrap, M. Leach, K. Kavukcuoglu, T. Graepel, and D. Hassabis, “Mastering the game of go with deep neural networks and tree search,” Nature, vol. 529, pp. 484–503, 2016.

• [3]

S. Wang, W. Chaovalitwongse, and R. Babuska, “Machine learning algorithms in bipedal robot control,” Trans. Sys. Man Cyber Part C, vol. 42, no. 5, pp. 728–743, Sep. 2012.

• [4]

P. Abbeel, A. Coates, M. Quigley, and A. Y. Ng, “An application of reinforcement learning to aerobatic helicopter flight,” in In Advances in Neural Information Processing Systems 19.   MIT Press, 2007, p. 2007.

• [5]

R. S. Sutton, D. McAllester, S. Singh, and Y. Mansour, “Policy gradient methods for reinforcement learning with function approximation,” in Proceedings of the 12th International Conference on Neural Information Processing Systems, ser. NIPS’99.   Cambridge, MA, USA: MIT Press, 1999, pp. 1057–1063.

• [6]

D. Silver, G. Lever, N. Heess, T. Degris, D. Wierstra, and M. Riedmiller, “Deterministic policy gradient algorithms,” in Proceedings of the 31st International Conference on Machine Learning, ser. ICML’14, 2014, pp. I–387–I–395.

• [7]

C. Watkins, “Learning from delayed rewards,” Ph.D. dissertation, King’s College, Cambridge, 1989.

• [8]

J. F. J. Garcia, “A comprehensive survey on safe reinforcement learning,” Journal of Machine Learning Research, vol. 16, pp. 1437–1480, 2013.

• [9]

G. Dalal, K. Dvijotham, M. Vecerik, T. Hester, C. Paduraru, and Y. Tassa, “Safe Exploration in Continuous Action Spaces,” 2018. [Online]. Available: https://arxiv.org/abs/1801.08757

• [10]

T. Pham, G. De Magistris, and R. Tachibana, “Optlayer - practical constrained optimization for deep reinforcement learning in the real world,” in 2018 IEEE International Conference on Robotics and Automation (ICRA), May 2018, pp. 6236–6243.

• [11]

S. Gros and M. Zanon, “Data-Driven Economic NMPC using Reinforcement Learning,” IEEE Transactions on Automatic Control, 2018, (in press).

• [12]

S. Gros, M. Zanon, and A. Bemporad, “Safe Reinforcement Learning via Projection on a Safe Set: How to Achieve Optimality?” in 21st IFAC World Congress, 2020, (submitted).

• [13]

T. Koller, F. Berkenkamp, M. Turchetta, and A. Krause, “Learning-based Model Predictive Control for Safe Exploration and Reinforcement Learning,” 2018, published on Arxiv.

• [14]

A. Aswani, H. Gonzalez, S. S. Sastry, and C. Tomlin, “Provably safe and robust learning-based model predictive control,” Automatica, vol. 49, no. 5, pp. 1216 – 1226, 2013.

• [15]

C. J. Ostafew, A. P. Schoellig, and T. D. Barfoot, “Robust Constrained Learning-based NMPC enabling reliable mobile robot path tracking,” The International Journal of Robotics Research, vol. 35, no. 13, pp. 1547–1563, 2016.

• [16]

F. Berkenkamp, M. Turchetta, A. Schoellig, and A. Krause, “Safe Model-based Reinforcement Learning with Stability Guarantees,” in Advances in Neural Information Processing Systems 30, I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, Eds.   Curran Associates, Inc., 2017, pp. 908–918.

• [17]

R. Murray and M. Palladino, “A model for system uncertainty in reinforcement learning,” Systems & Control Letters, vol. 122, pp. 24 – 31, 2018.

• [18]

F. L. Lewis and D. Vrabie, “Reinforcement learning and adaptive dynamic programming for feedback control,” IEEE Circuits and Systems Magazine, vol. 9, no. 3, pp. 32–50, 2009.

• [19]

F. L. Lewis, D. Vrabie, and K. G. Vamvoudakis, “Reinforcement learning and feedback control: Using natural decision methods to design optimal adaptive controllers,” IEEE Control Systems, vol. 32, no. 6, pp. 76–105, 2012.

• [20]

M. Zanon, S. Gros, and A. Bemporad, “Practical Reinforcement Learning of Stabilizing Economic MPC,” in Proceedings of the European Control Conference, 2019, (accepted).

• [21]

B. Amos, I. D. J. Rodriguez, J. Sacks, B. Boots, and J. Z. Kolter, “Differentiable mpc for end-to-end planning and control,” in Proceedings of the 32Nd International Conference on Neural Information Processing Systems, ser. NIPS’18.   USA: Curran Associates Inc., 2018, pp. 8299–8310.

• [22]

S. Dean, S. Tu, N. Matni, and B. Recht, “Safely Learning to Control the Constrained Linear Quadratic Regulator,” in 2019 American Control Conference (ACC), July 2019, pp. 5582–5588.

• [23]

S. Gros and M. Zanon, “Safe Reinforcement Learning Based on Robust MPC and Policy Gradient Methods,” IEEE Transactions on Automatic Control (submitted), 2019.

• [24]

J. B. Rawlings and D. Q. Mayne, Model Predictive Control: Theory and Design.   Nob Hill, 2012, ch. Postface to model predictive control: theory and design.

• [25]

L. Chisci, J. Rossiter, and G. Zappa, “Systems with persistent disturbances: predictive control with restricted constraints,” Automatica, vol. 37, pp. 1019–1028, 2001.

• [26]

D. Q. Mayne, “Model predictive control: Recent developments and future promise,” Automatica, vol. 50, no. 12, pp. 2967 – 2986, 2014.

• [27]

D. Mayne, “Robust and stochastic mpc: Are we going in the right direction?” IFAC-PapersOnLine, vol. 48, no. 23, pp. 1 – 8, 2015, 5th IFAC Conference on Nonlinear Model Predictive Control NMPC 2015.

• [28]

I. Kolmanovsky and E. Gilbert, “Theory and computation of disturbance invariant sets for discrete-time linear systems,” Math. Probl. Eng., vol. 4, no. 4, pp. 317–367, 1998.

• [29]

C. Büskens and H. Maurer, Online Optimization of Large Scale Systems.   Berlin, Heidelberg: Springer Berlin Heidelberg, 2001, ch. Sensitivity Analysis and Real-Time Optimization of Parametric Nonlinear Programming Problems, pp. 3–16.

• [30]

J. Nocedal and S. Wright, Numerical Optimization, 2nd ed., ser. Springer Series in Operations Research and Financial Engineering.   Springer, 2006.

• [31]

P. Scokaert and J. Rawlings, “Feasibility Issues in Linear Model Predictive Control,” AIChE Journal, vol. 45, no. 8, pp. 1649–1659, 1999.

• [32]

R. Fletcher, Practical Methods of Optimization, 2nd ed.   Chichester: Wiley, 1987.

• [33]

H. Bock, “Recent advances in parameter identification techniques for ODE,” in Numerical Treatment of Inverse Problems in Differential and Integral Equations, P. Deuflhard and E. Hairer, Eds.   Boston: Birkhäuser, 1983, pp. 95–121.

• [34]

D. Bertsekas and I. Rhodes, “Recursive state estimation for a set-membership description of uncertainty,” IEEE Transactions on Automatic Control, vol. 16, pp. 117–128, 1971.

• [35]

M. Zanon, S. Gros, and M. Diehl, “Indefinite Linear MPC and Approximated Economic MPC for Nonlinear Systems,” Journal of Process Control, vol. 24, pp. 1273–1281, 2014.

• [36]

M. Diehl, R. Amrit, and J. Rawlings, “A Lyapunov Function for Economic Optimizing Model Predictive Control,” IEEE Trans. of Automatic Control, vol. 56, no. 3, pp. 703–707, March 2011.

• [37]

R. Amrit, J. Rawlings, and D. Angeli, “Economic optimization using model predictive control with a terminal cost,” Annual Reviews in Control, vol. 35, pp. 178–186, 2011.

• [38]

M. A. Müller, D. Angeli, and F. Allgöwer, “On necessity and robustness of dissipativity in economic model predictive control,” IEEE Transactions on Automatic Control, vol. 60, no. 6, pp. 1671–1676, 2015.

• [39]

M. Zanon, S. Gros, and M. Diehl, “A Tracking MPC Formulation that is Locally Equivalent to Economic MPC,” Journal of Process Control, 2016.

• [40]

M. Zanon and T. Faulwasser, “Economic MPC without terminal constraints: Gradient-correcting end penalties enforce asymptotic stability,” Journal of Process Control, vol. 63, pp. 1 – 14, 2018.

• [41]

T. Faulwasser and M. Zanon, “Asymptotic Stability of Economic NMPC: The Importance of Adjoints,” in Proceedings of the IFAC Nonlinear Model Predictive Control Conference, 2018.

• [42]

L. Grüne and R. Guglielmi, “Turnpike properties and strict dissipativity for discrete time linear quadratic optimal control problems,” SIAM Journal on Control and Optimization, vol. 56, no. 2, pp. 1282–1302, 2018.

• [43]

L. Grüne, “Economic receding horizon control without terminal constraints,” Automatica, vol. 49, pp. 725–734, 2013.

• [44]

F. Y. Wang and I. T. Cameron, “Control studies on a model evaporation process — constrained state driving with conventional and higher relative degree systems,” Journal of Process Control, vol. 4, pp. 59–75, 1994.

• [45]

C. Sonntag, O. Stursberg, and S. Engell, “Dynamic Optimization of an Industrial Evaporator using Graph Search with Embedded Nonlinear Programming,” in Proc. 2nd IFAC Conf. on Analysis and Design of Hybrid Systems (ADHS), 2006, pp. 211–216.

• [46]

R. Amrit, J. B. Rawlings, and L. T. Biegler, “Optimizing process economics online using model predictive control,” Computers & Chemical Engineering, vol. 58, pp. 334 – 343, 2013.