Economic Machine-Learning-Based Predictive Control of Nonlinear Systems

Abstract

In this work, a Lyapunov-based economic model predictive control (LEMPC) method is developed to address economic optimality and closed-loop stability of nonlinear systems using machine learning-based models to make predictions. Specifically, an ensemble of recurrent neural network (RNN) models via a k-fold cross validation is first developed to capture process dynamics in an operating region. Then, the LEMPC using an RNN ensemble is designed to maintain the closed-loop state in a stability region and optimize process economic benefits simultaneously. Parallel computing is employed to improve computational efficiency of real-time implementation of LEMPC with an RNN ensemble. The proposed machine-learning-based LEMPC method is demonstrated using a nonlinear chemical process example.

1. Introduction

Economic model predictive control (EMPC) that integrates optimization of process economics with process control has attracted significant attention over the last decade (e.g., [1,2,3]). Unlike the traditional two-layered real-time optimization (RTO) that continuously optimizes process set-points in the upper level and applies a tracking model predictive controller or a proportional-integral-derivative controller in the lower layer, EMPC operates the system in a time-varying fashion (off steady-state) by dynamically optimizing an economic cost function accounting for stability constraints in one layer. Many research works have been done on EMPC in terms of closed-loop stability, economic optimality, model uncertainty and process operational safety (e.g., [4,5,6]). However, to achieve a desired closed-loop performance under EMPC, it is acknowledged that a reliable process model that can predict the future state evolution accurately is a key requirement. To that end, multiple linear state-space models have been used within EMPC to capture the nonlinear dynamics in the operating region in [7]. Additionally, a well-conditioned polynomial nonlinear state-space model has been used in [8] to approximate the actual (first-principles) nonlinear process model (normally substituting of the actual nonlinear process) for EMPC. Given the recent interest in machine learning modeling and particularly its usage in solving nonlinear regression problems ([9,10]), in this work, we will derive a nonlinear data-driven process model for EMPC by taking advantage of machine learning techniques.

Machine learning has been widely used in a variety of engineering problems such as classification, pattern recognition, and numerical prediction. Among machine learning methods, neural networks have shown great success in regression problems. With the rapid development of neural network libraries and application programming interfaces e.g., Tensorflow [11], Caffe [12], Keras [13], etc.), building a neural network has presently become intuitive and straightforward. Additionally, neural network-based modeling outperforms other data-driven modeling methods due to its ability to implicitly detect complex nonlinearities and the availability of multiple training algorithms ([14]). Neural network has undergone several phases in which different theoretical approaches led to increasing understanding of its power and limitations. For example, in the 1980s, the first RNN (i.e., Hopfield networks) was created for pattern recognition ([15]). At that time, implementing a neural network with many layers for complex systems was a technical challenge [16]. Subsequently, deep learning was proposed to alleviate the problems of vanishing gradient in RNNs, which provided a new perspective to the machine learning field [17]. Compared to feedfoward neural networks, RNN is able to memorize previous information in the network and thus, can be applied to model a nonlinear dynamical system from time-series data. However, since the majority of research is focused on validating the identified system, it is not clear how RNN models can be incorporated in EMPC with guaranteed closed-loop stability.

In addition, ensemble learning, a machine learning paradigm that trains multiple models for the same problem, is gaining increasing attention in recent years ([18,19]). By taking advantage of multiple learning algorithms or various training datasets, ensemble learning provides better predictive performance than a single machine learning model. However, the consequent challenge for implementations of ensemble regression models in process control is how to solve the ensemble of RNNs efficiently and in a timely fashion. Although reducing the computation time of neural network training process is always an important research topic in machine learning community, process engineers are more concerned about the real-time implementation of neural networks in process control systems. To address computational implementation issues, parallel computing ([20]) is employed to solve multiple RNN models simultaneously in EMPC.

Motivated by the above, a Lyapunov-based EMPC method using an ensemble of RNN models is developed in this work. Specifically, in Section 2, the notation, the class of nonlinear systems considered and the stabilizability assumptions are introduced. In Section 3, the class of RNN models and the learning algorithm are provided. Then, ensemble learning is used to obtain multiple RNN models based on a k-fold cross validation, while parallel computing is used to improve computational efficiency. Subsequently, in Section 4, the Lyapunov-based EMPC that incorporates an RNN ensemble is formulated to guarantee the boundedness of the state in the closed-loop stability region accounting for sufficiently small modeling error, sampling time and bounded disturbances. In the last section, the proposed LEMPC method is applied to a chemical process example to demonstrate its guaranteed closed-loop stability and economic optimality properties.

2. Preliminaries

2.1. Notation

Throughout the manuscript, the operator $\left|·\right|$ denotes the Euclidean norm of a vector. $x^T$ is used to denote the transpose of x. The function $f(·)$ is of class ${\mathcal{C}}^1$ if it is continuously differentiable in its domain. The notation $L_{f}V(x)$ denotes the standard Lie derivative $L_{f}V(x):=\frac{\partial V(x)}{\partial x}f(x)$. “\” is used to denote the set subtraction, i.e., $A\backslash B:=\{x\in {\mathbf{R}}^n|x\in A,x\notin B\}$. A continuous function $\alpha :[0,a)\to [0,\infty )$ is said to belong to class $\mathcal{K}$ if $\alpha (0)=0$ and it is strictly increasing.

2.2. Class of Systems

We consider the class of continuous-time nonlinear systems that is described by the following system of first-order nonlinear ordinary differential equations:

$$\dot{x}=F(x,u,w):=f(x)+g(x)u+h(x)w,x(t_0)=x_0$$

(1)

where the state vector is $x\in {\mathbf{R}}^n$, the manipulated input vector is $u\in {\mathbf{R}}^m$, and the disturbance vector is $w\in W$ with $W:=\{w\in {\mathbf{R}}^q||w|\le w_m,w_m\ge 0\}$. The control action constraint is defined by $u\in U:=\{u_i^{min}\le u_i\le u_i^{max},i=1,\dots ,m\}\subset {\mathbf{R}}^m$. It is assumed that $f(·):{\mathbf{R}}^n\to {\mathbf{R}}^n$, $g(·):{\mathbf{R}}^n\to {\mathbf{R}}^{n\times m}$, and $h(·):{\mathbf{R}}^n\to {\mathbf{R}}^{n\times q}$ are sufficiently smooth vector and matrix functions. Additionally, it is assumed that the initial time $t_0$ is zero (i.e., $t_0=0$) and $f(0)=0$. Therefore, a steady-state of the nominal system of Equation (1) (i.e., $w(t)\equiv 0$) is at the origin.

A stabilizing control law $u=\mathsf{\Phi }(x)\in U$ that renders the origin of the nominal system of Equation (1) (i.e., $w(t)\equiv 0$) exponentially stable is assumed to exist, in the sense that there exists a ${\mathcal{C}}^1$ Control Lyapunov function $V(x)$ such that the following inequalities hold for all x in an open neighborhood D around the origin:

$$c_1{|x|}^2\le V(x)\le c_2{|x|}^2,$$

(2a)

$$\frac{\partial V(x)}{\partial x}F(x,\mathsf{\Phi }(x),0)\le -c_3{|x|}^2,$$

(2b)

$$\left|\frac{\partial V(x)}{\partial x}\right|\le c_4|x|$$

(2c)

where $c_1$, $c_2$, $c_3$ and $c_4$ are positive constants. $F(x,u,w)$ is the nonlinear system of Equation (1). Based on the universal Sontag control law ([21]), a candidate controller $\mathsf{\Phi }(x)$ is given by the saturated control law accounting for the input constraint $u\in U$.

The stability region for the closed-loop system of Equation (1) is characterized as a level set of Lyapunov function $V(x)$ within D where the time-derivative of the Lyapunov function V is negative. Specifically, based on Equation (2), we can first find a region ${\varphi }_u$ where $\dot{V}(x)$ is rendered negative under the controller $u=\mathsf{\Phi }(x)\in U$: ${\varphi }_u:=\{x\in {\mathbf{R}}^n|\dot{V}(x)=L_{f}V+L_{g}Vu<-kV(x),u=\mathsf{\Phi }(x)\in U\}\cup \{0\}$, with $k>0$. Then the closed-loop stability region ${\Omega }_{\rho }$ is characterized as a subset in ${\varphi }_u$: ${\Omega }_{\rho }:=\{x\in {\varphi }_u|V(x)\le \rho \}$, where $\rho >0$. Since ${\Omega }_{\rho }\subset {\varphi }_u$, it follows that ${\Omega }_{\rho }$ is a forward invariant set. Additionally, the Lipschitz property of $F(x,u,w)$ combined with the bounds on u and w implies that there exist positive constants M, $L_x,L_w,L_x^{\prime },L_w^{\prime }$ such that the following inequalities hold for all $x,x^{\prime }\in D,u\in U$, and $w\in W$:

$$|F(x,u,w)|\le M$$

(3a)

$$|F(x,u,w)-F(x^{\prime },u,0)|\le L_x|x-x^{\prime }|+L_w|w|$$

(3b)

$$\begin{array}{c}\left|\frac{\partial V(x)}{\partial x}F(x,u,w)-\frac{\partial V(x^{\prime })}{\partial x}F(x^{\prime },u,0)\right|\\ \le L_x^{{}^{\prime }}|x-x^{\prime }|+L_w^{{}^{\prime }}|w|\end{array}$$

(3c)

Remark 1.

Since closed-loop stability is achieved for the nonlinear system of Equation (1) under the controller $u=\mathsf{\Phi }(x)\in U$ for all initial conditions $x_0\in {\Omega }_{\rho }$ ([4]), the dataset for training recurrent neural networks will be generated within the closed-loop stability region ${\Omega }_{\rho }$ such that there exist feasible control actions $u\in U$ for the resulting recurrent neural network that captures the nonlinear dynamics in ${\Omega }_{\rho }$.

3. Recurrent Neural Network

In this work, we develop a nonlinear recurrent neural network (RNN) model with the following form:

$$\dot{\widehat{x}}=F_{nn}(\widehat{x},u):=A\widehat{x}+{\mathsf{\Theta }}^{T}y$$

(4)

where $\widehat{x}\in {\mathbf{R}}^n$ is the RNN state vector and $u\in {\mathbf{R}}^m$ is the manipulated input vector. $y=[y_1,\dots ,y_n,y_{n+1},\dots ,y_{m+n}]=[\sigma ({\widehat{x}}_1),\dots ,\sigma ({\widehat{x}}_n),u_1,\dots ,u_m]\in {\mathbf{R}}^{n+m}$ is a vector of both the network state $\widehat{x}$ and the input u, where $\sigma (·)$ is the nonlinear activation function (e.g., a sigmoid function $\sigma (x)=1/(1+e^{-x})$). A is a diagonal matrix, i.e., $A=diag\{-a_1,\dots ,-a_n\}\in {\mathbf{R}}^{n\times n}$, and $\mathsf{\Theta }$ is a coefficient matrix $\mathsf{\Theta }=[{\theta }_1,\dots ,{\theta }_n]\in {\mathbf{R}}^{(m+n)\times n}$ with ${\theta }_i=b_i[w_{i1},\dots ,w_{i(m+n)}]$, $i=1,\dots ,n$. $a_i$ and $b_i$ are constants, where $a_i$ is assumed to be positive such that each state ${\widehat{x}}_i$ is bounded-input bounded-state stable. $w_{ij}$ is the weight connecting the jth input to the ith neuron where $i=1,\dots ,n$ and $j=1,\dots ,(m+n)$. Since the RNN model of Equation (4) is affine in the manipulated inputs u, the RNN model can also be represented in a similar form of Equation (1) by $\dot{\widehat{x}}=\widehat{f}(\widehat{x})+\widehat{g}(\widehat{x})u$, where $\widehat{f}(·)$ and $\widehat{g}(·)$ are derived from coefficients A and $\mathsf{\Theta }$ in Equation (4). Similarly, $\widehat{f}(·)$ and $\widehat{g}(·)$ are assumed to be sufficiently smooth vector and matrix functions of dimensions $n\times 1$ and $n\times m$, respectively.

It is shown in Figure 1 that RNNs possess a certain type of memory by introducing feedback loops in the hidden layer such that the states information derived from earlier inputs remains in the network. Additionally, if we unfold one of the RNN states in time as shown in Figure 1, the manipulated input vector $u_t$ and the previous state vector $x_{t-1}$ are fed into the current state $x_t$ with the weight matrix $\mathsf{\Theta }$, where $x_{t-1}$, $x_t$, and $x_{t+1}$ are all internal states of an RNN model.

Due to the capability of capturing nonlinear dynamic behavior and the universal approximation theorem ([22,23]) demonstrating that the RNN model with sufficient neurons is able to approximate any dynamic nonlinear system on compact subsets of the state-space for finite time, the RNN model of Equation (4) is used to approximate the continuous-time nonlinear system of Equation (1) in this work.

3.1. RNN Learning Algorithm

The RNN learning algorithm is developed to obtain the optimal weight matrix ${\mathsf{\Theta }}^{*}$, under which the modeling error between the nominal system of Equation (1) (i.e., $w(t)\equiv 0$) and the the RNN model of Equation (4) is minimized. The modeling error $\nu :=F(x,u,0)-F_{nn}(\widehat{x},u)$ under $\mathsf{\Theta }={\mathsf{\Theta }}^{*}$ is assumed to be nonzero since the obtained RNN model may not perfectly approximate the nominal system of Equation (1) due to insufficient number of nodes or layers. Additionally, the modeling error $\nu (t)$ is bounded (i.e., $|\nu (t)|\le {\nu }_m$, ${\nu }_m>0$) since the RNN model of Equation (4) is developed based on the dataset generated within ${\Omega }_{\rho }$ where $x(t)$ and $u(t)\in U$ are bounded. The weight vector ${\theta }_i$ of the RNN model of Equation (4) is also bounded by $|{\theta }_i|\le {\theta }_m$, with ${\theta }_m>0$ to prevent the weights from drifting to infinity. Throughout the manuscript, we will use $x(t)$ to represent the state of the nonlinear system of Equation (1) and $\widehat{x}(t)$ to represent the state of the RNN model of Equation (4).

According to the definition of the modeling error $\nu$, the nominal system of Equation (1) (i.e., $w(t)\equiv 0$) can be written in the following form:

$${\dot{x}}_i=-a_i{x}_i+{\theta }_i^{*T}y+{\nu }_i,i=1,\dots ,n$$

(5)

where the optimal weight vector ${\theta }_i^{*}$ is obtained as follows:

$${\theta }_i^{*}:=arg\underset{|{\theta }_i|\le {\theta }_m}{min}\{{\displaystyle \sum _{k=1}^{N_d}}|F_i(x_k,u_k,0)+a_i{x}_k-{\theta }_i^T{y}_k|\}$$

(6)

where $N_d$ is the number of data samples used for training. Consider the state error $e=\widehat{x}-x\in {\mathbf{R}}^n$. Based on Equations (4) and (5), the time-derivative of the state error is derived as follows:

$${\dot{e}}_i={\dot{\widehat{x}}}_i-{\dot{x}}_i=-a_i{e}_i+{\zeta }_i^{T}y-{\nu }_i,i=1,\dots ,n$$

(7)

where ${\zeta }_i={\theta }_i-{\theta }_i^{*}$ is the error between the current weight vector ${\theta }_i$ and the unknown optimal weight vector ${\theta }_i^{*}$. The weight vector $\theta$ is updated by the following equation:

$${\dot{\theta }}_i=-{\eta }_{i}y{e}_i,i=1,\dots ,n$$

(8)

where the learning rate ${\eta }_i$ is a positive real number. Based on the learning law of Equation (8), the following theorem is established to demonstrate that the state error e remains bounded and is constrained by the modeling error $|\nu |$.

Theorem 1.

Consider the RNN model of Equation (4) with the weights updating law of Equation (8). In the presence of a nonzero modeling error ν, the state error $e_i$ and the weight error ${\zeta }_i$ are bounded, and there exist $\lambda \in \mathbf{R}$ and $\mu >0$ such that the following inequality holds:

$${\int }_0^t{|e(\tau )|}^{2}d\tau \le \lambda +\mu {\int }_0^t{|\nu (\tau )|}^{2}d\tau$$

(9)

Proof. 

We first define a Lyapunov function $\tilde{V}=\frac{1}{2}{\sum }_{i=1}^n(e_i^2+{\zeta }_i^T{\eta }_i^{-1}{\zeta }_i)$ and prove the boundedness of $e_i$ and ${\zeta }_i$ by showing that $\tilde{V}$ is bounded. Specifically, using Equations (7) and (8) and the fact that $\dot{{\zeta }_i}=\dot{{\theta }_i}$, the time-derivative of $\tilde{V}$ is first derived as follows:

$$\begin{array}{cc}\hfill \dot{\tilde{V}}& ={\displaystyle \sum _{i=1}^n}(e_i\dot{e_i}+{\eta }_i^{-1}{\zeta }_i\dot{{\zeta }_i})\hfill \\ \hfill & ={\displaystyle \sum _{i=1}^n}(-a_i{e}_i^2-e_i{\nu }_i)\hfill \end{array}$$

(10)

In the presence of modeling error ${\nu }_i\ne 0$, it is observed that $\dot{\tilde{V}}\le 0$ does not hold for all times. Therefore, let $\alpha :=min\{a_i,1/({\lambda }_m),i=1,\dots ,n\}$ and $\beta :={\sum }_{i=1}^n(2{\theta }_m^2+{\nu }_m^2/2a_i)$, where ${\lambda }_m$ represents the maximum eigenvalue of ${\eta }_i^{-1}$. The following inequality is further derived based on Equation (10):

$$\begin{array}{cc}\hfill \dot{\tilde{V}}& ={\displaystyle \sum _{i=1}^n}(-\frac{a_i}{2}{e}_i^2-\frac{1}{2}|{\zeta }_i{|}^2)+(\frac{1}{2}|{\zeta }_i{|}^2-\frac{a_i}{2}{e}_i^2-e_i{\nu }_i)\hfill \\ \hfill & \le -\alpha \tilde{V}+{\displaystyle \sum _{i=1}^n}(\frac{1}{2}|{\zeta }_i{|}^2+\frac{1}{2a_i}{\nu }_i^2)\hfill \end{array}$$

(11)

Since the weight vector is bounded by $|{\theta }_i|\le {\theta }_m$, it is derived that $\frac{1}{2}{|{\zeta }_i|}^2\le 2{\theta }_m^2$. According to the definitions of $\alpha$ and $\beta$, it follows that $\dot{\tilde{V}}\le -\alpha \tilde{V}+\beta$. Therefore, for all $\tilde{V}\ge V_0=\beta /\alpha$, $\dot{\tilde{V}}\le 0$ holds. This implies that $\tilde{V}$ is bounded, and thus, $e_i$ and ${\zeta }_i$ are bounded as well. Additionally, from Equation (10), it is readily shown that the following equations hold:

$$\begin{array}{cc}\hfill \dot{\tilde{V}}& ={\displaystyle \sum _{i=1}^n}(-\frac{a_i}{2}{e}_i^2-\frac{a_i}{2}{e}_i^2-e_i{\nu }_i)\hfill \\ \hfill & ={\displaystyle \sum _{i=1}^n}(-\frac{a_i}{2}{e}_i^2-(\frac{a_i}{2}{e}_i^2+e_i{\nu }_i+\frac{1}{2a_i}{\nu }_i^2)+\frac{1}{2a_i}{\nu }_i^2)\hfill \\ \hfill & \le {\displaystyle \sum _{i=1}^n}(-\frac{a_i}{2}{e}_i^2+\frac{1}{2a_i}{\nu }_i^2)\hfill \end{array}$$

(12)

From Equation (12), $\tilde{V}(t)$ is obtained as follows:

$$\begin{array}{cc}\hfill \tilde{V}(t)& \le \tilde{V}(0)+{\displaystyle \sum _{i=1}^n}(-\frac{a_i}{2}{\int }_0^t{e}_i{(\tau )}^{2}d\tau +\frac{1}{2a_i}{\int }_0^t{\nu }_i{(\tau )}^{2}d\tau )\hfill \\ \hfill & \le \tilde{V}(0)-\frac{a_{min}}{2}{\int }_0^t{|e(\tau )|}^{2}d\tau +\frac{1}{2a_{min}}{\int }_0^t{|\nu (\tau )|}^{2}d\tau \hfill \end{array}$$

(13)

where $a_{min}$ is the minimum value of $a_i$, $i=1,\dots ,n$. Let $\lambda =\frac{2}{a_{min}}{sup}_{t\ge 0}(\tilde{V}(0)-\tilde{V}(t))$ and $\mu =1/a_{min}^2$. The following inequality is derived:

$$\begin{array}{cc}\hfill {\int }_0^t{|e(\tau )|}^{2}d\tau & \le \frac{2}{a_{min}}(\tilde{V}(0)-\tilde{V}(t))+\frac{1}{a_{min}^2}{\int }_0^t{|\nu (\tau )|}^{2}d\tau \hfill \\ \hfill & \le \lambda +\mu {\int }_0^t{|\nu (\tau )|}^{2}d\tau \hfill \end{array}$$

(14)

The above inequality shows that the state error e is bounded and is proportional to the modeling error $|\nu |$ between the RNN model of Equation (4) and the nominal system of Equation (1) (i.e., $w(t)\equiv 0$). Furthermore, if there exists a constant $C>0$ such that ${\int }_0^{\infty }{|\nu (t)|}^{2}dt$ is bounded (i.e., ${\int }_0^{\infty }{|\nu (t)|}^{2}dt=C<\infty$), then it follows that ${\int }_0^{\infty }{|e(t)|}^{2}dt\le \lambda +\mu C<\infty$. Therefore, $e(t)$ converges to zero asymptotically since $e(t)$ is a uniformly continuous function implied by the boundedness of its derivative $\dot{e}$. □

3.2. Ensemble of RNN Models

In this section, an RNN model is initially developed to approximate the nonlinear system of Equation (1) within an operating region. Based on that, the ensemble learning using the aforementioned RNN learning algorithm is used with a k-fold cross validation to improve prediction accuracy of RNN models.

3.2.1. Development of RNN Models

The development of RNN models consists of two steps: data generation and training process. In general, the dataset for training an RNN model is a collection of sampled data from industrial processes, laboratory experiments, or simulations of a given system. There are many high-quality training datasets, for example, the collection of databases from the UCI Machine Learning Repository [24] that is widely used by the researchers in machine learning community. In this work, we use the simulation data of the nominal system of Equation (1) (i.e., $w(t)\equiv 0$) as the training and validation datasets for the RNN model of Equation (4). Specifically, to generate the dataset that captures the system dynamics for $x\in {\Omega }_{\rho }$ and $u\in U$, open-loop simulations are first conducted for initial conditions $x_0\in {\Omega }_{\rho }$ with the constrained input vector $u\in U$. It should be mentioned that the continuous system of Equation (1) is simulated in a sample-and-hold fashion (i.e., the input is a piecewise constant function, i.e., $u(t)=u(t_k)$, $\forall t\in [t_k,t_{k+1})$, where $t_{k+1}:=t_k+\Delta$ and $\Delta$ is the sampling period). Then, the nominal system of Equation (1) is integrated via explicit Euler method with a sufficiently small integration time step $h_c<\Delta$. The dataset is generated with the following steps: (1) The targeted region ${\Omega }_{\rho }$ and the range of inputs $u\in U$ are first discretized with sufficiently small intervals. (2) Open-loop simulations of the nominal system of Equation (1) are performed in a sample-and-hold fashion with various initial states $x_0\in {\Omega }_{\rho }$ and inputs $u\in U$ to obtain a large number of state trajectories for a fixed length of time to sweep over all the values that $(x,u)$ can take. It is noted that data points are collected at every integration time step $h_c$ instead of every sampling step $\Delta$ to obtain intermediate states as many as possible in order to approximate a continuous system of Equation (1). (3) Then, the obtained state trajectories are separated into many time-series data samples with the length of $P_{nn}$, which represents the prediction period for RNN models. In other words, the RNN model is developed to predict future states over $t\in [t_k,t_k+P_{nn}]$, where all the data points between $t_k$ and $t_k+P_{nn}$ are treated as the internal states of the RNN model with the time interval of $h_c$ (i.e., the time interval between two consecutive states $x_{t-1}$ and $x_t$ for the unfolded RNN shown in Figure 1 is equal to $h_c$). (4) The dataset is partitioned into training, validation and testing datasets.

Next, the RNN model is developed based on the dataset generated from open-loop simulations via an application program interface (API), i.e., Keras ([13]). To guarantee that the RNN model of Equation (4) achieves good performance in a neighborhood around the origin and has the same equilibrium point as the nonlinear system of Equation (1), the modeling error is further required to be bounded by $|\nu |\le \gamma |x|\le {\nu }_m$ during the training process. Before starting an RNN training process, the number of layers and nodes used, the optimizer, and the type of loss function need to be specified at the beginning. Specifically, in our work, the number of layers and nodes for RNNs are determined through a grid search and the adaptive moment estimation method (i.e., Adam in Keras) is used as RNN optimizer. The loss function is the mean absolute percentage error between predicted states $\widehat{x}$ and actual states x from training data. Additionally, to prevent RNNs from over-training, the early stopping criteria is employed to terminate the training process before the maximum number of epochs is reached if the following two conditions are satisfied: (1) The modeling error is below its threshold. (2) The performance on the validation dataset stops improving.

Remark 2.

It is noted that original data samples from open-loop simulations are shortened to a fixed length, $P_{nn}$, in the step 3 of data generation above. Through this data preprocessing step, the dataset is applicable for developing an RNN with the prediction period $P_{nn}$, which is chosen to be an integer multiple of the sampling period Δ such that the RNN model of Equation (4) is able to predict future states for at least one sampling period. One of the advantages of introducing step 2 and step 3 is that a new dataset can be conveniently generated for an RNN with a different prediction period $P_{nn}$ while the original dataset in step 2 only needs to be created once.

Remark 3.

Unlike feedforward neural networks (FNN) that take state measurement $x(t_k)$ and manipulated inputs $u(t_k)$ at time $t=t_k$ to predict the state at $t=t_k+P_{nn}$ only, RNNs use all the intermediate states between $x(t_k)$ and $x(t_k+P_{nn})$ to predict the state trajectory over this prediction period. In other words, all the internal states with the time step $h_c$ are used as one of the inputs (the other one is the manipulated input u) to the RNN (shown in Figure 1) over $t\in [t_k,t_k+P_{nn}]$ during the training process. Therefore, besides the final predicted state $x(t_k+P_{nn})$, all the internal states before can be predicted and exported (if needed) by RNNs. Additionally, it is important to mention that since the time interval for RNN internal states is chosen to be the sufficiently small integration time step $h_c$ that is used to simulate the continuous-time nonlinear system of Equation (1), the obtained RNN model can be regarded as a good representation for the continuous-time RNN model of Equation (4).

Remark 4.

Data quality plays an important role in data-driven modeling techniques including RNN modeling. In case a good RNN model is initially not available due to poor data quality, the RNN model can be improved via on-line update through real-time process state measurements. Additionally, if the real-time process measurements are of poor quality, for example, noisy measurements, a data preprocessing step can be employed before training an RNN. Overall, creating a good dataset remains an important issue in RNN modeling; in the chemical reactor example in Section 5 the data is generated from open-loop simulations.

Remark 5.

The RNN models considered in our work involve the use of nonlinear terms in the dynamic model (this can be thought of as a nonlinear process gain $k(x)$, multiplying the input), which allows capturing gain sign change in the different regions in the operating state-space if it occurs. This issue does not occur in the operating region of the chemical reactor example considered in Section 5. The reader may refer to [25] for an approach to neural network modeling that captures gain change.

3.2.2. Ensemble Learning

Ensemble learning is a machine learning process that combines multiple learning algorithms to obtain better prediction performance ([26]) in terms of reduced variability and improved generalization performance. The reasons that ensemble regression models are able to improve prediction performance are given in [10,26]. In this work, homogeneous ensemble regression models are constructed following the RNN learning algorithm in Section 3.1 and a k-fold cross validation. Specifically, as shown in Figure 2, the dataset is initially separated into k folds with the same size, where $k-1$ folds are used for training and the remaining one is used for validation. By choosing a different subset as the validation dataset every time, a total number of k different RNNs can be obtained. Subsequently, the final prediction results of an ensemble of RNN models is calculated using the stacking method that averages the predictions of multiple RNN models. It is noted that besides the stacking method, other ensemble-based algorithms such as Boosting and Bagging [19] can also be applied to improve the overall prediction performance within the targeted operating region.

3.2.3. Parallel Computing

Parallel computing is a type of computation in which the execution of multiple processes or of the computational tasks are carried out simultaneously. Since the final prediction results in ensemble learning are obtained by combining the results of individual RNN models that are independent from each other, the computational tasks for all RNN models can be operated in parallel to improve computational efficiency. The algorithm for parallel computation of an ensemble of k RNN models is shown in Figure 3. The blocks on the left demonstrate the state prediction using parallel computing and the dashed blocks on the right represent the actual nonlinear system of Equation (1) and the model-based controller that will be discussed in the next section. Specifically, we first reserve $k+1$ nodes, where node 0 is the host node and the rest are worker nodes. The host node is used to exchange information while the worker nodes are mainly used for computation. Each worker node is assigned with a corresponding RNN model for prediction. At time $t=t_k$, the model-based controller receives the state measurement $x(t_k)$ and sends it to the host node 0 along with the manipulated inputs $u^g(t_k)$ (i.e., the guess of control actions over the prediction horizon). The host node broadcasts $x(t_k)$ and $u^g(t_k)$ to all worker nodes such that all RNN models share the same input vector. Prediction of RNN models are carried out in all worker nodes simultaneously. It is noted that a synchronization barrier is placed before the ensemble combiner in node 0 to ensure that all the computation tasks in worker nodes are completed. Lastly, the host node combines prediction results from multiple RNN models to obtain predicted states over $t\in [t_k,t_k+Pnn]$ and sends them to the model-based controller. The optimal control actions $u^{*}(t_k)$ is calculated by the optimizer solver in model-based controller by taking various guesses of control actions $u^g(t_k)$, and will be sent to the nonlinear process of Equation (1) to be applied for the next sampling period. The above process is repeated as the horizon of the nonlinear process of Equation (1) is rolled one sampling period forward each time.

4. Lyapunov-Based EMPC Using Ensemble RNN Models

In this section, the Lyapunov-based economic MPC (LEMPC) that incorporates an ensemble of RNN models to predict future states is developed for the nonlinear system of Equation (1). In the following subsections, a Lyapunov-based controller using the RNN model of Equation (4) is first used to characterize the closed-loop stability region in which the origin of the nonlinear system of Equation (1) can be rendered exponentially stable. Then, the formulation of the LEMPC using an ensemble of RNN models is given and Theorem 2 is established to demonstrate closed-loop stability for all initial conditions in the closed-loop stability region.

4.1. Lyapunov-Based Control Using RNN Models

We assume that there exists a stabilizing feedback controller $u={\mathsf{\Phi }}_{nn}(x)\in U$ that can render the origin of the RNN model of Equation (4) exponentially stable in an open neighborhood $\widehat{D}$ in the sense that there exists a ${\mathcal{C}}^1$ Lyapunov function $\widehat{V}(x)$ such that the following inequalities hold for all x in $\widehat{D}$:

$${\widehat{c}}_1{|x|}^2\le \widehat{V}(x)\le {\widehat{c}}_2{|x|}^2,$$

(15a)

$$\frac{\partial \widehat{V}(x)}{\partial x}{F}_{nn}(x,{\mathsf{\Phi }}_{nn}(x))\le -{\widehat{c}}_3{|x|}^2,$$

(15b)

$$\left|\frac{\partial \widehat{V}(x)}{\partial x}\right|\le {\widehat{c}}_4|x|$$

(15c)

where ${\widehat{c}}_1$, ${\widehat{c}}_2$, ${\widehat{c}}_3$ and ${\widehat{c}}_4$ are positive constants. Although the RNN model of Equation (4) is obtained based on the dataset within the closed-loop stability region ${\Omega }_{\rho }$ for the nonlinear system of Equation (1), it is noted that the existence a feasible control action $u\in U$ that can render the origin of the RNN model of Equation (4) exponentially stable is not guaranteed due to model mismatch. Therefore, based on Equation (15), the closed-loop stability region is re-characterized for the RNN model of Equation (4) accounting for the input constraints. Specifically, a Lyapunov-based controller based on the RNN model is first developed using the universal Sontag control law ([21]) since the RNN model of Equation (4) can be represented by $\dot{\widehat{x}}=\widehat{f}(\widehat{x})+\widehat{g}(\widehat{x})u$. Then, following the characterization method of ${\Omega }_{\rho }$ for the nonlinear system of Equation (1), an open set from which the origin of the RNN model of Equation (4) can be rendered exponentially stable under the controller $u={\mathsf{\Phi }}_{nn}(x)\in U$ is characterized as follows: ${\widehat{\varphi }}_u:=\{x\in {\mathbf{R}}^n|\dot{\widehat{V}}(x)<-k\widehat{V}(x),u={\mathsf{\Phi }}_{nn}(x)\in U\}\cup \{0\}$, where $\dot{\widehat{V}}(x)=\frac{\partial \widehat{V}(x)}{\partial x}{F}_{nn}(x,{\mathsf{\Phi }}_{nn}(x))$ and $k>0$. The closed-loop stability region for the RNN model of Equation (4) is defined as the largest level of $\widehat{V}$ inside ${\widehat{\varphi }}_u$: ${\Omega }_{\widehat{\rho }}:=\{x\in {\widehat{\varphi }}_u|\widehat{V}(x)\le \widehat{\rho }\}$, where $\widehat{\rho }>0$. From now on, ${\Omega }_{\widehat{\rho }}$ will be taken as the new operating region for the closed-loop system of Equation (1) in the following sections. Additionally, the Lipschitz property of $F_{nn}(x,u)$ and the boundedness of u imply that there exist positive constants $M_{nn}$ and $L_{nn}$ such that the following inequalities hold for all $x,x^{\prime }\in {\Omega }_{\widehat{\rho }},u\in U$:

$$|F_{nn}(x,u)|\le M_{nn}$$

(16a)

$$\begin{array}{c}\left|\frac{\partial \widehat{V}(x)}{\partial x}{F}_{nn}(x,u)-\frac{\partial \widehat{V}(x^{\prime })}{\partial x}{F}_{nn}(x^{\prime },u)\right|\hfill \\ \le L_{nn}|x-x^{\prime }|\hfill \end{array}$$

(16b)

The following proposition is developed to demonstrate that the state of the nonlinear system of Equation (1) can be driven towards the origin exponentially for all $x_0\in {\Omega }_{\widehat{\rho }}$ under the stabilizing controller $u={\mathsf{\Phi }}_{nn}(x)\in U$ designed for the RNN model of Equation (4) provided that the modeling error is sufficiently.

Proposition 1.

Consider the stabilizing controller $u={\mathsf{\Phi }}_{nn}(x)\in U$ that renders the origin of the RNN system of Equation (4) exponentially stable for all $x_0\in {\Omega }_{\widehat{\rho }}$. If there exists a positive real number γ that constrains the modeling error $|\nu |=|F(x,u,0)-F_{nn}(x,u)|\le \gamma |x|\le {\nu }_m$, $\forall x\in {\Omega }_{\widehat{\rho }}$, $u\in U$ and satisfies the following inequality:

$$\gamma <{\widehat{c}}_3/{\widehat{c}}_4$$

(17)

then the origin of the nominal closed-loop system of Equation (1) under $u={\mathsf{\Phi }}_{nn}(x)\in U$ is also exponentially stable for all $x\in {\Omega }_{\widehat{\rho }}$.

Proof. 

We demonstrate the exponential stability of the origin of the nominal system of Equation (1) under $u={\mathsf{\Phi }}_{nn}(x)\in U$, $\forall x\in {\Omega }_{\widehat{\rho }}$ by showing that there exists a positive real number $k^{\prime }$ such that $\dot{\widehat{V}}(x)\le -k^{\prime }{|x|}^2$ holds for all $\forall x\in {\Omega }_{\widehat{\rho }}$, where $\dot{\widehat{V}}(x)=\frac{\partial \widehat{V}(x)}{\partial x}F(x,{\mathsf{\Phi }}_{nn}(x),0)$. Based on Equations (15b) and (15c), the time-derivative of $\widehat{V}$ is derived as follows:

$$\begin{array}{cc}\hfill \dot{\widehat{V}}=& \frac{\partial \widehat{V}(x)}{\partial x}F(x,{\mathsf{\Phi }}_{nn}(x),0)\hfill \\ \hfill =& \frac{\partial \widehat{V}(x)}{\partial x}(F_{nn}(x,{\mathsf{\Phi }}_{nn}(x))+F(x,{\mathsf{\Phi }}_{nn}(x),0)\hfill \\ \hfill & -F_{nn}(x,{\mathsf{\Phi }}_{nn}(x)))\hfill \\ \hfill \le & -{\widehat{c}}_3{|x|}^2+{\widehat{c}}_4|x|(F(x,{\mathsf{\Phi }}_{nn}(x),0)-F_{nn}(x,{\mathsf{\Phi }}_{nn}(x)))\hfill \\ \hfill \le & -{\widehat{c}}_3{|x|}^2+{\widehat{c}}_4\gamma {|x|}^2\hfill \end{array}$$

(18)

If $\gamma$ satisfies Equation (17), then there exists a positive real number ${\tilde{c}}_3=-{\widehat{c}}_3+{\widehat{c}}_4\gamma >0$ such that $\dot{\widehat{V}}\le -{\tilde{c}}_3{|x|}^2\le 0$ holds. Let $k^{\prime }={\tilde{c}}_3$ and this completes the proof that the closed-loop state of the nominal system of Equation (1) is driven to the origin in an exponential rate of convergence under $u={\mathsf{\Phi }}_{nn}(x)\in U$ for all $x_0\in {\Omega }_{\widehat{\rho }}$. □

4.2. LEMPC Using an Ensemble of RNN Models

The Lyapunov-based economic MPC (LEMPC) design using an ensemble of RNN models is represented by the following optimization problem:

$$\mathcal{J}=\underset{u\in S(\Delta )}{\max }{\int }_{t_k}^{t_{k+N}}{l}_e(\tilde{x}(t),u(t))dt$$

(19a)

$$\mathrm{s}.\mathrm{t}.\dot{\tilde{x}}(t)=\frac{1}{N_e}{\displaystyle \sum _{j=1}^{N_e}}{F}_{nn}^j(\tilde{x}(t),u(t))$$

(19b)

$$u(t)\in U,\forall t\in [t_k,t_{k+N})$$

(19c)

$$\tilde{x}(t_k)=x(t_k)$$

(19d)

$$\begin{array}{c}\widehat{V}(\tilde{x}(t))\le {\widehat{\rho }}_e,\forall t\in [t_k,t_{k+N}),\hfill \\ \mathrm{if}x(t_k)\in {\Omega }_{{\widehat{\rho }}_e}\hfill \end{array}$$

(19e)

$$\begin{array}{c}\dot{\widehat{V}}(x(t_k),u)\le \dot{\widehat{V}}(x(t_k),{\mathsf{\Phi }}_{nn}(x(t_k)),\hfill \\ \mathrm{if}x(t_k)\in {\Omega }_{\widehat{\rho }}\backslash {\Omega }_{{\widehat{\rho }}_e}\hfill \end{array}$$

(19f)

where $S(\Delta )$ is the set of piecewise constant functions with period $\Delta$, $\tilde{x}$ is the predicted state trajectory, $N_e$ is the number of RNN models used in predicting future states, and N is the number of sampling periods in the prediction horizon. We use $\dot{\widehat{V}}(x,u)$ to represent the time-derivative of $\widehat{V}$, i.e., $\dot{\widehat{V}}(x,u)=\frac{\partial \widehat{V}(x)}{\partial x}(F_{nn}(x,u))$. The optimization problem of Equation (19) is solved by optimizing the time integral of the stage cost function $l_e(\tilde{x}(t),u(t))$ of Equation (19a) over the prediction horizon accounting for the constraints of Equations (19b)–(19f). The constraint of Equation (19b) is the ensemble of RNN models of Equation (4) that is used as the prediction model. Equation (19c) defines the input constraints that are applied over the entire prediction horizon. Equation (19d) defines the initial condition $\tilde{x}(t_k)$ of Equation (19b). The constraint of Equation (19e) maintains the closed-loop state predicted by Equation (19b) in ${\Omega }_{{\widehat{\rho }}_e}$ over the prediction horizon if the state $x(t_k)$ is inside ${\Omega }_{{\widehat{\rho }}_e}$. However, if $x(t_k)$ leaves ${\Omega }_{{\widehat{\rho }}_e}$ but still remains in ${\Omega }_{\widehat{\rho }}$, the contractive constraint of Equation (19f) drives the state towards the origin for the next sampling period such that the state will eventually enter ${\Omega }_{{\widehat{\rho }}_e}$ within finite sampling periods. By introducing a conservative region ${\Omega }_{{\widehat{\rho }}_e}$ with ${\widehat{\rho }}_e<\widehat{\rho }$ and the stability constraints of Equations (19e) and (19f), it will be proved in Theorem 2 that for any initial condition $x_0\in {\Omega }_{\widehat{\rho }}$, the closed-loop state of the nonlinear system of Equation (1) is bounded in ${\Omega }_{\widehat{\rho }}$ for all times. A schematic of ${\Omega }_{\widehat{\rho }}$ and ${\Omega }_{{\widehat{\rho }}_e}$ is given by Figure 4, in which it is shown that the state trajectory starting inside ${\Omega }_{\widehat{\rho }}$ is bounded in ${\Omega }_{\widehat{\rho }}$ under LEMPC.

The optimal input trajectory is calculated over the entire prediction horizon $t\in [t_k,t_{k+N})$ and is denoted by $u^{*}(t)$. However, only the control action $u^{*}(t_k)$ for the first sampling period in $t\in [t_k,t_{k+N})$ is sent to the nonlinear system of Equation (1) to be applied over the next sampling period. Before we demonstrate closed-loop stability for the nonlinear system of Equation (1) under the LEMPC of Equation (19), the next two propositions are first derived to provide useful tools for proving closed-loop stability in Theorem 2. Specifically, the following proposition provides an upper bound for the state error in the presence of sufficiently small modeling error and bounded disturbances.

Proposition 2.

Consider the solution $x(t)$ of the nonlinear system $\dot{x}=F(x,u,w)$ of Equation (1) in the presence of bounded disturbances $|w(t)|\le w_m$, and the solution $\widehat{x}(t)$ of the RNN model $\dot{\widehat{x}}=F_{nn}(\widehat{x},u)$ of Equation (4) with the same initial condition $x_0={\widehat{x}}_0\in {\Omega }_{\widehat{\rho }}$. If $x(t)$, $\widehat{x}(t)\in {\Omega }_{\widehat{\rho }}$, and the modeling error is bounded (i.e., $|\nu (t)|=|\dot{x}-\dot{\widehat{x}}|\le {\nu }_m$) for all times, then there exists a positive constant κ and a class $\mathcal{K}$ function $f_w(·)$ such that the following inequalities hold $\forall x,\widehat{x}\in {\Omega }_{\widehat{\rho }}$ and $w(t)\in W$:

$$|x(t)-\widehat{x}(t)|\le f_w(t):=\frac{L_w{w}_m+{\nu }_m}{L_x}(e^{L_{x}t}-1)$$

(20a)

$$\widehat{V}(x)\le \widehat{V}(\widehat{x})+\frac{{\widehat{c}}_4\sqrt{\widehat{\rho }}}{\sqrt{{\widehat{c}}_1}}|x-\widehat{x}|+\kappa |x-\widehat{x}{|}^2$$

(20b)

Proof. 

Following the proof in [5,7], we define the state error vector by $e(t)=x(t)-\widehat{x}(t)$. The time-derivative of $e(t)$ is obtained $\forall x,\widehat{x}\in {\Omega }_{\widehat{\rho }}$, $u\in U$ and $w(t)\in W$ as follows:

$$\begin{array}{cc}\hfill |\dot{e}|& =|F(x,u,w)-F_{nn}(\widehat{x},u)|\hfill \\ \hfill & \le |F(x,u,w)-F(\widehat{x},u,0)|+|F(\widehat{x},u,0)-F_{nn}(\widehat{x},u)|\hfill \\ \hfill & \le L_x|e(t)|+L_w{w}_m+{\nu }_m\hfill \end{array}$$

(21)

where the last inequality is derived from Equation (3b) and the fact that $|\nu |\le {\nu }_m$. Since $x(t)$ and $\widehat{x}(t)$ share the same initial condition, (i.e., $e(0)=0$), the upper bound for $|e(t)|$ is derived for all $x(t),\widehat{x}(t)\in {\Omega }_{\widehat{\rho }}$ and $|w(t)|\le w_m$ as follows:

$$|e(t)|=|x(t)-\widehat{x}(t)|\le \frac{L_w{w}_m+{\nu }_m}{L_x}(e^{L_{x}t}-1)$$

(22)

Subsequently, $\forall x,\widehat{x}\in {\Omega }_{\rho }$, Equation (20b) is derived based on the Taylor series expansion of $\widehat{V}(x)$ around $\widehat{x}$ as follows:

$$\begin{array}{cc}\hfill \widehat{V}(x)& \le \widehat{V}(\widehat{x})+\frac{\partial \widehat{V}(\widehat{x})}{\partial x}|x-\widehat{x}|+\kappa |x-\widehat{x}{|}^2\hfill \\ \hfill & \le \widehat{V}(\widehat{x})+\frac{{\widehat{c}}_4\sqrt{\widehat{\rho }}}{\sqrt{{\widehat{c}}_1}}|x-\widehat{x}|+\kappa |x-\widehat{x}{|}^2\hfill \end{array}$$

(23)

where $\kappa$ is a positive real number, and the last inequality is derived using Equations (15a) and (15c). □

The above proposition demonstrates that the error of state trajectories of the nonlinear system of Equation (1) and the RNN model of Equation (4), starting from the same initial condition, is bounded for finite time. The next proposition is developed to demonstrate that if $x\in {\Omega }_{\widehat{\rho }}\backslash {\Omega }_{{\rho }_s}$, the stabilizing controller $u={\mathsf{\Phi }}_{nn}(x)\in U$ implemented in a sample-and-hold fashion is able to maintain $\dot{\widehat{V}}(x)$ negative such that the state will be driven towards the origin and ultimately enters a small neighborhood around the origin (i.e., ${\Omega }_{{\rho }_s}$) in finite sampling steps.

Proposition 3.

Consider the system of Equation (1) under the controller $u={\mathsf{\Phi }}_{nn}(\widehat{x})\in U$ that meets the conditions of Equation (15) and is implemented in a sample-and-hold fashion, i.e., $u(t)={\mathsf{\Phi }}_{nn}(\widehat{x}(t_k))$, $\forall t\in [t_k,t_{k+1})$, where $t_{k+1}:=t_k+\Delta$. Let ${ϵ}_w,{ϵ}_s>0$, $\Delta >0$ and $\widehat{\rho }>{\widehat{\rho }}_e>{\rho }_s>0$ satisfy

$$-\frac{{\widehat{c}}_3}{{\widehat{c}}_2}{\rho }_s+L_{nn}{M}_{nn}\Delta \le -{ϵ}_s$$

(24a)

$$-\frac{{\tilde{c}}_3}{{\widehat{c}}_2}{\rho }_s+L_x^{{}^{\prime }}M\Delta +L_w^{{}^{\prime }}{w}_m\le -{ϵ}_w$$

(24b)

$${\widehat{\rho }}_e>max\{\widehat{V}(\widehat{x}(t_k+\Delta ))|\widehat{x}(t_k)\in {\Omega }_{{\rho }_s},u\in U,w\in W\}$$

(24c)

Then, for any $x(t_k)\in {\Omega }_{\widehat{\rho }}\backslash {\Omega }_{{\rho }_s}$, the following inequalities holds:

$$\widehat{V}(\widehat{x}(t))\le \widehat{V}(\widehat{x}(t_k)),\forall t\in [t_k,t_{k+1})$$

(25a)

$$\widehat{V}(x(t))\le \widehat{V}(x(t_k)),\forall t\in [t_k,t_{k+1})$$

(25b)

Proof. 

To show that the value of $\widehat{V}$ decreases along the trajectory $\widehat{x}(t)$ of the RNN model of Equation (4) over $t\in [t_k,t_{k+1})$, we calculate the time-derivative of $\widehat{V}(\widehat{x})$ based on $\widehat{x}(t)$ as follows:

$$\begin{array}{cc}\hfill \dot{\widehat{V}}(\widehat{x}(t))=& \frac{\partial \widehat{V}(\widehat{x}(t))}{\partial \widehat{x}}{F}_{nn}(\widehat{x}(t),{\mathsf{\Phi }}_{nn}(\widehat{x}(t_k)))\hfill \\ \hfill =& \frac{\partial \widehat{V}(\widehat{x}(t_k))}{\partial \widehat{x}}{F}_{nn}(\widehat{x}(t_k),{\mathsf{\Phi }}_{nn}(\widehat{x}(t_k)))\hfill \\ \hfill & +\frac{\partial \widehat{V}(\widehat{x}(t))}{\partial \widehat{x}}{F}_{nn}(\widehat{x}(t),{\mathsf{\Phi }}_{nn}(\widehat{x}(t_k)))\hfill \\ \hfill & -\frac{\partial \widehat{V}(\widehat{x}(t_k))}{\partial \widehat{x}}{F}_{nn}(\widehat{x}(t_k),{\mathsf{\Phi }}_{nn}(\widehat{x}(t_k)))\hfill \end{array}$$

(26)

Using Equations (15a) and (15b) and the Lipschitz condition of Equation (16), the following inequalities are obtained:

$$\begin{array}{cc}\hfill \dot{\widehat{V}}(\widehat{x}(t))\le & -\frac{{\widehat{c}}_3}{{\widehat{c}}_2}{\rho }_s+\frac{\partial \widehat{V}(\widehat{x}(t))}{\partial \widehat{x}}{F}_{nn}(\widehat{x}(t),{\mathsf{\Phi }}_{nn}(\widehat{x}(t_k)))\hfill \\ \hfill & -\frac{\partial \widehat{V}(\widehat{x}(t_k))}{\partial \widehat{x}}{F}_{nn}(\widehat{x}(t_k),{\mathsf{\Phi }}_{nn}(\widehat{x}(t_k)))\hfill \\ \hfill \le & -\frac{{\widehat{c}}_3}{{\widehat{c}}_2}{\rho }_s+L_{nn}{M}_{nn}\Delta \hfill \end{array}$$

(27)

Therefore, $\dot{\widehat{V}}(\widehat{x}(t))\le -{ϵ}_s$ holds $\forall \widehat{x}(t_k)\in {\Omega }_{\widehat{\rho }}\backslash {\Omega }_{{\rho }_s}$ and $t\in [t_k,t_{k+1})$ if Equation (24a) is satisfied. By integrating the above inequality, it follows that $\widehat{V}(\widehat{x}(t))\le \widehat{V}(\widehat{x}(t_k))-\Delta {ϵ}_s$, $\forall \widehat{x}(t_k)\in {\Omega }_{\widehat{\rho }}\backslash {\Omega }_{{\rho }_s}$, $t\in [t_k,t_{k+1})$ (i.e., Equation (25a).

Next, to show that $\widehat{V}(x(t))\le \widehat{V}(x(t_k))$ holds for all $t\in [t_k,t_{k+1})$, the time-derivative of $\widehat{V}(x)$ for the nonlinear system of Equation (1) (i.e., $\dot{x}=F(x,u,w)$) in the presence of bounded disturbances (i.e., $|w(t)|\le w_m$) is derived as follows:

$$\begin{array}{cc}\hfill \dot{\widehat{V}}(x(t))& =\frac{\partial \widehat{V}(x(t))}{\partial x}F(x(t),{\mathsf{\Phi }}_{nn}(x(t_k)),w)\hfill \\ \hfill & =\frac{\partial \widehat{V}(x(t_k))}{\partial x}F(x(t_k),{\mathsf{\Phi }}_{nn}(x(t_k)),0)\hfill \\ \hfill & +\frac{\partial \widehat{V}(x(t))}{\partial x}F(x(t),{\mathsf{\Phi }}_{nn}(x(t_k)),w)\hfill \\ \hfill & -\frac{\partial \widehat{V}(x(t_k))}{\partial x}F(x(t_k),{\mathsf{\Phi }}_{nn}(x(t_k)),0)\hfill \end{array}$$

(28)

Since Equation (18) shows that $\frac{\partial \widehat{V}(x(t_k))}{\partial x}F(x(t_k),{\mathsf{\Phi }}_{nn}(x(t_k)),0)\le -{\tilde{c}}_3{|x(t_k)|}^2$ holds for all $x\in {\Omega }_{\widehat{\rho }}$, the following inequality is obtained for all $x(t_k)\in {\Omega }_{\widehat{\rho }}\backslash {\Omega }_{{\rho }_s}$, $t\in [t_k,t_{k+1})$ using Equation (2a) and the Lipschitz condition in Equation (3):

$$\begin{array}{cc}\hfill \dot{\widehat{V}}(x(t))\le & -\frac{{\tilde{c}}_3}{{\widehat{c}}_2}{\rho }_s+L_x^{{}^{\prime }}|x(t)-x(t_k)|+L_w^{{}^{\prime }}|w|\hfill \\ \hfill \le & -\frac{{\tilde{c}}_3}{{\widehat{c}}_2}{\rho }_s+L_x^{{}^{\prime }}M\Delta +L_w^{{}^{\prime }}{w}_m\hfill \end{array}$$

(29)

Therefore, the following inequality is further derived for all $x(t_k)\in {\Omega }_{\widehat{\rho }}\backslash {\Omega }_{{\rho }_s}$ and $t\in [t_k,t_{k+1})$ if the condition of Equation (24b) is satisfied:

$$\dot{\widehat{V}}(x(t))\le -{ϵ}_w$$

(30)

Similarly, this implies that $\widehat{V}(x(t))\le \widehat{V}(x(t_k))-\Delta {ϵ}_w$, $\forall x(t_k)\in {\Omega }_{\widehat{\rho }}\backslash {\Omega }_{{\rho }_s}$, $t\in [t_k,t_{k+1})$. Therefore, the state of the nonlinear system of Equation (1) will enter ${\Omega }_{{\rho }_s}$ in finite sampling periods if $x(t_k)\in {\Omega }_{\widehat{\rho }}\backslash {\Omega }_{{\rho }_s}$. Additionally, if $x(t_k)\in {\Omega }_{{\rho }_s}$, in which Equations (29) and (30) do not hold, Equation (24c) guarantees that the state will not leave ${\Omega }_{{\widehat{\rho }}_e}$ in one sampling period for all $u\in U$ and $w\in W$. If the state $x(t_{k+1})$ leaves ${\Omega }_{{\rho }_s}$ but remains in ${\Omega }_{{\widehat{\rho }}_e}$, at the next sampling period $t\in [t_{k+1},t_{k+2})$, Equation (30) is satisfied again and the state will be driven towards the origin. Therefore, the state of the nonlinear system of Equation (1) is bounded in ${\Omega }_{\widehat{\rho }}$ for all times. □

Based on Propositions 2 and 3, the following theorem is established to show recursive feasibility of the LEMPC optimization problem of Equation (19), and the boundedness of the closed-loop state within ${\Omega }_{\widehat{\rho }}$ under the sample-and-hold implementation of the LEMPC.

Theorem 2.

Consider the closed-loop system of Equation (1) under the sample-and-hold implementation of the LEMPC of Equation (19) with the stabilizing controller ${\mathsf{\Phi }}_{nn}(x)$ that satisfies Equation (15). Let $\Delta >0$, ${ϵ}_w>0$ and $\widehat{\rho }>{\widehat{\rho }}_e>0$ satisfy Equation (24) and the following inequality:

$${\widehat{\rho }}_e\le \widehat{\rho }-\frac{{\widehat{c}}_4\sqrt{\widehat{\rho }}}{\sqrt{{\widehat{c}}_1}}{f}_w(\Delta )-\kappa {(f_w(\Delta ))}^2$$

(31)

If $x_0\in {\Omega }_{\widehat{\rho }}$ and the conditions of Propositions 2 and 3 are satisfied, then there always exists a feasible solution for the optimization problem of Equation (19), and the closed-loop state $x(t)$ is bounded in the closed-loop stability region ${\Omega }_{\widehat{\rho }}$, $\forall t\ge 0$.

Proof. 

The LEMPC of Equation (19) predicts future states by averaging the results of an ensemble of RNN models in Equation (19b). Since each RNN model is trained to satisfy the modeling error constraint and Equation (16), it is readily shown that the averaged results of an ensemble of multiple RNN models also satisfy the above conditions, and thus, the results derived in Propositions 1–3 for a single RNN model can be generalized to an RNN ensemble.

We first prove that the optimization problem of Equation (19) can be solved with recursive feasibility for all $x\in {\Omega }_{\widehat{\rho }}$. Specifically, if $x(t_k)\in {\Omega }_{{\widehat{\rho }}_e}$, the control actions ${\mathsf{\Phi }}_{nn}(x(t_{k+i}))$, $i=0,1,\dots ,N-1$ satisfy the input constraint of Equation (19c) and the Lyapunov-based constraint of Equation (19e) since it is shown in Equation (25a) that the states predicted by the RNN model of Equation (19b) remain inside ${\Omega }_{{\widehat{\rho }}_e}$ under the controller ${\mathsf{\Phi }}_{nn}(x)$. Additionally, if $x(t_k)\in {\Omega }_{\widehat{\rho }}\backslash {\Omega }_{{\widehat{\rho }}_e}$, the control action $u(t)={\mathsf{\Phi }}_{nn}(x(t_k))\in U$, $t\in [t_k,t_{k+1})$ satisfies the input constraint of Equation (19c) and the Lyapunov-based constraint of Equation (19f) such that the state can be driven towards the origin during the next sampling period. Therefore, the stabilizing controller $u={\mathsf{\Phi }}_{nn}(x)\in U$ provides a feasible solution that satisfies all the constraints of the LEMPC optimization problem of Equation (19) if $x(t)\in {\Omega }_{\widehat{\rho }}$ for all times.

Next, we prove that for $x_0\in {\Omega }_{\widehat{\rho }}$, the state of the closed-loop system of Equation (1) is bounded in ${\Omega }_{\widehat{\rho }}$ for all times. Specifically, if $x(t_k)\in {\Omega }_{{\widehat{\rho }}_e}$, the predicted states $\widehat{x}(t)$ of the RNN model of Equation (19b) are maintained in ${\Omega }_{{\widehat{\rho }}_e}$ under the constraint of Equation (19e). According to Proposition 2, the actual state $x(t)$, $t\in [t_k,t_{k+1})$ of the nonlinear system of Equation (1) is bounded by the following inequality:

$$\begin{array}{cc}\hfill \widehat{V}(x)& \le \widehat{V}(\widehat{x})+\frac{{\widehat{c}}_4\sqrt{\widehat{\rho }}}{\sqrt{{\widehat{c}}_1}}|x-\widehat{x}|+\kappa |x-\widehat{x}{|}^2\hfill \\ \hfill & \le \widehat{V}(\widehat{x})+\frac{{\widehat{c}}_4\sqrt{\widehat{\rho }}}{\sqrt{{\widehat{c}}_1}}{f}_w(\Delta )+\kappa {(f_w(\Delta ))}^2\hfill \end{array}$$

(32)

Therefore, if ${\Omega }_{{\widehat{\rho }}_e}$ is chosen as a level set of $\widehat{V}$ that satisfies Equation (31), $\widehat{V}(x)$ based on the actual state $x(t)$ is bounded in ${\Omega }_{\widehat{\rho }}$ for all $t\in [t_k,t_{k+1})$. However, if $x(t_k)\in {\Omega }_{\widehat{\rho }}\backslash {\Omega }_{{\widehat{\rho }}_e}$, the constraint of Equation (19f) is activated such that the control action u decreases the value of $\widehat{V}(\widehat{x})$ based on the states predicted by the RNN model of Equation (19b) within the next sampling period. According to Equation (25b) in Proposition 3, the value of $\widehat{V}$ also decreases along the state trajectory of the actual nonlinear system of Equation (1) over $t\in [t_k,t_{k+1})$. Therefore, it is concluded that for any initial condition in ${\Omega }_{\widehat{\rho }}$, the closed-loop state of the system of Equation (1) is bounded in ${\Omega }_{\widehat{\rho }}$ for all times under the LEMPC of Equation (19). □

5. Application to a Chemical Process Example

A chemical process example is used to illustrate the application of LEMPC using an ensemble of RNN models, under which the state of the closed-loop system is maintained within the stability region in state-space for all times. We consider a well-mixed, non-isothermal continuous stirred tank reactor (CSTR), in which an irreversible second-order exothermic reaction takes place. Specifically, the reaction transforms a reactant A to a product B ($A\to B$). The inlet temperature, the inlet concentration of A, and the feed volumetric flow rate of the reactor are $T_0$, $C_{A0}$ and F, respectively. Additionally, the CSTR is equipped with a heating jacket to supply or remove heat at a rate Q. The following material and energy balance equations are used to represent the CSTR dynamic model:

$${\displaystyle \frac{d{C}_A}{dt}}={\displaystyle \frac{F}{V}}(C_{A0}-C_A)-k_0{e}^{\frac{-E}{RT}}{C}_A^2$$

(33a)

$${\displaystyle \frac{dT}{dt}}={\displaystyle \frac{F}{V}}(T_0-T)+{\displaystyle \frac{-\Delta H}{{\rho }_L{C}_p}}{k}_0{e}^{\frac{-E}{RT}}{C}_A^2+{\displaystyle \frac{Q}{{\rho }_L{C}_{p}V}}$$

(33b)

where V is the volume of the reacting liquid in the reactor, $C_A$ is the reactant A concentration in the reactor, Q is the heat input rate and T is the temperature of the reactor. The feed volumetric flow rate and temperature are F and $T_0$, respectively. The feed concentration of reactant A is $C_{A0}$. The reacting liquid has a constant heat capacity of $C_p$ and density of ${\rho }_L$. $k_0$ is the pre-exponential constant, $\Delta H$ is the enthalpy of reaction, R is the ideal gas constant, and E is the activation energy. Process parameter values are listed in

Table 1. Parameter values of the CSTR.

$T_0$ = 300 K	$F=5{\mathrm{m}}^3/\mathrm{h}$
$V=1{\mathrm{m}}^3$	$E=5\times {10}^4\mathrm{kJ}/\mathrm{kmol}$
$k_0=8.46\times {10}^6{\mathrm{m}}^3/\mathrm{kmol}\mathrm{h}$	$\Delta H=-1.15\times {10}^4\mathrm{kJ}/\mathrm{kmol}$
$C_p=0.231\mathrm{kJ}/\mathrm{kg}\mathrm{K}$	$R=8.314\mathrm{kJ}/\mathrm{kmol}\mathrm{K}$
${\rho }_L=1000{\mathrm{kg}/\mathrm{m}}^3$	$C_{A{0}_s}=4{\mathrm{kmol}/\mathrm{m}}^3$
$Q_s=0.0\mathrm{kJ}/\mathrm{h}$	$C_{A_s}=1.2{\mathrm{kmol}/\mathrm{m}}^3$
$T_s=438\mathrm{K}$

.

The CSTR is initially operated at the unstable steady-state $(C_{As},T_s)=(1.95{\mathrm{kmol}/\mathrm{m}}^3,402\mathrm{K})$, and $(C_{A{0}_s}{Q}_s)=(4{\mathrm{kmol}/\mathrm{m}}^3,0\mathrm{kJ}/\mathrm{h})$. In this example, the two manipulated inputs are the heat input rate, $\Delta Q=Q-Q_s$, and the inlet concentration of species A, $\Delta C_{A0}=C_{A0}-C_{A{0}_s}$, respectively. Additionally, the manipulated inputs are subject to the following constraints: $|\Delta Q|\le 5\times {10}^5\mathrm{kJ}/\mathrm{h}$ and $|\Delta C_{A0}|\le 3.5{\mathrm{kmol}/\mathrm{m}}^3$. The states and the inputs of the system of Equation (33) are represented by $x^T=[C_A-C_{As}T-T_s]$ and $u^T=[\Delta C_{A0}\Delta Q]$, respectively. Therefore, the unstable steady-state of the system of Equation (33) is at the origin of the state-space, (i.e., $(x_s^{*},u_s^{*})=(0,0)$).

The control objective of the LEMPC of Equation (19) is to maximize the production rate of B in the CSTR by manipulating the inlet concentration $\Delta C_{A0}$ and the heat input rate $\Delta Q$, while maintaining the closed-loop state trajectories in the stability region ${\Omega }_{\widehat{\rho }}$ for all times under LEMPC. The objective function of the LEMPC optimizes the production rate of B as follows:

$$l_e(\tilde{x},u)=k_0{e}^{-E/RT}{C}_A^2$$

(34)

To simulate the dynamic model of Equation (33) numerically, the explicit Euler method is used with an integration time step of $h_c={10}^{-4}\mathrm{h}$. The Python module of the IPOPT software package [27], named PyIpopt, is used to solve the nonlinear optimization problem of the LEMPC of Equation (19) with a sampling period $\Delta ={10}^{-2}\mathrm{h}$. Additionally, a Message Passing Interface (MPI) for the Python programming language, named MPI4Py ([28]), is incorporated in the LEMPC optimization problem to execute multiple RNN predictions of Equation (19b) concurrently in independent computing processes.

Extensive open-loop simulations for the CSTR system of Equation (33) are initially conducted to generate the dataset for RNN models. The control Lyapunov function $\widehat{V}(x)=x^{T}Px$ is designed with the following positive definite P matrix:

$$P=\left[\begin{array}{cc}1060& 22\\ 22& 0.52\end{array}\right]$$

(35)

Two hidden recurrent layers consisting of 96 and 64 recurrent units, respectively, are used in the RNN model. Additionally, a 10-fold cross validation are used to construct an ensemble of RNN models for the LEMPC of Equation (19). The closed-loop stability region ${\Omega }_{\widehat{\rho }}$ with $\widehat{\rho }=368$ and a subset ${\Omega }_{{\widehat{\rho }}_e}$ with ${\widehat{\rho }}_e=320$ for the CSTR system are characterized based on the obtained RNN models and the stabilizing controller $u={\mathsf{\Phi }}_{nn}(x)\in U$. To demonstrate the effectiveness of the proposed LEMPC of Equation (19), a linear state-space model (i.e., $\dot{x}=Ax+Bu$ with $A=100\times \left[-0.154-0.003;5.190.138\right]$ and $B=\left[4.030;1.230.004\right]$) is also identified following the system identification method in [29], and is added into the following performance comparison under the ensemble of RNN models and the first-principles model of Equation (33).

The open-loop simulation studies are first carried out for the first-principles model of Equation (33), an RNN model ensemble, and the linear state-space model, respectively. The optimal number of RNN models used in the ensemble is determined to be five from extensive simulation tests due to its best performance in terms of closeness to the trajectories under the first-principles model. As shown in Figure 5, starting from various initial conditions in ${\Omega }_{\widehat{\rho }}$, the state trajectories under the RNN model ensemble almost overlap those under the first-principles model using the same input sequences, while the state trajectories under the linear state-space model show considerable deviation, especially in the top and bottom regions of ${\Omega }_{\widehat{\rho }}$. From Figure 5, it is demonstrated that the RNN model ensemble is a good representation for the first-principles model of Equation (33) in that it successfully captures the CSTR nonlinear dynamics everywhere in ${\Omega }_{\widehat{\rho }}$.

Subsequently, the closed-loop simulation of the nominal CSTR system (i.e., $w(t)\equiv 0$) under the first-principles model of Equation (33), the ensemble of RNN models, and the linear state-space model, respectively, are shown in Figure 6. The following material constraint is used in the LEMPC of Equation (19) to make the averaged reactant material available within the entire operating period $t_p$ to be its steady-state value, $C_{A0s}$ (i.e., the averaged reactant material in deviation form, $u_1$, is equal to 0).

$$\frac{1}{t_p}{\int }_0^{t_p}{u}_1(\tau )d\tau =0{\mathrm{kmol}/\mathrm{m}}^3$$

(36)

In Figure 6, it is shown that for the initial condition (0,0), the state trajectory under the ensemble of RNN models is close to that under the first-principles model, and both state trajectories are bounded in ${\Omega }_{{\widehat{\rho }}_e}$ in the absence of disturbances. However, it is observed that the state trajectory under the linear state-space model ultimately leaves ${\Omega }_{\widehat{\rho }}$ due to its large model mismatch. This implies that a more conservative set ${\Omega }_{{\widehat{\rho }}_e}$ needs to be characterized for the state spaced model to guarantee the boundedness of state in the stability region ${\Omega }_{\widehat{\rho }}$.

Additionally, Figure 7 and Figure 8 show the manipulated inputs profiles in deviation from the steady-state values. It is shown that the manipulated inputs are bounded within the input constraints for all times. Moreover, the material constraint of Equation (36) is satisfied as it is shown in Figure 7 that the LEMPC using the first-principles model and the RNN model ensemble both consume the maximum allowable $\Delta C_{A0}$ at the first few sampling steps, and thus have to lower the consumption of $\Delta C_{A0}$ near the end of operating period. The LEMPC using the linear state-space model also satisfies the material constraint, but shows persistent oscillation due to model mismatch. Additionally, in Figure 8, it is shown that the consumption of $\Delta Q$ under the RNN ensemble is close to that under the first-principles model (both correspond to left y-axis) since their closed-loop trajectories shown in Figure 6 are similar. However, $\Delta Q$ shows large oscillation under the LEMPC using the linear state-space model (right y-axis) since the closed-loop trajectory in Figure 6 leaves ${\Omega }_{{\widehat{\rho }}_e}$, and thus, the contractive constraint of Equation (19f) is activated frequently.

To demonstrate that the closed-loop system under LEMPC achieves high process economics than the steady-state operation using the same amount reactant $C_{A0}$, we calculate the accumulated economic benefits $L_E={\int }_0^{t_p}{l}_e(x,u)dt$ over the entire operating period $t_p$ = 0.2 h, which are $2.04$, $4.14$, and $4.22$ for the steady-state operation, the LEMPC under an ensemble of RNN models, and the LEMPC under the first-principles model, respectively. Therefore, it is demonstrated that both closed-loop stability and economic optimality are achieved under the proposed LEMPC using an ensemble of RNN models.

6. Conclusions

This work presented a new class of economic model predictive controllers that use an ensemble of recurrent neural network models as the prediction model for nonlinear systems. Specifically, an RNN ensemble was first developed based on the dataset of extensive open-loop simulations within the operating region. Parallel computing was used to reduce the computation time of prediction by multiple RNN models. Then, the LEMPC that incorporates RNN models and Lyapunov-based stability constraints was formulated to stabilize a nonlinear process within the closed-loop stability region while optimizing process economic benefits simultaneously. The proposed LEMPC using an RNN ensemble was applied to a nonlinear chemical process example to demonstrate its economic optimality and closed-loop stability.