Machine Learning‐based Distributed Model Predictive Control Of .

3 of 18CHEN ET AL.such that the origin of the nominal system of Equation (1) with3LS TM NETWO RK w(t) 0 is rendered exponentially stable in the sense that thereexists a C1 Control Lyapunov function V(x) such that the followingIn this work, we develop an LSTM network model with the follow-inequalities hold for all x in an open neighborhood D around theing form:origin: x Fnn ð x, u1 , u2 Þ A x ΘT y22c1 jxj V ðxÞ c2 jxj ,ð5Þð2aÞwhere x Rn is the predicted state vector, and u1 Rm1 and u2 Rm2 are V ðxÞFðx, Φ1 ðxÞ, Φ2 ðxÞ, 0Þ c3 jxj2 , xð2bÞ V ðxÞ x c4 j x j yhT y1 , , yn , yn 1 , , yn m1 , yn m1 1 , yni m1 m2 , yn m1 m2 1 Hð x1 Þ, , Hð xn Þ,u11 , , u1m1 , u21 , , u2m2 , 1 Rn m1 m2 1 is a vec-ð2cÞtor of the network state x, where H( ) represents a series of tors.nonlinear activation functions in each LSTM unit, the inputs u1 and u2,and the constant 1 which accounts for the bias term. A is a diagonalthe nonlinear system of Equation (1). A set of candidate controllerscoefficient matrix, that is, A diag{ α1, , αn} Rn n, and Θ ½θ1 , , θn Rðn m1 m2 1Þ n with θi βi ωi1 , , ωiðn m1 m2 Þ , bi , i 1,Φ1 ðxÞ Rm1 and Φ2 ðxÞ Rm2 , both denoted by Φk(x) where k 1, 2, is , n. αi and βi are constants, and ωij is the weight connecting the jthgiven in the following form:input to the ith neuron where i 1, , n and j 1, , (n m1 m2), andwhere c1, c2, c3, and c4 are positive constants. F(x, u1, u2, w) representsϕki ðxÞ 8 : p pffiffiffiffiffiffiffiffiffiffiffiffiffiffip2 q4qTq q08 min uk iΦki ðxÞ ϕki ðxÞ : maxuk ibi is the bias term for i 1, , n. We use x to represent the state ofif q 0ð3aÞif q 0actual nonlinear system of Equation (1) and use x for the state of theLSTM model of Equation (5). Here, αi is assumed to be positive suchthat each state xi is bounded-input bounded-state stable.Instead of having one-way information flow from the input layer toif ϕki ðxÞ uminkið3bÞmaxif uminki ϕki ðxÞ ukiif ϕki ðxÞ umaxkithe output layer in a FNN, RNNs introduce feedback loops into the network and allow information exchange in both directions between modules. Unlike FNNs, RNNs take advantage of the feedback signals tostore outputs derived from past inputs, and together with the currentwhere k 1, 2 represents the two candidate controllers, p denotes LfV T(x) and q denotes Lgki V ðxÞ, f [f1 fn]T, gki gki1 , , gkin , (i 1, 2, ,input information, give a more accurate prediction of the current out-m1 for k 1 corresponding to the vector of control actions Φ1(x), andresenting dynamic behaviors of time-series samples, therefore it is anput. By having access to information of the past, RNN is capable of rep-i 1, 2, , m2 for k 2 corresponding to the vector of control actionseffective method used to model nonlinear processes. Based on the uni-Φ2(x)). ϕki ðxÞ of Equation (3a) represents the i th component of theversal approximation theorem, it can be shown that the RNN modelcontrol law ϕk(x). Φki ðxÞ of Equation (3) represents the ith componentwith sufficient number of neurons is able to approximate any nonlinearof the saturated control law Φk(x) that accounts for the input con-dynamic system on compact subsets of the state-space for finitestraints uk Uk.time.16,17 However, in a standard RNN model, the problem of vanishingBased on Equation (2), we can first characterize a region wheregradient phenomena often arises due to the network's difficulty to cap-the time-derivative of V is rendered negative under the controllernΦ1(x) U1, Φ2(x) U2 as D x Rn j V ðxÞ Lf V Lg1 Vu1 Lg2 Vu2 ture long term dependencies; this is because of multiplicative gradients2that can be exponentially decaying with respect to the number of c3 jxj ,u1 Φ1 ðxÞ U1 , u2 Φ2 ðxÞ U2 g [ f0g . Then the closed-looplayers. Therefore, the stored information over extended time intervals isstability region Ωρ for the nonlinear system of Equation (1) isvery limited in a short term memory manner. Due to these consider-defined as a level set of the Lyapunov function, which is inside D:ations, Hochreiter and Schmidhuber18 proposed the LSTM network,Ωρ {x D j V(x) ρ}, where ρ 0 and Ωρ D. Also, the Lipschitzwhich is a type of RNN that uses three gated units (the forget gate, theproperty of F(x, u 1, u2 , w) combined with the bounds on u1 , u 2, andinput gate, and the output gate) to protect and control the memory cell00w implies that there exist positive constants M, Lx , Lw , Lx , Lw such0state, c(k), where k 1, , T, such that information will be stored andthat the following inequalities hold 8x, x Ωρ, u1 U1, u2 U2remembered for long periods of time.18 For these reasons, LSTM net-and w W:works may perform better when modeling processes where inputstoward the beginning of a long time-series sequence are crucial to thej Fðx, u1 , u2 ,wÞ j Mð4aÞj F ðx, u1 ,u2 , wÞ F ðx0 , u1 ,u2 , 0Þ j Lx j x x0 j Lw j w jð4bÞprediction of outputs toward the end of the sequence. This may bemore prevalent in large-scale systems where there may exist inherent V ðxÞ 00 V ðx0 Þ 00 x F ðx, u1 , u2 , wÞ x Fðx , u1 , u2 ,0Þ Lx j x x j Lw j w jtime delays between subsystems, causing discrepancies in the speed ofthe dynamics between the subsystems. The basic architecture of anð4cÞLSTM network is illustrated in Figure 1a. We develop an LSTM networkmodel to approximate the class of continuous-time nonlinear processes

4 of 18CHEN ET AL.F I G U R E 1 Structure of (a) theunfolded LSTM network and (b) theinternal setup of an LSTM unit [Colorfigure can be viewed atwileyonlinelibrary.com](a)(b)of Equation (1). We use m Rðn m1 m2 Þ T to denote the matrix ofinput sequences to the LSTM network, and x Rn Tto denote thematrix of network output sequences. The output from each repeatingmodule that is passed onto the next repeating module in theunfolded sequence is the hidden state, and the vector of hidden hf ðkÞ σ ωmf mðk Þ ωf hðk 1Þ bfð6bÞ hcðkÞ f ðkÞcðk 1Þ iðkÞtanh ωmc mðkÞ ωc hðk 1Þ bcð6cÞ hoðkÞ σ ωmo mðk Þ ωo hðk 1Þ boð6dÞhðkÞ oðkÞtanhðcðkÞÞð6eÞ xðkÞ ωy hðkÞ byð6fÞstates is denoted by h. The network output x at the end of the prediction period is dependent on all internal states h(1), , h(T), wherethe number of internal states T (i.e., the number of repeating modules)corresponds to the length of the time-series input sample. The LSTMnetwork calculates a mapping from the input sequence m to the output sequence x by calculating the following equations iteratively fromk 1 to k T:where σ( ) is the sigmoid function, tanh( ) is the hyperbolic tangentfunction; both of which are activation functions. h(k) is the internal hiðkÞ σ ωmi mðk Þ ωi hðk 1Þ bið6aÞstate, and xðkÞ is the output from the repeating LSTM module with ωyand by denoting the weight matrix and bias vector for the output,

5 of 18CHEN ET AL.respectively. The outputs from the input gate, the forget gate, and thethrough numerical approximation using the forward finite differenceoutput gate are represented by i(k), f(k), o(k), respectively; correspond-method. Given that the time interval between internal states of thehωmi , ωi ,are the weight matrices for the input vec-LSTM model is a multiple of the integration time step qnn hc, the timetor m and the hidden state vectors h within the input gate, the forgetgate, and the output gate, respectively, and bi, bf, bo represent the biasderivative of the LSTM predicted state xðtÞ at t tk can be approxi mated by x ðtk Þ F nn ðxðtk Þ, u1 , u2 Þ xðtk qnn hc Þ xðtk Þ. The time derivative ofvectors within each of the three gates, respectively. Furthermore, c(k)the actual state x(t) at t tk can be approximated byis the cell state which stores information to be passed down the net-x ðtk Þ Fðxðtk Þ, u1 , u2 , 0Þ xðtk qqnnnnhhc cÞ xðtk Þ . At time t tk, xðtk Þ xðtk Þ , theingly,hωmf , ωf ,work units, withωmc ,hωmo , ωoωhc ,and bc representing the weight matrices forqnn hcconstraint jν j γ j xj can be written as follows:the input and hidden state vectors, and the bias vector in the cell statej ν j j F ðxðtk Þ,u1 , u2 , 0Þ F nn ðxðtk Þ, u1 , u2 Þ jactivation function, respectively. The series of interacting nonlinearð7aÞfunctions carried out in each LSTM unit, outlined in Equation (6), canbe represented by Hð xÞ. The internal structure of a repeating module jwithin an LSTM network where the iterative calculations ofEquation (6) are carried out is shown in Figure 1b.xðtk qnn hc Þ xðtk qnn hc Þjqnn hc γ j xðtk Þ jThe closed-loop simulation of the continuous-time nonlinear sys-ð7bÞð7cÞtem of Equation (1) is carried out in a sample-and-hold manner, wherethe feedback measurement of the closed-loop state x is received byxðtk qnn hc Þwhich will be satisfied if j xðtk qnn hcxÞð j γqnn hc . Therefore, thetk Þthe controller every sampling period Δ. Furthermore, state informa-mean absolute percentage error between the predicted states x andtion of the simulated nonlinear process is obtained via numerical inte-the targeted states x in the training data will be used as a metric togration methods, for example, explicit Euler, using an integration timeassess the modeling error of the LSTM model. While the error boundsstep of hc. Since the objective of developing the LSTM model is itsthat the LSTM network model and the actual process should satisfy toeventual utilization in a controller, the prediction period of the LSTMensure closed-loop stability are difficult to calculate explicitly and are,model is set to be the same as the sampling period Δ of the modelin general, conservative, they provide insight into the key networkpredictive controller. The time interval between two consecutiveparameters that will need to be tuned to reduce the error betweeninternal states within the LSTM can be chosen to be a multiple qnn ofthe two models as well as the amount of data needed to build a suit-the integration time step hc used in numerical integration of theable LSTM model.nonlinear process, with the minimum time interval being qnn 1, thatIn order to gather adequate training data to develop the LSTMis, 1 hc. Therefore, depending on the choice of qnn, the number ofmodel for the nonlinear process, we first discretize the desired operat-internal states, T, will follow T qnnΔ hc . Given that the input sequencesing region in state-space with sufficiently small intervals as well as dis-fed to the LSTM network are taken at time t tk, the future statescretize the range of manipulated inputs based on the control actuatorpredicted by the LSTM network, xðtÞ, at t tk Δ, would be the net-limits. We run open-loop simulations for the nonlinear process ofwork output vector at k T, i.e., xðtk ΔÞ xðT Þ . The LSTM learningEquation (1) starting from different initial conditions inside the desiredalgorithm is developed to obtain the optimal parameter matrix Γ*,operating region, that is, x0 Ωρ, for finite time using combinations ofwhich includes the network parameters ωi, ωf, ωc, ωo, ωy, bi, bf, bc, bo,the different manipulated inputs u1 U1, u2 U2 applied in a sample-by. Under this optimal parameter matrix, the error between the actualand-hold manner, and the evolving state trajectories are recorded atstate x(t) of the nominal system of Equation (1) (i.e., w(t) 0) and thetime intervals of size qnn hc. We obtain enough samples of such tra-modeled states xðtÞ of the LSTM model of Equation (5) is minimized.jectories to sweep over all the values that the states and the manipu-The LSTM model is developed using a state-of-the-art applicationlated inputs (x, u1, u2) could take to capture the dynamics of theprogram interface, that is, Keras, which contains open-source neuralprocess. These time-series data can be separated into samples with anetwork libraries. The mean absolute percentage error between x(t)fixed length T, which corresponds to the prediction period of theand xðtÞ is minimized using the adaptive moment estimation optimizer,LSTM model, where Δ T qnn hc. The time interval between twothat is, Adam in Keras, in which the gradient of the error cost functiontime-series data points in the sample qnn hc corresponds to the timeis evaluated using back-propagation. Furthermore, in order to ensureinterval between two consecutive memory units in the LSTM net-that the trained LSTM model can sufficiently represent the nonlinearwork. The generated dataset is then divided into training and valida-process of Equation (1), which in turn ascertains that the LSTM modeltion sets.can be used in a model-based controller to stabilize the actualnonlinear process at its steady-state with guaranteed stability proper-Remark 1 The actual nonlinear process is a continuous-time modelties, a constraint on the modeling error is also imposed during training,that can be represented using Equation (1); therefore, to char-where jν j j F(x, u1, u2, 0) Fnn(x, u1, u2) j γ j xj, with γ 0. Addition-acterize the modeling error ν between the LSTM network andally, to avoid over-fitting of the LSTM model, the training process isthe nonlinear process of Equation (1), the LSTM network isterminated once the modeling error falls below the desired thresholdrepresented as a continuous-time model of Equation (5). How-and the error on the validation set stops decreasing. One way toever, the series of interacting nonlinear operations in the LSTMassess the modeling error ν F(x(tk), u1, u2, 0) Fnn(x(tk), u1, u2) ismemory unit is carried out recursively akin to a discrete-time

6 of 18CHEN ET AL.model. The time interval qnn hc between two LSTM memoryset of control actions u1(tk) U1\{0} and u2(tk) U2\{0} that are appliedunits is given by the time interval between two consecutivetime-series data points in the training samples. Since the LSTMin a sample-and-hold fashion for the time interval t [tk, tk hc) (hc isthe integration time step), f and g can be numerically approximated asnetwork provides a predicted state at each time intervalfollows:qnn hc calculated by each LSTM memory unit, similarly tohow we can use numerical integration methods to obtain the f ðxðtk ÞÞ state at the same time instance using the continuous-timemodel, we can use the predicted states from the LSTM network to compare with the predicted states from the nonlinearmodel of Equation (1) to assess the modeling error. The model-Ð tk hcing error is subject to the constraint of Equation (7) to ensurethat the LSTM model can be used in the model-based controller with guaranteed stability properties.tk g1 ðxðtk ÞÞ tkF nn ðx, 0,0Þdt xðtk Þð9aÞhcF nn ðx, u1 ðtk Þ, 0Þdt Ð tk hctkF nn ðx, 0,0Þdthc u1 ðtk ÞÐ tk hc g2 ðxðtk ÞÞ The integralÐ tk hctkÐ tk hctkF nn ðx, 0,u2 ðtk ÞÞdt Ð tk hctkF nn ðx, 0,0Þdthc u2 ðtk Þð9bÞð9cÞF nn ðx, u1 , u2 Þdt gives the predicted state xðtÞ att tk hc under the sample-and-hold implementation of the inputs3.1 Lyapunov-based control using LSTM modelsu1(tk) and u2(tk); xðtk hc Þ is the first internal state of the LSTM network, given that the time interval between consecutive internal statesOnce we obtain an LSTM model with a sufficiently small olleru1 Φnn1 ðxÞ U1 and u2 Φnn2 ðxÞ U2 that can render the origin of theLSTM model of Equation (5) exponentially stable in an open neighbor around the origin in the sense that there exists a C1 Controlhood D ðxÞ such that the following inequalities hold forLyapunov function V all x in D:of the LSTM network is chosen as the integration time step hc. Afterobtaining f , g1 , and g2 , the stabilizing control law Φnn ðxÞ and Φnn ðxÞ12can be computed similarly as in Equation (3), where f, g1, and g2 are can also bereplaced by f , g1 , and g2 , respectively. Subsequently, Vcomputed using the approximated f , g1 , and g2 . The assumptions ofEquations (2) and (8) are the stabilizability requirements of the firstprinciples model of Equation (1) and the LSTM network model ofEquation (5), respectively. Since the dataset for developing the LSTM ðxÞ c2 jxj2 , c1 jxj V2ð8aÞnetwork model is generated from open-loop simulations for x Ωρ,u1 U1, and u2 U2, the closed-loop stability region of the LSTM sys- ðxÞ VF nn ðx, Φnn1 ðxÞ, Φnn2 ðxÞÞ c3 jxj2 , x V ðxÞ c jxj x 4ð8bÞtem is a subset of the closed-loop stability region of the actualnonlinear system, Ω ρ Ωρ . Additionally, there exist positive constantsMnn and Lnn such that the following inequalities hold for all x, x0 Ω ρ ,ð8cÞu1 U1 and u2 U2:j F nn ðx, u1 , u2 Þ j Mnnwhere c1 , c2 , c3 , c4 are positive constants, and Fnn(x, u1, u2) representsthe LSTM network model of Equation (5). Similar to the characterization method of the closed-loop stability region Ωρ for the nonlinear V ðx0 Þ V ðxÞ Fnn ðx, u1 ,u2 Þ Fnn ðx0 , u1 , u2 Þ Lnn j x x0 j x xð10aÞð10bÞsystem of Equation (1), we first search the entire state-space to char where the following inequality holds:acterize a set of states D ðxÞ V V ðxÞ x F nn ðx,u1 , u2 Þ c3 jxj2 , u1 Φnn1 ðxÞ U1 , u2 Φnn2 ðxÞ U2 .The closed-loop stability region for the LSTM network model of Equa :tion (5) is defined as a level set of Lyapunov function inside Dno Ωρ x D j V ðxÞ ρ , where ρ 0. Starting from Ω ρ , the origin of the4 DISTRIBUTED LMPC USING LSTMNETWORK MODELSLSTM network model of Equation (5) can be rendered exponentiallyTo achieve better closed-loop control performance, some level ofstable under the controller u1 Φnn1 ðxÞ U1 , and u2 Φnn2 ðxÞ U2 . It iscommunication may be established between the different controllers.noted that the above assumption of Equation (8) is the same as theIn a distributed LMPC framework, we design two separate LMPCs—assumption of Equation (2) for the general class of nonlinear systemsLMPC 1 and LMPC 2—to compute control actions u1 and u2, respec-of Equation (1) since the LSTM network model of Equation (5) can bewritten in the form of Equation (1) (i.e., x f ð xÞ gð xÞu, where f ð Þ andtively; the trajectories of control actions computed by LMPC 1 and gð Þ are obtained from coefficient matrices A and Θ in Equation (5)).types of distributed control architectures: sequential and iterative dis-However, due to the complexity of the LSTM structure and the interactions of the nonlinear activation functions, f and g may be hard totributed MPCs. Having only one-way communication, the sequentialcompute explicitly. For computational convenience, at t tk, given atransmitting inter-controller signals. However, it must assume theLMPC 2 are denoted by ud1 and ud2 , respectively. We consider twodistributed MPC architecture carries less computational load with

7 of 18CHEN ET AL.F I G U R E 2 Schematic diagrams showing theflow of information of (a) the sequentialdistributed LMPC and (b) the iterative distributedLMPC systems with the overall process [Colorfigure can be viewed at wileyonlinelibrary.com](a)(b)input trajectories along the prediction horizon of the other controllerssampling period u d2 ðtk Þ to the corresponding actuators, and sends thedownstream of itself in order to make a decision. The iterative distrib-entire optimal trajectory to LMPC 1.uted MPC system allows signal exchanges between all controllers,3. LMPC 1 receives the entire optimal input trajectory of ud2 fromthereby allowing each controller to have full knowledge of theLMPC 2, and eva

and applied to distributed MPC. As distributed MPC systems also depend on an accurate process model, the development and implemen-tation of RNN models in distributed MPCs is an important area yet to be explored. In the present work, we introduce distributed control frameworks that employ a LSTM network, which is a particular type of RNN. The .