# Robust MPC with Zone Tracking

Zhiyinan Huanga, Jinfeng Liua,, Biao Huanga

aDepartment of Chemical & Materials Engineering, University of Alberta,

Edmonton, AB T6G 1H9, Canada

Corresponding author: J. Liu. Email<EMAIL_ADDRESS>
###### Abstract

We propose a robust nonlinear model predictive control design with generalized
zone tracking (ZMPC) in this work. The proposed ZMPC has guaranteed
convergence into the target zone in the presence of bounded disturbance. The
proposed approach achieves this by modifying the actual target zone such that
the effect of disturbances is rejected. A control invariant set (CIS) inside
the modified target zone is used as the terminal set, which ensures the
closed-loop stability of the proposed controller. Detailed closed-loop
stability analysis is presented. Simulation studies based on a continuous
stirred tank reactor (CSTR) are performed to validate the effectiveness of the
proposed ZMPC.

##  1 Introduction

Due to the development in computational hardware, nonlinear model predictive
control (MPC) is considered to be a promising advanced control strategy and
has attracted high research interests over the past decades. Inspired by the
linear quadratic regulator (LQR), MPC has the ability to solve nonlinear
optimization problems while considering constraints simultaneously. However
unlike LQR problems that can be solved explicitly for infinite horizon, the
same cannot be achieved for MPC due to the presence of constraints. Instead,
MPC is solved in a receding horizon manner [1, 2], where the optimization
problem is solved for a finite control horizon at each sampling instant with
only the very first control policy applied to the system. Conventional MPC
tracks a referencing trajectory, which is oftentimes an optimal steady-state
operating point [2]. This approach, however, may limit the flexibility of the
controller and sacrifice the process performance during transient operations.
One popular extension of the conventional MPC is economic MPC (EMPC), which
considers a general economic objective directly in dynamic optimizations [3].
Improved transient economic performance has been observed in EMPC compared to
MPC [4, 5].

Another approach that helps to improve the flexibility of conventional MPC is
zone MPC (ZMPC). ZMPC relaxes the set-point-tracking objective to a zone-
tracking one, which aims to drive the system into a bounded set [6, 7]. More
degrees of freedom are provided by ZMPC, which are beneficial for handling
multiple control objectives (e.g. tracking and economic) simultaneously.
Furthermore, ZMPC is potentially more robust in the presence of noise or
disturbance. Zone-tracking is also a natural objective that rises in many
real-world problems. For example, in agricultural practice, maintaining the
soil moisture in a certain range is often sufficient and more practical rather
than attempting to maintain it at a particular level. Various applications of
ZMPC have been reported in the literature, for example, treatment of diabetes
in [8, 9], control of heat system inside a building in [10], control of coal-
fired boiler-turbine generating systems in [11], and control of irrigation
systems in [12] and [13].

Closed-loop stability of the system is always an essential aspect to be
considered in controller design, as unstable controllers can lead to drastic
safety risks in process operation. The closed-loop stability theories of MPC
and EMPC are well-developed and widely used in literature based on the
Lyapunov stability theory [14, 15, 16, 17, 18]. However, most of these
stability theories rely on the assumption of the existence of an optimal
operating point for both MPC and EMPC. On the other hand, ZMPC is often
considered to be a relaxation to conventional MPC with not many investigations
done in terms of the stability property. The stability of ZMPC with a
secondary economic objective is investigated in [19]. An extension of [19]
considered the presence of disturbance, of which a ZMPC that tracks an
alternative economic zone is proposed [20]. In both works, the stability
theories rely on the assumption that there exists an optimal steady-state
operating condition. On the other hand, a generalized ZMPC was proposed with
no assumption on steady state optimal operating point in [21], however without
considering any system disturbance or noise.

In this work, we propose a robust ZMPC with guaranteed convergence in the
presence of bounded disturbance, which is an extension of the results obtained
in [21]. The proposed design rejects the effect of disturbance by modifying
the target set. The modified target set is a subset of the actual target set,
which the shrinkage is proved to be upper-bounded. A terminal constraint is
employed to ensure the closed-loop stability, which forces the system state to
converge to a forward control invariant set (CIS) inside the modified target
set. Based on the stability theory, a practical guideline for determining the
modified target zone is proposed.

On top of the generalized robust ZMPC, we consider a secondary economic
objective. Instead of introducing economic zones as proposed in [20], we
prioritize the zone tracking objective and establish stability based on the
ZMPC alone. This implies that the economic objective is optimized only if the
zone-tracking objective is satisfied, which provides a larger feasible region
and more flexibility, as the economic objective can be modified without
affecting the stability property of the controller.

The performance of the proposed ZMPC is validated through simulations based on
a continuous stirred tank reactor (CSTR). The performance of the generalized
ZMPC proposed in [21] is employed as the benchmark reference. To show the
effectiveness of the proposed design, simulations based on alternative
controller formulations are performed and the results are compared with those
obtained based on the proposed design. The effectiveness of the proposed
practical guideline is tested through simulations with different values of the
tuning parameter. A remark regarding the selection of the tuning parameter is
summarized.

The remainder of this script is organized as follows: Section 2 introduces the
preliminaries and the system of interest, Section 3 presents the proposed ZMPC
formulation, of which the stability proof is presented in Section 4. Based on
the stability theory, an algorithm for estimating the modified target set is
designed in Section 5. Section 6 introduces the proposed ZMPC with additional
economic objective. The simulation results obtained based on the CSTR are
presented in 7. Finally, concluding remarks and ideas for future works are
summarized in Section 8.

##  2 Preliminaries

###  2.1 Notation

Throughout this work, $|\cdot|$ denotes the Euclidean norm of a scalar or a
vector. The operator $||\cdot||_{n}$ denotes the $n$-norm of a scalar or a
vector. $\mathbb{I}_{M}^{N}$ represents the set of integers from $M$ to $N$
$\\{M,M+1,\cdots,N\\}$. $\mathbb{I}_{\geq 0}=\\{0,1,2,...\\}$ denotes the set
of non-negative integers. For two sets $A$ and $B$, the set subtraction is
defined as $A\backslash B=\\{a\in A|a\notin B\\}$. The operators
$\max{\\{\cdot\\}}$ and $\min{\\{\cdot\\}}$ find the maximum and minimum of
each element in a vector variable respectively.

The operator $a\odot b$ represents element-wise multiplication between vectors
$a$ and $b$.

###  2.2 Description and problem formulation

The following discrete-time nonlinear system is considered in this work:

$x(n+1)=f(x(n),u(n),w(n))$ (1)

where $x(\cdot)\in\mathbb{X}\subseteq\mathbb{R}^{n_{x}}$ denotes the state
vector, $u(\cdot)\in\mathbb{U}\subseteq\mathbb{R}^{n_{u}}$ denotes the control
input vector, and $w(\cdot)\in\mathbb{W}\subseteq\mathbb{R}^{n_{w}}$ denotes
the disturbance vector. $n\in\mathbb{I}_{\geq 0}$ denotes the time step.

A few assumptions are enforced on the system of interest throughout this work
and are listed below:

###### Assumption 1

The constraint sets $\mathbb{X}$ and $\mathbb{U}$ are compact and coupled such
that:

$(x(n),u(n))\in\mathbb{Z}\subseteq\mathbb{X}\times\mathbb{U}$

Furthermore, the sets can be expressed as following element-wise constraints:

$\mathbb{X}:=\left\\{x\Big{|}x_{lb}\leq x\leq
x_{ub}\right\\},\>\mathbb{U}:=\left\\{u\Big{|}u_{lb}\leq u\leq
u_{ub}\right\\}$

where $x_{lb}$, $x_{ub}$, $u_{lb}$, and $u_{ub}$ are constant vectors, with
each elements being the element-wise lower and upper bound of the state and
input. The expressions imply that the constraint sets are polyhedrons in the
corresponding space bounded by hyperplanes.

###### Assumption 2

The disturbance vector $w$ is bounded and contains the origin in its domain:

$\mathbb{W}:=\\{w\in\mathbb{R}^{n_{w}}:||w||_{\infty}\leq\theta,\theta>0\\}$
(2)

such that the infinity norm of $w$ is bounded by a positive constant $\theta$.

###### Assumption 3

The function
$f:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\times\mathbb{R}^{n_{w}}\rightarrow\mathbb{R}^{n_{x}}$
is locally Lipschitz with respect to $x$ and $w$ for all $x\in\mathbb{X}$,
$u\in\mathbb{U},w\in\mathbb{W}$. This implies there exist positive constants
$L_{x}$ and $L_{w}$ such that:

$|f(x,u,w)-f(z,u,0)|\leq L_{w}|w|+L_{x}|x-z|$ (3)

which indicates that the function $f$ cannot change drastically as $x$ and $w$
changes.

The control objective of this work is to drive the state of system (1) into a
tracking target set $\mathbb{X}_{t}$ and keeps it inside $\mathbb{X}_{t}$
thereafter in the presence of process uncertainties. The target set
$\mathbb{X}_{t}$ can be expressed as the following element-wise constraints:

$\mathbb{X}_{t}:=\left\\{x\Big{|}x^{t}_{lb}\leq x\leq x^{t}_{ub}\right\\}$ (4)

where $x^{t}_{lb}$ and $x^{t}_{ub}$ are constant vectors, with each element
being the element-wise lower and upper bound on the state.

To achieve the control objective, we propose a robust ZMPC with guaranteed
closed-loop stability in the Lyapunov sense, which is introduced in Section 3.
The proposed robust ZMPC is an extension of the generalized ZMPC proposed in
[21], which is introduced first in Section 2.3 for the completeness.

###  2.3 A generalized ZMPC formulation

The generalized ZMPC proposed in [21] is presented in this section. This
generalized ZMPC assumes perfect knowledge regarding the system dynamics (i.e.
$x(n+1)=f(x(n),u(n),0)$). Throughout this work, we refer to this generalized
ZMPC as the nominal ZMPC, which is used as the referencing benchmark for the
performance validation of our proposed ZMPC. The generalized ZMPC formulation
for time instant $n$ is defined as follows:

$\displaystyle V_{N}^{0}(x(n))=$
$\displaystyle\min\limits_{u(0),\cdots,u(N-1)}\sum_{i=0}^{N-1}\ell_{z}(\hat{x}(i))\vspace{2mm}$
(5a) $\displaystyle{\rm s.t.~{}}$
$\displaystyle\hat{x}(i+1)=f(\hat{x}(i),u(i),0),~{}~{}i\in\mathbb{I}_{0}^{N-1}\vspace{2mm}$
(5b) $\displaystyle\hat{x}(0)=x(n)\vspace{2mm}$ (5c)
$\displaystyle\hat{x}(i)\in\mathbb{X},~{}~{}i\in\mathbb{I}_{0}^{N-1}\vspace{2mm}$
(5d) $\displaystyle u(i)\in\mathbb{U},~{}~{}i\in\mathbb{I}_{0}^{N-1}$ (5e)
$\displaystyle\hat{x}(N)\in\mathbb{X}_{f}$ (5f)

where (5b) is the nominal process model constraint, (5c) is the initial
condition constraint, and (5d) and (5e) are the state and input constraints,
respectively. (5f) denotes the terminal constraint, where $\mathbb{X}_{f}$
denotes the terminal set. $N$ is a positive constant that represents the
control horizon of the controller. $\hat{x}(i)$ denotes the predicted states
over the prediction horizon for $i\in\mathbb{I}_{0}^{N-1}$. $\ell_{z}(\cdot)$
represents the zone-tracking objective and is defined as follows:

$\displaystyle\ell_{z}(z)=\min\limits_{z_{z}}\>$ $\displaystyle
c_{1}||z-z_{z}||_{1}+c_{2}||z-z_{z}||_{2}^{2}$ (6a) $\displaystyle s.t.\>$
$\displaystyle z_{z}\in\mathbb{X}_{t}$ (6b)

where $c_{1}$ and $c_{2}$ are non-negative weighting factors,
$\mathbb{X}_{t}\subset\mathbb{X}$ denotes the target zone, and $z_{z}$ is a
slack variable that is forced to stay inside the target zone by (6b). The
zone-tracking objective is a weighted summation of the 1-norm and squared
2-norm of the minimum difference between the actual state and the slack
variable, which is a representation of the distance between the system state
$x$ and the target set $\mathbb{X}_{t}$. When the system state is outside the
target set, $\ell_{z}$ is positive. When the system state converges to the
target set, $\ell_{z}$ equals to zero.

Taking the definition of (6) into consideration, the objective of the
controller (5) is to find the optimal input sequence for $N$ future steps such
that the system state is driven towards and kept inside the target zone while
constraints (5b) - (5e) are satisfied. In addition, the terminal state $x_{N}$
is forced to converge to
$\mathbb{X}_{f}\subseteq\mathbb{X}_{t}^{M}\subseteq\mathbb{X}_{t}$, where
$\mathbb{X}_{t}^{M}$ denotes the largest forward control invariant set (CIS)
inside $\mathbb{X}_{t}$. The definition of a forward CIS is provided as
follows:

###### Definition 1

(Forward control invariant set [22]) A set $\mathbb{X}_{r}\subseteq\mathbb{X}$
is said to be a forward control invariant set (CIS) for the system
$x(n+1)=f(x(n),u(n),0)$ if there exists a control policy $u(n)=\mu(x(n))$ for
every $x(n)\in\mathbb{X}_{r}$ such that
$x(n+1)=f(x(n),\mu(x(n)),0)\in\mathbb{X}_{r}$.

We refer to the forward CIS as CIS hereafter for simplicity. The nominal ZMPC
is proved to be closed-loop stable in the Laypunov sense when the system model
is exactly known. In this work, we extend the design of the nominal ZMPC to
consider nonlinear systems with model-plant mismatch in the form of process
disturbance as shown in (1). A robust ZMPC with guaranteed closed-loop
stability is proposed in the following section.

##  3 The proposed robust ZMPC

The proposed robust ZMPC is presented in this section. The effect of process
uncertainty is encountered through modifying the tracking target set and the
terminal set. The proposed robust ZMPC is designed based on the nominal system
$\tilde{x}(n+1)=f(x(n),u(n),0)$, as the value of the disturbance $w(n)$ is
unknown at any given time $n\in\mathbb{I}_{\geq 0}$. The proposed ZMPC design
for time instant $n$ is presented as follows:

$\displaystyle V_{N}^{0}(x(n))=$
$\displaystyle\min\limits_{u(0),\cdots,u(N-1)}\sum_{i=0}^{N-1}\tilde{\ell}_{z}(\tilde{x}(i))\vspace{2mm}$
(7a) $\displaystyle{\rm s.t.~{}}$
$\displaystyle\tilde{x}(i+1)=f(\tilde{x}(i),u(i),0),~{}~{}i\in\mathbb{I}_{0}^{N-1}\vspace{2mm}$
(7b) $\displaystyle\tilde{x}(0)=x(n)\vspace{2mm}$ (7c)
$\displaystyle\tilde{x}(i)\in\mathbb{X},~{}~{}i\in\mathbb{I}_{0}^{N-1}\vspace{2mm}$
(7d) $\displaystyle
u(i)\in\mathbb{U},~{}~{}i\in\mathbb{I}_{0}^{N-1}\vspace{2mm}$ (7e)
$\displaystyle\tilde{x}(N)\in\tilde{\mathbb{X}}_{f}$ (7f)

where $\tilde{\ell}_{z}(\cdot)$ is defined as follows:

$\displaystyle\tilde{\ell}_{z}(z)=\min\limits_{z_{z}}\>$ $\displaystyle
c_{1}||z-z_{z}||_{1}+c_{2}||z-z_{z}||_{2}^{2}\vspace{2mm}$ (8a) $\displaystyle
s.t.\>$ $\displaystyle z_{z}\in\tilde{\mathbb{X}}_{t}$ (8b)

Similar notations employed in (5) and (6) are adopted in the proposed
controller (7) and (8) except that $\tilde{x}$ is used to represent the
nominal state predicted by the nominal model. Compare to the nominal ZMPC, the
key changes took place in the proposed design are the modified target set
$\tilde{\mathbb{X}}_{t}$ and the terminal set $\tilde{\mathbb{X}}_{f}$ in (8b)
and (7f), respectively. The modified target set is a subset of the original
target set ($\tilde{\mathbb{X}}_{t}\subset\mathbb{X}_{t}$). The terminal
constraint is needed to ensure the closed-loop stability with a CIS inside the
modified target zone $\tilde{\mathbb{X}}_{t}$.

The proposed ZMPC is applied to system (1) in the receding horizon fashion. At
each time instant $n\in\mathbb{I}_{\leq 0}$, the optimal input
$u(i|n)^{*},~{}i\in\mathbb{I}_{0}^{N-1}$ for $N$ future steps is obtained by
solving the optimization problem (7). Only the optimal input for $i=0$ is
applied to the system at time $n$. At time $n+1$, the optimization problem (7)
is evaluated again with updated initial condition.

Figure 1: Potential nominal state trajectory (dotted line) and actual state
trajectory (dashed line) of the closed-loop system under the proposed ZMPC.

Figure 1 presents potential trajectories of the closed-loop nominal state and
the actual state based on the proposed robust ZMPC. The solid box represents
the state space $\mathbb{X}$, while the rectangles bounded by the dashed and
the dotted lines represent the actual target set $\mathbb{X}_{t}$ and the
modified target set $\tilde{\mathbb{X}}_{t}$, respectively. The ellipsoid
bounded by the dotted line denotes the largest CIS
$\tilde{\mathbb{X}}_{t}^{M}$ inside the modified target set, which is used as
the terminal set $\tilde{\mathbb{X}}_{f}$ of the proposed controller.
Oftentimes, a nominal ZMPC tends to drive the system state to the boundary of
the target set, which likely leads to target set violation in the presence of
process uncertainties. By ensuring the closed-loop convergence of the nominal
state to the modified target set, which is a subset of the actual target set,
a buffer zone is provided to reject the effect of process uncertainties. The
terminal constraint (7f) ensures that the nominal state will stay inside the
control invariant terminal set $\tilde{\mathbb{X}}_{f}$ and thus
$\tilde{\mathbb{X}}_{t}$ once it has first converged. More details regarding
the closed-loop stability of the proposed controller are discussed in the
following section.

###### Remark 1

Throughout the manuscript, we use $x(\cdot)$ to denote the true system state
at a given time step. Contrarily, $\tilde{x}(\cdot)$ denotes the system state
predicted based on the knowledge available to the controller. Following
similar ideas, the set $\mathbb{X}_{t}$ reflects the true zone control
objectives we have on the system of interest, which is referred to as the
actual target set. For example, the objective is to keep the reactor
temperature in a certain range. The set $\tilde{\mathbb{X}}_{t}$ in
comparison, is mathematically determined based on the actual target set with
no physical significance and is referred to as the modified target set.
Furthermore, we define the set $\mathbb{X}_{t}^{M}$ and
$\tilde{\mathbb{X}}_{t}^{M}$ as the largest CIS inside $\mathbb{X}_{t}$ and
$\tilde{\mathbb{X}}_{t}$, respectively.

##  4 Stability analysis

In this section, we investigate the closed-loop stability of the proposed
robust ZMPC (7). First, we introduce the definition of the $N$-step reachable
set, which is essential in characterizing the feasibility of the proposed
robust ZMPC.

###### Definition 2

($N$-step reachable set [21]) We use
$\mathbb{X}_{N}(\mathbb{X},\mathbb{X}_{f})$ to denote the set of states in
$\mathbb{X}$ that can be steered to $\mathbb{X}_{f}$ in $N$ steps while
satisfying the state and input constraints $(x,u)\in\mathbb{Z}$. That is,

$\mathbb{X}_{N}(\mathbb{X},\mathbb{X}_{f})=\big{\\{}x(0)\in\mathbb{X}\;|\;\exists\>(x(n),u(n))\in\mathbb{Z},n\in\mathbb{I}_{0}^{N-1},x(N)\in\mathbb{X}_{f}\big{\\}}$

The following assumption is employed to ensure the feasibility of the proposed
robust ZMPC.

###### Assumption 4

The terminal set $\tilde{\mathbb{X}}_{f}$ and the $N$-step reachable set
$\mathbb{X}_{N}(\mathbb{X},\tilde{\mathbb{X}}_{f})$ are compact and nonempty.

For the completion of this work, here we summarize the essential results
obtained in [21] regarding the nominal ZMPC as Theorem 1.

###### Theorem 1

(Stability of the nominal ZMPC [21]) If Assumptions 1 and 4 hold,
$x_{0}\in\mathbb{X}_{N}(\mathbb{X},\mathbb{X}_{f})$, $w(n)=0$, and in
addition, the optimal value function $V_{N}^{0}(x(n))$ is locally continuous
on $\mathbb{X}_{t}$, then: (i) the zone MPC is recursively feasible with
$x(n)\rightarrow\mathbb{X}_{t}^{M}$, (ii) the terminal set
$\mathbb{X}_{t}^{M}$ is asymptotically stable, (iii) the following inequality
holds:

$V_{N}^{0}(\tilde{x}(n+1))-V_{N}^{0}(\tilde{x}(n))\leq-\ell_{z}(\tilde{x}(n)),\;\;n\in\mathbb{I}_{\geq
0}$ (9)

where $V_{N}^{0}$ is the optimum control objective function and is a Lyapunov
function with respect to the set $\mathbb{X}_{t}^{M}$, and $\ell_{z}(\cdot)$
denotes the zone-tracking objective function defined in (6).

As the proposed ZMPC is applied to the system of interest in the receding
horizon fashion, the optimization problem is solved at each time step with the
updated information regarding the actual state at the initial time step over
the prediction horizon (i.e. $\tilde{x}(n)=x(n)$). This implies only the one-
step-ahead deviation caused by the plant-model mismatch between the nominal
model and the actual model needs to be considered. The following Proposition
provides an upper bound on the one-step-ahead deviation of the predicted state
from the actual state:

###### Proposition 1

(c.f. [20]) Consider the actual system (1) and the nominal system with
$w(n)=0$. If Assumption 2 is satisfied, starting from a known initial
condition $x(n)$, the deviation of the actual state $x(n+1)$ from the
predicted state $\tilde{x}(n+1)$ in one time step is bounded:

$|x(n+1)-\tilde{x}(n+1)|\leq\sqrt{n_{x}}L_{w}\theta,~{}\forall
x(n+1),~{}\tilde{x}(n+1)\in\mathbb{X}$ (10)

Based on the upper bound of the one-step-ahead deviation in the state,
Proposition 2 provides an upper bound on the progression of the Lyapunov
function $V_{N}^{0}(\cdot)$ in one time step from $x(n)$ to $x(n+1)$.

###### Proposition 2

Consider the Lyapunov function $V_{N}^{0}(x)$. There exist positive constants
$K_{V}$ and $H$ such that:

$V_{N}^{0}(x(n+1))-V_{N}^{0}(x(n))\leq-\tilde{\ell}_{z}(\tilde{x}(n))+f_{V}(\sqrt{n_{x}}L_{w}\theta)$
(11)

where $n_{x}$ and $\theta$ are known parameters, $\tilde{\ell}_{z}(\cdot)$
denotes the zone-tracking objective function defined in (8), and $f_{V}$ is a
function defined as follows:

$f_{V}(x)=K_{V}x+Hx^{2}$

where $K_{V}$ and $H$ are positive constants.

Proof. Upon applying Taylor expansion to $V_{N}^{0}(x)$ around the predicted
state $\tilde{x}(n+1)$, the following relation can be obtained:

$\begin{split}&V_{N}^{0}(x(n+1))=V_{N}^{0}(\tilde{x}(n+1))+\\\
&\left.\frac{\partial V_{N}^{0}}{\partial
x}\right|_{\tilde{x}(n+1)}|x(n+1)-\tilde{x}(n+1)|+H.O.T.\end{split}$ (12)

where $H.O.T.$ includes the higher order terms of the Taylor expansion. A
positive constant $H$ can be found for $x\in\mathbb{X}$ such that the
following holds:

$H.O.T.\leq H|x(n+1)-\tilde{x}(n+1)|^{2}$ (13)

Consider (12) and (13) correspondingly, we can obtain the following inequality
regarding the difference in $V_{N}^{0}$ due to the effect of disturbance in
one step:

$\begin{split}V_{N}^{0}(x&(n+1))\leq V_{N}^{0}(\tilde{x}(n+1))\\\
&+\left.\frac{\partial V_{N}^{0}}{\partial
x}\right|_{\tilde{x}(n+1)}|x(n+1)-\tilde{x}(n+1)|\\\
&+H|x(n+1)-\tilde{x}(n+1)|^{2}\end{split}$ (14)

If Assumptions 3 and 4 are satisfied, there exists a positive constant $K_{V}$
such that the magnitude of the partial derivative $\frac{\partial
V_{N}^{0}}{\partial x}$ is bounded by it:

$\left|\frac{\partial V_{N}^{0}}{\partial x}\right|\leq K_{V}$ (15)

Define $f_{V}(x)=K_{V}x+Hx^{2}$. Based on (14) and (15), applying Proposition
1, the following relationship can be derived:

$V_{N}^{0}(x(n+1))\leq
V_{N}^{0}(\tilde{x}(n+1))+f_{V}(\sqrt{n_{x}}L_{w}\theta)$ (16)

Taking into account of (9) and (16), plus the fact that $x(n)=\tilde{x}(n)$,
the relationship in (11) can be obtained and this proves Proposition 2.
$\blacksquare$

To derive the amount of shrinkage required in the target set that ensures the
closed-loop convergence of the system in the presence of disturbance, we
introduce the following definition of Lyapunov function level set:

###### Definition 3

(Level set of Lyapunov function) The set $\Omega_{\rho}$ is defined to be the
level set of the Lyapunov function $V_{N}^{0}(\cdot)$:

$\Omega_{\rho}:=\\{x\in\mathbb{X}:V_{N}^{0}(x)\leq\rho\\}$ (17)

Let $\Omega_{\rho_{max}}$ be the largest level set of $V_{N}^{0}$ inside
$\mathbb{X}_{t}^{M}$:

$\rho_{max}:=max\\{\rho:x\in\mathbb{X}_{t}^{M},V_{N}^{0}(x)\leq\rho\\}$ (18)

Define $\Omega_{\rho_{e}}\subset\Omega_{\rho_{max}}$ that satisfies:

$\rho_{e}=\rho_{max}-f_{V}(\sqrt{n_{x}}L_{w}\theta)$ (19)

###### Assumption 5

$\Omega_{\rho_{e}}$ is nonempty.

###### Theorem 2

Consider system (1) under the control of the proposed ZMPC (7). If
$x(n)\in\Omega_{\rho_{e}}$ and Assumptions 3 \- 5 are satisfied, then it is
guaranteed that $x(n+1)\in\Omega_{\rho_{max}}$.

Proof. At a given time instant $n$, with the initial condition
$x(n)\in\Omega_{\rho_{e}}$, it is guaranteed that there exists a feasible
control policy such that the nominal state at the next time step
$\tilde{x}(n+1)\in\Omega_{\rho_{e}}$, as $\Omega_{\rho_{e}}$ is control
invariant. Take the inequality expression (16) into consideration, the
deviation in the value of Lyapunov function $V^{0}_{N}$ from $\tilde{x}(n+1)$
to $x(n+1)$ is bounded by $f_{V}(\sqrt{n_{x}}L_{w}\theta)$. Recall that
$\Omega_{\rho_{e}}$ is a shrunk Lyapunov function level set with respect to
$\Omega_{\rho_{max}}$, with the shrinkage defined in the Lyapunov function
value by $f_{V}(\sqrt{n_{x}}L_{w}\theta)$, which proves Theorem 2.
$\blacksquare$

Theorem 2 indicates the convergence of the actual state into a control
invariant subset of the actual target set $\mathbb{X}_{t}$ if the initial
state $x(n)\in\Omega_{\rho_{e}}$. The following theorem provides statements on
the upper bound of the shrinkage amount in the target set required in the
proposed ZMPC design such that the state is guaranteed to be driven into
$\Omega_{\rho_{e}}$.

###### Definition 4

(The smallest polyhedron set around a known set) We define the smallest
polyhedron set bounded by orthogonal hyperplanes $\mathbb{P}(\mathbb{S})$
around a given set $\mathbb{S}\subseteq\mathbb{R}^{n_{s}}$ as follows:

$\mathbb{P}(\mathbb{S})=\\{s|s_{lb}\leq s\leq s_{ub}\\}$ (20)

where $s_{lb},s_{ub}\in\mathbb{R}^{n_{s}}$ are vectors with each element
$s_{lb}^{i},~{}i=1,2,\cdots,{n_{s}}$ and $s_{ub}^{i},~{}i=1,2,\cdots,{n_{s}}$
being the lower and upper bound of each dimension of the set $\mathbb{S}$:

$\displaystyle s_{lb}^{i}$
$\displaystyle=\min\\{s^{i}|s\in\mathbb{S}\\},~{}i=1,2,\cdots,{n_{s}}$ (21a)
$\displaystyle s_{ub}^{i}$
$\displaystyle=\max\\{s^{i}|s\in\mathbb{S}\\},~{}i=1,2,\cdots,{n_{s}}$ (21b)

where $s^{i}$ is the $i$-th element of the vector $s$.

###### Theorem 3

Consider system (1) under the control of the proposed ZMPC (7). If Assumptions
1 \- 5 are satisfied,
$x(0)\in\mathbb{X}_{N}(\mathbb{X},\tilde{\mathbb{X}}_{f})$, and the following
relationship holds:

$-\tilde{\ell}_{z}(x)+f_{V}(\sqrt{n_{x}}L_{w}\theta)<0,~{}\forall
x\in\mathbb{X}_{N}(\mathbb{X},\tilde{\mathbb{X}}_{f})\backslash\Omega_{\rho_{e}}$
(22)

then there exists
$\overline{\varepsilon},\underline{\varepsilon}\in\mathbb{R}^{n_{x}}$ with
non-negative elements that define the modified target zone
$\tilde{\mathbb{X}}_{t}$ as follows:

$\tilde{\mathbb{X}}_{t}=\\{x:\tilde{x}_{lb}^{t}\leq
x\leq\tilde{x}_{ub}^{t}|\tilde{x}_{lb}^{t}-x_{lb}^{t}\leq\underline{\varepsilon},x_{ub}^{t}-\tilde{x}_{ub}^{t}\leq\overline{\varepsilon}\\}$
(23)

such that the proposed ZMPC controller is recursively feasible and guarantees
the closed-loop convergence of the system state into the actual target zone
$\mathbb{X}_{t}$. Note that $x_{lb}^{t}$ and $x_{ub}^{t}$ are known vectors
containing the element-wise lower and upper bound of the actual target set
$\mathbb{X}_{t}$ as defined in Assumption 1.

Proof. We define the modified target set as the smallest polyhedron set around
$\Omega_{\rho_{e}}$ (i.e.
$\tilde{\mathbb{X}}_{t}=\mathbb{P}(\Omega_{\rho_{e}})$) and the corresponding
terminal set $\tilde{\mathbb{X}}_{f}$ as the largest CIS inside
$\tilde{\mathbb{X}}_{t}$. If Assumption 5 holds, it means there exists a valid
set of $\tilde{x}_{lb}^{t}$ and $\tilde{x}_{ub}^{t}$ such that
$\tilde{\mathbb{X}}_{t}\subset\mathbb{X}_{t}$ is bounded and nonempty. Recall
that the actual target set $\mathbb{X}_{t}$ is bounded, meaning that the
elements in vectors
$\overline{\varepsilon},\underline{\varepsilon}\in\mathbb{R}^{n_{x}}$ are
bounded and non-negative.

Take the relationship defined by (22) into consideration. The value of the
Lyapunov function $V_{N}^{0}$ is guaranteed to be decreasing until the state
converges to $\Omega_{\rho_{e}}$, thus ensuring the convergence of the nominal
state $\tilde{x}$ into $\Omega_{\rho_{e}}$.

In the presence of disturbance, the actual state $x$ may deviate from
$\Omega_{\rho_{e}}$ after the first convergence of the nominal state
$\tilde{x}$. Taking Theorem 2 into consideration, it is guaranteed that with
$x(n)\in\Omega_{\rho_{e}}$, $x(n+1)\in\Omega_{\rho_{max}}$, meaning the actual
state is kept inside
$\Omega_{\rho_{max}}\subset\mathbb{X}_{t}^{M}\subset\mathbb{X}_{t}$ in one
step. As
$\Omega_{\rho_{e}}\subset\Omega_{\rho_{max}}\subset\mathbb{X}_{N}(\mathbb{X},\tilde{\mathbb{X}}_{f})$,
the actual state that has deviated in one step due to the effect of
disturbance is guaranteed to be driven back to $\Omega_{\rho_{e}}$ in finite
steps as (22) holds. This indicates that the proposed ZMPC controller is able
to drive the system state into the actual target set $\mathbb{X}_{t}$ and
keeps it inside. $\blacksquare$

###### Remark 2

By definition, any CIS inside the modified target set $\tilde{\mathbb{X}}_{t}$
can be used as the terminal set $\tilde{\mathbb{X}}_{f}$. However, in order to
gain a controller with a larger feasibility region and more degrees of
freedom, we define $\tilde{\mathbb{X}}_{f}=\tilde{\mathbb{X}}_{t}^{M}$. The
largest CIS inside the modified target set is approximated using the algorithm
proposed in [23]. The algorithm provides an inner approximation of the largest
CIS in the form of a polyhedron bounded by hyperplanes.

##  5 An algorithm for estimating the modified target set
$\tilde{\mathbb{X}}_{t}$

The stability proof in the previous section ensures the existence of an upper
bound of the shrinkage amount on the actual target set, however, it is very
challenging to obtain explicit expressions of the Lyapunov function and
numerical values of the parameters. In order to make the proposed controller
more practical in applications, we propose an algorithm for estimating and
tuning the modified target set.

The proposed algorithm is inspired by a study on the problem in the context of
linear systems. Consider a linear system as follows:

$x(n+1)=Ax(n)+Bu(n)+Ew(n)$ (24)

where $A$, $B$, and $E$ are matrices of the appropriate dimensions. Let us
assume that the linear system is controllable and there exists the following
quadratic Lyapunov function for the system:

$V^{0}_{N}(x)=x^{T}Px$ (25)

where $P$ is a symmetric positive definite matrix. Recall the Lipschitz
continuity Assumption 3. If we substitute the linear system model (24) into
the left hand side of the inequality (3), the following expression can be
obtained:

$|f(x,u,w)-f(x,u,0)|=|Ew|\leq L_{w}|w|$ (26)

By the property of norm, the following expression holds:

$|Ew|\leq|E||w|$ (27)

Thus, we may define $L_{w}=|E|$ such that Assumption 3 holds. Recall the
inequality relation (10):

$|x(n+1)-\tilde{x}(n+1)|\leq\sqrt{n_{x}}L_{w}\theta$

where $n_{x}$ is the dimension of the state, $\theta$ is the upper bound of
the one-norm of the disturbance and is assumed to be known. Thus for the
linear system, the right-hand side of the inequality can be determined
explicitly. By calculating the term $\sqrt{n_{x}}L_{w}\theta$, a conservative
approximation of the norm of the one-step-ahead difference between the actual
state $x(n+1)$ and the nominal state $\tilde{x}(n+1)$ in the presence of
disturbance is obtained. For the nonlinear system of interest, the matrix $E$
is essentially the sensitivity matrix of the state with respect to the
disturbance (i.e. $\frac{\partial x}{\partial w}$). The value of the
sensitivity matrix is a function of the state, input, and disturbance and can
be calculated explicitly with a given combination of the variables.

Practically, based on inequality (10), we propose the following definition of
the shrinkage amount $s$ in the state boundaries of the target zone for the
nonlinear system:

$s=\gamma x^{d}_{max}$ (28)

where $\gamma\in\mathbb{R}^{+}$ is a risk factor that helps to tune how
conservative the proposed controller is, and $x^{d}_{max}$ essentially
estimates the largest potential effect the disturbance can have on the actual
state compared to the nominal state in one-step-ahead prediction, of which the
definition is inspired by (10). $x^{d}_{max}$ can be calculated using
Algorithm 1 based on sensitivity analysis.

Algorithm 1 Estimation of the state deviation due to the presence of
disturbance

$x\in\mathbb{X}^{M}_{t}$, $u\in\mathbb{U}$, $w\in\mathbb{W}$

$I_{t}\leftarrow[\>]$

for $i=0:n_{x}-1$ do

if $x_{lb}<x^{t}_{lb}$ or $x_{ub}>x^{t}_{ub}$ then

$I_{t}\leftarrow[I_{t},1]$

else $I_{t}\leftarrow[I_{t},0]$

end if

end for

Define $z\in\mathbb{R}^{n_{z}}\leftarrow[x,u]$, $n_{z}=n_{x}+n_{u}$

Determine the sensitivity function $\frac{\partial x}{\partial
w}:\mathbb{R}^{n_{z}+n_{w}}\rightarrow\mathbb{R}^{n_{x}\times n_{w}}$

Find and store the maximum and minimum of each variables
$z_{max}\in\mathbb{R}^{n_{z}}\leftarrow\max{\\{z\\}},\>z_{min}\in\mathbb{R}^{n_{z}}\leftarrow\min{\\{z\\}}$

$w_{max}\in\mathbb{R}^{n_{w}}\leftarrow\max{\\{w\\}},\>w_{min}\in\mathbb{R}^{n_{w}}\leftarrow\min{\\{w\\}}$

Calculate the set of the accumulated sensitivity $\mathbb{X}^{d}$ for the
state variables with zone-tracking objective at each extreme point of $z$ and
$w$:  $\mathbb{X}^{d}:=\left\\{\left[\frac{\partial x}{\partial
w}\Big{|}_{z^{i},w^{j}}\cdot w\right]\odot
I_{t}\;\Big{|}z^{i}\in[z_{max}^{i},z_{min}^{i}],i\in\mathbb{I}_{0}^{n_{z}-1},w^{j}\in[w_{max}^{j},w_{min}^{j}],j\in\mathbb{I}_{0}^{n_{w}-1}\right\\}$

Find the element $x^{d}_{max}$ with the largest norm in $\mathbb{X}^{d}$
$x^{d}_{max}\leftarrow\left\\{x^{d}_{max}\;\Big{|}x^{d}_{max}\in\mathbb{X}^{d},|x^{d}_{max}|\geq|x^{d}|,\forall
x^{d}\in\mathbb{X}^{d}\right\\}$

return $x^{d}_{max}$

In Algorithm 1, the set $\mathbb{X}^{d}$ contains $2^{n_{z}+n_{w}}$ number of
vectors. Each vector represents the amount of accumulated effects of the
disturbance $w$ have on the state $x$ at one of the extremes of the joint
space $\mathbb{X}^{M}_{t}\times\mathbb{U}\times\mathbb{W}$. It is worth
pointing out that the sensitivity matrix is computed for states inside the
maximum CIS $\mathbb{X}^{M}_{t}$, that is, inside the actual target set. This
is because it has been proved in [21] that for the nominal system, guaranteed
convergence into the actual target set $\mathbb{X}_{t}$ is equivalent to the
convergence into its maximum control invariant subset $\mathbb{X}^{M}_{t}$.
Since the purpose of modifying the target set is to avoid zone violation once
it has converged, we focus on rejecting the effect of disturbance when
$x\in\mathbb{X}^{M}_{t}$. Furthermore, the effect of disturbance is only
investigated for the individual state variables that have a zone-tracking
objective on them, as it is common that the zone-tracking objective is not
enforced on all state variables.

The modified target zone can be calculated based on the shrinkage amount $s$
as follows:

$\tilde{\mathbb{X}}_{t}:=\left\\{x\Big{|}x^{t}_{lb}+s\leq x\leq
x^{t}_{ub}-s\right\\}$ (29)

###### Remark 3

The inequality relationship (10) is defined in terms of the norm of vectors
and matrix. However, in practice, the shrinkage amount $s$ should be a vector
that has the same dimension of the state, as the shrinkage required for each
state needs to be determined individually. Thus, $x^{d}_{max}$ is selected to
be the accumulated effect of $w$ on $x$ that has the largest norm, which
mimics the idea of (10), however in the state space.

###### Remark 4

In Algorithm 1, the value of the sensitivity matrix needs to be calculated
$2^{n_{z}+n_{w}}$ times, as every combination of the maximum and minimum of
the state, input, and disturbance are considered. For higher-dimension
systems, the computational effort may be significant. One way to reduce the
computational cost is to utilize available physical knowledge in the
implementation. For example, the optimal operating condition may be estimated
based on the nominal model, or it is known to be on one edge of the zone. In
these cases, the sensitivity matrix may be calculated only at the estimated
optimal operating point, or at the extremes on one edge.

Furthermore, (29) assumes that the target set shrinks in the same amount on
the upper and lower bound. If the optimal operating point is known to be on
one edge of the target zone, then the shrinkage can be applied to that edge
only. This again helps to preserve a larger modified target zone and terminal
zone, which helps to reduce the economic performance loss and feasibility
issue.

##  6 Consideration of economic objectives

A natural extension to the proposed ZMPC (7) is to include an economic
objective, as the proposed design provides the system with more degrees of
freedom to optimize additional objectives. The proposed ZMPC design that
handles additional economic objective is presented as follows:

$\begin{split}\min\limits_{u_{0},\cdots,u_{N-1}}&\sum_{i=0}^{N-1}\tilde{\ell}_{z}(\tilde{x}_{i})+\ell_{e}(\tilde{x}_{i},u_{i})\vspace{2mm}\\\
{\rm s.t.~{}}&\text{(\ref{eqn:ZMPC2b})}-\text{(\ref{eqn:ZMPC2f})}\end{split}$
(30)

With weighting factors being large enough on the zone control objective,
namely $c_{1}$ in (8), the controller will prioritize the zone-tracking
objective over the economic objective such that the economic objective will be
optimized only after the system enters the target zone [21]. Thus the
stability properties of the controller (30) remain unchanged from that of (7).
In the following section, simulation results based on (30) are presented.

###### Remark 5

By tuning the parameter $\gamma$ in (28), a trade-off between the size of the
initial feasible region, the economic performance, and zone-tracking
performance can be achieved due to changing target zone shrinkage $s$. A
larger $\gamma$ leads to a smaller modified target zone and correspondingly, a
smaller terminal zone. This enhances convergence to the actual target zone
while potentially reducing the flexibility and feasibility of the controller,
as the system is forced to converge into a smaller region at the end of the
control horizon. Higher sacrifice in the economic performance may be observed
as well. On the other hand, smaller $\gamma$ poses less effect on the
controller feasibility and economic performance but the closed-loop system has
a greater chance to violate the actual target zone.

##  7 Simulations

In this section, we use a continuous stirred tank reactor (CSTR) as the
benchmark process to demonstrate the effectiveness of the proposed design.
After introducing the CSTR process (Section 7.1) and the detailed setting of
the proposed controller (Section 7.2), the performance of the proposed
controller is first validated without and with a secondary economic objective
with respect to the nominal ZMPC in Section 7.3. The importance of shrinking
the target set in the proposed controller design is verified in Section 7.4.
The significance of the terminal constraint in ensuring stability is
investigated in Section 7.5. Finally, the impact of the tuning parameter
$\gamma$ is investigated in Section 7.6.

###  7.1 Process description

The benchmark CSTR process is introduced in this section. An irreversible
first-order exothermic reaction $A\rightarrow B$ takes place inside the CSTR
reactor, where A and B are the reactants and the desired product,
respectively. A cooling jacket is used to control the temperature of the CSTR.
Assuming perfect mixing and constant mixture volume inside the reactor, the
system can be described by the following ordinary differential equations:

$\displaystyle\frac{dC_{A}}{dt}=\frac{q}{V}(C_{A_{f}}-C_{A})-k_{0}\exp{\Big{(}-\frac{E}{RT}\Big{)}}C_{A}$
(31a) $\displaystyle\frac{dT}{dt}=\frac{q}{V}(T_{f}-T)+\frac{UA}{V\rho
C_{p}}(T_{c}-T)+\frac{-\triangle H}{\rho
C_{p}}k_{0}\exp{\Big{(}-\frac{E}{RT}\Big{)}}C_{A}$ (31b)

where $C_{A}$ $[mol/L]$ denotes the molar concentration of the reactant, $T$
$[K]$ represents the temperature inside the reactor, $T_{c}$ $[K]$ represents
the temperature of the cooling stream, and $C_{A_{f}}$ $[mol/L]$ and $T_{f}$
$[K]$ denote the molar concentration of the reactant and the temperature of
the feed stream, respectively. $q$ $[L/min]$ represents the volumetric flow
rate of the streams entering and leaving the reactor, $V$ $[L]$ and $\rho$
$[g/L]$ denote the volume and the density of the mixture inside the reactor
respectively, and $k_{0}$, $E$ and $R$ are reaction-related parameters, namely
the pre-exponential factor of the reaction rate, the activation energy, and
the universal gas constant. $C_{p}$ denotes the specific heat capacity of the
mixture, and $\triangle H$ and $UA$ denote the heat of reaction and the heat
transfer coefficient between the reactor and the cooling jacket, respectively.
The values of the parameters are adopted from [20].

Recall the system defined in (1). For the CSTR example specifically, we define
the system state vector to be $x=[C_{A},T]^{T}$ and the control input to be
$u=T_{c}$. The disturbance vector is defined as
$w=[C_{A_{f}}-\bar{C}_{A_{f}},T_{f}-\bar{T}_{f}]^{T}$, where
$\bar{C}_{A_{f}}=1.0\,mol/L$ and $\bar{T}_{f}=350\,K$ are the nominal values
of the corresponding variables. The following bounds are employed for the
state, input, and disturbance, respectively:

$\displaystyle[0.0,345.0]^{T}$ $\displaystyle\leq x\leq[1.0,355.0]^{T}$ (32)
$\displaystyle 285.0$ $\displaystyle\leq u\leq 315.0$ (33)
$\displaystyle[-0.1,-2.0]^{T}$ $\displaystyle\leq w\leq[0.1,2.0]^{T}$ (34)

###  7.2 Controller setup and implementation

The proposed ZMPC controller is implemented on the CSTR system without and
with a secondary economic objective. Following the design proposed in (7) and
(30), we define the zone control objective and economic objective. The zone
control objective is to hold the reactor temperature $T$ between $348.0\,K$ to
$352.0\,K$. Recall the definition of $\ell_{z}$ in (6), which implies the
following:

$\mathbb{X}_{t}=\\{x:[0.0,348.0]^{T}\leq x\leq[1.0,352.0]^{T}\\}$

The largest CIS inside the original and modified target zones are approximated
using the algorithm proposed in [23].

We computed the value of $x_{max}^{d}$ to be 0.511 for the states inside the
largest CIS with the input and disturbances in their stated boundaries. It was
found that the disturbance has the strongest effect on the state at the
following point:

$x=[0.754,352.0],\>u=315.0,\>w=[0.1,2]$

Unless otherwise mentioned, $\gamma=1$ is used, which means
$s=x^{d}_{max}=0.511$. The economic objective is to minimize the concentration
of the reactant inside the reactor, i.e. $\ell_{e}=C_{A}$. For all the
simulations, the control horizon $N=5$ unless otherwise mentioned.

The value of the disturbance $w(n)$ at any time step is selected randomly
inside the boundaries specified by (34). Except for the CIS approximation that
is carried out in Julia, all other simulations are performed in Python. The
optimization problems are solved using the CasADi toolbox developed by [24].

###  7.3 Control performance validation

Figure 2: Proposed ZMPC without economic objectives. Figure 3: Proposed ZMPC
v.s. the nominal ZMPC with economic objective, $x_{0}=[0.12,355]$. Figure 4:
Proposed ZMPC v.s. the nominal ZMPC with economic objective,
$x_{0}=[0.9,345]$.

In this section, we validate the performance of the proposed controller.
Figure 2 presents the state trajectories under the proposed ZMPC without any
economic objective (7) with and without the presence of disturbance. Figures 3
and 4 present the performance of the proposed ZMPC with economic objective
(30) in comparison to the nominal ZMPC based on different initial conditions.
In all figures, the solid black box represents the overall state space, and
the dashed blue and red lines denote the boundaries of the original target set
$\mathbb{X}_{t}$ and the modified target set $\tilde{\mathbb{X}}_{t}$
respectively. The shaded areas represent the largest CIS inside
$\mathbb{X}_{t}$ and $\tilde{\mathbb{X}}_{t}$. The larger shaded area
represents $\mathbb{X}_{t}^{M}$, while the smaller one is
$\tilde{\mathbb{X}}_{t}^{M}$.

It can be observed from the blue trajectory in Figure 2 that without any
economic objective, the system converges to a steady-state in the center of
the target zone. In the presence of disturbance, as represented by the red
trajectory, the state oscillates around the optimal steady state but is able
to stay inside the terminal CIS once entered. Thus, for this particular setup,
the benefit of the proposed controller is not significant.

To validate the benefits of the proposed approach, the secondary economic
objective is added to the system. With the economic objective considered, the
optimal operating condition shifts to the boundary of the target zone. In
Figures 3 and 4, starting from different initial conditions, the state
progressions under the control of the nominal ZMPC are shown in blue, while
the state trajectories controlled by the proposed ZMPC are red. Under the
nominal ZMPC, the state leaves the terminal CIS after converging to the
optimal operating point due to the fluctuation caused by the presence of
disturbance. These violations on the zone-tracking objective are unfavored and
should be eliminated. As presented by the red trajectories, violations are
avoided by the proposed ZMPC. Although the state trajectories violate the
upper bound of the modified target zone $\tilde{\mathbb{X}}_{t}$, they are
still enclosed in the actual target zone $\mathbb{X}_{t}$. Furthermore,
Figures 3 and 4 verify that the proposed controller is able to drive the
system state to the same economically-optimal neighborhood even under
different initial conditions.

###  7.4 The significance of modifying the target set

Figure 5: ZMPC with economic objective that tracks the original target set
($\mathbb{X}_{t}$) with modified terminal set
($\mathbb{X}_{f}=\mathbb{X}_{t}^{M}$).

One natural question that rises is that instead of modifying just the terminal
set, why both the tracking target set and the terminal set are modified in the
proposed approach. Figure 5 investigates the controller performance if the
original target set $\mathbb{X}_{t}$ is tracked and only the terminal set is
modified. The state trajectories under the control of the nominal ZMPC (blue
solid trajectory with hexagon markers), the proposed ZMPC (red dashed
trajectory with dot markers), and the ZMPC that tracks the original target set
$\mathbb{X}_{t}$ and the modified terminal set $\tilde{\mathbb{X}}_{t}^{M}$
(cyan dotted trajectory with cross markers) are presented, respectively. It
can be observed that if only the terminal set is modified, the optimal state
trajectory is identical to that obtained based on the nominal ZMPC. This
implies that the control performance is not affected if only the terminal set
is modified. This matches the theory proposed in [19], which indicates the
controller performance remains equivalent in terms of stability as far as
$\tilde{\mathbb{X}}_{f}\subseteq\mathbb{X}_{t}^{M}$. Note that the state
trajectory obtained based on the proposed ZMPC is presented in the figure for
reference, which verifies the effectiveness of modifying the tracking target
set.

###  7.5 The significance of the terminal constraint

Figure 6: ZMPC without the terminal constraint v.s. the proposed ZMPC with
economic objective.

The significance of the terminal constraint in ensuring closed-loop stability
is investigated in this section. With the control horizon $N=3$, Figure 6
shows the state trajectory under the ZMPC without any terminal constraint in
cyan with cross markers, while the trajectory obtained based on the proposed
controller is presented in red with dot markers for reference. Without the
terminal constraint, the controller is able to drive the system to the optimal
operating point, however with aggressive violations of the hard state
constraint represented by the solid black box. These violations are avoided by
the proposed controller with the terminal constraint. Furthermore, it can be
observed that under the controller without terminal constraint, the system
takes significantly more steps to converge compared to that taken based on the
proposed approach, indicating the effectiveness of the terminal constraint in
ensuring closed-loop stability and enhancing convergence.

###  7.6 The impact of $\gamma$

Figure 7: The modified target sets and the corresponding maximum CISs under
different $\gamma$. Figure 8: The zone-tracking performance and economic
performance under different $\gamma$.

The impact of the tuning parameter $\gamma$ is investigated in this section.
Figure 7 presents the modified target sets and the corresponding maximum CIS
inside the modified target sets with different values of $\gamma$. Figure 8
displays the accumulated zone-tracking cost and the economic cost as $\gamma$
increases. The accumulated zone-tracking costs are presented in blue with
pentagon markers, while the accumulated economic costs are shown in red with
circle markers. Note that the accumulated zone-tracking cost is calculated
based on the actual target set for a fair comparison, as the modified target
zone provided to the controller changes with different $\gamma$. Table 1 acts
as a supplement of Figure 8 and presents the number of violations with respect
to the actual target set after the system converges to the optimal operating
neighborhood for the first time, plus the average magnitude of violation out
of all the violations with different $\gamma$.

Table 1: The number of violations and the average violation value after the first convergence to the actual target set with different $\gamma$. $\gamma$ | # of Vio. | Avg. Vio.
---|---|---
0.3 | 9 | 0.097
0.5 | 4 | 0.045
0.6 | 2 | 0.027
0.7 | 0 | -
1 | 0 | -
3 | 0 | -

The actual target set is bounded by solid black lines in Figure 7. The sets
with $\gamma$ equal to 0.5, 1, and 3 are bounded by green dash-dotted, red
dashed, and blue dotted curves, respectively. The information regarding the
optimal operating condition is not considered, both the upper bound and the
lower bound of the target zone are shrunk equally. Figure 7 indicates the size
of the target set and the corresponding largest terminal set available to the
controller reduce as the value of $\gamma$ increases. This implies the
controller becomes more conservative with a smaller feasible region, which
leads to a sacrifice in the economic performance. The economic performance is
affected in two ways. First, since a smaller feasible zone leads to a shifting
in the feasible economically-optimal operating condition. Second, with a more
conservative target and terminal zone, the controller tends to drive the
system into the zone more aggressively and potentially sacrifice the transit
economic performance. This can be seen from Figure 8 that the accumulated
economic cost function continuously increases as $\gamma$ increases.

It is however noteworthy that a more conservative controller does not
necessarily leads to better zone-tracking performance. It is observed that the
accumulated zone-tracking cost first decreases and then increases after
reaching a minimum at $\gamma=0.6$. From Table 1, it can be observed that at
$\gamma=0.6$, two violations are observed with a very small magnitude of the
average violation. At $\gamma=0.7$, no violation is observed after
convergence, however, the accumulative zone-tracking cost is slightly higher
than that obtained with $\gamma=0.6$. The reason for the increase in zone-
tracking cost as $\gamma$ increase is due to the shrinkage in the feasible
operating range and the more aggressive control action chosen by the
controller, which sacrifices the transit performance of the system. A larger
$\gamma$ helps to reduce the zone-violation and thus the accumulative tracking
cost; however, this reduction is saturated once the violations are
insignificant ($\gamma=0.6$) or are fully eliminated ($\gamma=0.7$). The
sacrifice in the transit performance on the other hand, continuously increases
as $\gamma$ increase, which slowly overtakes the cost reduced by shrinking the
target zone.

###### Remark 6

It is observed from Figure 8 and Table 1 that the optimal zone-tracking
performance is not achieved at $\gamma=1$. This indicates that directly
shrinking the target zone by the $x^{d}_{max}$ only provides a conservative
estimation of the amount of shrinkage required. It is essential to tune
$\gamma$ to obtain a better zone-tracking performance. On the other hand, it
is noteworthy that the economic performance always reduces as $\gamma$
increases. Thus, it is recommended to use the smallest $\gamma$ that provides
reasonable zone-tracking performance in applications.

##  8 Conclusions

In this work, we proposed a robust ZMPC formulation with guaranteed
convergence to the actual target set in the presence of bounded disturbance.
This is achieved by modifying the actual target set, which helps to reject the
effect of the disturbance. A terminal constraint is utilized to ensure the
closed-loop stability of the system. The system state is forced to converge to
a CIS inside the modified target set at the end of the control horizon.
Without assuming the existence of an optimal steady-state operating condition,
this generalized approach provides more degrees of freedom. Furthermore, the
proposed design is able to handle a secondary economic objective without
affecting closed-loop stability. Apart from the theoretical stability proof, a
practical guideline for determining the modified target set is provided. The
proposed design is applied to a CSTR system and is proved to be effective from
various perspectives. The effect of the tuning parameter in the proposed
practical guideline is investigated with the rule of thumb for parameter
selection provided.

##  9 Acknowledgement

Financial support from Natural Sciences and Engineering Research Council of
Canada is gratefully acknowledged.

## References

  * [1] E. S. Meadows, M. A. Henson, J. W. Eaton, and J. B. Rawlings, “Receding horizon control and discontinuous state feedback stabilization,” International Journal of Control, vol. 62, pp. 1217–1229, Nov. 1995.
  * [2] J. B. Rawlings, “Tutorial: model predictive control technology,” in Proceedings of the 1999 American Control Conference (Cat. No. 99CH36251), (San Diego, CA, USA), pp. 662–676 vol.1, IEEE, 1999.
  * [3] J. B. Rawlings, D. Angeli, and C. N. Bates, “Fundamentals of economic model predictive control,” in 2012 IEEE 51st IEEE Conference on Decision and Control (CDC), (Maui, HI, USA), pp. 3851–3861, IEEE, Dec. 2012.
  * [4] M. Diehl, R. Amrit, and J. B. Rawlings, “A Lyapunov Function for Economic Optimizing Model Predictive Control,” IEEE Transactions on Automatic Control, vol. 56, pp. 703–707, Mar. 2011.
  * [5] M. Ellis, H. Durand, and P. D. Christofides, “A tutorial review of economic model predictive control methods,” Journal of Process Control, vol. 24, pp. 1156–1178, Aug. 2014.
  * [6] P. O. M. Scokaert and J. B. Rawlings, “Feasibility issues in linear model predictive control,” AIChE Journal, vol. 45, pp. 1649–1659, Aug. 1999.
  * [7] M. Askari, M. Moghavvemi, H. A. F. Almurib, and A. M. A. Haidar, “Stability of Soft-Constrained Finite Horizon Model Predictive Control,” IEEE Transactions on Industry Applications, vol. 53, Nov. 2017.
  * [8] B. Grosman, E. Dassau, H. C. Zisser, L. Jovanovič, and F. J. Doyle, “Zone Model Predictive Control: A Strategy to Minimize Hyper- and Hypoglycemic Events,” Journal of Diabetes Science and Technology, vol. 4, pp. 961–975, July 2010.
  * [9] A. H. González, P. S. Rivadeneira, A. Ferramosca, N. Magdelaine, and C. H. Moog, “Stable impulsive zone model predictive control for type 1 diabetic patients based on a long‐term model,” Optimal Control Applications and Methods, vol. 41, pp. 2115–2136, Nov. 2020.
  * [10] S. Prívara, J. Široký, L. Ferkl, and J. Cigler, “Model predictive control of a building heating system: The first experience,” Energy and Buildings, vol. 43, pp. 564–572, Feb. 2011.
  * [11] Y. Zhang, B. Decardi-Nelson, J. Liu, J. Shen, and J. Liu, “Zone economic model predictive control of a coal-fired boiler-turbine generating system,” Chemical Engineering Research and Design, vol. 153, pp. 246–256, Jan. 2020.
  * [12] Y. Mao, S. Liu, J. Nahar, J. Liu, and F. Ding, “Soil moisture regulation of agro-hydrological systems using zone model predictive control,” Computers and Electronics in Agriculture, vol. 154, pp. 239–247, Nov. 2018.
  * [13] Z. Huang, J. Liu, and B. Huang, “Model predictive control of agro-hydrological systems based on a two-layer neural network modeling framework,” International Journal of Adaptive Control and Signal Processing, pp. 1–29, 2023.
  * [14] D. Angeli, R. Amrit, and J. B. Rawlings, “On Average Performance and Stability of Economic Model Predictive Control,” IEEE Transactions on Automatic Control, vol. 57, pp. 1615–1626, July 2012.
  * [15] M. Heidarinejad, J. Liu, and P. D. Christofides, “Economic model predictive control of nonlinear process systems using Lyapunov techniques,” AIChE Journal, vol. 58, pp. 855–870, Mar. 2012.
  * [16] S. Liu and J. Liu, “Economic model predictive control with extended horizon,” Automatica, vol. 73, pp. 180–192, Nov. 2016.
  * [17] C. Liu, J. Gao, and D. Xu, “Lyapunov-based model predictive control for tracking of nonholonomic mobile robots under input constraints,” International Journal of Control, Automation and Systems, vol. 15, pp. 2313–2319, Oct. 2017.
  * [18] H. Makhamreh, M. Sleiman, O. Kukrer, and K. Al-Haddad, “Lyapunov-Based Model Predictive Control of a PUC7 Grid-Connected Multilevel Inverter,” IEEE Transactions on Industrial Electronics, vol. 66, pp. 7012–7021, Sept. 2019.
  * [19] S. Liu and J. Liu, “Economic Model Predictive Control with Zone Tracking,” Mathematics, vol. 6, p. 65, Apr. 2018.
  * [20] B. Decardi-Nelson and J. Liu, “Robust Economic Model Predictive Control with Zone Control,” IFAC-PapersOnLine, vol. 54, no. 3, pp. 237–242, 2021.
  * [21] S. Liu, Y. Mao, and J. Liu, “Model-Predictive Control With Generalized Zone Tracking,” IEEE Transactions on Automatic Control, vol. 64, pp. 4698–4704, Nov. 2019.
  * [22] F. Blanchini, “Set invariance in control,” Automatica, vol. 35, pp. 1747 – 1767, 1999.
  * [23] B. Decardi-Nelson and J. Liu, “Computing robust control invariant sets of constrained nonlinear systems: A graph algorithm approach,” Computers & Chemical Engineering, vol. 145, p. 107177, Feb. 2021.
  * [24] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, “CasADi – A software framework for nonlinear optimization and optimal control,” Mathematical Programming Computation, In Press, 2018.