**1. Introduction**

The buzz around artificial intelligence (AI), machine learning, and data in recent years has sparked both excitement and skepticism from the process systems engineering community [1,2]. Some of the most prevalent uses of data in the process systems field have included its use in developing models of various processes (e.g., Reference [3]) with potential applications in model-based control [4], in learning control laws [5,6], and in process monitoring [7,8]. Control engineers have debated about whether control itself should be considered to be artificial intelligence, particularly as control laws become more advanced. For example, a particularly intelligent form of control (known as economic model predictive control (EMPC) [9–12]) is an optimization-based control strategy that determines the optimal manner in which to operate a chemical process in the sense that the control actions optimize a profit metric for the process over a prediction horizon, subject to process constraints. The significant potential benefits of this control law for next-generation manufacturing have prompted a wide range of investigations in the context of EMPC, including how it may be used for building temperature regulation [13], wastewater treatment [14], microgrid dispatch [15], and gas pipeline networks [16]. Though chemical processes have traditionally been operated at steady-state, EMPC does not necessarily enforce steady-state operation in its efforts to optimize process economic performance. This has raised key questions for this control design regarding

important properties of intelligent systems such as interpretability of its operating strategy and verification that it will work correctly for the real environment that it will need to control and interact with.

Interpretability is a desirable property for artificially intelligent systems. It has been considered in a variety of contexts; for example, the issue of building interpretable data-driven models has been considered to be enhanced by sparse regression, where a model with a small number of available possible terms which could be utilized to build it is derived (with an underlying assumption being that simpler models are more physically realistic and therefore should be more interpretable) [17]. Models identified via sparse regression techniques have been utilized in model predictive control for hydraulic fracturing [18]. Interpretability of other model-building strategies has also been a consideration; for example, for neural networks, where interpretability may be considered to be multidimensional, but to generally constitute whether a human can trace how a neural network obtained its conclusions via how the input information was processed [19], recurrent neural networks with long short-term memory were analyzed for how their cells processed different aspects of character-level language models [20].

It is recognized that interpretability of the control actions computed by an EMPC will be a major determining factor in the adoption of EMPC in the process industries (because, if operators and engineers do not know if the process is in an upset condition, they will likely disable features of the controller that make it difficult to understand due to the need to be sure that safety is maintained at all times). Interpretability for EMPC has not ye<sup>t</sup> received significant focus in the literature. The subset of EMPC formulations which track a steady-state [21] possess a form of interpretability in that the reference behavior is understood by engineers and operators. Reference [22] developed an EMPC formulation in which the desired closed-loop process response specified or restricted by an operator or engineer is tracked by the controller. However, developing the best means for ensuring interpretability for EMPC to appropriately trade off end user understanding with economic optimality remains a largely open question. This work provides new perspectives on this important issue, suggesting that a controller formulation that bridges the human–machine interface by allowing the adjustment of constraints in response to human opinions about the process behavior under the EMPC may provide new avenues of both democratizing advanced control and allowing end users to adjust the response to their liking from an interpretability standpoint.

Another important topic for intelligent control systems is enabling their verification (i.e., certifying that they will perform in practice as intended). Verification can take a significant amount of engineering time and expense, and methods for reducing the time required to validate the controller's performance could reduce the cost of advanced control, could promote operational safety, and could make the controller more straightforward to implement (a lack of ability to verify can prevent an intelligent system from being placed online at all). In the control community, a traditional approach to verification is to design controllers with guaranteed robustness to bounded uncertainty and to use this as a certificate that the controller will be able to maintain closed-loop stability in practice (e.g., References [23–25]). This requires some knowledge of the disturbance characteristics (e.g., upper bounds), which may be difficult to fully determine a priori but is important for EMPC, as the controller could drive the closed-loop state to operate at boundaries of safe operating regions to optimize profits, where the uncertainty in the disturbance characteristics could lead to unsafe conditions. Additional conservatism to account for the uncertainty could lead to over-conservatism that could decrease profits. Other methods for handling disturbances in EMPC have been developed, including methods that account for disturbances probabilistically (making assumptions on their distribution) [26] or adapting models used by the predictive controller online (e.g., References [27–29]). Results on the use of adapting models in EMPC have even included closed-loop stability guarantees when a recurrent neural network that is updated via error triggering is used as the process model [30]. An example of an adaptive control strategy which handles uncertain dynamics in batch processing is that in Reference [31], which uses model predictive control equipped with a probabilistic recursive least squares model parameter update algorithm with a forgetting factor to capture batch process

dynamics. In addition, Reference [32] analyzed a learning-based MPC strategy with a terminal constraint for systems with unmodeled dynamics, where performance is enhanced by using a learned model in the MPC but safety goals are met by ensuring that control actions computed via the MPC are stabilizing.

Another direction that has received attention for handling uncertainty is fault tolerance in the sense of controller reconfigurations upon detection of an actuator fault/anomaly (e.g., Reference [33]) or anomaly response cast in a framework of fault-tolerant control handled via fault/anomaly detection followed by updating the model used by a model-based controller [34]. In Reference [35], fault-tolerant control for nonlinear switched systems was analyzed in the context of safe parking for model predictive control with a steady-state tracking objective function for actuator faults. For EMPC, Reference [36] handled faults through error-triggered data-driven model updates in the controller, and the uniting of EMPC with driving the state into safety-based regions in state-space (e.g., References [37,38]) also constitutes a form of fault-handling. Despite these advances in handling anomalies and uncertainty, which are critical for addressing moving toward a verification paradigm for EMPC, verifying the controller today would still be expected to be time-consuming; additional work is needed to explore further ways of considering and establishing verification for the control design.

Another approach in verification of controllers has been online verification via data-driven models complemented by detection algorithms for problematic controller behavior leading to bounds on the time that would elapse before detection of problematic controller behavior [39]. A feature of this direction in verification, therefore, is the combination of data-driven modeling for control (to address model uncertainty) with guarantees that problematic behavior due to model inaccuracies can be flagged within a given time period. In the present work, we take a conceptually similar approach to verification for EMPC using online anomaly handling with a conservative Lyapunov-based EMPC (LEMPC) [24] design approach that facilitates guaranteed detection of significant plant/model mismatch under sufficient conditions and allows upper bounds on the amount of time available until the mismatch would need to be compensated via model updates without compromising closed-loop stability (as well as the characteristics of the resulting control law after model reidentification required to obtain these theoretical results) to be presented. The development of theoretical guarantees on closed-loop stability with data-driven models that can be updated online in LEMPC has some similarities to References [30,40] but is pursued from a different angle that allows the underlying process dynamics to suddenly change and also allows for more general nonlinear data-driven models to be considered (i.e., we do not restrict the modeling methodology to neural networks as in References [30,40]). It also has similarities to the framework for accounting for faults in LEMPC via model updates in Reference [41] but considers a theoretical treatment of anomaly conditions with data-driven LEMPC, which was not explored in that work.

Motivated by the above considerations, this work focuses on advancing both interpretability and verification for EMPC. These are important considerations for human–machine interaction and can be viewed as different aspects of a "responsive" control design in the sense that the controller is made responsive to changing or unexpected conditions like a human would be. We first address the interpretability concept suggested above in an LEMPC framework in which we elucidate conditions under which an LEMPC could be made responsive to potentially inaccurate metrics reflecting the reactions of end users to the LEMPC's behavior without loss of closed-loop stability. We subsequently move in the direction of addressing verification considerations for LEMPC by developing theoretical guarantees which can be made for the controller in the presence of process dynamics anomalies/changes when potentially adapting data-driven models are used in the controller. We evaluate the conditions under which closed-loop stability may be lost in such circumstances, with exploration of bounds on times before which detection and accommodation of the anomaly could be stabilized to avoid potential plant shutdown. Numerical examples utilizing continuous stirred tank reactors (CSTRs) are presented to illustrate major concepts. Throughout, we highlight cases where the proposed methods could interface with other artificial intelligence techniques (e.g., sentiment analysis or image-based sensing) without compromising closed-loop stability, highlighting the range of intelligent techniques which can be used to enhance next-generation control within an appropriate theoretical framework.

This work is organized as follows: in Section 2, preliminaries are presented. These are followed by the main results in Section 3, which consist of controller formulations and implementation strategies, with demonstration via numerical examples, where (1) the controller constraints can be adjusted online in response to potentially inaccurate stimuli without closed-loop stability being lost (Section 3.1) and (2) the control strategy has characterizable properties in the presence of process anomalies resulting in unanticipated changes in the underlying process dynamics (Section 3.2). Section 4 concludes and provides an outlook on the presented results. Proofs for theoretical results associated with the second control strategy noted above are provided in the Appendix. This manuscript is an extended version of Reference [42].

#### **2. Preliminaries**

#### *2.1. Notation*

The operator |·| denotes the vector Euclidean norm. A function *α* : [0, *a*) → [0, ∞) is in class K if it is continuous, if it strictly increases, and if *α*(0) = 0. The notation <sup>Ω</sup>*ρ* defines a level set of a scalar-valued function *V* (i.e., <sup>Ω</sup>*ρ* := {*x* ∈ *Rn* : *<sup>V</sup>*(*x*) ≤ *ρ*}). The operator / signifies set subtraction (i.e., *A*/*B* := {*x* ∈ *Rn* : *x* ∈ *A*, *x* ∈/ *B*}). *x<sup>T</sup>* represents the transpose of the vector *x*. We define a sampling time with the notation *tk* := *k*Δ, *k* = 0, 1, . . ..

#### *2.2. Class of Systems*

This work considers switched nonlinear systems of the following form:

$$\dot{x}\_{a,i} = f\_i(x\_{a,i}(t), u(t), w\_i(t)) \tag{1}$$

where *xa*,*<sup>i</sup>* ∈ *X* ⊂ *Rn* denotes the state vector, *u* ∈ *U* ⊂ *Rm* denotes the input vector (*u* = [*<sup>u</sup>*1, ... , *um*] *T*), and *wi* ∈ *Wi* ⊂ *Rz* denotes the disturbance vector, where *Wi* := {*wi* ∈ *Rz* : |*wi*| ≤ *θi*, *θi* > <sup>0</sup>}, for *i* = 1, 2, .... In this notation, the *i*th model is used for *t* ∈ [*ts*,*i*, *ts*,*i*+<sup>1</sup>), where *xa*,*<sup>i</sup>*(*ts*,*i*+<sup>1</sup>) = *xa*,*i*+<sup>1</sup>(*ts*,*i*+<sup>1</sup>) and *ts*,<sup>1</sup> = *t*0. The vector function *fi* is assumed to be a locally Lipschitz function of its arguments with *f*1(0, 0, 0) = 0 and *fi*(*xa*,*i*,*s*, *ui*,*s*, 0) = 0 for *i* > 1 (i.e., the steady-state of the updated models when *wi* = 0 is at *xa*,*<sup>i</sup>* = *xa*,*i*,*s*, *u* = *ui*,*<sup>s</sup>*). The system of Equation (1) with *wi* ≡ 0 is known as the nominal system. Synchronous measurement sampling is assumed, with measurements available at every *tk* = *k*Δ, *k* = 0, 1, .... It is noted that *ts*,*i*, *i* = 1, 2, ..., is not required to be an integer multiple of *tk*. We define *<sup>x</sup>*¯*a*,*<sup>i</sup>* = *xa*,*<sup>i</sup>* − *xa*,*i*,*<sup>s</sup>* and *u*¯*i* = *u* − *ui*,*<sup>s</sup>* and define ¯ *fi* as *fi* rewritten to have its origin at *<sup>x</sup>*¯*a*,*<sup>i</sup>* = 0, *u*¯*i* = 0, *wi* = 0. Similarly, we define *Ui* to be the set *U* in deviation variable form from *ui*,*<sup>s</sup>* and *Xi* to be the set *X* in deviation variable form from *xa*,*i*,*s*.

We assume that there exists an explicit stabilizing (Lyapunov-based) control law *hi*(*x*¯*a*,*<sup>i</sup>*) = [*hi*,<sup>1</sup>(*x*¯*a*,*<sup>i</sup>*) ... *hi*,*<sup>m</sup>*(*x*¯*a*,*<sup>i</sup>*)]*<sup>T</sup>* that renders the origin of the nominal system of Equation (1) asymptotically stable in the sense that the following inequalities hold:

$$\alpha\_{1,i}(|\mathcal{R}\_{a,i}|) \le V\_i(\mathcal{R}\_{a,i}) \le \alpha\_{2,i}(|\mathcal{R}\_{a,i}|) \tag{2}$$

$$\frac{\partial V(\pounds\_{4,i})}{\partial \mathbb{1}\_{4,i}} f\_i(\pounds\_{a,i}, h\_i(\pounds\_{a,i}), 0) \le -a\_{3,i}(|\pounds\_{a,i}|) \tag{3}$$

$$\left| \left| \frac{\partial V\_l(\vec{x}\_{d,i})}{\partial \vec{x}\_{d,i}} \right| \right| \leq \alpha\_{4,i} (\left| \mathcal{S}\_{d,i} \right|) \tag{4}$$

$$h\_i(\mathfrak{x}\_{a,i}) \in \mathcal{U}\_i \tag{5}$$

for all *<sup>x</sup>*¯*a*,*<sup>i</sup>* ∈ *Di* ⊆ *Rn* and *i* = 1, 2, ..., where *Di* is an open neighborhood of the origin of ¯ *fi*, and for a positive definite, sufficiently smooth Lyapunov function *Vi*. The functions *<sup>α</sup>*1,*i*, *<sup>α</sup>*2,*i*, *<sup>α</sup>*3,*i*, and *<sup>α</sup>*4,*i* are of class K. A level set of *Vi* denoted by <sup>Ω</sup>*ρi* ⊂ *Di* is referred to as the stability region of the system of Equation (1) under the controller *hi*(*x*¯ *<sup>a</sup>*,*<sup>i</sup>*). We consider that <sup>Ω</sup>*ρi* is selected to be contained within *X*. The Lyapunov-based controller is assumed to be Lipschitz continuous such that the following inequalities hold:

$$|h\_{i,j}(\mathbf{x}) - h\_{i,j}(\mathbf{x'})| \le L\_{h,l}|\mathbf{x} - \mathbf{x'}|\tag{6}$$

for a positive constant *Lh*,*<sup>i</sup>* for all *x*, *x* ∈ <sup>Ω</sup>*ρi*, and *i* = 1, 2, . . ., with *j* = 1, . . . , *m*.

Lipschitz continuity of *fi* and sufficient smoothness of *Vi* provide the following inequalities, for positive constants *Mi*, *Lx*,*i*, *Lw*,*i*, *<sup>L</sup>x*,*i*, and *<sup>L</sup>w*,*i*:

$$|\langle \bar{f}(x, u, w\_l) \rangle \subseteq \mathcal{M}\_l \tag{7}$$

$$|\bar{f}\_l(\mathbf{x}, u, w\_l) - \bar{f}\_l(\mathbf{x'}, u, \mathbf{0})| \le L\_{\mathbf{x}, l} |\mathbf{x} - \mathbf{x'}| + L\_{\mathbf{w}, l} |w\_l| \tag{8}$$

$$\left|\frac{\partial V\_{l}(\mathbf{x})}{\partial \mathbf{x}}\bar{f}\_{l}(\mathbf{x},\boldsymbol{\mu},\boldsymbol{w}\_{l}) - \frac{\partial V\_{l}(\mathbf{x}')}{\partial \mathbf{x}}\bar{f}\_{l}(\mathbf{x}',\boldsymbol{\mu},\mathbf{0})\right| \leq L\_{\mathbf{x},i}'|\mathbf{x}-\mathbf{x}'| + L\_{\mathbf{w},i}'|\boldsymbol{w}\_{i}|\tag{9}$$

for all *x*, *x* ∈ <sup>Ω</sup>*ρi*, *u* ∈ *Ui*, and *wi* ∈ *Wi*.

As this work considers responses to unexpected conditions, we consider that there may be cases in which the nonlinear model of Equation (1) may not be available, though an empirical model with the following form may be available:

$$\alpha\_{b,q}(t) = f\_{NL,q}(x\_{b,q}(t), u(t))\tag{10}$$

where *fNL*,*<sup>q</sup>* is a locally Lipschitz nonlinear vector function in *xb*,*<sup>q</sup>* ∈ R*n* and in the input *u* ∈ R*m* with *fNL*,<sup>1</sup>(0, 0) = 0 and *fNL*,*<sup>q</sup>*(*xb*,*q*,*s*, *uq*,*<sup>s</sup>*) = 0 for *q* > 1 (i.e., the steady-state of the updated models is at *xb*,*<sup>q</sup>* = *xb*,*q*,*s*, *u* = *uq*,*<sup>s</sup>*). Here, *q* = 1, 2, ..., to allow for the possibility that, as the underlying process dynamics change (i.e., the value of *i* increases in Equation (1)), it may be desirable to switch the empirical model used to describe the system. However, we utilize the index *q* instead of *i* for the empirical model to signify that we do not assume that the empirical model must switch with the same frequency as the process dynamics. When the model of Equation (10) does switch, we assume that the switch occurs at a time *ts*,*NL*,*q*+<sup>1</sup> in a manner where *xb*,*<sup>q</sup>*(*ts*,*NL*,*q*+<sup>1</sup>) = *xb*,*q*+<sup>1</sup>(*ts*,*NL*,*q*+<sup>1</sup>). We define *<sup>x</sup>*¯*b*,*<sup>q</sup>* = *xb*,*<sup>q</sup>* − *xb*,*q*,*<sup>s</sup>* and *u* ¯ *q* = *u* − *uq*,*<sup>s</sup>* and define ¯ *fNL*,*<sup>q</sup>* as *fNL*,*q*, rewritten to have its origin at *<sup>x</sup>*¯*b*,*<sup>q</sup>* = 0, *u*¯*q* = 0, as follows:

$$\mathfrak{a}\_{b,q}(t) = \check{f}\_{NL,q}(\mathfrak{x}\_{b,q}(t), \mathfrak{a}\_q(t)) \tag{11}$$

Similarly, we define *Uq* to be the set *U* in deviation variable form from *uq*,*<sup>s</sup>* and *Xq* to be the set *X* in deviation variable form from *xb*,*q*,*s*.

We consider that, for the empirical models in Equation (10), there exists a locally Lipschitz explicit stabilizing controller *hNL*,*<sup>q</sup>*(*x*¯*b*,*<sup>q</sup>*) that can render the origin asymptotically stable in the sense that:

$$
\hat{w}\_{1,q}(|\bar{\mathbf{x}}\_{b,q}|) \le \hat{V}\_q(\bar{\mathbf{x}}\_{b,q}) \le \hat{a}\_{2,q}(|\bar{\mathbf{x}}\_{b,q}|) \tag{12a}
$$

$$\frac{\partial \hat{V}\_{\emptyset}(\bar{x}\_{b,\emptyset})}{\partial \mathfrak{x}\_{b,\emptyset}} \, \, \breve{f}\_{NL,\emptyset}(\mathfrak{x}\_{b,\emptyset}, h\_{NL,\emptyset}(\mathfrak{x}\_{b,\emptyset})) \le -\mathfrak{k}\_{\emptyset,\emptyset}(|\mathfrak{x}\_{b,\emptyset}|) \tag{12b}$$

$$\left|\frac{\partial \hat{V}\_{q}(\mathfrak{x}\_{b,q})}{\partial \mathfrak{x}\_{b,q}}\right| \leq \hat{\mathfrak{a}}\_{4,q}(|\bar{\mathfrak{x}}\_{b,q}|) \tag{12c}$$

$$h\_{NL,q}(\mathfrak{x}\_{b,q}) \in \mathcal{U}\_q \tag{12d}$$

for all *<sup>x</sup>*¯*b*,*<sup>q</sup>* ∈ *DNL*,*<sup>q</sup>* (where *DNL*,*<sup>q</sup>* is a neighborhood of the origin of ¯ *fb*,*<sup>q</sup>* contained in *X*), where *V*ˆ*q* : R*n* → R+ is a sufficiently smooth Lyapunov function, *<sup>α</sup>*<sup>ˆ</sup>*i*,*q*, *i* = 1, 2, 3, 4, are class K functions, and *q* = 1, 2, .... We define <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* ⊂ *DNL*,*<sup>q</sup>* as the stability region of the system of Equation (10) under *hNL*,*<sup>q</sup>* and <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q* as a superset of <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* contained in *DNL*,*<sup>q</sup>* and *X*. Lipschitz continuity of *fNL*,*<sup>q</sup>* and sufficient smoothness of *V*ˆ*q* imply that there exist *ML*,*<sup>q</sup>* > 0 and *LL*,*<sup>q</sup>* > 0 such that

$$|f\_{NL,q}(\mathbf{x},\mu)| \le M\_{L,q} \tag{13a}$$

$$\left| \frac{\partial \hat{V}\_{\eta}(\mathbf{x}\_{1})}{\partial \mathbf{x}} \, f\_{\mathrm{NL},\eta}(\mathbf{x}\_{1}, \mathbf{u}) - \frac{\partial \hat{V}\_{\eta}(\mathbf{x}\_{2})}{\partial \mathbf{x}} \, f\_{\mathrm{NL},\eta}(\mathbf{x}\_{2}, \mathbf{u}) \right| \leq L\_{\mathrm{L},\eta} |\mathbf{x}\_{1} - \mathbf{x}\_{2}| \tag{13b}$$

∀*<sup>x</sup>*, *x*1, *x*2 ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* , *u* ∈ *Uq*, and *q* = 1, 2, . . ..

Furthermore, we define *<sup>x</sup>*¯*a*,*i*,*<sup>q</sup>* = *xa*,*<sup>i</sup>* − *xb*,*q*,*<sup>s</sup>* as the variable representing the deviation of each *xa*,*<sup>i</sup>* from the steady-state of the *q*th empirical model of Equation (10) and ¯ *fi*,*<sup>q</sup>* as the right-hand side of Equation (1) when the model is rewritten in terms of the deviation variables *<sup>x</sup>*¯*a*,*i*,*<sup>q</sup>* and *<sup>u</sup>*¯*q*, as follows:

$$\dot{\mathfrak{x}}\_{a,i,\boldsymbol{\uprho}} = \bar{f}\_{i,\boldsymbol{\uprho}}(\mathfrak{x}\_{a,i,\boldsymbol{\uprho}}(t), \mathfrak{u}\_{\boldsymbol{\uprho}}(t), \boldsymbol{w}\_{i}(t)) \tag{14}$$

We assume that the following holds:

$$|\bar{f}\_{i,q}(\mathbf{x}, \boldsymbol{u}', \boldsymbol{w}) - \bar{f}\_{i,q}(\mathbf{x}', \boldsymbol{u}', \mathbf{0})| \le L\_{\mathbf{x},i,q}|\mathbf{x} - \mathbf{x}'| + L\_{\mathbf{w},i,q}|\boldsymbol{w}|\tag{15}$$

$$\left| \left| \frac{\partial \mathcal{V}\_{\boldsymbol{q}}(\mathbf{x})}{\mathrm{d}\mathbf{x}} \right| \bar{f}\_{i,\boldsymbol{q}}(\mathbf{x}, \boldsymbol{u}', \boldsymbol{w}) - \frac{\partial \mathcal{V}\_{\boldsymbol{q}}(\mathbf{x}')}{\mathrm{d}\mathbf{x}} \bar{f}\_{i,\boldsymbol{q}}(\mathbf{x}', \boldsymbol{u}'', \mathbf{0}) \right| \leq L\_{\mathbf{x},i,\boldsymbol{q}}' |\mathbf{x} - \mathbf{x}'| + L\_{\boldsymbol{w},i,\boldsymbol{q}}' |\boldsymbol{w}| \tag{16}$$

for all *x*, *<sup>x</sup>*, *<sup>u</sup>*, *u* and *w* such that *x* + *xb*,*q*,*<sup>s</sup>* − *xa*,*i*,*<sup>s</sup>* ∈ <sup>Ω</sup>*ρi* , *x* + *xb*,*q*,*<sup>s</sup>* − *xa*,*i*,*<sup>s</sup>* ∈ <sup>Ω</sup>*ρi* , *u* + *uq* ∈ *U*, *u* + *uq* ∈ *U*, and *w* ∈ *Wi*. We define a level set of *V* ˆ *q* contained in <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q* that is also contained in <sup>Ω</sup>*ρi* to be <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*,*<sup>i</sup>* , and *Lx*,*i*,*q*, *Lw*,*i*,*q*, *<sup>L</sup>x*,*i*,*q*, *<sup>L</sup>w*,*i*,*<sup>q</sup>* > 0

#### *2.3. Economic Model Predictive Control*

Economic model predictive control (EMPC) [12] is an optimization-based control design formulated as follows:

$$\min\_{\mathcal{A}\_{i}\in S(\Delta)} \int\_{t\_{k}}^{t\_{k+N}} L\_{\mathfrak{e}}(\mathfrak{X}\_{\mathfrak{a},i}(\mathbf{r}), \mathfrak{a}\_{i}(\mathbf{r})) d\mathbf{r} \tag{17}$$

$$\text{s.t. } \dot{\mathbb{1}}\_{a,i}(t) = f\_i(\mathbb{1}\_{a,i}(t), \vec{u}\_i(t), 0) \tag{18}$$

$$\pounds\_{\mathfrak{a}\_{\mathfrak{i}}}(t\_{\mathbb{k}}) = \mathfrak{x}(t\_{\mathbb{k}}) \tag{19}$$

$$
\overline{u}\_l(t) \in \mathcal{U}\_l, \ \forall \ t \in [t\_k, t\_{k+N}) \tag{20}
$$

$$\mathcal{A}\_{d,i}(t) \in \mathcal{X}\_{i\star} \; \forall \; t \in [t\_{\mathbb{k}\star} t\_{\mathbb{k}+\cdot N}) \tag{21}$$

where *Le*(·, ·) represents the stage cost of the EMPC, which can be a general scalar-valued function that is optimized in Equation (17). The notation *u* ∈ *S*(Δ) signifies that *u* is a piecewise-constant input trajectory with period Δ. The prediction horizon is denoted by *N*. Equation (18) represents the nominal process model, with predicted state *<sup>x</sup>*˜¯*a*,*<sup>i</sup>* for the *i*th model. Equations (20) and (21) represent the input and state constraints, respectively. We denote the optimal solution of an EMPC at *tk* by *<sup>u</sup>*<sup>∗</sup>*p*(*tj*|*tk*), *p* = 1, ... , *m*, *j* = *k*, ... , *k* + *N* − 1, where each *<sup>u</sup>*<sup>∗</sup>*p*(*tj*|*tk*) holds for *t* ∈ [*tj*, *tj*+<sup>1</sup>) within the prediction horizon. *x*(*tk*) in Equation (19) signifies that the state measurement represents the actual system state at *tk* placed in deviation variable form with respect to *<sup>x</sup>*¯*a*,*i*,*s*. Due to the potential switching of the underlying process dynamics before the model in Equation (18) is updated, the measurement may come from a dynamic system different than the *i*th model used in Equation (18).

#### *2.4. Lyapunov-Based Economic Model Predictive Control*

A variety of variations on the general EMPC formulation in Equations (17)–(21) have been developed. One such variation which will receive focus in this paper is Lyapunov-based EMPC (LEMPC) [24], which is formulated as in Equations (17)–(21) but with the following Lyapunov-based constraints added as well:

$$V\_l(\tilde{\mathbf{x}}\_{\mu\_l}(t)) \le \rho\_{\epsilon^j s}, \forall \ t \in [t\_k, t\_{k+N}), \text{ if } t\_k \le t^l \text{ and } V\_l(\mathbf{x}(t\_k)) \le \rho\_{\epsilon j} \tag{22}$$

$$\begin{split} \frac{\partial V\_{i}(\mathbf{x}(t\_{k}))}{\partial \mathbf{x}} \bar{f}\_{i}(\mathbf{x}(t\_{k}), u(t\_{k}), 0) &\leq \frac{\partial V\_{i}(\mathbf{x}(t\_{k}))}{\partial \mathbf{x}} \bar{f}\_{i}(\mathbf{x}(t\_{k}), h\_{i}(\mathbf{x}(t\_{k})), 0), \\ \text{if } t\_{k} > t' \text{ or } V\_{i}(\mathbf{x}(t\_{k})) &> \rho\_{t'i} \end{split} \tag{23}$$

where <sup>Ω</sup>*ρ<sup>e</sup>*,*<sup>i</sup>* ⊂ <sup>Ω</sup>*ρi* is selected such that the closed-loop state is maintained within <sup>Ω</sup>*ρi* over time when the process of Equation (1) is operated under the LEMPC of Equations (17)–(23). *t* is a time after which the constraint of Equation (23) is always applied, regardless of the value of *Vi*(*x*(*tk*)). The activation conditions of the LEMPC with respect to *Vi*(*x*(*tk*)) ensure that the LEMPC can maintain closed-loop stability within <sup>Ω</sup>*ρi* as well as recursive feasibility.

#### *2.5. Lyapunov-Based Economic Model Predictive Control with an Empirical Model*

Several prior works have developed LEMPC formulations including empirical models [43,44] when the model of Equation (1) is either unknown or undesirable for use (e.g., more computationally intensive than an empirical model). They have the following form:

$$\min\_{\mathcal{A}\_q(t)\in\mathcal{S}(\Lambda)} \int\_{t\_k}^{t\_{k+N}} [L\_t(\mathfrak{x}\_{b,q}(\tau), \mathfrak{d}\_q(\tau))] d\tau \tag{24a}$$

$$\text{s.t.} \quad \mathring{\mathbf{x}}\_{b,\emptyset} = \mathring{f}\_{NL,\emptyset}(\ddot{\mathbf{x}}\_{b,\emptyset}(t), \ddot{\mathbf{u}}\_{\emptyset}(t)) \tag{24b}$$

$$\mathbf{x}\_{b,q}(t\_k) = \mathbf{x}(t\_k) \tag{24c}$$

$$\exists \bar{x}\_{b,q}(t) \in X\_{q,} \,\forall t \in [t\_k, t\_{k+N}) \tag{24d}$$

$$m\_q(t) \in \mathcal{U}\_{q\_{\prime}}, \forall \; t \in [t\_{k\prime}, t\_{k+N}) \tag{24e}$$

$$\mathcal{V}\_{q}(\mathfrak{x}\_{b,q}(t)) \le \mathfrak{p}\_{\mathfrak{e},q}, \quad \forall t \in [t\_k, t\_{k+N}) \text{ if } \mathfrak{x}(t\_k) \in \Omega\_{\mathfrak{p}\_{\mathfrak{e},q}} \tag{24f}$$

$$\begin{aligned} \frac{\partial \mathcal{V}\_{\boldsymbol{q}}(\mathbf{x}(t\_{k}))}{\partial \mathbf{x}}(\boldsymbol{f}\_{\text{NL},\boldsymbol{q}}(\mathbf{x}(t\_{k}),\boldsymbol{u}(t\_{k}))) &\leq \\ \frac{\partial \mathcal{V}\_{\boldsymbol{q}}(\mathbf{x}(t\_{k}))}{\partial \mathbf{x}}(\boldsymbol{f}\_{\text{NL},\boldsymbol{q}}(\mathbf{x}(t\_{k}),h\_{\text{NL},\boldsymbol{q}}(\mathbf{x}(t\_{k})))) &\text{ if } \mathbf{x}(t\_{k}) \notin \Omega\_{\boldsymbol{\hat{p}}\_{\text{eq}}} \\ \text{or } \quad t\_{k} \geq t' \end{aligned} \tag{24g}$$

where the notation follows that found in Equations (17)–(23) except that the predictions from the nonlinear empirical model are denoted by *<sup>x</sup>*¯*b*,*<sup>q</sup>* (Equation (24b)) and are initialized from a measurement of the state of the *i*th system of Equation (1) (i.e., from the state measurement of whichever model describes the process dynamics at *tk*). Regardless of which dynamic model describes the underlying process dynamics, the *q*th empirical model along with the state (Equation (24d)) and Lyapunov-based stability constraints corresponding to that model are used.

#### **3. Responsive Economic Model Predictive Control Design**

The next sections present two concepts for moving toward interpretability and verifiability goals for EMPC, cast within a framework of making EMPC more responsive to "unexpected" behavior.

#### *3.1. Automated Control Law Redesign*

In this section, we focus on a case in which the process model used does not change over time (i.e., the *i* = 1 process model in Equation (1) is used for all time) and consider the problem that, despite the pushes toward next-generation manufacturing, many companies that may benefit from automation can have difficulty implementing the appropriate advances if they do not have a knowledgeable control engineer on site due to both a lack of knowledge of advanced control as well as a lack of interpretability of the controller's actions. We present one idea for making an LEMPC easier to work with by giving it a "self-design" capability that allows the controller to update its formulation in a manner that satisfies end-user requirements without requiring understanding of the control laws on the part of the end users. Critically, closed-loop stability and recursive feasibility guarantees are retained. This can be considered to be a case in which the human response to the operating strategy is "unexpected" (in the sense that it is not easily predictable by the control designer), but the controller must have the ability to adjust its control law in response to the human reaction.

The first step toward designing an appropriate controller for this scenario is to recognize that the human response to the process behavior is some function of the pattern observed in the state and input data and that the pattern is dictated by the control formulation. For EMPC, for example, it is dictated by the constraints and objective function (though the process model of Equation (18) also plays a role in determining the response, we consider that the model must represent the process at hand and that therefore it cannot be tuned to impact the state/input behavior). Conceptually, then, the solution to handling the "unexpected" response of the end user of the controller is to learn the mapping between the end user's satisfaction with the response and the constraint/objective function formulation and then to use that mapping to find the constraint/objective function formulation that provides "optimal" satisfaction to the end user.

An open question is how to do this and, in particular, how to do it in a manner that provides theoretical guarantees on feasibility/closed-loop stability. To demonstrate this challenge, consider the LEMPC of Equations (17)–(23). The theoretical results for LEMPC which guarantee closed-loop stability and recursive feasibility under sufficient conditions when no changes occur in the underlying process dynamics rely on the constraints of Equations (22) and (23) being present in the control design [24]. Therefore, ad hoc constraint development in an attempt to optimize end-user "satisfaction" with the process response would not be a means for providing closed-loop stability and recursive feasibility guarantees. Instead, any modification of constraints must take place in a more rigorously defined manner.

One approach would be to develop constraints for EMPC which allow "tuning" of the process response but impact neither closed-loop stability nor feasibility as the tuning parameter in these constraints is adjusted. They thus offer some flexibility to the end user in modifying the response but also ensure that the end user's power to adjust the control law is appropriately restricted for feasibility/stability purposes.

An example of constraints which meet this requirement is the input rate of change constraints added to LEMPC in Reference [45]. In the following section, we will discuss in detail how these constraints may be incorporated within the proposed framework for providing an end user with a restricted flexibility in adjusting the process response without losing theoretical properties of LEMPC.

**Remark 1.** *The question of how the human response may be accurately sensed is outside the scope of the present manuscript. A process example will be provided below in which the end user is assumed to take time to rank his or her "satisfaction" with the process behavior under a number of different controllers to develop a mapping between satisfaction and the tuning parameter of the control law. However, human responses could also be considered to be obtained through other machine learning/artificial intelligence methods, such as sentiment analysis [46].*

**Remark 2.** *Potential benefits of an approach that adjusts the controller's behavior based on the end user's response (rather than assuming that some type of standard metric for evaluating control performance (e.g., settling time, rise time, or overshoot of the steady-state) is able to capture the desired response) are that (1) EMPC may operate processes in a potentially time-varying fashion, meaning that the closed-loop state may not be driven to a steady-state and that the behavior of the process under the EMPC may not be easily predictable a priori (e.g., without running closed-loop simulations). Therefore, determining what metrics to use to state whether performance under EMPC is acceptable or not may not be intuitive or easily generalizable, unlike in the case where steady-state operation is desired. (2) Again, unlike the steady-state case, not all end users of a given EMPC formulation may have the same definition of "good" behavior. Ideally, the "best" behavior is the one computed by the EMPC when it optimizes the process economics over the prediction horizon in whatever manner is necessary to ensure that the constraints are met but profit is maximized. However, an end user may not find this to constitute the "best" behavior due to other considerations that are perhaps difficult or costly to include in the control law (for example, the most profitable input trajectories from the perspective of the profit metric being used in Equation (17) may be expected to lead to more actuator wear than is desirable, which will be the subject of the example below). Therefore, it may be difficult to set a general metric on "good" behavior under EMPC, as the additional considerations defining "goodness" that are not directly included in the control law may vary between processes. (3) The concept of designing a controller that is responsive to unexpected evaluations of its behavior could have broader implications, if appropriately developed, than the initial goal of achieving desired process behavior for a given control law. Ideally, developments in this direction would serve as a springboard for reducing a priori control design efforts while increasing flexibility for next-generation manufacturing such that end users are able to achieve many goals during production that they may conceive over time as being important to their operation but without needing to interface extensively with vendors or even needing to update their software to achieve these updated process responses. The vision is one where modifications for manufacturing could become as flexible and safe through new responsive and intelligent controller formulations as modifications to codes are for computer scientists who do not work with physical processes and therefore can readily test and evaluate new protocols to advance the field quickly.*

#### 3.1.1. LEMPC with Self-Designing Input Rate of Change Constraints

In Reference [45], an LEMPC formulation with input rate of change constraints was designed with the form in Equations (17)–(23) but with the following rate of change constraints added on the inputs:

$$|u\_{\boldsymbol{\theta}}(t\_k) - h\_{1,\boldsymbol{\theta}}(\mathbf{x}(t\_k))| \le e\_{\boldsymbol{\ell}}, \ p = 1, \ldots, m \tag{25}$$

$$|u\_p(t\_j) - h\_{1,p}(\mathbb{F}\_{u,i}(t\_j))| \le \mathbf{e}\_{t\_r}, \ p = 1, \dots, m, \ j = k+1, \dots, k+N-1 \tag{26}$$

where *r* ≥ 0. This formulation is demonstrated in Reference [45] to maintain closed-loop stability and recursive feasibility under sufficient conditions and to cause the following constraints to be met:

$$|u\_p^\*(t\_k|t\_k) - u\_p^\*(t\_{k-1}|t\_{k-1})| \le e\_{\text{desired}}, \forall \, p = 1, \dots, m \tag{27}$$

$$|u\_p^\*(t\_j|t\_k) - u\_p^\*(t\_{j-1}|t\_k)| \le e\_{\text{desired}}, \forall \, p = 1, \dots, m, \, j = k+1, \dots, k+N-1 \tag{28}$$

where desired > 0. The goal of this formulation of LEMPC is to utilize input rate of change constraints to attempt to reduce variations in the inputs between sampling periods that have the potential to cause actuator wear.

However, as noted in Reference [47], despite the intent of the method to prevent actuator wear, there is no explicit relationship between desired or *r* and the amount of actuator wear. Therefore, a control engineer seeking to prevent actuator wear for a given process under the LEMPC of Equations (17)–(23), (25), and (26) might design the value of *r* by performing closed-loop simulations of the process under various values of *r* and then by selecting the one that gives the response that the engineer judges to present a sufficient tradeoff between optimizing economic performance and reducing actuator wear. A company with little control expertise on hand, however, may have difficulties with tuning *r* without vendor assistance. The fact that controllers today cannot readily "fix" their response if engineers who do not have control expertise would like the response to have different characteristics presents a hurdle to the adoption of even simple control laws, let alone the more complex designs which we would like to move into widespread use as part of the next-generation manufacturing paradigm.

These potential negative responses to a lack of on-site control expertise might be prevented by allowing the controller itself to be responsive to end-user preferences. For example, the value of *r* might be designed by allowing a short period of operation under the control law of Equations (17)–(23), (25), and (26) with different values of *r*. The engineers at the plant could then look at time periods in the plant data during which each of the values of *r* were used and could evaluate the performance of the plant through some metric that can be recorded. Then, the value of *r* that is predicted to provide the highest rate of satisfaction (based on some relationship between the value of *r* and the evaluation metrics which can be derived through techniques for fitting appropriate models to the kind of data generated, such as regression or other techniques of machine learning) could be selected for use (and further updated over time through a similar mechanism as necessary).

**Remark 3.** *One could argue that the algorithm by which a control engineer judges whether a given value of r is preferable could be represented mathematically (e.g., as an optimization problem with an objective function representing a tradeoff between penalties on input variation and loss of profit). However, for the reasons noted in Remark 2 above and also with the goal of developing an algorithm which may facilitate interpretability of LEMPC by allowing its control law to be self-adjusted based on how end users feel about the response of the process under the controller, we handle this within the general case of "unexpected" scenarios to which we would like to make EMPC responsive.*

LEMPC with Self-Designing Input Rate of Change Constraints: Theoretical Guarantees

The methodology proposed above incorporates human judgments on the process response for different values of *r* for setting *r* in Equations (17)–(23), (25), and (26). Despite the fact that human judgment is imprecise, the LEMPC formulations of Equations (17)–(23), (25), and (26), by design, maintains closed-loop stability and recursive feasibility under sufficient conditions (proven in Reference [45]) that are unrelated to the value of *r*, demonstrating that the combination of control theory and data-driven models for "unexpected" behavior or human intuition may be possible to achieve with theoretical guarantees.

When the proposed strategy for evaluating *r* online via human responses to different values of the parameter *r* is used, closed-loop stability and feasibility still hold; however, it may not be guaranteed that Equations (27) and (28) hold. Since desired is arbitrary in many respects since it is indirectly tied to actuator wear (primarily though human evaluation), the satisfaction of Equations (27) and (28) may not be significant during the time period that an operator or engineer is evaluating *r*.

There is no guarantee that the proposed method will produce a value of *r* that gives "optimal satisfaction" to the end user. However, this is not considered a limitation of the method, as the end user's satisfaction is subjective and various methods for modeling the relationship between *r* and the end user's satisfaction could be examined if one is found to produce an inadequate result. The value of *r* can also be adjusted further over time if the response after an initial value of *r* is chosen is determined not to be preferable. Reference [45] does guarantee however that, throughout all of the time of operation (both when various values of *r* are tested and when a single value of *r* is selected), closed-loop stability and recursive feasibility can be guaranteed. This is because the value of *r* only impacts whether Equations (27) and (28) are satisfied under the LEMPC of Equations (17)–(23), (25), and (26), and Equations (27) and (28) are only of potential concern for actuator wear and not closed-loop stability or feasibility. Furthermore, because Reference [45] demonstrates that *hi*(*x*˜¯*a*,*<sup>i</sup>*(*tq*)), ∀ *t* ∈ [*tq*, *tq*+<sup>1</sup>), *q* = *k*, ... , *k* + *N* − 1 is a feasible solution to Equations (17)–(23), (25), and (26) at every sampling time regardless of the value of *r* because Equations (25) and (26) can be satisfied by *hi*(*x*˜¯*a*,*<sup>i</sup>*(*tq*)), *t* ∈ [*tq*, *tq*+<sup>1</sup>), *q* = *k*, ... , *k* + *N* − 1 for any *r* ≥ 0, the value of *r* can change between two sampling periods as *r* is being evaluated and recursive feasibility (and therefore closed-loop stability, since closed-loop stability depends on Equations (22) and (23) and not on Equations (25) and (26)) will be maintained. Finally, though when *r* is being evaluated, the process profit or actuator wear level may not be the same as they would be after the value of *r* is selected, this is not expected to pose significant problems for many processes if it is performed over a short period of time. Furthermore, if there are hard process constraints defined by *Xi* that must be met in order to ensure that the product produced during the time when *r* is evaluated can be sold, these can be met even as various values of *r* are tried because *<sup>x</sup>*¯*a*,*<sup>i</sup>*(*t*) ∈ <sup>Ω</sup>*ρi* ⊆ *Xi* according to Reference [45] for any value of *r*. Furthermore, Reference [45] also guarantees that, even as the values of *r* are adjusted, the closed-loop state can be driven to a neighborhood of a steady-state to avoid production volume losses as *r* is adjusted if necessary.

**Remark 4.** *The fact that the above stability analysis holds regardless of the value of r indicates that the accuracy of the method used in obtaining r does not impact closed-loop stability. This is particularly important if the method used in obtaining r involves, for example, performing sentiment analysis of human speech data to determine how well humans like a given value of that parameter. We overcome the limitation of interfacing humans with machines by ensuring that the only parameter of the control law design which is modified in response to the algorithm that carries uncertainty is one which, deterministically, does not impact closed-loop stability.*

**Remark 5.** *Though this section on automated control law redesign has explored only input rate of change constraints, other online redesigns may also be possible in control. For example, in the LEMPC formulation of Equations (17)–(23), the value ρ<sup>e</sup>*,*<sup>i</sup> could be modified over time if an appropriate implementation strategy was developed. Specifically, there exist bounds on ρ<sup>e</sup>*,*<sup>i</sup> given in Reference [24] which are required for closed-loop stability to be maintained for the process of Equation (1) operated under the LEMPC of Equations (17)–(23). Given this, a similar strategy to that presented for the selection of r could be utilized to adjust the value of ρ<sup>e</sup>*,*<sup>i</sup> within its bounds online without impacting closed-loop stability. This holds because a value of ρ<sup>e</sup>*,*<sup>i</sup> between the minimum and maximum at a given time would always be utilized. According to Reference [24], the consequence of this is that, at the next sampling time, <sup>x</sup>*¯*a*,*<sup>i</sup>*(*tk*) ∈ <sup>Ω</sup>*ρi . If <sup>x</sup>*¯*a*,*<sup>i</sup>*(*t*) ∈ <sup>Ω</sup>*ρi at the end of every sampling period for any ρ<sup>e</sup>*,*<sup>i</sup> between its minimum and maximum, <sup>x</sup>*¯*a*,*<sup>i</sup>*(*t*) ∈ <sup>Ω</sup>*ρi at all times. If both r and ρ<sup>e</sup>*,*<sup>i</sup> were to be simultaneously varied, for example, closed-loop*

*stability would again hold, as the value of r does not impact closed-loop stability for the reasons noted above and the value of ρ<sup>e</sup>*,*<sup>i</sup> can vary between its minimum and maximum value as just described without impacting closed-loop stability. Recursive feasibility would also not be impacted. This suggests that it may be possible to design more complex control laws with multiple self-tuning parameters that are simultaneously optimized based on human response to develop control laws that behave in a desirable manner online without posing a safety concern due to loss of closed-loop stability.*

EMPC with Self-Designing Input Rate of Change Constraints: Application to a Chemical Process Example

In this section, we employ a process example that demonstrates the concept of self-designing input rate of change constraints. For simplicity, in this example, we do not employ the Lyapunov-based stability constraints of Equations (22) and (23); therefore, no theoretical stability guarantees can be made for this example. However, this does not present problems for illustrating the core concepts of the method of integrating human responses to operating conditions with EMPC.

The process under consideration is an ethylene oxidation process in a continuous stirred tank reactor (CSTR) from Reference [48] with reaction rates from Reference [49]. The following three reactions are considered to occur in the CSTR:

$$\rm C\_2H\_4 + \frac{1}{2}O\_2 \to C\_2H\_4O \tag{29}$$

$$\rm C\_2H\_4 + 3O\_2 \to 2CO\_2 + 2H\_2O \tag{30}$$

$$2\,\mathrm{C\_2H\_4O} + \,\frac{5}{2}\mathrm{O\_2} \to 2\,\mathrm{CO\_2} + 2\,\mathrm{H\_2O} \tag{31}$$

Mass and energy balances for the reactor, in dimensionless form, are as follows:

$$
\dot{\mathfrak{K}}\_1 = \mathfrak{H}\_1 (1 - \mathfrak{K}\_1 \mathfrak{K}\_4) \tag{32}
$$

$$\dot{\mathfrak{R}}\_2 = \mathfrak{A}\_1(\mathfrak{A}\_2 - \mathfrak{A}\_2 \mathfrak{A}\_4) - A\_1 e^{\gamma\_1/\mathfrak{A}\_4} (\mathfrak{A}\_2 \mathfrak{A}\_4)^{0.5} - A\_2 e^{\gamma\_2/\mathfrak{A}\_4} (\mathfrak{A}\_2 \mathfrak{A}\_4)^{0.25} \tag{33}$$

$$\pounds\_3 = -\-\-\-\-\pounds\_3 \pounds\_4 + A\_1 e^{\circ\_1/\pounds\_4} (\pounds\_2 \pounds\_4)^{0.5} - A\_3 e^{\circ\_1/\pounds\_4} (\pounds\_3 \pounds\_4)^{0.5} \tag{34}$$

$$\mathbb{1}\_{4} = \frac{a\_{1}}{\mathbb{1}\_{1}}(1-\mathbb{1}\_{4}) + \frac{B\_{1}}{\mathbb{1}\_{1}}e^{\varepsilon\_{1}/\beta\_{4}}(\mathbb{1}\_{2}\mathbb{1}\_{4})^{0.5} + \frac{B\_{2}}{\mathbb{1}\_{1}}e^{\varepsilon\_{2}/\beta\_{4}}(\mathbb{1}\_{2}\mathbb{1}\_{4})^{0.25} + \frac{B\_{3}}{\mathbb{1}\_{1}}e^{\varepsilon\_{3}/\beta\_{4}}(\mathbb{1}\_{3}\mathbb{1}\_{4})^{0.5} - \frac{B\_{4}}{\mathbb{1}\_{1}}(\mathbb{1}\_{4}-T\_{\varepsilon})\tag{35}$$

where the process model parameters are listed in Table 1; the state vector components *x*¯1, *x*¯2, *x*¯3, and *x*¯4 (i.e., *x*¯ = [*x*¯1 *x*¯2 *x*¯3 *x*¯4] *T*) are dimensionless quantities corresponding to the gas density, ethylene concentration, ethylene oxide concentration, and temperature in the CSTR, respectively; and the input vector components *u*¯1 and *u*¯2 are dimensionless quantities corresponding to the feed volumetric flow rate and the feed ethylene concentration. The process of Equations (32)–(35) has a steady-state at *x*¯1 = 0.998, *x*¯2 = 0.424, *x*¯3 = 0.032, *x*¯4 = 1.002, *u*¯1 = 0.35, and *u*¯2 = 0.5.

**Table 1.** Parameters for the continuous stirred tank reactor (CSTR) of Equations (32)–(35).


An EMPC is designed to control this process by maximizing the yield of ethylene oxide, which is defined by the following equation over a time interval from the initial time (*<sup>t</sup>*0 = 0) to the final time of operation *tf* :

$$Y(t\_f) = \frac{\int\_0^{t\_f} a\_1(\mathbf{r}) x\_3(\mathbf{r}) x\_4(\mathbf{r}) d\mathbf{r}}{\int\_0^{t\_f} a\_1(\mathbf{r}) a\_2(\mathbf{r}) d\mathbf{r}} \tag{36}$$

However, it is assumed that, in addition to the following bounds on the inputs,

$$0.0704 \le \bar{u}\_1 \le 0.7042\tag{37}$$

$$0.2465 \le \text{\u0} \le \text{\u0} \le 2.4648 \tag{38}$$

there is also a constraint on the total amount of material which can be fed to the CSTR over time:

$$\int\_{0}^{t\_f} \vec{n}\_1(\mathbf{r}) \vec{n}\_2(\mathbf{r}) d\mathbf{r} = 0.175 t\_f \tag{39}$$

As Equation (39) fixes the denominator of Equation (36), the stage cost to be minimized using the EMPC is as follows:

$$L\_{\mathbf{r}}(\mathbf{x}, \boldsymbol{\mu}) = -\boldsymbol{\mu}\_1(t)\boldsymbol{\p}\_3(t)\boldsymbol{\p}\_4(t) \tag{40}$$

To attempt to avoid actuator wear, input rate of change constraints will also be considered. The general form of the EMPC for this example is therefore as follows:

$$\min\_{\mathcal{A}\_1, \mathcal{A}\_2 \in \mathcal{S}(\Lambda)} \int\_{t\_k}^{t\_{k+N\_k}} -\vec{\alpha}\_1(\tau) \mathbb{1}\_3(\tau) \mathbb{1}\_4(\tau) d\tau \tag{41}$$

$$\text{3. s.t. Equations (32)-(35)}\tag{42}$$

$$\mathfrak{x}(t\_k) = \mathfrak{x}(t\_k) \tag{43}$$

$$0.0704 \le \bar{u}\_1(t) \le 0.7042, \forall \ t \in [t\_k, t\_{k+N\_k}) \tag{44}$$

$$\exists \ 0.246 \mathcal{S} \le \pi\_2(t) \le 2.4648, \forall \ t \in [t\_k, t\_{k+N\_k}) \tag{45}$$

$$\frac{1}{r\_v} \int\_{rt\_0}^{t\_k} \vec{u}\_1^\*(\mathbf{r}) \vec{u}\_2^\*(\mathbf{r}) d\mathbf{r} + \frac{1}{l\_v} \int\_{t\_k}^{t\_{k+N\_k}} \vec{u}\_1(\mathbf{r}) \vec{u}\_2(\mathbf{r}) d\mathbf{r} = 0.175 \tag{46}$$

$$|\bar{u}\_p(t\_j) - \bar{u}\_p(t\_{j-1})| \le e, \ j = k, \dots, k + N\_k - 1, \ p = 1, 2\tag{47}$$

In this formulation, no Lyapunov-based stability constraints are employed and no closed-loop stability issues arose in the simulations (i.e., the closed-loop state always remained within a bounded region of state-space). Furthermore, due to the lack of Lyapunov-based stability constraints, the input rate of change constraints of Equations (27) and (28) are enforced directly on input differences (i.e., they have the form of Equations Equations (27) and (28) rather than the form of Equations (25) and (26)). *x*˜¯ represents the predicted value of the process state according to the model of Equation (42). *u*¯ ∗ 1 and *u*¯ ∗ 2 represent the optimal values of *u*¯ 1 and *u*¯ 2 that have been applied in past sampling periods (i.e., *u*¯ ∗ 1 = *u*¯ <sup>1</sup>(*tk*−<sup>1</sup>), and *u*¯ ∗ 2 = *u*¯ <sup>2</sup>(*tk*−<sup>1</sup>)). The values of *u*¯ <sup>1</sup>(*tk*−<sup>1</sup>) and *u*¯ <sup>2</sup>(*tk*−<sup>1</sup>) for *k* = 0 are assumed to be the steady-state values of these inputs. *Nk* is a shrinking prediction horizon in the sense that, at the beginning of every operating period of length *tv* = 46.8, the value of *Nk* is reset to 5 but is then reduced by 1 at each subsequent sampling time of the operating period. This shrinking horizon allows the constraint of Equation (39) to be enforced within every operating period to ensure that, by the end of the time of operation, Equation (39) is met. In Equation (46), *r* signifies the operating periods completed since the beginning of the time of operation (e.g., in the first *tv* time units, *r* = 0 because no operating periods have been completed yet).

We assume that the engineers and operators do not know the value of that they would like to impose in the EMPC of Equations (41)–(46) but plan to determine an appropriate value by assessing the process behavior from the same initial condition under EMPC's with different values of and by selecting a value that they expect will give the optimal tradeoff between economic performance and actuator wear reduction. To represent the process behavior as is varied in these experiments, we performed eight closed-loop simulations of the process of Equations (32)–(35) under the EMPC of Equations (41)–(46) from the same initial condition *<sup>x</sup>*¯*I* = [*x*¯1*I <sup>x</sup>*¯2*I <sup>x</sup>*¯3*I <sup>x</sup>*¯4*I*]*<sup>T</sup>* = [0.997 1.264 0.209 1.004]*<sup>T</sup>* using eight different input-rate-of-change constraint formulations (the simulations were performed both with no input rate of change constraints and with values of 0.01, 0.05, 0.1, 0.3, 0.5, 1, and 3). The simulations lasted for 10 operating periods and used a sampling period of Δ = 9.36, an integration step for the model of Equation (42) (i.e., the model used by the controller) of 10−<sup>4</sup> and an integration step for the model of Equations (32)–(35) (i.e., the model of the plant) of 10−5. The open-source interior point solver Ipopt [50] was used to solve all optimization problems. Figures 1 and 2 show the state and input trajectories for each of the values of chosen. Table 2 shows how the yield varies with the choice of . To express the engineer's or operator's judgment of the relative "goodness" of the response that they see when both profit and input variations are considered, the engineers and operators are considered to have ranked the response for a given on a scale of 1 to 10 as shown in Table 2, with 1 being the worst and 10 being the best.

Figure 3 shows the rankings as a function of as solid blue circles. From this figure, we postulate that a model that may fit this data has the following form:

$$\text{Rankimg} = c\_1 e^{(-c\_2 \omega)} \epsilon^{c\_3} + c\_4 \tag{48}$$

Using the MATLAB function lsqcurvefit, the data from Table 2 for the various values of reported was fit to the function in Equation (48), resulting in *c*1 = 68.8901, *c*2 = 3.8356, *c*3 = 0.8480, and *c*4 = 0.7933. The plot of the function fit to the data is shown as the red curve on Figure 3. A more rigorous method could have been utilized to fit the model and the data (involving, for example, more samples and an evaluation of the deviation of the model from the data), but the present method is sufficient for demonstrating the concepts developed in this work.

The utility of the function in Equation (48) is that it provides a mathematical representation of the model that an engineer or operator is using within his or her mind to determine the best value of to utilize when this engineer or operator is not aware of the model himself or herself. This makes the advanced control design more tractable for the operator or engineer to utilize without advanced control knowledge by fitting the "mind of the human" to a function that can then be utilized in optimizing the control design automatically. To demonstrate this, we determine the "optimal" value of based on the model of Equation (48) by differentiating the equation with respect to and by setting it to 0. This gives an "optimal" value of of *<sup>c</sup>*3/*<sup>c</sup>*2 or 0.22. Simulations were performed for 10 operating periods of the process of Equations (32)–(35) under the EMPC of Equations (41)–(46) with this value of and initialized from *<sup>x</sup>*¯*I*, and the resulting state and input trajectories are shown in Figures 4 and 5. The yield is 8.33%.


**Table 2.** Yield variation with .

**Figure 1.** *x*¯1, *x*¯2, *x*¯3, and *x*¯4 trajectories under economic model predictive controllers (EMPCs) with different values of specified in the legend (the gray trajectory labeled "None" corresponds to no input rate of change constraint applied).

**Figure 2.** *u*¯1 and *u*¯2 trajectories under EMPCs with different values of specified in the legend (the gray trajectory labeled "None" corresponds to no input rate of change constraint applied).

**Figure 3.** Scatter plot reflecting rankings in Table 2 (solid blue circles) and the curve fit using lsqcurvefit (solid red line).

**Figure 4.** State trajectories under EMPC with = 0.22.

**Remark 6.** *The rankings in Table 2 are fabricated to demonstrate the concept that a human judgment could be translated to a modification of an EMPC formulation parameter. They were contrived to display a form to which a reasonable model could be readily fit using lsqcurvefit and, furthermore, are highly simplified (e.g., only a single ranking is provided for each value of rather than an average ranking with additional information such as standard deviation that might be expected if more than one individual was to rank the response). For an actual process,*

*the transformation of human opinion on the response to a function of would therefore be expected to be more complex and to potentially involve statistics-based techniques or other methods for obtaining models from process data; however, an investigation of such methods is outside of the scope of this paper, and therefore, a simplified ranking model was used to demonstrate the concept that a control law parameter might be decided upon by evaluating characteristics of a response where there is a tradeoff between competing operating objectives where at least one of them (in this case, the actuator wear) is more difficult to quantify with a simple model such that the incorporation of human judgment can make the control law design potentially simpler (than if, for example, a detailed actuator wear model was to be developed to allow the controller to more accurately predict the wear itself to then prevent it through a constraint on wear rather than input rate of change).*

**Figure 5.** Input trajectories under EMPC with = 0.22.

#### *3.2. EMPC Response to Unexpected Scenarios via Model Updates*

A second case for which we will explore EMPC designs which are responsive to unexpected events considers these "unexpected" events to be defined by a change in the underlying process dynamics (i.e., the value of *i* increases in Equation (1)). This class of problems covers anomaly responses for EMPC, for which we will adopt the common anomaly-handling strategy (as described in the Introduction section) of updating the process model. Mathematically, we assume that the process model was known with reasonable accuracy before the anomaly (i.e., there is an upper bound on the error between the model used in the LEMPC and the model of Equation (1) with *i* = 1).

We make several points with respect to model updates in this section. First, if the underlying dynamics change, it is possible that the structure of the underlying dynamic model has fundamentally changed. When identifying a new model, it may therefore be preferable to identify the parameters of one with a revised structure; this is a case of seeking to identify a more physics-based model from process data [51]. In keeping with the prior section where the potential was shown for integrating machine learning algorithms known to not be guaranteed to provide accurate data with control, we here highlight that, if machine learning-based sensors (e.g., image-based sensors) are utilized with the process, they may aid in suggesting how to update a process model's structure over time to attempt to keep the structure physically relevant. Because such sensing techniques may not provide correct suggestions, however, a model with a structure suggested by such an algorithm does not need to be automatically implemented in model-based control; instead, engineers could consider multiple models after a machine learning-based algorithm suggests that an anomaly/change in the underlying process model has occurred, where one model to be evaluated is that used until this point and the second is a model that includes any updates implied by the sensing techniques. Subsequently, the prediction accuracy of the two models could be compared, and whichever is most accurate can be considered for use in the LEMPC [52]. Like the methodology in Section 3.1.1, this method limits the ability of any attempts to integrate machine learning (in the sensors) and control from impacting closed-loop stability by using it to complement a rigorous control design approach rather than to dictate it.

Second, at a chemical plant, anomalies may be considered to be either those which pose an immediate hazard to humans and the environment and are considered to require plant shutdown upon detection or those which do not. When the anomaly detected requires plant shutdown, generally the safety system is used to take extreme actions like cutting feeds to shut down the plant as quickly as possible; these generally have a prespecified nature (e.g., closing the feed valve). Anomalies that do not present immediate hazards to humans may either result in sufficiently small plant/model mismatch that the controller is robust against or the plant/model mismatch could cause subsequent control actions to drive the closed-loop state out of the expected region of process operation (at which point, the anomaly may be a hazard). We consider that characterizing conditions under which closed-loop stability is not lost in the second case may constitute steps in moving toward verification of EMPC for the process industries with adaptive model updates in the presence of changing process dynamics.

#### 3.2.1. Automated Response to Anomalies: Formulation and Implementation Strategy

In the next section, we will present theoretical results regarding conditions under which an LEMPC could be conservatively designed to handle anomalies of different types in the sense that closed-loop stability would not be lost upon the occurrence of an anomaly or that impending loss of closed-loop stability could be detected by defining a region <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* (a superset of <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* ) which the closed-loop state should not leave unless the anomaly has been significant and the model used by the LEMPC should be attempted to be reidentified to try to maintain closed-loop stability. If the closed-loop state leaves <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* , however, it has also left <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* , so that the LEMPC of Equation (24) may not be feasible. For this reason, the implementation strategy below suggests that, if the closed-loop state leaves <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* , *hNL*,*<sup>q</sup>* should be applied to the process so that a control law with no feasibility issues is used.

The implementation strategy proposed below relies on the existence of two controllers *hNL*,*<sup>q</sup>* and *hNL*,*q*+1, where *hNL*,*<sup>q</sup>* can stabilize the origin of the nominal closed-loop system of Equation (10) and *hNL*,*q*+<sup>1</sup> can stabilize the origin of the nominal closed-loop system of Equation (10) with respect to the *q* + 1th model. Specifically, before the change in the underlying process dynamics that occurs at *ts*,*i*+<sup>1</sup> is detected at *td*,*q*, the process is operated under the LEMPC with the *q*th empirical model. After the change is detected (in a worst case via the closed-loop state leaving <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* ), a worst-case bound *th*,*<sup>q</sup>* is placed on the time available until the model must be updated at time *tID*,*<sup>q</sup>* to the *q* + 1th empirical model to prevent the closed-loop state from leaving a characterizable operating region.

We consider the following implementation strategy for carrying out the above methodology:

1. At *t*0, the *i* = 1 first-principles model (Equation (1)) describes the dynamics of the process. The *q* = 1 empirical model (Equation (10)) is used to design the LEMPC of Equation (24). An index *ihx* is set to 0. An index *ζ* is set to 0. Go to step 2.

	- (a) If *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*+<sup>1</sup> , operate the process under the LEMPC of Equation (24) with *q* ← *q* + 1 and set *ihx*= 0. Else, apply *hNL*,*q*+<sup>1</sup>(*x*(*tk*)) to the process. Return to step 3. *tk*← *tk*+1.
	- (b) If (*tk*+<sup>1</sup> − *td*,*<sup>q</sup>*) < *th*,*q*, gather online data to develop an improved process model as well as updated functions *V* ˆ *q*+1 and *hNL*,*q*+<sup>1</sup>(*x*) and an updated stability region <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*+<sup>1</sup> around the steady-state of the new empirical model but do not ye<sup>t</sup> update the LEMPC and control the process using the prior LEMPC. Else, if (*tk*+<sup>1</sup> − *td*,*<sup>q</sup>*) ≥ *th*,*q*, set *ihx* = 1 and apply *hNL*,*q*+<sup>1</sup>(*x*(*tk*)). Return to step 3. *tk* ← *tk*+1.
	- (c) Operate the process under the LEMPC of Equation (24) that was used at the prior sampling time. Return to step 3. *tk* ← *tk*+1.

We note that we do not specify the detection method to be used in step 3, but the use of a sufficiently conservative <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* (in a sense to be clarified in the following section) allows a worst-case detection mechanism to be that the closed-loop state exits <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* in step 3. We consider that each *ts*,*i*+<sup>1</sup> and *ts*,*i*+<sup>2</sup> are separated by a sufficient period of time such that no second change in the underlying process dynamics occurs before the first change has resulted in an update in the dynamic model and the closed-loop state is within <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*+<sup>1</sup> .

**Remark 7.** *A significant difference between the proposed procedure and that in References [53,54], which also involves switched systems under LEMPC, is that Reference [53] assumes that the time at which the model is to be switched is known a priori. In handling of anomalies, this cannot be known; therefore, the proposed approach corresponds to LEMPC for switched systems with unknown switching times. We place bounds in the next section on a number of properties of the LEMPC of Equation (24) for this case to demonstrate the manner in which closed-loop stability guarantees depend on, for example, how large the possible changes in the process model could be when they occur. The goal is to provide a perspective on the timeframes available for detecting various anomalies without loss of closed-loop stability, which could aid in verification and self-design studies for EMPC.*

#### 3.2.2. Automated Response to Anomalies: Stability and Feasibility Analysis

According to the implementation strategy above, when an anomaly occurs that changes the underlying process dynamics, one of two things will happen: (1) the model used in Equation (24b) remains the same or (2) the change in the underlying process dynamics is detected and the model used in Equation (24b) is changed within a required timeframe to a new model (i.e., *q* is incremented by one in Equation (10)). In this section, we present the conditions under which closed-loop stability can be maintained in either case. For readability, proofs of theorems presented in this section are available in the Appendix.

We first present several propositions. The first defines the maximum difference between the process model of Equation (1) and that of Equation (10) over time when the two models are initialized from the same state, as long as the states of both systems are kept within a level set of *V* ˆ *q* which is also contained within the stability region around the steady-state for the model of Equation (1) and as long as there is no change in the underlying dynamics. The second sets an upper bound on the difference between the value

of *V* ˆ *q* at any two points in <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* . The third provides the closed-loop stability properties of the closed-loop system of Equation (10) under the controller *hNL*,*q*.

**Proposition 1** ([51])**.** *Consider the systems*

$$
\dot{\vec{x}}\_{a,i,q} = \vec{f}\_{i,q}(\vec{x}\_{a,i,q}(t), \vec{u}\_q(t), w\_i(t)) \tag{49a}
$$

$$\dot{\mathfrak{x}}\_{b,q} = \tilde{f}\_{NL,q}(\mathfrak{x}\_{b,q}(t), \mathfrak{u}\_q(t)) \tag{49b}$$

*with initial states <sup>x</sup>*¯*a*,*i*,*<sup>q</sup>*(*<sup>t</sup>*0) = *<sup>x</sup>*¯*b*,*<sup>q</sup>*(*<sup>t</sup>*0) = *<sup>x</sup>*¯(*<sup>t</sup>*0) *contained within* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*,*<sup>i</sup> , with t*0 = 0*, u*¯*q* ∈ *Uq, and wi* ∈ *Wi. If <sup>x</sup>*¯*a*,*i*,*<sup>q</sup>*(*t*) *and <sup>x</sup>*¯*b*,*<sup>q</sup>*(*t*) *remain within* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*,*<sup>i</sup>for t* ∈ [0, *<sup>T</sup>*]*, then there exists a function fW*,*i*,*<sup>q</sup>*(·) *such that:*

$$|\mathcal{A}\_{a,i,q}(t) - \mathcal{A}\_{b,q}(t)| \le f\_{W,i,\rho}(t) \tag{50}$$

*with:*

$$f\_{\mathcal{W},i,\boldsymbol{\eta}}(t) := \frac{L\_{\mathcal{W},i,\boldsymbol{\eta}}\theta\_{\boldsymbol{i}} + M\_{\text{err},i,\boldsymbol{\eta}}}{L\_{\text{x},i,\boldsymbol{\eta}}} (e^{L\_{\text{x},i,\boldsymbol{\eta}}t} - 1) \tag{51}$$

*where Merr*,*i*,*<sup>q</sup>* > 0 *is defined by:*

$$|\left|\tilde{f}\_{i,\emptyset}(\mathbf{x},\boldsymbol{\mu},\mathbf{0})-\tilde{f}\_{NL,\emptyset}(\mathbf{x},\boldsymbol{\mu})\right|\leq M\_{err,i,\emptyset}\tag{52}$$

*for all x contained in* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*,*<sup>i</sup> and u* ∈ *Uq.*

**Proposition 2** ([24,55])**.** *Consider the Lyapunov function <sup>V</sup>*<sup>ˆ</sup>*q*(·) *of the nominal system of Equation (10) under the controller hNL*,*<sup>q</sup>*(·) *that meets Equation (12). There exists a quadratic function fV*,*<sup>q</sup>*(·) *such that:*

$$
\hat{\mathcal{V}}\_{\boldsymbol{\theta}}(\mathbf{x}) \le \hat{\mathcal{V}}\_{\boldsymbol{\theta}}(\mathbf{x}') + f\_{V,\boldsymbol{\theta}}(|\mathbf{x} - \mathbf{x}'|) \tag{53}
$$

*for all x*, *x*¯ ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf<sup>e</sup>*,*qwith*

$$f\_{V, \emptyset}(s) := \mathbb{A}\_{4, \emptyset}(\mathbb{A}\_{1, \emptyset}^{-1}(\mathbb{A}\_{\emptyset}))s + M\_{\mathbb{P}, \emptyset}s^2 \tag{54}$$

*where Mv*,*<sup>q</sup> is a positive constant.*

**Proposition 3** ([51])**.** *Consider the closed-loop system of Equation (10) under hNL*,*<sup>q</sup>*(*x*¯*b*,*<sup>q</sup>*) *that satisfies the inequalities of Equation (12) in sample-and-hold. Let* Δ > 0*,* <sup>ˆ</sup>*W*,*<sup>q</sup>* > 0*, and <sup>ρ</sup>*<sup>ˆ</sup>*saf <sup>e</sup>*,*q* > *ρ*ˆ*q* > *ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>* > *<sup>ρ</sup>*ˆmin*q* > *ρ*<sup>ˆ</sup>*s*,*<sup>q</sup>* > 0 *satisfy the following:*

$$-\hbar\_{\Im,q}(\hbar\_{2,q}^{-1}(\not p\_{s,q})) + L\_{L,q}M\_{L,q}\Delta \leq -\hbar\_{W,q}/\Delta \tag{55}$$

$$
\hat{\rho}\_{\text{min}\_{\theta}} := \max \{ \hat{\mathcal{V}}\_{\theta} (\bar{\mathbf{x}}\_{b,\emptyset}(t+\Delta)) \quad : \quad \hat{\mathcal{V}}\_{\theta} (\bar{\mathbf{x}}\_{b,\emptyset}(t)) \le \hat{\rho}\_{\theta,\emptyset} \}. \tag{56}
$$

*If <sup>x</sup>*¯*b*,*<sup>q</sup>*(0) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q , then,*

$$
\hat{\mathcal{V}}\_{\emptyset}(\pounds\_{b,\emptyset}(t\_{k+1})) - \hat{\mathcal{V}}\_{\emptyset}(\pounds\_{b,\emptyset}(t\_k)) \le -\pounds\_{W,\emptyset} \tag{57}
$$

*for <sup>x</sup>*¯*b*,*<sup>q</sup>*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*<sup>q</sup>*/Ω*ρ*<sup>ˆ</sup>*s*,*<sup>q</sup> and the state trajectory <sup>x</sup>*¯*b*,*<sup>q</sup>*(*t*) *of the closed-loop system is always bounded in* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q for t* ≥ 0 *and is ultimately bounded in* <sup>Ω</sup>*ρ*ˆmin*q .*

The next proposition bounds the error between the actual process state and a prediction of the process state using an empirical model initialized from the same value of the process state over a period of time in which the underlying process dynamics change, but the empirical model is not updated. This requires overlap in stability regions for the *i*th and *i* + 1th models of Equation (1) and for the *q*th model of Equation (10) within <sup>Ω</sup>*ρ*<sup>ˆ</sup> *q*,*i* while the *q*th model is used. The proof of this proposition is available in Appendix A.

**Proposition 4.** *Consider the following systems:*

$$
\hat{\mathbf{x}}\_{d,i,\emptyset} = \bar{f}\_{i,\emptyset}(\bar{x}\_{d,i,\emptyset}(t), \bar{u}\_{\emptyset}(t), w\_{l}(t)) \tag{58}
$$

$$
\hat{\mathbf{x}}\_{b,\emptyset} = \bar{f}\_{\text{NL},\emptyset}(\bar{\mathbf{x}}\_{b,\emptyset}(t), \bar{u}\_{\emptyset}(t)) \tag{59}
$$

$$\pounds\_{a,i+1,q} = f\_{i+1,q}(\pounds\_{a,i+1,q}(t), \pounds\_q(t), w\_{l+1}(t))\tag{60}$$

*with initial states <sup>x</sup>*¯*a*,*i*,*<sup>q</sup>*(*<sup>t</sup>*0) = *<sup>x</sup>*¯*b*,*<sup>q</sup>*(*<sup>t</sup>*0) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*,*<sup>i</sup> with t*0 = 0*, u*¯*q* ∈ *Uq, wi* ∈ *Wi, and wi*+1 ∈ *Wi*+1*. Also, <sup>x</sup>*¯*a*,*i*,*<sup>q</sup>*(*ts*,*i*+<sup>1</sup>) = *<sup>x</sup>*¯*a*,*i*+1,*q*(*ts*,*i*+<sup>1</sup>)*. If <sup>x</sup>*¯*a*,*i*,*<sup>q</sup>*(*t*)*, <sup>x</sup>*¯*b*,*<sup>q</sup>*(*t*)*, <sup>x</sup>*¯*a*,*i*+1,*q*(*t*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*,*<sup>i</sup> for t* ∈ [0, *t*1] *and*

$$\left| \left| \left| \left| \left| \left. \left. \left. \_{a,i+1,q} (s) , \left. \left. \_{q} (s) \right. \_{i+1} (s) \right. \right) - \left. \left. \_{i,q} (\left. \left. \left. \_{a,i,q} (s) \right. \_{q} (s) , \left. \_{q} (s) \right. \right) \right. \right. \right. \right. \right. \right. \right|\_{1} \right| \leq M\_{\text{change},i,q} \tag{61}$$

*for all <sup>x</sup>*¯*a*,*i*,*q*, *<sup>x</sup>*¯*a*,*i*+1,*<sup>q</sup>* ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*,*<sup>i</sup> , u*¯*q* ∈ *Uq, wi* ∈ *Wi, and wi*+1 ∈ *Wi*+1*, then*

$$|\mathcal{K}\_{a,i,q}(t) - \mathcal{K}\_{b,q}(t)| \le f\_{W,i,q}(t) \tag{62}$$

*where fW*,*i*,*<sup>q</sup>*(*t*) *is defined in Equation (51) for t* ∈ [0, *ts*,*i*+<sup>1</sup>] *and*

$$|\mathbb{E}\_{\mathsf{s}\boldsymbol{\beta}+\mathbf{1},\boldsymbol{\eta}}(t) - \mathbb{E}\_{\mathsf{b},\boldsymbol{\eta}}(t)| \leq f\_{\mathsf{W},\boldsymbol{\beta},\boldsymbol{\eta}}(t\_{\mathsf{s},\boldsymbol{\beta}+1} - t\_{\mathsf{0}}) + (M\_{\text{charge},\boldsymbol{\beta},\boldsymbol{\eta}})(t - t\_{\mathsf{s},\boldsymbol{\beta}+1}) + \frac{L\_{\text{s},\boldsymbol{\beta}}\mathsf{s}\_{\boldsymbol{\theta}} + M\_{\text{air},\boldsymbol{\eta}}}{L\_{\text{s},\boldsymbol{\beta}}} (\boldsymbol{\varepsilon}^{\mathcal{L}\_{\boldsymbol{\xi},\boldsymbol{\beta}}} - \boldsymbol{\varepsilon}^{\mathcal{L}\_{\boldsymbol{\xi},\boldsymbol{\beta}}}{\boldsymbol{\varepsilon}^{\mathcal{L}\_{\boldsymbol{\xi}}}}) \tag{63}$$

*for t* ∈ [*ts*,*i*+1, *<sup>t</sup>*1]*.*

The following theorem provides the conditions under which, when no change in the underlying dynamic model occurs throughout the time of operation and *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* , the LEMPC of Equation (24) designed based on *hNL*,*<sup>q</sup>* and the *q*th empirical model of Equation (10) guarantees that the closed-loop state is maintained within <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* over time and is ultimately bounded in a neighborhood of the origin of the model of Equation (10).

**Theorem 1** ([51])**.** *Consider the closed-loop system of Equation (1) under the LEMPC of Equation (24) based on the controller hNL*,*<sup>q</sup>*(*x*) *that satisfies the inequalities in Equation (12). Let W*,*i*,*<sup>q</sup>* > 0*,* Δ > 0*, N* ≥ 1*, and ρ*ˆ*q* > *ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>* > *<sup>ρ</sup>*ˆmin,*<sup>i</sup>*,*<sup>q</sup>* > *ρ*<sup>ˆ</sup>*s*,*<sup>q</sup>* > 0 *satisfy the following:*

$$-\hbar \varepsilon\_{\vartheta} (\hbar\_{\mathfrak{L}\_{\neq}}^{-1} (\not p\_{\varepsilon,\varnothing})) + \hbar \mathfrak{u}\_{\neq} (\hbar\_{1,\neq}^{-1} (\not p\_{\varnothing})) M\_{\mathfrak{err},i,\varnothing} + L\_{\mathfrak{x},i,\neq}^{\prime} M\_{\mathfrak{l}} \Delta + L\_{\mathfrak{w},i,\neq}^{\prime} \theta\_{\mathfrak{l}} \leq -\epsilon \chi\_{\mathfrak{z},i,\neq} / \Delta \tag{64}$$

$$
\mathfrak{d}\_{\mathfrak{e},\mathfrak{q}} \le \mathfrak{p}\_{\mathfrak{q}} - f\_{V,\mathfrak{q}}(f\_{W,i,\mathfrak{q}}(\Lambda)) \tag{65}
$$

*If x*(0) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q and Proposition 3 is satisfied, then the state trajectory <sup>x</sup>*¯*a*,*i*,*<sup>q</sup>*(*t*) *of the closed-loop system is always bounded in* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q for t* ≥ 0*. Furthermore, if t* > *t and*

$$-\hbar\_{3,q}(\hbar\_{2,q}^{-1}(\not p\_{s,q})) + \hbar\_{4,q}(\not p\_{1,q}^{-1}(\not p\_q))M\_{err,i,q} + L\_{\ge i,j,q}'M\_i\Delta + L\_{w,i,q}'\theta\_i \le -\epsilon\_{W,i,q}/\Delta\tag{66}$$

*then the state trajectory xa*,*<sup>i</sup>*(*t*) *of the closed-loop system is ultimately bounded in* <sup>Ω</sup>*ρ*ˆmin,*<sup>i</sup>*,*<sup>q</sup> and defined as follows:*

$$\hat{\rho}\_{\min,i,\ell\_{\ell}} := \max \{ \hat{V}\_{\emptyset}(\bar{\mathbf{x}}\_{\theta,i,\ell}(t+\Delta)) \mid \hat{V}\_{\emptyset}(\bar{\mathbf{x}}\_{\theta,i,\ell}(t)) \le \hat{\rho}\_{\theta,\ell} \} \tag{67}$$

The prior theorem provided conditions under which the closed-loop state is maintained within <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* in the absence of changes in the dynamic model. In the following theorem, we provide sufficient conditions under which the closed-loop state is maintained in <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* after *ts*,*i*. The proof of this result is presented in Appendix B.

**Theorem 2.** *Consider the closed-loop system of Equation (1) under the LEMPC of Equation (24) with hNL*,*<sup>q</sup> meeting Equation (12), where the conditions of Propositions 3 and 4 hold and where* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q is contained in both* <sup>Ω</sup>*ρi and* <sup>Ω</sup>*ρi*+<sup>1</sup> *. If ts*,*i*+<sup>1</sup> ∈ [*tk*, *tk*+<sup>1</sup>)*, such that, after ts*,*i*+1*, the system of Equation (1) is controlled by the LEMPC of Equation (24), where xa*,*<sup>i</sup>*(*ts*,*i*+<sup>1</sup>) = *xa*,*i*+<sup>1</sup>(*ts*,*i*+<sup>1</sup>) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q , and if the following hold true,*

$$-\hbar \gamma\_{\vartheta} (\hbar\_{2,\rho}^{-1}(\not p\_{\vartheta,\emptyset})) + \hbar \iota\_{\vartheta} (\hbar\_{1,\emptyset}^{-1}(\not p\_{\vartheta})) M\_{\ell\tau r,\flat p,\emptyset} + L\_{x,\flat,\flat q}' M\_{\mathbb{P}} \Delta + L\_{w,\flat,\flat q}' \theta\_{\mathbb{P}} \leq -\epsilon \eta\_{\vartheta,\emptyset,\emptyset} / \Delta \tag{68}$$

$$
\rho\_{e,q} \le \rho\_q - f\_{V,q}(f\_{W,p,q}(\Delta))\tag{69}
$$

*for both p* = *i and p* = *i* + 1*, and*

$$\rho\_{t,q} + f\_{V,\boldsymbol{q}}(f\_{W,i,\boldsymbol{q}}\Delta + (M\_{\text{chang},i})\Delta + \frac{L\_{\text{tr},i,\boldsymbol{q}}\theta\_{i} + M\_{\text{trr},i,\boldsymbol{q}}}{L\_{\text{x},i,\boldsymbol{q}}}(e^{L\_{x,i,\boldsymbol{q}}\Delta} - e^{L\_{x,i,\boldsymbol{q}}t\_{s,i+1}})) \leq \rho\_{\boldsymbol{q}}\tag{70}$$

$$\begin{aligned} &-\hbar\_{3,q}(\mathfrak{h}\_{2,q}^{-1}(\mathfrak{\boldsymbol{\mu}}\_{\ell,q}))+\hbar\_{4,q}(\mathfrak{h}\_{1,q}^{-1}(\mathfrak{\boldsymbol{\mu}}\_{\ell}))M\_{\mathtt{err},i,\boldsymbol{\uprho}}+L\_{\mathtt{x},j,\boldsymbol{\uprho}}'M\_{i}\Delta+L\_{\mathtt{w},j,\boldsymbol{\uprho}}'\theta\_{i}+\hbar\_{4,q}(\mathfrak{h}\_{1,q}^{-1}(\mathfrak{\boldsymbol{\mu}}\_{\ell}))M\_{\mathtt{char},\boldsymbol{\uprho},i,\boldsymbol{\uprho}}+L\_{\mathtt{x},j+1,\boldsymbol{\uprho}}'M\_{i+1}\Delta \\ &+L\_{\mathtt{w},j+1,\boldsymbol{\uprho}}'\theta\_{i+1}\leq -\mathfrak{e}\_{\mathtt{W},i,\boldsymbol{\uprho}}'/\Delta \end{aligned} \tag{71}$$

*then the closed-loop state is bounded in* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*qfor all t* ≥ 0*.*

We highlight that these conditions are conservative and not intended to form the least conservative bounds possible. However, they do help to elucidate some of the factors which impact whether a model used in an LEMPC will need to be reidentified to continue to maintain closed-loop stability when the underlying dynamics change, such as the extent to which the dynamics change. The above theorem indicates that, if <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* is initially chosen in a sufficiently conservative fashion and the empirical model is sufficiently close to the underlying process dynamics before the model change, closed-loop stability may be maintained even after the underlying dynamics change if the model changes are such that the empirical model remains sufficiently close to the new dynamic model after the change. In general, anomalies may occur that could violate the conditions of Theorem 2. The result of this could be that the closed-loop state may leave <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* . In this case, it is helpful to characterize conditions under which changes in the underlying dynamics that could be destabilizing could be detected, triggering a model update and controller redesign for the new dynamic model to stabilize the closed-loop system. Therefore, the following theorem characterizes the length of time that the closed-loop state can remain in <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q* after a change in the underlying process dynamics occurs if the conditions of Theorem 2 are not met. This can be used in determining how quickly a model reidentification algorithm would need to successfully provide a new model for the LEMPC of Equation (24) for closed-loop stability to be maintained as a function of factors such as the extent that the new model deviates from the empirical model used in the LEMPC when the underlying dynamics change, the sampling period, and the conservatism in the selection of *ρ*<sup>ˆ</sup>*q*. The proof of this theorem is presented in Appendix C.

**Theorem 3.** *Consider the closed-loop system of Equation (1) under the LEMPC of Equation (24) with hNL*,*<sup>q</sup> meeting Equation (12) and Proposition 3, where* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q is contained in both* <sup>Ω</sup>*ρi and* <sup>Ω</sup>*ρi*+<sup>1</sup> *. If at t* = *ts*,*i*+1*, where ts*,*i*+<sup>1</sup> ∈ [*tk*, *tk*+<sup>1</sup>)*, such that, after ts*,*i*+1*, the system of Equation (1) is controlled by the LEMPC of Equation (24), where xa*,*<sup>i</sup>*(*ts*,*i*+<sup>1</sup>) = *xa*,*i*+<sup>1</sup>(*ts*,*i*+<sup>1</sup>) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q , then if the following hold true with <sup>ρ</sup>*<sup>ˆ</sup>*saf <sup>e</sup>*,*q* > *ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* > *ρ*ˆ*q* > *ρ*<sup>ˆ</sup>*q*,*e, ρ* ˆ *q*,*<sup>e</sup>* > *<sup>ρ</sup>*ˆmin,*<sup>q</sup>*,*<sup>i</sup>* > *ρ*<sup>ˆ</sup>*s*,*<sup>q</sup>* > 0*, and ρ*<sup>ˆ</sup>*q*,*<sup>e</sup>* > *<sup>ρ</sup>*ˆmin,*<sup>i</sup>*+1,*q* > *ρ*<sup>ˆ</sup>*s*,*<sup>q</sup>* > 0*:*

$$-\hbar\_{3,q}(\hbar\_{2,q}^{-1}(\not p\_{s,q})) + \hbar\_{4,q}(\hbar\_{1,q}^{-1}(\not p\_q))M\_{\text{err},i+1,q} + L\_{\text{x},i+1,q}'M\_{i+1}\Delta + L\_{\text{w},i+1,q}'\theta\_{i+1} \le e\_{W,i+1,q}/\Delta \tag{72}$$

$$\rho\_{\varepsilon,\boldsymbol{q}} + f\_{V,\boldsymbol{q}}(f\_{W,i,\boldsymbol{q}}\Delta + (M\_{\text{charge},i,\boldsymbol{q}})\Delta + \frac{L\_{\text{w},i,\boldsymbol{q}}\theta\_{i} + M\_{\text{err},i,\boldsymbol{q}}}{L\_{\text{x},i,\boldsymbol{q}}}(e^{L\_{x,i,\boldsymbol{q}}\Delta} - e^{L\_{x,i,\boldsymbol{q}}t\_{s,i}}) \leq \rho\_{\text{amp},\boldsymbol{q}}\tag{73}$$

$$
\rho\_{\theta} + f\_{V,\emptyset} (f\_{W,i,\emptyset} \Delta + (M\_{\text{charge},i,\emptyset}) \Delta + \frac{L\_{\text{ $\mu$ ,j $\emptyset}} \theta\_{\text{i}} + M\_{\text{err},j,\emptyset}}{L\_{\text{$ \mu $,j$ \emptyset}}} (e^{L\_{x,i,\emptyset} \Delta} - e^{L\_{x,i,\emptyset} t\_{i,\emptyset} + 1})) \le \rho\_{\text{samp},\emptyset} \tag{74}
$$

$$
\rho\_{\mathcal{A},\emptyset} + f\_{V,\emptyset}(f\_{V,i+1,\emptyset}(\Delta)) \le \rho\_{\mathcal{A}mp,\emptyset} \tag{75}
$$

$$
\rho\_{\Psi} + e\_{\text{IV},i+1,\emptyset} \stackrel{\sim}{\rightharpoonup} \rho\_{\text{start},\emptyset} \tag{76}
$$

*as well as Equations (65)–(67), then if <sup>x</sup>*(*ts*,*i*+<sup>1</sup>) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q and* <sup>Ω</sup>*ρ*ˆmin,*<sup>i</sup>*+1,*q* ⊂ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup> and the change to the model is not detected until a sampling time td*,*<sup>q</sup> with <sup>x</sup>*¯(*td*,*<sup>q</sup>*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*<sup>q</sup>*/Ω*ρ*<sup>ˆ</sup>*<sup>q</sup> (x*¯(*td*,*<sup>q</sup>*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* ⊂ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q ) after which hNL*,*<sup>q</sup> is used to control the system in sample-and-hold, then the number of sampling periods between tID*,*<sup>q</sup> and td*,*<sup>q</sup> within which the model in the LEMPC can be updated to a new model meeting Equation (65) with i replaced by i* + 1 *and q replaced by q* + 1 *without the closed-loop state exiting* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q is given by th*,*<sup>q</sup>* = *floor*((*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q*<sup>−</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>*) *W*,*i*,*<sup>q</sup>* )*, where floor represents the "floor" function that returns the largest integer less than the value of the argument. x*¯(*t*) *refers either to <sup>x</sup>*¯*a*,*i*+1,*q*(*t*) *or <sup>x</sup>*¯*a*,*i*,*<sup>q</sup>*(*t*)*, depending on whether ts*,*i*+<sup>1</sup> *is within the sampling period preceding the closed-loop state exiting* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q .*

The following theorem provides the conditions under which the closed-loop state is maintained within <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q*+1 for all times after *tID*,*<sup>q</sup>* and is driven into <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*+<sup>1</sup> after the model reidentification. The proof of the result is presented in Appendix D.

**Theorem 4.** *If* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q* ⊂ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q*+1 *and if both* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q and* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q*+1 *are contained in* <sup>Ω</sup>*ρi and* <sup>Ω</sup>*ρi*+<sup>1</sup> *, then if hNL*,*q*+<sup>1</sup> *is used to control the system after tID*,*<sup>q</sup> while x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q*+<sup>1</sup>/Ω*ρ*<sup>ˆ</sup>*q*+<sup>1</sup> *with the conditions of Equations (65) and (66) met for the q* + 1*th empirical model for the i* + 1*th dynamic system and the LEMPC of Equation (24) using the q* + 1*th empirical model of Equation (10) is used to control the system for all times after x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*+<sup>1</sup> *, then the closed-loop state is then maintained within* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q*+1 *until it enters* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*+<sup>1</sup> *and is then maintained in* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*+<sup>1</sup> *for all subsequent sampling times.*

**Remark 8.** *From a verification standpoint, the proofs above move toward addressing the question of what may happen if a controller is designed and even tested for certain conditions, but the process dynamics change. It provides a theoretical characterization of conditions under which action would subsequently need to be taken as well as indications of the time available to take the subsequent action. However, the results above may be difficult to utilize directly in developing an online monitoring scheme, as many of the theoretical conditions rely on knowing properties of the current and updated models that would likely not be characterizable or would not be known until after the anomaly occurred. However, these still may aid in gaining an understanding of different possibilities. For example, a conservative stability region* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q suggests that larger anomalies could still be detected and mitigated by a combined detection and reidentification procedure without loss of closed-loop stability. Earlier detection may provide more time for reidentification.*

**Remark 9.** *If there is an indication from detection methods that are not based on the closed-loop state leaving the stability region that the underlying dynamics may have changed but that the closed-loop state has not yet left* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q , then until the closed-loop state leaves* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q , online experiments (e.g., modifying the objective function as in Reference [51]) could be performed if they do not impact the constraint set to attempt to probe whether the dynamics are more consistent with the prior process model or the potential model postulated after the anomaly is suggested. This may be a method for attempting to detect the changes before the closed-loop state leaves* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q , which could allow larger changes in the process model to be handled practically than could be guaranteed to be handled in the theorems above, as the magnitude of the deviations in the dynamic model allowed above without loss of closed-loop stability depends on the distance between* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q and* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup> . However, it is also highlighted that the above is a conservative result, meaning that, in general, larger changes may be able to be handled without loss of closed-loop stability.*

**Remark 10.** *The above results can be used to comment on why giving greater flexibility to the process after an anomaly to handle it could introduce additional complexity. Specifically, consider the possibility that some actuators may not typically be used for control but could be considered for use after an anomaly (similar to how safety systems activate for chemical processes, but in this case, they would not act according to a prespecified logic but might be able to be manipulated in either an on-off or continuous manner to give the process additional capabilities for handling the anomaly). It is noted that this would constitute dynamics not previously considered. According to the proofs above, one way to guarantee closed-loop stability in the presence of sufficiently small disturbances is to cause the dynamics after they change to not differ too radically from those assumed before the change and used in the prior dynamic model in the EMPC. If additional flexibility is given to the system, this would be an additional model that would have to match up well.*

**Remark 11.** *The results above suggest that, if a model identification algorithm could be guaranteed to provide an accurate model with a small amount of data that could be gathered between when the closed-loop state leaves* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q but before it leaves* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q (where the amount of data available in that timeframe could be known a priori by the number of measurements available in a given sampling period), then the model could be reidentified and placed within the LEMPC in a manner that is stabilizing.*

**Remark 12.** *Instead of changes to the underlying dynamic model, anomalies may present changes in the constraint set (e.g., anomalies may change equipment material limitations (e.g., maximum shear stresses, which can change with temperature) used to place constraints on the state in an LEMPC). Because the above results assume that the stability region is fully contained within the state constraint set, the detection and response procedure above would need to ensure that there is no time at which the stability region is no longer fully included within the state constraint set under the new dynamic model. This may be handled by making* <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q sufficiently conservative such that the closed-loop state never exits a region where the state constraints can be met under different dynamic models.*

3.2.3. Automated Response to Unexpected Hazards: Application to a Chemical Process Example

In this section, we demonstrate concepts described above through a process example. This example considers a nonisothermal reactor in which an *A* → *B* reaction takes place, but the reactant inlet concentration *CA*0 and the heat rate *Q* supplied by a jacket are adjusted by an LEMPC. The process model is as follows:

$$\mathcal{C}\_A = \frac{\mathbb{E}}{\mathbb{P}} (\mathbb{C}\_{A0} - \mathbb{C}\_A) - k\_0 e^{-\frac{\mathbb{E}}{\mathbb{E}\_\ell^\mathrm{Tr}}} \mathbb{C}\_A^2 \tag{77}$$

$$\dot{T} = \frac{F}{\nabla}(T\_0 - T) - \frac{\Lambda H l\_0}{\rho\_L C\_p} e^{-\frac{h}{R\_d T}} C\_A^2 + \frac{Q}{\rho\_L C\_p V} \tag{78}$$

where the parameters are listed in Table 3 and include the reactor volume *V*, inlet reactant temperature *T*0, pre-exponential constant *k*0, solution heat capacity *Cp*, solution density *ρ<sup>L</sup>*, feed/outlet volumetric flow rate *F*, gas constant *Rg*, activation energy *E*, and heat of reaction Δ*H*. The state variables are the reactant concentration *CA* and temperature *T* in the reactor, which can be written in deviation form from the operating steady-state vector *CAs* = 1.22 kmol/m3, *Ts* = 438.2 K, *CA*<sup>0</sup>*s* = 4 kmol/m3, and *Qs* = 0 kJ/h as *x* = [*x*1 *<sup>x</sup>*2]*<sup>T</sup>* = [*CA* − *CAs T* − *Ts*]*<sup>T</sup>* and *u* = [*<sup>u</sup>*1 *<sup>u</sup>*2]*<sup>T</sup>* = [*CA*0 − *CA*<sup>0</sup>*s Q* − *Qs*]*<sup>T</sup>*. The model of Equations (77) and (78) has the following form:

$$\dot{x} = f(x) + g(x)u \tag{79}$$

where ˜ *f* represents a vector function derived from Equations (77) and (78) that is not multiplied by *u* and where *g*(*x*)=[*g*1 *<sup>g</sup>*2]*<sup>T</sup>* = [ *FV* 0; 0 1 *ρLCpV* ]*T* represents the vector function which multiplies *u* in these equations.


**Table 3.** Parameters for the CSTR model of Equations (77) and (78).

The EMPC utilized to adjust the manipulated inputs *CA*0 and *Q* utilizes the following stage cost (to maximize the production rate of the desired product) and physical bounds on the inputs:

$$L\_e = -k\_0 e^{-E/(R\_\mathcal{J}T(r))} \mathcal{C}\_A(r)^2 \tag{80}$$

$$0.5 \le C\_{A0} \le \mathcal{T}. 5 \text{ kmol/}m^3 \tag{81}$$

$$- \mathbb{S} \times 10^5 \le Q \le \mathbb{S} \times 10^5 \text{ kJ/h} \tag{82}$$

Lyapunov-based stability constraints are also enforced (where a constraint of the form of Equation (22) is enforced at the end of every sampling time if *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*e*, and the constraint of the form of Equation (23) is enforced at *tk* when *x*(*tk*) ∈ <sup>Ω</sup>*ρ*ˆ/Ω*ρ*<sup>ˆ</sup>*e* but then followed by a constraint of the form of Equation (22) at the end of all sampling periods after the first).

We will consider several simulations to demonstrate the developments above. In the first, we explore several aspects of the case in which there is a change in the underlying dynamics while the process is operated under LEMPC that is minor such that the closed-loop state does not leave <sup>Ω</sup>*ρ*<sup>ˆ</sup> after the change in the underlying dynamics. For this case, the Lyapunov function selected was *V*ˆ*q* = *<sup>x</sup>TPx*, with *P* given as follows:

$$P = \begin{bmatrix} 1200 & 5\\ 5 & 0.1 \end{bmatrix} \tag{83}$$

The Lyapunov-based controller *hNL*,<sup>1</sup>(*x*) was designed such that its first component *hNL*,1,1(*x*) = 0 kmol/m<sup>3</sup> and its second component *hNL*,1,2(*x*) is computed as follows (Sontag's formula [56]):

$$h\_{NL,1,2}(x) = \begin{cases} -\frac{L\_{\hat{f}}\hat{V}\_{\hat{q}} + \sqrt{L\_{\hat{f}}\hat{V}\_{\hat{q}}^2 + L\_{\hat{g}\_2}\hat{V}\_{\hat{q}}^4}}{L\_{\hat{g}\_2}\hat{V}\_{\hat{q}}}, & \text{if } L\_{\hat{g}\_2}\hat{V}\_{\hat{q}} \neq 0\\ 0, & \text{if } L\_{\hat{g}\_2}\hat{V}\_{\hat{q}} = 0 \end{cases} \tag{84}$$

Then, it is saturated at the input bounds of Equation (82) if they are met. *L* ˜ *f V*ˆ *q* and *Lg*˜2*V*<sup>ˆ</sup> *q* are Lie derivatives of *V*ˆ *q* with respect to the vector functions ˜ *f* and *g*˜2, respectively. *ρ*ˆ and *ρ*ˆ*e* were taken from Reference [57] to be 300 and 225, respectively. The process state was initialized at *xinit* = [ −0.4 kmol/m<sup>3</sup> 8 K] *T*, with controller parameters *N* = 10 and Δ = 0.01 h. The process model of Equations (77) and (78) was integrated with the explicit Euler numerical integration method using an integration step size of 10−<sup>4</sup> h within the LEMPC and of 10−<sup>5</sup> h to simulate the process.

For this first simulation, we assume that a change in the underlying process dynamics occurs at 0.5 h that does not compromise closed-loop stability. Specifically, at 0.5 h, it is assumed that an additional source of heat arises outside the reactor such that the right-hand side of Equation (78) is modified by the addition of another term *Qextra* = 300 K/h. Figures 6 and 7 show the process responses when the LEMPC is not aware of the change in the process dynamic model when it occurs and when it is aware of the change in the process dynamic model after it occurs such that it is fully compensated (i.e., an accurate process model is used in the LEMPC at all times, even after the dynamics change). In both cases, the closed-loop state was maintained within the stability region at all times. These simulations were carried out in MATLAB R2016b using fmincon with the default settings except for the increased iterations/function evaluations allowed, scaling *u*2 down by 10<sup>5</sup> and providing the steady-state input values as the initial guess for the optimization problem solution at each sampling time. No attempt was made to check whether the LEMPCs in the simulations located globally optimal solutions to the LEMPC optimization problems. However, the profit was higher than that at the steady-state around which the LEMPC was designed.

**Figure 6.** State trajectories under Lyapunov-based EMPC (LEMPC) with *Qextra* = 300 K/h starting at 0.5 h, where the LEMPC has not been made aware ("Unaware") and has been made aware ("Aware") of the change in the energy balance.

**Figure 7.** Input trajectories under LEMPC with *Qextra* = 300 K/h starting at 0.5 h, where the LEMPC has not been made aware ("Unaware") and has been made aware ("Aware") of the change in the energy balance.

The oscillatory behavior of the states before 0.5 h is caused by the fact that the profit is maximized for this process at the boundary of <sup>Ω</sup>*ρ*<sup>ˆ</sup>*e*. Without plant-model mismatch, the LEMPC is able to maintain the closed-loop state exactly on the boundary of <sup>Ω</sup>*ρ*<sup>ˆ</sup>*e* and therefore always operates the process using the constraint of Equation (22); however, when the plant-model mismatch occurs (induced by the use of different integration steps to simulate the process dynamic model within the LEMPC and for the simulation of the process under the computed control actions), the closed-loop state then exits <sup>Ω</sup>*ρ*<sup>ˆ</sup>*e* when the LEMPC predicts it will stay inside of it under the control actions computed by the controller. The result is that the constraint of Equation (23) is then activated until the closed-loop state reenters <sup>Ω</sup>*ρ*<sup>ˆ</sup>*e*. This process of entering <sup>Ω</sup>*ρ*<sup>ˆ</sup>*e*, attempting to operate at its boundary, and then being kicked out only to be driven back in is the cause of the oscillatory response of the states and inputs in Figures 6 and 7. It is noted, however, that though this behavior may be undesirable from, for example, an actuator wear perspective, it does not reflect a loss of closed-loop stability or a malfunction of the controller. The controller is in fact maintaining the closed-loop state within <sup>Ω</sup>*ρ*<sup>ˆ</sup> as it was designed to do; the fact that it does so in perhaps a visually unfamiliar fashion means that we have not specified in the control law that it should not do that, so it is not aware that an end user would find that behavior strange (if the oscillatory behavior is deemed undesirable, one could consider, for example, input rate of change constraints and potentially the benefits of the human response-based input rate of change strategy in the prior section for handling unexpected events).

In the case that the LEMPC is not aware of the change in the process dynamics, the profit is 32.7103, whereas when the LEMPC is aware of the change in the dynamics, the profit is 32.5833. Though these values are very close, an interesting note is that the profit when the LEMPC is not aware of the change in the underlying dynamics is slightly higher than when it is aware. Intuitively, one would expect an LEMPC with a more accurate process model to be able to locate a more economically optimal trajectory for the closed-loop state to follow than an LEMPC that cannot provide as accurate predictions. Part of the reason for the enhanced optimality in the case without knowledge of the change in the underlying dynamics, however, comes from the two-mode nature of LEMPC. In the case that the LEMPC is aware of the change in the underlying dynamics, it drives the closed-loop state to an operating condition that remains closer to the boundary of <sup>Ω</sup>*ρ*<sup>ˆ</sup>*e* after 0.5 h than when it is not aware of the change in the underlying dynamics due to the plant/model mismatch being different in the different cases. The result is that the process accesses regions of state-space that lead to higher profits when the LEMPC does not know about the change in the dynamics than if the LEMPC knows more about the process dynamics.

The remainder of this example focuses on elucidating the conservativeness of the proposed approach. Specifically, we now consider the Lyapunov function selected as *V*ˆ *q* = *<sup>x</sup>TPx*, with *P* given as follows:

$$P = \begin{bmatrix} 2000 & -10 \\ -10 & 3 \end{bmatrix} \tag{85}$$

Again, *hNL*,<sup>1</sup>(*x*) is designed such that *hNL*,1,1(*x*) = 0 kmol/m3, and *hNL*,1,2(*x*) is computed via Sontag's formula but saturated at the input bounds of Equation (82) if they are met. *ρ*ˆ and *ρ*ˆ*e* were taken to be 1300 and 975, respectively, and *<sup>ρ</sup>*<sup>ˆ</sup>*saf e* was set to 1800. The process state was initialized at *xinit* = [0 kmol/m<sup>3</sup> 0 K] *T*, with controller parameters *N* = 10 and Δ = 0.01 h. The process model of Equations (77) and (78) was integrated with the explicit Euler numerical integration method using an integration step size of 10−<sup>4</sup> h within the EMPC and with an integration step size of 10−<sup>5</sup> h to simulate the process. The constraint of the form of Equation (23) is enforced at *tk* when *x*(*tk*) ∈ <sup>Ω</sup>*ρ*ˆ/Ω*ρ*<sup>ˆ</sup>*e* but then followed by a constraint of the form of Equation (22) at the end of all sampling periods.

At 0.5 h, it is assumed that an additional source of heat arises outside the reactor such that the right-hand side of Equation (78) is modified by the addition of another heat term *Qextra* = 500 K/h. In this case, with no change in the process model used by the EMPC or even in the control law

(i.e., in contrast to the implementation strategy in Section 3.2.1, *hNL*,<sup>1</sup> is not employed when the closed-loop state exits <sup>Ω</sup>*ρ*ˆ), the behavior in Figure 8 results. Notably, the closed-loop state does not leave <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf e*, and no infeasibility issues occurred. In contrast, if we begin to utilize *hNL*,<sup>1</sup> when the closed-loop state leaves <sup>Ω</sup>*ρ*ˆ, the closed-loop state will eventually leave <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf e* (Figure 9). While we can obtain a new empirical model (in this case, we assume that the dynamics become fully known at 0.54 h and are accounted for completely to demonstrate the result) and can use that to update *hNL*,<sup>1</sup> to *hNL*,<sup>2</sup> (i.e., *hNL*,<sup>1</sup> but with modified saturation bounds to reflect design around the new steady-state of the system with *QAdded* = 500 K/h) before the closed-loop state leaves <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf e* as suggested in the implementation strategy in Section 3.2.1 (creating the profile shown in Figure 10 corresponding to 2 h of operation in which the closed-loop state is driven back to the origin under *hNL*,2), the fact that the closed-loop state would not have left the stability region if the controller had not been adjusted illustrates the conservativeness of the approach. We note that Figure 10 does not complete the implementation strategy in Section 3.2.1 (which would involve the use of a new LEMPC after the closed-loop state reenters <sup>Ω</sup>*ρ*<sup>ˆ</sup> for this example) because that part of the implementation strategy will be demonstrated in the discussion for a slightly different LEMPC presented below.

Finally, we provide a result where the LEMPC computes a time-varying input policy due to the desire to enforce a constraint on the amount of reactant available in the feed over an hour (i.e., a material/feedstock constraint) as follows:

$$\lim\_{\Gamma \to 0} \int\_{t=0}^{t=1} \ln\_1(\tau) d\tau = 0 \text{ kmol/m}^3 \tag{86}$$

This constraint is enforced via a soft constraint formulation by introducing slack variables *s*1 and *s*2 that are penalized in a modified objective function as follows:

$$\int\_{t\_k}^{t\_{k+N}} \left[ -k\_0 e^{-\frac{\Gamma}{k\_0 \Gamma(\tau)}} \mathbb{C}\_A(\tau)^2 \right] d\tau + 100(s\_1^2 + s\_2^2) \tag{87}$$

They are used in the following constraints:

$$\sum\_{i=0}^{k-1} \left( u\_1^\*(t\_i|t\_i) \right) + \sum\_{i=k}^{k+N\_k} \left( u\_1(t\_i|t\_k) \right) - 3.5\delta(100 - \frac{t\_k}{\Delta} - N) \le s\_1 \tag{88}$$

$$-\sum\_{i=0}^{k-1} \left( u\_1^\*(t\_i|t\_i\rangle) - \sum\_{i=k}^{k+N\_k} \left( u\_1(t\_i|t\_k\rangle) - 3.5\delta(100 - \frac{t\_k}{\Delta} - N) \right) \le s\_2 \tag{89}$$

where *Nk* = *N* and *δ* = 1 when *tk* < 0.9 h and where *δ* = 0 and *Nk* is the number of sampling periods left in a 1 h operating period when *tk* ≥ 0.9 h. These constraints are developed based on Reference [12]. *u*<sup>∗</sup> 1 (*ti*|*ti*) signifies the value of *u*1 applied to the process at a prior sampling time, and *<sup>u</sup>*1(*ti*|*tk*) reflects the value of *u*1 predicted at the current sampling time *tk* to be applied for *t* ∈ [*ti*, *ti*+<sup>1</sup>), *i* = *k*, ... , *k* + *Nk*. The upper and lower bounds on *s*1 and *s*2 were set to 2 × 10<sup>19</sup> and −2 × 1019, respectively, to allow them to be effectively unbounded. The initial guesses of the slack variables were set to 0 at each sampling time.

When the LEMPC with the above modifications is applied to the process with *QAdded* = 500 K/h starting at 0.5 h, the closed-loop state again exits <sup>Ω</sup>*ρ*<sup>ˆ</sup> for some time after 0.5 h but reenters it and also does not exit <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf e*, once again reflecting the conservatism from a closed-loop stability standpoint of a strategy that updates the process model whenever the closed-loop state leaves <sup>Ω</sup>*ρ*ˆ. Furthermore, if *hNL*,<sup>1</sup> is utilized after it is detected that the closed-loop state leaves <sup>Ω</sup>*ρ*<sup>ˆ</sup> (the first sampling time at which this occurs is 0.51 h), then it exits <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf e* by 0.52 h, showing that the length of the sampling period or the size of <sup>Ω</sup>*ρ*<sup>ˆ</sup> with respect to <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf e* is not sufficiently small enough to impose model updates before closed-loop stability is jeopardized because measurements are only available every sampling time. If instead, however, *ρ*ˆ is updated to be 1200 and *ρ*ˆ*e* is set to 900, then the closed-loop state remains in <sup>Ω</sup>*ρ*<sup>ˆ</sup> between 0.51 and 0.52 h. If at 0.52 h, we assume that the new dynamics (i.e., with *QAdded* = 500 K/h) become available and are used in designing *hNL*,<sup>2</sup> (used from 0.52 h until the first sampling time at which *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup> again) and that a second LEMPC designed based on the updated model is used after the closed-loop state has reentered <sup>Ω</sup>*ρ*ˆ, the state-space trajectory in Figure 11 results.

**Figure 8.** State-space plot under LEMPC with *Qextra* = 500 K/h starting at 0.5 h and no change in the control law or model in response.

**Figure 9.** State-space plot under LEMPC with *Qextra* = 500 K/h starting at 0.5 h and the control law switched to *hNL*,<sup>1</sup> in response to the closed-loop state leaving <sup>Ω</sup>*ρ*ˆ.

**Figure 10.** State-space plot under LEMPC with *Qextra* = 500 K/h starting at 0.5 h and the control law switched to *hNL*,<sup>1</sup> in response to the closed-loop state leaving <sup>Ω</sup>*ρ*<sup>ˆ</sup> and then switched to *hNL*,<sup>2</sup> at 0.54 h.

**Figure 11.** State-space plot under LEMPC with *Qextra* = 500 K/h starting at 0.5 h and the control law switched to *hNL*,<sup>1</sup> in response to the closed-loop state leaving <sup>Ω</sup>*ρ*ˆ, then switched to *hNL*,<sup>2</sup> at 0.52 h, and finally switched back to an LEMPC incorporating an updated process model after the closed-loop state reenters <sup>Ω</sup>*ρ*ˆ.

#### **4. Conclusions**

This work developed a Lyapunov-based EMPC framework for handling unexpected considerations of different types. One of the types of considerations handled was end-user response to how a control law operates a process, providing a controller self-update capability through input rate of change constraints that allows even uncertain or imprecise information about the end-user response to be used in optimizing the controller formulation without loss of closed-loop stability or feasibility. The second type of consideration was the occurrence of anomalies, where conditions which would guarantee that

the closed-loop state can be stabilized in the presence of an anomaly that changes the underlying process dynamics as long as a detection method identifies a new process model sufficiently quickly, were presented that uses the LEMPC stability properties in developing an anomaly detection mechanism. Chemical process examples were presented for both cases to demonstrate the proposed approach.

The work above provides insights into interpretability and verification considerations for EMPC from a theoretical perspective. However, these remain significant challenges for this control design. For example, there is no guarantee that adjusting a given constraint (e.g., adjusting the upper bound on an input rate of change constraint) will cause process behavior to appear interpretable to an end user before it approaches steady-state behavior, which may reduce the benefits of using EMPC. Furthermore, the results related to anomaly handling were demonstrated via process examples to be highly conservative. No methods were presented for practically ascertaining time (online) until an anomaly would result in the closed-loop state leaving a known region of state-space after detection to facilitate appropriate actions to be taken. Further work on these issues needs to be undertaken to develop practical EMPC designs with appropriate safety and interpretability properties with low time required to verify the designs before putting them into the field for different processes.

**Funding:** Financial support from the National Science Foundation CBET-1839675, from the Air Force Office of Scientific Research award number FA9550-19-1-0059, and from Wayne State University is gratefully acknowledged.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Appendix A. Proof of Proposition 4**

**Proof.** The result in Equation (62) is stated in Proposition 1; therefore, it remains to prove that Equation (63) holds. To derive the result of Equation (63), Equations (59) and (60) are integrated as follows:

$$\mathbb{1}\_{a,i+1,q}(t) = \mathbb{1}\_{a,i,q}(t\_{i,i+1}) + \int\_{t\_{i,l}+1}^{t} \vec{f}\_{i+1,q}(\mathbb{1}\_{a,i+1,q}(s), \vec{a}\_q(s), w\_{i+1}(s)) ds \tag{A1}$$

$$\mathfrak{a}\_{b,\emptyset}(t) = \mathfrak{x}\_{b,\emptyset}(t\_{s,i+1}) + \int\_{t\_{s,i+1}}^{t} \tilde{f}\_{\mathcal{N}L,\emptyset}(\mathfrak{x}\_{b,\emptyset}(s), \mathfrak{u}\_{\emptyset}(s))ds\tag{A2}$$

¯

for *t* ∈ [*ts*,*i*+1, *<sup>t</sup>*1]. Subtracting Equation (A2) from Equation (A1) and taking norms of both sides of the resulting equation gives the following:

¯


From Equations (15), (52), and (61), we have the following:


Using Equation (50) we ge<sup>t</sup> the following,


#### **Appendix B. Proof of Theorem 2**

**Proof.** To guarantee the results, recursive feasibility of the LEMPC must hold. Feasibility of the LEMPC of Equation (24) follows from Theorem 1 when *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* . Subsequently, closed-loop stability must be proven both when *ts*,*i*+<sup>1</sup> = *tk* and when *ts*,*i*+<sup>1</sup> ∈ (*tk*, *tk*+<sup>1</sup>).

Consider first the case that *ts*,*i*+<sup>1</sup> = *tk*. In this case, if Equation (68) holds with *p* = *i* + 1 and *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* , then *x*(*t*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* from Theorem 1 for *t* ≥ 0. Consider second the case that *ts*,*i*+<sup>1</sup> ∈ (*tk*, *tk*+<sup>1</sup>). In this case, until *ts*,*i*+1, if Equations (68) and (69) hold for *p* = *i*, the closed-loop state is maintained within <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* from Theorem 1. To guarantee that the closed-loop state is maintained in <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* after *ts*,*i*+<sup>1</sup> until *tk*+1, it is first noted that, if *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>* and *ts*,*i*+<sup>1</sup> ∈ (*tk*, *tk*+<sup>1</sup>), then from Proposition 2, we have the following:

$$\mathcal{V}\_{\boldsymbol{q}}(\pounds\_{\boldsymbol{a},i+1,\boldsymbol{q}}(t)) \leq \mathcal{V}\_{\boldsymbol{q}}(\pounds\_{\boldsymbol{b},\boldsymbol{q}}(t\_{k+1})) + f\_{V,\boldsymbol{q}}(|\pounds\_{\boldsymbol{a},i+1,\boldsymbol{q}}(t) - \pounds\_{\boldsymbol{b},\boldsymbol{q}}(t\_{k+1})|) \tag{A6}$$

if *<sup>x</sup>*¯*a*,*i*+1,*q*(*t*), *<sup>x</sup>*¯*b*,*<sup>q</sup>*(*t*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* for *t* ∈ [*tk*, *tk*+<sup>1</sup>]. If Proposition 4 holds, then from Equation (24f), we have the following:

$$\dot{V}\_{q}(\mathbb{1}\_{4\check{\rho}+1,q}(t)) \le \dot{\rho}\_{\varepsilon,\check{q}} + f\_{V,\emptyset}(f\_{W,l,q}(t\_{s,i+1}-t\_k) + (M\_{\text{clMRage},l,q})(t-t\_{s,i+1}) + \frac{L\_{\text{ol},\check{q}}\theta\_{l} + M\_{\text{ol},\check{q},i}}{L\_{\text{cl},\check{q}}}(\epsilon^{L\_{\text{cl},\check{q}}l} - \epsilon^{L\_{\text{ol},\check{q}}l,i\_{\hat{\phi}}+1})) \tag{A7}$$

If Equation (70) holds, then *V* ˆ *<sup>q</sup>*(*x*¯*a*,*i*+1,*q*(*t*)) ≤ *ρ*ˆ*q* for *t* ∈ [*ts*,*i*+1, *tk*+<sup>1</sup>].

If instead *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*/Ω*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>* and if Equations (68) and (69) hold, the closed-loop state is maintained within <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* from Theorem 1 until *ts*,*i*+1. To guarantee that the closed-loop state is maintained in <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* after *ts*,*i*+<sup>1</sup> until *tk*+1, it is first noted that the following is true:

$$\begin{split} & \frac{\partial \hat{V}\_{q}(\mathbf{x}(t\_{k}))}{\partial \mathbf{x}} \left( f\_{\text{NL},q}(\mathbf{x}(t\_{k}), \hat{u}\_{q}(t\_{k})) \right) \\ & \leq \frac{\partial \hat{V}\_{q}(\mathbf{x}(t\_{k}))}{\partial \mathbf{x}} \left( f\_{\text{NL},q}(\mathbf{x}(t\_{k}), h\_{\text{NL},q}(\mathbf{x}(t\_{k}))) \right) \leq -h\_{3,q}(|\mathbf{x}(t\_{k})|) \end{split} \tag{A8}$$

from Equation (12b) and Equation (24g). When *tk* ≤ *t* < *ts*,*i*+1, then from Reference [51], if Equation (68) and the conditions of Theorem 2 hold with *p* = *i*, the following is true:

*<sup>∂</sup>V*<sup>ˆ</sup>*q*(*x*¯*a*,*i*,*<sup>q</sup>*(*τ*)) *∂x* ( ¯ *fi*,*<sup>q</sup>*(*x*¯*a*,*i*,*<sup>q</sup>*(*τ*), *<sup>u</sup>*¯*q*(*tk*), *wi*(*τ*))) ≤ −*α*<sup>ˆ</sup> 3,*q*(*α*<sup>ˆ</sup> −1 2,*q* (*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>*)) + *α*ˆ 4,*q*(*α*<sup>ˆ</sup> −1 1,*q* (*ρ*<sup>ˆ</sup>*q*))*Merr*,*i*,*<sup>q</sup>* + *<sup>L</sup>x*,*i*,*<sup>q</sup>Mi*<sup>Δ</sup> + *<sup>L</sup>w*,*i*,*qθ<sup>i</sup>* (A9)

for *τ* ∈ [*tk*, *ts*,*i*+<sup>1</sup>), and

$$
\hat{\mathcal{V}}\_{\emptyset}(\pounds\_{a,i,\emptyset}(t\_{s,i+1})) \le \hat{\mathcal{V}}\_{\emptyset}(x(t\_k)) \tag{A10}
$$

> Given that *<sup>x</sup>*¯*a*,*i*,*<sup>q</sup>*(*ts*,*i*+<sup>1</sup>) = *<sup>x</sup>*¯*a*,*i*+1,*q*(*ts*,*i*+<sup>1</sup>), the following holds:

*<sup>∂</sup>V*<sup>ˆ</sup>*q* (*x*¯*a*,*i*+1,*<sup>q</sup>* (*ts*,*i*+<sup>1</sup>)) *∂x* ( ¯ *fi*+1,*q*(*x*¯*a*,*i*+1,*q*(*ts*,*i*+<sup>1</sup>), *<sup>u</sup>*¯*q*(*tk* ), 0)) = *<sup>∂</sup>V*<sup>ˆ</sup>*q* (*x*¯*a*,*i*+1,*<sup>q</sup>* (*ts*,*i*+<sup>1</sup>)) *∂x* ( ¯ *fi*+1,*q*(*x*¯*a*,*i*+1,*q*(*ts*,*i*+<sup>1</sup>), *<sup>u</sup>*¯*q*(*tk* ), 0)) + *<sup>∂</sup>V*<sup>ˆ</sup>*q* (*x*¯*a*,*i*,*<sup>q</sup>* (*ts*,*i*+<sup>1</sup>)) *∂x* ( ¯*fi*,*<sup>q</sup>*(*x*¯*a*,*i*,*<sup>q</sup>*(*ts*,*i*+<sup>1</sup>), *<sup>u</sup>*¯*q*(*tk* ), 0)) − *<sup>∂</sup>V*<sup>ˆ</sup>*q* (*x*¯*a*,*i*,*<sup>q</sup>* (*ts*,*i*+<sup>1</sup>)) *∂x* ( ¯ *fi*,*<sup>q</sup>*(*x*¯*a*,*i*,*<sup>q</sup>*(*ts*,*i*+<sup>1</sup>), *<sup>u</sup>*¯*q*(*tk* ), 0)) ≤ −*α*<sup>ˆ</sup> 3,*q*(*α*<sup>ˆ</sup> −1 2,*q* (*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>*)) + *α*ˆ 4,*q*(*α*<sup>ˆ</sup> −1 1,*q* (*ρ*<sup>ˆ</sup>*q*))*Merr*,*i*,*<sup>q</sup>* + *<sup>L</sup>x*,*i*,*<sup>q</sup>Mi*<sup>Δ</sup> + *<sup>L</sup>w*,*i*,*qθ<sup>i</sup>* + ---- *<sup>∂</sup>V*<sup>ˆ</sup>*q* (*x*¯*a*,*i*+1,*<sup>q</sup>* (*ts*,*i*+<sup>1</sup>)) *∂x* ( ¯*fi*+1,*q*(*x*¯*a*,*i*+1,*q*(*ts*,*i*+<sup>1</sup>), *<sup>u</sup>*¯*q*(*tk* ), 0)) − *<sup>∂</sup>V*<sup>ˆ</sup>*q* (*x*¯*a*,*i*,*<sup>q</sup>* (*ts*,*i*+<sup>1</sup>)) *∂x* ( ¯ *fi*,*<sup>q</sup>*(*x*¯*a*,*i*,*<sup>q</sup>*(*ts*,*i*+<sup>1</sup>), *<sup>u</sup>*¯*q*(*tk* ), 0)) ---- ≤ −*α*<sup>ˆ</sup> 3,*q*(*α*<sup>ˆ</sup> −1 2,*q* (*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>*)) + *α*ˆ 4,*q*(*α*<sup>ˆ</sup> −1 1,*q* (*ρ*<sup>ˆ</sup>*q*))*Merr*,*i*,*<sup>q</sup>* + *<sup>L</sup>x*,*i*,*<sup>q</sup>Mi*<sup>Δ</sup> + *<sup>L</sup>w*,*i*,*qθ<sup>i</sup>* + ---- *<sup>∂</sup>V*<sup>ˆ</sup>*q* (*x*¯*a*,*i*,*<sup>q</sup>* (*ts*,*i*+<sup>1</sup>)) *∂x* ---- -- ¯*fi*+1,*q*(*x*¯*a*,*i*+1,*q*(*ts*,*i*+<sup>1</sup>), *<sup>u</sup>*¯*q*(*tk* ), 0) − ¯*fi*,*<sup>q</sup>*(*x*¯*a*,*i*,*<sup>q</sup>*(*ts*,*i*+<sup>1</sup>), *<sup>u</sup>*¯*q*(*tk* ), 0)-- ≤ −*α*<sup>ˆ</sup> 3,*q*(*α*<sup>ˆ</sup> −1 2,*q* (*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>*)) + *α*ˆ 4,*q*(*α*<sup>ˆ</sup> −1 1,*q* (*ρ*<sup>ˆ</sup>*q*))*Merr*,*i*,*<sup>q</sup>* + *<sup>L</sup>x*,*i*,*<sup>q</sup>Mi*<sup>Δ</sup> + *<sup>L</sup>w*,*i*,*qθ<sup>i</sup>* + *α*ˆ 4,*q*(|*x*¯*a*,*i*,*<sup>q</sup>*(*ts*,*i*+<sup>1</sup>)|)*Mchange*,*i*,*<sup>q</sup>* ≤ −*α*<sup>ˆ</sup> 3,*q*(*α*<sup>ˆ</sup> −1 2,*q* (*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>*)) + *α*ˆ 4,*q*(*α*<sup>ˆ</sup> −1 1,*q* (*ρ*<sup>ˆ</sup>*q*))*Merr*,*i*,*<sup>q</sup>* + *<sup>L</sup>x*,*i*,*<sup>q</sup>Mi*<sup>Δ</sup> + *<sup>L</sup>w*,*i*,*qθ<sup>i</sup>* + *α*ˆ 4,*q*(*α*<sup>ˆ</sup> −1 1,*q* (*ρ*<sup>ˆ</sup>*q*))*Mchange*,*i*,*<sup>q</sup>* (A11)

where the last inequality follows from the fact that *<sup>x</sup>*¯*a*,*i*,*<sup>q</sup>*(*ts*,*i*+<sup>1</sup>) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* if *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* when Equations (68) and (69) hold according to Theorem 1.

Finally, for *τ* ∈ [*ts*,*i*+1, *tk*+<sup>1</sup>),

*<sup>∂</sup>V*<sup>ˆ</sup>*q* (*x*¯*a*,*i*+1,*<sup>q</sup>* (*τ*)) *∂x* ( ¯ *fi*+1,*q*(*x*¯*a*,*i*+1,*q*(*τ*), *<sup>u</sup>*¯*q*(*tk* ), *wi*+<sup>1</sup>(*τ*)) = *<sup>∂</sup>V*<sup>ˆ</sup>*q* (*x*¯*a*,*i*+1,*<sup>q</sup>* (*τ*)) *∂x* ( ¯ *fi*+1,*q*(*x*¯*a*,*i*+1,*q*(*τ*), *<sup>u</sup>*¯*q*(*tk* ), *wi*+<sup>1</sup>(*τ*)) + *<sup>∂</sup>V*<sup>ˆ</sup>*q* (*x*¯*a*,*i*+1,*<sup>q</sup>* (*ts*,*i*+<sup>1</sup>)) *∂x* ( ¯*fi*+1,*q*(*x*¯*a*,*i*+1,*q*(*ts*,*i*+<sup>1</sup>), *<sup>u</sup>*¯*q*(*tk* ), 0)) − *<sup>∂</sup>V*<sup>ˆ</sup>*q* (*x*¯*a*,*i*+1,*<sup>q</sup>* (*ts*,*i*+<sup>1</sup>)) *∂x* ( ¯ *fi*+1,*q*(*x*¯*a*,*i*+1,*q*(*ts*,*i*+<sup>1</sup>), *<sup>u</sup>*¯*q*(*tk* ), 0)) ≤ −*α*<sup>ˆ</sup> 3,*q*(*α*<sup>ˆ</sup> −1 2,*q* (*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>*)) + *α*ˆ 4,*q*(*α*<sup>ˆ</sup> −1 1,*q* (*ρ*<sup>ˆ</sup>*q*))*Merr*,*i*,*<sup>q</sup>* + *<sup>L</sup>x*,*i*,*<sup>q</sup>Mi*<sup>Δ</sup> + *<sup>L</sup>w*,*i*,*qθ<sup>i</sup>* + *α*ˆ 4,*q*(*α*<sup>ˆ</sup> −1 1,*q* (*ρ*<sup>ˆ</sup>*q*))*Mchange*,*i*,*<sup>q</sup>* + ---- *<sup>∂</sup>V*<sup>ˆ</sup>*q* (*x*¯*a*,*i*+1,*<sup>q</sup>* (*τ*)) *∂x* ( ¯*fi*+1,*q*(*x*¯*a*,*i*+1,*q*(*τ*), *<sup>u</sup>*¯*q*(*tk* ), *wi*+<sup>1</sup>(*τ*)) − *<sup>∂</sup>V*<sup>ˆ</sup>*q* (*x*¯*a*,*i*+1,*<sup>q</sup>* (*ts*,*i*+<sup>1</sup>)) *∂x* ( ¯*fi*+1,*q*(*x*¯*a*,*i*+1,*q*(*ts*,*i*+<sup>1</sup>), *<sup>u</sup>*¯*q*(*tk* ), 0))---- ≤ −*α*<sup>ˆ</sup> 3,*q*(*α*<sup>ˆ</sup> −1 2,*q* (*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>*)) + *α*ˆ 4,*q*(*α*<sup>ˆ</sup> −1 1,*q* (*ρ*<sup>ˆ</sup>*q*))*Merr*,*i*,*<sup>q</sup>* + *<sup>L</sup>x*,*i*,*<sup>q</sup>Mi*<sup>Δ</sup> + *<sup>L</sup>w*,*i*,*qθ<sup>i</sup>* + *α*ˆ 4,*q*(*α*<sup>ˆ</sup> −1 1,*q* (*ρ*<sup>ˆ</sup>*q*))*Mchange*,*i*,*<sup>q</sup>* <sup>+</sup>*Lx*,*i*+1,*q*|*x*¯*a*,*i*+1,*q*(*τ*) − *<sup>x</sup>*¯*a*,*i*+1,*q*(*ts*,*i*+<sup>1</sup>)| + *<sup>L</sup>w*,*i*+1,*qθi*+<sup>1</sup> ≤ −*α*<sup>ˆ</sup> 3,*q*(*α*<sup>ˆ</sup> −1 2,*q* (*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>*)) + *α*ˆ 4,*q*(*α*<sup>ˆ</sup> −1 1,*q* (*ρ*<sup>ˆ</sup>*q*))*Merr*,*i*,*<sup>q</sup>* + *<sup>L</sup>x*,*i*,*<sup>q</sup>Mi*<sup>Δ</sup> + *<sup>L</sup>w*,*i*,*qθ<sup>i</sup>* + *α*ˆ 4,*q*(*α*<sup>ˆ</sup> −1 1,*q* (*ρ*<sup>ˆ</sup>*q*))*Mchange*,*i*,*<sup>q</sup>* <sup>+</sup>*Lx*,*i*+1,*qMi*+1<sup>Δ</sup> + *<sup>L</sup>w*,*i*+1,*qθi*+<sup>1</sup> (A12)

If Equation (71) holds, then integrating Equation (A12) gives that *V* ˆ *<sup>q</sup>*(*x*¯*a*,*i*+1,*q*(*t*)) ≤ *<sup>V</sup>*<sup>ˆ</sup>*q*(*x*¯*a*,*i*,*<sup>q</sup>*(*ts*,*i*+<sup>1</sup>)), for all *t* ∈ [*ts*,*i*+1, *tk*+<sup>1</sup>]. Since *<sup>x</sup>*¯*a*,*i*+1,*q*(*ts*,*i*+<sup>1</sup>) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* , this guarantees that the closed-loop state remains in <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* even after the switch in the process model occurs, regardless of whether it occurs at a sampling time or throughout a sampling period, when the conditions of the theorem hold.

#### **Appendix C. Proof of Theorem 3**

**Proof.** This proof consists of several parts. First, recursive feasibility of the LEMPC of Equation (24) until *td*,*<sup>q</sup>* is presented. Second, it is demonstrated that, after *ts*,*i*+<sup>1</sup> and before *td*,*q*, the closed-loop state is maintained in <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* under the conditions of the theorem. Third, it is demonstrated that, after *td*,*q*, the closed-loop state will be maintained in <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*for a number of sampling periods given by *th*,*q*.

*Part 1.* Until *td*,*q*, each state measurement provided to the LEMPC of Equation (24) is within <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* . From Reference [51], under the conditions of Equations (65) and (66), this guarantees feasibility of the LEMPC of Equation (24). After *td*,*q*, when the closed-loop state exits <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* , feasibility is no longer guaranteed for the LEMPC of Equation (24) but *hNL*,*<sup>q</sup>* is then used instead according to the statement of the theorem so that a characterizable control law is always used.

*Part 2.* Until *ts*,*i*+1, closed-loop stability within <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* is guaranteed under the LEMPC of Equation (24) under the conditions in Equations (65) and (66) from Reference [51]. Subsequently, until *td*,*q*, it must be demonstrated that, if the state measurement is contained within <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* at *tk*, then *x*(*t*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* ⊂ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf<sup>e</sup>*,*q* ,

*t* ∈ [*tk*, *tk*+<sup>1</sup>]. Here, one of two cases holds: either *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>* or *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*/Ω*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>* . The state of the underlying model before *ts*,*i*+<sup>1</sup> is denoted by *<sup>x</sup>*¯*a*,*i*,*<sup>q</sup>* and, after, is *<sup>x</sup>*¯*a*,*i*,*q*+1.

If *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>* and if *ts*,*i*+<sup>1</sup> ∈ [*tk*, *tk*+<sup>1</sup>), from Propositions 1 and 2 and Equation (24f), we have the following:

$$\begin{split} \dot{V}\_{\boldsymbol{q}}(\mathbb{x}\_{\boldsymbol{a},i,\boldsymbol{q}}(t)) &\leq \dot{V}\_{\boldsymbol{q}}(\mathbb{x}\_{\boldsymbol{b},\boldsymbol{q}}(t)) + f\_{\boldsymbol{V},\boldsymbol{q}}(|\mathbb{x}\_{\boldsymbol{a},i,\boldsymbol{q}}(t) - \mathbb{x}\_{\boldsymbol{b},\boldsymbol{q}}(t)|) \\ &\leq \dot{\rho}\_{\boldsymbol{e},\boldsymbol{q}} + f\_{\boldsymbol{V},\boldsymbol{q}}(f\_{\boldsymbol{W},\boldsymbol{i},\boldsymbol{q}}(\Delta)) \leq \dot{\rho}\_{\boldsymbol{q}} \end{split} \tag{A13}$$

for *t* ∈ [*tk*, *ts*,*i*+<sup>1</sup>) when Equation (65) holds, and

$$\begin{split} \mathcal{V}\_{\boldsymbol{\theta}}(\mathbb{z}\_{\boldsymbol{d},i+1,\boldsymbol{q}}(t)) &\leq \mathcal{V}\_{\boldsymbol{q}}(\mathbb{z}\_{\boldsymbol{b},\boldsymbol{q}}(t)) + f\_{V,\boldsymbol{q}}(|\mathbb{z}\_{\boldsymbol{d},i+1}(t) - \mathbb{z}\_{\boldsymbol{b},\boldsymbol{q}}(t)|) \\ &\leq \widehat{\rho}\_{\boldsymbol{\varepsilon},\boldsymbol{q}} + f\nu\_{\boldsymbol{\varepsilon}|}(f\_{V,\boldsymbol{\delta}|}(t\_{\boldsymbol{\delta},i+1} - t\_{\boldsymbol{k}}) + (M\_{\text{charge},i,\boldsymbol{q}})(t - t\_{\boldsymbol{\delta},i+1}) + \frac{L\_{\text{w},i,\boldsymbol{q}}\theta\_{i} + M\_{\text{tr},i,\boldsymbol{q}}}{L\_{\text{x},i,\boldsymbol{q}}}(\boldsymbol{\varepsilon}^{\mathcal{L}\_{\boldsymbol{d},i},\boldsymbol{t}} - \boldsymbol{\varepsilon}^{\mathcal{L}\_{\boldsymbol{d},\boldsymbol{q}}\boldsymbol{t}\_{\boldsymbol{d}} + 1})) \end{split} \tag{A14}$$

for *t* ∈ [*ts*,*i*+1, *tk*+<sup>1</sup>) from Proposition 4. From the conditions in Equation (73), this gives that *<sup>V</sup>*<sup>ˆ</sup>*q*(*x*(*t*)) is maintained within <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* for all *t* ∈ [*tk*, *tk*+<sup>1</sup>).

If instead *ts*,*i*+<sup>1</sup> occurs before or at *tk*, then *<sup>x</sup>*¯*b*,*<sup>q</sup>*(*tk*) = *<sup>x</sup>*¯*a*,*i*+1,*q*(*tk*) and Propositions 1 and 2 and Equation (24f) give the following:

$$\begin{split} \mathcal{V}\_{\boldsymbol{q}}(\mathfrak{x}\_{\boldsymbol{d},i+1,\boldsymbol{q}}(t)) &\leq \mathcal{V}\_{\boldsymbol{q}}(\mathfrak{x}\_{\boldsymbol{b},\boldsymbol{q}}(t)) + f\_{\boldsymbol{V},\boldsymbol{q}}(f\_{\boldsymbol{W},i+1,\boldsymbol{q}}(\boldsymbol{\Delta})) \\ &\leq \rho\_{\boldsymbol{\varepsilon},\boldsymbol{q}} + f\_{\boldsymbol{V},\boldsymbol{q}}(f\_{\boldsymbol{W},i+1,\boldsymbol{q}}(\boldsymbol{\Delta})) \end{split} \tag{A15}$$

for all *t* ∈ [*tk*, *tk*+<sup>1</sup>). From the conditions in Equation (75), this gives that *<sup>V</sup>*<sup>ˆ</sup>*q*(*x*(*t*)) is maintained within <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* for all *t* ∈ [*tk*, *tk*+<sup>1</sup>).

If *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*/Ω*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>* , then the constraint of Equation (24g) is used. In this case, we consider the cases where *ts*,*i*+<sup>1</sup> ∈ [*tk*, *tk*+<sup>1</sup>) and the case where *ts*,*i*+<sup>1</sup> occurs before *tk*, separately.

When *ts*,*i*+<sup>1</sup> ∈ [*tk*, *tk*+<sup>1</sup>), then before *ts*,*i*+1, Equation (24g) holds. From Reference [51], Equation (66) with Equation (67) cause *<sup>x</sup>*¯*a*,*i*,*<sup>q</sup>*(*t*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* for *t* ∈ [*tk*, *ts*,*i*+<sup>1</sup>). Subsequently, this result no longer holds because the underlying dynamic model changed so that Equation (24g) no longer provides an indication of the conditions which the closed-loop state meets, and a worst-case scenario in which the closed-loop state could subsequently move out of <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* is considered. Specifically, the first inequality in Equation (A14) continues to hold. Equation (24f) does not necessarily hold but instead it is guaranteed [51] that *<sup>x</sup>*¯*b*,*<sup>q</sup>*(*t*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q* under Equations (66) and (67), so that *V* ˆ *<sup>q</sup>*(*x*¯*b*,*<sup>q</sup>*) ≤ *ρ*<sup>ˆ</sup>*q*. Then, if Equation (74) holds, extending the first inequality in Equation (A14) guarantees that *V* ˆ *<sup>q</sup>*(*x*¯*a*,*i*+1,*q*(*t*)) ≤ *ρ*<sup>ˆ</sup>*samp*,*q*, for *t* ∈ [*ts*,*i*+1, *tk*+<sup>1</sup>). Therefore, throughout a sampling period containing *ts*,*i*+1, the closed-loop state does not leave <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* . If instead *ts*,*i*+<sup>1</sup> is before *tk*, then Equation (24g) is activated at *tk* and when *<sup>x</sup>*¯*a*,*i*+1,*q*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*/Ω*ρ*<sup>ˆ</sup>*s*,*<sup>q</sup>* [51]:

$$\begin{split} &\frac{\partial \bar{V}\_{\emptyset}(\tilde{x}\_{d,i+1,\emptyset}(\mathbf{r}))}{\partial \mathbf{x}} \left( \bar{f}\_{i+1,\emptyset}(\mathbb{X}\_{\emptyset,i+1,\emptyset}(\mathbf{r}),\bar{u}\_{\emptyset}(\mathbf{r}\_{k}),w\_{i+1}(\mathbf{r})) \right. \\ &\leq -\hbar \mathfrak{z}\_{\emptyset,\emptyset}(\mathbb{A}^{-1}\_{2,\emptyset}(\mathfrak{\boldsymbol{\mu}}\_{\emptyset,\emptyset})) + \hbar \mathfrak{z}\_{\emptyset,\emptyset}(\mathbb{A}^{-1}\_{1,\emptyset}(\mathfrak{\boldsymbol{\mu}}\_{\emptyset})) M\_{\mathbf{z}\mathbf{r},i+1,\emptyset} + L'\_{\mathbf{x},i+1,\emptyset} M\_{i+1} \Delta + L'\_{\mathbf{u},i+1,\emptyset} \theta\_{i+1} \end{split} \tag{A16}$$

When Equation (72) is satisfied,

$$\frac{\partial \bar{V}\_{\boldsymbol{q}}(\boldsymbol{x}\_{\boldsymbol{a},i+1,\boldsymbol{q}}(\boldsymbol{\tau}))}{\partial \boldsymbol{x}} \; (\bar{f}\_{i+1,\boldsymbol{q}}(\bar{\boldsymbol{x}}\_{\boldsymbol{a},i+1,\boldsymbol{q}}(\boldsymbol{\tau}), \bar{u}\_{\boldsymbol{q}}(\boldsymbol{t}\_{k}), \boldsymbol{w}\_{i+1}(\boldsymbol{\tau})) \leq c\_{W,i+1,\boldsymbol{q}}/\Delta \tag{A17}$$

or

$$\hat{\mathcal{V}}\_q(\pounds\_{d,i+1,q}(t)) \le \hat{\mathcal{V}}\_q(\pounds\_{d,i+1,q}(t\_k)) + \frac{\epsilon\_{\emptyset,\ell,i+1,q}}{\Delta}(t - t\_k) \tag{A18}$$

This indicates that *V* ˆ *q* is guaranteed to increase at a worst-case rate along the closed-loop state trajectories under the control actions determined by the LEMPC of Equation (24) if the condition of Equation (72) is satisfied after an anomaly occurs. To ensure that, at the end of the sampling period, *V* ˆ *<sup>q</sup>*(*x*¯*a*,*i*+1,*q*(*t*)) ≤ *ρ*<sup>ˆ</sup>*samp*,*q*, given that *<sup>V</sup>*<sup>ˆ</sup>*q*(*x*¯*a*,*i*+1,*q*(*tk*)) ≤ *ρ*<sup>ˆ</sup>*q*, Equation (76) must hold. If *ts*,*i*+<sup>1</sup> is before *tk* but *x* ¯ *<sup>a</sup>*,*i*+1,*q*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*s*,*<sup>q</sup>* , then if *<sup>ρ</sup>*ˆmin,*<sup>i</sup>*+1,*q* ⊂ *ρ*<sup>ˆ</sup>*samp*,*q*, then *<sup>x</sup>*¯*a*,*i*+1,*q*(*t*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* from Equation (67).

Thus, whether *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>* or *x*(*tk*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*/Ω*ρ*<sup>ˆ</sup>*e*,*<sup>q</sup>* , *<sup>x</sup>*(*tk*+<sup>1</sup>) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* . Applying this recursively indicates that, from *ts*,*i*+<sup>1</sup> until *td*,*q*, the closed-loop state is maintained within <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* . This also indicates that *V* ˆ *<sup>q</sup>*(*x*¯*a*,*i*+1,*q*(*td*,*<sup>q</sup>*)) ≤ *ρ*<sup>ˆ</sup>*samp*,*q*. Because <sup>Ω</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* ⊂ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q* , *<sup>x</sup>*¯*a*,*i*+1,*q*(*td*,*<sup>q</sup>*) ∈ <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q* as well.

*Part 3.* At *td*,*q*, *hNL*,*<sup>q</sup>* in sample-and-hold begins to be used to control the process. Again, Equations (A16)–(A18) hold.

The time *tout*,*<sup>q</sup>* at which the closed-loop state reaches <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q* (i.e., when *V* ˆ *<sup>q</sup>*(*x*¯*a*,*i*+1,*q*(*tout*,*<sup>q</sup>*)) = *<sup>ρ</sup>*<sup>ˆ</sup>*saf <sup>e</sup>*,*<sup>q</sup>*) when initialized from *V* ˆ *<sup>q</sup>*(*x*¯*a*,*i*+1,*q*(*tk*)) = *ρ*<sup>ˆ</sup>*samp*,*q*, where *ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>* ≤ *<sup>ρ</sup>*<sup>ˆ</sup>*saf <sup>e</sup>*,*q*, is at least (*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q*<sup>−</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>*)<sup>Δ</sup> *W*,*i*+1,*<sup>q</sup>* + *tk*. To ensure that the time between *tk* and *tout*,*<sup>q</sup>* is no greater than (*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q*<sup>−</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>*)<sup>Δ</sup> *W*,*i*+1,*<sup>q</sup>* , the number of sampling periods available after *td*,*<sup>q</sup>* until the model needs to be updated with one which meets the conditions in Equation (66) with *i* set to *i* + 1 and *q* set to *q* + 1 is floor((*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q*<sup>−</sup>*ρ*<sup>ˆ</sup>*samp*,*<sup>q</sup>*) *W*,*i*+1,*<sup>q</sup>* ).

#### **Appendix D. Proof of Theorem 4**

**Proof.** If *hNL*,*q*+<sup>1</sup> is used to control the system after *tID*,*<sup>q</sup>* and the conditions of Theorem 4 are met, then *xa*,*i*+1,*q*(*tID*,*<sup>q</sup>*) = *xa*,*i*+1,*q*+<sup>1</sup>(*tID*,*<sup>q</sup>*), which lies in both <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q* and in <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q*+1 so that the closed-loop state has not left either region. From Reference [51], if Equation (66) is met for the *q* + 1/*i* + 1 model combination, then *hNL*,*q*+<sup>1</sup> causes *<sup>V</sup>*<sup>ˆ</sup>*q*+<sup>1</sup> to decrease so that it will not leave <sup>Ω</sup>*ρ*<sup>ˆ</sup>*saf <sup>e</sup>*,*q*+1 before the closed-loop state enters <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*+<sup>1</sup> . Once the closed-loop state enters <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*+<sup>1</sup> , then the LEMPC of Equation (24) is used with the *q* + 1 model, and if Equations (65) and (66) are met for the *q* + 1/*i* + 1 model combination, the closed-loop state is maintained in <sup>Ω</sup>*ρ*<sup>ˆ</sup>*q*+<sup>1</sup> from Reference [51].
