# A Stochastic Tube-Based MPC Framework with Hard Input Constraints

Carlo Karam, Matteo Tacchi, Mirko Fiacchini

► To cite this version:

Carlo Karam, Matteo Tacchi, Mirko Fiacchini. A Stochastic Tube-Based MPC Framework with Hard Input Constraints. 2026. hal-05523746

HAL Id: hal-05523746

https://hal.science/hal-05523746v1

Preprint submitted on 23 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A Stochastic Tube-Based MPC Framework with Hard Input Constraints

Carlo Karam, Matteo Tacchi-Bénard, and Mirko Fiacchini

## Abstract

This work presents a stochastic tube-based model predictive control framework that guarantees hard input constraint satisfaction for linear systems subject to unbounded additive disturbances. The approach relies on a structured design of probabilistic reachable sets that explicitly incorporates actuator saturation into the error dynamics and bounds the resulting nonlinearity within a convex embedding. The proposed controller retains the computational efficiency and structural advantages of stochastic tube-based approaches while ensuring state chance constraint satisfaction alongside hard input limits. Recursive feasibility and mean-square stability are established for our scheme, and a numerical example illustrates its effectiveness.

## I. INTRODUCTION

Model Predictive Control (MPC) has established itself as a reliable framework for the optimal control of constrained dynamical systems and is widely adopted in applications where performance must be achieved subject to state and input limitations. In the presence of uncertainty, Robust MPC (RMPC) extends the nominal formulation by guaranteeing recursive feasibility and constraint satisfaction for all admissible disturbance realizations through a worst-case formulation, most notably through the tube-based approach [1]–[3]. Although such formulations provide strong deterministic guarantees, they typically induce conservatism by shrinking the admissible state and input sets, thereby limiting achievable performance [4]. Stochastic MPC (SMPC) has emerged as an alternative formulation that exploits probabilistic descriptions of disturbances [5]. By incorporating distributional information, SMPC often replaces hard constraint satisfaction with chance constraints, allowing for rare violations with user-prescribed probability levels. This relaxation generally enlarges the admissible region and improves closed-loop performance compared to worst-case designs. Moreover, the stochastic framework enables the treatment of disturbances without known bounded support [5], such as Gaussian noise, for which deterministic robust formulations are inherently inapplicable.

Despite the substantial progress achieved in SMPC design, the treatment of input constraints remains an open problem. In contrast to state constraints, which can be meaningfully relaxed in a probabilistic sense, input constraints must be satisfied deterministically at all times as they often correspond to physical actuator saturation limits. Ensuring such guarantees in the presence of unbounded disturbances is nontrivial, as the stochastic feedback component may generate arbitrarily large control actions. Scenario and sample-based SMPC approaches can, in principle, accommodate state chance constraints together with hard input bounds by enforcing constraints over a sufficiently large set of disturbance realizations [6]. However, these methods require access to disturbance samples, may entail significant computational effort, and typically provide guarantees only with respect to the sampled scenarios, without establishing strong closed-loop properties [5]. Analytical approaches instead rely on distributional or moment information of the disturbance. A prominent class of methods is based on affine disturbance-feedback parameterizations of the control law, initially introduced for RMPC in [7] and adapted for chance-constrained systems in [8]. Such parameterizations enable convex reformulations of the stochastic optimal control problem by optimizing jointly over open-loop inputs and feedback gains [5]. To address hard input constraints under unbounded noise, [9], [10] introduced saturation within the affine parametrization, thereby ensuring deterministic input feasibility in the absence of state constraints; this framework was subsequently extended to output-feedback settings in [11]. More recently, [12] generalized the approach to handle hard input bounds together with joint state

Matteo Tacchi-Bénard's work is supported by the French National Research Agency (ANR) under the PIVOINE project grant ANR-25-CE48-1599-01.

The authors are with Univ. Grenoble Alpes, Grenoble INP, CNRS, GIPSA-lab (e-mail: carlo.karam [at] gipsa-lab [dot] fr; matteo.tacchi [at] gipsa-lab [dot] fr; mirko.fiacchini [at] gipsa-lab [dot] fr).
constraints. While effective, these formulations are often computationally demanding and do not scale efficiently with
increasing prediction horizon or system dimension, as the number of decision variables grows accordingly. Stochastic
tube-based approaches offer a structurally different alternative. Inspired by robust tube MPC, they decouple nominal
trajectory optimization from uncertainty propagation by computing probabilistic reachable sets (PRS) offline and
tightening constraints accordingly, resulting in a deterministic and low-dimensional MPC problem online [5]. Despite
being extensively investigated as in [13]–[17], most existing stochastic tube schemes do not guarantee hard input
constraints when disturbances are unbounded, as the ancillary feedback gain is fixed offline and may generate
arbitrarily large inputs. Existing attempts to enforce deterministic actuator limits without assuming known robust
disturbance bounds ( [18], [19]) remain limited. In particular, output-feedback-based methods such as [20], [21]
ensure that the applied input satisfies the constraints deterministically, but compromise recursive feasibility due to
unrealistic input trajectories, only ensuring feasibility with a user-specified probability. More recently, a constraint-
tightening SMPC framework for general nonlinear systems accommodating hard input constraints was proposed in
[22], though its application to linear systems is not fully detailed and may introduce significant conservatism.

Motivated by the above considerations, this work develops a stochastic tube-based SMPC framework for linear
systems that guarantees hard input constraint satisfaction under unbounded disturbances by explicitly incorporating
actuator saturation into the system dynamics. Although this modeling choice introduces a nonlinearity in the error
dynamics, we embed it within an ellipsoidal PRS construction. This enables the retention of linear nominal dynamics
and ancillary feedback while ensuring deterministic satisfaction of input constraints together with state chance
constraint guarantees. The preliminary PRS design underlying the proposed framework can be found in [23]. In
that work, we developed a constructive method for computing ellipsoidal PRS for saturated linear systems subject
to unbounded additive disturbances, based on convex bounds of the saturated error dynamics and an ensuing
contraction analysis. The present work builds upon that foundation and makes the following contributions:

1) Integration into a complete SMPC scheme: we embed the PRS construction of [23] into a stochastic tube-based predictive control framework with hard input constraints, yielding a tractable scheme that we term Saturation-Aware SMPC (SA-SMPC).

2) For the resulting SA-SMPC controller, we establish recursive feasibility and state chance constraint satisfaction under unbounded disturbances with known second moments. In addition, we provide rigorous closed-loop guarantees, including mean-square boundedness, both for the practically implemented nominal-cost formulation and for the ideal nonlinear formulation.

The proposed framework guarantees the standard closed-loop properties typically required for SMPC schemes with hard input constraints, while explicitly accounting for the closed-loop structure of the error under the pre-stabilizing feedback in the offline constraint-tightening design. This enables a less conservative enforcement of hard input constraints than that provided in [22], where the tightening is effectively based on open-loop error predictions. To the best of our knowledge, no such framework exists.

The remainder of this paper is organized as follows. Section II introduces the problem formulation and standing assumptions. Section III presents the convex PRS synthesis for the saturated error dynamics and details the conditions under which the base bound can be refined to reduce conservatism. Section IV develops the full SA-SMPC scheme and establishes recursive feasibility, chance constraint satisfaction, and closed-loop stability properties. Finally, Section V provides numerical results validating the theoretical developments and compares the proposed method with a nominal MPC scheme as well as with the affine disturbance-feedback approach with hard input constraints of [12]. Lengthy technical proofs are deferred to Appendices I and II for completeness and clarity of exposition.

A. Notation

The set of positive integers is denoted by $\mathbb{N}$ and the set of nonnegative integers by $\mathbb{N}_0$. For integers $a, b \in \mathbb{N}$ with $a \le b$, we define $\mathbb{N}_a^b := \{n \in \mathbb{N} : a \le n \le b\}$. For a matrix $A \in \mathbb{R}^{n \times m}$, the notation $A_{(j)} \in \mathbb{R}^n$ denotes its $j$-th column, while $A_j \in \mathbb{R}^m$ denotes its $j$-th row. For symmetric matrix $A \in \mathbb{R}^{n \times m}$, $\Lambda_{\max}(A)$ is the maximum eigenvalue of $A$. For $P > 0$ and $r > 0$, we define the ellipsoid $\mathcal{E}(P, r) := \{e \in \mathbb{R}^n : e^\top P e \le r\}$. The Pontryagin difference of two sets $S, T \subseteq \mathbb{R}^n$ is defined as $S \ominus T := \{s \in S \mid s + t \in S, \forall t \in T\}$. For a matrix $K \in \mathbb{R}^{m \times n}$ and a set $S \subseteq \mathbb{R}^n$, we denote the image of $S$ under the linear map $K$ by $KS := \{Ks \in \mathbb{R}^m \mid s \in S\}$. Finally, the measured state at time $k$ is denoted $x_k$, and the $i$-step-ahead state prediction made at time $k$ is denoted $x_{i|k}$, with $x_{0|k} = x_k$.
## II. PROBLEM FORMULATION

Consider an uncertain linear time-invariant (LTI) system whose discrete-time dynamics obey

$$x_{k+1} = Ax_k + Bu_k + w_k, \qquad (1)$$

where $x_k \in \mathbb{R}^n$, $u_k \in \mathbb{R}^m$, and $w_k \in \mathbb{R}^n$ represent the system state, control input, and unknown process noise, respectively. The matrices $A \in \mathbb{R}^{n \times n}$ and $B \in \mathbb{R}^{n \times m}$ are fully known.

**Assumption 1 (Noise distribution):** The disturbance $w_k$ is independently and identically distributed (i.i.d.) according to an unknown distribution with zero mean and known variance, i.e.

$$\mathbb{E}[w_k] = 0, \quad \mathbb{E}[w_k^\top w_k] = W > 0$$

for all $k \in \mathbb{N}_0$.

Note that these conditions may be met by disturbance distributions with infinite support. The state constraint set $\mathcal{X} \subseteq \mathbb{R}^n$ and the input constraint set $\mathcal{U} \subseteq \mathbb{R}^m$ are defined as

$$\mathcal{X} := \{x \in \mathbb{R}^n \mid Hx \leq h\}, \quad \mathcal{U} := \{u \in \mathbb{R}^m \mid \|u\|_{\infty} \leq 1\},$$

We impose state chance constraints, with user-defined violation probability $\varepsilon \in (0, 1)$, as

$$\Pr\{x_k \in \mathcal{X} \mid x_0\} \geq 1 - \varepsilon, \quad \forall k \in \mathbb{N}, \qquad (2)$$

where the probability is conditioned on the initial state $x_0$, and hard input constraints

$$u_k \in \mathcal{U}, \quad \forall k \in \mathbb{N}. \qquad (3)$$

As usual in the SMPC framework, the aim is to solve an optimization problem at every step $k$, minimizing over a sequence of $i$-steps ahead predicted states $x_{i|k}$ and inputs $u_{i|k}$ which relate to each other by

$$x_{i+1|k} = Ax_{i|k} + Bu_{i|k} + w_{i+k}, \quad x_{0|k} = x_k, \qquad (4)$$

where the realized state $x_k$ is measurable. We introduce the finite-horizon cost

$$\mathcal{J}_N(\mathbf{x}_k, \mathbf{u}_k) := \mathbb{E} \left[ \sum_{i=0}^{N-1} \ell(x_{i|k}, u_{i|k}) + V_f(x_{N|k}) \right] \qquad (5)$$

to be minimized at time-step $k$, given a user-defined stage cost $\ell(\cdot, \cdot)$, terminal cost $V_f(\cdot)$, and some horizon $N \in \mathbb{N}$.
The resulting finite-horizon stochastic optimal control problem can be thus stated as

$$\min_{\mathbf{x}_k, \mathbf{u}_k} \mathcal{J}_N(\mathbf{x}_k, \mathbf{u}_k), \qquad (6a)$$

$$\text{s.t.} \quad x_{i+1|k} = Ax_{i|k} + Bu_{i|k} + w_{i+k}, \qquad (6b)$$

$$x_{0|k} = x_k, \qquad (6c)$$

$$u_{i|k} \in \mathcal{U}, \quad i \in \mathbb{N}_0^N, \qquad (6d)$$

$$\Pr\{x_{i|k} \in \mathcal{X} \mid x_k\} \geq 1 - \varepsilon, \quad i \in \mathbb{N}_0^{N-1}, \qquad (6e)$$

$$\Pr\{x_{N|k} \in \mathcal{X}_f \mid x_k\} \geq 1 - \varepsilon, \qquad (6f)$$

where $\mathbf{x}_k = (x_{i|k})_{i=0}^N$, $\mathbf{u}_k = (u_{i|k})_{i=0}^N$, and $\mathcal{X}_f$ is an appropriately designed terminal constraint set.

The optimization problem above is generally intractable due to the nonconvexity and implicit nature of the chance constraints, which require evaluating multivariate probability distributions over future uncertain states. To obtain a tractable formulation, the problem is approximated through a constraint-tightening approach based on reachable sets, which ensures probabilistic feasibility while enabling online optimization of control inputs within the admissible set $\mathcal{U}$. To this end, we adopt a standard stochastic tube approach inspired by robust MPC schemes [1], splitting the predicted state into nominal and stochastic components

$$x_{i|k} = z_{i|k} + e_{i|k}, \qquad (7)$$
where $z_{i|k}$ denotes the nominal state, and $e_{i|k}$ represents its deviation with respect to the stochastic predicted state. The control input is parameterized as

$$u_{i|k} = v_{i|k} + Ke_{i|k}. \qquad (8)$$

where $v_{i|k}$ is the nominal input sequence optimized online (at every step $k$), and $K$ is a fixed stabilizing feedback gain computed offline.

Because the process noise may be unbounded (see Assumption 1), the error term $e_{i|k}$ may become arbitrarily large. Since $K$ is fixed, this consequently makes it difficult to bound the behavior of the pre-stabilizing control action $Ke_{i|k}$, and in turn ensure hard constraints $u_{i|k} \in \mathcal{U}$ for all $k \in \mathbb{N}$. To address this issue, we incorporate the input constraints (3) directly into the system dynamics by introducing a saturation on the control input. The predicted state now evolves according to,

$$x_{i+1|k} = Ax_{i|k} + B\varphi(u_{i|k}) + w_{i+k}, \qquad (9)$$

where $\varphi: \mathbb{R}^m \to \mathbb{R}^m$ is the element-wise saturation function

$$\varphi_j([u_{i|k}]) = \operatorname{sign}([u_{i|k}]_j) \min \{|[u_{i|k}]_j|, 1\}, \qquad (10)$$

for all $j \in \mathbb{N}_1^m$, with saturation bounds normalized to 1. For compactness, define

$$f(e_{i|k}, v_{i|k}) := Ae_{i|k} + B(\varphi(Ke_{i|k} + v_{i|k}) - v_{i|k}), \qquad (11)$$

yielding the nominal and error dynamics

$$z_{i+1|k} = Az_{i|k} + Bv_{i|k}, \qquad (12a)$$

$$e_{i+1|k} = f(e_{i|k}, v_{i|k}) + w_{i+k}. \qquad (12b)$$

Note that, if $z_{0|k} = x_k$ then $e_{0|k} = 0$, from (7). For our development, we require the following regularity conditions.

**Assumption 2 (Regularity conditions):**

(a) The matrix $A$ is Schur stable.

(b) $\|v_{i|k}\|_{\infty} \le 1$, for all $i, k \in \mathbb{N}_0$.

Condition b) is imposed in the eventual SMPC scheme, and is equivalent to enforcing (3) on the nominal system (given normalized saturation bounds).

The error dynamics in (12b) provide the basis for recasting the chance constraints into deterministic, tightened constraints on the nominal state vector. This process is detailed in the following section.

### III. PROBABILISTIC REACHABLE SET CONSTRUCTION

This section details the methodology used to construct the sets required to make the probabilistic constraints in (6e) and (6f) tractable. We begin by formalizing the core set-theoretic concepts employed in that analysis.

**Definition 1 (Probabilistic Reachable Set):** For a given $k \in \mathbb{N}$, a set $\mathcal{R}_k^\varepsilon$ is a horizon $k$ Probabilistic Reachable Set (PRS) with violation probability $\varepsilon \in [0, 1]$ for a stochastic process $\{e_i\}_{i \in \mathbb{N}} \in \mathbb{R}^n$ if

$$\Pr\{e_k \in \mathcal{R}_k^\varepsilon | e_0 = 0\} \ge 1 - \varepsilon.$$

**Definition 2 (Probabilistic Ultimate Bound, [24]):** A set $\mathcal{R}^\varepsilon$ is a Probabilistic Ultimate Bound (PUB) with violation probability $\varepsilon \in [0, 1]$ for a stochastic process $\{e_k\}_{k \in \mathbb{N}} \in \mathbb{R}^n$ if, for any initial condition $e_0$, there exists a finite time $t(e_0)$ such that

$$\Pr\{e_k \in \mathcal{R}^\varepsilon\} \ge 1 - \varepsilon, \quad \forall k \ge t(e_0).$$

While the preceding definitions are valid for sets of arbitrary geometry, the subsequent development focuses on their ellipsoidal representations. Specifically, we develop a recursive procedure to compute ellipsoidal PRS for the $i$-step ahead predicted error $e_{i|k}$ conditioned on $e_{0|k} = 0$.
A. Convex Bounds of Saturated Functions

Reachability analysis typically relies on set propagation methods that require explicit knowledge of the system's dynamics to compute reachable sets under certain favorable conditions [25], [26]. Such conditions break down for the nonlinear error dynamics given above due to the hard nonlinearity introduced by the saturation function $\varphi(\cdot)$ and the dependence on the nominal input $v$. As a result, standard recursive set computation techniques are hardly applicable. To address this issue, we aim to bound $f(e_{i|k}, v_{i|k})$ within a well-defined convex set (as in [27]), which allows for the computation of outer approximations of the true reachable sets of the error dynamics.

**Definition 3 (Support Function):** For a nonempty convex set $C \subseteq \mathbb{R}^n$, the support function evaluated at $\eta \in \mathbb{R}^n$ is defined as $\delta_C(\eta) = \sup_{x \in C} \eta^\top x$.

The properties inherent to support functions lead to the following set inclusion theorem.

**Theorem 1 ([28]):** For nonempty, convex, and closed sets $C_1$ and $C_2$, the inclusion $C_1 \subseteq C_2$ holds if and only if $\delta_{C_1}(\eta) \le \delta_{C_2}(\eta)$ for every $\eta \in \mathbb{R}^n$.

From this result, we present the subsequent theorem, valid for any $e \in \mathbb{R}^n$.

**Theorem 2:** For any $v$ satisfying $\|v\|_{\infty} \le 1$, the function $f(e, v) = Ae + B(\varphi(Ke + v) - v)$ is contained within the set $F(e)$, where

$$F(e) = \co \left( \left\{ \left( A + \sum_{j \in \mathcal{J}} B_{(j)} K_j \right) e \right\}_{\mathcal{J} \subseteq \mathbb{N}_1^m} \right). \quad (13)$$

*Proof:* Consider first the scalar case $m = 1$ and assume $Ke \ge 0$. Given the constraint $-1 \le v \le 1$, it follows that $-1 \le Ke + v \le 1 + Ke$, and consequently

$$\varphi(Ke + v) = \begin{cases} Ke + v & \text{for } Ke + v \le 1, \\ 1 & \text{for } Ke + v \ge 1. \end{cases}$$

This leads to

$$\varphi(Ke + v) - v = \begin{cases} Ke & \text{for } Ke + v \le 1, \\ 1 - v & \text{for } Ke + v \ge 1. \end{cases}$$

Since $Ke$ is unbounded above and $1 - v \ge 0$, the inequality $0 \le \varphi(Ke + v) - v \le Ke$ holds. A similar argument for $Ke < 0$ yields $Ke \le \varphi(Ke + v) - v \le 0$. Generalizing to arbitrary $m$, for any index subset $\mathcal{J} \subseteq \mathbb{N}_1^m$ and $j \in \mathcal{J}$, we obtain

$$\left\{ \begin{array}{ll} 0 \le \varphi(K_j e + v_j) - v_j \le K_j e & \text{when } K_j e \ge 0, \\ K_j e \le \varphi(K_j e + v_j) - v_j \le 0 & \text{when } K_j e < 0. \end{array} \right.$$

By leveraging the support function, we deduce that for every $\eta \in \mathbb{R}^n$ and every $j \in \mathcal{J}$,

$$\eta^{\top} B_{(j)} (\varphi(K_j e + v_j) - v_j) \in \co \{0, \eta^{\top} B_{(j)} K_j e\} \subseteq \mathbb{R}.$$

Equivalently, the bounds

$$\min \{0, \eta^{\top} B_{(j)} K_j e\} \le \eta^{\top} B_{(j)} (\varphi(K_j e + v_j) - v_j) \le \max \{0, \eta^{\top} B_{(j)} K_j e\},$$

are satisfied. Therefore, for any $\eta \in \mathbb{R}^n$, any $v$ with $\|v\|_{\infty} \le 1$, and any $e \in \mathbb{R}^n$, there exists an index set $\mathcal{J}(e, \eta) \subseteq \mathbb{N}_1^m$ such that

$$\begin{aligned} \eta^{\top} f(e, v) &= \eta^{\top} Ae + \sum_{j \in \mathbb{N}_1^m} \eta^{\top} B_{(j)} (\varphi(K_j e + v_j) - v_j) \\ &\le \eta^{\top} Ae + \sum_{j \in \mathcal{J}(e, \eta)} \eta^{\top} B_{(j)} (\varphi(K_j e + v_j) - v_j). \end{aligned}$$

Application of Theorem 1 confirms that $f(e, v) \in F(e)$. ■

Theorem 2 establishes that for every $e \in \mathbb{R}^n$ and every $\|v\|_{\infty} \le 1$, the value $f(e, v)$ lies within the polytope $F(e)$, whose vertices correspond to the linear combinations $(A + \sum_{j \in \mathcal{J}} B_{(j)} K_j) e$ for all $\mathcal{J} \subseteq \mathbb{N}_1^m$. These vertices correspond to realizations of different saturation scenarios: full, partial, and no saturation [29].
The following corollaries are direct consequences of Theorem 2, derived via convexity arguments. To simplify notation, for every $\mathcal{J} \subseteq \mathbb{N}_1^m$ define $A_K(\mathcal{J}) = (A + \sum_{j \in \mathcal{J}} B_{(j)}K_j)$.

**Corollary 1:** If Assumption 2 is satisfied, then the inequality

$$f(e, v)^{\top} P f(e, v) \leq \max_{\mathcal{J} \subseteq \mathbb{N}_1^m} e^{\top} A_K(\mathcal{J})^{\top} P A_K(\mathcal{J}) e, \quad (14)$$

holds for every $e \in \mathbb{R}^n$.

*Proof:* This assertion follows directly from Theorem 2, by convexity of $\mathbb{R}^n \ni f \mapsto f^{\top} P f$.

<figure><img src="image_2.png" /></figure>

These results indicate that the nonlinear dynamics $f(e, v)$ are bounded by the behavior of the linear systems $A_K(\mathcal{J})e$ for all subsets $\mathcal{J} \subseteq \mathbb{N}_1^m$ and all admissible $v$. Consequently, bounds derived for $A_K(\mathcal{J})e$ are also valid for the saturated system.

## B. Base PRS Synthesis

Utilizing the results from Section III-A, we now construct a PRS for the error state $e_{i|k} \in \mathbb{R}^n$.

**Proposition 1:** Let Assumptions 1 and 2 hold, and consider a symmetric positive definite matrix $P = P^{\top} \succ 0$ that satisfies

$$A_K(\mathcal{J})^{\top} P A_K(\mathcal{J}) \leq \lambda P, \quad \forall \mathcal{J} \subseteq \mathbb{N}_1^m \quad (15)$$

for some contraction rate $\lambda \in [0, 1)$. Then, for all $i, k \in \mathbb{N}$, the bound

$$\mathbb{E} \left[ e_{i+1|k}^{\top} P e_{i+1|k} \right] \leq \lambda \mathbb{E} \left[ e_{i|k}^{\top} P e_{i|k} \right] + \mathrm{Tr}(PW) \quad (16)$$

holds for system (12b) and

$$\mathbb{E} \left[ e_{i|k}^{\top} P e_{i|k} \right] \leq \frac{1 - \lambda^i}{1 - \lambda} \mathrm{Tr}(PW), \quad (17)$$

if $e_{0|k} = 0$.

*Proof:* If (15) is satisfied, Corollary 1 can be applied, yielding

$$f(e_{i|k}, v_{i|k})^{\top} P f(e_{i|k}, v_{i|k}) \leq \lambda e_{i|k}^{\top} P e_{i|k}.$$

for all $e_{i|k} \in \mathbb{R}^n$. This inequality and Assumption 1 lead to the following relationship:

$$\begin{aligned} \mathbb{E}[(f(e_{i|k}, v_{i|k}) + w_{i+k})^{\top} P (f(e_{i|k}, v_{i|k}) + w_{i+k})] &= \mathbb{E} \left[ f(e_{i|k}, v_{i|k})^{\top} P f(e_{i|k}, v_{i|k}) \right] + \mathbb{E} \left[ w_{i+k}^{\top} P w_{i+k} \right] \\ &= \mathbb{E} \left[ f(e_{i|k}, v_{i|k})^{\top} P f(e_{i|k}, v_{i|k}) \right] + \mathrm{Tr}(PW) \leq \lambda \mathbb{E} \left[ e_{i|k}^{\top} P e_{i|k} \right] + \mathrm{Tr}(PW). \end{aligned}$$

Considering the error dynamics (12b), we obtain (16). The remainder of the proof proceeds by induction. For $i = 0$, $\mathbb{E}[e_{0|k}^{\top} P e_{0|k}] = 0$ by assumption. For $i = 1$,

$$\mathbb{E} \left[ e_{1|k}^{\top} P e_{1|k} \right] \leq \lambda \mathbb{E} \left[ e_{0|k}^{\top} P e_{0|k} \right] + \mathrm{Tr}(PW) \leq \mathrm{Tr}(PW).$$

For $i = t$, one gets

$$\mathbb{E} \left[ e_{t|k}^{\top} P e_{t|k} \right] \leq \lambda \mathbb{E} \left[ e_{t-1|k}^{\top} P e_{t-1|k} \right] + \mathrm{Tr}(PW) \leq \sum_{j=0}^{t-1} \lambda^j \mathrm{Tr}(PW).$$

The sequence of upper bounds for $\mathbb{E}[e_{i|k}^{\top} P e_{i|k}]$ forms a geometric series with ratio $\lambda$. Thus, inequality (17) is valid for any contraction rate $0 \leq \lambda < 1$.

<figure><img src="image_3.png" /></figure>

**Remark 1:** In the preceding proof, Assumption 1 was used to nullify the cross-terms. Should the analysis be extended to incorporate biased or correlated disturbances, the terms $\mathbb{E}[w_{i+k}]^{\top} P \mathbb{E}[w_{i+k}]$ and $\mathbb{E}[f(e_{i|k}, v_{i|k})^{\top} P w_{i+k}]$ would need to be introduced, respectively. Accommodating such terms within the present framework is possible by assuming known bounds on the means and covariance matrices, and following the approach detailed in [30].

Combining the previous result with Markov's inequality yields the following lemma.
**Lemma 1:** The sequence of sets defined by

$$ \mathcal{R}_i^{\varepsilon}(\lambda) = \left\{ e \in \mathbb{R}^n \mid e^\top Pe \le \frac{1 - \lambda^i}{\varepsilon(1 - \lambda)} \mathrm{Tr}(PW) \right\}, \qquad (18) $$

forms a sequence of Markov-based ellipsoidal PRS for violation probability $\varepsilon$.

*Proof:* Applying Markov's inequality to the random variable $e_k^\top Pe_{i|k}$ gives

$$ \mathrm{Pr} \left\{ e_{i|k}^\top Pe_{i|k} \ge r \right\} \le \frac{\mathbb{E} \left[ e_{i|k}^\top Pe_{i|k} \right]}{r}. $$

Plugging in the bound in (17) into the right hand side yields for $r = \frac{1 - \lambda^i}{\varepsilon(1 - \lambda)} \mathrm{Tr}(PW)$

$$ \mathrm{Pr} \left\{ e_{i|k}^{\top} Pe_{i|k} \le r \right\} \ge 1 - \frac{1 - \lambda^i}{r(1 - \lambda)} \mathrm{Tr}(PW) = 1 - \varepsilon. \qquad \blacksquare $$

Taking the limit of the Markov-based PRS $\mathcal{R}_i^{\varepsilon}(\lambda)$ as $i \to \infty$ produces the following Markov-based ellipsoidal PUB set.

**Lemma 2:** The set defined by

$$ \mathcal{R}^{\varepsilon}(\lambda) = \left\{ e \in \mathbb{R}^n \mid e^\top Pe \le \frac{1}{\varepsilon(1-\lambda)} \mathrm{Tr}(PW) + \zeta \right\}, \qquad (19) $$

with $\zeta > 0$, is a Markov-based ellipsoidal PUB for violation probability $\varepsilon$.

*Proof:* For any initial condition $e_{0|k} \in \mathbb{R}^n$, the recursion in (16) implies

$$ \lim_{i \to \infty} \mathbb{E} \left[ e_{i|k}^{\top} Pe_{i|k} \right] \le \lim_{i \to \infty} \left( \lambda^i \mathbb{E} \left[ e_{0|k}^{\top} Pe_{0|k} \right] + \frac{1 - \lambda^i}{1 - \lambda} \mathrm{Tr}(PW) \right) = \frac{1}{1 - \lambda} \mathrm{Tr}(PW). $$

Hence, by definition of limits, for every $\varepsilon$, $\zeta > 0$ and every $e_{0|k} \in \mathbb{R}^n$, there exists $t(e_{0|k}) \in \mathbb{N}$ such that

$$ \mathbb{E} \left[ e_{i|k}^{\top} Pe_{i|k} \right] \le \frac{1}{1 - \lambda} \mathrm{Tr}(PW) + \varepsilon\zeta, \quad \forall i \ge t(e_{0|k}). $$

Applying Markov's inequality yields

$$ \mathrm{Pr} \left\{ e_{i|k}^{\top} Pe_{i|k} \le \frac{1}{\varepsilon} \left( \frac{1}{1-\lambda} \mathrm{Tr}(PW) + \varepsilon\zeta \right) \right\} \ge 1 - \varepsilon, \quad \forall i \ge t(e_{0|k}), $$

which proves the result.

Note that for any stabilizing gain $K$, condition (15) also requires that the quadratic function satisfies the Lyapunov condition with the same rate $\lambda$ i.e.

$$ A^\top PA \le \lambda P. \qquad (20) $$

The resulting contraction rate $\lambda$ therefore characterizes the open-loop or fully saturated dynamics. This bound, equivalent to the one in [22] for a linearized system, is inherently conservative, as it ignores the stronger contraction achievable near the origin under non-saturated, closed-loop operation. The following section introduces a refined bound to address this limitation.

### C. Conditional Refinement of Base PRS

The Markov-based PRS and PUB presented in Lemmas 1 and 2 exhibit conservatism, as they rely on a contraction rate $\lambda$ that does not fully characterize the system dynamics. In practice, the system often evolves in a region where the control input $Ke + v$ does not saturate, and the dynamics are linear. This region is defined below.

**Definition 4 (Saturation Region of Linearity [31]):** The region of linearity for system (12b), denoted $\mathcal{R}_L$, is the set of all error states $e$ for which $\varphi(Ke + v) = Ke + v$.

Inside $\mathcal{R}_L$, the predicted error dynamics reduce to $e_{i+1|k} = (A + BK)e_{i|k} + w_{i+k}$. The following corollary, which follows from Theorem 2, Proposition 1, and the preceding definition, is presented.
**Corollary 2:** Let Assumptions 1 and 2 hold and $P = P^\top > 0$ be a matrix satisfying the Lyapunov inequality (15) with contraction rate $\lambda \in [0, 1)$. Suppose the gain $K \in \mathbb{R}^{n \times m}$ and a scalar $\lambda_L \in \mathbb{R}$ satisfy

$$ \lambda_L \le \lambda, \qquad (21a) $$

$$ (A + BK)^\top P(A + BK) \le \lambda_L P, \qquad (21b) $$

Then, the inequality

$$ \mathbb{E} [e_{k+1}^\top P e_{k+1}] \le \lambda_L \mathbb{E} [e_k^\top P e_k] + \mathrm{Tr}(PW), \qquad (22) $$

holds for every $k \in \mathbb{N}_0$ such that $e_k \in \mathcal{R}_L$.

*Proof:* Condition (21b) establishes that the contraction rate $\lambda_L$, applicable within $\mathcal{R}_L$, is less than or equal to $\lambda$. Thus $e_{i|k}^\top P e_{i|k} = \mathbb{E} [e_k^\top P e_k]$ for all $e_{i|k} \in \mathcal{R}_L$ and inequality (16) holds with $\lambda_L$ in spite of $\lambda$, resulting in (22).

<figure><img src="image_4.png" /></figure>

To harmonize the open-loop contraction from (15) with the local closed-loop contraction valid in $\mathcal{R}_L$, we impose the following compatibility condition.

**Assumption 3:** The matrices $K$ and $P$, along with parameters $\lambda$ and $\lambda_L$, are such that conditions (15) and (21) are satisfied.

**Remark 2:** Assumption 3 requires the existence of a common Lyapunov function $P$ and a feedback gain $K$ that simultaneously satisfy the open-loop and closed-loop contraction conditions. In particular, the aim is to find a solution with $\lambda_L \le \lambda$, ensuring the closed-loop dynamics contract at least as fast as the open-loop dynamics. A feasible pair $(P, K)$ can be synthesized by solving a convex optimization problem subject to the constraints (15) and (21). Through the change of variables $X = P^{-1}$ and $Y = KX$ and a Schur complement formulation, we get

$$ \min_{X,Y,\lambda_L,\lambda} \quad \lambda + \lambda_L \qquad (23a) $$

$$ \text{s.t.} \quad \lambda_L \le \lambda \qquad (23b) $$

$$ \left[ \begin{array}{cc} \lambda X & (AX + B_{(j)}Y_j)^\top \\ AX + B_{(j)}Y_j & X \end{array} \right] \ge 0, \quad \forall j \in \mathcal{J}, \quad \forall \mathcal{J} \subseteq \mathbb{N}_1^m, \qquad (23c) $$

$$ \left[ \begin{array}{cc} \lambda_L X & (AX + BY)^\top \\ AX + BY & X \end{array} \right] \ge 0, \qquad (23d) $$

where (23c) with $\mathcal{J} = \{\emptyset\}$ is equivalent to (20). The controller gain is then recovered as $K = YX^{-1}$. Problem (23) can be solved to global optimality using the moment-SoS hierarchy for optimization with polynomial matrix inequality constraints [32], but such a development is deferred to a later work. Alternatively, the parameters $\lambda$ and $\lambda_L$ can be optimized via a 2-level bisection search to find the fastest achievable contraction rates. From a practical perspective, it is conceivable that incorporating additional performance criteria may lead to a better overall scheme; e.g. regularization terms involving $KP$, or objectives that promote smaller ellipsoids (such as minimizing $-\log \det Q$). Since we do not develop a dedicated convex formulation that guarantees optimality, we omit such extensions. Note that problem (23) is always feasible under Assumption 2.

We now have two distinct bounds on the evolution of $\mathbb{E} [e_{i|k}^\top P e_{i|k}]$: one that is globally valid but conservative, and another that is tighter, valid only locally within the region of linearity. It is expected that the true system behavior falls somewhere between these two extremes. To capture this behavior more effectively, we combine the worst-case and best-case bounds probabilistically, with weights reflecting the likelihood of the error state's position in the state space. The following theorem helps formalize this approach within the chosen ellipsoidal set framework.

**Theorem 3 (Ellipsoidal Region of Linearity [31]):** The largest ellipsoid with shape matrix $P$ contained within the region of linearity, $\mathcal{E}(P, r_L) \subseteq \mathcal{R}_L$, has a scaling parameter $r_L$ given by

$$ r_L = \min_{1 \le j \le m} \frac{(1 - v_j)^2}{K_j P^{-1} K_j^\top}. \qquad (24) $$

This ellipsoid will be used to determine whether error states reside in $\mathcal{R}_L$. Specifically, if $e^\top P e \le r_L$, then $\varphi(Ke + v) = Ke + v$.
**Proposition 2:** Let Assumptions 1, 2 and 3 hold. Then, the inequality

$$ \mathbb{E} [e_{i+1|k}^{\top} Pe_{i+1|k}] \le \lambda_L \mathbb{E} [e_{i|k}^{\top} Pe_{i|k}] + \frac{\lambda - \lambda_L}{r_L} \mathbb{E}^2 [e_{i|k}^{\top} Pe_{i|k}] + \mathrm{Tr}(PW), \quad (25) $$

is satisfied for all $i, k \in \mathbb{N}_1^m$ for system (12b) with $e_{0|k} = 0$.

*Proof:* Distinct contraction rates $\lambda_L$ and $\lambda$ apply when $e_{i|k} \in \mathcal{R}_L$ and $e_{i|k} \notin \mathcal{R}_L$, respectively. The state space $\mathbb{R}^n$ is partitioned into $S_1 = \{e \in \mathbb{R}^n : e^\top Pe \le r_L\}$ and $S_2 = \{e \in \mathbb{R}^n : e^\top Pe > r_L\}$ with $r_L$ as in (24), where, by construction, $S_1 \subseteq \mathcal{R}_L$. The law of total expectation gives

$$ \mathbb{E} [e_{i+1|k}^{\top} Pe_{i+1|k}] = \sum_{j=1}^{2} \mathbb{E} [e_{i+1|k}^{\top} Pe_{i+1|k} | e_{i|k} \in S_j] \cdot \mathrm{Pr} \{e_{i|k} \in S_j\}, $$

Substituting with the bounds (16) and (22), we get

$$ \mathbb{E} [e_{i+1|k}^{\top} Pe_{i+1|k}] \le (\lambda_L \mathbb{E} [e_{i+1|k}^{\top} Pe_{i+1|k}] + \mathrm{Tr}(PW)) \cdot \mathrm{Pr} \{e_{i|k} \in S_1\} \\ + (\lambda \mathbb{E} [e_{i+1|k}^{\top} Pe_{i+1|k}] + \mathrm{Tr}(PW)) \cdot \mathrm{Pr} \{e_{i|k} \in S_2\}. $$

Since $\mathrm{Pr} \{e_{i|k} \in S_2\} = 1 - \mathrm{Pr} \{e_{i|k} \in S_1\}$, we have

$$ \mathbb{E} [e_{i+1|k}^{\top} Pe_{i+1|k}] \le \lambda \mathbb{E} [e_{i+1|k}^{\top} Pe_{i+1|k}] - (\lambda - \lambda_L) \mathbb{E} [e_{i+1|k}^{\top} Pe_{i+1|k}] \cdot \mathrm{Pr} \{e_{i|k} \in S_1\} + \mathrm{Tr}(PW). \quad (26) $$

Finally, substituting with the Markov bound on $\mathrm{Pr} \{e_{i|k} \in S_1\}$, yields

$$ \begin{aligned} \mathbb{E} [e_{i+1|k}^{\top} Pe_{i+1|k}] &\le \lambda \mathbb{E} [e_{i|k}^{\top} Pe_{i|k}] + (\lambda_L - \lambda) \mathbb{E} [e_{i|k}^{\top} Pe_{i|k}] \left(1 - \frac{\mathbb{E} [e_{i|k}^{\top} Pe_{i|k}]}{r_L}\right) + \mathrm{Tr}(PW), \\ &\le \lambda_L \mathbb{E} [e_{i|k}^{\top} Pe_{i|k}] + \frac{\lambda - \lambda_L}{r_L} \mathbb{E}^2 [e_{i|k}^{\top} Pe_{i|k}] + \mathrm{Tr}(PW). \end{aligned} $$

Therefore, inequality (25) holds for any pair of contraction rates $(\lambda_L, \lambda)$ satisfying $0 \le \lambda_L \le \lambda < 1$. ■

As $r_L$ increases, the bound approaches the linear case. This is consistent with the intuition that when $\mathcal{R}_L$ expands infinitely, nearly all error terms should fall within it. The relation (25) above further indicates that the true contraction rate of the system varies between its linear and saturated (or open-loop) bounds, depending on the magnitude of $r_L$. The following analysis aims at exploring this transition and derives a bound that captures this effective contraction behavior in a more interpretable form.

**Proposition 3:** Let Assumptions 1, 2 and 3 hold, and consider the bound given by (25). If the condition

$$ \frac{1}{1 - \lambda} \mathrm{Tr}(PW) < r_L, \quad (28) $$

is satisfied, then

$$ \mathbb{E} [e_{i|k}^{\top} Pe_{i|k}] \le \frac{1 - (\bar{\lambda}^*)^i}{1 - \bar{\lambda}} \mathrm{Tr}(PW), \quad (29) $$

holds, where $\bar{\lambda}^*$ is the unique value in $[\lambda_L, \lambda]$ satisfying

$$ \frac{1}{1 - \bar{\lambda}^*} \mathrm{Tr}(PW) = \frac{\bar{\lambda}^* - \lambda_L}{\lambda - \lambda_L} r_L. \quad (30) $$

*Proof:* The proof is given in Appendix I.

The result below follows directly from Proposition 3, see its proof.

**Corollary 3:** Let Assumptions 1, 2 and 3 hold, and consider the bound given by (25). If condition (28) holds with $\bar{\lambda}^* \in [\lambda_L, \lambda]$ satisfying (30), then

$$ \mathbb{E} [e_{i|k}^{\top} Pe_{i|k}] \le \frac{1}{1 - \bar{\lambda}^*} \mathrm{Tr}(PW) \quad \implies \quad \mathbb{E} [e_{i+j|k}^{\top} Pe_{i+j|k}] \le \frac{1}{1 - \bar{\lambda}^*} \mathrm{Tr}(PW), \quad \forall j \in \mathbb{N} \quad (31) $$

for any $i \in \mathbb{N}$.
Proposition 3 and Corollary 3 establish the existence of an effective contraction rate $\bar{\lambda}^*$ provided the set $\mathcal{E}(P, r_L)$ is sufficiently large. Such an effective rate enables the derivation of a deterministic bound on $\mathbb{E}[e_{i|k}^\top Pe_{i|k}]$ tighter than the one provided in (22). It can be shown that as $r_L$ grows infinitely large, $\bar{\lambda}^* \to \lambda_L$. Conversely, when $r_L$ falls below the critical threshold $(1 - \lambda)^{-1}\text{Tr}(PW)$, the existence of such a rate can neither be confirmed nor excluded, but it becomes computationally inaccessible by using the Markov inequality. In such a case, the loose Markov inequality

$$ \mathrm{Pr} \left\{ e_{i|k}^{\top} Pe_{i|k} \le r_L \right\} \ge 1 - \frac{\mathbb{E} \left[ e_{i|k}^{\top} Pe_{i|k} \right]}{r_L}, \quad (32) $$

fails to provide a useful bound on the error state's magnitude. Consequently, if (28) does not hold, one must resort to the more conservative bound in (17). Formally, we define

$$ \hat{\lambda} = \begin{cases} \bar{\lambda}^* & \text{if } \frac{1}{1-\bar{\lambda}} \mathrm{Tr}(PW) < r_L, \\ \lambda & \text{if } \frac{1}{1-\lambda} \mathrm{Tr}(PW) \ge r_L, \end{cases} \quad (33) $$

where $\bar{\lambda}^* < \lambda$ is the solution to (30).

**Lemma 3: The sequence of sets defined by**

$$ \mathcal{R}_i^\varepsilon(\hat{\lambda}) = \left\{ e \in \mathbb{R}^n \mid e^\top Pe \le \frac{1 - \hat{\lambda}^i}{\varepsilon(1 - \hat{\lambda})} \mathrm{Tr}(PW) \right\}, \quad (34) $$

forms a sequence of Markov-based ellipsoidal PRS for violation probability $\varepsilon$, with $\hat{\lambda}$ defined in (33).

*Proof:* This follows from Proposition 3 and application of Markov's inequality to the random variable $e_{i|k}^\top Pe_{i|k}$. ■

Taking the limit of the Markov-based PRS $\mathcal{R}_i^\varepsilon(\hat{\lambda})$ as $i \to \infty$ yields the following Markov-based PUB.

**Lemma 4: The set defined by**

$$ \mathcal{R}^\varepsilon(\hat{\lambda}) = \left\{ e \in \mathbb{R}^n \mid e^\top Pe \le \frac{1}{\varepsilon(1 - \hat{\lambda})} \mathrm{Tr}(PW) \right\}, \quad (35) $$

is a Markov-based ellipsoidal PUB for violation probability $\varepsilon$, with $\hat{\lambda}$ as defined in (33).

These sets are analogous to those in Lemmas 1 and 2, using a bound between (29) and the conservative bound (17) when the effective contraction rate $\bar{\lambda}^*$ is inaccessible. These sets will be used for constraint tightening in the full SMPC scheme detailed in the next section.

## IV. SATURATION-AWARE STOCHASTIC MPC

The different ingredients of the proposed SMPC scheme, which we name Saturation-Aware SMPC (SA-SMPC) will be detailed in this section. For the development below, we will be operating under the following assumption.

**Assumption 4:** The nominal control input $v_{i|k}$ is such that

$$ v_{i|k} \in \mathcal{V} := \{v \in \mathbb{R}^m \mid \|v\|_{\infty} \le v_{ss} < 1\}, $$

for all $i, k$, for which we can compute $r_L = \min_{1 \le j \le m} \frac{(1 - v_{ss})^2}{K_j P^{-1} K_j^\top}$ and corresponding $\bar{\lambda}^*$ satisfying Proposition 3, valid for all $i, k$.

Under this assumption, we can reformulate the state chance constraints in (2) as deterministic under the standard tightening below.

a) **State Constraints:** we impose $z_{i|k} \in \mathcal{Z}_i = \mathcal{X} \ominus \mathcal{R}_i^\varepsilon(\hat{\lambda})$ for all $i \in \mathbb{N}_0^{N-1}$.

b) **Terminal Constraint:** we impose $z_{N|k} \in \mathcal{Z}_f \subseteq (\mathcal{X} \ominus \mathcal{R}^\varepsilon(\hat{\lambda})) \cap \mathcal{X}_K$, where $\mathcal{Z}_f$ is a forward invariant set for the autonomous system $z^+ = (A + BK_f)z$ where $K_f$ defines the terminal control law, and $\mathcal{X}_K$ is defined as $\mathcal{X}_K = \{z \in \mathcal{X} \mid K_f z \in \mathcal{V}\}$.
## A. Cost Function and State Initialization

Recall that the finite-horizon stochastic optimal control problem introduced in (6) seeks to minimize the expected value of a cumulative stage cost and terminal cost evaluated along the closed-loop trajectories. For practical implementation, we consider quadratic stage and terminal costs typical in MPC formulations

$$ \ell(x, u) := \|x\|_Q^2 + \|u\|_R^2, \quad V_f(x) := \|x\|_S^2, \quad (36) $$

with $Q \succeq 0$, $R \succ 0$, and $S \succ 0$ satisfying

$$ (A + BK_f)^\top S (A + BK_f) - S \preceq -Q - K_f^\top R K_f, \quad (37) $$

Naturally, a possible choice for the pair $(S, K_f)$ is the solution of the LQR problem for which (37) is satisfied with equality.

A theoretically appealing SMPC formulation is the so-called indirect-feedback framework, in which the real predicted state and input trajectories enter the cost function, while the optimization is performed over constrained nominal variables [14]. In this case, the objective to minimize would be

$$ \mathcal{J}_N = \mathbb{E} \left[ \sum_{i=0}^{N-1} \ell(x_{i|k}, \varphi(u_{i|k})) \right] + \mathbb{E} \left[ V_f(x_{N|k}) \right]. \quad (38) $$

For this formulation, assuming that a minimizer exists at each time step and that appropriate design constraints are imposed on the terminal weight $S$, it is possible to establish closed-loop stability and convergence properties using standard MPC arguments based on cost decrease along shifted trajectories. The properties are formalized and proven in Appendix II.

However, minimizing the indirect-feedback cost in the present saturation-aware framework leads to a generally nonconvex optimization problem. Indeed, the saturation nonlinearity has already been embedded into the predictive model and the probabilistic reachable set construction; reintroducing it inside the objective function would undermine the tractability gains achieved through this reformulation and render the resulting problem unsuitable for online optimization.

As such, we adopt a tractable surrogate objective defined purely in terms of the nominal state and input trajectories. Specifically, we minimize a quadratic cost of the form

$$ \mathcal{J}_k^n = \sum_{i=0}^{N-1} \ell(z_{i|k}, v_{i|k}) + V_f(z_{N|k}) + \rho_k \xi_k^2. \quad (39) $$

where $\rho_k$ is a variable positive scalar weight and $\xi_k \in [0, 1]$ for all $k$ is an interpolating variable. This cost introduces feedback with respect to the measured state $x_k$ through the initialization strategy

$$ z_{0|k} = (1 - \xi_k) x_k + \xi_k z_{1|k-1}. \quad (40) $$

The variable $\xi_k$ allows a continuous interpolation between measured-state initialization and a shifted nominal trajectory, which is instrumental in preserving feasibility in the presence of unbounded disturbances. This initialization approach was originally suggested in [33], [34], [35]. Note that other feasibility-preserving approaches may be adopted, such as that presented in [17].

**Proposition 4 (Closed-loop PRS):** Let Assumptions 2 and 4 hold. Given the initialization rule (40), we have

$$ \Pr\{e_{i|k} \in \mathcal{R}_{i+k}^{\varepsilon}(\hat{\lambda})\} \geq 1 - \varepsilon, \quad (41) $$

for all $i, k \in \mathbb{N}_0$.

*Proof:* At time $k=0$, we have $e_{0|0} = 0$. By construction, the sets $\mathcal{R}_i^\varepsilon(\hat{\lambda})$ are Markov-based PRS for $e_{i|0}$, i.e., $\Pr\{e_{i|0} \in \mathcal{R}_i^\varepsilon(\hat{\lambda})\} \geq 1 - \varepsilon$ for all $i \in \mathbb{N}_0$. Assume now that $\mathcal{R}_{i+k}^\varepsilon(\hat{\lambda})$ is a PRS for $e_{i|k}$ at some time $k \geq 0$. Under the initialization rule (40), the error at the next time step satisfies $e_{0|k+1} = x_{k+1} - z_{0|k+1} = \xi_{k+1}(x_{k+1} - z_{1|k}) = \xi_{k+1}e_{1|k}$, which implies $\mathbb{E}[e_{0|k+1}^\top P e_{0|k+1}] = \xi_{k+1}^2 \mathbb{E}[e_{1|k}^\top P e_{1|k}] \leq \mathbb{E}[e_{1|k}^\top P e_{1|k}]$, since $\xi_k \in [0, 1]$ for all $k \in \mathbb{N}$. As such, the same moment bound constructed for $e_{1|k}$ holds for $e_{0|k+1}$, meaning that $\mathcal{R}_{k+1}^\varepsilon$ is a valid PRS for $e_{0|k+1}$ for all $k$ via the same Markov-based construction. Using the same recursive argument as in the proof of
Proposition 1, it follows that the shifted sets $\mathcal{R}_{i+k+1}^{\epsilon}(\hat{\lambda})$ are valid PRS for $e_{i|k+1}$ for all $i \ge 0$. By induction on $k$, the claim holds for all $i, k \in \mathbb{N}_0$.

The following section establishes the main closed-loop properties of the proposed SMPC scheme under this cost function and initialization strategy.

## B. Properties of SA-SMPC

Following our development thus far, the tractable SMPC formulation for the linear discrete-time system (9) is given by the following optimization problem at each time step $k$,

$$ \min_{\mathbf{v}_k, \mathbf{z}_k, \xi_k} \sum_{i=0}^{N-1} \ell(z_{i|k}, v_{i|k}) + V_f(z_{N|k}) + \rho_k \xi_k^2 \quad (42a) $$

$$ \text{s. t.} \quad z_{i+1|k} = A z_{i|k} + B v_{i|k}, \quad (42b) $$

$$ z_{0|0} = x_0 \quad (42c) $$

$$ z_{0|k} = (1 - \xi_k)x_k + \xi_k z_{1|k-1}, \quad k \in \mathbb{N}, \quad (42d) $$

$$ v_{i|k} \in \mathcal{V}, \quad i \in \mathbb{N}_0^{N-1}, \quad (42e) $$

$$ z_{i|k} \in \mathcal{Z}_{i+k} = \mathcal{X} \ominus \mathcal{R}_{i+k}^{\epsilon}(\hat{\lambda}), \quad i \in \mathbb{N}_0^{N-1}, \quad (42f) $$

$$ z_{N|k} \in \mathcal{Z}_f \subseteq (\mathcal{X} \ominus \mathcal{R}^{\epsilon}(\hat{\lambda})) \cap \mathcal{X}_K, \quad (42g) $$

$$ \xi_k \in [0, 1]. \quad (42h) $$

Note that constraint (42c) is only meaningful for $k = 0$, it can be dropped afterwards. Moreover, notice that the reachable sets in (42f) are shifted with $k$, to maintain the feasibility (see Proposition 5 below.)

**Assumption 5 (Initial Feasibility):** We assume that the optimization problem (42) is feasible at time $k = 0$. At this initial step, the system state is exactly known, and the nominal state is initialized as $z_{0|0} = x_0$, which implies $e_0 = 0$.

We first establish that the proposed SMPC problem remains feasible over time, which is a prerequisite for all subsequent closed-loop guarantees.

**Proposition 5 (Recursive Feasibility):** If Assumption 5 is satisfied, then problem (42) is recursively feasible.

*Proof:* Assume that at time step $k$, an optimal input sequence $\mathbf{v}_k^* = \{v_{0|k}^*, \dots, v_{N-1|k}^*\}$ exists that satisfies all input constraints (42e), state constraints (42f) for $i = 0, \dots, N-1$, and the terminal constraint (42g) for the predicted terminal state $z_{N|k}$. Moreover, suppose that there exists some $\xi_k$ satisfying (42h) and (42d). Now, consider the subsequent time step $k+1$. A candidate solution can always be constructed by choosing $\xi_{k+1} = 1$, shifting the previous optimal sequence, and appending the terminal control law, i.e., $\mathbf{v}_{k+1} = \{v_{1|k}^*, \dots, v_{N-1|k}^*, K_f z_{N|k}\}$. The feasibility of this candidate sequence at $k+1$ follows from the optimality of $\mathbf{v}_k^*$, the shifting of the constraint sets as shown in (42e) and (42f), and the positive invariance of the terminal set $\mathcal{Z}_f$ under the terminal control law. Specifically, the first $N-1$ elements of $\mathbf{v}_{k+1}$ were already part of the feasible solution at time $k$ and therefore trivially satisfy the input and state constraints at time $k+1$. Furthermore, since $z_{N|k} \in \mathcal{Z}_f$ and the set $\mathcal{Z}_f$ is positively invariant under the control law $u = K_f z$, it follows that the successor state remains in $\mathcal{Z}_f$ while the associated input $K_f z_{N|k}$ satisfies the input constraints for all future time steps. Therefore, $\mathbf{v}_{k+1}$ is a suboptimal but feasible solution at time $k+1$ for all $k \in \mathbb{N}$.

Having established recursive feasibility, we now show that the proposed tightening strategy ensures satisfaction of the original chance constraints.

**Proposition 6 (Chance Constraint Satisfaction):** Given Assumptions 4 and 5, the closed-loop PRS and recursive feasibility properties of Propositions 4 and 5, the proposed SA-SMPC scheme guarantees satisfaction of the chance constraints (6e) and (6f) conditioned on the known initial state $x_0$.

*Proof:* At time $k = 0$, the initial state is exactly known. Proposition 4 guarantees that the error process satisfies

$$ \Pr \{e_{i|k} \in \mathcal{R}_{i+k}^{\epsilon}(\hat{\lambda}) | x_0\} \ge 1 - \epsilon, \quad \forall i, k \in \mathbb{N}. $$

From Proposition 5, constraint $z_{i|k} \in \mathcal{X} \ominus \mathcal{R}_{i+k}^{\epsilon}(\hat{\lambda})$, which implies that $\Pr\{x_{i|k} \in \mathcal{X} | x_0\} \ge 1 - \epsilon$, is satisfied for every $i, k \in \mathbb{N}$.
We next analyze the closed-loop performance induced by the nominal cost and initialization strategy.

**Proposition 7 (Nominal Cost Bound):** Let Assumptions 4 and 5 hold. For problem (42) defined above, the optimal cost function $\mathcal{J}_k^*$ satisfies

$$ \mathcal{J}_{k+1}^* - \mathcal{J}_k^* \leq -\ell(z_{0|k}^*, v_{0|k}^*) + \rho_{k+1}. \quad (43) $$

for all $k \geq 0$. Consequently, for any integer $\bar{N} \geq 1$,

$$ \sum_{k=0}^{\bar{N}-1} \ell(z_{0|k}^*, v_{0|k}^*) \leq \mathcal{J}_0^* + \sum_{k=1}^{\bar{N}} \rho_k. \quad (44) $$

*Proof:* At time $k+1$, consider the shifted candidate sequences and interpolating variable

$$ z_{i|k+1} = z_{i+1|k}^*, \quad v_{i|k+1} = v_{i+1|k}^*, \quad \xi_{k+1} = 1, $$

with terminal input $v_{N-1|k+1} = K_f z_{N|k}^*$ and state $z_{N|k+1} = (A + BK_f)z_{N|k}^*$. These candidate sequences are feasible. The associated cost satisfies

$$ \tilde{\mathcal{J}}_{k+1} \leq \mathcal{J}_k^* - \ell(z_{0|k}^*, v_{0|k}^*) - \rho_k \xi_k^2 + \rho_{k+1}, $$

where the inequality

$$ \ell(z_{N-1|k+1}, v_{N-1|k+1}) + V_f(z_{N|k+1}) = \|z_{N|k}^*\|_Q^2 + \|K_f z_{N|k}^*\|_R^2 + \|(A + BK_f)z_{N|k}^*\|_S^2 \leq \|z_{N|k}^*\|_S^2, $$

which follows from (37), has been used. Since $\mathcal{J}_{k+1}^* \leq \tilde{\mathcal{J}}_{k+1}$ and $\rho_k \xi_k^2 \geq 0$, the cost decrease relationship in (43) follows. Summing the result from $k=0$ to $\bar{N}-1$ and using $\mathcal{J}_{\bar{N}}^* \geq 0$ yields the cost bound in (44). $\blacksquare$

The following result shows that the nominal cost decrease implies convergence of the nominal state and input.

**Proposition 8 (Nominal State and Input Convergence):** Let Assumptions 4 and 5 hold. If $\sum_{k=1}^{\infty} \rho_k < \infty$, then $\lim_{k \to \infty} z_{0|k}^* = 0$ and $\lim_{k \to \infty} v_{0|k}^* = 0$.

*Proof:* From Proposition 7, and $\sum_{k=1}^{\infty} \rho_k < \infty$, the right-hand side of (44) is uniformly bounded, implying $\sum_{k=0}^{\infty} \ell(z_{0|k}^*, v_{0|k}^*) < \infty$. Since $\ell(\cdot, \cdot)$ is positive definite, the result holds. $\blacksquare$

Finally, we relate the convergence of the nominal system to mean-square stability of the true closed-loop state.

**Proposition 9 (Mean-square Stability of $x_k$):** Let Assumptions 4 and 5 hold. Given $x_k = z_k + e_k$, and $z_k \to 0$, the true state $x_k$ is mean-square stable, satisfying

$$ \lim_{k \to \infty} \mathbb{E} \left[ \|x_k\|_Q^2 \right] \leq \Lambda_{\max} \left( P^{-1/2} Q P^{-1/2} \right) \frac{1}{1 - \hat{\lambda}} \mathrm{Tr}(PW), \quad (45) $$

*Proof:* From Proposition 3 and $Q \leq \Lambda_{\max} (P^{-1/2} Q P^{-1/2}) P$, we can derive the following bound

$$ \mathbb{E} \left[ \|e_k\|_Q^2 \right] \leq \Lambda_{\max} \left( P^{-1/2} Q P^{-1/2} \right) \frac{1 - (\hat{\lambda})^k}{1 - \hat{\lambda}} \mathrm{Tr}(PW). $$

Using $x_k = z_k + e_k$, we have $\|x_k\|_Q^2 = \|z_k\|_Q^2 + \|e_k\|_Q^2 + 2z_k^\top Q e_k$. Taking expectations and applying the Cauchy-Schwarz inequality yields

$$ \mathbb{E} \left[ \|x_k\|_Q^2 \right] \leq \|z_k\|_Q^2 + \mathbb{E} \left[ \|e_k\|_Q^2 \right] + 2 \sqrt{\|z_k\|_Q^2} \mathbb{E} \left[ \|e_k\|_Q^2 \right]. $$

Hence,

$$ \mathbb{E} \left[ \|x_k\|_Q^2 \right] \leq \left( \sqrt{\|z_k\|_Q^2} + \sqrt{\mathbb{E} \left[ \|e_k\|_Q^2 \right]} \right)^2. $$

Using $z_k \to 0$ from Proposition 8 and taking the limit as $k \to \infty$ gives the stated result. $\blacksquare$

The algorithm resulting from the optimization problem in (42) and the preceding development is presented below.
**Algorithm 1 Offline Design and Online Execution**

**Require:** System matrices *A*, *B*, constraints *X*, *U*, noise covariance *W*, cost matrices *Q*, *R*, violation level ε, parameter *v*<sub>ss</sub>.

**Offline Design**

1. Solve (23) for contraction rates $\lambda, \lambda_L$, Lyapunov matrix $P$, and feedback gain $K$.

2. Determine effective contraction rate $\hat{\lambda}$ corresponding to $v_{ss}$.

3. Compute tightened state constraints $Z_i$, and terminal set $Z_f$.

**Online Execution**

4. Initialize $z_{0|0} = z_0 \leftarrow x_0$ (measured-state initialization).

5. **for** $k = 0, 1, \dots$ **do**

6. Measure real state $x_k$.

7. Solve (42) for variable $\xi_k^*$ and corresponding sequences $z_k^*, v_k^*$.

8. Apply $u_k \leftarrow \varphi(v_{0|k}^* + K(x_k - z_{0|k}^*))$.

9. Store $z_{1|k}^*$ for next iteration.

10. **end for**

Enforcing (42e) ensures the nominal input $v_{i|k}$ is admissible, but cannot guarantee that the total control input $u_k = v_{0|k}^* + Ke_k$ satisfies the original input constraints (6d). However, by accounting for the saturation nonlinearity into the predictive model, the constructed PRS reflect the true closed-loop behavior of the system, thus enabling the controller to anticipate the potential saturation and plan against its effects. This means that the derived closed-loop properties are valid under the physically implementable control law $\varphi(u_k)$. In contrast, standard constraint-tightening SMPC formulations such as [13], establish guarantees under the idealized linear control law $u_k = v_k + Ke_k$, which becomes physically unrealizable when saturation occurs, thereby invalidating theoretical assurances of state constraint satisfaction and stability. Our proposed method replaces this weak assumption with a strong, verifiable guarantee of robust state safety and stability under real operating conditions.

**Remark 3 (Input Constraints Handling):** Note that $r_L$ depends, in general, on the nominal input $v$, see (24), which directly influences $\bar{\lambda}^*$. This introduces a circular dependency: large nominal inputs shrink the region of linearity and lead to conservative PRS, yet computing tight PRS requires prior knowledge of $v$, which itself depends on the resulting tightened constraints. A naïve solution is to rely solely on the open-loop induced PRS $\mathcal{R}_i^\varepsilon(\lambda)$ for the constraint tightening. While tractable and independent of $v$, this approach may be highly conservative and cause infeasibility. Similarly, although the strategy employed, by imposing Assumption 4, supports the theoretical development by ensuring the existence of an effective contraction rate, it may become overly restrictive in practice. In particular, for moderate or large disturbance levels, enforcing a constant input bound $v_{ss}$ along the prediction horizon can excessively limit the nominal input and render the nominal trajectory infeasible. An alternative strategy would consist in adopting a structured heuristic in which the nominal input is parameterized offline as a monotonically decreasing sequence $\bar{v}_i$, ranging from near-saturation toward the steady-state bound $v_{ss}$. The decay rate may be tuned according to design preferences. For instance, following a robust-style argument, $\bar{v}_i$ can be chosen to decrease sufficiently fast to reserve adequate control authority for the feedback action associated with the worst probabilistic error evolution characterized by $\mathcal{R}_i^\varepsilon(\lambda)$, while maintaining $\bar{v}_i \ge v_{ss}$. From a set-theoretic perspective, this corresponds to tightening the input constraints by bounded scalings $K\mathcal{R}_i^\varepsilon(\lambda)$. Related ideas appear in the literature; for instance, [12] enforces a saturation on the disturbance to preserve a non-empty set of admissible nominal inputs.

V. NUMERICAL EXAMPLE

This section illustrates the practical applicability of the proposed SA-SMPC framework on a representative system. All simulations were carried out in Python using the `CVXPY` modeling interface [36], [37] in conjunction with the `CLARABEL` solver [38]. The computations were executed on a Linux workstation equipped with an AMD Ryzen 7700 CPU (16 cores) and 32GB of RAM.

For comparison purposes, we evaluate the proposed method with $\mathcal{R}_i^\varepsilon(\bar{\lambda}^*)$-based constraint tightening against two alternatives: the conservative tightening using $\mathcal{R}_i^\varepsilon(\lambda)$, and the affine disturbance-feedback parameterization approach proposed in [12] with fixed uniform risk allocation. The latter is particularly relevant, as it represents an alternative
analytical SMPC framework capable of enforcing hard input constraints under unbounded disturbances. Since it results in a second-order cone program (SOCP) solved online at each time step, we refer to it as SOCP-SMPC. The three approaches will thus be denoted as $\bar{\lambda}^*$-SA-SMPC, lambda-SA-SMPC, and SOCP-SMPC, respectively.

**Remark 4 (Uniform Risk Allocation):** It is worth noting that [12] also discusses the possibility of optimizing the risk allocation across constraints. However, the joint optimization of the feedback gain and the risk allocation is inherently nonconvex. Although convex approximations have been proposed, they substantially increase computational effort and do not guarantee global optimality. The iterative procedure suggested in [12] alternates between optimizing the risk allocation for a fixed feedback gain and optimizing the feedback gain for a fixed risk allocation, repeating this process for a prescribed number of iterations. In practice, each iteration effectively entails solving an additional SOCP of comparable complexity, leading to a significant increase in computational burden. For this reason, we adopt a fixed uniform risk allocation in the SOCP-SMPC implementation.

## A. Setup

We consider an isothermal continuous-stirred-tank reactor (CSTR) as in [39], where the states $x_1$ and $x_2$ denote the concentrations of species A and B, respectively, and the control input $u$ is the dilution rate. The control objective is to regulate $x_2$ to the setpoint $x_2^* = 1$, which corresponds to the steady-state operating point $(x_1^*, u^*) = (2.5, 25)$. Linearizing the nonlinear CSTR dynamics about this equilibrium and discretizing with sampling time $T_s = 0.002$ s yields the discrete-time deviation model

$$x_{k+1} = Ax_k + Bu_k + w_k,$$

with

$$A = \begin{bmatrix} 0.95123 & 0 \\ 0.08833 & 0.81873 \end{bmatrix}, \quad B = \begin{bmatrix} -0.0048771 \\ -0.0020429 \end{bmatrix}.$$

This model coincides with [39]; in the present numerical study, the cost function weights, disturbance statistics, state constraints and initial conditions are adapted as specified below. Process disturbances are modeled as i.i.d. Gaussian noise $w_k \sim N(0, W)$ with $W = 0.003^2 I_2$. We fix $Q = \text{diag}(20, 100)$ and $R = 0.1$, and enforce hard input bounds $-25 \le \tilde{u}_k \le 25$. We will treat 4 different variations of the problem with changing state constraints and initial conditions (both specified in deviation variables). The considered variations are intentionally selected to challenge both feasibility and performance margins of the different control strategies. Specifically,

* Scenario 1: $\text{Pr}\{\tilde{x}_2 \le 0.25\} \ge 0.8, \tilde{x}_0 = [0.5, 0.2]^\top,$

* Scenario 2: $\text{Pr}\{\tilde{x} \le [0.75, 0.25]^\top\} \ge 0.8, \text{Pr}\{2\tilde{x}_1 + \tilde{x}_2 \le 1.5\} \ge 0.8, \tilde{x}_0 = [0.5, 0.2]^\top,$

* Scenario 3: $\text{Pr}\{[-0.5, -0.25]^\top \le \tilde{x} \le [0.75, 0.25]^\top\} \ge 0.8, \text{Pr}\{2\tilde{x}_1 + \tilde{x}_2 \le 1.5\} \ge 0.8, \tilde{x}_0 = [0.5, 0.2]^\top,$

* Scenario 4: $\text{Pr}\{\tilde{x}_2 \le 0.15\} \ge 0.8, \tilde{x}_0 = [0.5, 0.12]^\top.$

As such, we set a $\varepsilon = 0.2$ for the SA-SMPC methods, and a risk allocation $\varepsilon_i = 0.2/|\mathbb{N}_1^{nc}|$, where nc denotes the number of constraints, for the SOCP-SMPC strategy. Note that while the constraints on the combined states introduced in variations 2 and 3 may not always correspond to physically meaningful scenarios, they provide a useful benchmark for assessing the feasibility limits of any given SMPC approach. All scenarios are run with a prediction horizon $N = 15$.

## B. Offline Design for SA-SMPC

For all variations of the problem, the offline design process of the SA-SMPC scheme is identical, and yields the same results, with only the terminal set $\mathcal{Z}_f$ changing due to the addition of constraints. Solving the offline optimization problem in (23), we obtain an open-loop contraction $\lambda \approx 0.9049$, a linear contraction rate $\lambda_L \approx 0.67035$, along with the feedback gain $K$ and matrix $P$ as specified below

$$K = \begin{bmatrix} 27.1573 & 0.09622 \end{bmatrix}, \quad P = \begin{bmatrix} 71.2230 & -0.4498 \\ -0.4498 & 1.01712 \end{bmatrix}.$$

Given that the noise magnitude is relatively low, we set $v_{i|k} \le v_{ss} = 24$ for all $i, k$ following Assumption 4, for which we compute an effective contraction rate $\bar{\lambda}^* \approx 0.6752$. The constraint tightening that follows is standard, as described in Section IV.
<table><caption>TABLE I
COMPARISON OF $\bar{\lambda}^*$-SA-SMPC, $\lambda$-SA-SMPC, AND SOCP-SMPC METHODS ON THE BASIS OF AVERAGE COST AND COMPUTATIONAL TIME FOR EACH OF THE FOUR STUDIED SCENARIOS.</caption><thead><tr><th rowspan="2">Scenario</th><th colspan="2">Novel $\bar{\lambda}^*$-SA-SMPC</th><th colspan="2">$\lambda$-SA-SMPC</th><th colspan="2">SOCP-SMPC</th></tr><tr><th>E<sub>w</sub> [$\mathcal{J}_{MPC}$] ± $\sigma_{MPC}$</th><th>$\bar{\tau}$ (ms/iter)</th><th>E<sub>w</sub> [$\mathcal{J}_{MPC}$] ± $\sigma_{MPC}$</th><th>$\bar{\tau}$ (ms/iter)</th><th>E<sub>w</sub> [$\mathcal{J}_{MPC}$] ± $\sigma_{MPC}$</th><th>$\bar{\tau}$ (ms/iter)</th></tr></thead><tbody><tr><td>Scenario 1</td><td>77.1 ± 1.95</td><td>1.62</td><td>79.2 ± 2.11</td><td>1.65</td><td>67.0 ± 1.84</td><td>52</td></tr><tr><td>Scenario 2</td><td>77.1 ± 1.95</td><td>1.57</td><td>79.2 ± 2.11</td><td>1.63</td><td>88.8 ± 1.55</td><td>146</td></tr><tr><td>Scenario 3</td><td>77.1 ± 1.95</td><td>1.64</td><td>79.2 ± 2.11</td><td>1.71</td><td>Infeasible</td><td>Infeasible</td></tr><tr><td>Scenario 4</td><td>114.0 ± 2.27</td><td>1.59</td><td>Infeasible</td><td>Infeasible</td><td>57.5 ± 1.77</td><td>57</td></tr></tbody></table>

### C. Comparison and Discussion

The strategies are compared over 1000 trajectories on the basis of feasibility, average computational time $\bar{\tau}$ per iteration, and performance in terms of the index $\mathcal{J}_{MPC} = \sum_{k=0}^{T} (\|\tilde{x}_k\|_Q^2 + \|\tilde{u}_k\|_R^2)$. The main results are summarized in table I.

For the two variations in which all strategies are feasible, several trends emerge. In scenario 1, the SOCP-SMPC scheme exhibits superior performance, primarily due to its ability to optimize the feedback gain online and thereby accurately capture the effect of low noise magnitudes. This performance gain, however, comes at a substantial computational expense, as the SOCP-SMPC formulation is approximately 35 times more demanding than either SA-SMPC scheme. We also note that the reduced conservatism introduced by the refined tightening enables $\bar{\lambda}^*$-SA-SMPC to outperform its open-loop counterpart.

Scenario 2 augments the constraint set of scenario 1 with two additional constraints, leading to a significant degradation in the performance of SOCP-SMPC. This deterioration stems from the treatment of constraints as disjoint events via Boole's inequality (see [12]), which effectively assigns each of the three constraints a violation probability of $0.2/3 \approx 0.667$, thereby inducing substantially more conservative behavior. As mentioned in Remark 4, [12] proposes an online iterative procedure for risk allocation optimization per individual constraint that may mitigate this effect, such an approach further increases the computational effort and is unlikely to remain effective for higher-dimensional systems. Consistent with this observation, the computational time of SOCP-SMPC increases nearly three-fold to accommodate the additional constraints, and thus becomes around 100 times more computationally demanding than any of the SA-SMPC approaches. In contrast, both SA-SMPC methods display identical performance to that observed in Scenario 1, outperforming the SOCP-SMPC scheme, with no change in computational effort. The performance deterioration trend for the SOCP-SMPC scheme devolves into infeasibility in Scenario 3, where the addition of more constraints render the risk allocation too tight, as reported in table I.

In scenario 4, the $\lambda$-SA-SMPC scheme becomes infeasible due to the tighter constraint on $x_2$. Once again, SOCP-SMPC achieves superior performance. By comparison, the conservatism of the $\bar{\lambda}^*$-SA-SMPC formulation forces the mean trajectory to remain excessively backed off from the constraints, thereby requiring a greater dilution rate which increases the cost function value. This can be seen more clearly in figure 1 displaying the resulting input profiles of the two strategies.

It is worth noting that, although not explicitly illustrated in the presented results, the SOCP-SMPC scheme is expected to remain feasible under higher noise magnitudes and for initial conditions arbitrarily close to the constraint boundaries, whereas the SA-SMPC schemes lose feasibility in such scenarios.

Overall, our proposed SA-SMPC framework is computationally attractive, exhibiting low complexity and favorable scalability with respect to both system dimension and the number of constraints. The effective tightening of $\bar{\lambda}^*$-SA-SMPC further reduces conservatism, resulting in improved performance and enlarged feasibility regions relative to $\lambda$-SA-SMPC. However, this computational efficiency is achieved at the expense of increased conservatism, in no small part due to the use of the loose Markov inequality and fixed reachable set geometries, which leads to overly cautious constraint back-offs and reduced performance and feasibility in certain scenarios. The principal advantage of SOCP-SMPC lies in its nonconservative back-off strategy, which supports reduced costs, larger feasibility regions under tighter constraints and elevated noise levels, and satisfactory constraint violation frequencies. Nevertheless, this benefit is offset by poor computational scalability and increased sensitivity to the number of constraints, suggesting limited suitability for high-dimensional applications.
<figure><img src="image_5.png" /><figcaption>Fig. 1. Trajectory of state $x_2$ in scenario 4 under (a) $\bar{\lambda}^*$-SA-SMPC and (b) SOCP-SMPC. The inset in (a) shows the corresponding input profile for the $\bar{\lambda}^*$-SA-SMPC strategy, which exploits the full admissible input range. The inset in (b) reports the input profile for SOCP-SMPC, which operates over a restricted subset of the admissible inputs.</figcaption></figure>

## VI. CONCLUSION

We presented a tractable stochastic tube-based SMPC framework that guarantees hard input-constraint satisfaction for linear systems under unbounded disturbances. By incorporating actuator saturation directly into the prediction model and constructing PRS via sector bounds and convex embeddings, we obtain a control scheme that preserves a linear nominal optimization while appropriately characterizing the resulting nonlinear error dynamics. This saturation-based framework provides a systematic probabilistic link between linear and saturated regimes, enabling tighter PRS and improved performance. As demonstrated in the numerical study, the proposed approach outperforms schemes based on naïve open-loop tightenings, while remaining significantly more computationally efficient than its feedback-parameterized counterpart, albeit at the cost of increased conservatism. Future work aims to further reduce said conservatism by using the moment sum-of-squares hierarchy and polynomial optimization tools to incorporate higher-order disturbance moments into the PRS design, enabling the construction of semialgebraic PRS certified by moment-based probabilistic certificates that are less conservative than Markov-type bounds.

# APPENDIX I
## PROOF OF PROPOSITION 3

Given $\lambda - \lambda_L \ge 0$, inequality (26) implies that the globally valid upper bound is

$$ \mathbb{E}[e_{i+1|k}^{\top} Pe_{i+1|k}] \le \lambda \mathbb{E}[e_{i|k}^{\top} Pe_{i|k}] + \mathrm{Tr}(PW), $$

which holds when $\mathrm{Pr}\{e_{i|k} \in S_1\} = 0$. Iterating this recursion yields

$$ \mathbb{E}[e_{i|k}^{\top} Pe_{i|k}] \le \frac{1-\lambda^i}{1-\lambda} \mathrm{Tr}(PW) \le \frac{1}{1-\lambda} \mathrm{Tr}(PW), $$

for all $i, k \in \mathbb{N}$ and every $e_{i|k} \in \mathbb{R}^n$. Condition (28) guarantees the existence of a parameter $\mu \in [\lambda_L, \lambda]$ such that

$$ \frac{\mathbb{E}[e_{i|k}^{\top} Pe_{i|k}]}{r_L} \le \frac{\mu - \lambda_L}{\lambda - \lambda_L} \le 1. \quad (46) $$

Substituting this parametrized bound into (25) produces

$$ \mathbb{E}[e_{i+1|k}^{\top} Pe_{i+1|k}] \le \mu \mathbb{E}[e_{i|k}^{\top} Pe_{i|k}] + \mathrm{Tr}(PW), \quad (47) $$
which mirrors (16) with parameter $\mu$. Iteration then gives

$$ \mathbb{E}[e_{i|k}^{\top} Pe_{i|k}] \le \frac{1-\mu^i}{1-\mu} \mathrm{Tr}(PW) \le \frac{1}{1-\mu} \mathrm{Tr}(PW), \quad (48) $$

provided (46) holds for all $i, k \in \mathbb{N}$, thereby ensuring (48). We now identify the values of $\mu$ for which (46) holds universally, thus validating (48). Define the function

$$ g(\mu) := \frac{\mu - \lambda_L}{\lambda - \lambda_L} r_L - \frac{1}{1-\mu} \mathrm{Tr}(PW), \quad (49) $$

which is continuous and strictly concave on the interval $[\lambda_L, \lambda]$. From (28), $g(\lambda) > 0$. At $\mu = \bar{\lambda}^*$, defined by (30), $g(\bar{\lambda}^*) = 0$. The strict concavity of $g$ implies $g(\mu) > 0$ for all $\mu \in (\bar{\lambda}^*, \lambda]$. Equivalently,

$$ \frac{1}{1-\mu} \mathrm{Tr}(PW) < \frac{\mu - \lambda_L}{\lambda - \lambda_L} r_L, \quad \forall \mu \in (\bar{\lambda}^*, \lambda]. \quad (50) $$

Now, fix $\mu \in (\bar{\lambda}^*, \lambda]$. Three potential cases arise for $\mathbb{E}[e_{i|k}^{\top} Pe_{i|k}]$:

$$ \text{1) } \mathbb{E}[e_{i|k}^{\top} Pe_{i|k}] \le \frac{1}{1-\mu} \mathrm{Tr}(PW), $$

$$ \text{2) } \frac{1}{1-\mu} \mathrm{Tr}(PW) < \mathbb{E}[e_{i|k}^{\top} Pe_{i|k}] \le \frac{\mu - \lambda_L}{\lambda - \lambda_L} r_L, $$

$$ \text{3) } \frac{\mu - \lambda_L}{\lambda - \lambda_L} r_L < \mathbb{E}[e_{i|k}^{\top} Pe_{i|k}]. $$

Case 2 leads to a contradiction, as it violates the implication that (46) leads to (48). Case 3 is also impossible: by the concavity of $g$, if $\mathbb{E}[e_{i|k}^{\top} Pe_{i|k}]$ exceeded $\frac{\mu-\lambda_L}{\lambda-\lambda_L} r_L$, there would exist a $\nu \in (\mu, \lambda]$ such that

$$ \frac{1}{1-\nu} \mathrm{Tr}(PW) < \mathbb{E}[e_{i|k}^{\top} Pe_{i|k}] = \frac{\nu - \lambda_L}{\lambda - \lambda_L} r_L, $$

where the left-hand inequality follows from (50). This is equivalent to Case 2; a contradiction. Therefore, only Case 1 is possible. Similar reasoning applies for $\mu = \bar{\lambda}^*$. This demonstrates that any $\mu \in [\bar{\lambda}^*, \lambda]$ ensures (47) and (48) for all $i, k \in \mathbb{N}$. The smallest such $\mu$, i.e., that for which $g(\mu) = 0$, is $\mu = \bar{\lambda}^*$ and yields the tightest valid upper bound

$$ \mathbb{E}[e_{i|k}^{\top} Pe_{i|k}] \le \frac{1 - (\bar{\lambda}^*)^i}{1 - \bar{\lambda}^*} \mathrm{Tr}(PW). $$

This establishes (29) and concludes the proof.

■

## APPENDIX II
CONVERGENCE AND MEAN-SQUARE STABILITY UNDER THE REAL COST

**Proposition 10 (Average Cost Bound.):** Let Assumptions 2 and 4 hold, and let the MPC cost function be defined as in (38). For $\delta > 0$ and terminal weight $S = \alpha P$ with $\alpha \ge 1$ and $\delta \le (1-\hat{\lambda})/\hat{\lambda}$ such that

$$ \alpha (1 - (1 + \delta) \cdot \hat{\lambda}) P \succeq Q + K^{\top} R K, \quad (51) $$

the optimal cost function $\mathcal{J}_k^*$ satisfies

$$ \mathbb{E}[\mathcal{J}_{k+1}^*] - \mathbb{E}[\mathcal{J}_k^*] \le -\ell(x_{0|k}^*, \varphi(u_{0|k}^*)) + \left(1 + \frac{1}{\delta}\right) \frac{4}{1-\hat{\lambda}} \mathrm{Tr}(SW), \quad (52) $$

for all $k \ge 0$, and consequently

$$ \lim_{\bar{N} \to \infty} \frac{1}{\bar{N}} \sum_{k=0}^{\bar{N}-1} \mathbb{E} \left[ \ell \left( x_{0|k}^*, \varphi(u_{0|k}^*) \right) \right] \le \left( 1 + \frac{1}{\delta} \right) \frac{4}{1-\hat{\lambda}} \mathrm{Tr}(SW). \quad (53) $$
*Proof:* Consider the shifted candidate, whose cost is denoted by $\tilde{\mathcal{J}}_{k+1}$, obtained by applying the standard shift, and defined as: $\tilde{x}_{i|k+1} = x_{i+1|k}^*$ for $i \in \mathbb{N}_0^{N-1}$; $\tilde{u}_{i|k+1} = u_{i+1|k}^*$ for $i \in \mathbb{N}_0^{N-2}$; $\tilde{x}_{N|k+1} = A_K z_{N|k}^* + e_{N+1|k}$ and $\tilde{u}_{N-1|k+1} = K x_{N|k}^*$, where the dependence on the realization of $w_k$ is left implicit. Then

$$ 
\begin{aligned}
\mathbb{E}[\tilde{\mathcal{J}}_{k+1}] &= \mathbb{E}\left[\sum_{i=0}^{N-1} \|\tilde{x}_{i|k+1}\|_Q^2 + \|\varphi(\tilde{u}_{i|k+1})\|_R^2\right] + \mathbb{E}\left[\|\tilde{x}_{N|k+1}\|_S^2\right] \\
&= \mathbb{E}\left[\sum_{i=1}^{N-1} \|x_{i|k}^*\|_Q^2 + \|\varphi(u_{i|k}^*)\|_R^2\right] + \mathbb{E}\left[\|x_{N|k}^*\|_Q^2 + \|\varphi(Kx_{N|k}^*)\|_R^2\right] + \mathbb{E}\left[\|A_K z_{N|k}^* + e_{N+1|k}\|_S^2\right].
\end{aligned}
$$

The optimal cost at time $k$ satisfies

$$ 
\begin{aligned}
\mathbb{E}[\mathcal{J}_k^*] &= \mathbb{E}\left[\sum_{i=0}^{N-1} \|x_{i|k}^*\|_Q^2 + \|\varphi(u_{i|k}^*)\|_R^2\right] + \mathbb{E}\left[\|x_{N|k}^*\|_S^2\right] \\
&= \|x_{0|k}^*\|_Q^2 + \|\varphi(u_{0|k}^*)\|_R^2 + \mathbb{E}\left[\sum_{i=1}^{N-1} \|x_{i|k}^*\|_Q^2 + \|\varphi(u_{i|k}^*)\|_R^2\right] + \mathbb{E}\left[\|x_{N|k}^*\|_S^2\right],
\end{aligned}
$$

Subtracting the two expressions yields the cost difference

$$ 
\begin{aligned}
\mathbb{E}[\tilde{\mathcal{J}}_{k+1}] - \mathbb{E}[\mathcal{J}_k^*] &= - \|x_{0|k}^*\|_Q^2 - \|\varphi(u_{0|k}^*)\|_R^2 - \mathbb{E}[\|x_{N|k}^*\|_S^2] + \mathbb{E}[\|x_{N|k}^*\|_Q^2 + \|\varphi(Kx_{N|k}^*)\|_R^2] \\
&\quad + \mathbb{E}[\|A_K z_{N|k}^* + e_{N+1|k}\|_S^2] \\
&\leq - \|x_{0|k}^*\|_Q^2 - \|\varphi(u_{0|k}^*)\|_R^2 - \mathbb{E}[\|x_{N|k}^*\|_S^2] + \mathbb{E}[\|x_{N|k}^*\|_Q^2 + \|Kx_{N|k}^*\|_R^2] \\
&\quad + \mathbb{E}[\|A_K x_{N|k}^* + e_{N+1|k} - A_K e_{N|k}\|_S^2], \\
&\leq - \|x_{0|k}^*\|_Q^2 - \|\varphi(u_{0|k}^*)\|_R^2 + \mathbb{E}\left[(x_{N|k}^*)^\top (Q + K^\top R K + (1+\delta)A_K^\top S A_K - S)(x_{N|k}^*)\right], \\
&\quad + \left(1 + \frac{1}{\delta}\right) \mathbb{E}\left[\|e_{N+1|k} - A_K e_{N|k}\|_S^2\right],
\end{aligned}
$$

where the first inequality follows from the fact that $\|\varphi(u)\|_R^2 \le \|u\|_S^2$ for all $u$ and the fact that $z_{N|k}^* = x_{N|k}^* - e_{N|k}$. The second inequality is obtained by applying Young's inequality, i.e. $2xy \le \delta x^2 + (1/\delta)y^2$ for all $x, y \in \mathbb{R}$, in fact

$$ \|a+b\|_S^2 \le \|a\|_S^2 + \|b\|_S^2 + 2 \|a\|_S \|b\|_S \le (1+\delta) \|a\|_S^2 + \left(1+\frac{1}{\delta}\right) \|b\|_S^2, $$

with $\delta > 0$, to the term $\mathbb{E}\left[\|A_K x_{N|k}^* + e_{N+1|k} - A_K e_{N|k}\|_S^2\right]$.

Note that for every $\delta \le (1 - \hat{\lambda})/\hat{\lambda}$, equivalent to $1 - (1 + \delta)\hat{\lambda} \ge 0$, condition (51) holds for an $\alpha$ big enough. From this and $A_K^\top P A_K \le \hat{\lambda} P$, it follows that

$$ S - (1 + \delta) A_K^\top S A_K = \alpha \left( P - (1 + \delta) \cdot A_K^\top P A_K \right) \succeq \alpha \left( 1 - (1 + \delta) \cdot \hat{\lambda} \right) P \succeq Q + K^\top R K $$

which allows us to drop the terms in $x_{N|k}^*$ and get

$$ \mathbb{E}[\tilde{\mathcal{J}}_{k+1}] - \mathbb{E}[\mathcal{J}_k^*] \le - \|x_{0|k}^*\|_Q^2 - \|\varphi(u_{0|k}^*)\|_R^2 + \left(1 + \frac{1}{\delta}\right) \mathbb{E}\left[\|e_{N+1|k} - A_K e_{N|k}\|_S^2\right]. $$
From the parallelogram law, i.e. $\|a+b\|^2 \le 2\|a\|^2 + 2\|b\|^2$, we obtain

$$
\begin{aligned}
\mathbb{E} \left[ \|e_{N+1|k} - A_K e_{N|k}\|_S^2 \right] &\le 2 \mathbb{E} \left[ \|e_{N+1|k}\|_S^2 \right] + 2 \mathbb{E} \left[ \|e_{N|k}\|_{A_K^\top S A_K}^2 \right] \\
&\le 2 \left( \hat{\lambda} \mathbb{E} \left[ \|e_{N|k}\|_S^2 \right] + \text{Tr}(SW) \right) + 2 \mathbb{E} \left[ \|e_{N|k}\|_{A_K^\top S A_K}^2 \right] \\
&\le 2 \left( (1+\hat{\lambda}) \mathbb{E} \left[ \|e_{N|k}\|_S^2 \right] + \text{Tr}(SW) \right) \le 2 \left( \frac{1+\hat{\lambda}}{1-\hat{\lambda}} \text{Tr}(SW) + \text{Tr}(SW) \right) \\
&= \frac{4}{1-\hat{\lambda}} \text{Tr}(SW).
\end{aligned}
$$

Using the above result we obtain

$$
\begin{aligned}
\mathbb{E}[\tilde{\mathcal{J}}_{k+1}] - \mathbb{E}[\mathcal{J}_k^*] &\le - \|x_{0|k}^*\|_Q^2 - \|\varphi(u_{0|k}^*)\|_R^2 + \left(1+\frac{1}{\delta}\right) \frac{4}{1-\hat{\lambda}} \text{Tr}(SW) \\
&= -\ell(x_{0|k}^*, \varphi(u_{0|k}^*)) + \left(1+\frac{1}{\delta}\right) \frac{4}{1-\hat{\lambda}} \text{Tr}(SW).
\end{aligned}
$$

Finally, by optimality we have $\mathbb{E}[\mathcal{J}_{k+1}^*] \le \mathbb{E}[\tilde{\mathcal{J}}_{k+1}]$, hence the same inequality holds with $\mathcal{J}_{k+1}^*$ on the left-hand side. Taking expectations and summing from $k=0$ to $\bar{N}-1$ yields

$$
\mathbb{E}[\mathcal{J}_{\bar{N}}^*] - \mathbb{E}[\mathcal{J}_0^*] \le - \sum_{k=0}^{\bar{N}-1} \mathbb{E} \left[ \ell \left( x_{0|k}^*, \varphi \left( u_{0|k}^* \right) \right) \right] + \bar{N} \left( 1 + \frac{1}{\delta} \right) \frac{4}{1-\hat{\lambda}} \mathrm{Tr}(SW),
$$

and since $\mathbb{E}[\mathcal{J}_{\bar{N}}^*] \ge 0$, dividing by $\bar{N}$ and letting $\bar{N} \to \infty$ gives the result in (53). ■

**Remark 5:** The choice of $S = \alpha P$ with $\alpha \ge 1$ is made only to obtain a clean bound proportional to $\text{Tr}(SW)$. The argument remains valid for a generic terminal weight $S \succ 0$ satisfying a contraction condition of the form $S \succeq (1+\delta) \cdot A_K^\top S A_K + Q + K^\top R K$, which would introduce scaling constants into the final result. We restrict ourselves to $S = \alpha P$ for clarity and compactness of the final bound.

**Remark 6:** A more interpretable, albeit mathematically moot, way of looking at the cost bound above may be obtained by close inspection of the term $\mathbb{E}[\|e_{N+1|k} - A_K e_{N|k}\|_S^2]$. Under the terminal law $v_{N|k} = K z_{N|k}$ (hence $u_{N|k} = K x_{N|k}$), one has

$$
\begin{aligned}
\mathbb{E} \left[ \| e_{N+1|k} - A_K e_{N|k} \|_S^2 \right] &= \mathbb{E} \left[ \| A e_{N|k} + B \left( \varphi (K x_{N|k}) - K z_{N|k} \right) + w_{N+k} - A e_{N|k} - B K e_{N|k} \|_S^2 \right], \\
&= \mathbb{E} \left[ \| B \left( \varphi (K x_{N|k}) - K x_{N|k} \right) + w_{N+k} \|_S^2 \right].
\end{aligned}
$$

Since $w_{N+k}$ is zero mean and independent of $x_{N|k}$ by Assumption 1,

$$
\begin{aligned}
\mathbb{E} \left[ \| e_{N+1|k} - A_K e_{N|k} \|_S^2 \right] &= \mathbb{E} \left[ \| B \left( \varphi (K x_{N|k}) - K x_{N|k} \right) \|_S^2 \right] + \mathbb{E} \left[ \| w_{N+k} \|_S^2 \right], \\
&= \mathbb{E} \left[ \| B \left( \varphi (K x_{N|k}) - K x_{N|k} \right) \|_S^2 \right] + \text{Tr}(SW).
\end{aligned}
$$

Hence, the cost bound can be interpreted as being driven by the saturation mismatch and the disturbance covariance, which aligns with our physical intuition and the established results on nonlinear SMPC stability [40].

## REFERENCES

[1] D. Mayne, M. Seron, and S. Raković, “Robust model predictive control of constrained linear systems with bounded disturbances,” *Automatica*, vol. 41, no. 2, p. 219–224, 2005.

[2] M. Cannon, J. Buerger, B. Kouvaritakis, and S. Raković, “Robust tubes in nonlinear model predictive control,” *IFAC Proceedings Volumes*, vol. 43, no. 14, pp. 208–213, 2010.

[3] S. V. Raković, B. Kouvaritakis, R. Findeisen, and M. Cannon, “Homothetic tube model predictive control,” *Automatica*, vol. 48, no. 8, p. 1631–1638, 2012.

[4] M. Farina, L. Giulioni, and R. Scattolini, “Stochastic linear model predictive control with chance constraints – a review,” *Journal of Process Control*, vol. 44, p. 53–67, Aug. 2016.
[5] A. Mesbah, "Stochastic model predictive control: An overview and perspectives for future research," *IEEE Control Systems Magazine*, vol. 36, no. 6, pp. 30–44, 2016.

[6] L. Hewing and M. N. Zeilinger, "Scenario-based probabilistic reachable sets for recursively feasible stochastic model predictive control," *IEEE Control Systems Letters*, vol. 4, no. 2, pp. 450–455, 2020.

[7] P. J. Goulart, E. C. Kerrigan, and J. M. Maciejowski, "Optimization over state feedback policies for robust control with constraints," *Automatica*, vol. 42, no. 4, p. 523–533, 2006.

[8] F. Oldewurtel, C. N. Jones, and M. Morari, "A tractable approximation of chance constrained stochastic MPC based on affine disturbance feedback," in *2008 47th IEEE Conference on Decision and Control*. IEEE, 2008, p. 4731–4736.

[9] P. Hokayem, D. Chatterjee, and J. Lygeros, "On stochastic receding horizon control with bounded control inputs," in *Proceedings of the 48th IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference*, 2009, pp. 6359–6364.

[10] D. Chatterjee, P. Hokayem, and J. Lygeros, "Stochastic receding horizon control with bounded control inputs: A vector space approach," *IEEE Transactions on Automatic Control*, vol. 56, no. 11, p. 2704–2710, 2011.

[11] P. Hokayem, E. Cinquemani, D. Chatterjee, F. Ramponi, and J. Lygeros, "Stochastic receding horizon control with output feedback and bounded controls," *Automatica*, vol. 48, no. 1, p. 77–88, 2012.

[12] J. A. Paulson, E. A. Buehler, R. D. Braatz, and A. Mesbah, "Stochastic model predictive control with joint chance constraints," *International Journal of Control*, vol. 93, no. 1, pp. 126–139, 2020.

[13] L. Hewing and M. N. Zeilinger, "Stochastic model predictive control for linear systems using probabilistic reachable sets," in *2018 57th IEEE Conference on Decision and Control (CDC)*, 2018, pp. 5182–5188.

[14] L. Hewing, K. P. Wabersich, and M. N. Zeilinger, "Recursively feasible stochastic model predictive control using indirect feedback," *Automatica*, vol. 119, 2020, Art. no. 109095.

[15] L. M. Chaouach, M. Fiacchini, and T. Alamo, "Stochastic model predictive control for linear systems affected by correlated disturbances," *IFAC-PapersOnLine*, vol. 55, no. 25, p. 133–138, 2022.

[16] Y. Ao, J. Köhler, M. Prajapat, Y. As, M. N. Zeilinger, P. Fürnstahl, and A. Krause, "Stochastic model predictive control for sub-gaussian noise," 2025, arXiv: 2503.08795.

[17] M. Fiacchini, M. Mammarella, and F. Dabbene, "Recursive feasibility for stochastic MPC and the rationale behind fixing flat tires," 2025. [Online]. Available: https://arxiv.org/abs/2504.17718

[18] A. D. Bonzanini, T. L. Santos, and A. Mesbah, "Tube-based stochastic nonlinear model predictive control: A comparative study on constraint tightening," *IFAC-PapersOnLine*, vol. 52, no. 1, pp. 598–603, 2019.

[19] H. Schlüter and F. Allgöwer, "A constraint-tightening approach to nonlinear stochastic model predictive control under general bounded disturbances," *IFAC-PapersOnLine*, vol. 53, no. 2, pp. 7130–7135, 2020.

[20] E. Joa, M. Bujarbaruah, and F. Borrelli, "Output feedback stochastic MPC with hard input constraints," in *2023 American Control Conference (ACC)*, 2023, pp. 2034–2039.

[21] I. Perez-Salesa, D. V. Dimarogonas, C. Sagúés, and R. Aldana-López, "Cooperative stochastic MPC under hard input constraints and event-triggered communication," *IEEE Transactions on Control of Network Systems*, p. 1–12, 2025.

[22] J. Köhler and M. N. Zeilinger, "Predictive control for nonlinear stochastic systems: Closed-loop guarantees with unbounded noise," *IEEE Transactions on Automatic Control*, pp. 1–16, 2025.

[23] C. Karam, M. Tacchi-Bénard, and M. Fiacchini, "Probabilistic reachable set estimation for saturated systems with unbounded additive disturbances," in *2025 IEEE 64th Conference on Decision and Control (CDC)*. IEEE, Dec. 2025, p. 4547–4553.

[24] E. Kofman, J. A. De Doná, and M. M. Seron, "Probabilistic set invariance and ultimate boundedness," *Automatica*, vol. 48, no. 10, pp. 2670–2676, 2012.

[25] F. Blanchini, S. Miani et al., *Set-theoretic methods in control*. Springer, 2008, vol. 78.

[26] I. Kolmanovsky and E. G. Gilbert, "Theory and computation of disturbance invariant sets for discrete-time linear systems," *Mathematical problems in engineering*, vol. 4, no. 4, pp. 317–367, 1998.

[27] M. Fiacchini, S. Tarbouriech, and C. Prieur, "Quadratic stability for hybrid systems with nested saturations," *IEEE Transactions on Automatic Control*, vol. 57, no. 7, pp. 1832–1838, 2012.

[28] R. T. Rockafellar, *Convex Analysis*. Princeton University Press, 1970.

[29] T. Alamo, A. Cepeda, D. Limon, and E. Camacho, "A new concept of invariance for saturated systems," *Automatica*, vol. 42, no. 9, p. 1515–1521, 2006.

[30] M. Fiacchini and T. Alamo, "Probabilistic reachable and invariant sets for linear systems with correlated disturbance," *Automatica*, vol. 132, 2021, Art. no. 109808.

[31] S. Tarbouriech, G. Garcia, J. M. Gomes Da Silva, and I. Queinnec, *Stability and Stabilization of Linear Systems with Saturating Actuators*. Springer, 2011.

[32] D. Henrion and J.-B. Lasserre, "Convergent relaxations of polynomial matrix inequalities and static output feedback," *IEEE Transactions on Automatic Control*, vol. 51, no. 2, p. 192–202, 2006.

[33] J. Kohler and M. N. Zeilinger, "Recursively feasible stochastic predictive control using an interpolating initial state constraint," *IEEE Control Systems Letters*, vol. 6, p. 2743–2748, 2022.

[34] H. Schlüter and F. Allgöwer, "Stochastic model predictive control using initial state optimization," *IFAC-PapersOnLine*, vol. 55, no. 30, p. 454–459, 2022.

[35] ——, "Stochastic model predictive control using initial state and variance interpolation," in *2023 62nd IEEE Conference on Decision and Control (CDC)*. IEEE, 2023, p. 6700–6706.

[36] S. Diamond and S. Boyd, "CVXPY: A Python-embedded modeling language for convex optimization," *Journal of Machine Learning Research*, vol. 17, no. 83, pp. 1–5, 2016.

[37] A. Agrawal, R. Verschueren, S. Diamond, and S. Boyd, "A rewriting system for convex optimization problems," *Journal of Control and Decision*, vol. 5, no. 1, pp. 42–60, 2018.

[38] P. J. Goulart and Y. Chen, "Clarabel: An interior-point solver for conic programs with quadratic objectives," 2024. [Online]. Available: https://arxiv.org/abs/2405.12762
[39] T. A. N. Heirung, J. A. Paulson, J. O'Leary, and A. Mesbah, "Stochastic model predictive control — how does it work?" *Computers & Chemical Engineering*, vol. 114, pp. 158-170, 2018.

[40] R. D. McAllister and J. B. Rawlings, "Nonlinear stochastic model predictive control: Existence, measurability, and stochastic asymptotic stability," *IEEE Transactions on Automatic Control*, vol. 68, no. 3, p. 1524-1536, 2023.