PARODIS: One MPC framework to control them all. Almost.

—We introduce the M ATLAB framework PARODIS, the Pareto optimal Model Predictive Control framework for distributed Systems. It is a general-purpose, ﬂ exible and easy-to-use framework for discrete state space models. Special features are the support of distributed (hierarchical) systems, scenario-based optimization and built-in methods for determination of the Pareto front and selection of a solution. It uses the popular M ATLAB framework YALMIP for the symbolic formulation of optimization problems and models.


I. INTRODUCTION
Model Predictive Control (MPC) has seen rapid developments over the last decade, both in theory and applications. To validate new MPC algorithms in extensive simulation studies, researchers face the challenge of implementing both control algorithms and system models. Therefore, they need a framework which is i) general enough to cover their problem class and common control approaches, ii) flexible enough to also allow implementation of new specialized algorithms and iii) easy enough in use to maintain a reasonable time investment.
However, current MPC frameworks have a different use case, are very specific or lack important features to satisfy these requirements. Multiple frameworks have been developed for real-time control, e. g. ACADO [1] or GRAMPC [2]. Usually based on C/C++ and designed to run on embedded hardware, they are very strong in controlling highly dynamical nonlinear systems, even with sampling times in the (sub)millisecond range. While they are very apt to perform this task, they do not generalize well and are not suitable for prototyping more complicated algorithms such as economic or distributed MPC.
Other frameworks are specialized in a specific field of application, e. g. robotics [3] or, commonly, building control [4], [5]. In this case, they usually employ more sophisticated MPC algorithms and might also feature co-simulation with simulation tools such as Modelica or EnergyPlus. However, their use is limited to the respective application. They do not support a more general form of system classes and are not flexible enough for implementation of new control algorithms. The official MATLAB MPC toolbox [6] is powerful and can be used for more than just tracking MPC, but the implementation complexity rises significantly with deviation from this standard case, e. g. distributed systems. With the popular MATLAB toolbox MPT3 [7], regulation or tracking MPC for linear, (piecewise) affine and mixed logical dynamic systems can be set up extremely easily. However, despite some impressive features like automatic conversion to explicit MPC, it does not support economic MPC in general and is not designed for implementation of more complex algorithms, i. e. not flexible. MPCTools [8] is an Octave interface to the popular optimization framework CasADi [9] with means to simplify the definition of nonlinear MPC systems. While there are no significant explicit limitations due to the use of CasADi, it lacks general features such as direct support for parametrization with external data or distributed systems. Another framework that makes use of CasADi is do-mpc [10], based on Python. It is more capable than MPCTools, as it allows for full parameterization of the system dynamics. It focuses on nonlinear time-continuous systems, but does support discrete system formulations. It has a very modular system structure and offers various features. Overall, flexibility and ease of use are given. However, there are no pre-defined interfaces for distributed systems and usage of large data sources.
Here, we present an alternative MATLAB framework, called PARODIS -Pareto optimal MPC framework for distributed systems. It provides generality to a great extent while maintaining flexibility and easy usage and distinguishes it from other MPC frameworks by i) supporting distributed (hierarchical) systems, ii) providing convenient interfaces to use large-scale datasets for predicted and actual disturbances, iii) respecting scenarios for uncertain disturbances and parameters in a customizable fashion and iv) as its main innovation, the integrated support for multiple objectives by automated generation of the Pareto front and selection of solutions. Note that the Pareto optimization functionality is an optional feature. The complete functionality and its concept, showing its wide generality while maintaining flexibility, is explored in detail in Section II. The basic usage is explained in Section III and an exemplary case study is given in Section IV. We end with a conclusion in Section V.

A. Functionality
PARODIS allows fully parameterized (and thus time-variant) descriptions of discrete state space systems in the form  where x, u and d are the state, input and disturbance vectors, respectively. Multiple systems, in PARODIS called agents, can be defined and controlled in a distributed (or hierarchical) fashion. Model predictive controllers are defined for every agent by default, whereby the cost functions can be defined arbitrarily, i. e. as in economic MPC. However, controllers could easily be exchanged. There are two specialties of PARODIS. First, it distinguishes not only between predicted and real disturbances d, but it also supports the consideration of multiple scenarios in the prediction of disturbances and parameters, thus advocating it for the use of scenario-based stochastic MPC. As sources, either function handles or csvfiles can be defined, thus allowing for easy import of large data sets for long-term simulations. In contrast to do-mpc [10], the scenarios can be chosen freely (not only as a full scenario-tree of all possible combinations of distinct values). Furthermore, the incorporation of scenarios in the optimization problem is customizable. Second, it ships with algorithms for multiobjective optimization, i. e. both methods for determining the Pareto front for multiple objectives and methods for automatic selection of a solution, which is repeated at every time step. PARODIS utilizes YALMIP [11] and its symbolic variables to formulate the optimization problem, which makes the large selection of solvers supported by YALMIP accessible. 1

B. Concept
PARODIS is object-oriented. All necessary information is wrapped in an instance of the Simulation-class, which may contain multiple instances of the Agent-and TimeSeries-classes (figures), as illustrated in Figure 1.
A simplified flow chart diagram is shown in Figure 2. In 1 See https://yalmip.github.io/allsolvers/ for an overview. distributed or hierarchical MPC approaches, agents usually interact with each other and negotiate what actions they should take in the current time step. In PARODIS, this interaction is called the negotiation and is separate from the agents' actual step. There are two options for this negotiation: The first one is for the agents to execute their negotiation step in a fixed order. Here, their optimal control problem is solved as usual and a user-defined callback is called afterwards. Using this callback, agents can exchange information freely. The second option is a user-defined interaction using simulation.negotiationHandle: This is a user-defined function, where any arbitrary algorithm can be implemented. In this function, all agents, controllers and simulation data can be accessed. After the negotiation, the actual simulation step is executed.

Simulation
Step Agent's Sampling Time Note that multiple optional callback-functions exist and enable the user to adjust the standard simulation flow, e. g. agent.callbackMeasureState() to update an agent's state before optimization.
PARODIS supports different agent sampling times, as long as the first steps of the agents' horizons are multiples of each other. Thus, an execution order is determined at the beginning. Figure 3 illustrates this. Within the iteration of a simulation step, every agent's doStep()-function is called following this execution order. Thereby, all disturbances and parameters are set. Then, the optimal trajectory is determined, plots are updated and the first input is applied.
In PARODIS, the Controller-class handles the formulation and solving of the optimization problem. There are currently three different controller classes available in PARODIS, the SymbolicController, the ExplicitController and the ParetoController. All controllers are easily interchangeable, as they all provide the same interface.
Both the symbolic and explicit controller formulate a standard optimization problem, where the objective function is the (weighted) sum of the configured cost functions. The difference between the two controllers is the representation of the optimization problem: The symbolic controller pre-compiles a fully parameterized optimization problem using YALMIP's optimizer feature, while the explicit controller re-builds the problem at each time step and replaces all parameters with explicit values. The SymbolicController is to be considered the default controller, as it is very efficient for smaller problems, while the ExplicitController is suited for larger problems and debugging purposes as it scales much better and is more transparent in its model representation.

C. Pareto Optimization
Pareto optimization refers to the optimization of multiple (competing) objectives. Usually, the set of Pareto-optimal (or non-dominated) solutions is derived. A solution is nondominated, if there exists no other solution which is better in objective i without being worse in at least one other objective j [12]. This set is also called Pareto front, from which a compromise is then selected. The ParetoController provides such functionality. In the same way as any other controller, it derives an optimal input to be applied to the system. It does so in three steps. First, an extreme point function calculates the extreme points, i. e. the Pareto optimal points minimizing single objectives. Second, a front determination scheme systematically samples points on the Pareto front, using YALMIP optimizer objects for increased efficiency. Finally, a Pareto optimal solution with its corresponding trajectories is selected. While metric functions select a solution automatically, the implemented interactivity tool enables the user to instead select a point manually (for n obj = 2 or 3 objectives). For each category, PARODIS provides a set of predefined functions, which can be selected within the ParetoController's configuration. Due to the modular implementation, users can alternatively define their own methods. Table I summarizes the functions currently implemented.

III. USAGE
In this section, we give a brief explanation of how problems are implemented in PARODIS. For a detailed explanation, the user is referred to the documentation.

A. Installation
Installing and setting up PARODIS is very straightforward. Simply download the latest release from the GitHub repository, 2 unpack the archive and add the directory to the MATLAB path. As PARODIS is built on YALMIP, it has to be installed separately in the same way. PARODIS itself does not depend on any MATLAB toolboxes.

B. Creating Agents
The agent is the core of PARODIS. It represents a system and fully determines its optimal control problem. To create an agent, a model and a controller need to be defined first. For the model, the user has to implement a function which receives the sampling time T s and returns the corresponding ordinary difference equations, the number of states n x , inputs n u , and disturbances n d . For the controller, one of the three available types presented in Section II-B can be used, e. g.  Lexicographic Approach Optimization in lexicographic order, i. e. minimizing a single objective J i with optimal values for J j =i as constraints [12]. Miniscule Weight Approx.
Approximation of the extreme points with weights w i = 1 and w j =i = 10 −5 .

Normalized Miniscule Weight Approximation
Miniscule Weight Approximation with prior normalization of the objectives to a range of 0 to 1 for a better conditioning of the optimization problem.

Front Determination Scheme
Normal Boundary Intersection The boundary plane, i. e. the hyperplane connecting all extreme points, is evenly sampled. These plane points are then projected onto the Pareto front in the boundary planes' normal vector direction [12]. Focus Point Boundary Intersection Self-developed. Similar scalarization method like NBI, but with different plane sampling and projection direction.

Adaptive Weight Determination Scheme
Geometric interpretation of the weighted sum method. From n obj parents, weights are determined for which the new solution lies between the parents [13].

Metric Function
Closest-to-Utopia-point Herein lies one of the strengths of PARODIS: switching to Pareto optimization is as easy as swapping the controller.  To each subplot, either line plots or heatmap plots can be added, which may display information from different agents. The first argument defines from which agent data should be displayed. Using var, the user can plot any states, inputs, disturbances, eval-or cost functions. All plots can be configured to be plotted live during the simulation and/or at the end. This can be defined either individually for each figure, or globally for all figures.

D. Running a Simulation
With agents and plots defined, creating and running a simulation is very straightforward. Note that agents and figures have to be added to the simulation instance.

IV. EXEMPLARY HIERARCHICAL SYSTEM
To illustrate the most important features of PARODIS, we present the model of an energy management system of an office building, which is controlled using a hierarchical MPC approach.

A. Higher Level
The main model, in the following called the higher level (HL), contains the energy E of a stationary battery and the building's overall temperature ϑ b as states. The system's inputs are the electrical power obtained from the power grid P grid , the electrical power P chp from the building's combined heat and power plant (CHP), the thermal power from a gas heating unitQ rad and the thermal cooling power from a heating, ventilation and air conditioning systemQ cool . As uncontrollable disturbances acting on the system, the electrical power from a photo-voltaic plant P ren , the electrical power demand P dem and the air temperature ϑ air are considered. With these, the system's dynamic is given by Hair . T s is the sampling time, C th,b the building's total thermal capacity, H air,b the total heat transfer coefficient, ε c the energy efficiency ratio of the cooling machine and c CHP the CHP's current constant. For the higher level, Pareto optimization is applied with the monetary costs and the quadratic temperature deviation from 21°C as two objectives. For a detailed insight into both modeling and the Pareto approach, the reader is referred to [14], [15].

B. Lower Level
As a second (sub)model, in the following called the lower level (LL), the building's thermal energy system is split up into 9 different temperature zones (neglecting the electrical energy system). In the time-continuous case, the temperature ϑ b,i of zone i is described bẏ where C th,i is the thermal capacity from zone i, H air,i is the heat transfer coeffienct between zone i and the outside air, β ij is the heat transfer by coefficient between zones i and j, andQ heat,i andQ cool,i are the respective heating and cooling powers allocated to zone i. For brevity, we omit stating the (discretized) state space model used for implementation. Note that 9 i=1 C th,i = C th,b and 9 i=1 H air,i = H air,b . The lower level's task is to distribute the heating and cooling power produced on the higher level between the 9 temperature zones. Thus, its only objective (and cost function) is to minimize the temperature deviation of all zones from the desired 21°C, Note that the deviation for every zone is scaled by its capacity to prevent the preferred regulation of smaller zones. w_th is a weighting vector with elements w th,i = C th,i C th,b . Afterwards, the step is performed as usual.

C. Communication
The lower level has to use the total produced heating and cooling powers as constraints for its own input variables, i. e.
This is realized by adding 2 parameterized constraints to the controller in its definition before simulation:

D. Results
For both systems a prediction horizon of 24 h is used. The first 8 h are split in 15 min steps, the second 8 h in 30 min steps and the last 8 h in 1 h steps, i. e. a total of N pred = 56 steps. On the higher level, the resulting optimization problem (including, e. g., epigraph reformulations [14]) consists of 338 decision variables and 1066 constraints. The lower level results in 1521 decision variables and 3649 constraints. With GUROBI as the solver, the simulation of 2 days with realistic data from April 2020 (and Pareto optimization on the higher level) takes 247 s on an Intel i7-8550U notebook CPU. Without Pareto optimization, it reduces to 58 s. Figure 4 shows the resulting trajectories for the higher level, Figure 5 the Pareto fronts for its two objectives (monetary and comfort costs). Figure 6 shows the temperature and corresponding heating power progressions of the lower level agent's zones.

V. CONCLUSION
We present the MATLAB MPC framework PARODIS of which we hope that it will support researchers in both implementation and testing of newly developed MPC algorithms as well as empirical verification of applications by data-based simulation studies. It is the first of its kind with built-in Pareto optimization methods. Even if permanent multi-objective optimization is not intended, they can be used as an automatic setup for finding appropriate weights for multiple objectives. The modular object-oriented structure of PARODIS allows for easy adaption for both new models and control algorithms. Thus, PARODIS is not limited to its present scope.