A renewable heat plant Modelica library for dynamic optimization with Optimica

Almost half of the energy consumed globally is under the form of heat, produced mainly through fossil fuels. Switching to using renewable energy instead is a real challenge. Combining renewable thermal energy with thermal storage is a complex system to operate. To harness the full potential of thermal plants, advanced control strategies need to be implemented. Dynamic real-time optimization (DRTO) seems promising to fine tune controller setpoints of plants. The goal of our study is to ultimately enable DRTO by using Optimica because of its ease of use and Modelica’s modularity. This paper presents a Modelica library developed to first perform offline dynamic optimization with Optimica, and would ultimately be used in a DRTO strategy. The library enables to model a renewable thermal plant composed of solar thermals, heat pumps and thermal storages. The model of each subcomponent has been validated. Initial dynamic optimizations of plant operation give promising results.


Introduction
Completely replacing fossil fuels for heat production by renewable energy requires combining several heat sources.The integration of thermal renewable energy in conversion/storage/distribution systems, in particular, solar thermal energy, faces several obstacles.Renewable thermal systems have varying thermal inertias and are prone to environmental disturbances.These mixed sources have to be well controlled to maximize performance and competitiveness of the global system.
Combining both a mix of renewable thermal energy sources and large thermal energy storage is quite innovative and has never been implemented in France in large scale heat plants.Such innovative systems have successfully been built and operated in a few Northern European countries.The main learning highlighted by those first plants is the need to fine tune system setpoints in real time to both maximize energy output and minimize operational cost.In China, a study demonstrated that in specific cases of large-scale solar heating systems integrated with water-to-water heat pumps and pit storage, a 16% decrease of the leverage cost of heat (LCOH) can be achieved by operation optimization (Zhang et al. 2023).
Since renewable heat plants are sensitive to external variation such as changing weather, electricity cost, and heat demand, real-time optimization (RTO) is a method that seems perfectly adapted to them.Although it is widely used in the field of chemical engineering, using RTO in the field of energy is quite recent.Adding dynamic optimization would allow to handle the thermal inertias of the system.Then, Dynamic Real-Time Optimization (DRTO) seems particularly adapted to the management of thermal installations combining several renewable production sources.
In our global study, DRTO combined with Non-Linear Programming (NLP) will be used to optimize the operation of a multi-energy heat plant.The optimization will be performed using Modelica language and the Optimica Compiler Toolkit (OCT) from Modelon.
The objective of this optimization is to maximize the heat produced by the plant, while minimizing the operational cost (mainly electric consumption).Weather, electricity, and heat demand forecasts are included in the model.The optimization variables are temperature and power setpoints for each heat source.
The modelling work and initial dynamic optimizations are presented in this paper.DRTO theory and results won´t be discussed in detail in this paper.

Literature review
JModelica is an open-source platform (Åkesson et al. 2009) for numerically solving large-scale dynamic optimization problems of Modelica models.This tool evolved during the years, in particular by including CasADi.The actual framework JModelica.org(Magnusson and Åkesson 2015) became Optimica Compiler Toolkit (OCT) under a Modelon license in 2020.
JModelica.org is used since its development in many research works, with some in heat production optimization.Runvik et al. (2015) developed a short-term planning optimization for a district heating system solved in two steps, one MILP (Mixed Integer Linear Programming) and one NLP using JModelica.More recently, Rohde et al. (2020) used JModelica to dynamically optimize control setpoints for an integrated heating and cooling system, including heat pump and solar thermal, with thermal energy storages.Other works used Modelica language to model the heat system but chose other ways to optimize the heat production.Liu et al. (2018) used Dymola to model the system composed of CO2 heat pump coupled with hot and cold thermal storages, and genetic algorithm to maximize the efficiency of the system.
Dynamic optimization of heat plants including solar thermal and thermal storage using NLP has been also performed by Scolan et al. (2020).They optimized the design and the control of the plant using the optimization software GAMS for modeling and optimization.More recently, Untrau et al. (2023) proposed a DRTO of solar thermal plant including thermal storage also using GAMS, combined with a detailed simulation model of the real plant in Matlab.

Purpose of this work
Dynamic optimization with Optimica requires dynamic models of the thermal plant formulated with continuous equations; discontinuities would be hardly interpreted by the time discretization method (collocation method) used by OCT.The objective of this research work is to develop a library in Modelica compatible with Optimica to model a thermal plant composed of a solar thermal field, a heat pump, and a thermal energy storage as shown in Figure 1.Modelica libraries already exist and model some of the components of thermal plants, but either are not compatible with Optimica or have too complex fluid modeling, which hinders convergence of the optimization.
This library allows to model and dynamically optimize several plant layouts using the modularity of Modelica.
First, we will present the library and the models used for each component of the thermal plant, then we will discuss the simulation results, and finally present the optimization methodology and the first optimization results.

Library
The library is developed in Modelica language using Modelon Impact software.
Each main component of the library and associated controller will be detailed in the following sections.The controllers will be used to initialize the optimization problem as detailed in section 4.

Hydraulic representation and interfaces
Two different fluids will be used in this library, water on storage and process side, and a glycol-based water solution on solar side.Glycol-based water solution is used to avoid frosting inside the solar field during winter (in European countries).
To simplify modeling, the only hydraulic parameters considered in this library are temperature and mass flow with constant fluid properties.Density  and specific heat capacity  only vary of few percents over the operating temperature range (40 to 90°C).The main impact of those hypotheses is on the viscosity of the glycol-based water solution, which is varying a lot over the operating temperature range.To further simplify and ensure optimization convergence, pressure loss is not modeled here.This hypothesis leads to inexact electricity consumption of the pump on solar side.However, as the main electricity consumption is from the heat pump, around 15 times more than the solar pump in the plant studied (based on Newheat's proprietary plant design tool), it will only lead to a small error in the overall electricity consumption.
The Modelica connectors used in this library are propagating only temperature and mass flow of the fluid.The mass flow is generated by the pump and is then propagated through each component of the hydraulic loop.

Solar Field (SF)
The solar thermal collector is the equipment used to transform solar radiation into heat.The collectors modeled here are Flat Plate Collectors (FPC).They provide heat at low temperature (below 100°C).Figure 2 shows the main components of the FPC and how it works.The solar field is usually composed by several parallel loops of collectors, and each loop is composed by several collectors in series.This layout is described in Figure 3.

Model
The model follows the standard ISO/FDIS 9806:2017 (International Standard 2017), using the quasi-dynamic equation of solar thermal collectors.This standard determines characteristic parameters for each collector ( ,  ,  and  in our model).The whole solar field is modeled as an equivalent panel with  equal to the whole solar field area, assuming a uniform distribution between each panel loop of the solar field, whose power is  ̇.
The model in Figure 4 is described by equations (1)(2)(3), with  the global irradiance received by the collectors, ̇ the mass flow inside the solar field,  the ambient temperature, and  ,  ,  the inlet, outlet and mean temperatures of the solar field.

Control
In most solar thermal plants, the outlet temperature of the solar field is controlled to follow a temperature setpoint which could vary within the year.This temperature is controlled thanks to the mass flow inside the panels.Thus, the mass flow to be provided by the solar pump is calculated in the simulation to get the temperature setpoint at the outlet of the solar field by solving equations (1)(2)(3) without the differential term, as standardly done in solar thermal plants.

Plate Heat Exchanger (PHEX)
A heat exchanger is needed to transfer heat from the solar field loop (filled with glycol-based water) to the storage and supply loop (filled with water).PHEX are used in most solar thermal plants operating at low temperature (below 100°C).PHEX in Figure 5 is modelled with constant efficiency ( between 0 and 1, a value of 0.9 is used to match Newheat's plants data) in equation ( 4) and is used to have the heat capacity flow ratio equals to 1 in equation( 5).

Model
The energy conservation gives those latest equations to get the outlet temperatures: ̇= ̇  ( , −  , )

Control
The mass flow on side 2 will be computed to always respect the equality of the two heat capacity flows as in equation ( 5).

Tank thermal energy storage (TTES)
TTES is an essential element of solar thermal plants, it allows desynchronizing solar heat production and process heat demand, as the heat production depends on a fluctuating solar irradiance.During the charge phase the solar field provides heat to the tank (cold water from the bottom is warmed up by the solar field before coming back to the top of the tank).Conversely, during the discharge phase the tank provides heat to the process (hot water from the top of the tank transfers heat to cold water coming from the process).
The tank used in this work is an insulated water tank.

Model
Different tank models are available in the literature (Dumont et al. 2016).The model used in this work is a  The equation governing the stratified thermal model for conduction and convection is the energy equation ( 8) (Hawlader et al. 1988).
Where (, ) is the storage fluid temperature dependent on height  and time ,  the flow velocity,  the thermal diffusivity,  the diffusion coefficient due to turbulent mixing caused by buoyancy effect,  the overall heat transfer coefficient,  the tank cross-sectional area, and associated  perimeter.
The tank is discretized into n nodes from equation (8) using a finite difference method (Scolan et al., 2020).Each layer is exchanging heat with its neighboring, spatial derivatives being then expressed by second order finite differences.
In the studied heat plant, solar heat and heat pump could lead to a temperature inversion in the tank (an upper layer which is colder than a lower one).It generates a mixing flows called buoyancy effect or plume entrainment (Hawlader et al. 1988).A model of buoyancy was implemented in Modelica IBPSA library (IBPSA 2013).Based on IBPSA model we implemented an equivalent model computing ̇ , as a mixing mass flow between layer i and layer i-1 following equations (9)(10)(11).
where  is a time constant for mixing.

Water-to-water Heat Pump
The first heat pump used in this work is a water-towater compression heat pump.
Heat pump could be used in several ways in a thermal plant.In this work, as shown Figure 1, the solar heat stored in the tank is used as cold source (evaporator side); the condenser side is warming up the cold water coming from the industrial process.To model the heat pump, we need a simplified model because detailed one would be too complex for optimization and prevent convergence.We decided to neither model the refrigerant fluid inside the heat pump nor its dynamics.

Model
The heat pump is modeled using Coefficient Of Performance (COP) modeling with evaporator and condenser temperatures mapped on a datasheet.The heat pump can be used at partial load.We first assume that COP does not vary with the load of the heat pump.
• Energy conservation: • COP modeling: The COP is defined as: We model the COP as a polynomial expression of temperature lift between evaporator and condenser.This method is used by Fischer et al. (2017) and Ruhnau et al. (2019).
Coefficients  ,  and  are determined using the datasheet of the heat pump at full load.

• Compressor modeling:
The compressor is modeled with a constant efficiency  and its available power  ̇ , depends on the temperature of the fluid inside him.Fischer et al. (2017) propose to consider only evaporator temperature to determine the available power.
Coefficients  and  are determined using the datasheet of the heat pump at full load.
Partial load y is defined as below: The electric power of the compressor is considered as an expenditure in the heat plant operation.

Control
To model the COP, we used specific working conditions.All operating points are usually given with a constant temperature difference between the inlet and outlet of evaporator and condenser.We also want to be able to control the evaporator and condenser temperatures.
To control the temperature difference, we need to add a pump, and to control the temperature, we need to add an ideal three-way valve, on each side as shown Figure 8.The valve position and the mass flow are ideally controlled to get the temperature setpoints (this means specification equations are added to the system to be solved).The heat pump is turned on when the cold source is warm enough (above the minimal temperature accepted at evaporator) and the warm source is cold enough (below the maximal temperature accepted at condenser). , and y will be computed to follow heat demand temperature and mass flow.

Pump
The electricity consumption of the pumps will be considered as an expenditure in the heat plant operation.Pumps and heat pump are the main operating cost of the plant.

Model
The electricity consumption of a pump is determined by equation ( 21) where  , is the electric power measured on the real pump at its maximal mass flow ̇ .

Control
The pump works as a mass flow generator, it directly provides the mass flow desired.It does not need specific control.

Simulation results
In previous section, all the models have been detailed.They need now to be validated with experimental results or datasheet, before simulating a full plant.

Validation
All models (except heat pump) are compared to measurement from solar thermal plants operated by Newheat.For solar field and heat exchangers we used data from Solthermalt, a plant owned by Kyotherm providing heat to a malthouse (Newheat 2020) in Issoudun, France.Tank model has been compared to the one build in Lactosol project, a solar plant providing heat to a whey powder production site (Newheat 2023) in Verdun, France.
Following sections show the result of comparison between the model and the measurements.
In figures below, time of 0.0 day corresponds to 12AM, and 0.5 day to 12PM.The model fits quite well to the measurements except during night periods.The inlet and outlet temperature sensors are in an insulated pipe inside a building; therefore, they do not decrease to the ambient temperature.In the model, the solar field mean temperature decreases to ambient temperature when there is no irradiation.As ambient temperature is around 10°C and measured inlet around 30°C (instead of 10°C if pipes were not insulated), the outlet is logically computed to -10°C to keep the mean of the inlet and outlet temperatures at 10°C.Since there is no mass flow during this period, heat produced by the solar field is not impacted.Same as previous section, the model fits well to the measurements except when the solar and storage pumps are off.

Tank thermal energy storage
Measurements are coming from Lactosol plant.A 3000 m 3 insulated water tank is installed in this solar plant to buffer the heat provided by the solar field and the consumption of the whey powder production site.This 12-meter-high tank is instrumented with 12 temperature sensors (about one every meter).
Two models are simulated and presented, one with 12 layers in Figure 12 and one with 60 layers in Figure 13.Temperature sensors are compared to the temperature at equivalent height in the simulation.
The measurement period is composed of a 9-hour charge phase and an 8-hour discharge phase separated by a 7-hour standby phase.
On the one hand, simulation fits better to measurement with 60 layers than with 12.A higher number of layers in the tank allows indeed higher temperature gradients which are needed to represent the thermocline zone inside the tank.But on the other hand, each layer is adding a state variable to the model and then increasing the size of the optimization problem.The number of layers will have to be selected carefully to model with a satisfactory accuracy without adding too much complexity to the optimization problem.For optimization we decided to set the number of layers at 10 to reduce the size of the optimization problem.

Heat pump
Figure 14 below compares the COP at full load between the model and the datasheet of an industrial heat pump (WWHS ER3b) made by Ochsner Energie Technik.Temperatures are given for the secondary medium (water) and not the primary medium (refrigerant).
The model fits well the datasheet for mean temperatures but leads to a 6.6% error for cold evaporator temperature.The model overestimates a bit the COP, which results in an underestimation of the electricity consumption.

Full plant simulation
The simulated plant is a virtual plant composed of a solar field of 14248 m², a plate heat exchanger with a constant efficiency of 0.9, a 3000 m 3 tank, and a 3.1 MW heat pump.The layout of the plant is described in Figure 1.The heat pump is connected to the tank on the evaporator side while the condenser side is connected to the process.
The heat consumer is represented with a constant heat demand (constant return temperature, and constant mass flow and temperature setpoints).The plant is controlled by the expert rules defined in section 2.
The simulation runs over two spring days, the first one is cloudy and the second one sunny; Figure 15 represents the global irradiation of those two days.An optimization of the control of this plant will be performed on the same days in section 4.3.Inlet and outlet temperatures of the solar field are presented in Figure 16.The setpoint given to the solar field is 53°C.The outlet of the solar field is following well the setpoint except for the second day where the irradiation is too strong to limit the outlet temperature (because the solar pump is at its maximum speed as it can be seen in Figure 17).The tank is discretized in 10 layers (layer 1 is the top layer and layer 10 is the bottom layer).The first day is too cloudy to fill the tank completely, while the second day allows filling the tank at higher temperature (Figure 18).Evaporator outlet temperature is warmer than tank bottom temperature during the second day, which homogenizes the temperatures of layers from 3 to 10.Heat pump temperatures are shown in Figure 19.Heat pump is turned on when the top tank layer temperature exceeds the minimal temperature accepted at evaporator inlet.We can see in Figure 17 when the heat pump is on (supplied mass flow at 15kg/s).Finally, the evaporator temperature is controlled to be as high as possible while staying below its maximum (55°C).
The plant behaves as expected with validated models.We developed and validated a library to simulate renewable thermal plants.The next objective is now to optimize control variables and see how behave the models with optimization solver.

Optimization
Initial optimization presented in this paper is an offline optimization.It means that the optimizer is not yet connected to the real plant (or a highly detailed model).The goal is to try to optimize the plant in typical days and to define the ideal optimization sequences.Once the offline optimization is working well enough, it will be plugged to the real plant to try real-time optimization (optimization launched every hour considering the changes in the system states and the updated forecasts).

Tools
The tool used for optimization in this research work is Optimica Compiler Toolkit (OCT) under Modelon license for academic and commercial use.This tool is coming from JModelica.orgwhich became OCT since 2020.Magnusson and Åkesson (2015) presented how JModelica.org is working.OCT workflow is described in Figure 20.
The first step is to describe the continuous models of the system to optimize in Modelica language.Optimica language (which is an extension of Modelica language) will be then used to describe the constraints and objectives of the optimization problem.
Modelica and Optimica models are transformed into optimization problem in the form of DAEs (Differential-Algebraic system of Equations), which is sent to CasADi, before being discretized via orthogonal collocation.Finally, the discretized optimization problem is sent to the solver (IPOPT) which will optimize the control variables of the system.

Methodology
The objective of optimization is to maximize heat supplied by the thermal plant, while minimizing operational costs.In this first optimization, this equates to maximizing global profit of the plant, electricity price assumed to be fixed.The control variables to be optimized are the outlet temperature of the solar field, evaporator outlet temperature, and the activation of the heat pump.
The optimization is composed of several stages detailed in Figure 21: • An initial simulation of the plant is done with a control strategy, which would be implemented in a real plant (expert rules detailed in section 2).This simulation provides the optimizer an acceptable solution, which is the starting point of the optimizer.
• Several optimizations are done successively releasing degrees of freedom.Each result is initializing the next optimization.Last optimization result (with all degrees of freedom together) is then used to get optimal trajectories for several setpoints of the plant.This iterative approach was chosen to improve convergence of the optimization, otherwise too complex to be solved directly.Linear solver (MA57)

Result
• Finally, to perform DRTO, the controller of the real plant (or a highly detailed model) will use the optimal setpoints to operate the plant exposed to real disturbances.This part is not described in this paper.

First optimization results
In this section is presented the first optimization results using the developed Modelica library.The plant optimized is the one presented in section 3.2.Each figure below compares the standard control strategy used in section 3.2 (named sim) and optimized control (named optim).
The optimization variables are the outlet temperature of the solar field, evaporator temperature, and the activation of the heat pump.
One important thing to remind is that the efficiency of the heat pump (COP) is increasing with evaporator temperature (Figure 14). Figure 22 shows the mass flow supplied to the process (coming from the heat pump) shifting from the first day (around 0.5 day) to the first night (around 1.0 day), knowing that heat demand of the process is considered constant.This suggests that the heat pump does not turn on as soon as the tank is filled with enough hot water but rather is led to wait for the night.It allows the tank to rise in temperature during the first day (Figure 23). Figure 24 shows that the first hours of the night the heat pump is turning on with a higher evaporator temperature (around 42°C) than the rest of the night (around 34°C), which allows to have a better COP the first hours of the night.Figure 25 shows that the solar field outlet temperature is different for the two days.The first day, it is lower than the simulation, it leads to less heat losses in the solar field.But the second day this temperature is much higher.This could be explained by the fact that at the end of the second day the tank is full enough to provide heat continuously to the heat pump in both cases: simulation and optimization (Figure 23).Then, the higher is the solar field temperature, the higher the temperature could be at the evaporator, and thus the COP.It also means that, the optimizer does not care about what could happen after the optimization time horizon: we could imagine a third cloudy day which would need to have a tank more filled (with lower temperature) at the end of the second day than the optimization result.This example shows that the result of the optimization may depend on the considered time horizon.Further, to prevent the optimizer to empty the tank at the end of optimization and closer emulate the behavior of a real system, tank state could be constrained to approach the final value obtained through the initial simulation.

Conclusion and outlook
To conclude, we developed and validated a library which is compatible with Optimica and models a renewable thermal plant.Then we obtained interesting dynamic optimization results of the overall system.
However, the dynamic optimization result is not yet satisfactory since it does not consider the necessity to be able to provide heat after the optimization end time.It could lead to non-optimal results in the real plant, even if only the first hours of the optimization result would be sent to the real plant controller as the optimization will be updated every hour.It also points out the need to perform an offline dynamic optimization based on forecasts on a long enough time horizon before starting DRTO which will correct the control variables trajectories taking into account real disturbances.The final state of the tank could be considered in the optimization objective to get a better result.
We could also consider stratification indicators into the optimization objective, because stratification inside the tank is affecting a lot the efficiency of the storages.Some research works developed indicators to quantify the quality of stratification in thermal storages.
Equipment modeling could also be improved.Efficiency of the PHEX is considered constant.An operating point tabular efficiency could make PHEX model more accurate, with minimal added complexity.A slightly more complex model of COP could decrease error of the COP observed of the heat pump.The current storage model loses accuracy due to low discretization.Other storage models could be considered to either reduce the number of state variables or improve accuracy.
Once the offline optimization will be robust enough, this work will be extended to real-time optimization of the thermal plant.A highly detailed model describing the plant will be needed, so that the optimizer can be run and finetune the plant setpoints in real-time.

Figure 5 .
Figure 5. Plate Heat Exchanger model (icon from DLR ThermoFluid Stream Library) Session 1-C: Applications of Modelica for optimization and optimal control 1 DOI 10.3384/ecp20495 Proceedings of the Modelica Conference 2023 October 9-11, 2023, Aachen, Germany one-dimensional model called Multi-Node.Tank is vertically discretized in n layers (or nodes) with uniform temperature, considering a constant and incompressible volume of water in the tank (Figure 6).

Figure 8 .
Figure 8.Heat pump with recirculation loop model

Figure 9 .
Figure 9. Solar field validation Figure 9 shows the solar global tilted irradiation and the outlet temperature of a 5000 m² solar field during three

Figure 14 .
Figure 14.Heat pump COP validation All component models are now validated separately.The next step is to simulate a full plant with controllers.

Figure 16 .
Figure 16.Full plant simulation -Solar field temperature

Figure
Figure 21.Optimization methodology

Figure 22 .
Figure 22.Mass flows of each hydraulic loop.The time scale reads as follow (likewise for all following graphs): 0.5 day = 12PM, 1.0 day = 12AM.

Figure 24 .
Figure 24.Heat pump evaporator temperature.Tin evap sim and T from tank overlap across the whole time range.

Figure 25 .
Figure 25.Solar Field TemperatureFinally, Figure26shows the instant profit of the plant, depending on the heat supplied and electricity consumption.We can also see the shift of the heat production time from the first day to the first night.

Figure 26 .
Figure 26.Instant value of the plant

Table 1
gives the improvement provided by the optimizer.The gain seams huge (+21.4% of profit), but it is still not considering what is happening after the two days.

Table 1 .
Plant economic results for 2 simulated days