Nonlinear model predictive control for thermal balance in solar trough plants

European Journal of Control 67 (2022) 100717 Contents lists available at ScienceDirect European Journal of Control journal homepage: www.elsevier.com/locate/ejcon Nonlinear model predictive control for thermal balance in solar trough plants Antonio J. Gallego a , ∗, Adolfo J. Sánchez b , J.M. Escaño a , Eduardo F. Camacho a , 1 a Departamento de Ingeniería de Sistemas y Automática, Universidad de Sevilla, Camino de los Descubrimientos s/n., 41092 Sevilla, Spain b Department of Mechanical, Biomedical and Manufacturing Engineering, Munster Technological University, Cork T12 P928, Ireland a r t i c l e i n f o Article history: Received 15 September 2021 Revised 22 June 2022 Accepted 2 July 2022 Available online 8 July 2022 Recommended by Prof. T Parisini Keywords: Solar energy Hydraulic model Non-linear optimization algorithm Thermal balance a b s t r a c t The size of existing commercial solar trough plants poses new challenges in applying advanced control strategies to optimize operation. One of these challenges is to obtain a better thermal balance of the loops’ temperature. Since plants are made up of many loops, the efficiency of the loops can vary substantially if a group has been cleaned or affected by dust. This leads to the need to defocus the collectors of the most efficient loops to avoid overheating problems, thus producing energy losses. In order to minimize these energy losses, the input valves have to be manipulated to reduce the temperature difference of the loops. However, when the pipes connecting the loops are very long, the pressure drop and energy losses in those pipes become notorious, affecting the flow distribution. The need to consider the hydraulic model becomes very important. In this paper, a non-linear model predictive algorithm is presented that uses a hydraulic model of the solar field to compute the aperture of the input valves. The results show that when the length of the pipes is increased, the algorithm proposed in this paper obtains better results than other algorithms proposed in the literature. ©2022 The Author(s). Published by Elsevier Ltd on behalf of European Control Association. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) 1. Introduction The great impulse experienced by renewable energy systems is driven by the need to reduce the negative environmental impact that fossil fuels produce. In particular, solar energy is the most abundant renewable energy source. The only limitation of the use of solar thermal energy is its efficient and economical use [4] . Other solar technologies, such as photovoltaics (PV) are a very low-cost alternative, but the difficulty is to store electrical energy, which is more difficult and expensive than to store thermal energy [11] . Since 1980, a large number of commercial solar power plants have been constructed and commissioned around the world. The first commercial solar trough plants were the 30 MW SEGS plants in California (USA), commissioned in the 80s. From 20 0 0 onward, the commissioning of commercial solar thermal plants increased at a very high rate. As examples, the four 50 MW trough plants Solaben I,II, III and VI in Extremadura (Spain) owned by Atlantica ∗Corresponding author. E-mail addresses: agalleg[email protected] (A.J. Gallego), [email protected] (A.J. Sánchez), [email protected] (J.M. Escaño), [email protected] (E.F. Camacho) . 1 +34 954487347 Yield. The 50 MW solar trough plants Andasol I, II and III owned by Cobra/ACS group were constructed in Guadix (South of Spain) [40] . More solar projects were carried out in the USA: the Solana power plant, including thermal storage systems, and the Mojave solar parabolic trough plants, each of them producing 280 MW of electrical power [26,27] . Currently, there are many concentrating solar plants projects around the world as can be seen in [38] . One of the great challenges of the century identified by the National Academy and the European Commission is making solar energy economical and competitive [12,25] . Advanced control strategies and optimization algorithms can play a decisive role in improving the overall efficiency of solar energy plants [2,6,17] . Developing new control algorithms for large-scale solar trough plants is a very challenging issue. The Advanced Grant Optimal Control of Solar Energy Systems (OCONTSOLAR) funded by the European Research Council aims to contribute to these problems [13] . One of the main control objectives in solar trough plants is to regulate the average temperature of the solar field around a set point [1,18,19] . To fulfill this goal, the level of radiation that affects the whole solar field and the optical efficiency of all the loops are considered to be the same. For small solar fields such as the ACUREX plant, this assumption can be considered reasonable [5,9] . However, current commercial solar troughs cover vast expanses of https://doi.org/10.1016/j.ejcon.2022.100717 0947-3580/© 2022 The Author(s). Published by Elsevier Ltd on behalf of European Control Association. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) A.J. Gallego, A.J. Sánchez, J.M. Escaño et al. European Journal of Control 67 (2022) 100717 land. For example, the two 140 MW solar trough plants of Mojave Alpha and Beta are composed of 282 loops each and cover about 700 hectares of land. The SOLANA solar trough plant is even more extensive: It covers approximately 780 hectares of land and consists of 808 loops [27] . In such large plants, the efficiency of loops can vary substantially if a group has been cleaned or affected by dust [33] . This fact leads to the need to defocus the collectors of the most efficient loops to avoid overheating problems, thus producing energy losses, as explained in [34,36] . The valves of the most efficient loops must be opened, if possible, to increase the flow of heat transfer fluid (HTF) to reduce the energy losses. Loop valves are used only in existing commercial plants for steady-state flow balancing. Several preliminary works have been published on this topic. In [33] an optimization algorithm is presented that manipulates the input valves of each loop every 30 minutes. The algorithm tried to compensate for the difference in optical efficiency. In [35] and [8] , a similar optimization algorithm was proposed and tested in a 50 MW large-scale solar trough plant. The algorithm calculated the aperture of the valves and was compared to the case when the input valves were not manipulated. Thermal energy losses due to defocussing actions were significantly reduced, thus improving the efficiency of the plant. In [14] a distributed algorithm is proposed to calculate the aperture of the input valve of the loops to maximize the thermal power of the solar field. The performance obtained is close to the centralized approach but requires less computational time. The previous algorithms consider that the flow distribution is proportional to the aperture of the valve without considering the energy losses in the pipes; that is: the farther loops receive less flow due to energy losses in the pipe. This simplification is correct for small plants. Nevertheless, when the plant is very large and the pipes are very long, as happens with existing commercial plants, the energy losses in the pipes play a significant role in calculating the flow distribution. In [23] a mathematical hydraulic model is developed for a 1 MW solar trough plant in China. The paper showed the importance of considering the hydraulic model in these types of plants by implementing two feedforward control approaches tested on a simulation environment. In this paper, a predictive optimization algorithm for the model is proposed that uses a hydraulic model to compute the energy losses of the plant. A hydraulic model is developed for the ACUREX field. It is shown that when the pipe length increases, the effect of the energy losses becomes perceptible, and the thermal balance worsens when the hydraulic model is not considered. The differences between the algorithm proposed here and other control strategies proposed in [23,33,35] are pointed out below. 1. The algorithm presented here uses a heuristic to improve performance when strong transients affect the plant. This heuristic avoids the computation of the input valves when solar radiation transients affect the field. If the solar radiation is close to 0, the outlet temperature decreases even with the minimum flow. If the optimization algorithm is computed in those conditions, all the valves will be set to the minimum value with undesirable effects. The com putation of input valve apertures is carried out only if the plant is working under stable conditions. If not, the aperture remains unchanged. 2. In this algorithm, a hydraulic model is used to compute the flow distribution depending on the opening valves and how far the loop is from the main pumps. Previous works related to the thermal balance of solar trough plants [33,35] considered that the flow distribution is proportional to the aperture of the input valves. As shown in Section 5 , the algorithm allows working close to the maximum temperature without any defocus actions. The best working temperature depends on the solar radiation reaching the field: on summer days with high levels of radiation, working at higher temperatures may produce higher electrical production [6] . Minimizing the number of defocus actions reduces not only the energy losses but the degradation of the actuators [36] . The hydraulic model used is explained in Section 3 and uses equations similar to those proposed in [23] . 3. In [23] , two control strategies that use a hydraulic model were developed for a 1 MW pilot plant in China. In this case, the optical efficiency was considered the same for the 3 loops, unlike the algorithm proposed here, where the optical efficiency is estimated. Furthermore, the aperture of the valves is computed using a feedforward control strategy, whereas, in this paper, the aperture of the valves is computed employing a nonlinear optimization algorithm. The maximum difference between the maximum and minimum temperature was about 5 ◦C where using the algorithm presented here, the thermal discrepancy between the maximum and minimum outlet temperatures is smaller than 2 ◦C. The hydraulic model used here is solved using a simple iterative algorithm that avoids the need to compute the Jacobian matrix, unlike the Newton-Raphson numerical method proposed in [23] . The hydraulic model used is explained in Section 3 and uses equations similar to those proposed in [23] . 4. The optimization algorithm proposed in this paper uses an unscented Kalman filter (UKF) [37] to estimate the metal-fluid temperature profile instead of a CART algorithm used in [33] . It was found that this approach offered an adequate trade-off between precision and computational time. This approach avoids the need to solve a nonlinear optimization problem for estimating the optical efficiency of the loops since the UKF can estimate both the optical efficiency and the metal-fluid temperature profiles with one algorithm as done in [35] . Similar results are obtained when the pipe length is short, but a better thermal balance is obtained when the pipe length increases. The paper is organized as follows: Section 2 describes the mathematical model of the loop. Section 3 presents the hydraulic model used in this document and the effect of pipe length in flow distribution. Section 4 develops the optimization algorithm. Section 5 presents the simulation results. Finally, the paper comes to a close with concluding remarks. 2. Mathematical modeling of the ACUREX solar field This section describes the equations that govern the dynamics of the solar collector field. The equations were presented and described in [9] . The ACUREX solar field was one of the first experimental solar trough plants commissioned in the 1980s at the Plataforma Solar de Almería (see Fig. 1 ). The field was formed by 10 North-west oriented loops of parabolic trough collectors which provided about 1.2 MW of thermal power for an incident solar radiation of 900 W/m 2 . The length of each loop was 172 m long: 142 m of active parts where solar radiation reaches the metal tube to heat up the heat transfer fluid and 30 m of passive parts composed of joints and pipes connecting the active parts. These passive parts are not reached by solar radiation. There are two ways to model this kind of system: the concentrated parameter model and the distributed parameter model [3] . Since the distributed parameter model is a more precise description of the solar field dynamics, it is used for simulation purposes in this work. The model of a solar collector loop is described by the following system of partial differential equations (PDE) describing the energy balance ( [7] ): 2 A.J. Gallego, A.J. Sánchez, J.M. Escaño et al. European Journal of Control 67 (2022) 100717 Fig. 1. ACUREX collector field (courtesy of PSA). Table 1 List of abbreviations. HTF Heat Transfer Fluid MPC Model Predictive Control NMPC Nonlinear Model Predictive Controller PDE Partial Differential Equation PV Photovoltaics UKF Unscented Kalman Filter Table 2 Parameter Description. Symbol description Units tTime s lSpace m ρDensity kgm −3 CSpecific heat capacity JK −1 kg −1 A Cross sectional area m 2 T (l, t) Temperature K, ◦C q (t) Oil flow rate m 3 s −1 I(t) Direct Solar Radiation Wm −2 n o geometric efficiency Unitless K opt Optical efficiency Unitless G a Collector aperture m T a (t) Ambient temperature K, ◦C H l Global coefficient of thermal loss Wm −2 ◦C −1 H t Coefficient of heat transmission metal-fluid Wm −2 ◦C −1 L wetted perimeter m Re Reynolds number Unitless fBarr friction coefficient Unitless μdynamic viscosity Pa ˙ s r relative rugosity Unitless ggravitational acceleration m / s 2 hp l pressure drop m v velocity m/s L pipe Length of the pipe m ρm C m A m ∂T m ∂t = IK opt n o G −H l G (T m −T a ) −LH t (T m −T f ) ρf C f A f ∂T f ∂t + ρf C f q ∂T f ∂l = LH t (T m −T f ) (1) Where the subindex m refers to metal and frefers to the fluid. The model parameters and their units are shown in Table 2 . The whole plant can be modeled by adding loops in parallel. The PDE system is solved by dividing the metal and fluid into 172 segments of 1 m long each. The integration step is chosen to be 0.5 seconds, and the integration technique is a Euler forward method. The same system of equations is used to model the fluid dynamics of the pipes connecting the loops. The difference is that the Fig. 2. Acurex solar field: flow distribution scheme. incident radiation is nil and the thermal loss coefficient is smaller because they are better isolated than the metal tube of the loops. 3. Hydraulic model of the ACUREX solar collector field In this section, the hydraulic model used in this paper is described. It should be noted that the hydraulic model used in this paper is not intended to be very precise. It uses the well-known Darcy equations [24] to calculate the energy losses in the pipes, and thus the oil flow circulating through each loop can be obtained. Fig. 2 shows a scheme of the ACUREX solar field. The loops are connected to the oil pump through a pipe where the cold fluid is circulating (blue pipe). Then the fluid passes through the loops and is heated. The heated fluid returns to the hot pipe (red line) to the steam generator, where the accumulated heat is consumed, and the cold HTF returns to the pump. In most previous work, it was usually considered that the flow that passes through every loop was the same when all the valves were fully open. However, this simplification is only valid if the pipes connecting the loops are short. The main reason is a drop in pressure in the pipe due to energy losses. These energy losses depend on the length of the pipe. The farther loops receive a lesser amount of flow. Furthermore, there are also energy losses in valves, sudden enlargements, sudden contractions, [10] . 3.1. Computation of the pressure drop in pipes The well-known Darcy equations are used to calculate the pressure drop and energy losses inside the pipes [20] . First, the Reynolds number Re and the Barr friction coefficient fare calculated using Eqs. (2) and (3) as indicated in [28] : Re = ρf q d A f μ(2) 1  f = −2 log 10 (r / 3 . 7 + 5 . 1286 /Re 0 . 89 ) (3) Here, μis the dynamic viscosity and r is the relative rugosity. The Darcy equation for computing the pressure drop is given by: h pl = f L pipe c v 2 2 g (4) 3 A.J. Gallego, A.J. Sánchez, J.M. Escaño et al. European Journal of Control 67 (2022) 100717 Fig. 3. Acurex solar field: closed circuit formed by two loops. Here, L pipe is the length of the loop, and dis the diameter of the pipe. The symbol gstands for gravitational acceleration, and the value used in this article is 9.81 m/s 2 , and v is the velocity of the fluid. The form used in this paper for equation (4) is the following: h pl = 8 fq 2 gπ2 L pipe d 5 (5) It is obtained by substituting v = q A = 4 q πd 2 into equation (4) . The advantage of this form is that it is directly dependent on the flow that flows through the pipe. The energy losses produced in a valve can be modeled as follows. h plv = K v v 2 2 g (6) K v can be found in the manufacturer’s datasheet, and v is the average velocity of the fluid. In this paper, for simplicity, a value of K v of 0.46 is considered. 3.2. Computation of the flow distribution through the loops To obtain the flows that circulate through each loop, the fact that the sum of the pressure drop in a closed circuit has to be equal to zero is used. A closed circuit comprises two consecutive loops and the pipe segments that connect them. For example, closed-circuit 1 consists of loop 1, segment 2 of the hot and cold main pipes, and loop 2 ( Fig. 3 ). Another constraint to be satisfied is that the sum of the flows circulating for all loops must be equal to the flow delivered by the main pump. The conditions mentioned above give rise to a set of nonlinear equations that can be solved using iterative methods to obtain the flows that circulate through every loop (see Fig. 3 ). Let H cold N+1 be the pressure drop in the N+1 segment of the cold tube. The variable H hot N+1 denotes the pressure drop produced in the N+1 segment of the hot pipe. The flow circulating through the N + 1 segment of the cold and hot pipes is denoted by Q N+1 . Let Hp loop N be the pressure drop produced in loop N, including the energy losses caused by the input valve and other hydraulic elements. The flow that circulates through loop N is denoted by q N . For every closed circuit formed by loops N and N+1 (see Fig. 3 ), the pressure drop for each segment can be calculated using the formulas (5) and (6) . It depends on the flow, the temperature of the segment, and its geometric dimensions (diameter and length). The dependence of the segment temperature is because the density and the dynamic viscosity of the oil depend on the working temperature [5] . ρf = 903 −0 . 672 T (7) μ= 1 . 41 ·10 −2 −1 . 6 ·10 −4 ·T + 6 . 41 ·10 −7 ·T 2 −8 . 66 ·10 −10 ·T 3 (8) The following set of equations must be fulfilled (system of Eqs. (9) ): Hp loop N = H p loop N+1 + H cold N+1 + H hot N+1 Q cold N+1 = Q cold N −q N Q hot N = Q hot N+1 + q N Hp N+1 = H p loop N −Hp loop N+1 −H cold N+1 −H hot N+1 (9) Notice that the flow of the first loop is imposed: the sum of all the loops’ flow is equal to the flow delivered by the pump. The complete system for the entire plant comprises 9 closed circuits. The equations can be solved iteratively using the iterative Algorithm 1 . Algorithm 1: Algorithm to solve the flow distribution. All the flows for the loops q N are considered to be distributed proportionally. T ol ← 10 6 . k ← 1 . while T ol > T ol max and k < k max do Set q 1 = Q pump − N=10 N=2 q N for N=1 to 9 do Solve problem 9 for the closed circuit formed by loops N and N+1. Save the value Hp N+1 to compute the flow q N+1 in the next iteration. end Compute the flow increment for the next iteration as q N+1 = K alg ·Hp N+1 . Update: q k +1 N+1 = q k N+1 + q k N+1 . T ol =  N=10 N=2 | q N | , k = k + 1 . end The constant K alg > 0 must be chosen carefully to ensure that the algorithm converges. Based on several tests carried out considering different conditions, it has been found that for the flow levels used in the plant (2–12 l/s), the algorithm converges to a correct solution if K alg ≤0 . 1 . A value K alg = 0 . 05 was chosen. For the rest of this work, it is assumed that the oil pump can provide the necessary pressure drop to function properly. The pump is capable of providing the pressure drop needed for the solar field, the steam generator, and the auxiliary pipes. P pump = P sol ar fiel d + P steamgenerator + P pipes (10) The following tests consider proportional valves (the flow entering the loop is proportional to its aperture). The total flow delivered by the pump is 9 l/s. The diameter of the cold and hot pipes is considered to be 0.09 m. Fig. 4 shows the flow distribution when a 40 m pipe is considered. The loops are distributed every 4 m. As can be seen, the flow distribution is very close to that obtained without considering the hydraulic model. The maximum error is about 1%. However, when the length of the pipe increases to 280 m, the flow distribution changes substantially, as seen in Fig. 5 . The separation between loops is 28 m. This problem becomes very important in commercial solar plants since their size is considerable and covers many hectares. The maximum error is about 6%, which is very important. Fig. 6 shows the flow distribution when the input valves of the loops (the valve apertures are given in percentage) are considered as follows: Ap v alv es (%) = [100 100 90 70 85 75 90 95 70 100] (11) As can be seen, the estimation error between considering the hydraulic model and not considering it can be up to 5% approximately. 4 A.J. Gallego, A.J. Sánchez, J.M. Escaño et al. European Journal of Control 67 (2022) 100717 Fig. 4. Flow distribution: hydraulic vs Non hydraulic model. Cold and hot pipes of 40 m. Fig. 5. Flow distribution: hydraulic vs Non hydraulic model. Cold and hot pipes of 280 m. Fig. 6. Flow distribution: hydraulic vs Non hydraulic model. Cold and hot pipes of 280 m. Valves not fully opened. 4. Control strategy for the solar plant In this section, the proposed Model Predictive Control (MPC) strategy is presented. The control objective of the plant is two-fold accomplished by two different algorithms with two different sampling times: 1 First, regulate the average temperature of the solar field around a set point. A non-linear model predictive control (NMPC) algorithm is used to achieve this goal. This algorithm computes the flow set point for the main pumps with a sampling time of 30 s. The non-linear model needs an estimation of the metalfluid temperature profiles and the optical efficiency. The entire solar field is modeled as an equivalent loop with average optical efficiency. An unscented Kalman filter (UKF) is used as an estimation algorithm. 2 The second objective is to minimize the temperature discrepancy between the loops by manipulating the input valves. The control algorithm designed to carry out this task aims at opening the input valves of the most efficient loops and closing the valves of the least efficient ones. The contribution of this paper is that the algorithm uses a hydraulic model of the solar field to consider the pressure drop in pipes. The algorithm is computed every 20 minutes. 4.1. Nonlinear MPC for control of the average temperature of the loops In this section, the non-linear MPC control strategy is presented. In general, the mathematical expression of the Model Predictive Control problem can be posed as follows (optimization problem (12) ). min u J = N p  t=1 y k + t| k −y ref k + t  y k + t| k −y ref k + t  + λ N c −1  t=0 u  k + t| k u k + t| k s.t. x k + t| k = F (u, x k + t−1 , x k + t−2 , . . . ) y k + t| k = C ·(x k + t| k ) u k + t| k = u k + t−1 | k + u k + t| k u min ≤u k + t| k ≤u max t = 0 , . . . , N p −1 (12) Where N p and N c stand for the prediction and control horizons, respectively. The parameter λpenalizes the control effort. Only u k ≡u k | k is applied to the system, recalculating the sequence every sampling time that corresponds to a receding horizon technique [29] . The difference between MPC strategies is the model used to predict the evolution of the output. In the case of a solar collector field, the state vector x corresponds to the metal-fluid temperature profile segments. The output y is the output temperature of the field, which corresponds to the last segment of the fluid temperature profile. The input is the manipulated variable, that is, the oil flow of the main pump q [22] . The vector Cis equal to [0 0 ... 1] . In the actual plant, only the inlet and outlet temperatures are measurable. The rest of the segment temperatures must be estimated by the UKF algorithm. Since the solar plant dynamics is highly non-linear, the use of conventional linear control strategies does not perform optimally in the entire range of plant operation [19,30,32] . In this paper, a non-linear model is used to predict the future evolution of the average temperature of the solar field. The function F used in the 5 A.J. Gallego, A.J. Sánchez, J.M. Escaño et al. European Journal of Control 67 (2022) 100717 MPC control problem to predict the evolution of the average temperature and the metal-fluid temperature profiles is a simplification of the PDE systems of Eqs. (1) . The model considers only 16 segments for metal-fluid temperatures instead of 172 segments. This simplification is essential to reduce the computational burden of the MPC problem. The integration step used to integrate the model is the same as that used for the full non-linear distributed parameter model. The model is integrated for a time of Np ∗30 seconds, but the output used in the cost function is picked every 30 seconds to compute its value. The main reason is that the control signal changes only every 30 seconds, and the MPC problem formulation is discrete. The non-linear MPC is solved using the fmincon function belonging to the Matlab optimization toolbox®. The parameters used are N p = 10 , N c = 4 and λ= 50 . The main reason to choose the control horizon smaller than the prediction horizon is that the computational time invested in the non-linear optimization problem depends on the number of decision variables. The values used here allow the optimization problem to be solved in less than 10 seconds, which is enough for the sampling time of 30 seconds. Only constraints in the maximum and minimum flow are considered. The possible values for the flow are within the range 2–12 l/s. 4.2. Model based optimization algorithm for the input valves This subsection describes the non-linear model-based optimization algorithm for manipulating the input valves. The algorithm presented here is based on that described in [33] and [35] . The difference and contributions of the algorithm presented here are listed below: 1. In this algorithm, a hydraulic model is used to compute the flow distribution depending not only on the opening valves but how far is the loop from the main pumps. The model is explained in Section 3 . 2. The optimization algorithm uses an unscented Kalman filter (UKF) to estimate the metal-fluid temperature profile instead of a CART algorithm used in [33] . The algorithm uses 16 segments for the fluid and the metal as done in [33] . It was found that this approach offered an adequate trade-off between precision and computational time. This approach avoids the need to solve a nonlinear optimization problem for estimating the optical efficiency of the loops since the UKF can perform the two tasks with one algorithm computed every 30 seconds. 3. The manipulation of the input valves aims at correcting the difference between the loop efficiencies in steady-state. Thus, the algorithm is computed every 20 minutes. If strong transients in solar radiation are affecting the plant, the temperature of the loops may not be regulated around a set-point. Furthermore, the outlet temperature of the loops decreases even with the minimum flow if solar radiation levels are low. These conditions produce that the model-based optimization algorithm sets all the valves to the minimum value, intending to increase the temperature. A significant discrepancy between loops can be produced when solar radiation levels recover. A heuristic is used to address this problem. The heuristic works as follows: the input valve aperture is only computed in stable operating conditions when the average temperature of the field is close to the desired reference, the derivative of the temperature is low, and the measured solar radiation is high enough. This heuristic is based on the experience obtained in commercial solar trough plants. The values for the heuristic have been tuned up by running several simulations under different operating conditions. The algorithm is only computed when the average temperature of all loops is close to the reference, the derivative of that temperature is less than 2 ◦C/min (stable operation), and the measured solar radiation is higher than 300 W / m 2 . If those conditions are not presented, input valve apertures are not computed. The mathematical formulation of the non-linear model-based optimization problem is given by the set of Eqs. (13) . The cost function penalizes the deviation of each loop temperature, in a steady-state, concerning the average value. Since the NMPC controls the average value, all loop temperatures will tend to the temperature reference. Where T Np n is the outlet temperature of loop n in steady-state. Ap(n ) is the aperture of the input valve of loop n, and T a v g is the average outlet temperature of the loops. F is the function that implements the evolution of the temperature of the loops. It corresponds to the system of Eqs. (1) but considers 16 segments to the metal-fluid profiles instead of 172. His the function that computes the flow distribution taking into account the aperture of the input valves, the length of the pipes, and the fluid temperatures. The function His the algorithm 1 that computes the flow distribution of the loops taking into account the hydraulic model explained in Subsection 3.2 . Ap min and Ap max are the minimum and maximum apertures allowed, respectively. T in is the HTF temperature at the main pump. min Ap J = N loops  n =1 T Np n −T a v g  T Np n −T a v g  s.t q n = H(Ap(n ) , T n (k −1) , L pipe , . . . ) T Np n = F (T n (k −1) , q n , I, T in , . . . ) T a v g = 1 N loops N loops  n =1 T Np n Ap min ≤Ap(n ) ≤Ap max (13) Notice that although the variables to be computed are the input valve apertures, the cost function uses only the temperature distribution in steady-state. The input valve apertures A p are used to calculate the flow distribution q n using the Hfunction. This flow distribution is then used to obtain the steady-state outlet temperature for the loops using the function F . These two non-linear models are used in the optimization problem so that the fmincon can minimize the cost function value. A simulation is carried out to show how the heuristic works. Fig. 7 shows a simulation in which a strong transient in solar radiation affects the field 13 to 13.4 h. If no heuristic is applied, the optimization algorithm reduces the input valves to a minimum in several loops to increase the temperature. When solar radiation recovers, a significant discrepancy in loop temperatures appear. It is because some loops receive less flow than others. At 14 h, the algorithm is computed, and the discrepancy is reduced. Fig. 8 shows the result when the heuristic is used. As can be seen, when the transient is produced, the algorithm is not computed, and the aperture of the valve remains constant, thus avoiding the undesirable effect produced in the previous case. Fig. 9 shows the evolution of solar radiation and oil flow during operation. The fall of solar radiation to 0 is produced from 13 to 13.4 h. 4.3. Estimation algorithm: The unscented kalman filter As mentioned above, the nonlinear optimization algorithm and the NMPC need both the metal-fluid temperature profiles and an estimation of the optical efficiency of the loops. Instead of using a CART-based algorithm to estimate the metal-fluid temperature profile and a nonlinear optimization problem to estimate the optical 6 A.J. Gallego, A.J. Sánchez, J.M. Escaño et al. European Journal of Control 67 (2022) 100717 Fig. 7. Temperature distribution where a strong transient affect from 13 to 13.4 h using the algorithm proposed in [33] . Top part: Average temperature (black solid line) and loops temperature (1–10). Bottom part: inlet valves aperture (1–10). Fig. 8. Temperature distribution where a strong transient affect from 13 to 13.4 h using the algorithm proposed in this paper where the heuristic is used. Top part: Average temperature (black solid line) and loops temperature (1–10). Bottom part: inlet valves aperture (1–10). Fig. 9. Test changing the average temperature reference with a strong transient: Evolution of flow and solar radiation. efficiency, a UKF is used as an estimator in this paper. The main reason is that this nonlinear estimator can estimate both metalfluid temperature profiles and optical efficiency with one algorithm, offering a more straightforward approach with similar performance. In this subsection, the actual values for the temperature profiles are those obtained from the full nonlinear distributed parameter model. The real value of the optical efficiency is that used in the model. The UKF is based on the unscented transformation, which represents a method to calculate the mean and covariance of a random variable that undergoes a nonlinear transformation [31,39] . The main advantage of the UKF algorithm concerning other approaches, such as the extended Kalman filter, is the direct use of the nonlinear model instead of a first- or second-order linear approximation. Several variations and improvements of the UKF estimation algorithm have been presented and proposed in the literature [21,41,42] , some of them related to the application of the UKF algorithm to solar power plants to estimate states and parameters [15] . Several research works related to the use of the UKF estimator with the model predictive control algorithm have been developed in the literature with good performance [18] . Fig. 10 shows a comparison between the optical efficiencies for loops 1,4, and 8 estimated by the UKF algorithm and those used by the full distributed parameter model (considered the real ones). The optical efficiencies estimated by the UKF converge to the actual values, even though the initial estimate is far from the real one. Fig. 11 shows the metal-fluid temperature profiles for loops 1,4, and 8 obtained by the UKF. As seen, the estimation obtained by the UKF is very close to the actual values obtained by the full nonlinear model. 4.4. Final control scheme Fig. 12 shows the overall control scheme. It works as follows: Every 30 seconds, the NMPC control algorithm receives the plant variables: effective solar radiation (product IK opt n o G ), ambient temperature, average metal-fluid temperature profiles, and average optical efficiency. Then it computes an oil flow value for the main pump. The UKF estimator algorithm estimates every 30 seconds the metal-fluid temperature profiles and the optical efficiency (which 7 A.J. Gallego, A.J. Sánchez, J.M. Escaño et al. European Journal of Control 67 (2022) 100717 Fig. 10. UKF efficiency estimation for loops 1,4 and 8. is unknown a priori) for every loop. It receives the plant variables and the outlet temperature of all loops. Every 20 minutes, the nonlinear optimization algorithm receives the plant variables: the metal-fluid temperature profiles and the optical efficiency of all loops. Using these estimates, it computes the value of the input valve apertures using the simplified distributed parameter model and the hydraulic model of the plant. The algorithm aims to reduce the discrepancy in the outlet temperature of all loops. 5. Simulation results In this section, simulation results that compare the proposed algorithm with that described in [33] are presented. Simulations show the importance of considering the hydraulic model in the optimization problem as the length of pipes increases. Fig. 11. UKF temperature profile estimation for loops 1,4 and 8. The minimum variation allowed for valve apertures is considered 1%, and valve apertures are limited to be within the range 70–100%. The inlet temperature of the solar field is considered to be the average temperature of the field minus a thermal jump of 80 ◦C as done in [6] . The set of optical efficiencies (loops 1 to 10) used in this paper is as follows: K opt = [0 . 513 0 . 543 0 . 543 0 . 57 0 . 57 0 . 582 0 . 582 0 . 61 0 . 61 0 . 623]  For the sake of clarity, the algorithm proposed in this paper is called HMBA (Hydraulic Model-Based Algorithm), and the algorithm proposed in [33] is called NHMBA (NonHydraulic Model- Based Algorithm). Fig. 12. Overall Control Scheme. 8 A.J. Gallego, A.J. Sánchez, J.M. Escaño et al. European Journal of Control 67 (2022) 100717 Fig. 13. Sim1: Test changing the average temperature reference. NHMBA performance.Top part: Average temperature (black solid line) and loops temperature (1– 10). Bottom part: inlet valves aperture (1–10). Fig. 14. Sim1: Test changing the average temperature reference. HMBA performance.Top part: Average temperature (black solid line) and loops temperature (1– 10). Bottom part: inlet valves aperture (1–10). Fig. 13 shows a test in which the temperature reference for the average temperature of the solar field changes throughout the day. The distance between the loops is 4 m; that is, the pipe connecting all the loops measures 40 m. The results correspond to the NHMBA algorithm. Fig. 15. Sim2: NHMBA performance when considering that the pipes length is 300 m. Top part: Average temperature (black solid line) and loops temperature (1– 10). Bottom part: inlet valves aperture (1–10). Fig. 14 shows the results using HMBA in this paper. As can be seen, both algorithms perform very similarly. They managed to steer all loop temperatures within a 1 ◦C range. When the pipe length increases, the difference between the two algorithms’ performance becomes more notorious. Fig. 15 shows the results of the NHMBA method when the length of the pipe is considered 300 m long. The distance between loops is 30 m long. The test conditions are the same as the previous one. The proportional flow model is less accurate than the hydraulic model due to energy losses in pipes. As can be seen, the difference between the maximum and minimum temperatures is approximately 4 ◦C when using the NHMBA algorithm. In current commercial solar trough plants, the pipes connecting the loops can have a km length. Considering a hydraulic model to predict the flow distribution becomes more important. Fig. 16 shows the result using the HMBA algorithm. The difference between the maximum and minimum loop temperatures is less than 1 ◦C. This effect is significant when the average temperature of the plant is close to the maximum allowable temperature. In the ACUREX field, when a loop outlet temperature reaches 300 ◦C, the collectors are defocused to avoid overheating situations. If there are much hotter loops than others, the collectors have to be defocused, thus losing thermal energy, which is undesirable. The defocusing procedure works by increasing the tracking angle of the collectors concerning the solar beam in the anti-clockwise direction. If the collector perfectly tracks the Sun, the defocus angle is nil. If the defocus angle is greater than 0, the incidence angle between the solar beam and the collector increases, thus decreasing the geometric efficiency [16] . The algorithm to defocus the mirrors implemented in the simulation works as follows: if the outlet temperature of the loop reaches 300 ◦C, the defocus angle [36] is gradually increased to ensure that the outlet temperature does not exceed 300 ◦C. The increased defocus angle decreases the geometric efficiency of the 9