Data Driven Economic Model Predictive Control

This manuscript addresses the problem of data driven model based economic model predictive control (MPC) design. To this end, first, a data-driven Lyapunov-based MPC is designed, and shown to be capable of stabilizing a system at an unstable equilibrium point. The data driven Lyapunov-based MPC utilizes a linear time invariant (LTI) model cognizant of the fact that the training data, owing to the unstable nature of the equilibrium point, has to be obtained from closed-loop operation or experiments. Simulation results are first presented demonstrating closed-loop stability under the proposed data-driven Lyapunov-based MPC. The underlying data-driven model is then utilized as the basis to design an economic MPC. The economic improvements yielded by the proposed method are illustrated through simulations on a nonlinear chemical process system example.


Introduction
Control systems designed to manage chemical process operations often face numerous challenges such as inherent nonlinearity, process constraints and uncertainty.Model predictive control (MPC) is a well-established control method that can handle these challenges.In MPC, the control action is computed by solving an open-loop optimal control problem at each sampling instance over a time horizon, subject to the model that captures the dynamic response of the plant, and constraints [1].In early MPC designs, the objective function was often utilized as a parameter to ensure closed-loop stability.In subsequent contributions, Lyapunov-based MPC was proposed where feasibility and stability from a well characterized region was built into the MPC [2,3].
With increasing recognition (and ability) of MPC designs to focus on economic objectives, the notion of Economic MPC (EMPC) was developed for linear and nonlinear systems [4][5][6], and several important issues (such as input rate-of-change constraint and uncertainty) addressed.The key idea with the EMPC designs is the fact that the controller is directly given the economic objective to work with, and the controller internally determines the process operation (including, if needed, a set point) [7].
Most of the existing MPC formulations, economic or otherwise, have been illustrated using first principles models.With growing availability of data, there exists the possibility of enhancing MPC implementation for situations where a first principles model may not be available, and simple 'step-test', transfer-function based model identification approaches may not suffice.One of the widely utilized approaches in the general direction of model identification are latent variable methods, where the correlation between subsequent measurements is used to model and predict the process evolution [8,9].In one direction, Dynamic Mode Decomposition with control (DMDc) has been utilized to extract low-order models from high-dimensional, complex systems [10,11].In another direction, subspace-based system identification methods have been adapted for the purpose of model identification, where state-space model from measured data are identified using projection methods [12][13][14].To handle the resultant plant model mismatch with data-driven model based approaches, monitoring of the model validity becomes especially important.
One approach to monitor the process is to focus on control performance [15], where the control performance is monitored and compared against a benchmark control design.To focus more explicitly on the model behavior, in a recent result [16], an adaptive data-driven MPC was proposed to evaluate model prediction performance and trigger model identification in case of poor model prediction.In another direction, an EMPC using empirical model was proposed [17].The approach relies on a linearization approach, resulting in closed-loop stability guarantees for regions where the plant-model mismatch is sufficiently small, and illustrate results on stabilization around nominally stable equilibrium points.In summary, data driven MPC or EMPC approaches, which utilize appropriate modeling techniques to identify data from closed-loop tests to handle operation around nominally unstable equilibrium points, remain to be addressed.
Motivated by the above considerations, in this work, we address the problem of data driven model based predictive control at an unstable equilibrium point.In order to identify a model around an unstable equilibrium point, the system is perturbed under closed-loop operation.Having identified a model, a Lyapunov-based MPC is designed to achieve local and practical stability.The Lyapunov-based design is then used as the basis for a data driven Lyapunov-based EMPC design to achieve economical goals while ensuring boundedness.The rest of the manuscript is organized as follows: first, the general mathematical description for the systems considered in this work and a representative formulation for Lyapunov-based model predictive control are presented.Then, the proposed approach for closed-loop model identification is explained.Subsequently, a Lyapunov-based MPC is designed and illustrated through a simulation example.Finally, an economic MPC is designed to consider economical objectives.The efficacy of the proposed method is illustrated through implementation on a nonlinear continuous stirred-tank reactor (CSTR) with input rate of change constraints.Finally, concluding remarks are presented.

Preliminaries
This section presents a brief description of the general class of processes that are considered in this manuscript, followed by closed-loop subspace identification and Lyapunov based MPC formulation.

System Description
We consider a multi-input multi-output (MIMO) controllable systems where u ∈ R n u denotes the vector of constrained manipulated variables, taking values in a nonempty convex subset U ⊂ R n u , where U = u ∈ R n u | u min ≤ u ≤ u max , u min ∈ R n u and u max ∈ R n u denote the lower and upper bounds of the input variables, and y ∈ R n y denotes the vector of measured output variables.In keeping with the discrete implementation of MPC, u is piecewise constant and defined over an arbitrary sampling instance k as: where ∆t is the sampling time and x k and y k denote state and output at the kth sample time.The central problem that the present manuscript addresses is that of designing a data driven modeling and control design for economic MPC.

System Identification
In this section, a brief review of a conventional subspace-based state space system identification methods is presented [16,18,19].These methods are used to identify the system matrices for a discrete-time linear time invariant (LTI) system of the following form: where x ∈ R n x and y ∈ R n y denote the vectors of state variables and measured outputs, and w ∈ R n x and v ∈ R n y are zero mean, white vectors of process noise and measurement noise with the following covariance matrices: where Q ∈ R n x ×n x , S ∈ R n x ×n y and R ∈ R n y ×n y are covariance matrices, and δ ij is the Kronecker delta function.The subspace-based system identification techniques utilize Hankel matrices constructed by stacking the output measurements and manipulated variables as follows: where i is a user-specified parameter that limits the maximum order of the system (n), and, j is determined by the number of sample times of data.By using Equation ( 4), the past and future Hankel matrices for input and output are defined: Similar block-Hankel matrices are made for process and measurement noises V p , V f ∈ R in y ×j and W p , W f ∈ R in x ×j are defined in the similar way.The state sequences are defined as follows: Furthermore, these matrices are used in the algorithm: By recursive substitution into the state space model equations Equations ( 1) and (2), it is straightforward to show: where: Equation ( 9) can be rewritten in the following form to have the input and output data at the left hand side of the equation [20]: In open loop identification methods, in the next step, by orthogonal projecting of Equation ( 15) onto Ψ p : Note that, the last two terms in RHS of Equation ( 15) are eliminated since the noise terms are independent, or othogonal to the future inputs.Equation (16) indicates that: Therefore, Γ i and H d i can be calculated using Equation ( 17) by decomposition methods.These can in turn be utilized to determine the system matrices (some of these details are deferred to Section 3.1).For further discussion on system matrix extraction, the readers are referred to references [18,19].

Lyapunov-Based MPC
The Lyapunov-based MPC (LMPC) for linear system has the following form: min ũk ,..., ũk+P subject to: where xk+j , ỹk+j , y SP k+j and ũk+j denote predicted state and output, output set-point and calculated manipulated input variables j time steps ahead computed at time step k, and xl is the current estimation of state, and 0 < α < 1 is a user defined parameter.The operator ||.|| 2 Q denotes the weighted Euclidean norm defined for an arbitrary vector x and weighting matrix W as ||x|| 2 W = x T Wx.Furthermore, Q y > 0 and R du ≥ 0 denote the positive definite and positive semi-definite weighting matrices for penalizing deviations in the output predictions and for the rate of change of the manipulated inputs, respectively.Moreover, N y and N u denote the prediction and control horizons, respectively, and the input rate of change, given by ∆ ũk+j = ũk+j − ũk+j−1 , takes values in a nonempty convex subset U • ⊂ R m , where U • = ∆u ∈ R n u | ∆u min ≤ ∆u ≤ ∆u max .Note finally that, while the system dynamics are described in continuous time, the objective function and constraints are defined in discrete time to be consistent with the discrete implementation of the control action.
Equations ( 23) and ( 24) are representatives of Lyapunov-based stability constraint [21,22], where V(x k ) is a suitable control Lyapunov function, and α, * > 0 are user-specified parameters.In the presented formulation, * > 0 enables practical stabilization to account for the discrete nature of the control implementation.
Remark 1. Existing Lyapunov-based MPC approaches exploit the fact that the feasibility (and stability) region can be pre-determined.The feasibility region, among other things, depends on the choice of the parameter α, the requested decay factor in the value of the Lyapunov function at each time step.If (reasonably) good first principles models are available, then these features of the MPC formulation provide excellent confidence over the operating region under closed-loop.In contrast, in the presence of significant plant-model mismatch (as is possibly the case with data driven models), the imposition of such decay constraints could result in unnecessary infeasibility issues.In designing the LMPC formulation with a data driven model, this possible lack of feasibility must be accounted for (as is done in Section 3.2).

Integrating Lyapunov-Based MPC with Data Driven Models
In this section, we first utilize an identification approach necessary to identify good models for operation around an unstable equilibrium point.The data driven Lyapunov-based MPC design is presented next.

Closed-Loop Model Identification
Note that, when interested in identifying the system around an unstable equilibrium point, open-loop data would not suffice.To begin with, nominal open-loop operation around an unstable equilibrium point is not possible.If the nominal operation is under closed-loop, but the loop is opened to perform step tests, the system would move to the stable equilibrium point corresponding to the new input value, thereby not providing dynamic information around the desired operating point.The training data, therefore, has to be obtained using closed-loop step tests, and an appropriate closed-loop model identification method employed.Such a method is described next.
In employing closed-loop data, note that the assumption of future inputs being independent of future disturbances no longer holds, and, if not recognized, can cause biased results in system identification [18].In order to handle this issue, the closed-loop identification approach in the projection utilizes a different variable Ψ pr instead of Ψ p .The new instrument variable, which satisfies the independence requirement, is used to project both sides of Equation ( 15) and the result is used to determine LTI model matrices.For further details, refer to [16,18,23].
By projecting Equation ( 15) onto Ψ pr we get: Since the future process and measurement noises are independent of the past input/output and future setpoint in Equation ( 25), the noise terms cancel, resulting in: By multiplying Equation ( 26) by the extended orthogonal observability Γ ⊥ i , the state term is eliminated: Therefore, the column space of Ψ f /Ψ pr is orthogonal to the row space of (Γ ⊥ i ) T −(Γ ⊥ i ) T Φ d i .By performing singular value decomposition (SVD) of Ψ f /Ψ pr : where Σ 1 contains dominant singular values of Ψ f /Ψ pr and, theoretically, it has the order n u i + n [18,23].Therefore, the order of the system can be determined by the number of the dominant singular values of the Ψ f /Ψ pr [20].
The orthogonal column space of where M ∈ R (n y −n)i×(n y −n)i is any constant nonsingular matrix and is typically chosen as an identity matrix [18,23].One approach to determine the LTI model is as follows [18]: From Equation ( 29), Γ i and Φ d i can be estimated: which results in (using MATLAB (2017a, MathWorks, Natick, MA, USA) matrix index notation): The past state sequence can be calculated as follows: The future state sequence can be calculated by changing data Hankel matrices as follows [18]: where Γi is obtained by eliminating the last n y rows of Γ i , and H d i is obtained by eliminating the last n y rows and the last n u columns of H d i .Then, the model matrices can be estimated using least squares: Note that the difference between the proposed method in [18] and described method is that, in order to ensure that the observer is stable (eigenvalues of A − KC are inside unit circle), instead of innovation form of LTI model, Equations ( 1) and ( 2) are used [16] to derive extended state space equations.The system matrices can be calculated as follows: With the proposed approach, process and measurement noise Hankel matrices can be calculated as the residual of the least square solution of Equation (39): Then, the covariances of plant noises can be estimated as follows: Model identification using closed-loop data has a positive impact on the predictive capability of the model (see the simulation section for a comparison with a model identified using open-loop data).

Control Design and Implementation
Having identified an LTI model for the system (with its associated states), the MPC implementation first requires a determination of the state estimates.To this end, an appropriate state estimator needs to be utilized.In the present manuscript, a Luenberger observer is utilized for the purpose of illustration.Thus, at the time of control implementation, state estimates xk are generated as follows: where L is the observer gain and is computed using pole placement method, and y k is the vector of measured variables (in deviation form, from the set point).
In order to stabilize the system at an unstable equilibrium point, a Lyapunov-based MPC is designed.The control calculation is achieved using a two-tier approach (to decouple the problem of stability enforcement and objective function tuning).The first layer calculates the minimum value of Lyapunov function that can be reached subject to the constraints.This tier is formulated as follows: where x, ỹ are predicted state and output and ũ1 is the candidate input computed in the first tier.x SP is underlying state setpoint (in deviation form from the nominal equilibrium point), which here is the desired unstable equilibrium point (and therefore zero in terms of deviation variables).For setpoint tracking, this value can be calculated using the target calculation method; readers are referred to [24] for further details.
Note that the first tier has a prediction horizon of 1 because the objective is to only compute the immediate control action that would minimize the value of the Lyapunov function at the next time step.V is chosen as a quadratic Lyapunov function with the following form: where P is a positive definite matrix computed by solving the Riccati equation with the LTI model matrices as follows: where Q ∈ R n x ×n x and R ∈ R n u ×n u are positive definite matrices.Then, in the second tier, this minimum value is used as a constraint (upper bound for Lyapunov function value at the next time step).The second tier is formulated as follows: subject to: where N p is the prediction horizon and ũ2 denotes the control action computed by the second tier.
In essence, in the second tier, the controller calculates a control action sequence that can take the process to the setpoint in an optimal fashion optimally while ensuring that the system reaches the minimum achievable Lyapunov function value at the next time step.Note that, in both of the tiers, the input sequence is a decision variable in the optimization problem, but only the first value of the input sequence of the second tier is implemented on the process.The solution of the first tier, however, is used to ensure and generate a feasible initial guess for the second tier.The two-tiered control structure is schematically presented in Figure 1.Remark 2. Note that Tiers 1 and 2 are executed in series and at the same time, and the implementation does not require a time scale separation.The overall optimization is split into two tiers to guarantee feasibility of the optimization problem.In particular, the first tier computes an input move with the objective function only focusing on minimizing the Lyapunov function value at the next time step.Notice that the constraints in the first tier are such that the optimization problem is guaranteed to be feasible.With this feasible solution, the second tier is used to determine the input trajectory that achieves the best performance, while requiring the Lyapunov function to decay.Again, since the second tier optimization problem uses the solution from Tier 1 to impose the stability constraint, feasibility of the second tier optimization problem, and, hence, of the MPC optimization problem, is guaranteed.In contrast, if one were to require the Lyapunov function to decay by an arbitrary chosen factor, determination of that factor in a way that guarantees feasibility of the optimization problem would be a non-trivial task.
Remark 3. It is important to recognize that, in the present formulation, feasibility of the optimization problem does not guarantee closed-loop stability.A superfluous (and incorrect) reason is as follows: the first tier computes the control action that minimizes the value of the Lyapunov function at the next step, but does not require that it be smaller than the previous time step, leading to potential destabilizing control action.The key point to realize here, however, is that if such a control action were to exist (that would lower the value of the Lyapunov function at the next time step), the optimization problem would determine that value by virtue of the Lyapunov function being the objective function, and lead to closed-loop stability.The reasons closed-loop stability may not be achieved are two: (1) the current state might be such that closed-loop stability is not achievable for the system dynamics and constraints; and ( 2) due to plant model mismatch, where the control action that causes the Lyapunov function to decay for the identified model does not do so for the system in question.The first reason points to a fundamental limitation due to the presence of input constraints, while the second is due to the lack of availability of the 'correct' system dynamics, and as such will be true in general for data driven MPC formulations.Note that inclusion of a noise/plant model mismatch term in the model may help with the predictive capability of the model, however, unless a bound on the uncertainty can be assumed, closed-loop stability can not be guaranteed.

Remark 4.
Along similar lines, consider the scenario where, based on the model, and constraints, an input value exists for which V(x(k)) <= V(x(k − 1)) is achievable.It can be readily shown that any solution computed by the first tier of the optimization problem would also result in V(x(k)) <= V(x(k − 1)) by virtue of the objective function being the Lyapunov function at the next time step.Thus, in such a case, the explicit incorporation of the constraint V(x(k)) <= V(x(k − 1)) (as is traditionally done in Lyapunov based MPC) does not help, and is not required.On the other hand, for the scenario where such an input does not exist, the inclusion of the constraint will cause the optimization problem to be infeasible.In contrast, in the proposed formulation, the MPC will compute a control action where the value of the Lyapunov function might be greater than the previous value, but greater by the smallest margin possible.The real impact of this phenomenon is in making the MPC formulation more pliable, especially when dealing with plant-model mismatch.In such scenarios, the proposed MPC continues to compute feasible (best possible, in terms of stabilizing behavior) solutions, and, should the process move into a region from where stabilization is possible, smoothly transits to computing stabilizing control action.
Remark 5.In the current manuscript, we focus on the cases where a first principal model is not available.If a good first principles model was available, it could be utilized directly in a nonlinear MPC design, or linearized if one were to implement a linear MPC.In the case of linearization, the applicability would be limited by the region over which the linearization holds.In contrast, note that the model utilized in the present manuscript does not result from a linearization of a nonlinear model.Instead, it is a linear model, possibly with a higher number of states than the original nonlinear model, albeit identified, and applicable, over a 'larger' region of operation, compared to a linearized model.

Remark 6.
To account for possible plant-model mismatch, model validity can be monitored with model monitoring methods [16], resulting in appropriately triggering re-identification in case of poor model prediction.
In another direction, in line with control performance monitoring approaches, the Lyapunov function value could be utilized.Thus, unacceptable increases in Lyapunov function value could be utilized as a means of triggering re-identification.

Remark 7.
As mentioned previously, in order to create rich training data around unstable operating points, closed-loop data must be generated.In turn, since open-loop methods result in biased estimation [25,26] in model identification, a suitable closed-loop identification method is utilized, and adapted to ensure that the model accurately captures the key dynamics.

Simulation Results
We next illustrate the proposed approach using a nonlinear CSTR example [27].To this end, consider a CSTR where a first-order, exothermic and irreversible reaction of the form A k − → B takes place.The mass and energy conservation laws results in the following mathematical model: The description of the process variables and the values of the system parameters are presented in Table 1.The control objective is to stabilize the system at an unstable equilibrium point using inlet concentration, C A 0 , and the rate of heat input, Q, while the manipulated inputs are constrained to be within the limits |C A 0 | ≤ 1 kmol/m 3 and |Q| ≤ 9 × 10 3 KJ/min, and the input rate is constrained as |∆C A 0 | ≤ 0.1 Kmol/m 3 and |∆Q| ≤ 9 × 200 KJ/min.We assume that both of the states are measured.The system has an unstable equilibrium point at C A = 0.573 Kmol/m 3 and T = 395.3K.The goal is to stabilize the system at this equilibrium point.To this end, first an LTI model is identified using closed-loop data; then, an MPC is designed to stabilize the system at the unstable equilibrium point.For system identification of the CSTR model, proportional-integral (PI) controllers (pairing C A with C A,in and T with Q) are implemented in the process.In particular, pseudo-random binary signals are used as set-points for PI controllers.The identified LTI model order is selected as n = 4 and i = 12, in order to achieve the best fit in model prediction (using cross-validation).Note that these four states are the states of the identified LTI model.When dealing with setpoint tracking, these states can be augmented with additional states and utilized as part of an offset-free MPC design.Model validation results under a different set of set-point changes from training data are presented in Figures 2 and 3.
The identified system is unstable with absolute eigenvalues 0.9311 0.9311 0.9998 1.0002 , which has an eigenvalue outside unit circle.The unstable nature of the identified model is consistent with the operation of the system around the unstable equilibrium point.For the model validation, initially, a steady state Kalman filter (gain calculated by the identification method) is utilized to update state estimate until t = 0.8 min and after convergence of the states (gaged via convergence of the outputs), the model and the input trajectory (without the state estimator) are used to predict the future output.Figure 2 illustrates the results of the model validation, and compares against model obtained from open-loop step pseudo-random binary sequence (PRBS) on the input.As expected, the model identified using closed-loop data predicts better.
Next, closed-loop simulation results for proposed controller and conventional MPC (i.e., MPC without Lyapunov constraint) with horizons 1 and 10 are presented in Figures 4-7.The controllers parameters are presented in Table 2.As can be seen, the LMPC has the best performance in stabilizing the system at the unstable equilibrium point.The MPC with a horizon of 1 is not capable of stabilizing the system, and the controller with a horizon of 10 reaches the set-point later compared to the LMPC.In addition, the evolution of the subspace states indicates better performance under the proposed LMPC.

Data-Driven EMPC Design and Illustration
Having illustrated the ability of the LMPC to achieve stabilization, it is next utilized to achieve economical objectives while ensuring stability.The Lyapunov based EMPC formulation is as follows: max ũk ,..., ũk+P where the value of ρ dictates the neighborhood that the process states are allowed to evolve within.c y and c u indicate output and input cost vectors.Other variables have the same definition as Equation (47).
Remark 8.In recent contributions [17,28], a Lyapunov-Based EMPC is proposed that utilizes data-driven methods to identify an empirical model for the system where the number of empirical model states is equal to the order of the plant model.In contrast, in the present work, the order of the model is selected based on the ability of the model to fit and predict dynamic behavior over a suitable range of operation, in turn allowing for an EMPC design that can reliably operate over a larger region.

Remark 9.
The EMPC formulation in the present manuscript utilizes a linear form of the cost function for the purpose of illustration.The proposed approach is not limited by this particular choice.Any other form of the cost function, including those where the costs could be time dependent, could be readily utilized within the proposed formulation.In such scenarios, the presence of the stability constraints provide the safeguards that allow the EMPC to move the process as needed to achieve economical goals.
Remark 10.The use of linear models in the control design opens up the possibility of utilizing MPC formulations [3,29] that enable stabilization from the entire null controllable region (the region from which stabilization is achievable subject to input constraints).The use of the NCR can, in turn, be utilized to maximize the region over which the EMPC can be implemented, thereby maximizing the potential economic benefit.Such an implementation, however, needs to account for potential plant model mismatch owing to the use of the linear model, and remains the subject of future work.
Next, the proposed Lyapunov-based EMPC (LEMPC) is implemented on the CSTR simulation example and compared to the LMPC implementation.The closed-loop results are presented in  Exploiting the flexibility of operation within a neighborhood of the origin, the LEMPC drives the system to a point on the border of that neighborhood, which happens to be the optimal operating point, instead of the nominal operating point.Figure 12 shows the comparison of the LEMPC and LMPC.As expected, the LEMPC achieves improved economic returns compared to the conventional MPC.Note that the LEMPC drives the system to a point within the acceptable neighborhood of the origin.

Conclusions
In this study, a novel data-driven MPC is developed that enables stabilization at nominally unstable equilibrium points.This LMPC is then utilized within an economic MPC formulation to yield a data driven EMPC formulation.The proposed approach is described and compared against a representative MPC, and shown to be able to provide improved closed-loop performance.

Figure 2 .
Figure 2. Data driven model validation results: measured outputs (dash-dotted line), state and output estimates using the (linear time invariant) LTI model model from closed-loop data and identification (dashed line), state and output estimates using the LTI model from open-loop data and identification (dotted line), observer stopping point (vertical dashed line).

Figure 4 .
Figure 4. Closed-loop profiles of the measured variables obtained from the proposed Lyapunov-based MPC (continuous line), MPC with horizon 1 (dash-dotted line), MPC with horizon 10 (dashed line), and MPC with horizon 1 and open-loop identification (narrow dash-dotted line) and set-point (dashed line).

Figure 5 . 4 Figure 6 .
Figure 5. Closed-loop profiles of the manipulated variables obtained from the proposed LMPC (continuous line), MPC with horizon 1 (dash-dotted line), MPC with horizon 1 and open-loop identification (narrow dash-dotted line) and MPC with horizon 10 (dashed line).

Figure 8 .
Figure 8. Closed-loop profiles of the measured variables obtained from the proposed Lyapunov-based economic MPC (continuous line) and the nominal equilibrium point (dashed line).

Figure 12 .
Figure 12.A comparison of the economic cost between the LEMPC (continuous line) and LMPC (dotted line).

Table 1 .
Variable and parameter description and values for the continuous stirred-tank reactor (CSTR) example.

Table 2 .
List of controllers parameters for the CSTR reactor.