Design, Implementation and Simulation of an MPC Algorithm for Switched Nonlinear Systems under Combinatorial Constraints

Within this work, we present a warm-started algorithm for Model Predictive Control (MPC) of switched nonlinear systems under combinatorial constraints based on Combi-natorial Integral Approximation (CIA). To facilitate high-speed solutions, we introduce a preprocessing step for complexity reduction of CIA problems, and include this approach within a new toolbox for solution of CIA problems with special focus on MPC. The proposed algorithm is implemented and utilized within an MPC simulation study for a solar thermal climate system with nonlinear system behavior and uncertain operation conditions. The results are analyzed in terms of solution quality, constraint satisfaction and runtime of the solution steps, showing the applicability of the proposed algorithm and implementations


Motivation
The term hybrid systems refers to dynamic systems whose development is influenced by continuous as well as discrete variables (Mendes et al., 2017). A particular class of hybrid systems entitled switched systems refers to systems that are described by continuous states but which are subject to state-dependent discontinuities (Bock et al., 2018). 1 Such discontinuities can either be invoked internally depending on current system states, or by external influences which can be considered as controls of a system. Within this work, we will focus on the latter case. 2 Switched systems can be described by a finite number of subsystems where the current activation status of a subsystem depends on the current system states and controls (Kreiss et al., 2017). They can be found, e. g., in the process industry, the energy sector and in Heating, Ventilation and Air Conditioning (HVAC) of buildings, where systems often depend on both continuous controls (e. g. mass flow rates, valve positions) and discrete controls (e. g. operation status and modes of machinery). According systems and their components often show nonlinear behavior and are constrained to certain operation conditions, while additional constraints such as minimum operation time of machinery and start up costs have to be considered.
For control of systems that are generally subject to constraints and/or strongly depend on external factors (e. g. renewable energy systems), MPC yields great potential in terms of constraint-satisfactory and efficient system operation. However, application of MPC techniques leads to optimization problems that need to be solved on a time scale suitable for real-time control of these systems.

Relevant literature
If the number and sequence of switches on the considered control horizon is determined a priori for a system, switching time optimization can be applied in order to determine the optimal time points when the subsystems of the considered system needs to be exchanged in order to follow a certain objective, without needing to actively consider discrete controls within the optimization problem. Applications and variants of this approach can be found, e. g., with Gerdts (2006), Heidarinejad et al. (2013), Bau et al. (2015) and Farshidian et al. (2017). However, application of this approach is problematic if the sequence and number of switches are not known a priori but have to be determined, cf. Sager (2009).
If active consideration of discrete controls cannot be avoided within the problem formulation, this leads to Mixed-Integer Optimal Control Problems (MIOCPs) that need to be solved on a suitable time scale for application within MPC of the system. For solving MIOCPs and Optimal Control Problems (OCPs) in general, the use of direct methods, especially direct multiple shooting (Bock and Plitt, 1984) and direct collocation (Tsang et al., 1975), is favorable (Sager, 2009;Binder et al., 2001). For MIOCPs, this results either in Mixed-Integer Linear Programs (MILPs) or in Mixed-Integer Non-Linear Programs (MINLPs), depending, i. a., on the applied method, whether the considered dynamic system is linear(ized) or not, and on the type of constraints present within the problem formulation.
In the literature, several approaches for solution of mixed-integer programs within MPC applications are presented. For some applications, direct use of general solvers can be sufficient, as shown in Kobayashi and Imura (2012) and Khakimova et al. (2017) where corresponding MILPs are solvable within suitable time scales using CPLEX (IBM Corp., 2018) and Gurobi (Gurobi Optimization, Inc., 2016), respectively. Bemporad and Morari (1999) present an approach for modeling and predictive control of hybrid systems referred to as Mixed Logical Dynamical (MLD) systems (linear systems with both continuous and discrete inputs and states and subject to linear constraints) and solve the resulting MIOCPs using an available, general Mixed-Integer Quadratic Program (MIQP) solver.
Today, tailored algorithms that aim to exploit special aspects of MIOCPs for MLD systems have been presented, e. g., by Frick et al. (2015). This is exemplary for several approaches that can be found in the literature that aim towards decreasing solution times of mixed integer problems for real-time applications by tailored solution methods. Mayer et al. (2016) present a tailored Branch-and-Bound (BnB) method incorporating a search space reduction to save computation time. Luchini et al. (2017) solve a mixed-integer MPC problem by application of a linear model-predictive controller with a longer prediction horizon in combination with a mixed-integer MPC problem with shorter horizon.
Compared to MILP, less work can yet be found regarding the utilization of MINLPs within real-time MPC. Since MINLPs are comparatively hard to solve, nonlinearities are in fact often avoided, e. g., by application of piecewise linearization for nonlinear functions, cf. Rawlings et al. (2018). However, such approximations might not always be applicable or sufficient. Apart from that, approaches can be found that aim towards faster solution of MINLPs by decomposing an original problem into a sequence of super-and subordinate optimization problems along different layers of abstraction of a problem, cf. Dias et al. (2018), Schweiger et al. (2017) and Manaf et al. (2017). Possibilities for such separations are problem-specific and application can be problematic if abstraction layers cannot be sufficiently separated, e. g., if the allowance of operation of scheduled machinery depends on the actual system states that are subject to subordinate control levels.
A more general decomposition approach for fast solution of MINLP problems is presented by Sager (2009) where it is proposed to first solve a relaxed version of an MINLP, which is an Non-Linear Program (NLP), and obtain an integer trajectory for the discrete controls by application of a special rounding scheme entitled Sum-Up-Rounding (SUR). Though this might result in a suboptimal solution for the original problem, the error introduced by this approach is bounded and the method allows for solution of MINLPs within a time-scale that renders those real-time usable. Successful applications of this approach can be found in Kirches (2011) and Ebrahim et al. (2018).
In Sager et al. (2011) and more generalized in Zeile et al. (2018), an enhanced method for generating binary trajectories 3 named Combinatorial Integral Approximation (CIA) is proposed which can directly consider constraints such as SOS1-, maximum switching and dwell time constraints, but is more computationally demanding than SUR. A CIA problem is an MILP that can either be solved using a general MILP solver or by application of a fast, tailored BnB method. Additionally, Jung (2013) proposes to solve the MINLP once more subsequent to the binary approximation and with the obtained binary controls fixed, which again is an NLP, in order to adjust continuous controls and states to the obtained binary trajectories. This procedure yields the possibility to solve MIOCPs for nonlinear systems under quite arbitrary constraints, such as soft constraints for state boundaries, constraints on maximum switching and dwell times for discrete controls, as well as for specification of machinery operation conditions.
In Bürger et al. (2018), the applicability of the decomposition method described above for optimization-based control of a complex renewable energy system in the form of a Solar Thermal Climate System (STCS) is investigated by solving individual MIOCPs for a nonlinear model of the system which incorporates a single binary control switch and combinatorial constraints. It is shown that the solution time is tremendously decreased in comparison to using the general MINLP solver Bonmin (Bonami et al., 2008) while comparable solution quality is achieved.

Contribution of this work
In this work, we build upon Bürger et al. (2018) and utilize the MINLP decomposition solution method within a warm-started algorithm for MPC of switched nonlinear systems under combinatorial constraints.
For this, we first present a new Python toolbox for solving CIA problems with special focus on application within MPC to facilitate easier use of the presented methods. Since the runtime of the CIA method rises with the number of binary variables treated, we introduce preprocessing heuristics for complexity reduction of CIA problems and show that those can decrease the solution time to a high extent for an exemplary benchmark problem.
Afterwards, we present a software implementation of the entire MPC algorithm and apply it within a simulated MPC loop for an extended version of the STCS model incorporating an increased number of states and continuous controls as well as two mutually exclusive binary switches under dwell-time-and maximum-switching-constraints and whose operation conditions are subject to uncertainty. The performance of the algorithm and its implementation is assessed by an analysis of the obtained results in terms of solution quality, constraint satisfaction and runtime of the particular solution steps.

Notation
In the following, we denote by a i the i-th element of a vector a. Further, we denote by A i the i-th column of a matrix A and by A (i,...,j) a submatrix of A that consists of all columns from i to j of this matrix. Indexing of vector and matrix elements starts at 0.

Used hardware and software
All computations are conducted on a Fujitsu P920 Desktop PC with an Intel Core i5-4570 3.20 GHz CPU and 32 GB RAM running Xubuntu 16.04 and with Python 2.7, gcc 5.4, CasADi 3.4.4, Ipopt 3.12.3. and SCIP 5.0.1.

Solution of MIOCPs for MPC applications
Within this section, we describe the decomposition approach for solution of MIOCPs utilized within this work as well as the CIA procedure in more detail. Further, we highlight selected aspects of problem formulation for MIOCPs within MPC.

Decomposition approach for approximate solution of MIOCPs
Let us consider a generic MIOCP with differential states x(t) ∈ R nx , continuous controls u(t) ∈ R nu , binary controls b(t) ∈ {0, 1} n b , slack variables s(t) ∈ R ns and time varying parameters c(t) ∈ R nc with n x , n u , n b , n s , n c ∈ N on a given time horizon The objective of the (MIOCP) is to minimize the sum of a continuous-time Lagrangian cost functional L(·) and a Mayer term M (·). The Lagrangian term measures the "running" cost during the time horizon for given differential states, control and slack variables and is assumed to be integrable on the whole time horizon. The Mayer term, on the other hand, offers the option of including a desired differential state cost at the end of the time horizon in the objective. Note that we assume in (1b) a nonlinear system of explicit Ordinary Differential Equations (ODEs) 4 , where the binary controls appear affinely so that we can write the right hand side as a sum of functions f 0 (·) and f i (·), i = 1, . . . , n b . The function r(·) contains path constraints with possibly time-dependent lower bound r lb (t) and upper bound r ub (t). Within r(·), different types of constraints, such as specifications of operation conditions for machinery, maximum switching constraints, dwell time constraints etc., can be specified. If desired, path constraints can also be formulated as soft constraints using elements of s. We assume that the functions f 0 (·), f i (·), i = 1, . . . , n b and r(·) are continuously differentiable in all arguments within the domain of interest. By X and U we denote the feasible domains of the states and continuous controls, respectively.
For solution of OCPs, the use of direct methods (first discretize, then optimize approach) is favorable, especially direct multiple shooting (Bock and Plitt, 1984) and direct collocation

Algorithm 1: Decomposition algorithm for solution of MIOCPs
Input : Discretized (MIOCP) instance with time grid T N , initial guesses for x, u, b. Output: (Local) Optimal variables x * , u * , b * , s * with objective L * = L(x * , u * , b * , s * ). 1 Initialize s = 0 and solve NLP rel → x, u, b rel , s, L rel . (Tsang et al., 1975), which for MIOCPs results in MINLPs. For solving MINLPs within an MPC setting, we apply the decomposition approach proposed by Sager (2009) and solve first the relaxed MINLP with dropped binary control constraint, which is an NLP and therefore called NLP rel . Next, we approximate the obtained relaxed binary controls b rel ∈ [0, 1] n b with binary controls b bin ∈ {0, 1} n b using CIA, which as described by Sager et al. (2011) is an MILP. The final step is to solve the MINLP again with binary controls fixed, which again is an NLP and referred to as NLP bin , in order to adjust the states x, continuous controls u and slack variables s to the obtained binary solution, cf. Jung (2013). Algorithm 1 summarizes the steps of the decomposition approach. Note that this approach does not solve the exact MINLP but an approximation given by a sequence of subproblems of which each one is less hard to solve than the original MINLP. Nevertheless, the approach is favorable for application here due to its real-time capabilities which we prefer due to our focus on MPC, while asymptotic objective error bounds for the approach depending on the grid step size have been proven (Sager et al., 2011). In the upcoming section, the CIA method is described in more detail.

Combinatorial integral approximation under combinatorial constraints
The idea of CIA is to approximate the relaxed binary controls b rel ∈ [0, 1] n b with binary variables b bin ∈ {0, 1} n b in the sense that the maximum integrated difference of b rel and b bin for each discretization grid point is minimized. This distance is commonly measured with the maximum norm, cf. Jung (2013), and can be reformulated into an MILP which provides the opportunity to intuitively include combinatorial constraints. In the following, we are going to formally introduce this CIA problem.
Let the ordered set G N = {t 0 < . . . < t N = t f } denote a time grid with ∆t i = t i+1 − t i used for discretization of the controls in (MIOCP) meaning that the controls can only change values on these grid points i = 0, . . . , N − 1. Hence, the relaxed binary solution resulting from solving NLP rel is in matrix form b rel ∈ [0, 1] n b ×N . For formulating the dwell time constraints with minimum dwell times L k ∈ R + , k = 1, . . . , n b , we introduce the grid index subsets G i,L k . By further introducing slack variables η ∈ R, s σmax ∈ [−1, 1] n b ×N −1 , we define (CIA) in the presence of constraints on the maximum amount of switching actions per control σ k,max , k = 1, . . . , n b on the given time horizon, and additional dwell time constraints, as follows: Note that we minimize the accumulated difference between relaxed and binary control variables, as illustrated later in, e.g., Figure 3, by adding constraint (2c). Since we are interested in the absolute value of this sum, we add the constraint both with plus and minus signs, indicated by "±". The inequality constraints (2d) to (2h) represent the combinatorial constraints described above, while (2i) is the so-called Special Order Set 1 (SOS1) condition, which causes exactly one mode or the off mode to be active at any time.
Overall, the solution of this problem yields a binary approximation b bin ∈ {0, 1} n b ×N for b rel with an approximation error that is bounded depending on the discretization grid. Further details on CIA can be found in Sager et al. (2011) andJung (2013).

Time transformation of ODEs
If a time grid G N used within a moving-horizon MPC application consists of nonequidistant time points, the duration of a certain grid interval ∆t i can change in between two MPC steps. To avoid regeneration of the problem discretization and the according derivatives needed by an NLP solver, a time transformation to the ODE right hand side can be applied, so that at time point t i , instead of is used for formulation of the problem discretization, while the corresponding values of the time points t i and t i+1 enter NLP rel and NLP bin as parameters. According formulations are also used within switching time optimization, where the durations of time steps enter an OCP as optimization variables, cf. Gerdts (2006).

Formulation of combinatorial constraints as smoothened vanishing constraints
If the allowance for activation of binary switches b within a mixed-integer optimization problem depends on conditions h(·) as in which is a particular case of the path constraints (1c), there exists several ways to formulate such conditions, while the choice of formulation can have different effects on the optimization problem. An overview and comparison of different formulation techniques in the context of this work is given by Jung et al. (2013). According to the results given there, we will in the following use smoothened vanishing constraints for formulation of combinatorial constraints to prevent, e. g., violation of Linear Independence Constraint Qualification (LICQ) conditions, cf. Jung et al. (2013).

A new software package for solution of CIA problems
In the following, we introduce a new software package for solution of CIA problems including preprocessing heuristics for complexity reduction of this problem class.

Introduction of the software package pycombina
According to Sager et al. (2011), (CIA) can either be solved using a general MILP solver such as Gurobi or SCIP (Gleixner et al., 2017), or a BnB algorithm tailored to the structure of (CIA) can be applied which was found to be more efficient in terms of computation time. Jung (2013) presents several variants of the BnB algorithm further referred to as (BnB).
To simplify the use of methods and algorithms presented within this work, we introduce the open-source software package pycombina for solution of (CIA). With pycombina, we provide a Python module for solving (CIA) under consideration of combinatorial constraints. pycombina is licensed under the GNU Lesser General Public License (LGPL) and contains a fast implementation of (BnB) presented in Sager et al. (2011) and Jung (2013) using a bestfirst search strategy, i. e., the first solution found by the method is also a global optimum of (CIA).
While major parts of pycombina such as the problem setup and preprocessing heuristics (see Section 3.2) are written in Python to favor flexible applicability and extensibility, computationally heavy parts of (BnB) are implemented in C++ and interfaced using pybind11 (Wenzel et al., 2016) to increase computational performance. Maximum-switching-, SOS1-and dwell-time-constraints are, if present in a problem, actively considered within the implementation of (BnB) to allow for taking off bigger steps and shortcuts, cf. Jung (2013). In addition to (BnB), (CIA) can be solved using Gurobi and SCIP, which is mainly intended for benchmarking purposes. 5

A complexity reduction heuristic for CIA problems
Since the effort to solve (CIA) can increase strongly with the number of variables considered for the problem, we propose preprocessing heuristics for complexity reduction of CIA problems. The main part of the procedure is based on the idea that, if already in the relaxed solution a binary variable is constantly 0 or 1 (or sufficiently close to one of these values) over several time steps, then this is likely to persist for the corresponding time points in the binary solution. Therefore, the CIA solution method does not need to process every time point in this case, but can handle several time points of the grid within bigger time steps.
Since it is likely that many relaxed binary values will already be close to either 0 or 1 due to the typical bang-bang-type solution of OCPs (Sager, 2009), grouping of such sequences can be expected to reduce the number of time points that need to be processed to a high extent.
A detailed description of the procedure is given in Algorithm 2. Iterating over the relaxed binary solution given in matrix form, submatrices of a column span δ are checked whether all elements are sufficiently close to either 0 or 1 depending on a threshold ε and whether all columns of such a submatrix are identical. If this is the case, the center column of this submatrix is marked not to be considered within the reduced problem. The smaller the column span parameter δ is chosen, the more columns will potentially be removed, but also the more degrees of freedom might be taken away from the CIA solution method.
In addition to that, further measures have been identified and included in pycombina that can reduce the effort for solving (CIA) to a greater extent: 1. binaries whose relaxed solution is constantly 0 over the whole time horizon are directly set to 0 and not considered within the problem formulation; 2. in case a control is constantly 1 over the whole horizon, the solution are directly given since all other controls need to be equal to 0 due to the SOS1 constraints; 3. in case the required number of switches indicated by the relaxed solution for a binary control is less than the number of maximum switches defined by the user, it can be automatically reduced to this number.

Performance benchmark for pycombina
Within a brief benchmark, we compare the performance of the different solution methods incorporated in pycombina and the impact of the preprocessing heuristics on the solution time and quality of a problem.
Algorithm 2: Complexity reduction heuristic for CIA problems

Benchmarking setup
The input data used within this benchmark is the relaxed solution of the binary variables from the Lotka Volterra Multimode fishing problem taken from the MIOCP benchmark library mintoc.de (Sager, 2012) as shown in Figure 1, which contains n b = 3 binary control inputs that depict different fishing options. To investigate runtime development for an increasing number of optimization variables, we solve a CIA problem for data sets representing the trajectories in Figure 1 by N 1 = 60, N 2 = 120, N 3 = 240 and N 4 = 480 data points on an equidistant time grid from t 0 = 0 to t f = 12, once using (BnB) and once using SCIP. For each solver, we solve the problem once in original size and once after application of the preprocessing heuristics. For this study, we determine a maximum of σ max = 2 switching actions per control input and that at most one control can be active at a time.   Figure 2 shows a comparison of the solution times obtained by application of the BnB algorithm and SCIP as well as the impact of the preprocessing heuristics. In the top graph of Figure 2, it can be seen that (BnB) is capable of solving the same CIA problem up to several orders of magnitude faster than SCIP. While SCIP takes ≈ 170.6 s to solve the biggest problem instance N 4 , (BnB) is able to solve the problem in ≈ 3.4 s.

Runtime of solution methods and impact of the preprocessing heuristics
The bottom graph of Figure 2 shows that the preprocessing heuristic is able to reduce the number time steps of (CIA) that need to be processed within this benchmark to a high extent, from 43.3 % for N 1 up to 65.8 % for N 4 . By doing this, the solution time for a problem can be reduced by a minimum of approximately one order of magnitude for this example for both solution methods, so that for the biggest problem instance N 4 , the solution time of (BnB) is further reduced to ≈ 0.4 s.

Uniqueness of solutions and optimality after preprocessing
Exemplary, the solution of (CIA) for the data set represented by N 3 = 240 data points obtained by (BnB) is shown in Figure 3. Both maximum-switching-and SOS1-constraints are fulfilled. However, it can be observed that the solutions of the original and the reduced problem with N 3,red = 92 slightly differ.
On the one hand, this situation can occur if the global optimum of a CIA problem is not unique, i. e., there can exist several solutions that are able to yield the same objective value. In that case, the result is still a global optimum, while different binary trajectories are obtained due to different execution orders of the solution methods. On the other hand however, since the preprocessing heuristics remove degrees of freedom from (CIA), situations might also occur where the global optimum of the reduced problem is in fact different from the solution of the original problem, which cannot be obtained anymore due a reduction of degrees of freedom in the modified problem setup.
Within the scenario shown in Figure 3 the latter holds true, and the obtained objective of the original problem is η 3 = 0.4171 and η 3,red = 0.534 for the reduced problem. This reflects a trade-off between runtime reduction and optimality possibly introduced by the preprocessing heuristics and which can be influenced, i. a., by adjustment of the column span parameter δ shown in Algorithm 2.

Design and implementation of the MPC algorithm
Within this section, we describe the methods and software used for implementation of the MPC algorithm and highlight specific aspects. The overall procedure and the main stages of the implemented algorithm are illustrated in Algorithm 3.

Discretization of the MIOCP and utilized solvers
The discretization of the MIOCP is conducted using direct collocation with Lagrange polynomials with Radau collocation points of order 3 (Biegler, 2010) and implemented using the open-source dynamic optimization framework CasADi (Andersson et al., 2018) in Python. For solution of NLPs, we use Ipopt (Wächter and Biegler, 2006) with linear solver MA27 (HSL, 2017). 6 The CIA problem is solved using the pycombina software package presented in the previous section.

Initial guess generation and warm-starting of subsequent MPC steps
To facilitate high performance in terms of solution quality and runtime of NLP solution methods, careful initialization of optimization variables is required.
For initialization of the first NLP within an MPC loop, initial guesses can be obtained, e. g., by simulation of a system f C (·) which, in addition to f (·), contains a conventional control scheme based on PID-and set-point-based control. We use Modelica to implement f C (·) and OpenModelica (The Open Source Modelica Consortium, 2017) for simulation of f C (·) via the OpenModelica Python interface OMPython (Lie et al., 2016).
For initialization of the next MPC solution step, it is particularly useful to re-use the results of a previous optimization cycle. If applied reasonably, warm starting of solvers can save valuable computation time within MPC applications. For warm starting of Interior Point Methods (IPMs), a methodology is presented by Zanelli et al. (2017) which we partially adopt for warm-starting of Ipopt.
However, the solution of the second NLP with binaries fixed is not necessarily a good initialization for the first NLP of the next MPC step where binaries are both free and relaxed again, since the solution of the second NLP heavily relies on the results of the approximation of the binary controls. To overcome this issue, we first divide our MPC algorithm into an online part that needs to be executed as fast as possible right after new data for the system state has been obtained, and a subsequent offline part where more computation time is available. Within the online part of the algorithm, optimized system controls are obtained, while we use the offline part to determine an initialization for the next upcoming MPC step. To achieve that, once the online part has finished and optimized controls for the binary problem have been obtained, we start to solve the same NLP again in its relaxed form, but only until a certain barrier parameter is reached. Once new data for the system state and for the time-varying parameters has been obtained, we update x 0 and c and afterwards continue to solve the relaxed NLP, which is now parametrized for the new MPC step.

Algorithm 3: Main stages of the implemented MPC algorithm
Input : Discretized (MIOCP) instance with initial time grid G N , initial guess generator G, forecast source F , field interface I, user-defined stopping criterion S (optional) Parameters: warm-start IPM barrier parameter µ ws , final IPM barrier parameter µ f 1 Initialize timer τ ← G 0 N and start timer; 2 Obtain c through F and current system state x 0 through I; 3 Obtain initial guess x init , u init , b init given c, x 0 using G; 4 while not S fulfilled do

Simulation study for a solar thermal climate system
Within this section, the framework introduced within the previous sections is applied for the model of a Solar Thermal Climate System (STCS) within a simulation study. First, we give a detailed description of the considered system modeling and its dynamics, components, dependencies and operation conditions. Then we present the MIOCP for the system and the implementation and conduction of the simulation study. Finally, we analyze the obtained results in terms of solution quality, constraint satisfaction and runtime.

System description
The system considered for this study is a modified version of a climate system installed at the Faculty of Management Science and Engineering at Karlsruhe University of Applied Sciences. A schematic depiction of the system is given in Fig. 4.
The core component of the system is a silica-gel/water Adsorption Cooling Machine (ACM) which, when switched on (b acm = 1), utilizes the heat from a stratified High Temperature (HT) water storage of volume V hts = 2.0 m 3 to absorb heat from a stratified Low Temperature (LT) water storage with volume V lts = 1.0 m 3 and emits the combined heat through a recooling tower located on the roof of the building. In times of low ambient temperature and if the ACM is not in operation, the recooler can also be used autonomously in a free cooling or Heat Dissipation (HD) mode (b hd = 1) to cool down the LT storage directly at the ambient without use of the ACM.
The HT storage is supported by an array of horizontally placed solar thermal collectors which are connected to the system through a heat exchanger that separates the outdoor, frost-proof solar circuit from the indoor water circuits. The solar collectors have a total collector surface of A sc = 35.0 m 2 , an optical efficiency of η sc = 0.8 and are exposed to the solar irradiation I g . Regarding the mass flow rateṁ sc of the solar fluid through the collectors, we assume that the pump that produces the mass flow, if turned on, can be set on a continuous scaleṁ sc ∈ [0.1, 0.5] kg/s. 7 The pump on the water side of the heat exchanger behaves identically withṁ hx ∈ [0.1, 0.5] kg/s. The heat exchanger has a volume of V hx = 3.77 · 10 −3 m 3 on each side with heat exchange surface A hx = 4.02 m 2 .
Cooled-down water can be extracted from the bottom of the LT storage to support a set of fan coils operated for heat exchange with the air mass m r,a = 2.16 · 10 3 kg of a room. The water mass flow rateṁ fc,w through the fan coils can be set on a continuous scalė m fc,w ∈ [0.1, 0.3] kg/s if the corresponding pump is turned on. For simplicity, we assume the air mass flowing through the fan coils to be constant atṁ fc,a = 0.43 kg/s. 8 In addition to the fan coils, the room temperature T r,a is influenced by a time-varying heat loadQ load that is usually due to influences from heating through solar irradiation on 7 In the following, we will treat mass flow rates produced by pumps as continuous controls with a minimum flow of 0.0 kg/s and not include the switching behavior into the problem formulation, as explained in the upcoming sections.
8 For real applications, constant operation of fan coils would be inefficient and could instead be, e. g., set depending on pump operation, so that the fans are activated only when the pump is working, which would result in a similar system behavior but consume less electrical energy for fan operation. the building, air exchange, internal heat loads from machines and humans etc. The room air further exchanges heat with the concrete wall mass m r,c = 2.376 · 10 5 kg of temperature T r,c , while the concrete wall mass exchanges heat with the environment depending on the ambient temperature T amb .
The main objective for the system control is to keep the room's air temperature T r,a within a defined comfort range.

System modeling
For modeling of the system components, nonlinear gray-box models based on mass and energy balances are used that fulfill all necessary conditions regarding differentiability for use within derivative-based optimization methods, cf. Biegler (2010). If not stated explicitly, energy losses of the components to their surroundings are neglected. Material and media are assumed incompressible and with constant material properties as given in Table 1.

Solar collector model
To calculate the solar collector temperature T sc , the single-node solar collector model is used, where C sc = 2.6 MJ/K is the lumped heat capacity of the solar collectors including the contained medium and α sc = 1.4 W/(m 2 K) is the heat transfer coefficient for the heat losses of the collector to the environment, cf. Wesselak et al. (2013).

Heat exchanger
Using a lumped model for the heat exchanger and a heat transfer coefficient α hx = 3150.0 W/(m 2 K), the temperatures on the solar side T hx,sc and the water side T hx,hts of the heat exchanger can be calculated aṡ where m hx,sc = ρ sf V hx = 3.77 kg and m hx,hts = ρ w V hx = 3.77 kg.

HT water storage
To facilitate depiction of the temperature distribution in the stratified HT storage tank, we discretize the volume V hts of the storage into L = 4 equally big layers so that the water mass of a layer is m hts = (ρ w V hts )/L = 500 kg and calculate an energy balance for each layer to determine its temperature, cf. Streckiene et al. (2011), Eicker (2003. Magnitude and direction of mass flows between these layers depend on the current mass flow through the heat exchangerṁ hx and ACM operation status b acm ∈ {0, 1}, while we assume thaṫ m hx (t) <ṁ acm,ht at all times. The energy balance that determines the temperature T hts,1 of the top layer can then be calculated according to (9), the temperatures for the the middle layers are determined by (10) and for the bottom layer by (11).

Recooling tower
We assume the temperature of the water 9 returning from the sufficiently sized recooling tower to be of a constant difference ∆T rc = 2 K above the current ambient temperature T amb . Depending on the current operation mode, the medium is either supported for recooling of the ACM or for direct cooling of the LT storage, as shown in the following sections.

Adsorption cooling machine
For operation of an ACM, upper and lower limits on a machine's inlet temperatures must be considered, while low Medium Temperature (MT) input temperatures and high HT and LT input temperatures are favorable conditions for efficient operation (Chang et al., 2007). Within these limits, both the mean cooling powerQ acm,lt and the mean Coefficient Of Performance (COP) of an ACM during one adsorption cycle are determined and can be depicted by, e. g., suitable (polynomial) curve fittings fQ acm,lt (·) and f cop (·). The energy balances for the relevant machine circuits can then be calculated as iṅ Q acm,lt (t) = fQ acm,lt (T hts,1 (t), T lts,1 (t), T amb (t) + ∆T rc ) Q acm,ht (t) =Q acm,lt (t) f cop (T hts,1 (t), T lts,1 (t), T amb (t) + ∆T rc ) (12) and using the mass flow ratesṁ acm,lt = 48 kg/min andṁ acm,ht = 41.8 kg/min, the corresponding output temperatures of the circuits as in The typical cyclic output temperature behavior of an ACM is not depicted by this model. However, the model becomes applicable due to the storages connected to the HT and LT side of the machine and the comparatively high mass flow rate of the MT circuitṁ acm,mt = 85 kg/min which dampens the effect of fluctuating output temperatures, see Sawant and Pfafferott (2017), Sawant and Doan (2017) and Bürger et al. (2017).

LT water storage
The model for the LT storage is formulated analogously to the model of the HT storage. In addition to the current ACM status b acm and water mass flow rate through the can foilṡ m fc,w , the mass flows in between the M = 3 modeled storage layers depend on the current status of the free cooling mode b hd ∈ {0, 1}. The temperature T lts,1 of the top layer of the LT storage is calculated as in (14), the temperature for the the middle layer is determined by (15) and for the bottom layer by (16).

Fan coil and room models
The fan coils are modeled as one single gas-liquid heat exchanger and by a single-node model with heat transfer coefficient (Aα) fc = 475 W/K. With m fc,w = 3.6 kg the water mass and m fc,a = 0.79 kg the air mass inside the fan coil, the temperatures of the water side T fc,w and the air side T fc,a of the fan coils are given bẏ The temperature of the room's air mass T r,a and the temperature of the concrete mass of the wall T r,c are determined by single-node models given in (18). The air and concrete masses exchange heat depending on the wall surface A c = 540 m 2 and a heat transfer coefficient α a,c = 10.0 W/(m 2 K). The concrete mass additionally exchanges heat with the environment depending on the same factors and the ambient temperature T amb .

Model summary
To summarize, the system is described by n x = 14 differential states which are all temperatures of system components, and influenced by n u = 3 continuous controls, n b = 2 binary controls and n c = 3 time-varying parameters

Operation conditions and constraints
System temperatures should not exceed T ub = 105 • C or go below T lb = 5 • C. The limits on the inlet temperatures for operation of the ACM are given in Table 2. Since there is only one recooler present in the system, b acm and b hd cannot be active at the same time, i. e., b acm (t) + b hd (t) ≤ 1.
Both ACM and HD are considered subject to dwell-time-and maximum-switchingconstraints. The minimum operation time for the ACM should be 120 min and for the HD mode 60 min. The operation status of ACM and HD should not be changed more than σ acm,max = σ hd,max = 4 times within 24 h.

Temperature Lower bound
Upper bound The Lagrange term of the objective (22a) contains the sum of squares of the continuous controls u and the slack variables s, weighted by appropriate weighting matrices W u ∈ R nu×nu and W s ∈ R ns×ns . While minimization of u favors energy efficient pump use, minimization of s reduces the deviation ∆T r,a of the room temperature T r,a from an assumed comfort range [T r,a,lb , T r,a,ub ] = [20 • C, 22 • C] in (22c) as well as the use of s lb and s ub for relaxation of the temperature boundaries for ACM operation (22d) and (22e) and the use of s T for relaxation of the general system temperature limits (22f). Constraint (22i) ensures that at most one binary control is active at a time.
The Mayer term contains the storage layer temperatures T hts,1 (t f ) and T lts,3 (t f ) at the final time point of the control horizon t f , weighted by appropriate scalar weightings w T hts,1 and w T lts,3 . Since high temperatures for T hts,1 and low temperatures for T lts,3 are generally advantageous, it is favorable to drive the system towards such a configuration.
The dynamics of the system that must be fulfilled are given in (22b). Path constraints that enforce compliance of the temperature bounds for ACM operation are introduced in (22d) and (22e) as smoothened vanishing constraints (cf. Section 2.4). While the boundaries on the continuous controls are given in (22g) as hard constraints, the boundaries on the system states in (22f) are formulated as soft constraints to preserve feasibility of the problem in case of state constraint violations, cf. . In addition to that, the binary constraint (22h) must be fulfilled as well as the initial state constraint in (22j). Additionally, the system is subject to dwell time constraints and maximum switching constraints, but which are not considered within solution of the NLPs but only within (CIA), since such constraints increase the size of the NLP very much, while they in fact shorten the solution time for (BnB) in pycombina.

Setup and implementation
The simulation study is implemented according to Section 4.1. Since we consider a solar driven system, we can assume more need for control interaction during day than during night, so for a control horizon of t f = 24 h we set up the time grid using time steps of ∆t d = 5 min for day times in between 4:40 and 21:30 10 , and longer time steps of ∆t n = 20 min for night times in between 21:30 and 4:40 to reduce the number of optimization variables in the MINLP. The resulting time grid G N consists of N = 226 time points. To avoid regeneration of the discretization in between MPC steps, we apply the time transformation as described in Section 2.3 so that the current time steps for each grid point enter the problem as time-varying parameters. In total, the MINLP contains 18000 continuous and 450 binary optimization variables as well as 12600 equality and 4950 inequality constraints.
The simulation model that we implement in Modelica can either be run using a PID-and set-point-based, conventional control scheme that is contained directly within the Modelica model, or it can be run using the control inputs from the solution of the MINLP. Further, the Modelica model includes a safety monitoring of the solar collectors temperature T sc that assures that once T sc exceeds a critical temperature, the mass flow rate through the collectorsṁ sc is always set to its maximum value, e. g., when due to imprecise forecasts the temperature development inside the collectors is underestimated.
Mass flow rates requested by the controller that are below a pump's minimum flow rate introduced in Section 5.1 are treated by pulsing the corresponding mass transported over an interval for the requested flow rate within two shorter pulses at the corresponding minimum flow rate. The necessary pulsing time ∆t k,p within an interval k is then obtained from the simple relation while the pulsing of mass is initiated at time points t k and t k + ∆t k 2 . This approach becomes applicable only due to the short start-up-time of the pumps and is not suitable for, e. g., operation of the ACM.

Conduction of the simulation study
Using the implementations and methods described in the previous sections, the simulation study for the system is conducted as follows: 1. Starting from an initial state of the system at 2:00 am at night of an initial guess for the MINLP solution algorithm is obtained from simulation in OpenModelica using the included conventional controller and forecast values for solar irradiation I g,fc , ambient temperature T amb,fc and heat loadQ load,fc . 2. Then, the MINLP is solved using the initial guess obtained in the previous step and the forecasts I g,fc , T amb,fc andQ load,fc . 3. After that, a simulation is conducted in OpenModelica applying the control inputs obtained from the solution of the MINLP, but using measured values for solar irradiation I g,m , ambient temperature T amb,m and heat loadQ load,m that deviate from the forecast values at a varying extend. 4. The result of the simulation for the upcoming time point is then used as initial state vector x 1 for the next MPC steps and corresponding MINLP. 11 We repeat this process to simulate an MPC loop of 1000 steps, which corresponds to more than 100 h of system operation time. The profile of forecasted and corresponding measured solar irradiation I g , ambient temperature T amb and heat loadQ load used within this study are depicted in Figure 5. 12 The heat load is calculated for both the measured loadQ load,m and forecasted loadQ load,fc by a simple relation to a reference temperature of T (25)

Simulation results
The result of the simulation is shown in Figure 6. The upper two plots show the development of the HT and LT side of the system on selected system states. The second lowest plot shows the continuous controls, while the bottom plot depicts binary controls.
The graph shows that the MPC is able to drive the system in a desirable way. The controller heats up the solar collector temperature T sc at low pump operationṁ sc anḋ m hx and increases pump operation with rising temperature T sc to charge the HT storage, observable from increasing storage temperatures T hts,1 and T hts,4 . The ACM is operated to decrease the temperature of the LT storage depicted by the storage layer temperatures T lts,1 and T lts,3 . For doing this, the controller operates the ACM (b acm = 1) not only during occurrence of thermal loads, but also at different points in the early day or during night to benefit from the low ambient temperature T amb and resulting in a high COP for the ACM to predictively regenerate the LT storage in preparation for high loads of an upcoming day, cf. Bürger et al. (2017). Free cooling (b hd = 1) is thoroughly used during times of low ambient temperature to cool down room air and concrete mass especially during nights. Nevertheless, the system is not completely able keep the room temperature T r,a within its desired comfort region at all times, which is due to the (intentionally chosen) comparatively high heat load affecting the system compared to the system's cooling capacities.
Despite mismatches between forecasted and occurring solar irradiation, ambient temperature and heat load profiles, the results in Figure 6 show that the MPC is capable of controlling the system sufficiently. However, imprecise forecasts can cause suboptimal system operation, as shown, e. g., by the fast switching of the ACM around hour 80. Also, underestimation of the solar irradiation might lead to an unexpected increase in the solar collector temperature, and with this, to an intervention of the solar temperature safety monitoring or violation of temperature boundaries. Figure 6 shows that the dwell time constraints are always fulfilled. Regarding constraint satisfaction, violations of the operation temperature constraints T ht,lb and T lt,lb as in

Constraint satisfaction
as well as a violation of the maximum system temperature T ub for the solar collector temperature T sc result in utilization of the corresponding slack variables for relaxation of the constraints to preserve feasibility of the NLPs. Apart from these, no operations constraints are violated. The actual values of violations are shown in Figure 7. One reason for constraint violation can be found within the approximation of the binary controls within CIA, which might result in scheduling of machine operation in times where corresponding operation conditions are not (yet) completely met, cf. Bürger et al. (2018). Also, mismatches between forecasted and measured time-varying parameter values result in under-and/or overestimation of system temperatures. In order to avoid such violations, larger safety margins for the constraints within the OCP formulation could be introduced, or the quality of the utilized forecasts could be improved.
In case of larger constraint violations, internal safety measures similar to the solar temperature monitoring included in this study can be utilized to perform, e. g., a safety shutdown of a machine. For the components used within this study, however, we can presume constraint violations of the presented dimension, which are always smaller than 1 K, acceptable.  Figure 8 shows a box-plot for the individual runtimes for the several solution stages of the algorithm as well as the total duration of the online and offline parts and the overall solution time of one MINLP solution process. As customary, the box ranges from the 25 % to the 75 % quartile. The whiskers of the boxes extend these by 1.5 times the distance between the 75 % and 25 % quartile, which is the so-called Inter-Quartile Range (IQR). The remaining values outside of this range are regarded as outliers, denoted as dots within the plot. The median is denoted as a vertical line, however for NLP init which is evaluated only once at the start of the MPC loop the vertical line indicates the duration of this single run. Figure 9 shows an additional plot for the relation between online and total solution times. The corresponding minimum and maximum values shown in the plots as well as the values of the median and quartiles are depicted in Table 3.

Runtime of solution stages
From the plots and  Table 3: Numerical values of the properties shown in Figure 8.
MINLP takes less than 26 s, and less than 33 s for more than 75 %. In 50 % of the cases, the online part of the algorithm takes up less than 63 % of the total runtime, and less than 69 % in 75 % of the cases. Both the maxima of the online and total runtimes are far below the smallest element of the discretization time grid G N . All CIA solution times obtained using the branch-and-bound method of pycombina after preprocessing are below 0.5 s. The majority of the algorithm runtime is spent within solution of the NLP problems, while similar shares are spent within preparation of the next NLP solution step NLP prep and the solution of NLP rel . All NLPs solved fully within this study converged to a feasible optimal solution.
Since it has been shown in Bürger et al. (2018) that the solution time for a similar but less complex problem can already increase up to several hours when a general MINLP solver is used, we apply Bonmin with default settings for solution of the first MINLP of the simulation study for comparison, however neglecting the complex dwell-time-and maximum-switchingconstraints for simplicity. Still, it takes approx. 10 hrs for Bonmin to find an optimal integer solution for the problem, while our decomposition approach is able to find a solution of comparable quality in approx. 50 s.
Amongst others, such solution times for a problem of the given complexity renders the algorithm suitable for applications within process control and control of energy and HVAC systems. This field of applications could possibly be enlarged if the runtime of the algorithm was further reduced, e. g., by application of tailored NLP solution methods. Though Ipopt is tailored for solution of large-scale, sparse systems (Vigerske et al., 2016), further possible reductions in runtime could be achieved by application of suitable Sequential Quadratic Programming (SQP) solvers which are intrinsically more suitable for warm-starting, cf. . In addition to that, the C-code generation and parallel computation facilities within CasADi could be utilized to further speed up computations. 5.9. Comparison to a conventional control scheme Finally, a system simulation using the Conventional Controller (CC) as applied for initial guess generation of the first MPC step is conducted in OpenModelica and its performance and behavior is compared to the simulated MPC loop presented in the previous sections.
Since the CC is tuned for scenarios of high heat load occurrence as used within this study, behavior and performance of the CC shown in Figure 10 at first glance appear quite similar to the results of the MPC. However, a closer investigation reveals several aspects that illustrate the benefits of the MPC application, as described in the following.

Solar collector temperature
Since the CC has no information on the future development of the temperature T sc , the controller needs to set a high mass flow through the solar collectors in order to prevent possible overheating caused by fast temperature increases, as shown on the high values of the controlṁ sc .
The MPC on the other hand can directly take into account the estimated development of T sc , and therefore typically increases T sc earlier during the day and to a comparatively higher temperature. This can be observed, e. g., during the second and third day of the simulation study. Due to the quadratic contributions ofṁ sc to the MPC objective which approximates the superlinear increase of electricity consumption of pumps at increasing speed, collector operation is realized at comparatively lower mass flows and therefore in a more energy-efficient way.

Rules for ACM and HD activation
Within the CC, specific and hierarchical rules for activation of ACM and HD cooling based on current system temperatures need to be defined. As observable in Figure 10, this results in utilization of HD cooling during night and ACM cooling during day.
In contrast to that, the MPC can make flexible decisions based on the current situation, which can be observed, e. g., during the night between the second and third day of the study. There, since the comfort region of T r,a is already reached, the MPC favors ACM cooling for preparation of the LT storage for upcoming loads during efficient ACM operation conditions, contrary to the CC which utilizes HD cooling for further reducing the room temperature.

Reference temperature and utilization of the comfort range
For the CC, the comfort range for the room temperature allows for setting of a reference temperature and a corresponding hystereses at which cooling needs to be activated or deactivated, respectively. Once the cooling is active, the mass flowṁ fc,w is determined based on the current deviation of T r,a from its reference value. However, this behavior is fixed and independent of the expected heat loads of the upcoming hours or days.
Within the MPC, the comfort range introduces comparatively more flexibility since the controller can decide on its utilization in a situation-dependent way. This allows to decrease room temperature, e. g., in time of low heat loads and in this way prepare for expected higher loads during the upcoming day. Here, setting a lower reference temperature for the CC however could result in a room temperature that is generally lower than necessary, resulting in a unnecessary high energy usage.

Conclusions and outlook
In this work, we presented an algorithm for MPC of switched nonlinear systems under combinatorial constraints and a corresponding new toolbox for solving CIA problems with special focus on application within MPC including preprocessing heuristics for complexity reduction of CIA problems. Within a simulated MPC loop, the implemented algorithm was applied for the case study of an STCS with nonlinear system behavior and combinatorial constraints whose operation conditions are subject to uncertainty. It has been shown that the algorithm is capable to efficiently control the system while preserving sufficient constraint satisfaction and algorithm runtimes.
Future work should focus on further runtime reduction by applying more tailored NLP solution methods. Also, the application of the algorithm for other systems or MPC ap-plications with possibly increasing complexity may be investigated and finally the physical installation of the STCS presented within this work is desirable.