Solar Field Output Temperature Optimization Using a MILP Algorithm and a 0D Model in the Case of a Hybrid Concentrated Solar Thermal Power Plant for SHIP Applications

: Using solar power for industrial process heat is an increasing trend to ﬁght against climate change thanks to renewable heat. Process heat demand and solar ﬂux can both present intermittency issues in industrial systems, therefore solar systems with storage introduce a degree of freedom on which optimization, on a mathematical basis, can be performed. As the efﬁciency of solar thermal receivers varies as a function of temperature and solar ﬂux, it seems natural to consider an optimization on the operating temperature of the solar ﬁeld. In this paper, a Mixed Integer Linear Programming (MILP) algorithm is developed to optimize the operating temperature in a system consisting of a concentrated solar thermal ﬁeld with storage, hybridized with a boiler. The MILP algorithm optimizes the control trajectory on a time horizon of 48 h in order to minimize boiler use. Objective function corresponds to the boiler use, for completion of the heat from the solar ﬁeld, whereas the linear constraints are a simpliﬁed representation of the system. The solar ﬁeld mass ﬂow rate is the optimization variable which is directly linked to the outlet temperature of the solar ﬁeld. The control trajectory consists of the solar ﬁeld mass ﬂow rate and outlet temperature, along with the auxiliary mass ﬂow rate going directly to the boiler. The control trajectory is then injected in a 0D model of the plant which performs more detailed calculations. For the purpose of the study, a Linear Fresnel system is investigated, with generic heat demand curves and constant temperature demand. The value of the developed algorithm is compared with two other control approaches: one operating at the nominal solar ﬁeld output temperature, and the other one operating at the actual demand mass ﬂow rate. Finally, a case study and a sensitivity analysis are presented. The MILP’s control shows to be more performant, up to a relative increase of the annual solar fraction of 4% at 350 ◦ C process temperature. Novelty of this work resides in the MILP optimization of temperature levels presenting high non-linearities, applied to a solar thermal system with storage for process heat applications.


Introduction
SHIP (Solar for Industrial Process Heat) is gaining more and more attention in research and in industry worldwide. Recent European projects regarding Solar Thermal are oriented towards SHIP such as InSun, SHIP2Fair, ASTEP or FriendSHIP. Besides, a dedicated IEA task (task 64) is devoted to Solar Process Heat and entities like Solar Heat Europe are promoting the use of solar thermal technology for renewable heating and cooling.
Indeed, heat for industry accounts for 34% of the total energy need, and represents 73% of the total energy needed by industries [1]. Nowadays, less than 5% of this energy need is met by renewables. Solar thermal energy is one of the most efficient ways to provide renewable heat. As a result, there is a growing interest in solar thermal due to its potential to save thousands of tons of greenhouse gas emissions in the context of climate change.
Most of the actual SHIP plants are working at low temperature (below 150 • C) with solar fields based on flat plate collectors or evacuated tube collectors. Medium temperature (150-400 • C) is also investigated: this corresponds to 27% of the industrial heat need [2], and can be addressed by linear concentrators, i.e., linear Fresnel (LFR) and parabolic trough (PTC). Some of these systems were implemented in different industries or commercial and administrative places, as is presented per example in Pietrushka et al. [3]. The website SHIP plants [4], created in the framework of IEA's task 49, provides a non-exhaustive list of solar thermal plants through the world; 336 projects are referenced in 2021, 69 of which use concentrating collectors. [2] Solar system can almost always benefit from thermal storage to decorrelate the production and the use of the heat, allowing a broader solar fraction. It also often needs a complementary source of heat, because of the daily, weekly and seasonal variability of solar radiation. Moreover, the heat demand varies over time, depending on the season, the time of the week and the process schedule. Because of the complexity of those kinds of systems, optimization on a mathematical basis enables higher revenues and therefore should favor SHIP use through the world.
Little investigations were performed on operating optimizations of SHIP systems. Wallerand et al. [5] used a MILP algorithm (Mixed Integer Linear Programming) to optimize the design and control of a solar dish with a photovoltaic cells for Hybrid Concentrated Photovoltaics; a stratified thermal storage was also present. Optimization was performed on the thermal part, considering a constant demand; temperature levels were kept constant, and optimization was performed on a yearly basis: the objective was to compare with a TRNSYS precise formulation. The MILP algorithm proved to be faster in estimating annual revenues. They found a 5% difference, and considered that this MILP algorithm was suited for quick estimations. Omu et al. [6] performed design and control optimization on a solar domestic hot water system, through the use of two iterative MILP: one minimized the design cost, while another evaluated the energy outcome of the system for one year. Iteration between the sizing and the operation module was done as such: sizing module gave the size of the system components, while the operation module returned the efficiency of the collectors for each hour, in order to simplify both MILP algorithms and allow for quicker convergence and optimal design of the system. Tilahun et al. [7] and Fitsum et al. [8] performed an agent-based optimization on design and control of a dyeing textile industry in Ethiopia. In [7], control optimization is performed in order to keep constant the outlet temperature of the storage, in opposition to a constant mass flow rate in the solar collector. Three subprocesses are investigated. The constant outlet temperature strategy shows a decrease of 0.3 years for the payback period. Ghazouani et al. [9] performed an optimization based on genetic algorithms, in which the hourly mass flow rates in a small PTC and the different design parameters were optimized. No storage unit was considered, and the heat demand was assumed as constant over the year. Results showed that lowering the inlet temperature and increasing the mass flow rates increased the exergetic efficiency.
MILP algorithms are often used to continuously define the best scheduling strategies for a given time horizon. For example, many articles describe its use for electricity production in a Concentrated Solar Power (CSP) plant. Because of the varying cost of electricity in day-ahead, spot market and ancillary services, an optimization regarding storage and electricity price can be performed and lead to higher revenues, such as in Wagner et al. [10], Vasallo and Bravo [11] optimized the plant schedule through a MILP procedure and then run a simulation of the real power plant, with the MILP results, to estimate the real production profits-they obtain a relative increase in profits from 3% to 4%. Uncertainty upon future solar radiation is important for such works: different methods, such as robust formulation in He et al. [12] or chance-constrained programming were evaluated. Petrollese et al. [13] made a MILP model to investigate the differences between robust, stochastic and deterministic optimization. They found out that taking forecast uncertainty into account increases the revenue, but that choice between both methods ought to be done depending on the quality of the forecast. Dowling et al. [14] made a different approach: they used a MILP algorithm for choosing the on-off period of the plant, and then solved a NLP for accurate results; this was justified by the necessity to take fast dynamics into account when working on real time market. Discussion over the use of MINLP (Mixed Integer Non-Linear Programming) and meta-heuristic optimization such as particle swarm optimization or artificial neural network can be found in Omu et al. [6]. The main reason for the choice of a MILP algorithm is the quality of the calculated optimum and the low computation cost.
Other authors used MILP to investigate CSP dispatch with hybridization; a few differences arise, regarding the objective function implemented: solar fraction is put in evidence, in order to reduce fuel consumption; a multi-objective function can be defined. Pousinho et al. [15] applied a MILP to a fossil-fuel/CSP hybrid, regarding the same markets than above; they implemented robust optimization to take into account the forecast uncertainty. Pousinho et al. [16] investigated the day-ahead market and the reserve market for a wind/CSP hybrid, and compared the results with a genetic algorithm; they found out the MILP to be more efficient. Wind/CSP hybridization, along with electrical heaters, was investigated by Yang et al. [17] and Zhao et al. [18] through the use of MILP algorithm; the latter introduced chance-constrained algorithm for taking under consideration the wind and solar uncertainty. Hamilton et al. [19] investigated a CSP-TES-PV-batteries system, with dispatch optimization with a MILP algorithm.
The concept of flexible heat integration was investigated by Rashid et al. [20] and Ellingwood et al. [21] for CSP (electricity production). As heat losses get larger with higher temperature and/or lower Direct Normal Irradiance (DNI), operating at lower temperature leads to higher revenues because of a better efficiency and so a lower use of the complementary heat source (called boiler in the following). Rashid et al. [20] showed an increase of the solar fraction of the hybrid plant compared to a hybrid plant with optimum flexible heat integration from 13.3% to 15%.
The novelty of this paper resides in the use of a MILP algorithm for optimization of the boiler use by varying the temperature level and mass flow rates instead of energy as was done in previous works: to the authors' knowledge, such optimization on temperature levels with a MILP algorithm was only recently performed by Wirtz et al. [22] for district heating. This approach is applied to optimizing the solar field output temperature for SHIP applications considering a thermal storage along with a boiler.

Concentrated Solar Thermal System under Consideration
The heat integration point and the nature of storage are determining factors in the design of a SHIP system [2]. For the purpose of the study, it is therefore necessary to choose one system architecture on which to perform the optimization. Nonetheless, the proposed MILP algorithm can be applied to any kind of solar collector. Indeed, the temperature range and the collector-receiver performances underline the interest of optimization: a higher process temperature will increase the heat losses and give more range for this optimization. This section describes the envisaged system and its simulation, along with the main hypothesis.

Schematic and System Description
The considered concentrating solar system is constituted of a solar field (SF), with a two-tank thermal storage in series and a boiler. Figure 1 shows a solar field connected in series with a hot tank. Mass flow rate from the solar field ( hot tank (M st ) and its temperature (T st ). The cold tank is not modeled. It is considered that the mass fluctuations in the hot tank are compensated by the cold tank. hot tank ( ) and its temperature ( ). The cold tank is not modeled. It is considered that the mass fluctuations in the hot tank are compensated by the cold tank.
There is a third inlet to the boiler that permits to bypass the solar field and storage (̇).The three mass flow rates at the inlet of the boiler must be equal to the process' demand mass flow rate (̇). If necessary, the boiler heats up the total HTF's mass flow rate in order to reach the process' demand temperature.
The system's cold inlet mass flow rate of HTF equals the system's hot outlet mass flow rate of HTF, which ensures that the total HTF's mass remains constant in the considered system. The block separation between SF Model and Plant Model comes from the optimization performed on the solar field, as explained in Section 2.2.

Main Hypothesis
Different hypotheses were made in this first order 0D model, as the objective of this paper is to assess the interest of temperature optimization in a solar field for SHIP: 1. Some thermal and parasitic losses are neglected: heat losses from night's recirculation to avoid freezing, line pumping power... 2. Perfectly mixed hot tank: uniform temperature in the tank. 3. Heat losses' approximation in the SF: they are classically taken as the module's losses, with the average temperature of the SF. 4. Constant fluid properties with temperature. There is a third inlet to the boiler that permits to bypass the solar field and storage ( . m aux ).The three mass flow rates at the inlet of the boiler must be equal to the process' demand mass flow rate ( . m demand ). If necessary, the boiler heats up the total HTF's mass flow rate in order to reach the process' demand temperature.
The system's cold inlet mass flow rate of HTF equals the system's hot outlet mass flow rate of HTF, which ensures that the total HTF's mass remains constant in the considered system. The block separation between SF Model and Plant Model comes from the optimization performed on the solar field, as explained in Section 2.2.

Main Hypothesis
Different hypotheses were made in this first order 0D model, as the objective of this paper is to assess the interest of temperature optimization in a solar field for SHIP:

1.
Some thermal and parasitic losses are neglected: heat losses from night's recirculation to avoid freezing, line pumping power...

2.
Perfectly mixed hot tank: uniform temperature in the tank.

3.
Heat losses' approximation in the SF: they are classically taken as the module's losses, with the average temperature of the SF.

4.
Constant fluid properties with temperature.

5.
Constant boiler efficiency. No limitations on the power (maximal or minimal) are considered.

6.
No heat losses in the cold tank: no calculations were made on the cold tank, although heat losses do happen in it as well. 7.
Perfect solar forecast: this hypothesis is chosen in order to exclude the influence of the forecasting model.  Figure 2 describes the co-simulation approach performed in this work.

Algorithm Structure
5. Constant boiler efficiency. No limitations on the power (maximal or minimal) are considered. 6. No heat losses in the cold tank: no calculations were made on the cold tank, although heat losses do happen in it as well. 7. Perfect solar forecast: this hypothesis is chosen in order to exclude the influence of the forecasting model. Figure 2 describes the co-simulation approach performed in this work. Inlet data are the process demand curve (hourly mass flow rates and return temperature) and the process temperature, the meteorological data (hourly solar irradiation) and the design of the system, i.e., the solar field (number of modules per loop, number of loops, aperture area, nominal and maximal mass flow rates, optical efficiency and heat losses); as well as the design of the storage (maximal mass in the storage, heat losses). The temperature range is discretized; the MILP algorithm choses the best temperature in a given set. and are taken as the minimum return temperature and the process temperature respectively. The SF Model is described in Section 3.1. Solar production, i.e., heats and mass flow rates, are calculated for every hour on a time horizon of 48 h and for every temperature issued from the discretization, considering thermal inertia of the SF. Those calculations are performed upstream of the MILP algorithm, enabling linear optimization despite the high non-linearities of the SF's equations.

Algorithm Structure
The MILP algorithm is described in Section 4. This optimization algorithm chooses the optimal trajectory that minimizes the use of the boiler, through the optimization of a simplified linearized system: it is not possible to optimize on the real system as it contains numerous non-linearities. It is run on 48 h in order for the MILP algorithm to optimize the Inlet data are the process demand curve (hourly mass flow rates and return temperature) and the process temperature, the meteorological data (hourly solar irradiation) and the design of the system, i.e., the solar field (number of modules per loop, number of loops, aperture area, nominal and maximal mass flow rates, optical efficiency and heat losses); as well as the design of the storage (maximal mass in the storage, heat losses). The temperature range is discretized; the MILP algorithm choses the best temperature in a given set. T min and T max are taken as the minimum return temperature and the process temperature respectively. The SF Model is described in Section 3.1. Solar production, i.e., heats and mass flow rates, are calculated for every hour on a time horizon of 48 h and for every temperature issued from the discretization, considering thermal inertia of the SF. Those calculations are performed upstream of the MILP algorithm, enabling linear optimization despite the high non-linearities of the SF's equations.
The MILP algorithm is described in Section 4. This optimization algorithm chooses the optimal trajectory that minimizes the use of the boiler, through the optimization of a simplified linearized system: it is not possible to optimize on the real system as it contains numerous non-linearities. It is run on 48 h in order for the MILP algorithm to optimize the storage for the night and the following day; if run only for 24 h, the optimization would just discharge the storage.
The Plant Model is described in Section 3.2. This model deduces the overall system behaviour from the optimized control trajectory on 24 h. Indeed, the control trajectory consists of (  (24)), is then sent back as input for the algorithm in order to loop over the year, while interesting data are stored for future analysis.

Physical Model
In this section, the physical model is presented; it is physical in opposition to the MILP algorithm, which is a mathematical approximation of the physical model: MILP algorithms necessitate linearized equations. In Section 3.1, the SF model is presented, separately from the plant model, as the SF Model results are injected in the MILP algorithm. In Section 3.2, the Plant Model is described, as the MILP algorithm equations (Section 4) are derived from the Plant Model.

0D Solar Field Model
The chosen solar field is made of LF-11 Fresnel collectors from Industrial Solar [23] which includes a vacuum tube absorber. A North-South orientation is chosen.

• Sun position
The sun position is computed with the declination and hour angle related to the day and hour of the simulation, using a simplified version of the PSA algorithm for declination. The solar vector is computed as in Blanco et al. [24], considering (Ox) oriented towards east, (Oy) towards north and (Oz) towards the vertical.

•
Optical efficiency The optical efficiency is computed using an interpolation on the longitudinal (θ L ) and transversal angle(θ T ) of their respective incidence angle modifiers (IAM), shown in Figure 3. Those angles are defined as: θ L = arcsin s y and θ T = arcsin − s x cos(θ L ) , with s x the solar vector coordinates pointing to the east and s y to the north. This definition is adapted from Morin et al. [25] and can be easily retrieved through algebraic manipulations. The optical efficiency is then computed as: With ̇, the heat losses in W/m 2 of aperture area, , the outlet temperature of the solar field in °C, , the inlet temperature of the solar field in °C and the ambient temperature in °C, and 1 and 4 the heat losses correlation coefficients. For the LF-11 Fresnel collector, 1 = 0.032913 2 . and 4 = 1.48 * 10 −9 2 . 4 . Those coefficients are taken from a single module and extrapolated to the entire SF, assuming that the temperature evolves linearly in the SF. With η 0 the optical efficiency at normal incidence. For the LF-11 of Industrial Solar, η 0 = 0.686 [23].

• Heat losses
The SF's heat losses are calculated with polynomial heat losses: . q losses,SF (T SF,out , T SF,in , T amb ) = a 1 * T SF,out + T SF,in With . q losses,SF the heat losses in W/m 2 of aperture area, T SF,out the outlet temperature of the solar field in • C, T SF,in the inlet temperature of the solar field in • C and T amb the ambient temperature in • C, and a 1 and a 4 the heat losses correlation coefficients. For the LF-11 Fresnel collector, a 1 = 0.032913 W m 2 .K and a 4 = 1.48 * 10 −9 W m 2 .K 4 . Those coefficients are taken from a single module and extrapolated to the entire SF, assuming that the temperature evolves linearly in the SF.
The heat losses of the LF-11 Fresnel collector as a function of the outlet temperature are shown in Figure 4.

•
Heat from the solar field

Heat losses
The SF's heat losses are calculated with polynomial heat losses: ̇, ( , , , , ) With ̇, the heat losses in W/m 2 of aperture area, , the outlet temperature of the solar field in °C, , the inlet temperature of the solar field in °C and the ambient temperature in °C, and 1 and 4 the heat losses correlation coefficients. For the LF-11 Fresnel collector, 1 = 0.032913 2 . and 4 = 1.48 * 10 −9 2 . 4 . Those coefficients are taken from a single module and extrapolated to the entire SF, assuming that the temperature evolves linearly in the SF.
The heat losses of the LF-11 Fresnel collector as a function of the outlet temperature are shown in Figure 4.


Heat from the solar field The solar production is then expressed as: With . q abs the heat absorbed by the solar field in W, T SF,in and T SF,out respectively the inlet and outlet temperatures of the solar field, → s sun the solar vector, T amb the ambient temperature, DNI the direct normal irradiation in W/m 2 , η opt the optical efficiency, and . q losses the heat losses in W/m 2 .

•
Mass flow rate of the SF from outlet temperature, considering inertia.
Wagner et al. [26] describes the heat balance at a node as done in SAM (see Figure 5). In the present simple model, the node is the entire SF. U is the internal energy of the solar field, T is the average temperature of the solar field, m SF * C p * T out is the outlet heat of the SF.
With ̇ the heat absorbed by the solar field in W, , and , respectively the inlet and outlet temperatures of the solar field, ⃗ the solar vector, the ambient temperature, DNI the direct normal irradiation in W/m 2 , the optical efficiency, and ̇ the heat losses in W/m 2 .
 Mass flow rate of the SF from outlet temperature, considering inertia.
Wagner et al. [26] describes the heat balance at a node as done in SAM (see Figure 5). In the present simple model, the node is the entire SF. U is the internal energy of the solar field, ̅ is the average temperature of the solar field, ̇=̇ is the SF's mass flow rate, ̇=̇ * * is the inlet heat to the SF, while ̇=̇ * * is the outlet heat of the SF. The heat balance gives: As the outlet temperature is given by the control model, mass flow rate has to be deduced from this equation: The heat balance gives: .
Energies 2021, 14, 3731 As the outlet temperature is given by the control model, mass flow rate has to be deduced from this equation: In general, mass flow rate is limited by various technological or physical constraints. In this work, the maximal mass flow rate written . m SF,max , is considered as twice the nominal mass flow rate.
Furthermore, according to [26], the internal energy variation can be expressed as: with (mc) bal the inertia of the receiver in J/K, m HTF the HTF mass present in the receiver in kg, C p,HTF the specific heat capacity at constant pressure of the HTF in J/kg·K and the average temperature in the solar field. The inertia term is written as: I t = (mc) bal + m HTF * C p,HTF , in J/K.

• Validation
Validation of the solar field model was performed using SAM [27]: the Linear Fresnel Molten Salt Model was run with constant fluid properties with temperature. Inlet and outlet temperatures of SAM were reinjected in the 0D solar field model, which calculated the expected mass flow rate. The relative difference, defined in this part as m SAM the daily mass flow rate in SAM and . m 0D the daily mass flow rate from the proposed model, gave a value of 9.8%. A difference of 2 • could be observed on the longitudinal and transversal angles, which can come from a difference in the reference time for solar position evaluation-SAM defines the solar vector differently. Using the solar position of SAM gave a 3.6% relative difference. Remaining differences are expected to come from inertia considerations and modeling assumptions (mainly 1D vs. 0D).
In this section, the SF Model was presented. In the following section, the Plant Model is exposed: for this model, the control trajectory, i.e., ( . m SF , T SF , . m aux ) is known, for every hour h of the upcoming day, as well as the initialization state, which in our case can be summarized as the hot HTF mass present in the hot tank and its temperature at t 0 (midnight).

Mass and Mass Flow Rates Balance (Storage and Demand)
The demand mass flow rate has to be met, which gives: The dispatch of the mass flow rate from the solar field gives: The storage's mass balance is then written as: Reinjection of Equations (7) and (8) in Equation (9)  Discretization gives: Verification is then performed on the control trajectory: if the stored mass of hot HTF tends to exceed the maximum or minimum storage capacity, either a defocus is imposed or an increase of the auxiliary mass flow rate. m SF,h .

Temperature of the Storage
The storage is supposed perfectly mixed, which gives a unique storage temperature T st of the hot tank.
With h st the specific enthalpy of the fluid in the hot tank, h SF,out the specific enthalpy of the fluid at the outlet of the Solar Field and . q losses,st the heat losses of the storage. Heat losses of the storage are modeled through a constant heat transfer coefficient: With U Newton's heat transfer coefficient, A the external area of the tank, T st the temperature of the storage and T amb the ambient temperature.
The reference enthalpy is 0 at 0 K. After discretization of Equation (12), it can be written that: One can therefore obtain the temperature of the storage from Equation (14); if the storage is empty, equilibrium with ambient temperature is reached.

Boiler Model
The boiler heats up the HTF coming from the storage, the HTF coming from the SF and the auxiliary mass flow rate. Therefore the heat production of the boiler can be written as: Energies 2021, 14, 3731 10 of 22

Plant Metrics
In this section, the definition of the solar fraction and of the defocusing percentage in the solar field is presented. Solar Fraction is defined as the part of the heat not supplied by the boiler: Defocus is defined as the difference between the potential power at the outlet and the actual power. It is an approximation as defocus is not continuous but rather constant per pieces. The defocused energy is written as the difference of production between the possible mass flow rate (from Equation (5)) and the chosen mass flow rate, multiplied by C p ∆T: .
After presenting the SF and plant model, the following section presents the MILP model derived from those models. Two simpler and deterministic control approaches are also introduced for comparison.

General Description of MILP Algorithm and Solver
MILP algorithms are optimization algorithms consisting in a linear cost function (such as boiler use, operating costs, yearly revenues...) to be minimized or maximized, and linear constraints representing the model. Optimization variables can be integers or continuous. The general mathematical formulation is [28]: with f , b min , b max , x min , x max vectors of n coordinates, matrix and J a subset of {0, . . . , n}.
In this work, the solver used is the open-source MIP solver Coin-Or CBC v.2.10.5, downloaded from [29]; time step is one hour. Typical solve time is around 1.5 min for the entire year.

Methodology: How Is the Mathematical Problem Posed
As is described in Figure 2, considering the temperature discretization from the process data, the solar production is calculated for every temperature at every hour h, considering every temperature at h-1: because thermal inertia is taken into account, every possibility has to be calculated. Binary variables per temperature level at hour h and hour h-1 choose the optimal temperature level and mass flow rate, the one that minimizes the boiler use. As we do not go into details of the dispatch, the boiler use is taken as the addition of the heat up of the outlet SF's mass flow rate and the compensation of heat losses in the storage.
The optimization is performed considering a time horizon of 48 h, with a receding horizon of 24 h. This means that the optimization is performed on 48 h, and that the plant model is used to get accurate results on 24 h, and then optimization is restarted: this way, optimization takes the storage into account for the following day. Indeed, if optimization is run on 24 h, the storage just empties itself.

From Algorithm Parameters to MILP Algorithm's Parameters
In this chapter, the words "parameters" and "variables" will be used from the MILP point of view, and not from the overall algorithm point of view. For example, from the overall algorithm point of view, the different calculated mass flow rates are variables; but from the MILP algorithm's point of view, these are parameters, calculated and fixed beforehand. The algorithm parameters such as the solar field area, process temperatures and mass flow rate are noted straight. The MILP algorithm's parameters such as solar production and the temperature range are noted in italic. On the other side, the operating temperature and the auxiliary mass flow rate are optimization variable; they are noted in bold. t is the time set; hourly resolution is chosen, with a time horizon of 48 h.

•
Temperature range discretization: With N discr the number of discretization chosen by the user.
• Absorbed heat from the sun, from Equation (3): • Mass flow rate calculation: Discretization of the inertia term from Equation (6) is necessary: , with T h the temperature at hour h and ∆t = 3600 s the time discretization interval. Therefore, the mass flow rate at hour h, supposing T h = T i and T h−1 = T j can be calculated from Equation (5):

From MILP Parameters to MILP Formulation
In this section, variables are related to each other and to the MILP algorithm's parameters. This is done through the constraints. The cost function is the objective to be minimized by the algorithm.
In our approach, binary variables were implemented to choose the most interesting temperature; as inertia is taken into account, those binary variables are to be valued over a minimum of 2 time steps: The implementation of this definition as constraints is given by: z h,i,j integer, h ∈ T, (i, j) ∈ {0, . . . , N discr } One can observe that T h = ∑ Boundary constraints for those binary variables are also necessary: when no production is possible, it might be of better interest to preheat the solar field in order to produce more, and at higher temperature, as soon as the demand restarts. Nonetheless, it is necessary to force the binary variable to be 0 when its corresponding mass flow rate is also 0. If not, the algorithm would choose to work at the highest possible temperature in order to withdraw this energy at restart; which is not realistic. This gives: Following this definition of the binary variables, it is possible to implement the mass flow constraint: The optimization will naturally maximize . m SF,h ; specifying an inequality instead of an equality allows partial defocus, in case the storage tank is filled, or if it is of higher interest to store later during the day, for example.
The mass flow rate is also constrained by the maximal mass flow rate: The conservation of mass in the storage is discretized as in Equation (11), which gives: The storage is constrained by its capacity, which can be written as: (28) Heat losses in the storage have to be considered: in this work, they are integrated to the MILP in the cost function-as the boiler will have to compensate those losses. A 0D estimation of those losses can be written through the following equations: as V storage = M storage ρ storage = πR 2 h, the surface in contact with the hot fluid can be written as A = 2πRh + πR 2 , with V storage the hot fluid volume, R the radius of the tank and h the height of fluid in the tank. The constant losses through the bottom (πR 2 ) would not influence the cost function of the MILP algorithm, they are therefore ignored; heat losses in the MILP algorithm are written as: Finally, the boiler consumption is written as: .
It has to be noted that the boiler consumption might be slightly overestimated in case of defocus, as . m SF,h is inferior to the chosen . m SF,h,i,j as Equation (26) is an inequation. Nonetheless, defocus mostly happens when solar energy is too high, and therefore the working temperature is T process , which puts the term equal to zero. The difference with the Plant Model is quite visible here: the boiler use is estimated at solar production time and not at boiler use time. The cost function to be minimized is expressed as:

Control Models for Comparison
To evaluate the approach of MILP temperature optimization in the solar field, two other control approaches-i.e., for the choice of the working temperature-were implemented.

Control Approach 1 (CA1)
The CA1 approach tries to work at T process whenever it is possible; if it is not (i.e., too high heat losses), it will work at the highest possible temperature.
If the solar field produces more mass flow rate than necessary for the demand, the excess part is stored. If the storage is full, defocus is applied. If the solar field produces less mass flow rate than necessary for demand, storage is discharged; if there is not enough in storage, the auxiliary mass flow rate is used through the boiler.

Control Approach 2 (CA2)
The second control approach is more complex and the corresponding logic flowchart is proposed in Figure 6. If the solar field can produce more than the demand mass flow rate at the process temperature, the solar field works at the process' temperature, storing the excess. If storage is full, defocus is applied. If the solar field cannot produce the necessary heat at the demand temperature, the demand mass flow rate is sent into the solar field, heating it up as much as possible and the boiler only completes (no storage). When solar energy is not sufficient, then if the storage is not empty, it is used; else, it is completed by the boiler.

Case Study
In this subsection, a case study is detailed, i.e., the parameters for the overall algorithm. Table 1 resumes the different parameters.  If the solar field can produce more than the demand mass flow rate at the process temperature, the solar field works at the process' temperature, storing the excess. If storage is full, defocus is applied. If the solar field cannot produce the necessary heat at the demand temperature, the demand mass flow rate is sent into the solar field, heating it up as much as possible and the boiler only completes (no storage). When solar energy is not sufficient, then if the storage is not empty, it is used; else, it is completed by the boiler.

Case Study
In this subsection, a case study is detailed, i.e., the parameters for the overall algorithm. Table 1 resumes the different parameters.  [30]. Process data: a demand curve is specified; the same curve was taken every day, with a demand mass flow rate of 50 kg/s from 8 a.m. to 19 p.m., solar hour and 10 kg/s in the remaining time-lapse.

Hourly Results
In this section, results from co-simulations for the three control models are exposed: examples of summer and winter days are discussed. Figure 7 reports the main simulation's hourly results on two consecutive summer days, in June. The following variables are displayed: DNI (Figure 7a), the absorbed power on the solar field (Figure 7b), the ambient temperature (Figure 7c), the demand mass flow rate (Figure 7d), and then for the 3 dispatch strategies, the solar field temperature (Figure 7(e.x)), the solar field mass flow rate (Figure 7(f.x)), the storage level in percentage (Figure 7(g.x)), the storage temperature (Figure 7(h.x)), the power of the boiler (Figure 7(i.x)), the auxiliary mass flow rate (Figure 7(j.x)), the power losses in the SF (Figure 7(k.x)) and the energy defocused in the solar field (Figure 7(l.x)).
On the first day, all strategies present defocus (Figure 7(l.x)) and work at the process temperature of 350 • C (Figure 7(e.x)). The main difference between the three strategies on this day is when the defocus takes place: MILP's control defocuses in the morning (Figure 7(l.1)) and fills the storage during the afternoon (Figure 7(g.1)), while the two other strategies fill the storage in the morning (Figure 7(g.2,g.3)) and then defocus in the afternoon (Figure 7(l.2,l.3)). This dissymmetry is as well observed on the mass flow rates (Figure 7(f.x)). It is explained by the heat losses in the storage: MILP strategy intends to lower heat losses in the storage, by storing later, and therefore maintaining a higher temperature in the storage (Figure 7(h.1), as compared to Figure 7(h.2,h.3)). It can be observed that the MILP doesn't fill the storage at 100%, as 80% is enough to go through the night: no auxiliary mass flow rate is needed (Figure 7(j.1)). Figure 7(i.3) shows CA2 limits: there is a peak in boiler use at 19 p.m., because rather than using the storage, CA2 chooses to run the solar field at the demand mass flow rate although the outlet temperature is not the process temperature (Figure 7(e.3)). Figure 8 reports the main results obtained by using the dispatch strategies MILP, CA1 and CA2 in cosimulations for two consecutive winter days in February, in the same way as in Figure 7. The main differences are observed on the first reported day characterized by a variable DNI, which will be commented hereafter. Figure 8(h.1) shows that MILP's control uses the storage at lower temperature than process temperature. Indeed, the maximum storage temperature is 150 • C while the process temperature is set to 350 • C. CA1 and CA2 do not use the storage as no excess production is possible for those control approaches upon that day (Figure 8(g.2,g.3)). This allows lower heat losses in the solar field for the MILP strategy (maximum value of 0.4 MWh in Figure 8(k.1)) on the first day while CA1 show values twice as high (Figure 8(k.2)). CA2 shows low heat losses as well, compared to CA1: it works at higher mass flow rate and so lower solar field temperature. Indeed, CA1 operates at 350 • C (Figure 8(e.2)) while the maximum operating temperature of CA2 is 315 • C (Figure 8(e.3)) and 180 • C for MILP's control (Figure 8(e.1)). Therefore the solar fraction is improved thanks to the MILP's control; over this two-day range, the solar fraction reaches 35.2% for MILP, 33.2% for CA1 and 34.9% for CA2. MILP's control boiler usage (Figure 8(i.1)) is smoothed during the DNI peak because of the storing strategy that compensates the fluctuations, as compared to CA1 and CA2 which stopped the boiler during the DNI peak around 1 p.m. Figure 8(e.2) illustrates the CA1 s strategy of keeping a constant temperature, whereas Figure 8(f.3) illustrates the CA2 s strategy of keeping the solar field flow rate at the demand mass flow rate.  Figure 8 reports the main results obtained by using the dispatch strategies MILP, CA1 and CA2 in cosimulations for two consecutive winter days in February, in the same way as in Figure 7. The main differences are observed on the first reported day characterized by a variable DNI, which will be commented hereafter. Figure 8(h.1) shows that MILP's control uses the storage at lower temperature than process temperature. Indeed, the maximum storage temperature is 150 °C while the process temperature is set to 350 °C. CA1 and CA2 do not use the storage as no excess production is possible for those control approaches upon that day (Figure 8(g.2,g.3)). This allows lower heat losses in the solar field for the MILP strategy (maximum value of 0.4 MWh in Figure 8(k.1)) on the first day while CA1 show values twice as high (Figure 8(k.2)). CA2 shows low heat losses as well, com- boiler usage (Figure 8(i.1)) is smoothed during the DNI peak because of the storing strategy that compensates the fluctuations, as compared to CA1 and CA2 which stopped the boiler during the DNI peak around 1 p.m. Figure 8(e.2) illustrates the CA1′s strategy of keeping a constant temperature, whereas Figure 8(f.3) illustrates the CA2′s strategy of keeping the solar field flow rate at the demand mass flow rate.  Figure 9 presents the monthly results for the three control strategies. Figure 9a displays the monthly solar fraction, Figure 9b displays the average percentage of mass in the hot tank per month, Figure 9c presents the energy defocused in the solar field and Figure  9d displays a box repartition of the mass flow rate in the solar field when it is greater than 0. Figure 9a shows that MILP's control always reaches a higher solar fraction than the two CAX control; CA2 surpasses CA1 in winter, due to its lower temperature in the solar field. Nonetheless, CA2′s discharge strategy in summer implies greater defocusing (Figure 9c), as its storage is more filled on average (Figure 9b). Figure 9b highlights that MILP's control uses storage in winter, as compared to the two CAX control approaches; nonetheless, it is  Figure 9 presents the monthly results for the three control strategies. Figure 9a displays the monthly solar fraction, Figure 9b displays the average percentage of mass in the hot tank per month, Figure 9c presents the energy defocused in the solar field and Figure 9d displays a box repartition of the mass flow rate in the solar field when it is greater than 0. Figure 9a shows that MILP's control always reaches a higher solar fraction than the two CAX control; CA2 surpasses CA1 in winter, due to its lower temperature in the solar field. Nonetheless, CA2 s discharge strategy in summer implies greater defocusing (Figure 9c), as its storage is more filled on average (Figure 9b). Figure 9b highlights that MILP's control uses storage in winter, as compared to the two CAX control approaches; nonetheless, it is less used in summer for the reason discussed in Section 5.2. Figure 9d shows that the median mass flow rate in the solar field in the MILP is higher than the two other control approach; CA1 has the lowest median mass flow rate due to the high temperature aimed, and CA2 runs almost constantly at 50 kg/s. MILP shows higher maximal mass flow rate, as it doesn't necessarily work at the highest temperature whenever hourly demand is fulfilled. less used in summer for the reason discussed in Section 5.2. Figure 9d shows that the median mass flow rate in the solar field in the MILP is higher than the two other control approach; CA1 has the lowest median mass flow rate due to the high temperature aimed, and CA2 runs almost constantly at 50 kg/s. MILP shows higher maximal mass flow rate, as it doesn't necessarily work at the highest temperature whenever hourly demand is fulfilled. Figure 9. Monthly results.

Yearly Results
The annual solar fraction gives 50.9% for the MILP, 49.6% for CA1 and 48.6% for CA2. This highlights the better performance of the MILP approach for this base case. In the following, a sensitivity study is proposed to enlarge the analysis.

Sensitivity Analysis
In this section, discussion is driven through a sensitivity analysis on different parameters and over 4 different demand curves.

Methodology
The following demand curves are considered: Figure 10 describes the different demand mass flow rate investigated and their acronym. DC1 corresponds to a constant demand; DC2 corresponds to the case study demand, which is a higher demand during the day and lower one at night; DC3 corresponds to a demand occurring only at night; DC4 corresponds to a batch demand. The same process temperature is taken for each demand curve. All mass flow demand curves were chosen to give the same daily heat demand, in order to obtain the same system design from the following methodology. The annual solar fraction gives 50.9% for the MILP, 49.6% for CA1 and 48.6% for CA2. This highlights the better performance of the MILP approach for this base case. In the following, a sensitivity study is proposed to enlarge the analysis.

Sensitivity Analysis
In this section, discussion is driven through a sensitivity analysis on different parameters and over 4 different demand curves.

Methodology
The following demand curves are considered: Figure 10 describes the different demand mass flow rate investigated and their acronym. DC1 corresponds to a constant demand; DC2 corresponds to the case study demand, which is a higher demand during the day and lower one at night; DC3 corresponds to a demand occurring only at night; DC4 corresponds to a batch demand. The same process temperature is taken for each demand curve. All mass flow demand curves were chosen to give the same daily heat demand, in order to obtain the same system design from the following methodology.
The system is systematically designed as follows: the number of modules in one loop is determined with the nominal mass flow rate of the module, the inlet and outlet temperature at design, under a condition of 900 W/m 2 : m nominal * C p *

Process Temperature
In this section, the influence of process temperature on the Annual Solar Fraction is investigated. Inlet temperature is kept constant at 50 • C. Figure 11 presents the evolution of the annual solar fraction to the variation of the process temperature on the different demand curves for the three control approaches. On every figure, although not instinctive, the solar fraction starts increasing: this is due to the limitation in maximum mass flow rate, not adapted to such low temperatures. DC3 shows low values of solar fraction due to the size of storage and because of a night only demand. In DC3 and DC4, MILP's control does not show a significant decrease of solar fractions with temperature; this is because a better response is given from the MILP's control with higher temperature: storage is much more used, but at lower temperature, and those two demand curves have higher share of demand during the night. Indeed, every curve shows an increase in MILP's control's average mass in storage with the outlet temperature (not shown here), which can be explained by the higher number of loops, due to design methodology. Figure 12 presents the evolution of the annual solar fraction with the increase of heat losses in the solar field. It can be observed that the annual solar fraction is linearly related to this coefficient; nonetheless, MILP's control shows a lower slope than the two other control approaches. MILP algorithm's interest is therefore directly related to the heat losses in the solar field. every figure, although not instinctive, the solar fraction starts increasing: this is due to the limitation in maximum mass flow rate, not adapted to such low temperatures. DC3 shows low values of solar fraction due to the size of storage and because of a night only demand. In DC3 and DC4, MILP's control does not show a significant decrease of solar fractions with temperature; this is because a better response is given from the MILP's control with higher temperature: storage is much more used, but at lower temperature, and those two demand curves have higher share of demand during the night. Indeed, every curve shows an increase in MILP's control's average mass in storage with the outlet temperature (not shown here), which can be explained by the higher number of loops, due to design methodology. Figure 11. Evolution of Solar Fraction as a function of process temperature for the 4 demand curves and 3 control strategies. Figure 12 presents the evolution of the annual solar fraction with the increase of heat losses in the solar field. It can be observed that the annual solar fraction is linearly related to this coefficient; nonetheless, MILP's control shows a lower slope than the two other control approaches. MILP algorithm's interest is therefore directly related to the heat losses in the solar field.  (2)) for the 4 demand curves and 3 control strategies.

Other Sensibility Analyses
Other analyses were performed, and will be only commented here for the sake of  (2)) for the 4 demand curves and 3 control strategies.

Other Sensibility Analyses
Other analyses were performed, and will be only commented here for the sake of brevity: the number of discretization of the temperature range showed very small impact beyond N discr = 8. Solar Multiple would not influence strongly the difference between the control approaches, as well as the size of the storage-which mostly influenced the results of night demand DC3.

Discussion and Perspectives
Rashid et al. [20] performed flexible heat integration for electricity production. Applying optimization on the temperature level raised the solar fraction of the considered solar system from 0.133 to 0.150, which appears coherent with the present simulations.
MILP algorithms present themselves as interesting optimization for real-time optimization: calculation time is very low, which would allow re-calculation for every solar irradiation prediction update. Nonetheless, the complexity of the considered system will be limited, due to the need for linear equations. Future works should investigate the possibility of re-calculation during the day considering various probabilities in the prediction.
The solar field 0D model did not compare perfectly with the SAM's 1D model, due to inertia, the dimension difference and recirculation considerations mostly. Further works will involve 1D detailed solar field model, along with non-constant fluid properties with temperature or even steam properties.
Model of the boiler should incorporate various efficiencies considering load per example, as this varies depending on the technology and/or fuel.
Two-tank storage is not foreseen for SHIP applications, although it might remain interesting for medium and high temperature processes. Further work will involve dualphase thermocline storage due to its lower cost.
Finally, a new dimension of optimization has to be investigated: various processes, or varying temperature demand would allow more flexibility, and therefore open a higher range for optimization.

Conclusions
This paper presented a systematic approach for temperature optimization in the case of concentrated solar plant with storage that assists a boiler to deliver heat to a process. First, a 0D plant model was described, then the MILP algorithm for the optimization was detailed, and finally a case study and a sensitivity analysis were proposed. Results are highly dependent on the demand curve and other inlet parameters. The pertinence of this optimization was shown to be particularly relevant for solar technologies with high temperature range and/or presenting high heat losses. Indeed, at low temperature range and highly efficient system, the relative difference in solar fraction with the other control approaches could reach values as low as 0.05%, whereas it could reach 4% for high temperature demand and/or high heat losses. This approach could be applied to other technologies, and might promote less efficient design at lower costs.
Novelty of this work consists in the use of binary variables for a choice on calculations done upstream of a MILP algorithm. This seems an attractive solution for finding optimal behavior of systems which present high non-linearities in only one part of the system. Informed Consent Statement: Not applicable.

Data Availability Statement:
The data generated for this study as well as Jupyter Notebooks for analysis are openly available in Data for Solar Field Output Temperature Optimization Using a MILP Algorithm and a 0D Model in the Case of a Hybrid Concentrated Solar Thermal Power Plant for SHIP Applications at https://doi.org/10.5281/zenodo.5006230 [31].