A Framework for Time-Consistent, Risk-Sensitive Model Predictive Control: Theory and Algorithms

TL;DR

This paper introduces a novel risk-sensitive model predictive control framework for linear stochastic systems, ensuring time-consistency and stability, and providing a convex optimization-based algorithm for real-time implementation.

Contribution

It develops a unified, time-consistent risk evaluation framework for MPC, with a provably stabilizing online algorithm and convex reformulation for practical use.

Findings

The proposed MPC algorithm is provably stabilizing.

The control law can be computed via convex optimization.

Simulation results demonstrate effectiveness and stability.

Abstract

In this paper we present a framework for risk-sensitive model predictive control (MPC) of linear systems affected by stochastic multiplicative uncertainty. Our key innovation is to consider a time-consistent, dynamic risk evaluation of the cumulative cost as the objective function to be minimized. This framework is axiomatically justified in terms of time-consistency of risk assessments, is amenable to dynamic optimization, and is unifying in the sense that it captures a full range of risk preferences from risk-neutral (i.e., expectation) to worst case. Within this framework, we propose and analyze an online risk-sensitive MPC algorithm that is provably stabilizing. Furthermore, by exploiting the dual representation of time-consistent, dynamic risk measures, we cast the computation of the MPC control law as a convex optimization problem amenable to real-time implementation. Simulation results are presented and discussed.

Figures

Tables

Table 1. Table 1 : Solve times per MPC iteration for varying lookahead horizons N 𝑁 N for a system with N x = 5 subscript 𝑁 𝑥 5 N_{x}=5 , N u = 2 subscript 𝑁 𝑢 2 N_{u}=2 , L = 6 𝐿 6 L=6 , over 100 simulations. The # Scenarios column corresponds to the number of control nodes in each online MPC scenario tree.

$N$	# Scenarios	Mean (max) [s]
2	7	0.02 (0.035)
3	43	0.15 (0.211)
4	259	1.84 (2.612)
5	1555	54.73 (60.89)

Equations

$$\begin{array}{rl}{\rho }_{k,N}=Z_k& +{\rho }_k(Z_{k+1}+{\rho }_{k+1}(Z_{k+2}+\dots \\ & +{\rho }_{N-2}(Z_{N-1}+{\rho }_{N-1}(Z_N))\dots )),\end{array}$$

$$x_{k+1}=A(w_k)x_k+B(w_k)u_k,$$

$${\mathcal{U}}^{\mathrm{poly}}=\left\{\zeta \in \mathcal{B}\mid S^I\zeta \le T^I,S^E\zeta =T^E\right\}.$$

$$\rho (Z)=\underset{q\in {\mathcal{U}}^{\mathrm{poly}}}{\max }{\mathbb{E}}_q[Z],$$

$$\mathcal{U}^{\mathrm{poly}}=\Bigl{\{}q\in\Delta^{L}\mid S^{I}q\leq T^{I},\,\,S^{E}q=T^{E}\Bigr{\}},$$

$$\text{CVaR}_{\alpha}(Z):=\inf_{y\in}\left[y+\frac{1}{\alpha}\mbox{$\mathbb{E}\left[(Z-y)^{+}\right]$}\right],$$

$$\mathcal{U}^{\mathrm{poly}}=\Bigl{\{}q\in\Delta^{L}\mid 0\leq q(j)\leq\frac{p(j)}{\alpha},j\in\{1,\ldots,L\}\Bigr{\}}.$$

$${\rho }_k(Z(x_{k+1}))=\underset{q\in {\mathcal{U}}_k^{\mathrm{poly}}(x_k,p)}{\max }{\mathbb{E}}_q[Z(x_{k+1})]$$

$${\rho }_{0,k}(0,\dots ,0,x_k^{\top }{x}_k)\le c{\lambda }^k{x}_0^{\top }{x}_0,\text{for all }k\in \mathbb{N},$$

$$\mbox{$\mathbb{E}\left[x_{k}^{\top}x_{k}\right]$}\leq c\,\lambda^{k}\,x_{0}^{\top}x_{0},\quad\emph{for all }k\in{\mathbb{N}},$$

$$x_{k+1}=\{\begin{array}{ll}\sqrt{0.5}{x}_k& \text{with probability }0.2,\\ \sqrt{1.1}{x}_k& \text{with probability }\mathrm{0.8.}\end{array}$$

$$\mbox{$\mathbb{E}\left[x_{k+1}^{\top}Px_{k+1}\right]$}-x_{k}^{\top}Px_{k}\leq-x_{k}^{\top}Lx_{k},$$

$${\text{CVaR}}_{0.5}(Z(x_{k+1}))=\underset{q\in {\mathcal{U}}^{\mathrm{poly}}}{\max }{\mathbb{E}}_q[Z(x_{k+1})],$$

$${\text{CVaR}}_{0.5}(x_{k+1}^2)\ge 0.5x_k^2\frac{0.1}{1.1}+1.1x_k^2\frac{1}{1.1}=1.0455x_k^2.$$

$${\rho }_{0,k}(x_{k+1}^2)={\text{CVaR}}_{0.5}\circ \dots \circ {\text{CVaR}}_{0.5}(x_{k+1}^2)\ge a^{k+1}{x}_0^{\top }{x}_0,$$

$$J_{0,k}(x_{0},\pi):=\rho_{0,k}\Bigl{(}C(x_{0},\pi(x_{0})),\ldots,C(x_{k},\pi(x_{k}))\Bigr{)},$$

$$\begin{array}{rl}\mathbb{X}& :=\{x{\in }^{N_x}:\Vert T_{x}x{\Vert }_2\le x_{\max }\},\\ \mathbb{U}& :=\{u{\in }^{N_u}:\Vert T_{u}u{\Vert }_2\le u_{\max }\}.\end{array}$$

$$\underset{\pi \in \Pi }{\inf }$$

$$x_{k+1}=A(w_k)x_k+B(w_k)\pi (x_k)$$

$$x_k\in \mathbb{X},\pi (x_k)\in \mathbb{U}\forall k$$

$$\text{System is UGRSES}.$$

$$\rho_{0,k}\Bigl{(}C(x_{0},\pi(x_{0})),\ldots,C(x_{k},\pi(x_{k}))\Bigr{)}=C(x_{0},\pi(x_{0}))+\rho(C(x_{1},\pi(x_{1}))+\ldots+\rho(C(x_{k},\pi(x_{k})))\ldots),$$

$$\begin{array}{rl}& b_1\Vert x{\Vert }^2\le V(x)\le b_2\Vert x{\Vert }^2,\text{and}\\ & \rho (V(f(x,w)))-V(x)\le -b_3\Vert x{\Vert }^2.\end{array}$$

$$\begin{array}{rl}& {\rho }_{0,k+1}(0,\dots ,0,b_1\Vert x_{k+1}{\Vert }^2)\\ & \le {\rho }_{0,k+1}(0,\dots ,0,V(x_{k+1}))\\ & ={\rho }_{0,k}(0,\dots ,0,V(x_k)+\rho (V(x_{k+1})-V(x_k)))\\ & \le {\rho }_{0,k}(0,\dots ,0,V(x_k)-b_3\Vert x_k{\Vert }^2)\\ & \le {\rho }_{0,k}(0,\dots ,0,(b_2-b_3)\Vert x_k{\Vert }^2).\end{array}$$

$$\begin{array}{rl}{\rho }_{0,k+1}(0,\dots ,0,V(x_{k+1}))& \le {\rho }_{0,k}(0,\dots ,0,V(x_k)-b_3\Vert x_k{\Vert }^2)\\ & \le \left(1-\frac{b_3}{b_2}\right){\rho }_{0,k}\left(0,\dots ,0,V(x_k)\right).\end{array}$$

$$\begin{array}{rl}& {\rho }_{0,k+1}(0,\dots ,0,V(x_{k+1}))\\ & \le {\left(1-\frac{b_3}{b_2}\right)}^k{\rho }_{0,1}\left(V(x_1)\right)={\left(1-\frac{b_3}{b_2}\right)}^k\rho \left(V(x_1)\right)\\ & \le {\left(1-\frac{b_3}{b_2}\right)}^k\left(V(x_0)-b_3\Vert x_0{\Vert }^2\right)\le b_2{\left(1-\frac{b_3}{b_2}\right)}^{k+1}\Vert x_0{\Vert }^2.\end{array}$$

$${\rho }_{0,k+1}(0,\dots ,0,x_{k+1}^{\top }{x}_{k+1})\le \frac{b_2}{b_1}{\left(1-\frac{b_3}{b_2}\right)}^{k+1}{x}_0^{\top }{x}_0.$$

$${\mathcal{E}}_{\max }(W):=\{x{\in }^{N_x}\mid x^{\top }{W}^{-1}x\le 1\},$$

$$\underset{\begin{array}{c}W=W^{\top }\succ 0\\ P=P^{\top }\succ 0\\ F\end{array}}{\max }$$

$$F^{\top }\frac{T_u^{\top }{T}_u}{u_{\max }^2}F-W^{-1}⪯0$$

$$\begin{array}{rl}& \sum _{j=1}^L{q}_l(j)(A_j+B_{j}F{)}^{\top }P(A_j+B_{j}F)-P\\ & +(F^{\top }RF+Q)\prec 0,\forall q_l\in {\mathcal{U}}^{\mathrm{poly},V}\end{array}$$

$$\forall j\in \{1,\dots ,L\}:$$

$$(A_j+B_{j}F{)}^{\top }\frac{T_x^{\top }{T}_x}{x_{\max }^2}(A_j+B_{j}F)-W^{-1}⪯0$$

$$(A_j+B_{j}F{)}^{\top }{W}^{-1}(A_j+B_{j}F)-W^{-1}⪯0.$$

Full text

1 Introduction

Safety-critical control demands the consideration of uncertain events, and in particular of events with small probabilities that can nevertheless have catastrophic effects if realized. Accordingly, one of the current main research thrusts in Model Predictive Control (MPC) is to find techniques that can robustly address uncertainty  [1, 2, 3]. Techniques for handling uncertainty within the MPC framework fall broadly into three categories: (1) min-max formulations, where the performance indices to be minimized are computed with respect to the worst possible disturbance realization [4, 5, 1], (2) tube-based formulations, where classical (uncertainty-unaware) MPC is modified to use tightened constraints and augmented with a tracking ancillary controller to maintain the system within an invariant tube around the nominal MPC trajectory [6, 7, 8], and (3) stochastic formulations, where risk-neutral expected values of performance indices (and possibly constraints) are considered [9, 10, 11, 12, 13] (see also the recent reviews [1, 3]). The main drawback of the min-max approach is that the control law may be too conservative, since the performance index is being optimized under the worst-case disturbance realizations (which may have an arbitrarily small probability of occurring). The tightened constraints in tube-based formulations induce similar conservatism upon the optimized cost function. On the other hand, stochastic formulations, where the assessment of future random outcomes is accomplished through the expectation operator, may be unsuitable in scenarios where one desires to account for risk, i.e., increased awareness of events of small probability and detrimental consequences.

In general, there are three main challenges with incorporating risk-sensitivity into control and decision-making problems:

Rationality and consistency: The behavior of a control system using a certain risk measure (i.e., a function that maps an uncertain cost to a real number) should be consistent over time. Intuitively, time-consistency stipulates that if a given sequence of costs incurred by the system, when compared to another sequence, has the same current cost and lower risk in the future, then it should be considered less risky at the current time (see Section 2.2 for a formal statement). Examples of “irrational” behavior that can result from a time-inconsistent risk measure include: (1) a control system intentionally seeking to incur losses [14], or (2) deeming states to be dangerous when in fact they are favorable under any realization of the underlying uncertainty [15], or (3) declaring a decision-making problem to be feasible (e.g., satisfying a certain risk threshold) when in fact it is infeasible under any possible subsequent realization of the uncertainties [16]. Remarkably, some of the most common strategies for incorporating risk aversion in decision-making (discussed below) display such inconsistencies [17, 15].

Computational tractability: A risk measure generally adds a nonlinear structure to the optimization problem one must solve in order to compute optimal actions. Hence, it is important to ensure the computational tractability of the optimization problem induced by the choice of a risk measure, particularly in dynamic decision-making settings where the control system must plan and react to disturbances in real-time.

Modeling flexibility: One would like to calibrate the risk measure to the control application at hand by: (1) exploring the full spectrum of risk assessments from worst-case to risk-neutral, and (2) ensuring that the risk measure can be applied to a rich set of uncertainty models (e.g., beyond Gaussian models).

Most popular methods in the literature for assessing risks do not satisfy these three requirements. The Markowitz mean-variance criterion [18], which has dominated risk management for over $50$ years, leads to time-inconsistent assessments of risk in a multi-stage stochastic control setting and also generally leads to computationally intractable problems [14]. Moreover, it is rather limited in terms of modeling flexibility since it relies only on the first two moments of the distribution and there is only a single tuning parameter to trade off between these two moments. Thus, the mean-variance criterion is not well-suited to applications where the disturbance model is non-Gaussian and has been shown to drastically underestimate the effect of extreme events characterized by severe losses [19].

A popular alternative to the mean-variance criterion is the entropic risk measure: $\rho(X)=\log\left(\mathbb{E}[e^{\theta X}]\right)/\theta$, with $\theta\in(0,1)$. The entropic risk measure has been widely studied in the financial mathematics [20, 21] and sequential decision making [22, 23, 24] literatures, and for modeling risk aversion in LQG control problems [25, 26]. While the entropic risk is a more computationally tractable alternative to the mean-variance criterion and can also lead to time-consistent behavior [27], practical applications of the entropic risk measure have proven to be problematic. Notice that the first two terms of the Taylor series expansion of $\rho(X)$ form a weighted sum of mean and variance with regularizer $\theta$, i.e., $\rho(X)\approx\mathbb{E}(X)+(\theta/2)\mathbb{E}(X-\mathbb{E}[X])^{2}$. Consequently, the primary concerns are similar to those associated with the mean-variance measure of risk, with the added complication that the exponential term may induce numerical conditioning problems [28]. The entropic risk measure is a particular example of the general class of methods that model risk aversion by using concave utility functions (convex disutility functions in the cost minimization setting). While the expected (dis)utility framework captures the intuitive notion of diminishing marginal utility, it suffers from the issue that even very little risk aversion over moderate costs leads to unrealistically high degrees of risk aversion over large costs [29, 30] (note that this is a limitation of any concave utility function). Additionally, the linear treatment of stochasticity in expected (dis)utility theory is in general too restrictive to be able to account for more general statistics, particularly in the context of capturing extreme events [31, Chapter 6]. Specifically, within the expected (dis)utility model, such events would be simply averaged over along with the low cost (or “safe”) events. We refer the reader to [31, Chapter 6] where the limitations of purely using utility functions for representing risk are discussed from a foundational perspective.

In order to overcome such challenges, in this paper we incorporate risk sensitivity in MPC by leveraging recent strides in the theory of dynamic risk measures developed by the operations research community [17]. This allows us to propose a framework that satisfies the requirements outlined above with respect to rationality and consistency, computational tractability, and modeling flexibility. Specifically, the key property of dynamic risk measures is that, by reassessing risk at multiple points in time, one can guarantee time-consistency of risk preferences and the agent’s behavior [17]. In particular, it is proven in [17] that time-consistent risk measures can be represented as a composition of one-step coherent risk measures. Coherent risk measures [32, 33] have been thoroughly investigated and widely applied for static decision-making problems in operations research and finance. Coherent risk measures were originally conceived in [32] from an axiomatization of properties that any rational agent’s risk preferences should satisfy (see Section 2 for a formal statement of these axioms). In addition to being axiomatically justified, coherent risk measures capture a wide spectrum of risk assessments from risk neutral to worst-case and thus provide a unifying approach to static risk assessments. Since time-consistent dynamic risks are composed of one-step coherent risks, they inherit the same modeling flexibility.

Statement of Contributions: The contribution of this paper is threefold. First, we introduce a class of dynamic risk measures, referred to as Markov dynamic polytopic risk measures, that capture a full range of risk assessments and enjoy a geometrical structure that is particularly favorable from a computational standpoint. Second, we present and analyze a risk-sensitive MPC algorithm that minimizes in a receding-horizon fashion a Markov dynamic polytopic risk measure, under the assumption that the system’s model is linear and is affected by stochastic multiplicative uncertainty. Finally, by exploring the geometric structure of Markov dynamic polytopic risk measures, we present a convex programming formulation for risk-sensitive MPC that is amenable to a real-time implementation (for moderate horizon lengths). Our framework has three main advantages: (1) it is axiomatically justified, in the sense that risk, by construction, is assessed in a time-consistent fashion; (2) it is amenable to dynamic and convex optimization, primarily due to the compositional form of Markov dynamic polytopic risk measures and their geometry; and (3) it is general, in that it captures a full range of risk assessments from risk-neutral to worst-case. In this respect, our formulation represents a unifying approach for risk-sensitive MPC.

Our approach is inspired by the work in [17], whereby the authors consider a similar risk-sensitive objective function for controlled Markov processes within an infinite-horizon formulation. Differently from [17], in this paper we consider an MPC formulation, define and address notions of persistent feasibility and stability, and provide real-time algorithms for the solution of problems with continuous state- and control-spaces.

A preliminary version of this paper was presented in [34]. In this extended and revised version, we present the following key extensions: (1) the introduction of constraints on state and control variables, (2) a new offline/online MPC formulation for handling these constraints, and (3) additional numerical experimental results including (i) an illustration of the effects of varying levels of risk-sensitivity (thereby also providing a comparison with the standard risk-neutral MPC formulation), and (ii) a scalability study to assess the computational limitations of the proposed approach.

Organization: The rest of the paper is organized as follows. In Section 2 we provide a review of the theory of dynamic risk measures. In Section 3 we discuss the stochastic model we address in this paper. In Section 4 we introduce and discuss the notion of Markov dynamic polytopic risk measures. In Section 5 we state the infinite horizon optimal control problem we wish to address and in Section 6 we derive conditions for risk-sensitive closed-loop stability. In Section 7 we present the MPC adaptation of the infinite horizon problem and present various solution approaches in Section 8. Numerical experiments are presented and discussed in Section 9. Finally, in Section 10 we draw some conclusions and discuss directions for future work.

2 Review of Dynamic Risk Theory

In this section, we briefly review the theory of coherent and dynamic risk measures, on which we will rely extensively in this paper. The material presented in this section summarizes several novel results in risk theory achieved in the past ten years. Our presentation strives to present this material in an intuitive fashion and with a notation tailored to control applications.

2.1 Static, Coherent Measures of Risk

Consider a probability space (\Omega,\mathcal{F},\mbox{\mathbb{P}}), where $\Omega$ is the set of outcomes (sample space), $\mathcal{F}$ is a $\sigma$-algebra over $\Omega$ representing the set of events we are interested in, and $\mathbb{P}$ is a probability measure over $\mathcal{F}$. In this paper we will focus on disturbance models characterized by probability mass functions (pmfs), hence we restrict our attention to finite probability spaces (i.e., $\Omega$ has a finite number of elements or, equivalently, $\mathcal{F}$ is a finitely generated algebra). Denote with $\mathcal{Z}$ the space of random variables $Z:\Omega\mapsto(-\infty,\infty)$ defined over the probability space $(\Omega,\mathcal{F},\mathbb{P})$. In this paper a random variable $Z\in\mathcal{Z}$ is interpreted as a cost, i.e., the smaller the realization of $Z$, the better. For $Z,W$, we denote by $Z\leq W$ the point-wise partial order, i.e., $Z(\omega)\leq W(\omega)$ for all $\omega\in\Omega$.

By a risk measure we understand a function $\rho(Z)$ that maps an uncertain outcome $Z$ into the extended real line $\cup\{+\infty\}\cup\{-\infty\}$. In this paper we restrict our analysis to coherent risk measures, defined as follows:

Definition 2.1** (Coherent Risk Measures).**

A coherent risk measure is a mapping $\rho:\mathcal{Z}\rightarrow$, satisfying the following four axioms: for all $Z,W\in\mathcal{Z}$,

• A1

Monotonicity: $Z\leq W\Rightarrow\rho(Z)\leq\rho(W)$;

• A2

Translation invariance: $\forall a\in$, $\rho(Z+a)=\rho(Z)+a$;

• A3

Positive homogeneity: $\forall\lambda\geq 0$, $\rho(\lambda Z)=\lambda\rho(Z)$;

• A4

Subadditivity: $\rho(Z+W)\leq\rho(Z)+\rho(W)$.

These axioms were originally conceived in [32] and ensure the “rationality” of one-step risk assessments (we refer the reader to [32] for a detailed motivation of these axioms). One of the main properties for coherent risk measures is a universal representation theorem for coherent risk measures, which in the context of finite probability spaces takes the following form:

Theorem 2.2** (Representation Theorem for Finite Probability Spaces [32]).**

Consider the probability space \{\Omega,\mathcal{F},\mbox{\mathbb{P}}\} where $\Omega$ is finite with cardinality $L\in{\mathbb{N}}$, $\mathcal{F}=2^{\Omega}$, and \mbox{\mathbb{P}}=(p(1),\ldots,p(L)), with all probabilities positive. Let $\mathcal{B}$ be the set of probability density functions: \mathcal{B}:=\Bigl{\{}\zeta\in^{L}:\sum_{j=1}^{L}\,p(j)\zeta(j)=1,\zeta\geq 0\Bigr{\}}. The risk measure $\rho:\mathcal{Z}\rightarrow$ is a coherent risk measure if and only if there exists a convex bounded and closed set $\mathcal{U}\subset\mathcal{B}$ such that $\rho(Z)=\max_{\zeta\in\mathcal{U}}\mathbb{E}_{\zeta}[Z]$.

This result states that any coherent risk measure can be written as an expectation with respect to a worst-case density function $\zeta$, chosen adversarially from a suitable set of test density functions (referred to as the risk envelope).

2.2 Dynamic, Time-Consistent Measures of Risk

This section provides a multi-period generalization of the concepts presented in Section 2.1 and follows closely the discussion in [17]. Consider a probability space $(\Omega,\mathcal{F},\mathbb{P})$, a filtration $\mathcal{F}_{0}\subset\mathcal{F}_{1}\subset\mathcal{F}_{2}\cdots\subset\mathcal{F}_{N}\subset\mathcal{F}$, and an adapted sequence of real-valued random variables $Z_{k}$, $k\in\{0,\ldots,N\}$. We assume that $\mathcal{F}_{0}=\{\Omega,\emptyset\}$, i.e., $Z_{0}$ is deterministic. The variables $Z_{k}$ can be interpreted as stage-wise costs. For each $k\in\{0,\ldots,N\}$, denote with $\mathcal{Z}_{k}$ the space of random variables defined over the probability space $(\Omega,\mathcal{F}_{k},\mathbb{P})$; also, let $\mathcal{Z}_{k,N}:=\mathcal{Z}_{k}\times\cdots\times\mathcal{Z}_{N}$. Given sequences $Z=\{Z_{k},\ldots,Z_{N}\}\in\mathcal{Z}_{k,N}$ and $W=\{W_{k},\ldots,W_{N}\}\in\mathcal{Z}_{k,N}$, we interpret $Z\leq W$ component-wise, i.e., $Z_{j}\leq W_{j}$ for all $j\in\{k,\ldots,N\}$.

A dynamic risk measure is a sequence of mappings $\rho_{k,N}:\mathcal{Z}_{k,N}\rightarrow\mathcal{Z}_{k}$, $k\in\{0,\ldots,N\}$, obeying the following monotonicity property: $\rho_{k,N}(Z)\!\leq\!\rho_{k,N}(W)\text{ for all }Z,W\!\in\!\mathcal{Z}_{k,N}$ such that $Z\leq W$. This monotonicity property (analogous to axiom A2 in Definition 2.1) is a natural requirement for any meaningful dynamic risk measure. In this paper, we restrict our discussion to dynamic risk measures that ensure time-consistency [17]. Informally, this property states that if a certain ‘situation’ is considered less risky than another situation in all states of the world at stage $k+1$, then it should also be considered less risky at stage $k$. Arguably, this is a desirable property to enforce when designing controllers for automatic control systems (in contrast to, e.g., analyzing human’s behavior, which may or may not display such a property).

To define the functional form of time-consistent dynamic risk measures, one must first generalize static coherent risk measures as follows:

Definition 2.3** (Coherent One-step Conditional Risk Measures ([17])).**

A coherent one-step conditional risk measure is a mapping $\rho_{k}:\mathcal{Z}_{k+1}\rightarrow\mathcal{Z}_{k}$, $k\in\{0,\ldots,N-1\}$, with the following four properties: for all $W,W^{\prime}\in\mathcal{Z}_{k+1}$

• •

Monotonicity: $W\leq W^{\prime}\Rightarrow\rho_{k}(W)\leq\rho_{k}(W^{\prime})$;

• •

Translation invariance: $\forall Z\in\mathcal{Z}_{k}$, $\rho_{k}(Z+W)=Z+\rho_{k}(W)$;

• •

Positive homogeneity: $\forall\lambda\geq 0$, $\rho_{k}(\lambda W)=\lambda\rho_{k}(W)$;

• •

Subadditivity: $\rho_{k}(W+W^{\prime})\leq\rho_{k}(W)+\rho_{k}(W^{\prime})$.

We now state the main result of this section.

Theorem 2.4** (Dynamic, Time-consistent Risk Measures ([17])).**

Consider, for each $k\in\{0,\ldots,N\}$, the mappings $\rho_{k,N}:\mathcal{Z}_{k,N}\rightarrow\mathcal{Z}_{k}$ defined as

[TABLE]

where the $\rho_{k}$’s are coherent one-step conditional risk measures. Then, the ensemble of such mappings is a dynamic, time-consistent risk measure.

Remarkably, Theorem 1 in [17] shows (under weak assumptions) that the “multi-stage composition” in equation (1) is indeed necessary for time-consistency. Accordingly, in the remainder of this paper, we will focus on the dynamic, time-consistent risk measures characterized in Theorem 2.4.

3 Model Description

Consider the discrete time system:

[TABLE]

where $k\in{\mathbb{N}}$ is the time index, $x_{k}\in^{N_{x}}$ is the state, $u_{k}\in^{N_{u}}$ is the (unconstrained) control input, and $w_{k}\in\mathcal{W}$ is the process disturbance. We assume that the initial condition $x_{0}$ is deterministic and that $\mathcal{W}$ is a finite set of cardinality $L$, i.e., $\mathcal{W}=\{w^{[1]},\ldots,w^{[L]}\}$. Accordingly, denote $A_{j}:=A(w^{[j]})$ and $B_{j}:=B(w^{[j]})$, $j\in\{1,\ldots,L\}$.

For each stage $k$ and state-control pair $(x_{k},u_{k})$, the process disturbance $w_{k}$ is drawn from set $\mathcal{W}$ according to the pmf $p=[p(1),\,p(2),\ldots,\,p(L)]^{\top}$, where $p(j)=\mathbb{P}(w_{k}=w^{[j]})$, $j\in\{1,\ldots,L\}$. Without loss of generality, we assume that $p(j)>0$ for all $j$. Note that the pmf for the process disturbance is time-invariant, and that the process disturbance is independent of the process history and of the state-control pair $(x_{k},u_{k})$. Under these assumptions, the stochastic process $\{x_{k}\}$ is clearly a Markov process.

Remark 3.1**.**

The results presented in this paper can be immediately extended to the case where the process disturbance pmf is time-varying (for example, it is driven by a separate stationary Markov process, as in the popular Markov Jump model [35]) by defining an augmented state such as $(x_{k},w_{k-1})$. We omit this generalization in the interest of brevity and clarity.

4 Markov Polytopic Risk Measures

In this section we refine the notion of dynamic, time-consistent risk measures (as defined in Theorem 2.4) in two ways: (1) we add a polytopic structure to the dual representation of coherent risk measures, and (2) we add a Markovian structure. This will lead to the definition of Markov dynamic polytopic risk measures, which enjoy favorable computational properties and, at the same time, maintain most of the generality of dynamic, time-consistent risk measures.

4.1 Polytopic Risk Measures

According to the discussion in Section 3, the probability space for the process disturbance has a finite number of elements. Thus, one has $\mathbb{E}_{\zeta}[Z]=\sum_{j=1}^{L}\,Z(j)p(j)\zeta(j)$. In our framework (inspired by [36]), we consider coherent risk measures where the risk envelope $\mathcal{U}$ is a polytope, i.e., there exist matrices $S^{I}$, $S^{E}$ and vectors $T^{I}$, $T^{E}$ of appropriate dimensions such that

[TABLE]

We will refer to coherent risk measures representable with a polytopic risk envelope as polytopic risk measures. Consider the bijective map $q(j):=p(j)\zeta(j)$ (recall that, in our model, $p(j)>0$). Then, by applying such a map, one can easily rewrite a polytopic risk measure as

[TABLE]

where $q$ is a pmf belonging to a polytopic subset of the standard simplex, i.e.:

[TABLE]

where \Delta^{L}:=\bigl{\{}q\in^{L}:\sum_{j=1}^{L}\,q(j)=1,q\geq 0\bigr{\}}. Accordingly, one has $E_{q}[Z]=\sum_{j=1}^{L}\,Z(j)q(j)$ (note that, with a slight abuse of notation, we are using the same symbols as before for $\mathcal{U}^{\mathrm{poly}}$, $S^{I}$, and $S^{E}$). We will refer to the set of vertices of $\mathcal{U}^{\mathrm{poly}}$ as $\mathcal{U}^{\mathrm{poly},V}$.

The class of polytopic risk measures is large; common examples include the expected value (the polytope reduces to the singleton pmf $\{p\}$) and the Conditional Value-at-Risk ($\text{CVaR}_{\alpha}$), defined as:

[TABLE]

where $\alpha\in(0,1]$. Its polytopic risk envelope is given by:

[TABLE]

Additional examples include semi-deviation measures, comonotonic risk measures, spectral risk, and optimized certainty equivalent; see [37] for further examples. The key point is that polytopic risk measures cover a full gamut of risk assessments, ranging from risk-neutral to worst case.

4.2 Markov Dynamic Polytopic Risk Metrics

Note that in the definition of dynamic, time-consistent risk measures, since at stage $k$ the value of $\rho_{k}$ is $\mathcal{F}_{k}$-measurable, the evaluation of risk can depend on the whole past, see [17, Section IV]. For example, the level $\alpha$ in the definition of the $\mathrm{CVaR}_{\alpha}$ risk measure can be an $\mathcal{F}_{k}$-measurable random variable (see [17, Example 3]). This generality, which appears of little practical value in many cases, leads to optimization problems that are intractable. This motivates us to add a Markovian structure to dynamic, time-consistent risk measures (similarly as in [17]). In particular, we consider dynamic risk measures $\rho_{k,N}$ as defined in eq. (1) with coherent one-step conditional risk measures $\rho_{k}(\cdot)$ of the form:

[TABLE]

where $\mathcal{U}^{\mathrm{poly}}_{k}(x_{k},p)=\{q\in\Delta^{L}\mid S^{I}_{k}(x_{k},p)q\leq T^{I}_{k}(x_{k},p),S^{E}_{k}(x_{k},p)q=T^{E}_{k}(x_{k},p)\}$ is the polytopic risk envelope, dependent only on the state $x_{k}$. We term such dynamic, time-consistent risk measures as Markov Dynamic Polytopic Risk Measures.

Consistent with the stationarity assumption on the process disturbance pmf $p$, we will further assume that the polytopic risk envelopes $\mathcal{U}^{\mathrm{poly}}_{k}$ are independent of time $k$ and state $x_{k}$, i.e. $\mathcal{U}^{\mathrm{poly}}_{k}(x_{k},p)=\mathcal{U}^{\mathrm{poly}}(p)$, for all $k$. We stress that our results readily generalize to the more general case as per Remark 3.1.

5 Problem Formulation

In light of Sections 3 and 4, we are now in a position to state the risk-sensitive optimization problem we wish to solve in this paper. We start by introducing a notion of stability tailored to our risk-sensitive context.

Definition 5.1** (Uniform Global Risk-Sensitive Exponential Stabilty).**

System (2) is said to be Uniformly Globally Risk-Sensitive Exponentially Stable (UGRSES) if there exist constants $c\geq 0$ and $\lambda\in[0,1)$ such that for all initial conditions $x_{0}\in^{N_{x}}$,

[TABLE]

where $\{\rho_{0,k}\}$ is a Markov dynamic polytopic risk measure. If condition (6) only holds for initial conditions within some bounded neighborhood $\Omega$ of the origin, the system is said to be Uniformly Locally Risk-Sensitive Exponentially Stable (ULRSES) with domain $\Omega$.

Note that, in general, UGRSES is a more restrictive stability condition than mean-square stability, as illustrated by the following example:

Example 5.2** (Mean-Square Stability versus Risk-Sensitive Stability).**

System (2) is said to be Uniformly Globally Mean-Square Exponentially Stable (UGMSES) if there exist constants $c\geq 0$ and $\lambda\in[0,1)$ such that for all initial conditions $x_{0}\in^{N_{x}}$,

[TABLE]

see [38, Definition 1] and [10, Definition 1]. Consider the discrete time system

[TABLE]

A sufficient condition for system (7) to be UGMSES is that there exist positive definite matrices $P=P^{\top}\succ 0$ and $L=L^{\top}\succ 0$ such that

[TABLE]

for all $k\in{\mathbb{N}}$, see [10, Lemma 1]. One can easily check that with $P=100$ and $L=1$ the above inequality is satisfied, and, hence system (7) is UGMSES.

Assuming risk is assessed according to the Markov dynamic polytopic risk measure $\rho_{0,k}=CVaR_{0.5}\circ\ldots\circ CVaR_{0.5}$, we next show that system (7) is not UGRSES. In fact, using the dual representation of CVaR, one can write

[TABLE]

where $\mathcal{U}^{\mathrm{poly}}=\left\{q\in\Delta^{2}\mid 0\leq q_{1}\leq 0.4,\,\,0\leq q_{2}\leq 1.6\right\}$. Consider the pmf $\overline{q}=[0.1/1.1,\,1/1.1]^{\top}$. Since $\overline{q}\in\mathcal{U}^{\mathrm{poly}}$, one has

[TABLE]

By repeating this argument, one can then show that

[TABLE]

where $a=1.0455$. Hence, one cannot find constants $c$ and $\lambda$ that satisfy equation (6). Consequently, system (7) is UGMSES but not UGRSES.

Consider the MDP described in Section 3 and let $\Pi$ be the set of all stationary feedback control policies, i.e., \Pi:=\Bigl{\{}\pi:^{N_{x}}\rightarrow^{N_{u}}\}. Consider the quadratic cost function $C:^{N_{x}}\times^{N_{u}}\rightarrow_{\geq 0}$ defined as $C(x,u):=\|x\|_{Q}^{2}+\|u\|_{R}^{2}$, where $Q=Q^{\top}\succ 0$ and $R=R^{\top}\succ 0$ are given state and control penalties, and $\|x\|_{A}^{2}$ defines the weighted norm, i.e., $x^{T}Ax$. Define the multi-stage cost function:

[TABLE]

and the compact, convex constraint sets:

[TABLE]

The problem we wish to address is as follows.

Optimization Problem $\mathcal{OPT}$ — Given an initial state $x_{0}\in^{N_{x}}$, solve

[TABLE]

We denote the optimal cost function as $J^{\ast}_{0,\infty}(x_{0})$. Note that the risk measure in the definition of UGRSES is assumed to be identical to the risk measure used to evaluate the cost of a policy. Also, by time-invariance of the conditional risk measures, one can write

[TABLE]

where $\rho(\cdot)$ is a given Markov polytopic risk measure that models the “degree” of risk sensitivity. This paper addresses problem $\mathcal{OPT}$ along three main dimensions: (1) find sufficient conditions for risk-sensitive stability (i.e., for UGRSES); (2) design a convex MPC algorithm to efficiently compute a suboptimal state-feedback control policy; and (3) assess algorithm performance via numerical experiments.

6 Risk-Sensitive Stability

In this section we provide a sufficient condition for system (2) to be UGRSES, under the assumptions of Section 5. This condition relies on Lyapunov techniques and extends to the risk-sensitive setting the condition provided in [10] (specifically, Lemma 6.1 reduces to Lemma 1 in [10] when the risk measure is simply an expectation).

Lemma 6.1** (Sufficient Condition for UGRSES).**

Consider a policy $\pi\in\Pi$ and the corresponding closed-loop dynamics for system (2), denoted by $x_{k+1}=f(x_{k},w_{k})$. The closed-loop system is UGRSES if there exists a function $V(x):^{N_{x}}\rightarrow$ and scalars $b_{1},b_{2},b_{3}>0$, such that for all $x\in^{N_{x}}$,

[TABLE]

We refer to the function $V(x)$ as a risk-sensitive Lyapunov function.

Proof.

From the time-consistency, monotonicity, translational invariance, and positive homogeneity of Markov dynamic polytopic risk measures, condition (9) implies

[TABLE]

Also, since $\rho_{0,k+1}$ is monotonic, one has $b_{1}\rho_{0,k+1}(0,\ldots,0,\|x_{k+1}\|^{2})\geq 0$, which implies $b_{2}\geq b_{3}$ and in turn $(1-b_{3}/b_{2})\in[0,1)$. Since $V(x_{k})/b_{2}\leq\|x_{k}\|^{2}$, by using the previous inequalities one can write:

[TABLE]

Repeating this bounding process, one obtains:

[TABLE]

Again, by monotonicity, the above result implies

[TABLE]

By setting $c=b_{2}/b_{1}$ and $\lambda=(1-b_{3}/b_{2})\in[0,1)$, the claim is proven. ∎

The closed-loop system is ULRSES with domain $\Omega$ if (9) holds within the bounded set $\Omega$.

7 Model Predictive Control Problem

This section describes a MPC strategy that approximates the solution to $\mathcal{OPT}$. We note that while an exact solution to $\mathcal{OPT}$ would lead to time-consistent risk assessments, MPC is not guaranteed to be time consistent over an infinite horizon realization, due to its receding horizon nature. In this regard, the MPC strategy provides an efficiently implementable policy that approximately mimics the time-consistent nature of the optimal solution to $\mathcal{OPT}$.

Our receding horizon framework consists of two steps. First, offline, we search for the largest ellipsoidal set $\mathcal{E}_{\max}$ and accompanying local feedback control law $u(x)=Fx$ that renders $\mathcal{E}_{\max}$ control invariant and ensures satisfaction of the state and control constraints. Additionally, within the offline step, we search for a terminal cost matrix $P$ (for the online MPC problem) to ensure that the closed-loop dynamics under the model predictive controller are risk-sensitive exponentially stable. The online MPC optimization then constitutes the second step of our framework.

Consider first, the offline step. We parameterize $\mathcal{E}_{\max}$ as follows:

[TABLE]

where $W$ (and hence $W^{-1}$) is a positive definite matrix. The (offline) optimization problem to compute $W$, $F$, and $P$ is presented below.

Optimization Problem $\mathcal{PE}$ — Solve

[TABLE]

Note that inequality (13) (crucial to guaranteeing closed-loop risk-sensitive stability) is bi-linear in the decision variables. In Section 8, we will derive an equivalent Linear Matrix Inequality (LMI) characterization of (13) in order to derive efficient solution algorithms. We first state the implications of problem $\mathcal{PE}$.

Lemma 7.1** (Properties of $\mathcal{E}_{\max}$).**

Suppose problem $\mathcal{PE}$ is feasible and $x\in\mathbb{X}\cap\mathcal{E}_{\max}(W)$. Let $u(x)=Fx$. Then, the following statements are true:

$\|T_{u}u\|_{2}\leq u_{\max}$, i.e., the control constraint is satisfied. 2. 2.

$\|T_{x}\left(A(w)x+B(w)u\right)\|_{2}\leq x_{\max}$* surely, i.e., the state constraint is satisfied at the next step almost surely.* 3. 3.

$A(w)x+B(w)u\in\mathcal{E}_{\max}(W)$* surely, i.e., the set $\mathcal{E}_{\max}(W)$ is robust control invariant under the control law $u(x)$.*

Thus, $u(x)\in\mathbb{U}$ and $A(w)x+B(w)u\in\mathbb{X}\cap\mathcal{E}_{\max}(W)$ almost surely.

Proof.

See Appendix A. ∎

Lemma 7.1 establishes $\mathbb{X}\cap\mathcal{E}_{\max}(W)$ as a robust control invariant set under the feasible feedback control law $u(x)=Fx$. This result will be crucial to ascertain the persistent feasibility and closed-loop stability of the online optimization algorithm.

Remark 7.2**.**

One can extend problem $\mathcal{PE}$ to the case with time-varying disturbance pmf by leveraging the state-space augmentation proposed in Remark 3.1. Specifically, one would now require the existence of a state-dependent feedback law for the terminal set $\mathcal{E}_{\max}(W)$. That is, instead of the terminal control law $u(x)=Fx$ for $x\in\mathcal{E}_{\max}(W)$, one now must search for a set of gain matrices $\{F_{j}\}_{j=1}^{L}$ where $u=F_{j}x$ is used when the augmented state is $(x,w^{[j]})$. Thus, this would simply induce an additional $L$ copies of the constraints in problem $\mathcal{PE}$.

We are now ready to formalize the MPC problem. Suppose the feasible set of solutions in problem $\mathcal{PE}$ is non-empty and define $W=W^{\ast}$ and $P=P^{\ast}$, where $W^{\ast},P^{\ast}$ are the maximizers for problem $\mathcal{PE}$. Given a prediction horizon $N\geq 1$, define the MPC cost function:

[TABLE]

where $x_{h|k}$ is the state at time $h$ predicted at stage $k$, $\pi_{h|k}:\mathbb{X}\rightarrow\mathbb{U}$ is the control policy to be applied at time $h$ as determined at stage $k$, and $C_{P}(x):=x^{T}Px$ is the terminal cost function. Then, the online MPC problem is formalized as:

Optimization problem $\mathcal{MPC}$ — Given current state $x_{k|k}\in\mathbb{X}$ and prediction horizon $N\geq 1$, solve

[TABLE]

Note that a Markov policy is guaranteed to be optimal for problem $\mathcal{MPC}$ (see [17, Theorem 2]). The optimal cost function for problem $\mathcal{MPC}$ is denoted by $J^{*}_{k}(x_{k|k})$, and a minimizing policy is denoted by $\{\pi^{\ast}_{k|k},\ldots,\pi^{\ast}_{k+N-1|k}\}$. For each state $x_{k}$, we set $x_{k|k}=x_{k}$ and the (time-invariant) model predictive control law is then defined as

[TABLE]

Note that problem $\mathcal{MPC}$ involves an optimization over time-varying closed-loop policies, as opposed to the classical deterministic case where the optimization is over open-loop control inputs. We will show in Section 8 how to solve problem $\mathcal{MPC}$ efficiently. We now address the persistent feasibility and stability properties for problem $\mathcal{MPC}$.

Theorem 7.3** (Persistent Feasibility).**

Define $\mathcal{X}_{N}$ to be the set of initial states for which problem $\mathcal{MPC}$ is feasible. Assume $x_{k|k}\in\mathcal{X}_{N}$ and the control law is given by (21). Then, it follows that $x_{k+1|k+1}\in\mathcal{X}_{N}$ surely.

Proof.

See Appendix A. ∎

Theorem 7.4** (Stochastic Stability with MPC).**

Suppose the initial state $x_{0}$ lies within $\mathcal{X}_{N}$. Then, under the model predictive control law given in (21), the closed-loop system is ULRSES with domain $\mathcal{X}_{N}$.

Proof.

See Appendix B. ∎

8 Solution Algorithms

Prior to solving problem $\mathcal{MPC}$, one would first need to find a matrix $P$ that satisfies (13) – a bilinear semi-definite inequality in $(P,F)$. While checking feasibility of a bilinear semi-definite inequality is NP-hard [39], one can transform this inequality into an LMI by applying the Projection Lemma [40]. The next two results present LMI characterizations of conditions (12)–(15). The proofs are provided in Appendix C.

Theorem 8.1** (LMI Characterization of Stability Constraint).**

Let $\overline{A}:=\begin{bmatrix}{A}_{1}^{\top}&\ldots&A_{L}^{\top}\end{bmatrix}^{\top}$, $\overline{B}:=\begin{bmatrix}B^{\top}_{1}&\ldots&B_{L}^{\top}\end{bmatrix}^{\top}$, and for each $q_{l}\in\mathcal{U}^{\mathrm{poly},V}$, define $\Sigma_{l}:=\mathrm{diag}(q_{l}(1),\ldots,q_{l}(L))\succ 0$. Consider the following set of LMIs with decision variables $Y$, $G$, $\overline{Q}=\overline{Q}^{\top}\succ 0$:

[TABLE]

for all $l\in\{1,\ldots,\mathrm{card}(\mathcal{U}^{\mathrm{poly},V})\}$. The expression in (13) is equivalent to the set of LMIs in (22) by setting $F=YG^{-1}$ and $P=\overline{Q}^{-1}$.

Furthermore, by the application of the Projection Lemma to the expressions in (12), (14) and (15), we obtain the following corollary:

Corollary 8.2**.**

Suppose the following set of LMIs with decision variables $Y$, $G$, and $W=W^{\top}\succ 0$ are satisfied:

[TABLE]

Then, by setting $F=YG^{-1}$, the LMIs above represent sufficient conditions for the LMIs in (12), (14) and (15).

Note that in Corollary 8.2, strict inequalities are imposed only for the sake of analytical simplicity when applying the Projection Lemma. Using similar arguments as in [4], non-strict versions of the above LMIs may also be derived, for example, leveraging some additional technicalities [41].

A solution approach for the receding horizon adaptation of problem $\mathcal{OPT}$ is to first solve the LMIs in Theorem 8.1 and Corollary 8.2. If a solution for $(P,Y,G,W)$ is found, problem $\mathcal{MPC}$ can be solved via dynamic programming (see [17, Theorem 2]) after state and action discretization, see, e.g., [42, 43]. Note that the discretization process might yield a large-scale dynamic programming problem for which the computational complexity scales exponentially with the resolution of the discretization. This motivates the convex programing approach presented next.

8.1 Convex Programming Approach

While problem $\mathcal{MPC}$ is defined as an optimization over Markov control policies, in the convex programming approach, we re-define the problem as an optimization over history-dependent policies. One can show (with a virtually identical proof) that the stability Theorem 7.4 still holds when history-dependent policies are considered. Furthermore, since Markov policies are optimal in our setup, the value of the optimal cost stays the same. The key advantage of history-dependent policies is that their additional flexibility leads to a convex formulation of the online problem. Consider the following parameterization of history-dependent control policies. Let $j_{0},\ldots,j_{h}\in\{1,\ldots,L\}$ be the realized indices for the disturbances in the first $h+1$ steps of the $\mathcal{MPC}$ problem, where $h\in\{1,\ldots,N-1\}$. The control to be exerted at stage $h$ is denoted by ${U}_{h}(j_{0},\ldots,j_{h-1})$. Similarly, we refer to the state at stage $h$ as ${X}_{h}(j_{0},\ldots,j_{h-1})$. The dependence on $(j_{0},\ldots,j_{h-1})$ enables us to keep track of the growth of the scenario tree. In terms of this new notation, the system dynamics (2) can be rewritten as:

[TABLE]

where $X_{0}:=x_{k|k}$, and constraints (19) and (20) can be rewritten as:

[TABLE]

for all $j_{0},\ldots,j_{N-1}\in\{1,\ldots,L\}$. The final solution algorithm, termed convex MPC ($\mathcal{CMPC}$), is presented below.

Algorithm $\mathcal{CMPC}$ — Given an initial state $x_{0}\in\mathbb{X}$ and a prediction horizon $N\geq 1$, solve

[TABLE]

[TABLE]

Denote optimizers: $\{W^{\ast},\overline{Q}^{\ast}\}$.

[TABLE]

Let $(W,P)=(W^{\ast},(\overline{Q}^{\ast})^{-1})$. At each step $k\in\{0,1,\ldots\}$, solve:

[TABLE] 2. 2.

Set $\pi^{MPC}(x_{k|k})=\overline{U}_{0}$.

The cost function in the online problem can be expressed as a nested sequence of convex quadratic inequalities by iteratively applying an epigraph reformulation (see Appendix D for an illustrative example). This results in a convex quadratically-constrained quadratic program (QCQP) which may be solved very efficiently even for moderate (clarified in the next section) values of $N$. In particular, the epigraph reformulation introduces an extra $O(L^{N-1})$ variables and $O(ML^{N-1})$ quadratic inequalities (where $M=\mathrm{card}(\mathcal{U}^{\mathrm{poly},V})$) to the existing $O(N_{u}L^{N-1})$ control variables in the online MPC problem.

As a degenerate case, when we exclude all lookahead steps, problem $\mathcal{MPC}$ is reduced to an offline optimization. By trading off performance, one can compute the control policy offline and implement it directly online without further optimization:

Algorithm $\mathcal{MPC}^{0}$ — Given $x_{0}\in\mathbb{X}$, solve:

[TABLE]

Then, set $\pi^{MPC}(x_{k})=YG^{-1}x_{k}$.

The domain of feasibility for $\mathcal{MPC}^{0}$ is the control invariant set $\mathbb{X}\cap\mathcal{E}_{\max}(W)$. Showing ULRSES for algorithm $\mathcal{MPC}^{0}$ is more straightforward than the corresponding analysis for problem $\mathcal{MPC}$ and is summarized within the following corollary.

Corollary 8.3** (Quadratic Lyapunov Function).**

Suppose problem $\mathcal{MPC}^{0}$ is feasible. Then, system (2) under the offline MPC policy: $\pi^{MPC}(x_{k})=YG^{-1}x_{k}$ is ULRSES with domain $\mathbb{X}\cap\mathcal{E}_{\max}(W)$.

Proof.

From Theorem 8.1, we know that the set of LMIs in (22) is equivalent to the expression in (13) when $F=YG^{-1}$. Then since $x_{0}\in\mathbb{X}\cap\mathcal{E}_{\max}(W)$, a robust control invariant set under the local feedback control law $u(x)=YG^{-1}x$, exploiting the dual representation of Markov polytopic risk measures yields the inequality

[TABLE]

where $L=Q+\left(YG^{-1}\right)^{T}R\left(YG^{-1}\right)=L^{\top}\succ 0$. Define the Lyapunov function $V(x)=x^{\top}Px$. Set $b_{1}=\lambda_{\min}(P)>0$, $b_{2}=\lambda_{\max}(P)>0$ and $b_{3}=\lambda_{\min}(L)>0$. Then by Lemma 6.1, this stochastic system is ULRSES with domain $\mathbb{X}\cap\mathcal{E}_{\max}(W)$. ∎

Note that our algorithms require a vertex representation of the polytopic risk envelopes (rather then the hyperplane representation in eq. (3)). In our implementation, we use the vertex enumeration function included in the MPT toolbox [44], which relies on the simplex method.

9 Numerical Experiments

In this section we present several numerical experiments that were run on a 2.6 GHz Intel Core i7 laptop, using the MATLAB YALMIP Toolbox (version 3.0 [45]) with the Mosek solver [46].

9.1 Effects of Risk Aversion

In the first example, we consider the system studied in [10], a similarly motivated work with an identical dynamical model but restricted to a risk neutral formulation. We remove the conditional dependence in the Markov chain governing $w_{k}$, i.e., all rows of the transition matrix governing the $w_{k}$ Markov chain are set to be the same to be consistent with model (2). As per Remark 3.1, the extension to the general case is straightforward. The goal of the first experiment is to study the effects of using a risk-sensitive objective.

Specifically, consider the second-order system defined by the transition matrices:

[TABLE]

where $\bar{w}_{j}\in\{0.8,1.2,-0.4\}$ with pmf $p=[0.5,0.3,0.2]$. The state and control constraints are defined by $T_{x}=\mathrm{diag}(1/10,1/2)$, $x_{\max}=1$, $T_{u}=1$, $u_{\max}=1$. The cost matrices are $Q=\mathrm{diag}(1,5)$ and $R=1$. We choose the conditional Markov polytopic risk measure: $\mathrm{CVaR}_{\alpha}$. Note that $\alpha=1$ corresponds to the standard risk neutral objective while $\alpha\ll 1$ corresponds to a worst-case risk assessment. For each value of $\alpha$ within the set $\{0.001,0.5,1.0\}$, we ran $1000$ simulations starting at $x_{0}=(6,1)^{T}$ (a point lying in $\mathbb{X}\setminus\mathcal{E}_{\max}$), with 15 online MPC iteration steps and lookahead horizon $N=4$. Each MPC iteration took on average 24.82ms (for the case $\alpha=0.5$ corresponding with the largest vertex set $\mathcal{U}^{\mathrm{poly},V}$). In Figure 1, we plot the empirical cumulative density functions (cdfs) for the cumulative111Recall that by translational invariance, $\rho_{0,k-1}(C_{0},\ldots,C_{k})=\rho_{0}\circ\cdots\circ\rho_{k-1}(C_{0}+\cdots+C_{k})$. cost distribution at various time indexes.

Notice that as $\alpha$ decreases, the optimization further targets the high-cost tail of the cost distribution, at the expense of higher mean cost222One could additionally consider the convex combination $(1-\beta)\mathbb{E}[\cdot]+\beta\,\mathrm{CVaR}_{\alpha}(\cdot)$ for $\beta\in[0,1]$ to characterize the risk-sensitive Pareto trade-off curve.. This is clearly observed in Figure 1 that shows the tail quantile value decreasing as $\alpha$ decreases. Thus, using a single algorithm ($\mathcal{CMPC}$), we are able to generate tunable risk-sensitive policies from risk-neutral to worst-case.

9.2 Computational Limits of $\mathcal{CMPC}$

In this example, we randomly generate a set of $L=6$ systems using MATLAB’s $\mathrm{drss}$ function with $N_{x}=5$ and $N_{u}=2$ to investigate the computational limits of our algorithms. The constraints are defined by $(T_{x},x_{\max})=(2I_{5},5)$, $(T_{u},u_{\max})=(I_{2},1)$, and the cost is defined by the weighting matrices $Q=2I_{5}$ and $R=I_{2}$. The nominal pmf $p$ is randomly generated and experiments were performed with conditional risk measure $\mathrm{CVaR}_{\alpha}$ with $\alpha=0.2$ (yielding a vertex set of size 20), varying lookahead horizons, and 15 MPC iterations for each simulation. The results are summarized in Table 1.

The table illustrates the applicability of the algorithm on a fairly large (with respect to the number of problem variables and constraints) dimensional example, with appreciable lookahead. The exponential growth in computation time is an unavoidable feature of scenario-based optimization. Thus, as typical of jump dynamic systems, we envision the applicability of this work to model temporally-extended (i.e., mode-switching) dynamics as opposed to fast dynamical systems. Improving the runtime capabilities of this algorithm must undoubtedly rely on massively parallel sampling in concert with branch-and-bound techniques, and is left for future research.

10 Conclusion and Future Work

In this paper we presented a framework for risk-sensitive MPC by leveraging recent advances in the theory of dynamic risk measures developed by the operations research community. The proposed approach has the following advantages: (1) it is axiomatically justified and leads to time-consistent risk assessments; (2) it is amenable to dynamic and convex programming; and (3) it is general, in that it captures a full range of risk assessments from risk-neutral to worst case (due to the generality of Markov polytopic risk measures). Our framework thus provides a unifying perspective on risk-sensitive MPC.

This paper opens several directions for future research. First, we plan to extend our work to handle cases where the state and control constraints are required to hold only with a given probability threshold (in contrast to hard constraints) by exploiting techniques such as probabilistic invariance [47]. This relaxation has the potential to provide significantly improved performance at the risk of occasionally violating constraints. Second, we plan to combine our approach with methods for scenario tree optimization in order to reduce the online computation load. Third, while polytopic risk measures encompass a wide range of possible risk assessments, extending our work to non-polytopic risk measures and more general stage-wise costs can broaden the domain of application of the approach. Fourth, we plan to generalize this framework to allow for nonlinear dynamics and more expressive models of uncertainty (e.g., time-varying distributions). Fifth, an important consideration from a practical standpoint is the choice of risk measure appropriate for a given application. We plan to develop principled approaches for making this choice, e.g., by computing polytopic risk envelopes based on confidence regions for the disturbance model.

Appendix A Proof of Lemma 7.1 and Theorem 7.3

Proof of Lemma 7.1.

We first prove the first and second statements and thereby establish $u(x)$ as a feasible control law within the set $\mathcal{E}_{\max}(W)$. Notice that:

[TABLE]

From (10), applying the Schur complement, we know that $\|W^{-\frac{1}{2}}x\|_{2}\leq 1$ for any $x\in\mathcal{E}_{\max}(W)$. Thus, by the Cauchy Schwarz inequality, a sufficient condition for (29) is given by $\|T_{u}FW^{\frac{1}{2}}\|_{2}\leq u_{\max}$, which can be written as

[TABLE]

Re-arranging the inequality above yields the expression given in (12). The state constraint can be proved in an identical fashion by leveraging (10) and (14). It is omitted for brevity.

We now prove the third statement. By definition of a robust control invariant set, we are required to show that for any $x\in\mathcal{E}_{\max}(W)$, that is, for all $x$ satisfying the inequality: $x^{\top}W^{-1}x\leq 1$, application of the control law $u(x)$ yields the following inequality:

[TABLE]

Equivalently, by the S-procedure [48], we are required to show the existence of a $\lambda\geq 0$ such that the following condition holds:

[TABLE]

for all $j\in\{1,\ldots,L\}$. By setting $\lambda=1$, one obtains the largest feasibility set for $W$ and $F$. The expression in (15) corresponds to the (1,1) block in the matrix above. ∎

Proof of Theorem 7.3.

Given $x_{k|k}\in\mathcal{X}_{N}$, problem $\mathcal{MPC}$ may be solved to yield a closed-loop optimal control policy:

[TABLE]

such that $x_{k+N|k}\in\mathbb{X}\cap\mathcal{E}_{\max}(W)$ almost surely. Consider problem $\mathcal{MPC}$ at stage $k+1$ with initial condition $x_{k+1|k+1}$. From Lemma 7.1, we know that

[TABLE]

is a feasible control policy at stage $k+1$. Note that this is simply a concatenation of the optimal tail policy from the previous iteration $\{\pi^{*}_{k+h|k}(x_{k+h|k})\}_{h=1}^{N-1}$, with the state feedback law $Fx_{k+N|k}$ for the final step.

Since a feasible control policy exists at stage $k+1$, $x_{k+1|k+1}=A_{j}x_{k|k}+B_{j}\pi^{*}_{k|k}(x_{k|k})\in\mathcal{X}_{N}$ for any $j\in\{1,\ldots,L\}$, completing the proof. ∎

Appendix B Closed-loop Stability of MPC

Proof of Theorem 7.4.

Let $J^{*}_{k}(x_{k|k})$ denote the optimal value function for problem $\mathcal{MPC}$. We will show that $J^{*}_{k}$ is a risk-sensitive Lyapunov function (Lemma 6.1). Specifically, we first show that $J^{*}_{k}$ satisfies the two inequalities in equation (9). Consider the bottom inequality in equation (9). At time $k$ consider problem $\mathcal{MPC}$ with state $x_{k|k}$. The sequence of optimal control policies is given by $\{\pi^{*}_{k+h|k}\}_{h=0}^{N-1}$. Now, consider the sequence:

[TABLE]

which, as we know from Lemma 7.1, is a feasible solution to problem $\mathcal{MPC}$ at stage $k+1$. Thus, for problem $\mathcal{MPC}$ at stage $k+1$ with initial condition given by $x_{k+1|k+1}=A(w_{k})x_{k|k}+B(w_{k})\pi_{k|k}^{\ast}(x_{k|k})$, denote by $\overline{J}_{k+1}(x_{k+1|k+1})$, the problem $\mathcal{MPC}$ objective corresponding to the control policy sequence $\pi_{k+h|k+1}(x_{k+h|k})$. Note that $x_{k+1|k+1}$ (and therefore $\overline{J}_{k+1}(x_{k+1|k+1})$) is a random variable with $L$ possible realizations, given $x_{k|k}$. Define:

[TABLE]

By exploiting the dual representation of Markov polytopic risk metrics, one can write

[TABLE]

Combining the equation above with equation (13), one readily obtains the inequality

[TABLE]

One can then construc the following chain of inequalities:

[TABLE]

where the first equality follows from the definitions of $Z_{k+N}$ and of dynamic, time-consistent risk measures, the second inequality follows from equation (31) and the monotonicity property of Markov polytopic risk metrics (see also [17, Page 242]), the third equality follows from the definition of $\overline{J}_{k+1}(x_{k+1|k+1})$, and the fourth inequality follows from the definition of $J^{*}_{k+1}$ and the monotonicity of Markov polytopic risk metrics.

Consider now the top inequality in equation (9). One can easily bound $J^{*}_{k}(x_{k|k})$ from below according to:

[TABLE]

where $\lambda_{\min}(Q)>0$ by assumption. The upper bound for $J^{*}_{k}(x_{k|k})$ is derived in two steps. First, define

[TABLE]

Suppose $x_{k|k}\in\mathbb{X}\cap\mathcal{E}_{\max}(W)$. From Lemma 7.1, we know that the control policy $\pi_{k+h|k}(x_{k+h|k})=\{Fx_{k+h|k}\}_{h=0}^{N-1}$ is feasible and consequently, $\mathbb{X}\cap\mathcal{E}_{\max}(W)\subseteq\mathcal{X}_{N}$. Defining $\theta_{f}:=\|Q+F^{\top}RF\|_{2}$, we thus have

[TABLE]

for all $x_{k|k}\in\mathbb{X}\cap\mathcal{E}_{\max}(W)$. Exploiting the translational invariance and monotonicity property of Markov polytopic risk metrics, one obtains the upper bound for all $x_{k|k}\in\mathbb{X}\cap\mathcal{E}_{\max}(W)$:

[TABLE]

In order to derive an upper bound for $J^{*}_{k}(x_{k|k})$ with the above structure for all $x_{k|k}\in\mathcal{X}_{N}$, we draw inspiration from a similar proof in [49, Proposition 2.18]. By leveraging the finite cardinality of the disturbance set $\mathcal{W}$ and the set closure preservation property of the inverse of continuous functions, it is possible to show that $\mathcal{X}_{N}$ is closed. Then, since $\mathcal{X}_{N}$ is necessarily a subset of the bounded set $\mathbb{X}$, it follows that $\mathcal{X}_{N}$ is compact. Thus, there exists some constant $\Gamma>0$ such that $J^{*}_{k}(x_{k|k})\leq\Gamma$ for all $x_{k|k}\in\mathcal{X}_{N}$. That $\Gamma$ is finite follows from the fact that $\{\|x_{k+h|k}\|_{2}\}_{h=0}^{N}$ and $\{\|\pi_{k+h|k}(x_{k+h|k})\|_{2}\}_{h=0}^{N-1}$ are finitely bounded for all $x_{k|k}\in\mathcal{X}_{N}$. Now since $\mathcal{E}_{\max}(W)$ is compact and non-empty, there exists a $d>0$ such that $\mathcal{E}_{d}:=\{x\in^{N_{x}}\mid\|x\|_{2}\leq d\}\subset\mathcal{E}_{\max}(W)$. Let $\hat{\beta}=\max\{\beta\|x\|_{2}^{2}\mid\|x\|_{2}\leq d\}$.

Consider, now, the function: $(\Gamma/\hat{\beta})\beta\|x\|_{2}^{2}$. Then since $\beta\|x\|^{2}_{2}>\hat{\beta}$ for all $x\in\mathcal{X}_{N}\setminus\mathcal{E}_{d}$ and $\Gamma\geq\hat{\beta}$, it follows that

[TABLE]

as desired. Combining the results in equations (33), (34), (36), and given the time-invariance of our problem setup, one concludes that $J^{*}_{k}(x_{k|k})$ is a risk-sensitive Lyapunov function for the closed-loop system (2), in the sense of Lemma 6.1. This concludes the proof. ∎

Appendix C Proof of Theorem 8.1 and Corollary 8.2

We first present the Projection Lemma:

Lemma C.1** (Projection Lemma).**

For matrices $\Omega(X)$, $U(X)$, $V(X)$ of appropriate dimensions, where $X$ is a matrix variable, the following statements are equivalent:

There exists a matrix $W$ such that

[TABLE] 2. 2.

The following inequalities hold:

[TABLE]

where $A^{\perp}$ is the orthogonal complement of $A$.

Proof.

See Chapter 2 in [40]. ∎

We now give the proof for Theorem 8.1 by leveraging the Projection lemma:

Proof.

(Proof of Theorem 8.1) Using simple algebraic factorizations, for all $l\in\{1,\ldots,\text{cardinality}\left(\mathcal{U}^{\text{poly},V}(p)\right)\}$, inequality (13) can be be rewritten as

[TABLE]

By Schur complement, the above expression is equivalent to

[TABLE]

where $\overline{Q}=P^{-1}$. Now since $\overline{Q}=\overline{Q}^{\top}\succ 0$ and $R=R^{\top}\succ 0$, we also have the following identity:

[TABLE]

Next, notice that

[TABLE]

Now, set:

[TABLE]

Then by Lemma C.1, inequalities (37) and (38) are equivalent to the existence of a matrix $G$ that satisfies the following inequality for all $l\in\{1,\ldots,\text{cardinality}\left(\mathcal{U}^{\text{poly},V}(p)\right)\}$:

[TABLE]

Setting $F=YG^{-1}$ and pre-and post-multiplying the above inequality by $\text{diag}(I,R^{\frac{1}{2}},I,I)$ yields the LMI given in (22). Furthermore, from the inequality $-\overline{Q}+G+G^{\top}\succ 0$ where $\overline{Q}\succ 0$, we know that $G+G^{\top}\succ 0$. Thus, by the Lyapunov stability theorem, the linear time-invariant system $\dot{x}=-Gx$ with Lyapunov function $x^{\top}x$ is asymptotically stable (i.e., all eigenvalues of $G$ have positive real part). Therefore, $G$ is an invertible matrix and $F=YG^{-1}$ is well defined. ∎

Proof.

(Proof of Corollary 8.2) We will prove that the third inequality in (23) implies inequality (15). Details of the proofs on the implications of the first two inequalities in (23) follow from identical arguments and will be omitted fin the interest of brevity. Using simple algebraic factorizations, inequality (15) may be rewritten (in strict form) as:

[TABLE]

By Schur complement, the above expression is equivalent to

[TABLE]

Furthermore since $W\succ 0$, we also have the identity

[TABLE]

Now, notice that:

[TABLE]

Then by Lemma C.1, inequalities (40) and (41) are equivalent to the existence of a matrix $G$ such that the following inequality holds for all $j\in\{1,\ldots,L\}$:

[TABLE]

Note that Lemma C.1 provides an equivalence (necessary and sufficient) condition between (42) and (15) if $G$ is allowed to be any arbitrary LMI variable. However, in order to restrict $G$ to be the same variable as in Theorem 8.1, the equivalence relation reduces to sufficiency only. Setting $F=YG^{-1}$ in the above expression gives the claim. ∎

Appendix D Epigraph Reformulation of MPC Cost Function

Suppose $N=2$. We wish to solve

[TABLE]

where $C_{T}$ is the terminal state cost $x_{2}^{T}Px_{2}$. Using the polytopic dual representation of a coherent risk measure and the history-dependent parameterization, the MPC problem can be equivalently written as:

[TABLE]

This epigraph reformulation is justified by the fact that the maximum in (5) must occur at one of the vertices in $\mathcal{U}^{\mathrm{poly},V}$.