Safeness Index-Based Economic Model Predictive Control of Stochastic Nonlinear Systems

Process operational safety plays an important role in designing control systems for chemical processes. Motivated by this, in this work, we develop a process Safeness Index-based economic model predictive control system for a broad class of stochastic nonlinear systems with input constraints. A stochastic Lyapunov-based controller is first utilized to characterize a region of the state-space surrounding the origin, starting from which the origin is rendered asymptotically stable in probability. Using this stability region characterization and a process Safeness Index function that characterizes the region in state-space in which it is safe to operate the process, an economic model predictive control system is then developed using Lyapunov-based constraints to ensure economic optimality, as well as process operational safety and closed-loop stability in probability. A chemical process example is used to demonstrate the applicability and effectiveness of the proposed approach.


Introduction
Process operational safety has become crucially important in the chemical industry since the failure of process safety devices/human error often leads to disastrous incidents causing human and capital loss [1].Motivated by this, recently, a new class of economic model predictive control systems (EMPC), in which the cost function penalizes process economics instead of the distances from the steady-state in a general quadratic form, was utilized to account for process operational safety and economic optimality based on a function called the Safeness Index [2,3].These new EMPC methods complement previous efforts on economic model predictive control (e.g., [4][5][6][7]), which were not concerned explicitly with process operational safety.Specifically, in [2], a Safeness Index function that indicates the level of safety of a given state was utilized to characterize a safe operating region and used as a constraint in the EMPC design such that the closed-loop state of a nonlinear process is guaranteed to be driven into the safe operating region in finite time in the presence of sufficiently small bounded disturbances and, if the Safeness Index takes a special form related to a Lyapunov function used in the EMPC design, to never again exit that safe operating region while maximizing the economics of the process.However, in general, the Safeness Index does not have to take this special form and may therefore leave the safe operating region for finite periods of time (this may be acceptable depending on how the notion of a "safe" region of operation is selected; e.g., perhaps a "safe" region of operation means it is safe to operate in for all times, but that if the state is not in that region for short periods of time, there is not an immediate concern).Therefore, with a general form of the Safeness Index, the hard constraint on this function in the EMPC design of [2] with a Safeness Index-based constraint may not be feasible.Due to the potential infeasibility issue caused by the hard constraint, the potential for the state to leave the safe operating region unless the Safeness Index has a specific form and the fact that disturbances may not be sufficiently small to guarantee that the closed-loop state re-enters this safe operating region, the EMPC design with a Safeness Index-based constraint may be limited in terms of its applicability to stochastic nonlinear systems.
On the other hand, MPC and EMPC of stochastic nonlinear systems have received a lot of attention recently (e.g., [8,9]).Uncertainty in the process model may be considered to have a worst-case upper and lower bound, or it may be considered to have unbounded variation and therefore be treated in a probabilistic manner.Since the variation of disturbances is not bounded in a stochastic nonlinear system, the Lyapunov-based economic model predictive control (LEMPC) framework [4] developed for nonlinear systems with small bounded disturbances is unable to guarantee closed-loop stability (i.e., the state of the closed-loop system stays within a well-characterized region of the state-space); instead, probabilistic closed-loop stability results are expected in this case.To that end, in [10], the Markov-chain Monte Carlo technique was used to derive the probabilistic convergence to a near-optimal solution for a constrained stochastic optimization problem.In [9], a Lyapunov-based model predictive control (LMPC) method was proposed for stochastic nonlinear systems to drive the state to a steady-state within an explicitly characterized region of attraction in probability.Recently, the work [11] developed a Lyapunov-based EMPC method for stochastic nonlinear systems by utilizing the probability distribution of the disturbance term to derive closed-loop stability and recursive feasibility results in probability.
In the same direction, this work focuses on the design of Safeness Index-based economic model predictive control systems for a broad class of stochastic nonlinear systems with input constraints.Specifically, under the assumption of the stabilizability of the origin of the stochastic nonlinear system via a stochastic Lyapunov-based control law, a process Safeness Index function and the level sets of multiple Lyapunov functions are first utilized to characterize a safe operating region in state-space, starting from which recursive feasibility and process operational safety are derived in probability for the stochastic nonlinear system under an economic model predictive controller.This economic model predictive control method is then designed that utilizes stochastic Lyapunov-based constraints to achieve economic optimality, as well as feasibility and process operational safety in probability in the well-characterized safe operating region.
The rest of the manuscript is organized as follows: in the Preliminaries, the notation, the class of systems and the stabilizability assumptions are given.In the Main Results, the process Safeness Index and the Safeness Index-based LEMPC are introduced.Subsequently, the Safeness Index-based LEMPC using multiple level sets of Lyapunov functions (to broaden the state-space set for which it is recursively feasible) is developed for the nominal system.Based on this, the corresponding stochastic Safeness Index-based LEMPC and its probabilistic process operational safety and feasibility properties are developed for the nonlinear stochastic system.Finally, a nonlinear chemical process example is used to demonstrate the application of the proposed stochastic Safeness Index-based LEMPC.

Notations
Throughout the paper, we use the notation (Ω, F , P) to denote a probability space.The notation |•| is used to denote the Euclidean norm of a vector, and the notation |•| Q denotes the weighted Euclidean norm of a vector (i.e., |x| Q = x T Qx where Q is a positive definite matrix).x T denotes the transpose of x.R + denotes the set [0, ∞).The notation L f V(x) denotes the standard Lie derivative and it is strictly increasing.The function f (x) is said to be a class C k function if the i-th derivative of f exists and is continuous for all i = 1, 2, ..., k.Consider a stochastic process x(t, w) : [0, ∞) × Ω → R n on (Ω, F , P).For each w ∈ Ω, x(•, w) is a realization or trajectory of the stochastic process, and we abbreviate x(t, w) as x w (t).E(A), P(A), E(A | •) and P(A | •) are the expectation, the probability, the conditional expectation and the conditional probability of the occurrence of the event A, respectively.The hitting time τ X of a set X is the first time that the state trajectory hits the boundary of X.Additionally, we define τ X,T (t) = min{τ X , T, t}, where T is the operation time.

Class of Systems
Consider a class of continuous-time stochastic nonlinear systems described by the following system of stochastic differential equations: where x ∈ R n is the stochastic state vector and u ∈ R m is the input vector.The available control action is defined by .., m}.The disturbance w(t) is a standard q-dimensional independent Wiener process defined on the probability space (Ω, F , P). f (•), g(•), and h(•) are sufficiently smooth vector and matrix functions of dimensions n × 1, n × m and n × q, respectively.It is assumed that the steady-state of the system with w(t) ≡ 0 is (x * s , u * s ) = (0, 0).The initial time t 0 is defined as zero (t 0 = 0).We also assume that h(0) = 0 such that the disturbance term h(x(t))dw(t) of Equation ( 1) vanishes at the origin.Definition 1.Given a C 2 Lyapunov function V : R n → R + , the infinitesimal generator (LV) of the system of Equation ( 1) is defined as follows: We assume that L f V(x), L g V(x) and h(x) T ∂ 2 V ∂x 2 h(x) are locally Lipschitz throughout the work.

Definition 2.
Assuming that the equilibrium of the uncontrolled system dx(t) = f (x(t))dt + h(x(t))dw(t) is at the origin, then the origin is said to be asymptotically stable in probability, if for any > 0, the following conditions hold ( [12]): lim , and the origin of the uncontrolled system is asymptotically stable in probability ( [12]).

Stabilizability Assumptions
We assume there exists a stochastic stabilizing feedback control law u = Φ s (x) ∈ U (e.g., [13,14]) such that the origin of the system of Equation ( 1) can be rendered asymptotically stable in probability for all x ∈ D ⊂ R n , where D is an open neighborhood of the origin, in the sense that there exists a positive definite C 2 stochastic control Lyapunov function V that satisfies the following inequality: where α 1 (•) is a class K function.
Based on the controller Φ s (x), we characterize the set We also choose a level set Ω ρ := {x ∈ φ d | V(x) ≤ ρ} of V(x) inside φ d as the stability region for the system of Equation (1).Therefore, the origin of the system of Equation ( 1) is rendered asymptotically stable via the controller Φ s (x) in probability if x(0) = x 0 ∈ Ω ρ .
In this work, we develop an economic MPC design that takes advantage of the Safeness Index function [2] in its design to achieve probabilistic process operational safety in the following sense: Definition 3. Consider the system of Equation ( 1) with input constraints u ∈ U.If there exists a control law u = Φ ∈ U such that the state trajectories of the system for any initial state x(0) = x 0 ∈ S satisfy x(t) ∈ S, ∀ t ≥ 0 with the probability p, where S is a safe operating region in state-space that excludes the unsafe region D, we say that the control law Φ maintains the process state within a safe operating region S with probability p. Remark 1.In general, the safe operating region S is characterized as a subset of the stability region (because process operation is safe provided that the system is operated within a closed-loop stability region) for the closed-loop system of Equation (1) to account for the additional safety constraints.Therefore, if there exists a control law u = Φ(x) ∈ U that maintains the process state within S with the probability p, it also maintains the process state within the stability region at least with probability p.This implies that the probability of process operational safety of the system of Equation ( 1), which we will discuss in the following sections, also gives a lower bound on probabilistic closed-loop stability.

Main Results
In this section, the process Safeness Index and the optimization problem of Safeness Index-based LEMPC designed for the nominal system of Equation (1) with w(t) ≡ 0 are first presented.Based on that, the Safeness Index-based LEMPC using multiple level sets of Lyapunov functions is developed for the nominal system of Equation (1) to guarantee recursive feasibility and to guarantee that the closed-loop state does not enter an unsafe operating region D. Subsequently, the stochastic Safeness Index-based LEMPC is developed for the system of Equation (1) to account for the disturbances w(t) with unbounded variation.The stochastic safety and feasibility in probability of the closed-loop system of Equation (1) are finally investigated under the sample-and-hold implementation of the proposed stochastic Safeness Index-based LEMPC.

Process Safeness Index
In [2], the Safeness Index function S(x) was developed to indicate the level of safety of a given state, through which process operational safety was integrated with process control system design to account for the process operational safety considerations resulting from multivariable interactions or interactions between units.There are various methods of determining the functional form of S(x), for example by utilizing first-principles process models or using systematic safety analysis tools such as HAZOP and fault tree analysis.
Based on the functional form of S(x), the closed-loop state predictions are required to be maintained within a safe region S (where S(x) is below the threshold on the Safeness Index S TH ) by using the Safeness Index-based constraint within the process control design.Additionally, the safety systems (e.g., the alarm, emergency shut-down and relief systems) can be triggered if the threshold S TH is sufficiently exceeded, which implies that the process operation becomes unsafe and further actions are required.

Safeness Index-Based LEMPC
Safeness Index-based LEMPC optimizes an economic cost function L e (•, •) and maintains the closed-loop state of the nominal system of Equation (1) with w(t) ≡ 0 in a safe operating region by utilizing the Safeness Index function as a hard constraint within the LEMPC design.Specifically, the formulation of the Safeness Index-based LEMPC is as follows: max where x is the predicted state trajectory, ST(∆) is the set of piecewise constant functions with sampling period ∆, τ P is the number of sampling periods of the prediction horizon and is the stabilizing feedback control law designed for the nominal system of Equation ( 1) with w(t) ≡ 0 such that the origin of the system of Equation ( 1) can be rendered asymptotically stable.Under the controller Φ n (x), we first characterize the set where 0 < ρ e < ρ is further designed to make the region Ω ρ a forward invariant set in the presence of sufficiently small bounded disturbances.The constraint of Equation (5e) allows the cost function of Equation (5a) to be maximized while keeping the predicted closed-loop state within . The safety constraint of Equation (5f) is applied to maintain the predictions of the closed-loop state within the safe operating region or x(t k ) is outside of S, the constraint of Equation (5g) is activated to decrease V(x) such that x(t) will move towards the origin within the current sampling period.Remark 2. Since the safe operating region S is not necessarily a forward invariant set based on the formulation of the Safeness Index function, the threshold S TH set on the Safeness Index may define a region that is irregularly shaped, for example the grey region in Figure 1 [2] corresponding to a chemical reactor example similar to the one in the section "Application to a Chemical Process Example" of this manuscript.Therefore, the existence of feasible solutions (i.e., the satisfaction of the constraints of Equation ( 5)) of the Safeness Index-based LEMPC is not guaranteed in S due to the constraint of Equation (5f).Additionally, S(x(t)) may not even decrease under the constraint of Equation (5g) due to the same reason (that S is not an invariant set).Considering the above feasibility issue in the formulation of the Safeness Index-based LEMPC, a new Safeness Index-based LEMPC is developed in the following subsection by using multiple Lyapunov functions to characterize the safe operating region S.  5) for the initial condition (0, 0).

Safeness Index-Based LEMPC Using Multiple Level Sets
The improved Safeness Index-based LEMPC for the nominal system of Equation ( 1) with w(t) ≡ 0 is developed utilizing the level sets of two Lyapunov functions V 1 and V 2 to characterize the safe and unsafe operating regions.Throughout this work, we assume that the shape of the stability regions, D, and their intersection are amenable to the treatment in this work, such as the use of only two Lyapunov functions in the LEMPC design and also the types of overlap of the stability regions described.Specifically, as shown in Figure 2, we define two level sets: from which the origin of the nominal system of Equation ( 1) is rendered asymptotically stable.Ω ρ represents the stability region as it is in the Safeness Index-based LEMPC of Equation ( 5), and U s is designed to exclude the unsafe region D where S(x) > S TH .Therefore, the safe operating region S becomes the union of S 1 := Ω ρ ∩ U s and S 2 := Ω ρ \(S 1 ∪ D) in Figure 2.This new Safeness Index-based LEMPC design is formulated by the following optimization problem: where the notation follows that in Equation (5).Ω ρ e and U s e are again chosen as the level sets inside φ n to make Ω ρ and U s forward invariant sets, respectively.In the optimization problem of Equation ( 6 .The contractive constraint of Equation (6g) is activated to decrease both V 1 and V 2 such that the closed-loop state enters the smaller level sets of V 1 and V 2 (i.e., towards the interior of S 1 ).Therefore, under the Safeness Index-based LEMPC of Equation ( 6), if x(t k ) ∈ S 1 , the constraints of Equations (6e)-(6g) maintain the closed-loop state in S 1 .If x(t k ) ∈ S 2 , the constraints of Equations (6e) and (6g) are applied to maintain the closed-loop state in Ω ρ , under which x(t) will stay in S 2 or enter S 1 in some time.
Remark 3. Based on the Safeness Index function S(x) and its threshold S TH , the level set U s of the Lyapunov function V 2 is chosen to exclude the unsafe region D that is originally in the level set Ω ρ of the Lyapunov function V 1 as shown in Figure 2. Since U s and Ω ρ are both forward invariant sets for the nominal system (or the system with sufficiently small bounded disturbances) of Equation ( 1) under the controller Φ n (x) ∈ U that satisfies Vi + κV i (x) ≤ 0, i = 1, 2, κ > 0, it follows that under the corresponding constraint of Equation (6g), the overlapping region S 1 is also an invariant set.Therefore, the infeasibility problem caused by the Safeness Index constraint of Equation (5f) is solved by introducing the second level set U s into the LEMPC design.For the remaining part of the safe operating region S 2 , the constraints of Equations (6e) and (6g) are utilized to ensure that the closed-loop state stays in Ω ρ all the time, which is similar to closed-loop stability under the traditional LEMPC [4].Since the sampling period ∆ has to be sufficiently small in the sample-and-hold implementation of the Safeness Index-based LEMPC of Equation ( 6), we can utilize a sufficiently small ∆ such that x(t k+1 ) is unable to jump into D within one sampling period if x(t k ) ∈ S 2 .This implies that at the next sampling time, the state x(t k+1 ) either stays in S 2 or enters S 1 via the boundary between S 1 and S 2 .In both cases, it is considered that the state is maintained in the safe operating region according to Definition 3.

Remark 4.
Besides the above development of Safeness Index-based LEMPC using multiple Lyapunov functions, there are also other methods that can guarantee the feasibility of the Safeness Index-based constraint in the LEMPC design.For example, in the optimization problem of Equation ( 5), we can choose a more conservative level set of V(x) (i.e., a small level set inside Ω ρ that excludes D) as the safe operating region.However, if the unsafe region characterized by the Safeness Index function is a set of points inside the stability region and is difficult to exclude by a single level set like U s , we may want to use control Lyapunov barrier functions to design the constraints that account for the unsafe region in state-space [15] and overcome the infeasibility problem.

Stochastic Safeness Index-Based LEMPC
Inspired by the Safeness Index-based LEMPC design of Equation ( 6), the stochastic Safeness Index-based LEMPC design is given by the following optimization problem: where the notation follows that in Equation ( 6) except using ρ, ρ e , s, s e , Φ s (x) and LV to replace ρ , ρ e , s , s e , Φ n (x) and V, respectively.For the system of Equation ( 1) with multiple Lyapunov functions, and U s e are level sets of V 1 and V 2 inside φ d , where 0 < ρ e < ρ and 0 < s e < s.Similar to the LEMPC designs of Equations ( 5) and ( 6), the optimal input trajectory determined by the optimization problem of the stochastic Safeness Index-based LEMPC is denoted by u * (t), which is calculated over the entire prediction horizon t ∈ [t k , t k + τ P ∆).The control action computed for the first sampling period of the prediction horizon u * (t k ) is sent to the actuators to be applied over the sampling period, and the optimization problem of Equation ( 7) is re-solved at the next sampling time.
The constraint of Equation (7e) maintains the predicted state in Ω o ρ e when the current state x(t k ) ∈ Ω o ρ e and the constraint of Equation (7f) maintains the predicted state in U o s e when the current state Since there exists a disturbance w(t) with unbounded variation dw(t) in the system of Equation ( 1), process operational safety (i.e., the closed-loop state is bounded in the safe operating region S) can only be ensured in probability.Therefore, in the following sections, we will establish the probabilities of process operational safety of the system of Equation (1) under the stochastic Safeness Index-based LEMPC of Equation (7).

Sample-And-Hold Implementation
We first investigate the impact of the sample-and-hold implementation of Equation ( 7) on the stability of the closed-loop system of Equation ( 1) following similar arguments to those in [9,11].Specifically, the probabilities of the sets Ω ρ and U s remaining invariant under the sample-and-hold implementation of the Safeness Index-based LEMPC of Equation ( 7) with a sampling period ∆ are given as follows.
Theorem 1.Consider the system of Equation (1) with Ω ρ and U s inside φ d under the control actions u computed by the LEMPC of Equation (7).Let u(t) = u(t k ), ∀t ∈ [t k , t k + ∆).Then, given any probability λ ∈ (0, 1], there exist positive real numbers ρ s < ρ e < ρ and ρ s < s e < s where Ω ρ s and U ρ s are level sets of V 1 and V 2 , respectively, around the origin where LV i , i = 1, 2 are not required to remain negative for the nominal system of Equation ( 1) under the sample-and-hold implementation of u(t), and there also exists a sampling period ∆ * := ∆ * (λ), such that if ∆ ∈ (0, ∆ * ], then: Using the results for standard Brownian motion [16], given any probability λ ∈ (0, 1], there exists a sufficiently small B, s.t.P(A B ) = 1 − λ.For each realization x w (t) with x(0) ∈ Ω ρ ∪ U s and w ∈ A B , there almost surely exists a positive real number , where r < 1/2, according to the local Hölder continuity.Therefore, the probability of the event We first prove that the probabilities of Equations ( 8) and ( 9) hold for the first sampling period.It should be noted that the probabilities of Equations ( 8)-( 10) can be generalized to any sampling period t ∈ [t k , t k + ∆] with the measurement of x(t k ) playing the role of x(0) in Equations ( 8)- (10).
Since V i (x), i = 1, 2 satisfies the local Lipschitz condition, there exist positive real numbers We now prove the probability of Equation ( 10) by using the equation ) for the nominal system of Equation ( 1) based on the definition of the value of LV i in φ d .However, LV i (x(t)) < − only holds in probability for the system in the presence of the disturbances w(t).Based on the local Lipschitz 0)) + κρ s − (which follows from the application of the Lipschitz properties of the components of LV i with ∆ * < ∆ 3 ) and the fact that the probability that LV i (x(t)) < − is as follows: Finally, let ∆ * ≤ min{∆ 1 , ∆ 2 , ∆ 3 }, and the probabilities of Equations ( 8)-( 10) are all satisfied for ∆ ∈ (0, ∆ * ].

Stability in Probability
Based on the results from the above section, the probabilistic process operational safety of the closed-loop system of Equation ( 1) under the Safeness Index-based LEMPC of Equation ( 7) applied in a sample-and-hold fashion is established by the following theorem.
Theorem 2. Consider the system of Equation ( 1) under the stochastic Safeness Index-based LEMPC of Equation ( 7) applied in a sample-and-hold implementation (i.e., u(t) = u(i∆), ∀ i∆ ≤ t < (i + 1)∆, i = 0, 1, 2, ...).Then, given ρ e ∈ (0, ρ), s e ∈ (0, s) and probability λ ∈ (0, 1], there exist a sampling time ∆ ∈ (0, ∆ * (λ)] and probabilities β, β , γ, γ ∈ [0, 1]: such that the following probabilities hold: where S e := S 1e ∪ S 2e is a subset of S that subtracts the risk margins ρ − ρ e and s − s e .The relationship among the sets S 1e := Ω ρ e ∩ U s e and S 2e := S e \S 1e and the unsafe region D are shown in Figure 3. Proof.The proof consists of three parts.We first show that under the Safeness Index-based LEMPC of Equation ( 7), any state trajectory initiated from x(0) ∈ S 1e has the probability defined by Equation ( 13) of staying in S 1 := Ω ρ ∩ U s .However, if x(0) ∈ S 2e , we prove that under the Safeness Index-based LEMPC of Equation ( 7), there exists the probability of Equation ( 14) for the state of the closed-loop system to stay in Ω ρ and with a sufficiently small ∆ to stay in the part of Ω ρ that excludes D. Finally, if x(0) is inside S 1 \S o 1e , we can show that the closed-loop state trajectory reaches the boundary of S 1e first before it leaves S 1 (implying it does not enter D) with the probability of Equation (15).However, if x(0) ∈ Ω ρ \Ω o ρ e and x(0) / ∈ (U s ∪ D) (i.e., the white risk margin around S 2e in Figure 3), we show that it does not enter D, ∀t ∈ [0, ∆) in probability, as well.Additionally, for the sake of simplicity, we denote the probabilities and expectations conditional on the event of A W given in the section "Sample-And-Hold Implementation" as P * (•) and E * (•).
Part 1: To show that Equation ( 13) holds for all x(0) ∈ S 1e , we consider both the case that x(0) ∈ S o 1e and that x(0) ∈ ∂S 1e .The former case is handled by Equations ( 8) and ( 9).Specifically, if x(0) ∈ S o 1e , then both x(0) ∈ Ω o ρ e and also x(0) ∈ U o s e .Then, 8) and ( 9) for 13) holds when x(0) ∈ S o 1e .When x(0) ∈ ∂S 1e , Equation ( 13) is also satisfied.To show this, we first assume x(0) ∈ ∂Ω ρ e and prove that the probability of x(t) staying in Ω ρ within one sampling period conditioned on the event of A W is (1 − β).When x(0) ∈ ∂Ω ρ e , Equation (7g) will be utilized in the LEMPC of Equation (7).Under the constraint of Equation (7g), the optimization problem of Equation ( 7) is solved such that LV 1 is forced to be negative for any x(t k ) ∈ Ω ρ \Ω o ρ e , which implies that Equation ( 10) holds (i.e., LV 1 < − for t ∈ [0, ∆] with the probability of the event A W ). Using Dynkin's formula [17], the following equation can be derived: The following probability is derived using similar arguments as in [11], for all x(0) ∈ ∂Ω ρ e : Bounding Equation (17) with Equation (12a) and taking the complementary events, the following probability is obtained: Using the same steps as performed above, we can prove that ∀x(0) ∈ ∂U s e , the probability of x(t) staying in U s within one sampling period conditioned on the event of A W is as follows: Since the set of initial conditions x(0) ∈ S 1e is the intersection of Ω ρ e and U s e , by combining the probabilities of Equations ( 18) and ( 19) together and using Equation ( 10), the probability of Equation ( 13) is obtained via the definition of conditional probability.
Part 2: If x(0) ∈ S 2e ⊂ Ω ρ e , then either x(0) ∈ Ω o ρ e or x(0) ∈ ∂Ω ρ e .If x(0) ∈ S 2e and Ω o ρ e , then Equation ( 8) holds and P(sup and Equation ( 14) therefore holds.If instead x(0) ∈ S 2e and ∂Ω ρ e , then the results of Part 1 indicate that Equation ( 18) holds.Applying the definition of conditional probability, this also gives that Equation ( 14) holds.Moreover, we show that x(t) is maintained inside the safe operating region S within one sampling period with the probability of Equation ( 14) (i.e., ∀t ∈ [0, ∆), x(t) will not jump into D in probability).It is shown in the section "Sample-And-Hold Implementation" that ∀t ∈ [0, ∆), the change of V i (x) is limited (i.e., |V 1 (x(t)) with a sufficiently small sampling period ∆ * (maybe smaller than the one derived by ∆ * ≤ min{∆ 1 , ∆ 2 , ∆ 3 }).Therefore, if x(0) ∈ S 2e ⊂ Ω ρ e , x(t) cannot move across the entire level set U s and jump into D within a sufficiently small ∆ with the probability (1 − λ).Instead, the closed-loop state at the next sampling time either stays in S 2e or moves into S 1e in probability.If x(t) enters S 1e , the probability of Equation ( 13) will be used to estimate the probability of closed-loop process operational safety thereafter.Because (1 − λ) ≥ (1 − β)(1 − λ), for β, λ ∈ (0, 1], Equation ( 14) establishes the probability of x(t) staying in the safe operating region S within one sampling period ∀x(0) ∈ S 2e .
Part 3: If x(0) ∈ S 1 \S o 1e , we show that it is possible that the closed-loop state trajectory hits the boundary of S o 1e before it hits the boundary of S 1 .If both hitting times τ R n \S 1e (∆) and τ S 1 (∆) are longer than a sampling period ∆, Equation ( 15) is trivially satisfied.However, if one of them or both occur within one sampling period, we show that the probability of Equation ( 15) holds by first showing that the extreme case that x(0) ∈ (Ω ρ \Ω o ρ e ) ∩ (U s \U o s e ) (which are the corners where the risk margins ρ − ρ e and s − s e overlap in Figure 3) satisfies Equation (15).We first show that the probability of the event The event A T indicates that the state of the closed-loop system of Equation (1) reaches the boundary of Ω ρ before it reaches the boundary of Ω ρ e .The probability of Equation ( 20) is determined via Equation ( 17) and the fact that the event {τ , the following probability is derived by bounding Equation (20) by Equation (12c): Using the same steps as performed above, we can prove that ∀x(0) ∈ U s \U o s e , the probabilities similar to Equations ( 20) and ( 21) are derived as follows: where , Equation ( 15)) for the case where 1e is obtained by taking the complementary event of Equations ( 21) and (22b) and using the definition of conditional probability.We now address the other two possibilities for x(0 Consider the case where s e before it leaves Ω ρ \Ω o ρ e , for some t ∈ [0, ∆), then for sure it holds that τ R n \S 1e (∆) < τ S 1 (∆) because the closed-loop state trajectory crosses the boundary of S 1e first.Therefore, Equation ( 15) holds for both cases.The same analysis can be performed for the case where Additionally, since it is demonstrated in Part 2 that the change of V i (x) within one sampling period is limited in probability, it follows that ∀x(0) ∈ Ω ρ \Ω o ρ e and x(0) / ∈ (U s ∪ D), x(t) does not enter D in one sampling period with the probability of 1 − λ, which implies that the closed-loop state either stays in S 2 or moves into S 1 in probability.
Remark 5.The Safeness Index-based LEMPC of Equation ( 7) is unable to ensure process operational safety for the closed-loop system of Equation (1) because of stochastic disturbances with unbounded variation.Additionally, in order to achieve process operational safety with higher probability, we should characterize the safe operating region S well and design large enough risk margins (i.e., ρ − ρ e and s − s e ) to avoid frequent activations of backup safety systems.Specifically, in Theorem 2, it is shown that as ρ e and s e decrease, the probabilities of Equations ( 13)-( 15) become larger, which implies that if we want to improve process operational safety, the Safeness Index-based LEMPC design of Equation ( 7) should be designed with more conservatism (i.e., choosing smaller ρ e and s e ).However, an operating region with smaller ρ e and s e in turn leads to less economic benefits, which is undesired for the Safeness Index-based LEMPC of Equation (7).Therefore, the uncertain process operational safety caused by stochastic disturbances with unbounded variation is essentially a trade-off between economic benefits and probabilistic process operational safety (i.e., in practice, we will choose a conservative operating region to make the process sufficiently safe with respect to the unbounded disturbances, especially considering the other safety systems online and the risks involved, while also optimizing process economics).

Feasibility in Probability
Recursive feasibility for the nominal system of Equation (1) with w(t) ≡ 0 under the Safeness Index-based LEMPC of Equation ( 6) is guaranteed since there always exists a solution (e.g., the Lyapunov-based controller Φ n (x) in sample-and-hold) that satisfies all the constraints of Equation ( 6).Now, consider the system of Equation ( 1) that has disturbance w(t) with unbounded variation.Recursive feasibility under the stochastic Safeness Index-based LEMPC of Equation ( 7) can only be guaranteed in probability over the operation period t ∈ [0, τ N ∆).The probability is established as follows, from which it is shown that the probabilistic bounds on recursive feasibility for the remainder of the entire time of operation decrease as the operation period becomes longer (however, this does not necessarily mean the closed-loop system will not remain recursively feasible because at every sampling time, the remaining time of operation decreases and therefore the probability that the LEMPC will remain recursively feasible for the remaining time of operation increases at the next sampling time if the closed-loop state was maintained within S throughout the prior sampling period).Theorem 3. Consider the system of Equation ( 1) under the stochastic Safeness Index-based LEMPC of Equation (7) applied in a sample-and-hold fashion.Then, if x(0) ∈ S, let V 1 (x(t + i∆)) = ρ i < ρ, V 2 (x(t + i∆)) = s i < s, i = 0, 1, ..., τ N − 1, and let A F represent the event that the optimization problem of Equation ( 7) is solved with the satisfaction of recursive feasibility for time t ∈ [0, τ N ∆).The probability of A F can be calculated as follows: where β i and β i are given as follows: Proof.We can derive the probability of Equation ( 23) following similar arguments to those in [11].
Since the deterministic prediction model of Equation (7b) is used in the stochastic Safeness Index-based LEMPC of Equation ( 7), it follows that there always exists a solution u(t) = Φ s ( x(t q )) ∈ U, ∀t ∈ [t q , t q+1 ), q = k, . . ., k + τ P − 1 that satisfies the constraints of Equation (7d-g) over the prediction horizon provided that x(t k ), t k ≥ 0 is inside the safe operating region S. Therefore, this implies that the probability of recursive feasibility (i.e., Equation ( 23)) is equal to the probability of closed-loop process operational safety over t ∈ [0, τ N ∆), which can be obtained via the recursive application of Equation ( 13) with β i and β i of Equation ( 24) and the definition of conditional probability.Additionally, it should be noted that if x(0) ∈ S 2e , the state is not in ∂U s i in Equation (24).In this case, β i simply takes the value of β , and the probability of Equation ( 23) still holds since it is shown in the proof of Theorem 2 that the state either stays in S 2 or moves into S 1 with the probability of 1 − λ (i.e., Equation (23) gives a conservative result in this case).
Remark 6.In Theorem 3, probabilistic process operational safety and probabilistic recursive feasibility over the operation period t ∈ [0, τ N ∆) are established for the closed-loop system of Equation ( 1) under the Safeness Index-based LEMPC of Equation (7).Due to the disturbance w(t) with unbounded variation, the closed-loop state x(t) may leave S at any sampling step, and thus, closed-loop process operational safety and recursive feasibility of the Safeness Index-based LEMPC of Equation ( 7) can only be derived in a probabilistic manner (i.e., ∀t ∈ [0, τ N ∆), these properties hold with the probability of Equation ( 23)).Since the existence of a feasible control action is only guaranteed in the safe operating region S, backup safety systems should be designed to handle the process if the state exits the safe operating region.Additionally, since the probabilities of Equations ( 13)-( 15) are less than one if ρ e < ρ and s e < s, the probabilities of recursive feasibility and process operational safety for t ∈ [0, τ N ∆) decrease as the operation period τ N ∆ becomes longer.However, it should be noted that this dependence is not unique to the MPC, but to all control designs that try to keep the process state within a specific region in state-space in the presence of stochastic disturbances with unbounded variation (i.e., the probability to keep the closed-loop state in S for all the remaining time of operation goes to zero at t 0 as the process operation time τ N → ∞).

Application to a Chemical Process Example
A chemical process example is used to illustrate the application of the stochastic Safeness Index-based LEMPC of Equation (7) to maintain the closed-loop state within a safe operating region in state-space in probability.Specifically, a well-mixed, non-isothermal continuous stirred tank reactor (CSTR) where an irreversible second-order exothermic reaction takes place is considered.The reaction transforms a reactant A to a product B (A → B).The inlet concentration of A, the inlet temperature and the feed volumetric flow rate of the reactor are C A0 , T 0 and F, respectively.The CSTR is equipped with a heating jacket that supplies/removes heat at a rate Q.The CSTR dynamic model is described by the following material and energy balance equations: where C A is the concentration of reactant A in the reactor, V L is the volume of the reacting liquid in the reactor, T is the temperature of the reactor and Q denotes the heat input rate.The concentration of reactant A in the feed is C A0 .The feed temperature and the volumetric flow rate are T 0 and F, respectively.The reacting liquid has a constant density of ρ L and a heat capacity of C p .∆H, k 0 , E and R represent the enthalpy of reaction, pre-exponential constant, activation energy and ideal gas constant, respectively.Process parameter values are given in Table 1.The disturbance terms dw 1 and dw 2 in Equation ( 25) are independent standard Gaussian white noise with the standard deviations σ 1 = 2.5 × 10 −3 and σ 2 = 0.15, respectively.It is noted that the disturbance terms of Equation (25) vanish at the steady state.
Table 1.Parameter values of the continuous stirred tank reactor (CSTR).The control objective of the stochastic Safeness Index-based LEMPC of Equation ( 7) is to maximize the production rate of B, while maintaining the closed-loop state trajectories in the safe operating region S in probability.The objective function of Equation (7a) is the production rate of B: L e (x, u) = k 0 e −E/RT C 2 A .The Lyapunov functions are designed using the standard quadratic form V i (x) = x T P i x, i = 1, 2, where the positive definite matrices P 1 = 1060 22 22 0.52 and P 2 = 1060 10 10 5 are chosen to characterize the set φ d for the stochastic system of Equation (25).The nonlinear feedback controllers in [13,18] are utilized as Φ n (x) and Φ s (x), respectively.The level sets of the Lyapunov functions V 1 (x) and V 2 (x) are chosen as ρ = 368 and s = 8100 to create a safe operating region S. The explicit Euler method with an integration time step of h c = 10 −4 h is applied to numerically simulate the dynamic model of Equation ( 25).The nonlinear optimization problem of the stochastic Safeness Index-based LEMPC of Equation ( 7) is solved using the IPOPT software package [19] with the sampling period ∆ = 10 −2 h.With the fixed sampling period ∆ = 10 −2 h, ρ = 368 and s = 8100, we focus on the impact of ρ e and s e on probabilistic process operational safety in the following simulations.
It is first shown in Figure 4 that under the Safeness Index-based LEMPC of Equation ( 6) designed for the nominal system of Equation ( 25), the closed-loop state of the nominal system of Equation ( 25) stays in the safe operating region S within the entire operation period t s = 1 h.Additionally, the Safeness Index-based LEMPC of Equation ( 6) is solved successfully in each iteration to obtain a feasible control action u(t) that is applied in the next sampling period.6) for the initial condition (0, 0) (in deviation variable form) with the additional material constraint: 1 It follows that under the stochastic Safeness Index-based LEMPC of Equation ( 7), the state of the closed-loop system of Equation (25) stays in S with different probabilities for different ρ e and s e .To better understand the relationship between probabilistic process operational safety and the choices of ρ e and s e , we derived the experimental probabilities via 500 simulation runs for the same initial condition (∆C As , ∆T s ) = (0 kmol/m 3 , 0 K) and different choices of ρ e and s e (without the material constraint applied for the nominal system).Let A V denote the event that the closed-loop state stays in S over the operation period t s = 1 h.The results are reported in Table 2. From Table 2, it is observed that with fixed s e , P(A V ) becomes larger as ρ e decreases.Likewise, with fixed ρ e , P(A V ) increases as s e decreases.It is demonstrated that a higher probability of closed-loop process operational safety of the system of Equation ( 25) is achieved when ρ e and s e are more conservative.Let ρ e = 320 and s e = 6800.It is obtained that the probability of the states of the closed-loop system of Equation (25) remaining in the safe operating region S reaches 97.4%.Additionally, the averaged total economic benefit (i.e., the time integral of the stage cost L e over the operation period t s = 1 h) is 24.3 under the Safeness Index-based LEMPC of Equation (7), which has an improvement of 81% compared to 13.4 under steady-state operation.Therefore, in this example, the closed-loop system of Equation (25) under the Safeness Index-based LEMPC achieves a relatively high probability of process safety and a satisfactory process economic performance simultaneously with ρ e = 320 and s e = 6800.For an actual process, additional work should likely be performed, which can use techniques like those demonstrated here, to increase the probability of the states of the closed-loop system remaining within the safe operating region to higher values considered acceptable for the process at hand given its design, hazards and the backup measures (alarms/operators, safety systems, relief systems) in place.
On the other hand, it is observed from Table 2 that decreasing ρ e increases the probability P(A V ).By looking at unsafe closed-loop trajectories (i.e., trajectories that leave the safe operating region S under the Safeness Index-based LEMPC of Equation ( 7) during the operation period t s ) in 500 simulation runs (one of them is shown in Figure 5), it is observed that almost all of the unsafe trajectories leave S through the boundary of Ω ρ (i.e., the right edge of Ω ρ in Figure 5).The reason for this behavior is that the local optimum value of L e is calculated to be at the right edge of Ω ρ , which is shown as the yellow region in Figure 6.Therefore, under the Safeness Index-based LEMPC of Equation ( 7), the closed-loop trajectory is optimized to approach this high production rate region and begin circling back due to the disturbances, which leads to a higher probability of leaving the safe operating region S from Ω ρ .Additionally, it is observed in Figure 6 that the production rate decreases as the safe operating region shrinks (i.e., the color becomes darker), which is consistent with the fact that smaller ρ e and s e lead to safer process operation, at the cost of lower economic performance.

Conclusions
In this work, a Safeness Index-based LEMPC design was developed for stochastic nonlinear systems.Under the assumption of stabilizability of the origin of the stochastic nonlinear system via a stochastic Lyapunov-based control law, an economic model predictive controller was developed to account for process operational safety by utilizing Lyapunov-based constraints to maintain the closed-loop state in a safe operating region defined by a Safeness Index function.Under the stochastic Safeness Index-based LEMPC, economic optimality may be achieved with respect to the objective function and sampling period.Additionally, recursive feasibility and process operational safety of the closed-loop stochastic nonlinear system were derived in probability for a well-characterized safe operating region.A chemical reactor example was used to demonstrate the effectiveness of the proposed control method.
Given a set D, we denote the boundary of D by ∂D, the closure of D by D and the interior of D by D o .Set subtraction is denoted by "\", i.e., A\B :

Figure 1 .
Figure 1.A schematic representing the safe operating region S (the gray region) with an example closed-loop trajectory under the Safeness Index-based Lyapunov-based economic model predictive control (LEMPC) design of Equation (5) for the initial condition (0, 0).

Figure 2 .
Figure 2. A schematic representing the unsafe region D (dark gray) and the safe operating region S := S 1 ∪ S 2 (light gray).

Figure 3 .
Figure 3.A schematic representing the unsafe region D (dark gray) and the region S e := S 1e ∪ S 2e (light gray), which is the safe operating region S subtracting the risk margins ρ − ρ e and s − s e .

Figure 4 .
Figure 4. Closed-loop trajectory under the Safeness Index-based LEMPC of Equation (6) for the initial condition (0, 0) (in deviation variable form) with the additional material constraint:1

Figure 5 .
Figure 5.An example closed-loop trajectory under the Safeness Index-based LEMPC of Equation (7) for the initial condition (0, 0) that leaves the safe operating region S, in which ρ e = 320 and s e = 6800.

Figure 6 .
Figure 6.The production rate L e = k 0 e −E/RT C 2 A within the safe operating region S.

Table 2 .
Experimental probability for different values of ρ e and s e .