Experimental operation of a solar-driven climate system with thermal energy storages using mixed-integer nonlinear model predictive control

This work presents the results of experimental operation of a solar-driven climate system using mixed-integer nonlinear model predictive control (MPC). The system is installed in a university building and consists of two solar thermal collector fields, an adsorption cooling machine with different operation modes, a stratified hot water storage with multiple inlets and outlets as well as a cold water storage. The system and the applied modeling approach is described and a parallelized algorithm for mixed-integer nonlinear MPC and a corresponding implementation for the system are presented. Finally, we show and discuss the results of experimental operation of the system and highlight the advantages of the mixed-integer nonlinear MPC application.

In this work, we present the real-life experimental operation of such a system installed at Karlsruhe University of Applied Sciences using mixed-integer nonlinear MPC. The system consists of solar thermal collector fields, an adsorption cooling machine (ACM) with different operation modes, a stratified hot water storage with multiple inlets and outlets as well as a cold water storage. Utilized for air conditioning of a part of the university building, the system operates under real conditions. As the system shows highly nonlinear behavior, complex component interaction, delays, and strong dependencies on external factors such as solar irradiation and ambient temperature, the plant is suitable for experimental studies on mixed-integer nonlinear MPC of renewable energy systems.
The contributions of this work are as follows. First, we extend the model presented in our previous work 30 and provide insight on the control-oriented modeling of the solar thermal collector fields, which are connected in parallel and show a nonlinear mass flow distribution as a function of pump and mixing valve operation, as well as on the modeling of the stratified hot water storage with multiple inlets and outlets. Then, we introduce a parallelized and modified version of the algorithm for mixed-integer nonlinear MPC used in our previous simulation study 30 that enables faster updates of control decisions and apply strategies for further decrease of algorithm runtime. Finally, we present and discuss experimental results of real-time MPC operation for the physical system on a day with challenging operation conditions and highlight the advantages of mixed-integer nonlinear MPC application for the system.
The remainder of this work is organized as follows. In Section 2, the solar thermal climate system subject to this study is introduced, followed by a description of a control-oriented nonlinear model for the system in Section 3. The mixed-integer nonlinear optimal control problem (OCP) for the system and a corresponding moving horizon estimation (MHE) problem for state estimation are presented in Section 4, including a description of the applied solution methods, MPC algorithm, and a selection of important implementation details. Finally, in Section 5, the results of experimental operation of the system are shown and discussed. Section 6 concludes this article.

DESCRIPTION OF THE SOLAR THERMAL CLIMATE SYSTEM
The system considered for this study is a solar-thermally-driven climate system installed at the Faculty of Management Science and Engineering at Karlsruhe University of Applied Sciences. The system is operated for covering cooling loads of the building's atrium during summer time and for heating support during winter time. Within this study, we will focus on summer operation, which is the more comprehensive case as it involves operation of more system components compared to winter operation and in particular switching of the ACM. A more detailed description of the system and its components is given in the upcoming sections.

Adsorption cooling machine
The central component of the system is a silica-gel/water ACM which is located in the cellar of the building, see Figure 1. The ACM is connected to a recooling tower (RT) on the roof of the building which is operated in a frost-proof circuit, see Figure 2, via a heat exchanger in the cellar further entitled as recooling tower heat exchanger (RTHX). The ACM can operate in three different modes, which can only be switched on or off completely and whereof at most one mode can be active at a time. Operation of the machine internal pumps and valves as well as the RT pump and fan speeds for each mode are subject to a dedicated internal controller of the ACM. * Running in adsorption cooling (AC) mode (b ac = 1), the ACM utilizes the heat from a stratified high-temperature storage (HTS) with volume V hts = 2 m 3 to cool down a stratified low-temperature storage (LTS) with volume V lts = 1 m 3 , heat is emitted through the RT. In this mode, cooling power and COP of an ACM depend on the inlet temperatures of the three circuits of the machine, while higher high temperature (HT)-and low temperature (LT)-and lower recooling-temperatures are favorable. 33 In times of low ambient temperature, the machine can run in a free cooling (FC) mode(b fc = 1), which allows to cool down the LTS via the RT. Using the third mode, which is entitled emergency cooling (EC) mode (b ec = 1), the ACM can be utilized to cool down the HTS at the ambient via the RT to prevent overheating of the HTS, and with this, overheating of the solar collectors. * While we focus on the optimal switching of an ACM within this work, it should be noted that MPC can also yield large potential for optimization of the internal processes of an ACM. 32 F I G U R E 1 System components in the cellar of the building: (1) control cabinet, (2) low temperature storage, (3) adsorption cooling machine, (4) high temperature storage, (5) pump P ssc , (6) pump P psc , (7) pump P lc , (8) solar heat exchanger, (9) recooling tower heat exchanger, (10) inlet I hts,t , (11) outlet O hts,m , (12) outlet O hts,b , (13) mixing valve M ssc , (14) outlet O hts,t , (15) inlet I hts,m , (16)  The temperature at the inlets and outlets of the ACM and the current flows in each circuit are measured, as well as the temperatures at the inlets and outlets of the RTHX and the outlet temperature of the RT.

High-temperature circuits
The heat stored in the HTS is collected by two roof-top solar thermal collector arrays as shown in Figure 2. One array consists of horizontally placed vacuum tube solar collectors (VTSC) with total area A vtsc = 31.2 m 2 and the other consists of flat plate solar collectors (FPSC) with total area A fpsc = 53.3 m 2 mounted at a 45 • angle facing south. Both arrays are connected to a circuit entitled primary solar circuit (PSC) with approximate volume V psc = 0.4 m 3 . The PSC is filled with an anti-freeze brine and is separated from the indoor water circuits by a solar heat exchanger (SHX) located in the cellar, see Figure 1. The amount of liquid flowing through the PSC and the individual collector arrays can be controlled by modulating the speed v ppsc of the pump P psc and the position p mpsc of the mixing valve M psc which is located in the return line of the PSC. The water-filled circuit on the other side of the SHX is entitled the secondary solar circuit (SSC) and is directly connected to the HTS. The amount of liquid flowing through the SSC can be controlled by modulating the speed v pssc of the pump P ssc .
Water returning from the SHX to the HTS in the SSC is always supplied into the top inlet I hts,t of the HTS as shown in Figure 1. However, the water flowing toward the SHX can be extracted from the middle outlet O hts,m , the bottom outlet O hts,b , or a combination of both, depending on the position p mssc of the mixing valve M ssc . For hot water supply of the ACM the top outlet O hts,t of the HTS is used. The water returning from the ACM to the HTS can be inserted into the middle inlet I hts,m , the bottom inlet I hts,b , or a combination of both, depending on the position p macm of the mixing valve M acm .
The temperature sensors T vtsc and T fpsc for the collector arrays are installed directly at the array outlets in the return line of the PSC. Therefore, the array temperatures can not be measured directly, but only via the temperature of the medium returning from the arrays, that is, only when P psc is operating. Additionally, the ambient temperature and the current solar irradiation on the roof are measured using a temperature sensor T amb and a pyranometer I g , respectively. Further, the temperatures at each inlet and outlet of the SHX are measured using the sensors T shx,psc,i , T shx,psc,o , T shx,ssc,i F I G U R E 3 System components in the atrium of the building: (1) location of installation North-West, (2) location of installation North-East, (3) fan coil units, (4) sensor for measurement of concrete temperature, and (5) detail view of installation South-West [Colour figure can be viewed at wileyonlinelibrary.com] and T shx,ssc,o . Also, the inlet and outlet temperature of the HTS can be monitored, as well as the flows in each circuit. Each mixing valve is equipped with a potentiometer which allows to determine the valve's current position. Also, n hts = 4 temperature sensors T hts,k , k = 1, … , n hts are inserted into thermometer pockets installed at different heights of the HTS.

Low-temperature circuits and room installations
The atrium of the building, which is depicted in Figure 3, has a length of l atrium = 17.2 m and width of w atrium = 11 m. With a height of h atrium = 12.7 m, the room extends over three floors. The glass roof of the atrium can be covered using outdoor blinds, which are typically closed during summer and cover most of the glass surface.
Behind the plumbing wall covered by wooden panels, four fan coil units (FCUs) are installed as shown in Figure 3. These are supplied by chilled water from the bottom of the LTS, which is then returned from the FCUs into the top of the LTS. This circuit is entitled the load circuit (LC) and the medium flowing through this circuit can be controlled by modulating the speed v plc of the pump P lc . † On four of the depicted black steel beams, installations for temperature measurement are attached, see Figure 3. According to their locations, these are labeled North-West (NW), North-East (NE), South-West (SW), and South-East (SE). Each installation consists of 10 temperature sensors, of which seven are used to measure air temperature at different heights and the remaining three are attached to the several floors using thermal grease for measurement of concrete temperature. Further, the water temperature at inlet and outlet of the FCUs is measured via the temperature sensors T fcu,i and T fcu,o as well as the total flow in the LC, and n lts = 2 temperature sensors T lts,k , k = 1, … , n lts are inserted into thermometer pockets at the top and the bottom of the LTS.

Operation conditions and system boundaries
Water circuit temperatures should not exceed T w,max = 95 • C or go below T w,min = 5 • C. Temperatures in the PSC should not exceed 110 • C to prevent the solar liquid from boiling and changing to a gaseous state, since this would result in a shutdown of the system until all medium has returned to a liquid state, which can typically not be ensured before the upcoming day. To include a certain tolerance region, we therefore define that the temperature of the solar liquid should not exceed T sl,max = 98 • C. As a further safety measure, P psc must be in full operation as soon as one solar collector array temperature exceeds T sc,so = 65 • C. Operation of the ACM in FC mode is allowed only at ambient temperatures T amb ≥ 4 • C. For operation in AC mode, the temperature of the water entering the HT side of the machine must not exceed T ac,ht,max = 95 • C and should not go below T ac,ht,min = 55 • C, the temperature of the water entering the LT side of the ACM must not go below T ac,lt,min = 10 • C, and the ACM should only be active at ambient temperatures above T ac,amb,min = 15 • C and below T ac,amb,max = 38 • C. Further, minimum operation times of 60 minutes for the AC mode and 15 min for the FC mode are defined.

Programmable logic controller
The control signals are sent to the several components using a programmable logic controller (PLC). This controller can operate the system either using a conventional, proportional-integral-derivative (PID) and set-point-based control scheme that is implemented directly on the PLC itself, or it can operate the system using external control signals that the controller receives via a network interface. Independently from the currently used control strategy, the PLC monitors the system temperatures and detects violations of previously specified safety limits. In such a case, the PLC automatically takes over control of the system, and aims to rapidly restore valid operation conditions. For example, in case the temperatures in the PSC exceed 108 • C, the EC mode of the ACM is automatically activated.

CONTROL-ORIENTED MODELING OF THE SOLAR THERMAL CLIMATE SYSTEM
The system is modeled using nonlinear grey-box models based on mass and energy balances, resulting in a system of explicit ordinary differential equations (ODEs). The model fulfills all necessary conditions regarding differentiability for later application within derivative-based optimization methods. 34 Material and media are assumed to be incompressible and with constant material properties. A schematic depiction of the system model is given in Figure 4. Model parameters are either taken from manufacturer data sheets, if possible, or are the results of system identification and parameter estimation procedures based on measurement data. In the following, component-and circuit-models are described in more detail and information on the corresponding system identification processes is given.

Mass flows, pumps, and mixing valves
For simplicity, we do not model pumps and mixing valves explicitly, but assume that the mass flows in the circuits of the system can be controlled directly and on a continuous scale starting from zero. For this to be applicable, a priori measurements from the physical system that determine the volume flows resulting from different pump speeds and mixing valve positions have been obtained. For pump modeling, we also neglect their minimal part load behavior. Therefore, it must be ensured that if the controller requests a mass flow rate below the minimal rate that a pump can produce, the amount of mass transport resulting F I G U R E 5 Exemplary depiction of mass flow pulsing: in the depiction,ṁ request is the mass flow requested by the controller for the time interval [t 0 , t 1 ], which is smaller than the minimum mass flowṁ min that the concerned pump can produce. Therefore, the pump is operated several times at its minimum flow rate during the interval, such that ∫ from the request over a certain period of time is realized for the physical system by pulsing at minimal flow rate, cf.
For modeling the effect of the mixing valve M ssc , we assume that the flow separation between outlets O hts,m and O hts,b depends linearly on the position p mssc ∈ [0, 1], so that the mass flows through O hts,m and O hts,b are given bẏ However, to avoid the increase of model nonlinearity introduced by the multiplication of two control variables, we introduce the mass flowṁ o,hts,b through O hts,b as an individual optimization variable instead of p mssc , and formulate accordinglyṁ Variableṁ o,hts,b enters the model as a state of the system as and the change of the mass flow rate Δṁ o,hts,b as a control variable. Sinceṁ o,hts,b is introduced as a substitute for the mixing valve position p mssc , this formulation facilitates to penalize position changes for M ssc via inclusion of Δṁ o,hts,b in the mixed-integer optimal control problem (MIOCP) objective, and with this, to avoid frequent opening and closing of the valve. Analogously, for inlets I hts,m and I hts,b and their mass flowsṁ i,hts,m andṁ i,hts,b with ‡ This approach becomes applicable due to the assumption that the pumps are typically operated at higher speeds and their comparatively short startup time of approximately 1-2 s, and is therefore not suitable for modeling of the switching behavior of the ACM. m i,hts,m (t) =p macm (t)b ac (t)ṁ ac,ht , and p macm ∈ [0, 1], we introduceṁ i,hts,b as an individual optimization variable, and formulate accordinglẏ while the HT side mass flow induced by the ACM is constant atṁ ac,ht = 0.69 kg s when the machine is operating in AC mode. Variableṁ i,hts,b also enters the model as a state of the system as and the change of the mass flow rate Δṁ i,hts,b as a control variable. The corresponding mixing valve positions for the PLC to set are determined from the ratios of the requested mass flows. The flow sensor installed in the PSC solely measures the total flow in the circuit but not in the individual collector arrays. Therefore, the flow in each array depending on P psc and M psc was determined for different speeds v ppsc and valve positions p mpsc from the individual array outlet temperatures, the total flow, and the temperature resulting from mixing of both flows at the outlet of M psc . Since the individual flows can only be identified from their mixing temperature when the temperature difference between flows is sufficiently high, one collector array was covered using tarpaulins during the measurement process. Afterward, the process was repeated with the other array covered. Since density and viscosity of the solar liquid are influenced by its current temperature but not explicitly considered, the identification was carried out with liquid temperatures above 40 • C so that the identified flows match the real occurring flows better at higher temperature, which is the desired operation state and where inaccuracies are assumed to be more critical.
The resulting flows are depicted in Figure 6 and show a highly nonlinear correlation. To depict these correlations in a suitable way and respecting the necessity for differentiability of the model, suitable B-spline interpolation functions 35 vtsc (⋅) and fpsc (⋅) are used to depict the individual mass flowsṁ vtsc andṁ fpsc aṡ The valve position p mpsc is considered as a state of the system with dp mpsc (t) dt = Δp mpsc (t) (13) and the change of the position Δp mpsc as a control variable.

Influence of external factors
The system is influenced by external factors, which we distinguish between ambient conditions and indoor conditions. Since no facilities for measurement or forecasting of indoor conditions such as building usage and occupancy are available, indoor conditions are in the following considered as noise acting on the system which is estimated during operation.
As ambient conditions, we regard the ambient temperature T amb and several quantities of solar irradiation. More precisely, we consider the total radiation on the VTSC I vtsc and on the FPSC I fpsc as well as the direct irradiation I r,dir and diffuse irradiation I r,diff acting on the building, which are adapted according to the geometrical properties of their corresponding model components. For these quantities, forecast data can be obtained.
Since the quality of the utilized forecasts highly influences the overall MPC performance, forecast data should be adapted based on current measurements of the corresponding quantities. For this, additional states which in the following are entitled as auxiliary states are introduced to the model, whose values are estimated during system operation. Exemplified on the ambient temperature, the term T amb used within modeling of the system is where T amb,fc is the forecast of the ambient temperature and Δ T amb is the auxiliary state to compensate the difference between the forecast and its currently measured value. The development of Δ T amb is described by with Δ T amb an adequate time constant, which results in a decaying influence of the current measurement on future time points. The term w Δ T amb is an additive white process noise term which is used within MHE and considered zero within MPC.

High-temperature storage
The HTS is discretized into n hts = 4 volumes with a water mass m hts = w V hts ∕n hts each, where w is the density of water, to depict the temperature distribution in the stratified HTS, cf. Reference 36. The number of layers has been determined by a trial-and-error approach, where n hts = 4 showed sufficient quality while still incorporating only a reasonable amount of states. The direction and magnitude of mass flows in between those layers depend on their position along the storage height and the current machinery operation. While in simulation models and software alternating mass flows are typically depicted using min/max-or if/else-clauses, these introduce discontinuities to a model and must be avoided to preserve differentiability. Therefore, we formulate the energy balances for the HTS in a differentiable form by following an approach presented by Reference 37. Under the assumption that inlet I hts,t and outlet O hts,t only interact with the top layer of the HTS, we formulate the energy balance that determines the temperature T hts,1 of the top layer of the HTS as with c w the specific heat capacity of water, hts the thermal conductance for the heat loss of the HTS to its ambient § identified from storage temperature measurements, andṁ hts,t,s the sum of the mass flows entering and leaving the top layer. Variableṁ hts,t,s is a smooth approximation for the absolute value ofṁ hts,t,s , which is better the smaller hts is chosen. However, smaller values for hts result in increased nonlinearity of the model, which renders hts a tuning parameter to be chosen based on the actual quantities of the involved mass flows and the desired accuracy of the approximation. In the following, we use hts = 10 −2 . The state T shx,ssc,n shx is the temperature state at the outlet of the SSC-side of the SHX model, which is described in more detail later in this section. The bottom layer of the HTS contains the inlet I hts,b and outlet O hts,b and the energy balance that determines its temperature T hts,n hts is given by with T ac,ht the outlet temperatures on the HT side of the ACM in AC mode,ṁ hts,b,s the sum of the mass flows entering and leaving the bottom layer, andṁ hts,b,s a smooth approximation for the absolute value ofṁ hts,b,s . The HTS layer 3 is assumed to interact only with its neighboring layers 2 and n hts , so that its energy balance is given by Finally, layer 2 of the HTS is assumed to interact with inlet I hts,m and outlet O hts,m in addition to its neighboring layers, so that the energy balance that determines the temperature T hts,2 of this layer reads as

Solar circuits
For modeling of the solar collector arrays, system identification showed that the use of one state per array already provides sufficient quality in depiction of the collector temperatures. Therefore, the temperatures T vtsc and T fpsc of the VTSC and FPSC, respectively, are determined by § The HTS is located in a cellar room of the building whose temperature is influenced by the ambient. Since forecasts for the exact room temperature are not available, current heat losses of an HTS layer to its surrounding are approximated via a relation to the current ambient temperature.
with C vtsc , C fpsc the heat capacities of each collector array including the contained medium, vtsc , fpsc their optical efficiencies, and vtsc , fpsc the heat transfer coefficients for the heat losses to the ambient of each array. 38 For consideration of unmodeled influences on collector heat losses induced by, for example, wind, the auxiliary states Δ vtsc and Δ fpsc are introduced, whose values are estimated during system operation and whose development is described by with w Δ vtsc and w Δ fpsc additive white process noise terms used within MHE and considered zero within MPC. As described in Section 2, the temperature of the collector arrays cannot be measured directly, but only via the sensors T vtsc and T fpsc . Their temperatures T vtsc,s and T fpsc,s are determined by energy balances over a short piece of the PSC pipe as with V vtsc,s , V fpsc,s the pipe volumes for the sensors, vtsc,s , fpsc,s their thermal conductances for the heat losses to the ambient, and c sl and sl the specific heat capacity and density of the solar liquid, respectively. Since the amount of liquid in the PSC is rather large, the mass transport from the roof to the cellar introduces a delayed availability of the heat collected by the solar arrays. This behavior is modeled using two energy balances that depict the temperature of the medium in the feed line pipe T pscf and the return line pipe T pscr , respectively, as follows with C psc the heat capacity of one pipe including the contained medium, psc the thermal conductance for the heat losses of one pipe to the ambient, and T shx,psc,n shx the temperature state at the outlet of the PSC-side of the SHX model. The model for the SHX consists of n shx = 4 energy balances per side of the SHX where opposing volumes exchange heat depending on a heat transfer coefficient shx and heat exchange area A shx . The temperature states T shx,psc,1 and T shx,ssc,1 at the inlets of each side of the SHX are determined by and the remaining temperature states per side of the SHX are determined by

Cooling machine
The ACM is a static model describing the relations between inlet and outlet temperatures of the machine's circuits in its different operation modes, which in turn influence the temperatures of the HTS and LTS once a mode is active. While only the AC mode utilizes the ACM itself, the FC and EC modes only utilize the RTHX. Measurements at the physical system have shown that the LT outlet temperatures of the ACM in FC mode can be sufficiently depicted by a simple temperature difference to the current ambient temperature as The cooling powerQ ac,lt and COP of the machine COP ac in AC mode depend on the inlet temperatures of all three circuits. For the machine used in this study,Q ac,lt and COP ac data can be obtained from the data sheet provided by the manufacturer as averaged values over one complete cycle of the machine and for different combinations of inlet temperatures.
For the ACM model, the temperatures T hts,1 and T lts,1 are regarded as inlet temperatures of the HT and LT side, respectively. For the temperature at the medium temperature (MT) inlet, we assume that the RT returns medium at a constant difference to the ambient temperature analogously to the FC mode. Therefore,Q ac,lt and COP ac can be determined by 30,37Q ac,lt (t) =̇Q ac,lt (T hts,1 (t), T lts,1 (t), T amb (t)), COP ac (t) = COP ac (T hts,1 (t), T lts,1 (t), T amb (t)), witḣQ ac,lt (⋅) and COP ac (⋅) suitable polynomial fittings to machine data. Using these, the temperatures T ac,ht and T ac,lt at the HT and LT outlet of the ACM can be calculated as Both the AC mode and FC mode of the ACM are subject to startup and shutdown phases which can take several minutes to complete, while the phases of the FC mode typically take less time to finish than the AC mode phases. Since these phases are not considered in the developed model, frequent switching of modes, especially from and to the AC mode, must be avoided.
As the EC mode of the machine is considered more as a safety measure than a desirable operation mode, active consideration of this mode within MPC is neglected. However, the mode can still be activated by the PLC in case of violations of temperature safety limits.

Low-temperature storage and fan coil units
System identification showed that the LTS can be depicted by n lts = 2 volumes with water mass m lts = w V lts ∕n lts to depict the temperature distribution in the stratified LTS. Due to the circumstance thatṁ lc,max <ṁ fc,lt andṁ lc,max <ṁ ac,lt witḣ m fc,lt =ṁ ac,lt = 0.8 kg s , the direction of mass flows inside the LTS is directly determined by whether the ACM is currently operating or not. Therefore, the temperatures of the upper layer T lts,1 and the bottom layer T lts,n lts are determined by the energy balances dT lts,n lts (t) with T ac,lt and T fc,lt the LT outlet temperatures of the ACM in AC and FC mode, respectively. The temperatures of the water side T fcu,w and of the air side T fcu,a of the FCUs are determined by Q fcu,r (t) = fcu (T fcu,w (t) − T fcu,a (t)), with C fcu,w and C fcu,a the thermal capacities of the water and air side of the FCUs, respectively, fcu the thermal conductance between those sides, T r,a,1 the temperature of the bottom air layer of the room model, andQ fcu,r the total heat exchanged between FCUs and room air. The air mass flowṁ fcu,a is assumed constant within the model, however, for the physical system the FCUs are only operated ifṁ lc > 0 kg s for saving of energy.

Room model
Since the atrium extends over several floors, the model that depicts the room temperature is discretized into n r = 3 layers of air and concrete masses along the room's height with corresponding temperatures T r,a,k and T r,c,k , k = 1, … , n r . The heat flows acting on these temperatures depend on the ambient temperature T amb , the direct solar irradiation I r,dir and the diffuse solar irradiation I r,diff on the building, as well as the heat exchangeQ fcu,r with the FCUs. Under these assumptions, a grey-box resistance-capacitance (RC) model has been identified for the atrium within an extensive system identification campaign which reads as dT r,a,1 (t) dt = 1 C r,a (− 1 R r,a,c,1 (T r,a,1 (t) − T r,c,1 (t)) + 1 R r,a,amb (T amb (t) − T r,a,1 (t)) + 1 R r,a,a (T r,a,2 (t) − T r,a,1 (t)) − r fcu,r,1Qfcu,r (t) + r i,dir,a,1 I r,dir (t) + r i,diff,a,1 I r,diff (t) +Q n,a,1 (t)), dT r,a,2 (t) dt = 1 C r,a (− 1 R r,a,c,2 (T r,a,2 (t) − T r,c,2 (t)) + 1 R r,a,amb (T amb (t) − T r,a,2 (t)) − 1 R r,a,a (T r,a,2 (t) − T r,a,1 (t)) (T r,a,n r (t) − T r,a,2 (t)) − r fcu,r,2Qfcu,r (t) + r i,dir,a,2 I r,dir (t) + r i,diff,a,2 I r,diff (t) +Q n,a,2 (t)), dT r,a,n r (t) dt = 1 C r,a ( − 1 R r,a,c,n r (T r,a,n r (t) − T r,c,n r (t)) + 1 R r,a,amb (T amb (t) − T r,a,n r (t)) − 1 R r,a,a (T r,a,n r (t) − T r,a,2 (t)) − r fcu,r,n rQ fcu,r (t) + r i,dir,a,n r I r,dir (t) + r i,diff,a,n r I r,diff (t) +Q n,a,n r (t) The model identification was carried out in multiple steps. First, a reduced version of the model was identified that contained only the thermal capacities C r,a , C r,c and resistances R r,a,amb , R r,c,amb , R r,a,c,k , k = 1, … , n r . For this, measurement data obtained during night times at different ambient conditions was utilized, so that the influence of solar irradiation and occupancy on the room temperature was excluded. Also, the FCUs were not operated during that time.
Using the identified capacities and resistances, several ratios r fcu,r,k , k = 1, … , n r were introduced to the model that depict the ratio of the total FCU heat exchangeQ fcu,r for an air mass. Identification of these ratios was carried out using data obtained especially while operating the FCUs during night times.
Afterward, further ratios have been introduced to the model, whereof r i,dir,a,k , r i,diff,a,k , k = 1, … , n r depict the ratios of the direct and diffuse solar irradiation acting on an air mass and r i,dir,c,k , r i,diff,c,k , k = 1, … , n r the ratios acting on a concrete mass. These were identified using data obtained during day times especially on weekends and bank holidays to minimize the influence of occupancy on the measured room temperature.
Finally, heat noise termsQ n,a,k ,Q n,c,k , k = 1, ..., n r are introduced to the model to depict unknown influences from, for example, occupancy and air exchange, of which no data are available. Similar to (15), these are regarded as states which are estimated during operation, while their development is described by dQ n,c,k (t) dt = −Q n,c,k (t) Q n,c,k + ẇQ n,c,k (t), k = 1, ..., n r , witḣQ n,a ,̇Q n,c adequate time constants and ẇQ n,a , ẇQ n,c additive white process noise terms used within MHE and considered zero within MPC.

Operation conditions and constraints
For formulation of the operation conditions and constraints, we first introduce a vector s of n s = 35 slack variables with For control of HVAC systems, it is favorable to define a temperature comfort range instead of tracking a fixed set point temperature, as this provides flexibility for the control algorithm, and with this, further energy saving potential. 39 Using the slack variable ΔT r,a and the room air temperature T r,a,1 , the deviation of the room air temperature from such a comfort range can be formulated as T r,a,min − ΔT r,a (t) ≤ T r,a,1 (t) ≤ T r,a,max + ΔT r,a (t), and using the slack variables s ac,lb and s ac,ub , the operation conditions for the AC mode of the ACM stated in Section 2.4 can be formulated as smoothened vanishing constraints 30,40 as T ac,ht,min − T hts,1 (t) T ac,lt,min − T lts,1 (t) with chosen ac,lb = ac,ub = 0.1 wherein the slack variables can be utilized to relax these constraints if necessary to preserve feasibility when applied within an optimization problem. Further, the constraints (T fpsc (t) − T sc,so )(v ppsc,so − s ppsc,fpsc (t)) ≤ sc,so , with chosen sc,so = 0.1 implement the safety measure for P psc not to operate below v ppsc,so once one of the collector temperatures exceeds T sc,so . Throughṁ it is ensured that the mass flow through the several storage outlets and inlets does not exceed the total mass flow induced by the respective pumps. The upper and lower temperature bounds for the water states in x selected by the selection matrix Z w are defined as soft constraints 41 as and the upper bounds for the solar liquid states in the PSC selected by the selection matrix Z sl are defined as Different from that, the bounds for state p mpsc are given as hard constraints as in with p mpsc,min = 0.1 and p mpsc,max = 0.9, since due to limit switches of the mixing valve, constraint violation here is prevented physically. Finally, an additional upper bound for the temperature state T shx,psc,n shx is introduced as T shx,psc,n shx (t) ≤ T shx,psc,out,max + s shx,psc,out (t) with T shx,psc,out,ma = 85 • C and s shx,psc,out ∈ s sl for further regulation of the temperature of the solar liquid returned from the PSC-side outlet of the SHX to the solar collector array inlets.

Model summary
In total, the system consists of n x = 42 states, n u = 6 continuous controls, n b = 2 binary controls, n c = 5 time-varying parameters, n y = 22 measured quantities, and n w = 11 process noise terms, as given by c ⊤ = [T amb,fc , I vtsc,fc , I fpsc,fc , I r,dir,fc , I r,diff,fc ], and the system model f consisting of Equations (1), (2), (5), (6), and (9)-(51) is a system of explicit ODEs which is directly given in the so-called partial outer convexified form, wherein the system dynamics are given as the sum of a function f 0 and the functions f i , i = 1 … , n b whose activity is determined by the corresponding binary controls b i , of which at most one can be active at a time. Availability of this special structure is important for the solution methods applied in the scope of this work, and can, if not directly given, also be achieved via systematic reformulation, cf. Reference 42.
The system is subject to the constraints (53)-(64) and additional combinatorial constraints regarding the maximum amount of allowed switches and minimum dwell times for the binary controls summarized as In (71), the notation b(⋅) is chosen to indicate that such combinatorial constraints can be coupled over time. The treatment of combinatorial constraints in the scope of this work is explained in more detail later in Section 4.3.

DESCRIPTION AND IMPLEMENTATION OF THE MPC ALGORITHM
In this section, first the MIOCP used within MPC of the system is presented, followed by a description of the algorithm and software used for implementation and solution of the problem. Afterward, the formulation and implementation of the MHE problem utilized within state estimation for the system is presented. The section is concluded by a description of the utilized computational hardware, software, and data sources. The information flow between MPC, MHE, and the plant is illustrated in Figure 7, which shows that the control signals u and b computed by MPC are applied for plant operation, while the measurements y obtained from the plant are used for estimation of the current system statê x via MHE. For a demo implementation of the presented algorithms with reduced dependencies, we refer the reader to Reference 43.

Optimal control problem formulation
With the states x, continuous controls u, binary controls b, time-varying parameters c, and slack variables s introduced in the previous section, the MIOCP reads as F I G U R E 7 Information flow between MPC, MHE, and the plant minimize The objective (72a) of the MIOCP is a Lagrange-type economic objective and contains both the sum of squares and the sum of the slack variables s, as well as the sum of squares and the sum of the continuous controls u, weighted by appropriate diagonal weighting matrices W s ∈ R n s ×n s and W u ∈ R n u ×n u and vectors s ∈ R n s and u ∈ R n u , respectively. ¶ The system dynamics described in Section 3 enter the MIOCP in (72b). The Special Order Set 1 (SOS1) constraint in (72f) ensures that at most one binary control is active at a time. The inequality and combinatorial constraints described in Section 3 enter the MIOCP in (72c). The limits of the continuous controls u are given in (72d) and the binary constraints for b in (72e). Constraint (72g) ensures that the slack variables s only take positive values and constraint (72h) defines that process noise is considered zero within the MIOCP. The initial state constraint is given in (72i).

Description of the MPC algorithm
For solution of OCPs, application of direct methods (first discretize, then optimize approach) and especially direct multiple shooting 44 and direct collocation 45 is favorable. For discretization of the MIOCP (72), we utilize the fact that typically more control interaction is needed during daytime when operating a solar-driven system. Therefore, we set up a time grid for the control horizon with t f = 24 h with time steps of Δt d = 15 min in times between 4:30 and 21:30 and Δt n = 30 min between 21:30 and 4:30 to reduce the number of discrete time points, and with this, the size of the resulting MINLP. Since the step size Δt d might still be to large to successfully control the solar collector temperatures, as the solar irradiation can be subject to fast and large fluctuations (e. g. through movement of clouds) and create need for more frequent updates of P psc and M psc , the first discrete interval of the time grid is further discretized in n t,mpc,s = 10 shorter intervals. This results in a time grid with a total of n t,mpc = 91 intervals. The discretization is conducted using direct collocation with Lagrange polynomials with Radau collocation points of order 3, cf. Reference 34. To avoid the need for regeneration of the discretization in between subsequent rolling-horizon MPC steps due to changing step length in between two time points caused by the nonequidistant time grid, a time transformation is applied to the ODE right hand side so that the time steps can enter the resulting MINLP as parameters, cf. References 30 and 46. For controls and time-varying parameters, we assume that their values only change at the discrete time points. Also, different from the description in Section 3, the variables p mpsc ,ṁ o,hts,b , andṁ i,hts,b enter the discretized MIOCP not as states but directly as controls of the system, and the change of these controls is penalized by including the squared difference of the piece-wise constant control values for subsequent time steps in the MINLP objective.
For solution of the resulting MINLPs on real-time suitable time scales, we apply the combinatorial integral approximation (CIA) decomposition approach 26,27 which has shown good performance regarding solution time and quality in previous studies. 30,47 The idea of the approach is to first solve a relaxed version of the MINLP with dropped binary constraints, which is a nonlinear program (NLP) further referred to as NLP rel , then to approximate the obtained relaxed binary controls q ∈ [0, 1] n b with binary controls p ∈ {0, 1} n b by solution of a CIA problem, 31 and afterward to solve the MINLP again with binary controls b ∶= p fixed, which is again an NLP further referred to as NLP bin , to adjust the states x, continuous controls u, and the other remaining optimization variables to the obtained binary solution, cf. Reference 48.
As explained above, the continuous controls should be updated frequently within this application. However, such frequent updates are not necessary for the binary controls, which renders the solution of the complete MINLP for each short term step unnecessary. Also, despite the fact that the decomposition approach can have a drastically reduced solution time compared to general MINLP solvers, 47 the time spent within the individual solution steps can still be significant, especially for the NLPs. 30 Under these considerations, we introduce the algorithm depicted in Figure 8, which is a parallelized algorithm for mixed-integer nonlinear MPC. After generation of an initial guess using a suitable system simulation f init for the estimated initial statex init of the system, the algorithm is running in two parallel threads. While the thread depicted in the top solves the relaxed MINLP and the CIA problem to obtain updates for the relaxed and approximated binary solution q and p, respectively, the thread depicted on the bottom computes more frequent updates of continuous controls u, states x etc. by repeated solution of NLP bin . After m solutions of NLP bin , the previously used binary controls are replaced with the updated binary approximation p. For the first iteration k = 0 of the algorithm, the binary controls for use within NLP bin are provided by the initial guess simulation.
To compensate the delay introduced by the solution time of the NLPs, a method entitled delay compensation by prediction 41 is applied. Here, the NLPs are not formulated and solved to generate controls for the current discrete time point t k,i at which the computation started, but for an upcoming discrete time point t k,i+1 and applied once this time point is reached. For this, a predictionx + k,i+1 of the state for the next relevant discrete time point is required, which is generated by solving a corresponding initial value problem (IVP) f (⋅) for the ODE model f on the corresponding time interval, considering the current statex k,i of the system and the controls applied during the time interval. Accordingly, NLP rel and the instance of NLP bin solved at the time point t m correspond to the first discrete time point t k+1,0 of the next iteration k + 1 of the algorithm.
Due to the chosen time grid, the scheduled time for one iteration of the total algorithm is Δt d = 15 min during day and Δt n = 30 min during night. Accordingly, the available time for preparation and solution of one instance of NLP bin is Δt d ∕n t,mpc,s = 90 s during day and Δt n ∕n t,mpc,s = 180 s during night, while the available time for preparation and solution of one sequence of NLP rel and CIA is ((n t,mpc,s − 1)∕n t,mpc,s )Δt d = 810 s during day and ((n t,mpc,s − 1)∕n t,mpc,s )Δt n = 1620 s during night. Within this study, the preparation time reserved for each problem instance for, for example, data preparation and simulation, is approximately 6 s, the remaining time is available for solution of the problem instance.

Implementation of the MPC algorithm
The algorithm implementation is based on Python, while the separate threads are realized using the Python multiprocessing module 49 and data exchange and synchronization in between processes are realized using blocking queues. 49 The MINLP is implemented using CasADi 50 via its Python interface. For this application, preliminary tests showed that using CasADi's SX-type variables results in reduced NLP solution times in comparison to using MX-type variables, since less time is spent for evaluating NLP functions and derivatives. However, utilization of the spline interpolation functionalities 35 to depict the solar collector mass flows shown in Section 3.1 requires the use of MX-type variables. But nevertheless, the code generation functionality of CasADi can be used to a priori generate C code for the NLP that contains all functions required by the solver, which can be compiled prior to starting the MPC algorithm and the resulting object later be used (and reused) with the NLP solver. For this application, the generated code was compiled using Clang 51 with optimization level 1 and the resulting object was able to yield function evaluation times comparable to using SX-type variables.
The NLPs within the MPC algorithm are solved using Ipopt 52 with linear solver MA57. 53 Ipopt is configured to use its adaptive barrier parameter update strategy. The CIA problem is solved using pycombina. 30 The initial guess simulation is implemented in Modelica and conducted using OpenModelica 54 via its Python interface OMPython. 55 The maximum switching and dwell time constraints for the ACM operation modes are considered only in the CIA problem stage of the MPC algorithm and not in the NLP rel stage to avoid further increase of the size of NLP rel . In particular, the AC mode is required to remain active for at least 1 h after activation and inactive for at least 0.5 h after deactivation, while the FC mode is required to remain both active after activation and inactive after deactivation for at least 0.25 h. The maximum amount of allowed switches is determined heuristically from the relaxed solution of the preceding NLP rel stage, however, with a minimum number of three allowed switches per mode.

MHE problem formulation and implementation
The relation of states x and associated measurements y is depicted by the measurement function Φ(⋅). We assume that the states T hts of the HTS and T lts of the LTS are associated with the temperature sensors T hts and T lts , and the temperature states at the outlet of the solar collector arrays T vtsc,s and T fpsc,s with the sensors T vtsc and T fpsc , respectively. The states T shx,psc,1 , T shx,psc,n shx , T shx,ssc,1 , and T shx,ssc,n shx are determined via T shx,psc,i , T shx,psc,o , T shx,ssc,i , and T shx,ssc,o . The state of the water side of the FCUs is associated with T fcu,o . The air temperature states T r,a of the room are determined from a combination of air room sensors at different heights and locations, similar to the concrete temperature states T r,c which are determined from the sensors attached to the concrete masses.
The current value of the auxiliary state Δ T amb is determined via the corresponding forecast T amb,fc and the ambient temperature sensor T amb . Similarly, the states Δ I vtsc and Δ I fpsc are determined via the forecasts I vtsc,fc and I fpsc,fc and the pyranometer I g with measurements converted to the orientation of the corresponding collector array's surface.
Using a time grid with n t,mhe = 20 equidistant time steps over the past 20 minutes, the MHE problem for estimation of the current system statex = x n t,mhe is formulated as with Φ the measurement function and y the corresponding measurements at the discrete time points with inverse covariance matrix Σ −1 y , w the additive white process noise with inverse covariance matrix Σ −1 w , and ||x 0 − x arr || 2 Σ −1 x arr the arrival cost with inverse covariance matrix Σ −1 x arr . Function f (⋅) is the solution of an IVP for the ODE model f for a discrete time interval [t k , t k+1 ].
Similar to the MIOCP, the MHE problem is implemented using direct collocation and code generation in CasADi. As (73) contains no discrete optimization variables, an application of the decomposition approach is not required and the arising NLPs are directly solved using Ipopt and MA57. Arrival cost updates are conducted using a smoothened extended Kalman filter (EKF), cf. Reference 56. The problem is solved every 60 s in an additional third thread, separate from the MPC loop.

4.5
Computational hardware, software, and data sources

EXPERIMENTAL MIXED-INTEGER NONLINEAR MPC OPERATION
In the upcoming sections, results of experimental operation of the system using mixed-integer nonlinear MPC are presented. In the following, first the experimental setup is described, followed by a presentation and discussion of the experimental results. Finally, the solution times of the algorithm are presented.

Setup of the experimental study
To investigate the applicability of mixed-integer nonlinear MPC for the solar thermal climate system, experimental operation of the system has been realized. The results presented in the following were obtained from an operation test on August 1, 2019 from 05:00 to 19:10 UTC. For that day, weather forecasts indicated high solar irradiation and relatively high ambient temperature for a summer day. Especially due to the high availability of driving energy, it was expected that the day would provide beneficial conditions to evaluate the capabilities of the MPC for operating the system at high temperatures and for prevention of overheating. For running the experiment, the MPC was activated in the morning after a night where the climate system was not operating. Due to a longer hot weather period prior to the experiment, both the air and concrete temperature of the building were relatively high, partially exceeding 25 • C. Due to that and the high availability of driving energy, it was decided to set the room temperature comfort range for the experiment to a deliberately low range of [T r,a,min , T r,a,max ] = [21 • C, 22 • C], as even though this room temperatures would not be achieved due to limitations in cooling power, such setting would cause the MPC to aim at maximum cooling power generation. The outdoor blinds of the glass roof of the building were closed during the whole day. The heat input caused by people in the building was very limited as the experiment was conducted during the semester break.

5.2
Results and discussion of the experimental study Figure 9 shows the results of the experimental MPC operation. From top to bottom, the plot shows the measured ambient conditions, solar collector temperatures, HTS temperatures, LTS temperatures, and room air and concrete temperatures measured over the considered horizon and corresponding, influencing controls. In the upcoming sections, a detailed explanation and interpretation of the system behavior is given. The ambient conditions measured during the day of the experiment are shown in the top plot of Figure 9. The solar irradiation throughout the day has been measured using the pyranometer on the roof. While the irradiation profile I vtsc on the VTSC was measured directly, the irradiation I fpsc on the FPSC mounted at 45 • angle is computed from the pyranometer measurements using pysolar. 58 The depicted measurements indicate that, after a partly cloudy morning, the irradiation profile was very stable during noon with high irradiation on the collector surfaces, apart from a few disturbances in the afternoon, until the evening again was partly cloudy. This confirms the expected high availability of solar energy during the day at predominantly stable conditions.
The ambient temperature T amb measured on the roof of the building rose steadily from approximately 18 • C in the morning and peaked in the late afternoon at 35 • C. Later, T amb decreased again to approximately 25 • C at the end of the experiment. The profile shows the expected higher ambient temperature during the day, as well as free cooling potential that could be utilized by the system at the beginning of the day.

Solar circuits
The second plot from top in Figure 9 shows the measurements for the temperature sensors T vtsc and T fpsc at the outlet of the VTSC and FPSC, respectively. It can be observed that despite the partly fluctuating but high solar irradiation, the MPC is able to achieve high, however slightly fluctuating, collector temperatures close to their defined maximum temperature throughout the major part of the day, which is a desirable result. Prior to approximately 9:00, an interesting observation can be made where the controller utilizes the solar mixing valve to keep the mass flow through the FPSC low while increasing the flow through the VTSC. Even though T fpsc,s shows a higher temperature than T vtsc,s , it appears to be more favorable not to realize a higher flow through the FPSC and to avoid temperature alignment of the arrays, which would be the aim of a conventional control approach. As both the area and the heat loss coefficient of the FPSC are bigger compared to the VTSC, the controls identified by the MPC increase the heat dissipation of the system to the environment. One plausible interpretation is that the MPC aims to reduce heat storage during that time, and with this, to prevent overheating of the system later during the day, which would be particularly remarkable.
During a phase of active heating of the collector fields between 09:30 and 10:00, the temperature exceeds the specified maximum by approx. 5 K. However, this temperature is still below the limits at which the EC mode of the ACM is activated by the PLC. In fact, activation of the EC mode was avoided throughout the entire day, and the safety constraints (56)- (58) which ensure that P psc is always in full operation once one collector temperature exceeds 65 • C are always fulfilled.
The MPC performance could possibly be further improved by introduction of a more detailed model of the pipe that connects the SHX and the collector fields. Though the volume of one pipe distance from the cellar to the roof contains approximately 0.2 m 3 of solar liquid, currently only one temperature state is used for modeling of this distance. Therefore, a more fine-grained depiction of the temperature distribution along the tube using more temperature states could probably improve both temperature constraint satisfaction and reduce the fluctuation of collector temperatures when operating close to their maximum value. However, this increase in the number of states would also increase the solution time for the MINLP, while reliable estimation of the corresponding state values would probably require installation of further temperature sensors.

High-temperature storage
The middle plot in Figure 9 shows the measurements for the n hts = 4 temperature sensors T hts,k , k = 1, … , n hts of the stratified HTS. It can be observed how the MPC actively utilizes the several inlets and outlets of the storage within its control decisions, which is detailed in the following. Prior to 11:00, the MPC uses mainly the middle outlet O hts,m and middle inlet I hts,m of the HTS. Doing this, it utilizes only the upper part of the stratified HTS for heat exchange with the PSC via SHX, which results in a rapid heating of this part of the storage as illustrated by the sensors T hts,1 and T hts,2 . Due to that, high HT temperatures for ACM operation can be achieved earlier as if the total storage would be heated up, which results in both a longer and more efficient operation of the ACM.
Past 11:00, the MPC first starts to utilize the lower outlet O hts,b . By mixing in water from the bottom region of the HTS, the temperature of the solar liquid at the SHX outlet can be reduced more as if only water from the upper region of the HTS was supported to the SSC. This is actively considered to prevent the solar collectors from overheating.
Later during the day and after the peak of solar energy availability, the MPC often fully utilizes the lower outlet O hts,b and lower inlet I hts,b . As no storage capacity needs to be held available for reduction of the collector temperatures, the total HTS can be used at that part of the day so that more of the heat energy remaining in the PSC can be extracted.
These observations are particularly remarkable, as they exemplify the advantages of MPC for situational and individual control decisions. While it was not favorable in the morning to use the total HTS but to keep the lower part of the storage at lower temperatures, it is favorable in the afternoon to extract as much energy as possible from the PSC and store it in the HTS. As the optimality of this decision however depends strongly on the expected ambient conditions and differs for other scenarios, it is hard to achieve comparable behavior using a conventional controller based exclusively on current measurements.

Adsorption cooling machine
As shown in the second plot from the bottom in Figure 9, the ACM is operated in FC mode during the morning when ambient temperature is still comparatively low and it can be observed that the temperature of the LTS is reduced. At the same time, P lc and the FCUs are operated to decrease the room air temperature, as shown in the bottom plot of Figure 9. Once a sufficiently high temperature in the top of the HTS is reached, the ACM is started and operates in AC mode at comparatively high LT storage temperature and still low ambient temperature, which results in high COP and cooling power of the ACM. This results in a rapid temperature decrease in the LTS. During the day, operation in AC mode continues and the LTS temperature decreases further until approximately 12:00. After that, rising ambient temperatures and later also decreasing HTS temperatures result in a slight but steady increase of LTS temperatures until the machine is turned off at around 18:30. After that time, the MPC decides to use the remaining cooling energy stored in the LTS, before continuing in FC mode at the end of the experiment when the LTS temperature meets the ambient temperature level again.
Apart from its purpose to provide cooling power, operation of the ACM in AC mode is important for the system functionality itself, as without the machine acting as a heat consumer, emergency cooling of the HTS via the RT or a shutdown of the system due to gasification of the solar liquid might become necessary at some point. Accordingly, surplus of heat during the day and corresponding negative effects can be avoided if the activation time of the ACM is chosen based on the expected solar irradiation over the day instead of fixed HTS temperature set points. In contrast to using MPC, this can hardly be achieved using conventional, set-point and PID-based controllers.
This can be further exemplified on the controllers choice for ACM operation at the end of the day. Though the HTS temperature would still be sufficiently high to drive the ACM in AC mode after 18:30, the MPC decides to utilize FC also before all driving energy in the HTS is used up. In between, ACM operation is paused so that the cooling energy stored in the LTS can be utilized first, as the temperature level for FC utilization is rather high.
Doing this, driving energy is retained in the HTS for an earlier AC start in the upcoming day. Such an earlier start has been scheduled for approximately 5:00 of the upcoming day, as shown in Figure 10 which depicts the predicted system operation for the upcoming 24 h at the end time point of the experiment. While the controller aims to use FC during the night where ambient temperatures are lower, the AC is activated early on the upcoming day, which results in favorable operation conditions for the ACM while generating storage capacity for the solar energy expected during the day.

Low-temperature circuits and room installations
As expected from the setup of the experiment, the pump P lc and with it the FCUs are operated permanently, due to the relatively high temperatures in the building compared to T r,a,min and T r,a,max , which is depicted in the bottom plot in Figure 9.
In the beginning of the day, a major temperature drop of the room air temperatures T r,a,1 and T r,a,2 can be observed. This is, however, not solely caused by the operation of the climate system alone, but due to major ambient air exchange caused by multiple open windows in the building opened by the faculty staff in the beginning of the day, which were closed again once the ambient temperature achieved a higher level than the room air temperature. This illustrates the high influence of building usage on system operation.
As soon as the windows got closed at around 7:00, the room temperature started to increase again since the concrete temperatures T r,c,1 and T r,c,2 were still high, which illustrates the strong influence of the heat stored in the concrete mass of the building on the air temperature development. Due to that, the climate system is as expected not capable of driving the room temperatures to the specified comfort region, as the system setup is not capable of providing the required cooling power. However, the climate system is able to prevent further heating of the building and to compensate most of the occurring heat loads, illustrated by the stagnating room temperature development.

Runtime of the algorithm
An overview of the runtime of the individual steps of the MPC algorithm and the return status of the solution processes are given in Figure 11. It can be observed that solving NLP bin takes typically less time than NLP rel , which was assumed during algorithm design. More than 95 % of the 416 solutions of NLP bin finished prior to the timeout of approximately 84 s. For all of these finished instances and for all instances of NLP rel , often optimal or at least acceptable solutions have been obtained, so that no problem has been classified infeasible. While the runtime of the CIA is typically much lower than those of the NLPs, several outliers have been registered with an algorithm runtime of up to 75 s. However, also the maximum runtime for both NLP rel and CIA are still far below the maximum possible runtime of approximately 15 min.

CONCLUSIONS AND FUTURE WORK
In this work, we presented a real-life experiment in which we demonstrate the real-time applicability of nonlinear mixed-integer MPC for a solar thermal climate system. For this, we first described a control-oriented nonlinear model of the system. Then, we introduced a parallelized algorithm for mixed-integer nonlinear MPC including a corresponding implementation for the system. Afterward, we presented and discussed the results of real-time experimental operation of the physical system. We have shown that the controller was able to successfully operate the system and actively utilized ACM, pumps, and mixing valves within situational and individual control decisions. The advantages of mixed-integer nonlinear MPC for control of the solar collector fields, active utilization of storage temperature stratifications, and consideration of ACM operation efficiency were highlighted. Finally, we have shown that the solution times of the utilized algorithm and implementations were suitable for real-time control of the system. Future work should focus on extended tests for MPC operation, including operation on longer time periods and during different weather scenarios and ambient conditions. Improvement of the system model and systematic tuning of problem and solver parameters could yield further control performance gains. In this regard, comparisons also to other control strategies, such as conventional controllers, linear(ized) MPC strategies, or strategies based on meta-heuristics, could be considered. Also, an extension toward a stochastic MPC approach could be considered, where different forecast variants for upcoming operation conditions could be taken into account simultaneously. Further, if a continuous MPC operation of the system should be facilitated also during cold seasons, the winter operation mode of the system where collected solar heat is used for heating support of the building needs to be included in the system model. In this context, an extension of the current plant configuration by installation of a reversible heat pump could be considered, as an extension of the plant by such an additional switched component could not only yield the possibility for experimental nonlinear mixed-integer MPC applications also during cold seasons, but, however, also allow for a higher number of possible operation modes during summer operation and an increased impact of the cooling system on the building temperatures.