# Sensitivity-based dynamic performance assessment for model predictive
control with Gaussian noise

Jianbang Liua, Song Bob, Benjamin Decardi-Nelsonb, Jinfeng Liub,, Jingtao Huc,
Tao Zoud

aSchool of Electrical and Information Engineering, Jiangsu University,

Zhenjiang 212013, China

bDepartment of Chemical & Materials Engineering, University of Alberta,

Edmonton, AB T6G 1H9, Canada

cShenyang Institute of Automation, Chinese Academy of Sciences,

Shenyang 110016, China

dSchool of Mechanical and Electrical Engineering, Guangzhou University,

Guangzhou 510006, China
Corresponding author: J. Liu. Tel: +1-780-492-1317. Fax: +1-780-492-2881.
Email<EMAIL_ADDRESS>
###### Abstract

Economic model predictive control and tracking model predictive control are
two popular advanced process control strategies used in various of fields.
Nevertheless, which one should be chosen to achieve better performance in the
presence of noise is uncertain when designing a control system. To this end, a
sensitivity-based performance assessment approach is proposed to pre-evaluate
the dynamic economic and tracking performance of them in this work. First,
their controller gains around the optimal steady state are evaluated by
calculating the sensitivities of corresponding constrained dynamic programming
problems. Second, the controller gains are substituted into control loops to
derive the propagation of process and measurement noise. Subsequently, the
Taylor expansion is introduced to simplify the calculation of variance and
mean of each variable. Finally, the tracking and economic performance surfaces
are plotted and the performance indices are precisely calculated through
integrating the objective functions and the probability density functions.
Moreover, boundary moving (i.e., back off) and target moving can be pre-
configured to guarantee the stability of controlled processes using the
proposed approach. Extensive simulations under different cases illustrate the
proposed approach can provide useful guidance on performance assessment and
controller design.

Keywords: Dynamic performance assessment, sensitivity analysis, controller
gain, model predictive control, Gaussian noise.

##  1 Introduction

Over the past decade, economic model predictive control (EMPC), as a
generalized form of model predictive control (MPC), has received significant
attention [3, 25, 12, 9, 10]. It has been shown that the asymptotic average
performance of an EMPC is never worse than the optimal steady-state operation
in the nominal case [3]. In [21], an economic model predictive controller with
extended horizon was proposed and it was shown that in the nominal condition
the transient or asymptotic average performance of EMPC is no worse than the
traditional tracking MPC if the prediction horizon is sufficiently large.
However, in the presence of model-plant mismatch (e.g., measurement and
process noise, process disturbance, parameter uncertainty), whether EMPC can
achieve better economic performance than the traditional tracking MPC is
uncertain as demonstrated in the applications of EMPC to a post-combustion
carbon capture plant [8], a wind energy conversion system [7], and a coal-
fired boiler-turbine system [30]. There is a need of control performance
assessment methods for EMPC in the presence of model-plant mismatch.

In the context of tracking MPC, many studies on performance assessment have
been conducted to determine tracking MPC’s ideal and reachable performance and
to guide the design and implementation of MPC [14, 16, 4, 17, 6]. Majority of
these methods are based on process operation data. For example, in [27], Wang
and co-workers proposed a data-driven subspace approach for control
performance assessment of multi-rate control systems based on minimum variance
control. In [23], the idea was extended and applied to the cross-directional
control of paper machines. In [29], Xu and co-workers developed an economic
performance assessment (EPA) approach and provided constraint/variance tuning
guidelines using routine process operation data plus the process steady-state
gain matrix. In [28], a multivariate economic performance assessment approach
was developed and applied to a MPC controlled electric arc furnace. In [19],
practical and selective tuning guidelines for MPC was proposed based on
sensitivity analysis of process variables. However, the above performance
assessment approaches heavily rely on actual process operation data, which are
not unavailable before the deployment of a control system.

Based on the above considerations, it is clear that there is a lack of methods
that can characterize the control performance of EMPC and tracking MPC in the
presence of model-plant mismatch before their actual deployment. Such a method
can be very useful in determining which predictive control method (tracking
MPC or EMPC) is more appropriate for a specific process. Motivated by this, in
this work, we propose a method to evaluate the performance of EMPC and
tracking MPC in the presence of Gaussian process and measurement noise.

The remainder of this paper is organized as follows: Section 2 gives a brief
introduction of preliminaries, Section 3 illustrates the motivations of this
work, Section 4 introduces the implementation of the proposed sensitivity-
based dynamic performance assessment approach, Section 5 applies the proposed
approach into a chemical process to demonstrate its effectiveness, Section 6
presents the conclusions and future directions.

##  2 Preliminaries

###  2.1 System description and problem formulation

In the work, we consider a general class of discrete time-invariant nonlinear
systems described as follows:

$\displaystyle{\boldsymbol{x}}(k+1)$
$\displaystyle=F({\boldsymbol{x}}(k),{\boldsymbol{u}}(k))+{\boldsymbol{w}}(k)$
(1a) $\displaystyle{\boldsymbol{y}}(k)$
$\displaystyle=H({\boldsymbol{x}}(k))+{\boldsymbol{v}}(k)$ (1b)

where ${\boldsymbol{x}}(k)\in R^{n}$, ${\boldsymbol{u}}(k)\in R^{m}$, and
${\boldsymbol{y}}(k)\in R^{r}$ denote the state, input, and output vectors at
time $k$, respectively; ${\boldsymbol{w}}(k)\in R^{n}$ and
${\boldsymbol{v}}(k)\in R^{r}$ denote the process and measurement noise
vectors at time $k$, respectively; $F(\cdot)$ and $H(\cdot)$ denote the state
and output functions, respectively.

In the work, we focus on characterizing the control performance of economic
MPC and tracking MPC in the presence of Gaussian noise. The results can be
useful when we need to choose between economic MPC and tracking MPC for a
specific control application. The results may also be used to fine tune the
control target (operating point) or constraints for improved control
performance.

###  2.2 Formulation of tracking MPC

Conventional tracking MPC minimizes a quadratic cost to track an optimal
steady state. In the work, we consider tracking MPC that can be formulated as
follows:

$\displaystyle\min\limits_{{\hat{\boldsymbol{u}}(\cdot)},{\hat{\boldsymbol{x}}(\cdot)}}$
$\displaystyle
J(k)=\sum\limits_{i=k}^{k+N-1}{\left\|{{\hat{\boldsymbol{x}}}(i)-{{\boldsymbol{x}}_{\rm{s}}}}\right\|_{\boldsymbol{Q}}^{2}}+\sum\limits_{i=k}^{k+N-1}{\left\|{{\hat{\boldsymbol{u}}}(i)-{{\boldsymbol{u}}_{\rm{s}}}}\right\|_{\boldsymbol{R}}^{2}}$
(2a) s.t.
$\displaystyle{\hat{\boldsymbol{x}}}(i+1)=F({\hat{\boldsymbol{x}}}(i),{\hat{\boldsymbol{u}}}(i)),~{}\forall
i=k,\cdots,k+N-1$ (2b)
$\displaystyle{\hat{\boldsymbol{u}}}(i)\in{\boldsymbol{U}},~{}\forall
i=k,\cdots,k+N-1$ (2c)
$\displaystyle{\hat{\boldsymbol{x}}}(k)={\boldsymbol{x}}(k),~{}{\hat{\boldsymbol{x}}}(i)\in{\boldsymbol{X}},~{}\forall
i=k,\cdots,k+N-1$ (2d)

where $N$ denotes the prediction horizon; ${\hat{\boldsymbol{x}}}(\cdot)$ and
${\hat{\boldsymbol{u}}}(\cdot)$ denote the predicted state and input
trajectories over the prediction horizon, respectively; ${\boldsymbol{x}}(k)$
is the actual system state at instant $k$; ${\boldsymbol{x}}_{\rm s}$ and
${\boldsymbol{u}}_{\rm s}$ are the given steady-state targets to be tracked;
$\boldsymbol{X}$ and $\boldsymbol{U}$ denote the known state and input
constraint sets. In (2), (2a) is the cost function of the tracking MPC, in
which $\boldsymbol{Q}$ and $\boldsymbol{R}$ are positive definite weighting
matrices. (2b) is the system model while (2c) and (2d) being the known
constraints on states and inputs.

###  2.3 Formulation of EMPC

EMPC optimizes a general cost function instead of a quadratic cost function.
In the work, we consider EMPC that can be formulated as follows:

$\displaystyle\min\limits_{{\hat{\boldsymbol{u}}(\cdot)},{\hat{\boldsymbol{x}}(\cdot)}}$
$\displaystyle
J(k)=\sum\limits_{i=k}^{k+N-1}{E({\hat{\boldsymbol{x}}}(i),{\hat{\boldsymbol{u}}}(i))}$
(3a) s.t.
$\displaystyle{\hat{\boldsymbol{x}}}(i+1)=F({\hat{\boldsymbol{x}}}(i),{\hat{\boldsymbol{u}}}(i)),~{}\forall
i=k,\cdots,k+N-1$ (3b)
$\displaystyle{\hat{\boldsymbol{u}}}(i)\in{\boldsymbol{U}},~{}\forall
i=k,\cdots,k+N-1$ (3c)
$\displaystyle{\hat{\boldsymbol{x}}}(k)={\boldsymbol{x}}(k),~{}{\hat{\boldsymbol{x}}}(i)\in{\boldsymbol{X}},~{}\forall
i=k,\cdots,k+N-1$ (3d)

where (3a) is the cost function of the EMPC to minimize and $E(\cdot)$ denotes
the economic stage cost.

###  2.4 Performance indices

In this work, two different performance indices are used: the average tracking
performance index $\bar{J}_{\rm tp}$ and the average economic performance
index $\bar{J}_{\rm ep}$ [4, 9]:

$\displaystyle{\bar{J}_{{\rm{tp}}}}$
$\displaystyle=\frac{1}{{{k_{f}}-{k_{0}}+1}}\sum\limits_{k={k_{0}}}^{{k_{f}}}{{J_{{\rm{tp}}}}(k)}$
(4a) $\displaystyle{\bar{J}_{{\rm{ep}}}}$
$\displaystyle=\frac{1}{{{k_{f}}-{k_{0}}+1}}\sum\limits_{k={k_{0}}}^{{k_{f}}}{{J_{{\rm{ep}}}}(k)}$
(4b)

where ${J_{{\rm{tp}}}}(k)$ and ${J_{{\rm{ep}}}}(k)$ defined as follows
denoting the instantaneous tracking and economic indices at time $k$:

$\begin{array}[]{l}{J_{{\rm{tp}}}}(k)={\left\|{{\boldsymbol{x}}(k)-{{\boldsymbol{x}}_{\rm{s}}}}\right\|_{\boldsymbol{Q}}^{2}+\left\|{{\boldsymbol{u}}(k)-{{\boldsymbol{u}}_{\rm{s}}}}\right\|_{\boldsymbol{R}}^{2}}\\\
{J_{{\rm{ep}}}}(k)={E({\boldsymbol{x}}(k),{\boldsymbol{u}}(k))}\end{array}$

##  3 Motivating example

Figure 1: Schematic diagram of the CSTR

We consider a nonlinear non-isothermal continuous stirred tank reactor (CSTR)
as shown in Figure 1 to illustrate the motivations of this work. The dynamics
of the CSTR are described as follows [13]:

$\displaystyle\frac{{d{C_{A}}}}{{dt}}=\frac{F}{{{V_{R}}}}({C_{A0}}-{C_{A}})-{k_{0}}{e^{-\frac{E}{{RT}}}}C_{A}^{2}$
(5a)
$\displaystyle\frac{{dT}}{{dt}}=\frac{F}{{{V_{R}}}}({T_{0}}-T)-\frac{{\Delta
H{k_{0}}}}{{{\rho_{R}}{c_{p}}}}{e^{-\frac{E}{{RT}}}}C_{A}^{2}+\frac{Q}{{{\rho_{R}}{c_{p}}{V_{R}}}}$
(5b)

where ${C_{A}}$ is the molar concentration of the reactant $A$; $T$ is the
temperature of reactor contents; $F$ is the outlet flow rate of the reactor;
$Q$ is the heat rate. The remaining notation definitions and process parameter
values are given in Table 1. The continuous model is discretized using forward
finite difference with a sampling time $\Delta T=10$ sec. At each sampling
time, ${C_{A}}$ and $T$ are both measured. For this process, the input vector
${\boldsymbol{u}}=[F~{}Q]^{\rm T}$, the state vector
${\boldsymbol{x}}=[C_{A}~{}T]^{\rm T}$, and the output vector
${\boldsymbol{y}}=[C_{A}~{}T]^{\rm T}$.

Table 1: Process parameters of the CSTR Symbol | Description | Value
---|---|---
$C_{A0}$ | Feed concentration of A | 3.5 ${\rm{kmol/m^{3}}}$
$T_{0}$ | Feedstock temperature | 300 $\rm K$
$V_{R}$ | Reactor fluid volume | 1.0 $\rm{m^{3}}$
$E$ | Activation energy | 5.0e4 $\rm{kJ/kmol}$
$k_{0}$ | Pre-exponential rate factor | 8.46e6 $\rm{m^{3}/kmol\cdot h}$
$\Delta H$ | Reaction enthalpy change | -1.16e4 $\rm{kJ/kmol}$
$c_{p}$ | Heat capacity | 0.231 $\rm{kJ/kg\cdot K}$
$\rho_{R}$ | Density | 1000 $\rm{kg/m^{3}}$
$R$ | Gas constant | 8.314 $\rm{kJ/kmol\cdot K}$

In this section, we configure two group experiments to show the different
control performance of EMPC and tracking MPC in the presence of noise. In each
group, EMPC and tracking MPC finally arrive at the same optimal steady state.
In group one, the tracking and economic indices expressed in Eqs. (2), (3),
and (4) are set as:

$\displaystyle{J_{\rm{tp}}}$
$\displaystyle={\left\|{{\boldsymbol{x}}(k)-{{\boldsymbol{x}}_{\rm{s}}}}\right\|_{\boldsymbol{Q}}^{2}+\left\|{{\boldsymbol{u}}(k)-{{\boldsymbol{u}}_{\rm{s}}}}\right\|_{\boldsymbol{R}}^{2}}$
(6a) $\displaystyle{J_{\rm{ep}}}$
$\displaystyle={{x_{1}}{u_{1}}+10^{5}{x_{2}}+\frac{2}{{{{({x_{1}}-1)}^{2}}}}+\frac{10^{4}}{{{{({x_{2}}-600)}^{2}}}}}$
(6b)

where the last two terms in Eq. (6b) are known as barrier functions to push
the corresponding variables away from their boundaries, $\boldsymbol{Q}$ and
$\boldsymbol{R}$ are configured as follows to avoid the influence of parameter
tuning in tracking MPC.

$\begin{array}[]{l}{\boldsymbol{Q}}={\rm{diag}}\left\\{{1/x_{1,{\rm{s}}}^{2},1/x_{2,{\rm{s}}}^{2},\cdots,1/x_{n,{\rm{s}}}^{2}}\right\\}\\\
{\boldsymbol{R}}={\rm{diag}}\left\\{{1/u_{1,{\rm{s}}}^{2},1/u_{2,{\rm{s}}}^{2},\cdots,1/u_{m,{\rm{s}}}^{2}}\right\\}\end{array}$

The constraints considered in the experiments are shown as follows:

$1\leq{x_{1}}\leq 3,~{}600\leq{x_{2}}\leq 1000$ (7)

The process and measurement noise added on states and measurements is as
follows:

$\begin{array}[]{l}{w_{1}}\sim{\rm N}(0,0.00001),{v_{1}}\sim{\rm N}(0,0.01)\\\
{w_{2}}\sim{\rm N}(0,0.00009),{v_{2}}\sim{\rm N}(0,0.09)\end{array}$ (8)

where ${\rm N}(\mu,\sigma^{2})$ denotes the normal distribution function with
mean $\mu$ and variance $\sigma^{2}$.

Group two adopts the same constraint and noise configuration except for a
different economic index shown as below:

$\begin{array}[]{l}{J_{\rm{ep}}}=10^{2}x_{1}+10^{-4}x_{2}u_{1}-2\times
10^{2}\ln(\frac{x_{2}-600}{200})-10^{-3}\ln(\frac{1000-x_{2}}{200})\\\ \kern
125.0pt-10^{-3}\ln(3-x_{1})-1.5\times 10^{3}\ln(x_{1}-1)\end{array}$ (9)

Table 2 shows the average tracking and economic performance of tracking MPC
and EMPC. It is shown that sometimes EMPC achieves better economic performance
than tracking MPC, sometimes not. Therefore, an unified conclusion about which
one is better cannot be directly given.

Table 2: Average performance indices of EMPC and tracking MPC in the presence of noise Group | Controller | Steady state ${\boldsymbol{x}}_{\rm s}$ | $\bar{J}_{\rm{ep}}$ | $\bar{J}_{\rm{tp}}$
---|---|---|---|---
1 | Economic MPC | $[1.1601~{}615.7373]^{\rm T}$ | 453.5738 | 2.0905e-03
Tracking MPC | $[1.1601~{}615.7373]^{\rm T}$ | 448.9241 | 1.0275e-03
2 | Economic MPC | $[2.7284~{}631.7380]^{\rm T}$ | 199.8444 | 7.1578e-03
Tracking MPC | $[2.7284~{}631.7380]^{\rm T}$ | 202.5862 | 6.7702e-03

Therefore, in this work we aim at providing an approach to analyze the
actually accessible performance and giving guidance on which controller should
be employed for a given chemical process. We start from the Gaussian noise and
propose a sensitivity based dynamic performance assessment approach to pre-
evaluate the actual performance of tracking MPC and EMPC. The results are also
helpful to conduct boundary moving (i.e. back off) and target moving when
$\boldsymbol{x}$ is located on or close to its boundaries, to guarantee the
safety of controlled systems.

##  4 Sensitivity-based dynamic performance assessment

Sensitivity analysis is to evaluate the sensitivity of the optimal solution
with respect to given parameters in constrained optimization problems [11, 26,
15, 2]. For a MPC problem, we usually convert it to a standard nonlinear
programming formulation, then analyze its sensitivity information[26, 2]. This
technology has been applied in industrial process control to reduce the
computational burden as well as the computational time [24, 22, 5]. In this
work, we use it to analyze the sensitivity of the input $\boldsymbol{u}$ with
respect to the state ${\boldsymbol{x}}$ and use it as the gain of the MPC
controller. Using the calculated controller gain, we can derive the
propagation of noise in the closed-loop system and get a noise transfer model.
Based on the model, we can calculate the mean and variance of each state which
reflect the distribution information of state variables. Finally, we can
evaluate the economic and tracking performance of the controlled system by
integrating the closed-loop performance functions which include only states
and the calculated distribution functions of states. Figure 2 shows the
implementation procedure of the proposed performance assessment approach. The
algorithm proposed in [2] is employed to carry out sensitivity analysis using
CasADi with Python in the work [1].

Figure 2: Implementation procedure of the proposed performance assessment
approach

###  4.1 Sensitivity based controller gain

First, we calculate the controller gain of a MPC controller using sensitivity
analysis in the section. We treat the initial state ${\boldsymbol{x}}(0)$ as a
sensitivity parameter. Future state ${\boldsymbol{x}}(1:N)$ and input
${\boldsymbol{u}}(0:N-1)$ are decision variables.

Let us reform the MPC problems (2) and (3) to the following standard nonlinear
programming (NLP) formulation:

$\begin{array}[]{l}\mathop{\min}\limits_{{\boldsymbol{d}},{\boldsymbol{p}}}f({\boldsymbol{d}},{\boldsymbol{p}})\\\
{\rm{s}}{\rm{.t}}{\rm{.}}\\\
{\boldsymbol{p}}={{\boldsymbol{p}}_{0}},{{\boldsymbol{d}}^{{\rm{lb}}}}\leq{\boldsymbol{d}}\leq{{\boldsymbol{d}}^{{\rm{ub}}}},{{\boldsymbol{g}}^{{\rm{lb}}}}\leq
g({\boldsymbol{d}},{\boldsymbol{p}})\leq{{\boldsymbol{g}}^{{\rm{ub}}}}\end{array}$
(10)

where

$\begin{array}[]{l}{\boldsymbol{d}}:=\left[{{\boldsymbol{x}}(1)^{\rm
T}~{}{\boldsymbol{x}}(2)^{\rm T}~{}\cdots~{}{\boldsymbol{x}}(N)^{\rm
T}~{}{\boldsymbol{u}}(0)^{\rm T}~{}{\boldsymbol{u}}(1)^{\rm
T}~{}\cdots~{}{\boldsymbol{u}}(N-1)^{\rm T}}\right]^{\rm T}\\\
{\boldsymbol{p}}:=\left[{{\boldsymbol{x}}(0)^{\rm
T}~{}{\boldsymbol{\tilde{p}}}^{\rm T}}\right]^{\rm T}\end{array}$

and $N$ is the prediction horizon; ${\boldsymbol{x}}(k)$ and
${\boldsymbol{u}}(k)$ denote the state and input vectors at time $k$,
respectively; ${\boldsymbol{d}}$ denotes the decision variable, which is
formed by the state vector $\boldsymbol{x}(\cdot)$ and the input vector
$\boldsymbol{u}(\cdot)$; ${\boldsymbol{p}}$ includes the parameters which need
to conduct sensitivity analysis and in our work it only contains the initial
state ${\boldsymbol{x}}{(0)}$, note that if necessary other sensitivity
parameter ${\boldsymbol{\tilde{p}}}$ can also be included;
${\boldsymbol{p}}_{0}$ denotes the nominal value of $\boldsymbol{p}$;
$f(\cdot)$ and $g(\cdot)$ denote the objective and constraint functions,
respectively; ${\boldsymbol{d}}^{\rm{lb}}$, ${\boldsymbol{d}}^{\rm{ub}}$,
${\boldsymbol{g}}^{\rm{lb}}$, and ${\boldsymbol{g}}^{\rm{ub}}$ denote the
lower and upper boundaries of $\boldsymbol{d}$ and $g(\cdot)$, respectively.

To solve the above constrained optimization problem, we need to introduce
Lagrangian multipliers ${\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{lb}}}$,
${\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{ub}}}$,
${\boldsymbol{\lambda}}_{g}^{{\rm{lb}}}$,
${\boldsymbol{\lambda}}_{g}^{{\rm{ub}}}$, and
${{\boldsymbol{\lambda}}_{\boldsymbol{p}}}$ to convert it to an unconstrained
optimization problem as follows:

$\begin{array}[]{l}\tilde{L}({\boldsymbol{d}},{\boldsymbol{p}},{\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{lb}}},{\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{ub}}},{\boldsymbol{\lambda}}_{g}^{{\rm{lb}}},{\boldsymbol{\lambda}}_{g}^{{\rm{ub}}},{{\boldsymbol{\lambda}}_{\boldsymbol{p}}})={\mathop{f}\nolimits}({\boldsymbol{d}},{\boldsymbol{p}})+{{\boldsymbol{\lambda}}_{\boldsymbol{p}}}({\boldsymbol{p}}-{{\boldsymbol{p}}_{0}})-{\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{lb}}}({\boldsymbol{d}}-{{\boldsymbol{d}}^{{\rm{lb}}}})-{\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{ub}}}({{\boldsymbol{d}}^{{\rm{ub}}}}-{\boldsymbol{d}})\\\
\kern
194.0pt-{\boldsymbol{\lambda}}_{g}^{{\rm{lb}}}(g({\boldsymbol{d}},{\boldsymbol{p}})-{{\boldsymbol{g}}^{{\rm{lb}}}})-{\boldsymbol{\lambda}}_{g}^{{\rm{ub}}}({{\boldsymbol{g}}^{{\rm{ub}}}}-g({\boldsymbol{d}},{\boldsymbol{p}}))\end{array}$
(11)

The Karush-Kuhn-Tucker(KKT) conditions for the above unconstrained problem
(11) is as follows:

$\begin{array}[]{l}{\nabla_{\boldsymbol{d}}}\tilde{L}({\boldsymbol{d}},{\boldsymbol{p}},{\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{lb}}},{\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{ub}}},{\boldsymbol{\lambda}}_{g}^{{\rm{lb}}},{\boldsymbol{\lambda}}_{g}^{{\rm{ub}}},{{\boldsymbol{\lambda}}_{\boldsymbol{p}}})={\boldsymbol{0}}\\\
{\nabla_{\boldsymbol{p}}}\tilde{L}({\boldsymbol{d}},{\boldsymbol{p}},{\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{lb}}},{\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{ub}}},{\boldsymbol{\lambda}}_{g}^{{\rm{lb}}},{\boldsymbol{\lambda}}_{g}^{{\rm{ub}}},{{\boldsymbol{\lambda}}_{\boldsymbol{p}}})={\boldsymbol{0}}\\\
{\boldsymbol{p}}={{\boldsymbol{p}}_{0}},{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{\kern 1.0pt}{\kern
1.0pt}{\boldsymbol{d}}-{{\boldsymbol{d}}^{{\rm{lb}}}}\geq{\boldsymbol{0}},{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{{\boldsymbol{d}}^{{\rm{ub}}}}-{\boldsymbol{d}}\geq{\boldsymbol{0}}\\\
g({\boldsymbol{d}},{\boldsymbol{p}})-{{\boldsymbol{g}}^{{\rm{lb}}}}\geq{\boldsymbol{0}},{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{{\boldsymbol{g}}^{{\rm{ub}}}}-g({\boldsymbol{d}},{\boldsymbol{p}})\geq{\boldsymbol{0}}\\\
{\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{lb}}}\geq{\boldsymbol{0}},{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{ub}}}\geq{\boldsymbol{0}},{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{\boldsymbol{\lambda}}_{g}^{{\rm{lb}}}\geq{\boldsymbol{0}},{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{\boldsymbol{\lambda}}_{g}^{{\rm{ub}}}\geq{\boldsymbol{0}}\\\
{\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{lb}}}({\boldsymbol{d}}-{{\boldsymbol{d}}^{{\rm{lb}}}})={\boldsymbol{0}},{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{ub}}}({{\boldsymbol{d}}^{{\rm{ub}}}}-{\boldsymbol{d}})={\boldsymbol{0}}\\\
{\boldsymbol{\lambda}}_{g}^{{\rm{lb}}}(g({\boldsymbol{d}},{\boldsymbol{p}})-{{\boldsymbol{g}}^{{\rm{lb}}}})={\boldsymbol{0}},{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{\boldsymbol{\lambda}}_{g}^{{\rm{ub}}}({{\boldsymbol{g}}^{{\rm{ub}}}}-g({\boldsymbol{d}},{\boldsymbol{p}}))={\boldsymbol{0}}\end{array}$
(12)

To simplify the above KKT conditions (12), we introduce
${{\boldsymbol{\lambda}}_{\boldsymbol{d}}}={\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{ub}}}-{\boldsymbol{\lambda}}_{\boldsymbol{d}}^{{\rm{lb}}}$
and
${{\boldsymbol{\lambda}}_{g}}={\boldsymbol{\lambda}}_{g}^{{\rm{ub}}}-{\boldsymbol{\lambda}}_{g}^{{\rm{lb}}}$
and a new Lagrangian function with all linear terms eliminated can be obtained
as follows:

$L({\boldsymbol{d}},{\boldsymbol{p}},{{\boldsymbol{\lambda}}_{g}})=f({\boldsymbol{d}},{\boldsymbol{p}})+{{\boldsymbol{\lambda}}_{g}}g({\boldsymbol{d}},{\boldsymbol{p}})$
(13)

The new equivalent KKT conditions can be formulated as below:

$\displaystyle{\nabla_{\boldsymbol{d}}}L({\boldsymbol{d}},{\boldsymbol{p}},{{\boldsymbol{\lambda}}_{g}})+{{\boldsymbol{\lambda}}_{\boldsymbol{d}}}={\boldsymbol{0}}$
(14a)
$\displaystyle{\nabla_{\boldsymbol{p}}}L({\boldsymbol{d}},{\boldsymbol{p}},{{\boldsymbol{\lambda}}_{g}})+{{\boldsymbol{\lambda}}_{\boldsymbol{p}}}={\boldsymbol{0}}$
(14b)
$\displaystyle{{\boldsymbol{d}}^{{\rm{lb}}}}\leq{\boldsymbol{d}}\leq{{\boldsymbol{d}}^{{\rm{ub}}}},{{\boldsymbol{g}}^{{\rm{lb}}}}\leq
g({\boldsymbol{d}},{\boldsymbol{p}})\leq{{\boldsymbol{g}}^{{\rm{ub}}}}$ (14c)
$\displaystyle{\boldsymbol{I}}_{\boldsymbol{d}}^{0}{{\boldsymbol{\lambda}}_{\boldsymbol{d}}}+{\boldsymbol{I}}_{\boldsymbol{d}}^{{\rm{ub}}}{{\boldsymbol{d}}^{{\rm{ub}}}}+{\boldsymbol{I}}_{\boldsymbol{d}}^{{\rm{lb}}}{{\boldsymbol{d}}^{{\rm{lb}}}}-{{\boldsymbol{I}}_{\boldsymbol{d}}}{\boldsymbol{d}}={\boldsymbol{0}}$
(14d)
$\displaystyle{\boldsymbol{I}}_{g}^{0}{{\boldsymbol{\lambda}}_{g}}+{\boldsymbol{I}}_{g}^{{\rm{ub}}}{{\boldsymbol{g}}^{{\rm{ub}}}}+{\boldsymbol{I}}_{g}^{{\rm{lb}}}{{\boldsymbol{g}}^{{\rm{lb}}}}-{{\boldsymbol{I}}_{g}}g({\boldsymbol{d}},{\boldsymbol{p}})={\boldsymbol{0}}$
(14e)

where ${\boldsymbol{I}}_{\bullet}^{{\rm{lb}}}$,
${\boldsymbol{I}}_{\bullet}^{{\rm{ub}}}$, ${{\boldsymbol{I}}_{\bullet}}$, and
${\boldsymbol{I}}_{\bullet}^{0}$ are the active constraints flags that denote
whether the corresponding constraints are active or not [2].

Using (14a), we can eliminate ${{\boldsymbol{\lambda}}_{\boldsymbol{d}}}$ from
(14d) and get an implicit function $F({\boldsymbol{z}},{\boldsymbol{q}})$ as
follows:

$F({\boldsymbol{z}},{\boldsymbol{q}}):=\left[{\begin{array}[]{*{20}{c}}{{\boldsymbol{I}}_{\boldsymbol{d}}^{0}{\nabla_{\boldsymbol{d}}}L({\boldsymbol{d}},{\boldsymbol{p}},{{\boldsymbol{\lambda}}_{g}})+{{\boldsymbol{I}}_{\boldsymbol{d}}}{\boldsymbol{d}}-{\boldsymbol{I}}_{\boldsymbol{d}}^{{\rm{ub}}}{{\boldsymbol{d}}^{{\rm{ub}}}}-{\boldsymbol{I}}_{\boldsymbol{d}}^{{\rm{lb}}}{{\boldsymbol{d}}^{{\rm{lb}}}}}\\\
{{{\boldsymbol{I}}_{g}}g({\boldsymbol{d}},{\boldsymbol{p}})-{\boldsymbol{I}}_{g}^{0}{{\boldsymbol{\lambda}}_{g}}-{\boldsymbol{I}}_{g}^{{\rm{ub}}}{{\boldsymbol{g}}^{{\rm{ub}}}}-{\boldsymbol{I}}_{g}^{{\rm{lb}}}{{\boldsymbol{g}}^{{\rm{lb}}}}}\end{array}}\right]={\boldsymbol{0}}$
(15)

where

$\begin{array}[]{l}{\boldsymbol{z}}:=[{\boldsymbol{d}}^{\rm
T}~{}{{\boldsymbol{\lambda}}_{g}}^{\rm T}]^{\rm T}\\\
{\boldsymbol{q}}:=[{{\boldsymbol{d}}^{{\rm{lb}}}}^{\rm
T}~{}{{\boldsymbol{d}}^{{\rm{ub}}}}^{\rm
T}~{}{{\boldsymbol{g}}^{{\rm{lb}}}}^{\rm
T}~{}{{\boldsymbol{g}}^{{\rm{ub}}}}^{\rm T}~{}{\boldsymbol{p}}^{\rm T}]^{\rm
T}\end{array}$

If we define a corresponding explicit function
${\boldsymbol{z}}:=G({\boldsymbol{q}})$, we can get the sensitivity of
$\boldsymbol{z}$ with respect to $\boldsymbol{q}$,
$\frac{d\boldsymbol{z}}{d\boldsymbol{q}}$, using the implicit function theorem
[18]:

$\frac{{d\boldsymbol{z}}}{{d{\boldsymbol{q}}}}=-{\left[{\frac{{\partial
F}}{{\partial{\boldsymbol{z}}}}}\right]^{-1}}\frac{{\partial
F}}{{\partial{\boldsymbol{q}}}}$ (16)

where

$\frac{{\partial
F}}{{\partial{\boldsymbol{z}}}}=\left[{\begin{array}[]{*{20}{c}}{{\boldsymbol{I}}_{\boldsymbol{d}}^{0}\nabla_{\boldsymbol{d}}^{2}L({\boldsymbol{d}},{\boldsymbol{p}},{{\boldsymbol{\lambda}}_{g}})+{{\boldsymbol{I}}_{\boldsymbol{d}}}}&{{\boldsymbol{I}}_{\boldsymbol{d}}^{0}{{\left[{\frac{{\partial
g}}{{\partial{\boldsymbol{d}}}}({\boldsymbol{d}},{\boldsymbol{p}})}\right]}^{\rm{T}}}}\\\
{{{\boldsymbol{I}}_{g}}\frac{{\partial
g}}{{\partial{\boldsymbol{d}}}}({\boldsymbol{d}},{\boldsymbol{p}})}&{-{\boldsymbol{I}}_{g}^{0}}\end{array}}\right]$
$\frac{{\partial
F}}{{\partial{\boldsymbol{q}}}}=\left[{\begin{array}[]{*{20}{c}}{-{\boldsymbol{I}}_{\boldsymbol{d}}^{{\rm{lb}}}}&{-{\boldsymbol{I}}_{\boldsymbol{d}}^{{\rm{ub}}}}&{\boldsymbol{0}}&{\boldsymbol{0}}&{{\boldsymbol{I}}_{\boldsymbol{d}}^{0}\frac{\partial}{{\partial{\boldsymbol{p}}}}{\nabla_{\boldsymbol{d}}}L({\boldsymbol{d}},{\boldsymbol{p}},{{\boldsymbol{\lambda}}_{g}})}\\\
{\boldsymbol{0}}&{\boldsymbol{0}}&{-{\boldsymbol{I}}_{g}^{{\rm{lb}}}}&{-{\boldsymbol{I}}_{g}^{{\rm{ub}}}}&{{{\boldsymbol{I}}_{g}}\frac{\partial}{{\partial{\boldsymbol{p}}}}g({\boldsymbol{d}},{\boldsymbol{p}})}\end{array}}\right]$

In MPC, the controller gain ${{\boldsymbol{K}}^{{\rm{MPC}}}}$ can be regarded
as the sensitivity of ${\boldsymbol{u}}(0)$ with respect to
${\boldsymbol{x}}(0)$, which is a sub-block of the above calculated
sensitivity matrix $\frac{d\boldsymbol{z}}{d\boldsymbol{q}}$. Therefore, we
can get the following controller gain ${{\boldsymbol{K}}^{{\rm{MPC}}}}$ from
the calculated $\frac{d\boldsymbol{z}}{d\boldsymbol{q}}$.

${{\boldsymbol{K}}^{{\rm{MPC}}}}={\left.\frac{{\partial{\boldsymbol{u}}(0)}}{{\partial{\boldsymbol{x}}(0)}}\right|_{{{\boldsymbol{x}}_{0}}}}$
(17)

where ${{\boldsymbol{K}}^{{\rm{MPC}}}}$ is the controller gain of MPC and
${\boldsymbol{x}}_{0}$ denotes the nominal value of ${\boldsymbol{x}}(0)$.

Finally, the optimal feedback control law of MPC can be formulated as below:

${\boldsymbol{u}}(0)={{\boldsymbol{u}}_{\rm{0}}}+{{\boldsymbol{K}}^{\rm{MPC}}}({\boldsymbol{x}}(0)-{{\boldsymbol{x}}_{\rm{0}}})$
(18)

where ${\boldsymbol{u}}_{0}$ denotes the nominal optimal solution,
${\boldsymbol{x}}(0)$ is the actual state value at time 0, and
${\boldsymbol{u}}(0)$ is the actual optimal input at time 0.

###  4.2 The propagation of process and measurement noise

Based on the above calculated controller gain, we can derive the noise
propagation of noise in a closed-loop controlled system. Assume that Gaussian
process and measurement noise is added on states and measurements and all
states are directly measured, the complete closed-loop system around the given
steady state can be described as follows:

$\displaystyle{\boldsymbol{y}}(k)={{\boldsymbol{x}}}(k)+{\boldsymbol{v}}(k)$
(19a) $\displaystyle{{\boldsymbol{x}}_{\rm{m}}}(k)={\boldsymbol{y}}(k)$ (19b)
$\displaystyle{\boldsymbol{u}}(k)={{\boldsymbol{u}}_{\rm{s}}}+{{\boldsymbol{K}}^{{\rm{MPC}}}}({{\boldsymbol{x}}_{\rm{m}}}(k)-{{\boldsymbol{x}}_{\rm{s}}})$
(19c)
$\displaystyle{{\boldsymbol{x}}}(k+1)=F({{\boldsymbol{x}}}(k),{\boldsymbol{u}}(k))+{\boldsymbol{w}}(k)$
(19d)

where ${\boldsymbol{x}}(k)$ and ${\boldsymbol{x}}_{\rm m}(k)$ denote the
actual and measured state values at time $k$, respectively;
${\boldsymbol{x}}(k)$ is usually supposed to be unknown and
${\boldsymbol{x}}_{\rm m}(k)$ instead of ${\boldsymbol{x}}(k)$ is used to
compute the optimal input ${\boldsymbol{u}}(k)$.

Finally, Eq. (20) can be achieved by successively substituting (19c), (19b),
and (19a) into (19d):

$\begin{array}[]{l}{{\boldsymbol{x}}}(k+1)=F({{\boldsymbol{x}}}(k),{{\boldsymbol{u}}_{\rm{s}}}+{{\boldsymbol{K}}^{{\rm{MPC}}}}({{\boldsymbol{x}}}(k)+{\boldsymbol{v}}(k)-{{\boldsymbol{x}}_{\rm{s}}}))+{\boldsymbol{w}}(k)\end{array}$
(20)

It is clear that Eq. (20) only contains the actual state
${\boldsymbol{x}}(\cdot)$, process noise ${\boldsymbol{w}}(k)$, and
measurement noise ${\boldsymbol{v}}(k)$, which reveals the propagation of
process and measurement noise in a closed-loop system.

###  4.3 Calculation of variance and mean

Based on the above noise transfer function, we can calculate the mean and
variance of each variable, which can be used to get the approximate
distribution of states in a well-controlled system. However, although the
propagation of process and measurement noise in a closed-loop controlled
system has been revealed as Eq. (20), calculating the variances and means of
relevant variables is still rather challenging. To simplify the calculation,
Taylor expansion is introduced to approximate the nonlinear model (20) around
the given steady state $({\boldsymbol{u}}_{\rm s},{\boldsymbol{x}}_{\rm s})$.

Let us denote Eq. (20) by follows:

${{\boldsymbol{x}}}(k+1):=\varphi({{\boldsymbol{x}}}(k),{\boldsymbol{w}}(k),{\boldsymbol{v}}(k))$
(21)

The first-order Taylor expansion of (20) will be as follows:

${{\boldsymbol{x}}}(k+1)=\varphi({{\boldsymbol{x}}_{\rm{s}}},{{\boldsymbol{w}}_{\rm{s}}},{{\boldsymbol{v}}_{\rm{s}}})+{\left.\frac{{\partial\varphi}}{{\partial{{\boldsymbol{x}}}}}\right|_{{{\boldsymbol{x}}_{\rm{s}}},{{\boldsymbol{w}}_{\rm{s}}},{{\boldsymbol{v}}_{\rm{s}}}}}({{\boldsymbol{x}}}(k)-{{\boldsymbol{x}}_{\rm{s}}})+{\left.\frac{{\partial\varphi}}{{\partial{{\boldsymbol{w}}}}}\right|_{{{\boldsymbol{x}}_{\rm{s}}},{{\boldsymbol{w}}_{\rm{s}}},{{\boldsymbol{v}}_{\rm{s}}}}}({\boldsymbol{w}}(k)-{{\boldsymbol{w}}_{\rm{s}}})+{\left.\frac{{\partial\varphi}}{{\partial{{\boldsymbol{v}}}}}\right|_{{{\boldsymbol{x}}_{\rm{s}}},{{\boldsymbol{w}}_{\rm{s}}},{{\boldsymbol{v}}_{\rm{s}}}}}({\boldsymbol{v}}(k)-{{\boldsymbol{v}}_{\rm{s}}})$
(22)

If we define

$\begin{array}[]{l}{{\boldsymbol{K}}_{{\boldsymbol{x}}2{\boldsymbol{x}}}}:={\left.\frac{{\partial\varphi}}{{\partial{{\boldsymbol{x}}}}}\right|_{{{\boldsymbol{x}}_{\rm{s}}},{{\boldsymbol{w}}_{\rm{s}}},{{\boldsymbol{v}}_{\rm{s}}}}}\\\
{{\boldsymbol{K}}_{{\boldsymbol{x}}2{\boldsymbol{w}}}}:={\left.\frac{{\partial\varphi}}{{\partial{{\boldsymbol{w}}}}}\right|_{{{\boldsymbol{x}}_{\rm{s}}},{{\boldsymbol{w}}_{\rm{s}}},{{\boldsymbol{v}}_{\rm{s}}}}}\\\
{{\boldsymbol{K}}_{{\boldsymbol{x}}2{\boldsymbol{v}}}}:={\left.\frac{{\partial\varphi}}{{\partial{{\boldsymbol{v}}}}}\right|_{{{\boldsymbol{x}}_{\rm{s}}},{{\boldsymbol{w}}_{\rm{s}}},{{\boldsymbol{v}}_{\rm{s}}}}}\end{array}$
(23)

Eq. (22) can be expressed as follows:

${{\boldsymbol{x}}}(k+1)=\varphi({{\boldsymbol{x}}_{\rm{s}}},{{\boldsymbol{w}}_{\rm{s}}},{{\boldsymbol{v}}_{\rm{s}}})+{{\boldsymbol{K}}_{{\boldsymbol{x}}2{\boldsymbol{x}}}}({{\boldsymbol{x}}}(k)-{{\boldsymbol{x}}_{\rm{s}}})+{{\boldsymbol{K}}_{{\boldsymbol{x}}2{\boldsymbol{w}}}}({\boldsymbol{w}}(k)-{{\boldsymbol{w}}_{\rm{s}}})+{{\boldsymbol{K}}_{{\boldsymbol{x}}2{\boldsymbol{v}}}}({\boldsymbol{v}}(k)-{{\boldsymbol{v}}_{\rm{s}}})$
(24)

Based on Eq. (24), we can achieve the following covariance matrix equation:

${\boldsymbol{\Sigma}}_{{{\boldsymbol{x}}}}={{\boldsymbol{K}}_{{\boldsymbol{x}}2{\boldsymbol{x}}}}{\boldsymbol{\Sigma}}_{{{\boldsymbol{x}}}}{\boldsymbol{K}}_{{\boldsymbol{x}}2{\boldsymbol{x}}}^{\rm{T}}+{{\boldsymbol{K}}_{{\boldsymbol{x}}2{\boldsymbol{w}}}}{\boldsymbol{\Sigma}}_{\boldsymbol{w}}{\boldsymbol{K}}_{{\boldsymbol{x}}2{\boldsymbol{w}}}^{\rm{T}}+{{\boldsymbol{K}}_{{\boldsymbol{x}}2{\boldsymbol{v}}}}{\boldsymbol{\Sigma}}_{\boldsymbol{v}}{\boldsymbol{K}}_{{\boldsymbol{x}}2{\boldsymbol{v}}}^{\rm{T}}$
(25)

where ${\boldsymbol{\Sigma}}_{{{\boldsymbol{x}}}}$,
${\boldsymbol{\Sigma}}_{\boldsymbol{w}}$, and
${\boldsymbol{\Sigma}}_{\boldsymbol{v}}$ denote the covariance matrices of
${\boldsymbol{x}}$, $\boldsymbol{w}$, and $\boldsymbol{v}$, respectively;
${\boldsymbol{\Sigma}}_{\boldsymbol{w}}$ and
${\boldsymbol{\Sigma}}_{\boldsymbol{v}}$ are supposed to be known.

We can also achieve the following mean vector equation:

${{\boldsymbol{\mu}}_{{{\boldsymbol{x}}}}}=\varphi({{\boldsymbol{x}}_{\rm{s}}},{{\boldsymbol{w}}_{\rm{s}}},{{\boldsymbol{v}}_{\rm{s}}})={{\boldsymbol{x}}_{\rm{s}}}$
(26)

where ${\boldsymbol{\mu_{x}}}$ denote the mean vector of $\boldsymbol{x}$.

Similarly, the covariance matrices of the measured output $\boldsymbol{y}$,
state $\boldsymbol{x}_{\rm m}$, and input $\boldsymbol{u}$ can also be
obtained based on Eqs. (19a), (19b), and (19c) using the same approach:

$\begin{array}[]{l}{\boldsymbol{\Sigma}}_{\boldsymbol{y}}={{\boldsymbol{K}}_{{\boldsymbol{y}}2{\boldsymbol{x}}}}{\boldsymbol{\Sigma}}_{{{\boldsymbol{x}}}}{\boldsymbol{K}}_{{\boldsymbol{y}}2{\boldsymbol{x}}}^{\rm{T}}+{{\boldsymbol{K}}_{{\boldsymbol{y}}2{\boldsymbol{v}}}}{\boldsymbol{\Sigma}}_{\boldsymbol{v}}{\boldsymbol{K}}_{{\boldsymbol{y}}2{\boldsymbol{v}}}^{\rm{T}}\\\
{\boldsymbol{\Sigma}}_{{{\boldsymbol{x}}_{\rm{m}}}}={{\boldsymbol{K}}_{{\boldsymbol{x}}2{\boldsymbol{y}}}}{\boldsymbol{\Sigma}}_{\boldsymbol{y}}{\boldsymbol{K}}_{{\boldsymbol{x}}2{\boldsymbol{y}}}^{\rm{T}}\\\
{\boldsymbol{\Sigma}}_{\boldsymbol{u}}={{\boldsymbol{K}}^{{\rm{MPC}}}}{\boldsymbol{\Sigma}}_{{{\boldsymbol{x}}_{\rm{m}}}}{\left({{{\boldsymbol{K}}^{{\rm{MPC}}}}}\right)^{\rm{T}}}\end{array}$
(27)

where ${{\boldsymbol{K}}_{{\bullet}2{\star}}}$ denotes the corresponding
sensitivity matrix of $\bullet$ with respect to $\star$.

###  4.4 Dynamic performance assessment

In this section, we use the above results to obtain the distribution functions
of states and the closed-loop performance functions which contains only
states. Then, we integrate the distribution functions and the closed-loop
performance functions to evaluate the dynamic economic and tracking
performance of the controlled systems.

####  4.4.1 Distribution function of states

Each state can be supposed to satisfy normal distribution and all states
satisfy joint normal distribution as (28), if only Gaussian process and
measurement noise exists in a well controlled system. Note that the measured
${\boldsymbol{x}}_{\rm m}$ is be used to evaluate the performance in practical
engineering instead of $\boldsymbol{x}$ due to its unpredictability.

${\boldsymbol{x}}_{\rm m}\sim{\rm
N}({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$ (28)

where $\boldsymbol{\mu}$ is the mean vector and $\boldsymbol{\Sigma}$ is the
covariance matrix.

Therefore, the joint probability density function of states can be expressed
as follows:

$p({\boldsymbol{x}}_{\rm
m}):=\psi({\boldsymbol{\mu}},{\boldsymbol{\Sigma}},{\boldsymbol{\rho}})$ (29)

where $\boldsymbol{\rho}$ denotes the correlation coefficients between
different states and can be calculated based on the covariance matrix
$\boldsymbol{\Sigma}$, $\psi(\cdot)$ is the joint normal probability density
function.

####  4.4.2 Closed-loop performance function and surface

If we substitute the controller gain ${\boldsymbol{K}}^{\rm{MPC}}$ into the
performance indices $J_{\rm ep}(k)$ and $J_{\rm tp}(k)$ shown in (4), we can
get the closed-loop economic and tracking performance functions as follows
which eliminate inputs and contain states only.

$\displaystyle{J_{{\rm{ep}}}}({\boldsymbol{x}}_{\rm m})$
$\displaystyle:={J_{{\rm{ep}}}}({\boldsymbol{x}}_{\rm
m},{{\boldsymbol{u}}_{\rm{s}}}+{{\boldsymbol{K}}^{{\rm{MPC}}}}({\boldsymbol{x}}_{\rm
m}-{{\boldsymbol{x}}_{\rm{s}}}))$ (30a)
$\displaystyle{J_{{\rm{tp}}}}({\boldsymbol{x}}_{\rm m})$
$\displaystyle:={J_{{\rm{tp}}}}({\boldsymbol{x}}_{\rm
m},{{\boldsymbol{u}}_{\rm{s}}}+{{\boldsymbol{K}}^{{\rm{MPC}}}}({\boldsymbol{x}}_{\rm
m}-{{\boldsymbol{x}}_{\rm{s}}}))$ (30b)

The above closed-loop performance functions can be used to evaluate the
tracking and economic performance of controlled systems. For better understand
their influence, we can plot the dynamic tracking and economic performance
surfaces around the given steady state using them, which could intuitively
show the sensitivity degree of the performance indices with respect to states.
The calculated closed-loop performance functions and plotted surfaces are also
helpful for designing a better objective function with lower sensitivity with
respect to noise. Figure 3 shows a simple example of the dynamic performance
surfaces, where Figure 3(a) denotes the dynamic economic performance surface
and Figure 3(b) is the dynamic tracking performance surface.

(a) Economic performance surface

(b) Tracking performance surface

Figure 3: Example of dynamic economic and tracking performance surfaces

####  4.4.3 Average performance assessment

Finally, the approximate average tracking and economic performance can be
evaluated through integrating the joint normal probability density function
(29) and the closed-loop performance functions (30), which are shown as Eq.
(31).

$\begin{array}[]{l}{\bar{J}_{{\rm{ep}}}}\approx\int\limits_{{{\boldsymbol{x}}_{\rm{s}}}-(3\sim
5){\boldsymbol{\sigma}}}^{{{\boldsymbol{x}}_{\rm{s}}}+(3\sim
5){\boldsymbol{\sigma}}}{p({\boldsymbol{x}}_{\rm
m}){J_{{\rm{ep}}}}({\boldsymbol{x}}_{\rm m})}{\rm{d}}{\boldsymbol{x}_{\rm
m}}\\\
{\bar{J}_{{\rm{tp}}}}\approx\int\limits_{{{\boldsymbol{x}}_{\rm{s}}}-(3\sim
5){\boldsymbol{\sigma}}}^{{{\boldsymbol{x}}_{\rm{s}}}+(3\sim
5){\boldsymbol{\sigma}}}{p({\boldsymbol{x}}_{\rm
m}){J_{{\rm{tp}}}}({\boldsymbol{x}}_{\rm m})}{\rm{d}}{\boldsymbol{x}_{\rm
m}}\end{array}$ (31)

where ${\boldsymbol{\sigma}}$ denotes the standard deviation vector of
${\boldsymbol{x}}_{\rm m}$.

Theoretically, the performance assessment will be more accurate if the
integral region is larger. However, the integral region cannot be too large
because the calculated controller gain ${\boldsymbol{K}}^{\rm{MPC}}$ only
works well around the given steady state. Thereby, the integral region can be
set as $(3\sim 5){\boldsymbol{\sigma}}$, i.e.
$\left[{{{\boldsymbol{x}}_{\rm{s}}}-(3\sim
5){\boldsymbol{\sigma}},{{\boldsymbol{x}}_{\rm{s}}}+(3\sim
5){\boldsymbol{\sigma}}}\right]$, which should be enough to cover more than
99.75% states for a well-controlled steady-state operating system.

##  5 Application to a chemical process example

In this section, we apply the proposed procedure to the CSTR shown in Figure 1
to illustrate the applicability and effectiveness of the proposed approach.
The model expression has been carefully explained in Section 3. Four different
cases are considered including Case 1 - $\boldsymbol{x}$ and $\boldsymbol{u}$
are both away from their boundaries, Case 2 - $\boldsymbol{x}$ locates on its
boundaries and $\boldsymbol{u}$ is free (away from its boundaries), Case 3 -
$\boldsymbol{u}$ locates on its boundaries and $\boldsymbol{x}$ is free, and
Case 4 - other conditions.

###  5.1 Case 1: $\boldsymbol{x}$ and $\boldsymbol{u}$ are both away from
their boundaries

In Case 1, we configure one scenario that the optimal $\boldsymbol{x}$ and
$\boldsymbol{u}$ are both away from their boundaries, which is very common and
desired in practical engineering. We are going to pre-evaluate the dynamic
economic and tracking performance of EMPC and tracking MPC based on the
proposed approach, which is helpful for us to choose an appropriate controller
to achieve better performance. Simulation configuration has been described in
detail by the first group experiment in Section 3.

####  5.1.1 Economic MPC

Based on the above configuration, the optimal economic steady state is
${{\boldsymbol{x}}_{\rm{s}}}=[1.1601~{}615.7373]^{\rm T}$. The controller gain
of EMPC obtained from sensitivity analysis is shown as below:

${{\boldsymbol{K}}^{{\rm{EMPC}}}}=\left[{\begin{array}[]{*{20}{c}}4.4428e\\!+\\!02&4.4232e\\!+\\!00\\\
1.9369e\\!+\\!07&1.8381e\\!+\\!05\end{array}}\right]$ (32)

Using EMPC and sensitivity analysis, we can get two groups input trajectories
as Figure 4(b)based on the same state trajectory shown in Figure 4(a), i.e.
the optimal dynamic input trajectory generated by EMPC and the calculated
input trajectory using sensitivity analysis. In Figure 4, “$\rm x_{s}$” and
“$\rm u_{s}$” denote the optimal state and input steady states, respectively;
“$\rm x$” is the actual state trajectory; “$\rm u_{mpc}$” denotes the optimal
input trajectory generated by EMPC; and “$\rm u_{sen}$” denotes the calculated
input trajectory using sensitivity analysis. It is shown that the “$\rm
u_{mpc}$” generated by EMPC and the “$\rm u_{sen}$” from sensitivity analysis
are basically the same. We can conclude that the controller gain obtained from
sensitivity analysis is accurate enough.

(a) State trajectories

(b) Input trajectories

Figure 4: State and input trajectories of EMPC in Case 1

Using the proposed approach, we can get the mean and variance of each
variable. Table 3 shows the corresponding results, where $\mu$ denotes the
mean, $\sigma^{2}$ is the variance, “Sim” denotes the results from
simulations, and “Cal” denotes the results calculated by sensitivity
evaluation. It is shown that the calculated results and the simulated results
are almost exactly the same.

Table 3: Mean and variance of each variable of EMPC in Case 1 |  | $x_{1}$ | $x_{2}$ | $y_{1}$ | $y_{2}$ | $u_{1}$ | $u_{2}$
---|---|---|---|---|---|---|---
$\mu$ | Sim | 1.1603 | 615.7342 | 1.1603 | 615.7334 | 2.7894e+02 | 1.2774e+07
Cal | 1.1601 | 615.7373 | 1.1601 | 615.7373 | 2.7885e+02 | 1.2769e+07
$\sigma^{2}$ | Sim | 8.2736e-04 | 2.4467e+00 | 9.1822e-04 | 2.5323e+00 | 5.5399e+01 | 1.1233e+11
Cal | 8.3498e-04 | 2.4718e+00 | 9.2498e-04 | 2.5618e+00 | 5.5615e+01 | 1.1276e+11

Figure 5 shows the dynamic economic and tracking performance surfaces around
the optimal steady state. It is shown that the optimal dynamic performance
point locates around the optimal steady state. However, the dynamic
performance surfaces may be not symmetric.

(a) Economic performance surface

(b) Tracking performance surface

Figure 5: Dynamic performance surfaces of EMPC in Case 1

Suppose the two states satisfy joint normal distribution and their joint
probability density function is as below:

$p({\boldsymbol{x}}_{\rm
m})={\left({2\pi{\sigma_{x_{1}}}{\sigma_{x_{2}}}\sqrt{1-{\rho^{2}}}}\right)^{-1}}e^{{-\frac{1}{{2(1-{\rho^{2}})}}\left({\frac{{{{({x_{1}}-{\mu_{x_{1}}})}^{2}}}}{{\sigma_{x_{1}}^{2}}}-\frac{{2\rho({x_{1}}-{\mu_{x_{1}}})({x_{2}}-{\mu_{x_{2}}})}}{{{\sigma_{x_{1}}}{\sigma_{x_{2}}}}}+\frac{{{{({x_{2}}-{\mu_{x_{2}}})}^{2}}}}{{\sigma_{x_{2}}^{2}}}}\right)}}$
(33)

where $\mu_{x_{i}}$ denotes the mean of $x_{i}$, $\rho$ is the correlation
coefficient.

Table 4 shows the average performance indices from simulations and
calculations. To avoid the uncertainties of experiments, the simulation length
is set very large to be $N_{\rm sim}=100000$. The calculated average
performance indices using the proposed method includes the results of
integrating $3\boldsymbol{\sigma}$, $4\boldsymbol{\sigma}$, and
$5\boldsymbol{\sigma}$ zones. It is shown that the calculated economic and
tracking performance results are the same as the results of simulations, which
reveals that the proposed approach works well in this situation.

Table 4: Calculated and simulated performance indices of EMPC in Case 1 |  | 3$\boldsymbol{\sigma}$ zone | 4$\boldsymbol{\sigma}$ zone | 5$\boldsymbol{\sigma}$ zone
---|---|---|---|---
$\bar{J}_{\rm ep}$ | Cal | 451.4201 | 453.9029 | 454.0342
Sim |  | 453.5738 |
$\bar{J}_{\rm tp}$ | Cal | 2.0345e-03 | 2.0831e-03 | 2.0475e-03
Sim |  | 2.0905e-03 |

####  5.1.2 Tracking MPC

The controller gain of tracking MPC obtained from sensitivity analysis is
shown as below:

${{\boldsymbol{K}}^{{\rm{TMPC}}}}=\left[{\begin{array}[]{*{20}{c}}2.8538e\\!+\\!02&2.8488e\\!+\\!00\\\
1.4720e\\!+\\!07&1.3795e\\!+\\!05\end{array}}\right]$ (34)

The same results and conclusions as EMPC’s can be attained from Figure 6,
Table 5, Figure 7, and Table 6. They further illustrate the effectiveness of
the proposed sensitivity-based dynamic performance assessment approach.

(a) State trajectories

(b) Input trajectories

Figure 6: State and input trajectories of tracking MPC in Case 1 Table 5: Mean and variance of each variable of tracking MPC in Case 1 |  | $x_{1}$ | $x_{2}$ | $y_{1}$ | $y_{2}$ | $u_{1}$ | $u_{2}$
---|---|---|---|---|---|---|---
$\mu$ | Sim | 1.1601 | 615.7488 | 1.1601 | 615.7480 | 2.7890e+02 | 1.2773e+07
Cal | 1.1601 | 615.7373 | 1.1601 | 615.7373 | 2.7885e+02 | 1.2769e+07
$\sigma^{2}$ | Sim | 5.0130e-04 | 1.5292e+00 | 5.9141e-04 | 1.6164e+00 | 1.9444e+01 | 5.4403e+10
Cal | 5.0824e-04 | 1.5499e+00 | 5.9824e-04 | 1.6399e+00 | 1.9513e+01 | 5.4636e+10

(a) Economic performance surface

(b) Tracking performance surface

Figure 7: Dynamic performance surfaces of tracking MPC in Case 1 Table 6: Calculated and simulated performance indices of tracking MPC in Case 1 |  | 3$\boldsymbol{\sigma}$ zone | 4$\boldsymbol{\sigma}$ zone | 5$\boldsymbol{\sigma}$ zone
---|---|---|---|---
$\bar{J}_{\rm ep}$ | Cal | 446.6845 | 448.9983 | 449.0709
Sim |  | 448.9241 |
$\bar{J}_{\rm tp}$ | Cal | 1.0010e-03 | 1.0287e-03 | 1.0040e-03
Sim |  | 1.0275e-03 |

####  5.1.3 Performance comparison

(1) Performance analysis of group 1

Compare Figure 5 with Figure 7, the economic performance surfaces of EMPC and
tracking MPC are very similar, while the tracking performance surfaces being
different in consideration of the values of z-axis. The possible reason should
be that EMPC adopts a nonsymmetric economic objective while tracking MPC using
a symmetric quadratic tracking objective so that tracking MPC will be more
sensitive to the variations of states. The differences can also be verified by
Table 3 and Table 5. Meanwhile, Table 3 and Table 5 show that the variances of
variables in EMPC and tracking MPC are different, which makes their economic
performance be different. Table 4 and Table 6 reveal that tracking MPC has 1%
better economic performance than EMPC in this group experiments.

(2) Performance comparison of group 1 and group 2

In group 2 of Section 3, the optimal economic steady state is
${\boldsymbol{x}}_{\rm s}=[2.7284~{}631.7380]^{\rm T}$ and the average
performance indices are shown in Table 7 and Table 8. In this group
experiments, the economic performance of EMPC is better than tracking MPC
(about 1%). Compare them with Table 4 and Table 6, it is shown that which
controller can achieve better performance is uncertain in practical processes.
However, the proposed dynamic performance assessment approach can provide
accurate pre-evaluation on the average performance of EMPC and tracking MPC,
which can be used to choose the most appropriate controller when engineers
design controllers for given systems.

Table 7: Calculated and simulated performance indices of EMPC in group 2 of Case 1 |  | 3$\boldsymbol{\sigma}$ zone | 4$\boldsymbol{\sigma}$ zone | 5$\boldsymbol{\sigma}$ zone
---|---|---|---|---
$\bar{J}_{\rm ep}$ | Cal | 197.9753 | 198.9069 | 198.9312
Sim |  | 199.8444 |
$\bar{J}_{\rm tp}$ | Cal | 6.7851e-03 | 7.0437e-03 | 7.0418e-03
Sim |  | 7.1578e-03 |
Table 8: Calculated and simulated performance indices of tracking MPC in group 2 of Case 1 |  | 3$\boldsymbol{\sigma}$ zone | 4$\boldsymbol{\sigma}$ zone | 5$\boldsymbol{\sigma}$ zone
---|---|---|---|---
$\bar{J}_{\rm ep}$ | Cal | 200.6039 | 201.6109 | 201.6389
Sim |  | 202.5862 |
$\bar{J}_{\rm tp}$ | Cal | 6.5879e-03 | 6.7543e-03 | 6.7537e-03
Sim |  | 6.7702e-03 |

###  5.2 Case 2: $\boldsymbol{x}$ locates on its boundaries and
$\boldsymbol{u}$ is free

In Case 2, we configure one scenario that $\boldsymbol{x}$ locates on its
boundaries and $\boldsymbol{u}$ is free (away from its boundaries). In this
scenario, violation of boundaries and offset of steady states are very common
in tracking MPC and EMPC, which will be presented in the follows. Based on the
proposed performance assessment approach, we can pre-evaluate the violation
degree and pre-set boundary moving or target moving to guarantee the stability
of controlled processes.

The economic objective is set as below:

${J_{\rm{ep}}}={{x_{1}}+{x_{2}}}$ (35)

The constraints considered in the experiments are shown as follows:

$1\leq{x_{1}}\leq 3,~{}500\leq{x_{2}}\leq 1000$ (36)

The noise added on states and measurements is shown as follows:

$\begin{array}[]{l}{w_{1}}\sim{\rm N}(0,0.0001),{v_{1}}\sim{\rm N}(0,0.1)\\\
{w_{2}}\sim{\rm N}(0,0.0004),{v_{2}}\sim{\rm N}(0,0.4)\end{array}$ (37)

####  5.2.1 Violation and offset

(1) Problem description

In the EMPC and tracking MPC controlled systems, when $\boldsymbol{x}$ is
close to its boundaries, the actual state values will violate their boundaries
in the presence of noise, which is harmful to the safety and stability of
controlled systems. Figure 8(a) gives a simple example to show this problem,
in which “$\rm x_{lb}$” denotes the lower boundary of $\boldsymbol{x}$. Except
this, there is another problem called offset shown in Figure 8(b), which means
that the actual states will diverge from the given steady state
${{\boldsymbol{x}}_{\rm{s}}}$. It may be caused by the punishment of both
${{\boldsymbol{x}}-{{\boldsymbol{x}}_{\rm{s}}}}$ and
${{\boldsymbol{u}}-{{\boldsymbol{u}}_{\rm{s}}}}$ in tracking MPC. The
violation problem exists in both EMPC and tracking MPC, while the offset
problem existing in only tracking MPC for the reason that EMPC is one kind of
optimal dynamic performance searching not a given target tracking.

(a) Violation phenomenon

(b) Offset phenomenon

Figure 8: Violation and offset phenomenons of state variables

(2) Reasonable solutions

① Move the boundaries (Boundary moving)

If there is not offset, moving the boundaries of $\boldsymbol{x}$ (i.e.
backoff) will be a good method to avoid the violation while ensuring minimum
economy loss. Therefore, it is a suitable way to deal with the violation in
EMPC. However, this method works not so well with tracking MPC because the
existence of offset in tracking MPC will make the states diverge from the
adjusted boundaries again and lead to more economy loss.

② Move the optimal steady state (Target moving)

When the optimal steady state locates on the boundaries, there may co-exist
violation and offset problems in tracking MPC. In this situation, moving the
optimal steady state away from its boundaries will be a good method to solve
these problems. After target moving, the reformulated optimization problem
will be the same as Case 1 in Section 5.1. Target moving is a suitable way to
deal with the violation and offset in tracking MPC but not suitable for EMPC
due to the lack of given steady states.

####  5.2.2 Economic MPC

The optimal steady state locates on the original lower boundary shown as
below:

${{\boldsymbol{x}}_{\rm{s}}}={\boldsymbol{x}}_{{\rm{ori}}}^{{\rm{lb}}}=[1~{}500]^{\rm
T}$ (38)

The controller gain of EMPC obtained from sensitivity analysis is as below:

${{\boldsymbol{K}}^{{\rm{EMPC}}}}=\left[{\begin{array}[]{*{20}{c}}-9.5472e\\!+\\!01&4.8640e\\!-\\!01\\\
-5.5836e\\!+\\!06&-7.0123e\\!+\\!04\end{array}}\right]$ (39)

The variances of states calculated by the proposed method are shown as below:

${\boldsymbol{\sigma}}^{2}=[2.8039e\\!-\\!04~{}~{}~{}6.7546e\\!-\\!01]^{\rm
T}$ (40)

The new lower boundary obtained through moving $3\boldsymbol{\sigma}$ distance
is as below:

$\begin{array}[]{l}{\boldsymbol{x}}_{{\rm{new}}}^{{\rm{lb}}}={\boldsymbol{x}}_{{\rm{ori}}}^{{\rm{lb}}}+3{\boldsymbol{\sigma}}=[1.0502~{}502.4656]^{\rm
T}\end{array}$ (41)

The simulation results after boundary moving are shown in Figure 9, in which
“$\rm x_{lb,ori}$” and “$\rm x_{lb,new}$” denote the original and new lower
boundaries of $\boldsymbol{x}$ respectively. To show the violation of
$\boldsymbol{x}$, we plot more than two hours results in Figure 9(a). It is
shown the percentages of the data that violates the original lower boundary
are $[0.1036\%~{}0.1978\%]^{\rm T}$, which are close to the theoretical values
$[0.125\%~{}0.125\%]^{\rm T}$. Therefore, the sensitivity-based boundary
moving works well for reducing the risk of violation in EMPC.

(a) State trajectories

(b) Input trajectories

Figure 9: Simulation results after boundary moving of EMPC in Case 2

####  5.2.3 Tracking MPC

The original given steady state and lower boundary are shown as below:

${{\boldsymbol{x}}_{{\rm{s,ori}}}}={{\boldsymbol{x}}^{{\rm{lb}}}}=[1~{}500]^{\rm
T}$ (42)

The reformulated problem after target moving will be a “Case 1” problem shown
in Section 5.1. The actual variances of states at the new adjusted steady
state is unknown because the adjusted steady state cannot be determined before
knowing its variances. However, the variances at the original steady state can
be used to replace the variances at the new adjusted steady state because the
adjusted steady state will not be far away from the original steady state. It
is worthy noting that the active constraint flags should be set as
${{\boldsymbol{I}}_{\boldsymbol{x}}}={\boldsymbol{0}},{{\boldsymbol{I}}_{\boldsymbol{u}}}={\boldsymbol{0}}$
when we calculate the controller gain of tracking MPC.

The controller gain of tracking MPC obtained from sensitivity analysis is as
below:

${{\boldsymbol{K}}^{{\rm{TMPC}}}}=\left[{\begin{array}[]{*{20}{c}}1.3458e\\!-\\!01&1.5024e\\!-\\!01\\\
-4.1946e\\!+\\!04&-1.2903e\\!+\\!03\end{array}}\right]$ (43)

The variances of states calculated by sensitivity analysis are as follows:

${\boldsymbol{\sigma}}^{2}=[2.4150e\\!-\\!04~{}~{}~{}2.0429e\\!+\\!00]^{\rm
T}$ (44)

The new steady state obtained through moving $3\boldsymbol{\sigma}$ distance
is as below:

${{\boldsymbol{x}}_{{\rm{s,new}}}}={{\boldsymbol{x}}_{{\rm{s,ori}}}}+3{\boldsymbol{\sigma}}=[1.0466~{}504.2879]^{\rm
T}$ (45)

The simulation results after target moving are shown in Figure 10, in which
“$\rm x_{s,ori}$” and “$\rm x_{s,new}$” denote the original and new steady
states of $\boldsymbol{x}$ respectively. It is shown the percentages of the
data that violates the original lower boundary (${\boldsymbol{x}}^{\rm
lb}={\boldsymbol{x}}_{\rm s,ori}$) are $[0.0552\%~{}0.0126\%]^{\rm T}$, which
shows the proposed sensitivity-based target moving guarantees as little data
as possible crosses the original boundaries. The differences between
${\boldsymbol{u}}_{\rm mpc}$ and ${\boldsymbol{u}}_{\rm sen}$ may result from
the controller gain approximation for using ${\boldsymbol{x}}_{\rm s,ori}$
instead of ${\boldsymbol{x}}_{\rm s,new}$. In conclusion, the sensitivity-
based target moving works well for reducing the risk of violation and offset
in tracking MPC.

(a) State trajectories

(b) Input trajectories

Figure 10: Simulation results after target moving of TMPC in Case 2

####  5.2.4 Performance assessment

Table 9 shows the economic performance loss resulted from boundary moving and
target moving. The economy loss of EMPC is about 0.50%. Both the simulation
results and the proposed assessment method give the same conclusion. The
economy loss of tracking MPC is about 0.87%. The results illustrate that the
proposed method could provide accurate pre-evaluation about the possible
performance loss when boundary moving and target moving are applied in EMPC
and tracking MPC.

Table 9: Economic performance loss resulted from boundary moving and target moving |  | Ideal $\bar{J}_{\rm{ep}}$ | Actual $\bar{J}_{\rm{ep}}$ | Loss of $\bar{J}_{\rm{ep}}$
---|---|---|---|---
Economic MPC | Cal | 501 | 503.5158 | 2.5158 (0.50%)
Sim | 501 | 503.5214 | 2.5214 (0.50%)
Tracking MPC | Cal | 501 | 505.3345 | 4.3345 (0.87%)
Sim | 501 | 505.2858 | 4.2858 (0.86%)

Another group of experiments are conducted as a comparison, in which the
original lower boundary of $\boldsymbol{x}$ is changed to
${{\boldsymbol{x}}^{{\rm{lb}}}}=[1~{}600]^{\rm T}$. Table 10 shows the
economic performance loss in the new experiments. Compare Table 10 with Table
9, it is shown that which controller can achieve better performance is
uncertain after boundary moving and target moving. But the proposed approach
can tell us how much we need moving the boundaries or targets for stable
operation and which controller can achieve less performance loss.

Table 10: Economic performance loss resulted from boundary moving and target moving in group 2 of Case 2 |  | Ideal $\bar{J}_{\rm{ep}}$ | Actual $\bar{J}_{\rm{ep}}$ | Loss of $\bar{J}_{\rm{ep}}$
---|---|---|---|---
Economic MPC | Cal | 601 | 608.0289 | 7.0289 (1.17%)
Sim | 601 | 608.0316 | 7.0316 (1.17%)
Tracking MPC | Cal | 601 | 605.5987 | 4.5987 (0.77%)
Sim | 601 | 605.5677 | 4.5677 (0.76%)

###  5.3 Case 3: $\boldsymbol{x}$ is free and $\boldsymbol{u}$ locates on its
boundaries

In Case 3, we configure one scenario that $\boldsymbol{x}$ is free (away from
its boundaries) and $\boldsymbol{u}$ locates on its boundaries. In this
scenario, Gaussian noise will influence the controller gain of tracking MPC
due to the deficiency of control freedom caused by the active constraints of
$\boldsymbol{u}$. This is not desired and expected to be avoided in process
control.

The economic objective is shown as below [20]:

$J_{\rm ep}=-1700000{u_{1}}({C_{A0}}-{x_{1}})+{u_{2}}$ (46)

The constraints considered in the experiments are shown as follows:

$0\leq{u_{1}}\leq 10,~{}-200000\leq{u_{2}}\leq 200000$ (47)

The noise added on states and measurements is shown as follows:

$\begin{array}[]{l}{w_{1}}\sim{\rm N}(0,0.00001),{v_{1}}\sim{\rm N}(0,0.01)\\\
{w_{2}}\sim{\rm N}(0,0.00004),{v_{2}}\sim{\rm N}(0,0.04)\end{array}$ (48)

####  5.3.1 Economic MPC

The optimal steady states are
${{\boldsymbol{x}}_{\rm{s}}}=[0.5096~{}536.7485]^{\rm T}$ and
${{\boldsymbol{u}}_{\rm{s}}}=[10~{}200000]^{\rm T}$. The controller gain
obtained from sensitivity analysis is shown in Eq. (49), in which the elements
are all 0. It means the optimal inputs will not change and their constraints
are active as long as the process operates around the given steady state. It
can be explained by that the inputs are economic decision variables in the
optimization of EMPC.

${{\boldsymbol{K}}^{{\rm{EMPC}}}}=\left[{\begin{array}[]{*{20}{c}}0&0\\\
0&0\end{array}}\right]$ (49)

The mean and variance of each variable and the average performance indices are
shown in Table 11 and Table 12, respectively. They illustrate that the
proposed sensitivity-based dynamic performance assessment approach works well
in the EMPC controlled system.

Table 11: Mean and variance of each variable of EMPC in Case 3 |  | $x_{1}$ | $x_{2}$ | $y_{1}$ | $y_{2}$ | $u_{1}$ | $u_{2}$
---|---|---|---|---|---|---|---
$\mu$ | Sim | 5.1004e-01 | 5.3662e+02 | 5.1006e-01 | 5.3662e+02 | 1.e+01 | 2.e+05
Cal | 5.0955e-01 | 5.3675e+02 | 5.0955e-01 | 5.3675e+02 | 1.e+01 | 2.e+05
$\sigma^{2}$ | Sim | 3.7723e-05 | 9.0709e-01 | 7.6350e-05 | 9.4428e-01 | 0 | 0
Cal | 3.6951e-05 | 9.3944e-01 | 7.6951e-05 | 9.7944e-01 | 0 | 0
Table 12: Calculated and simulated performance indices of EMPC in Case 3 |  | 3$\boldsymbol{\sigma}$ zone | 4$\boldsymbol{\sigma}$ zone | 5$\boldsymbol{\sigma}$ zone
---|---|---|---|---
$\bar{J}_{\rm{ep}}$ | Cal | -5.0370e+07 | -5.0631e+07 | -5.0638e+07
Sim |  | -5.0639e+07 |
$\bar{J}_{\rm{tp}}$ | Cal | 2.9194e-04 | 2.9942e-04 | 2.9732e-04
Sim |  | 3.0030e-04 |

####  5.3.2 Tracking MPC

In this case, the controller gain is varying with respect to the positions of
states because ${\boldsymbol{u}}_{\rm s}$ locates on the boundaries. When the
optimal $\boldsymbol{u}$ needs acting toward the boundaries, the actual inputs
will be restricted inside the boundaries. On the contrary, there is no
problem. Therefore, we calculate the controller gains with active constraints
${\boldsymbol{I_{u}}}=\boldsymbol{I}$ and inactive constraints
${\boldsymbol{I_{u}}}=\boldsymbol{0}$, respectively:

${\boldsymbol{K}}_{{{\boldsymbol{I}}_{\boldsymbol{u}}}={\boldsymbol{I}}}^{{\rm{TMPC}}}=\left[{\begin{array}[]{*{20}{c}}0&0\\\
0&0\end{array}}\right]$ (50)
${\boldsymbol{K}}_{{{\boldsymbol{I}}_{\boldsymbol{u}}}={\boldsymbol{0}}}^{{\rm{TMPC}}}=\left[{\begin{array}[]{*{20}{c}}-2.0092e\\!+\\!00&5.6354e\\!-\\!02\\\
-2.1417e\\!+\\!04&-5.7478e\\!+\\!02\end{array}}\right]$ (51)

It is impossible to know exactly the controller gain in this case, so here we
just discuss and compare the actual actions of tracking MPC with the predicted
results using sensitivity analysis. Figure 11 shows the trajectories of states
and inputs, in which “$\rm u_{sen,k0}$” denotes the input trajectory
calculated using Eq. (50) and “$\rm u_{sen,k1}$” denotes the results from Eq.
(51). The possible trajectory and tendency of ${\boldsymbol{u}}_{\rm mpc}$ can
be roughly predicted based on the results of sensitivity analysis. Some
conclusions about the ideal economic performance, possible tendency of
$\boldsymbol{u}$, and approximate economic loss can be roughly analyzed. It
may be helpful for the engineers to understand the actual actions of inputs in
the scenario. Actually, the situation that $\boldsymbol{u}$ locates on its
boundaries should be avoided in tracking MPC controlled systems in
consideration of the deficiency of control freedom. If target moving is
applied, Case 3 will convert to be a “Case 1” problem and can be analyzed
using the discussion in Section 5.1.

(a) State trajectories

(b) Input trajectories

Figure 11: State and input trajectories of tracking MPC in Case 3

###  5.4 Case 4: Other conditions

One of the largest differences between tracking MPC and EMPC is that the
former is one kind passive target tracking which may lead to offset and other
tracking problems, while the latter is a dynamic autonomous economy searching
which will not result in offset problem for the overlook of given targets.
Therefore, the proposed sensitivity-based dynamic performance assessment
method works well with all EMPC controlled systems as long as the process and
measurement noise is not too large.

Figure 12: Different conditions of tracking MPC

However in tracking MPC controlled systems, different conditions need to be
distinguished to make the analysis more reliable. Figure 12 gives a complete
view on all possible conditions in tracking MPC. Case 1 - Case 3 have been
discussed in detail in this work. Case 4, “Mixed active constraints”, denotes
the condition that some states and inputs are away from their boundaries and
some locate on their boundaries. The case can be converted to be a new case
problem discussed above through target moving.

##  6 Conclusion

In this work, a sensitivity-based dynamic performance assessment method is
proposed to analyze the tracking and economic performance of tracking MPC and
EMPC in the presence of constraints and Gaussian noise. It can provide
effective guidance on which controller can achieve better performance for a
given system and is very helpful for controller design. Extensive experiments
on a CSTR process demonstrate the effectiveness of the proposed approach. In
the future work, more concentrations will be focused on the dynamic
performance assessment of MPC with colored noise and disturbance.

## References

  * [1] Joel AE Andersson, Joris Gillis, Greg Horn, James B Rawlings, and Moritz Diehl.  Casadi: a software framework for nonlinear optimization and optimal control.  Mathematical Programming Computation, 11(1):1–36, 2019.
  * [2] Joel AE Andersson and James B Rawlings.  Sensitivity analysis for nonlinear programming in casadi.  IFAC-PapersOnLine, 51(20):331–336, 2018.
  * [3] David Angeli, Rishi Amrit, and James B Rawlings.  On average performance and stability of economic model predictive control.  IEEE transactions on automatic control, 57(7):1615–1626, 2011.
  * [4] Margret Bauer and Ian K Craig.  Economic assessment of advanced process control–a survey and framework.  Journal of process control, 18(1):2–18, 2008.
  * [5] Lorenz T Biegler.  A survey on sensitivity-based nonlinear model predictive control.  IFAC Proceedings Volumes, 46(32):499–510, 2013.
  * [6] Viviane Botelho, Jorge Otavio Trierweiler, Marcelo Farenzena, and Ricardo Duraiski.  Perspectives and challenges in performance assessment of model predictive control.  The Canadian Journal of Chemical Engineering, 94(7):1225–1241, 2016\.
  * [7] Jinghan Cui, Su Liu, Jinfeng Liu, and Xiangjie Liu.  A comparative study of mpc and economic mpc of wind energy conversion systems.  Energies, 11(11):3127, 2018.
  * [8] Benjamin Decardi-Nelson, Su Liu, and Jinfeng Liu.  Improving flexibility and energy efficiency of post-combustion co2 capture plants using economic model predictive control.  Processes, 6(9):135, 2018.
  * [9] Matthew Ellis, Helen Durand, and Panagiotis D Christofides.  A tutorial review of economic model predictive control methods.  Journal of Process Control, 24(8):1156–1178, 2014.
  * [10] Matthew Ellis, Jinfeng Liu, and Panagiotis D Christofides.  Economic model predictive control, volume 2.  Springer, 2017.
  * [11] AH Evers.  Sensitivity analysis in dynamic optimization.  Journal of optimization theory and applications, 32(1):17–37, 1980\.
  * [12] Mohsen Heidarinejad, Jinfeng Liu, and Panagiotis D Christofides.  Economic model predictive control of nonlinear process systems using lyapunov techniques.  AIChE Journal, 58(3):855–870, 2012.
  * [13] Michael A Henson and Dale E Seborg.  Nonlinear process control.  Prentice Hall PTR Upper Saddle River, New Jersey, 1997.
  * [14] Biao Huang and Sirish L Shah.  Performance assessment of control loops: theory and applications.  Springer Science & Business Media, 1999.
  * [15] Rui Huang, Victor M Zavala, and Lorenz T Biegler.  Advanced step nonlinear model predictive control for air separation units.  Journal of Process Control, 19(4):678–685, 2009.
  * [16] Mohieddine Jelali.  An overview of control performance assessment technology and industrial applications.  Control engineering practice, 14(5):441–466, 2006.
  * [17] Mohieddine Jelali.  Control performance management in industrial automation: assessment, diagnosis and improvement of control loop performance.  Springer Science & Business Media, 2012.
  * [18] Steven G Krantz and Harold R Parks.  The implicit function theorem: history, theory, and applications.  Springer Science & Business Media, 2012.
  * [19] Kwan Ho Lee, Biao Huang, and Edgar C Tamayo.  Sensitivity analysis for selective constraint and variability tuning in performance assessment of industrial mpc.  Control engineering practice, 16(10):1195–1215, 2008.
  * [20] Xinchun Li, Liqin Zhang, Makoto Nakaya, and Azusa Takenaka.  Application of economic mpc to a cstr process.  In 2016 IEEE Advanced Information Management, Communicates, Electronic and Automation Control Conference (IMCEC), pages 685–690. IEEE, 2016\.
  * [21] Su Liu and Jinfeng Liu.  Economic model predictive control with extended horizon.  Automatica, 73:180–192, 2016.
  * [22] Rodrigo López-Negrete and Lorenz T Biegler.  A moving horizon estimator for processes with multi-rate measurements: A nonlinear programming sensitivity approach.  Journal of Process Control, 22(4):677–688, 2012.
  * [23] Qiugang Lu, Michael G Forbes, Ratna Bhushan Gopaluni, Philip D Loewen, Johan U Backström, and Guy A Dumont.  Performance assessment of cross-directional control for paper machines.  IEEE Transactions on Control Systems Technology, 25(1):208–221, 2016.
  * [24] Hans Pirnay, Rodrigo López-Negrete, and Lorenz T Biegler.  Optimal sensitivity based on ipopt.  Mathematical Programming Computation, 4(4):307–331, 2012.
  * [25] James B Rawlings, David Angeli, and Cuyler N Bates.  Fundamentals of economic model predictive control.  In 2012 IEEE 51st IEEE conference on decision and control (CDC), pages 3851–3861. IEEE, 2012.
  * [26] Andreas Wächter and Lorenz T Biegler.  On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming.  Mathematical programming, 106(1):25–57, 2006.
  * [27] Xiaorui Wang, Biao Huang, and Tongwen Chen.  Multirate minimum variance control design and control performance assessment: A data-driven subspace approach.  IEEE Transactions on Control Systems Technology, 15(1):65–74, 2006\.
  * [28] Donghui Wei, Ian K Craig, and Margret Bauer.  Multivariate economic performance assessment of an mpc controlled electric arc furnace.  ISA transactions, 46(3):429–436, 2007.
  * [29] Fangwei Xu, Biao Huang, and Seyi Akande.  Performance assessment of model pedictive control for variability and constraint tuning.  Industrial & engineering chemistry research, 46(4):1208–1219, 2007\.
  * [30] Yi Zhang, Benjamin Decardi-Nelson, Jianbang Liu, Jiong Shen, and Jinfeng Liu.  Zone economic model predictive control of a coal-fired boiler-turbine generating system.  Chemical Engineering Research and Design, 153:246–256, 2020.