Application of the OpenModelica-MATLAB Interface to Integrated Simulation and Successive Linearization Based Model Predictive Control

This paper presents the implementation of successive linearization based model predictive control (SLMPC) efforts through the interfacing of OpenModelica and MAT-LAB using the OMMatlab tool. The dynamic system (here a chemical process) and the model predictive control (MPC) algorithm are implemented in OpenModelica and MATLAB, respectively. The model linearization procedure is carried out through OMMatlab, which is highly optimized in terms of run-time by using a single executable file and adapting it at each sample time. Also, necessary theories for a continuous model discretization are discussed for both nonlinear Modelica and linearized continuous models. A procedure for constructing an Extended Kalman Filter (EKF) from a continuous Modelica model is also presented. The usability of the OpenModelica-MATLAB interface for SLMPC is demonstrated by con-trol of liquid levels in a tanks-in-series problem.


Introduction
In the recent years, the Modedica language has been widely used for modeling of systems described by differential algebraic equations (DAEs).The object-oriented nature of Modedica facilitates modular model construction and the use of powerful solvers as provided by commercial and open source Modedica-based simulators (e.g., Dymola, OpenModelica (Fritzson et al. 2020), Jmodelica).Some integrated features such as Optimica (Åkesson 2008) and CasAdi are also supplied for solving dynamic optimization problems.The reader is referred to (Ruge et al. 2014;Magnusson and Åkesson 2015) for a thorough discussion of these topics.Moreover, model predictive control (MPC) problems have been implemented in Open-Modelica; see e.g., (Bachmann et al. 2012) for MPC implementation of a batch reactor.
Conventional nonlinear model predictive control (NMPC) works by predicting the future behavior of a system and solving a dynamic optimization problem for minimizing a performance index, e.g., tracking error or operational cost; see, for example, (Ellis, J. Liu, and Christofides 2017;Heidarinejad, J. Liu, and Christofides 2013).Although using a rigorous nonlinear model increases the prediction accuracy, it leads to high computational expense (Zhakatayev et al. 2017) and potential convergence issues in the optimization routine.Successive linearization MPC (SLMPC) is an efficient alternative to NMPC (Seki, Ooyama, and Ogawa 2002;Cannon, Ng, and Kouvaritakis 2009;Cortinovis et al. 2014;C. Liu et al. 2015), especially in large-scale applications.In this method, the nonlinear model of the system is linearized successively at each sample time, and the prediction model is updated over time to preserve the prediction accuracy (Vrlić, Ritzberger, and Jakubek 2020).As a result of using a linear model, the dynamic optimization can be performed faster, favoring SLMPC over NMPC for large-scale systems in terms of computational expense, although the prediction accuracy may be undermined to some extent.
Many dynamic optimization or MPC efforts using Modelica models are reported in the literature.Franke (2002) employed Modelica to study optimal startup of a power plant, in which Dymola was used for generating S-Functions that would be imported into Simulink.Gräber et al. (2012) proposed a framework based on functional mockup interface generated from Modelica models for implementing NMPC on a vapor compression cycle.Also, L. Imsland, P. Kittilsen, and T. Schei (2010) integrated Dymola models into commercial software designed for NMPC and did a case study with an offshore oil and gas processing plant.Pandey et al. (2021a) employed OMJulia (Lie et al. 2019) for stochastic MPC of solar and hydropower plants described by a set of small-scale DAEs.Pandey et al. (2021b) implemented MPC for a power grid system, where a small-scale model was implemented in OpenModelica as the actual plant, and a separate model was developed in Julia for use in the control algorithm that included an unscented Kalman filter for the state estimation.OMJulia was used for interfacing OpenModelica with Julia.Also, when it comes to optimal control problems that require a linearized model, Jmodelica can offer a tool for linearizing the nonlinear model of the plant.Jmodelica supplies the interfacing capability with Python, enabling it to integrate Jmodelica simulation, optimization, and linearization features and design MPC problems that need linearization.It is possible to interact with the Modelica models and functional mockup units through IPython, a command shell for the programming environment (Andersson et al. 2018).The interested readers are referred to, e.g., ((Romero, Goldar, and Garone 2019;Pippia et al. 2021;Lars Imsland, Pål Kittilsen, and T. S. Schei 2009;Jorissen, Boydens, and Helsen 2019)) for other dynamic optimization studies involving Modelica models.
The use of engineering software such as MATLAB for Modelica models offers easy implementation for algorithm prototyping.However, a challenge in the application of MPC for Modelica models is the computational cost of interfacing the Modelica software with engineering software such as MATLAB.This is particularly true in case of SLMPC, where repeated model linearization through Jacobian calculations is required.The computational cost can grow quickly for large-scale models, making it prohibitively expensive to apply SLMPC or similar algorithms through interfacing of Modelica models with MATLAB or other engineering software.Therefore, the procedure needs to be carried out decently ensuring that the Jacobian calculations are undertaken in a time-efficient manner.
In this paper, efficient implementation of SLMPC for Modelica models in OpenModelica is addressed.In particular, the dynamic model of the system (serving as the virtual plant) is implemented in OpenModelica.Also, OMMatlab (OpenModelica 2021) is used to take control of the model simulation and make an interconnection between MATLAB and the OpenModelica model, which is transformed into an executable file.Moreover, the OM-Matlab linearization method is enhanced to operate considerably faster than its older version, making the new version more favorable for large-scale SLMPC.
These enhancements can also be implemented for other OpenModelica interfaces such as OMOctave, OMPython, and OMJulia so they can support efficient prototyping of advanced control algorithms requiring successive linearization.
The rest of the paper is structured as follows.In section 2, the general structure of the SLMPC problem is briefly introduced.In section 3, the method of assembling the discrete system model from its continuous form is presented.Section 4 discusses the details of the linear prediction model, the pertinent theories, and the enhancements made to OMMatlab to boost the linearization pro-cess.The state estimator design algorithm is discussed in section 5.In section 6, the dynamic optimization problem solved at each sample time is formulated.A case study of the SLMPC implementation is demonstrated in section 7. Finally, the paper is concluded in section 8.

Overview of Problem Blocks
As depicted in Figure 1, the problem under study includes an often nonlinear model that represents the actual physical system (here a chemical plant).Some of the plant variables are measured at every sample time.Since measuring all the variables is impractical, a state estimator is used to estimate the unmeasured state variables.These values are used to initialize the prediction model, which is a linearized form of the nonlinear model updated at each sample time.It is also possible to increase the frequency of prediction model updates by regenerating the linear model at each step over the prediction horizon instead of updating it only at the beginning of each sample time.The optimal control inputs are calculated over the prediction horizon by minimizing the objective function (e.g., the cumulative tracking error).Then, the first element of the optimal input trajectory is passed to the plant, and the procedure is repeated in the next sample time.The building blocks of the SLMPC are detailed in the following subsections.

System Model
The plant model is represented in a continuous, nonlinear, time-invariant state-space form as where x ∈ R n x is the set of n x system states, u ∈ R n u is the vector of n u control inputs, y ∈ R n y is the vector of n y system outputs, .., g n y ] are the sets of equations describing the state evolution and outputs of the system, respectively.OpenModelica automatically generates these functions from a high-level, object-oriented, equation-based model description, and can simulate the nonlinear model and generate an executable file that will be used for the state estimator and prediction model.The existing features of OMMatlab allow the user to take control of the executable file and overwrite the simulation setup and model parameters from MATLAB.Equation 1 is converted to a discrete state-space form with additive white Gaussian noise for process states and measurements.This is done by the methodology presented in (Brembeck 2019;Brembeck, Otter, and Zimmer 2011), the implementation of which is adapted to the MATLAB-OpenModelica environment.The discrete form reads x k = f k|k−1 + w k−1 (3) where w k and v k are additive process and measurement noises, respectively; and Q k and R k denote the process and measurement noise covariance matrices, respectively.The LHS of Equation 4 is computed by integration of the system over one sample time.In the continuous model implemented in OpenModelica, initial state values (defined in the initial equation section) are set parametrically so that the initial values can be overwritten and Equation 4 can be evaluated from MATLAB at each sample time.This is illustrated in Listing 1 for an arbitrary continuous Modelica model.The initial conditions x k−1 could be set from MATLAB by the setParameters() method of OMMatlab.Similarly, the system inputs u k−1 are specified by the setInputs() method.Finally, the integration in Equation 4 is performed over one sample time using built-in OpenModelica solvers such as DASSL and IDA.

Prediction Model
The system linearization required for the prediction model and its implementation are described in this section.

System linearization
Equation 1 and Equation 2 can be linearized around an arbitrary operating point (x op , u op ) by using the first-order Taylor expansion as follows. where Considering Equation 1 and Equation 2, the bias values are and getSolutions() methods.Equation 10and Equation 11 can then be expressed as where K xc = ( f (x op , u op ) − A c x op − B c u op ), and K yc = (g(x op , u op ) − C c x op − D c u op ).The continuous linearized model can be discretized for a given sample time T s as (see (Zhakatayev et al. 2017;Vrlić, Ritzberger, and Jakubek 2020) for a detailed proof).
in which and

Improving OMMatlab for successive linearization
Repeated linearization of a nonlinear dynamic model was computationally inefficient in the previous OMMatlab distributions.For example, a single invocation of the linearize() method for a system of DAEs with about 1100 equations took around 420 seconds on an Intel Core i5 7200U CPU.This would make implementation of SLMPC impractical as the linearization task should be performed at every sample time.The leading cause for this inefficiency was that the OpenModelica model had to be recompiled each time the linearize() command was called.The resulting linearized model had to be rebuilt to create an XML file so that the values of matrices could be extracted into MATLAB.
To solve this problem, OMMatlab is edited by the authors so that the initial executable file can be adapted and used in all invocations without being recompiled.Moreover, instead of using the time-consuming perturbation and finite differences for Jacobian approximation, the built-in automatic differentiation algorithm (invoked by the -generateSymbolicLinearization flag) is employed for exact linearization and construction of the model matrices (Braun, Ochel, and Bachmann 2011).The generated linearized model is populated in a .mfile, from which the matrix values are read.With these improvements, a single invocation of the linearize() command for the same DAE system now takes only a fraction of a second, which is a remarkable computational enhancement.This improvement can benefit any application requiring nonlinear Modelica model linearization through MATLAB, including SLMPC.

State Estimator
The extended Kalman filter (EKF) is used for state estimation.An EKF for Modelica models can be implemented in MATLAB based on the approach presented in (Brembeck 2019).The EKF uses linearization for state and measurement covariance propagation.The prediction and correction steps of the EKF for Equation 3 to Equation 5 proceed as follows.
• Prediction: • Correction: where the − and + superscripts respectively denote the predicted and corrected values.The EKF is initialized by setting x+ 0 and P + 0 to arbitrary or approximate values.It should be noted that the required Jacobian values F k−1 and G k can be calculated easily by a call to the linearize() method through OMMatlab.

Optimization Problem
MPC solves the following general optimization problem at the k-th sample time.
where J k is the objective function optimized over the prediction horizon H p .Note that N H p = H p T s .The decision variable set r is the trajectory of control inputs.Also, in order to reduce the computational cost, the control horizon is set as H c < H p , and r k+N Hc = r k+N Hc +1 = T s .The function h describes the discrete state evolution that is used for prediction.Also, M and S are the inequality and equality constraints, respectively, and X k is the estimated current state vector used to initialize the prediction model.
The MPC objective can be a tracking index, an economic index, or a combination of the two; see e.g., (Heidarinejad, J. Liu, and Christofides 2013; Ellis, J. Liu, and Christofides 2017).In this work, the dynamic optimization problem is solved using the sequential method (Chachuat 2009), where the control input trajectories are parameterized (usually in a piece-wise constant manner), and the dynamic model and the nonlinear optimization problem are solved in sequence as shown in Figure 2.This work uses the active-set optimization algorithm to solve the nonlinear optimization problem.

Case Study
A tracking SLMPC is implemented for a system of three interconnected cylindrical tanks as depicted in Figure 3.The flow rates of the three inlet streams are considered the control inputs.The outlet flow rates of the tanks are the measured outputs, and the tanks' liquid levels are regarded as the system states being estimated by the EKF.The tanks are empty at the initial time, and the control scenario is to move the tank levels to a specific set point by adjusting the inlet flow rates.The continuous dynamic model describing the system is as follows.
where A 1 , A 2 , and A 3 are the vessel cross-section areas, C v is the valve coefficient; q in1 , q in2 , and q in3 are the volumetric flow rates of the inlet streams, and q o1 , q o2 , and q o3 are the volumetric flow rates of the outlet streams.In case of reverse flow through the valves, the square root term in the valve equations becomes negative, leading to a complex number and solver failure.To avoid this situation and preserve local Lipschitz continuity, the valve equations are regularized as (Barton, Banga, and Galán 2000;Sahlodin 2022;Casella 1998) where ε > 0. Let A 1 = A 2 = A 3 = 1m 2 and C v = 0.5.The model nonlinearity comes from the outlet stream equations.Also, the desired set points for the tank levels are h sp = [0.56,0.52, 0.36] T .For the dynamic optimization given in Equation 29, the following quadratic tracking objective function with control move penalization is defined.
The weighting factors W 1 and W 2 are positive-definite matrices set as Also, the following path constraints are added to avoid reverse flow between the tanks.
where δ > 0 is a small regularization parameter (Chachuat 2009).The system is simulated for a time span of 50 minutes with sampling time of 6 seconds.The prediction and the control horizons are set to H p = 1.5 and H c = 1 min, respectively.The control inputs are bounded as 0 ≤ q in1 , q in2 , q in3 ≤ 0.3.Note that a discretized linear model is applied as the prediction model is updated at each sample time.
The results of the SLMPC are presented in the sequel.Figure 4 shows trajectories of the measured outlet flow rates that are used to estimate the tank levels.
Figure 5 depicts the actual and estimated trajectories of the tank levels.The initial guess for the state estimator is ĥ+ 0 = [0.1,0.1, 0.1] T , which is different than the actual values h 0 = [0, 0, 0] T .Despite this difference, the EKF is able to track the actual state trajectories successfully.It is also seen that the actual tank levels approach the set points shortly after the SLMPC is executed.The optimal October 9-11, 2023, Aachen, Germany Figure 3. Schematic of the system with the SLMPC and estimator blocks.control inputs are plotted in Figure 6.It is observed that the control inputs are slightly oscillatory as a result of the process noise.
It is worth noting that the total simulation takes only 25 minutes by the enhanced OMMatlab, while it would take around 7 hours when employing the previous OMMatlab versions.This proves that the linearization part is the main bottleneck in the SLMPC procedure.
The code for the case-study presented in this manuscript can be found on GitHub at https:// github.com/pseAUT/SLMPC_OMMatlab.

Conclusion
The implementation of the SLMPC algorithm using Open-Modelica and MATLAB has been demonstrated in this paper.The OMMatlab API is upgraded to avoid repeated compilation of the OpenModelica model into an executable file at each sample time.In this way, the system can be linearized efficiently and the model matrices can be obtained directly in a .mfile, thanks to the existing OpenModelica flags that make it possible to generate the linearized model in different formats.Therefore, the successive linearization runtime is optimized consid- ) Therefore, it is possible to construct continuous linearized models by having A c , B c , C c , D c , and the bias values.OMMatlab can provide all these matrices along with ( ẋ, y) at any specific simulation time by the linearize() Session 3-C: Applications of Modelica for optimization and optimal control
Application of the OpenModelica-Matlab Interface to Integrated Simulation and Successive Linearization Listing 1. Continuous system model in Modelica Application of the OpenModelica-Matlab Interface to Integrated Simulation and Successive Linearization Based Model Predictive Control