Market-based clustering of model predictive controllers for maximizing collected energy by parabolic-trough solar collector fields

Applied Energy 306 (2022) 117936 Available online 13 October 2021 0306-2619/© 2021 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Applied Energy journal homepage: www.elsevier.com/locate/apenergy Market-based clustering of model predictive controllers for maximizing collected energy by parabolic-trough solar collector fields Eva Masero∗, José M. Maestre, Eduardo F. Camacho Departamento de Ingeniería de Sistemas y Automática, Universidad de Sevilla, C/ Camino de los Descubrimientos, s/n., 41092 Sevilla, Spain ARTICLE INFO Keywords: Concentrated solar power Nonlinear model predictive control Coalitional control Distributed solar collector field Solar thermal applications Thermal energy efficiency ABSTRACT This article focuses on maximizing the thermal energy collected by parabolic-trough solar collector fields to increase the production of the plant. To this end, we propose a market-based clustering model predictive control strategy in which controllers of collector loops may offer and demand heat transfer fluid in a market. When a transaction is made between loop controllers, a coalition is formed, and the corresponding agents act as a single entity. The proposed hierarchical algorithm fosters the formation of coalitions dynamically to improve the overall control objective, increasing the thermal energy delivered by the field. Finally, the proposed controller is assessed via simulation with other control methods in two solar parabolic-trough fields. The results show that the energy efficiency with the clustering strategy outperforms by 12% that of traditional controllers, and the method is implementable in real-time to control large-scale solar collector fields, where significant gains in thermal collected energy can be obtained, due to its scalability. 1. Introduction The increasing trend of energy demand causes a significant impact on the environment. Efficient, clean, and secure renewable energy sources such as solar, hydropower, wind, biomass, and geothermal are essential to deal with this steep demand. Indeed, renewable energy accounted for nearly 28% of global electricity generation in the first quarter of 2020 and is projected to increase to almost 30% in 2021 [1, 2]. Since the most abundant and cost-effective renewable source in many countries is the sun, solar power has been widely researched in the last decades [3,4]. The most standard solar technologies are photovoltaic cells (PV), which directly generate electricity from sunlight, and concentrating solar power (CSP), which concentrates sun radiation in a heat transfer fluid (HTF) to produce steam and drive an electricity generator (see Fig. 1). For large-scale power generation, CSP presents further advantages due to storing thermal energy to produce electricity when there is low or no sunlight, e.g., cloudy days and nights. This work focuses on parabolic-trough collector fields, which belong to the family of CSP systems [5]. Other members are linear fresnel collectors [6], parabolic dishes [7], and solar power towers [8]. Our challenge here is to maximize the energy production despite the daylight discontinuity by manipulating the inlet valves of the collector loops. Numerous control techniques for solar collector fields are presented in the exhaustive surveys [9,10]. A remarkable control strategy with ∗Corresponding author. E-mail addresses: [email protected] (E. Masero), [email protected] (J.M. Maestre), [email protected] (E.F. Camacho). widespread use in the industry and academia is model predictive control (MPC) because it handles complexities such as delays, constraints, and unstable, non-minimum phase, and multivariable systems [11]. Many applications of MPC in solar collector fields can be found in the literature. For example, [12] proposes an adaptive MPC with a Kalman filter to estimate the solar radiation and temperature profiles, [13] details an MPC scheme for optimal scheduling, and [14] presents a dual MPC based on a linear model for tracking and disturbance rejection. The nonlinear behavior of distributed solar collector fields also promotes nonlinear model predictive strategies (NMPC). For instance, [15] proposes an NMPC with guaranteed stability, and studies as [16,17] present NMPC methods that deal with dead-times and nonlinear dynamics. Most of the studies mentioned above optimize the global flow while tracking a reference for the outlet temperature of the solar collector fields. Nonetheless, these centralized strategies may fail to manage large-size fields due to the significant computational times. For instance, in a parabolic-trough solar collector plant like Solana, which has around 800 loops, centralized control cannot be performed due to the significant number of decision variables. For that reason, distributed control methods are more suitable to deal with these issues [18] because they divide the system into subsystems controlled by local agents. In a distributed fixed cooperation structure, local controllers with partial system information cooperate to decide what control action https://doi.org/10.1016/j.apenergy.2021.117936 Received 9 June 2021; Received in revised form 21 September 2021; Accepted 22 September 2021 Applied Energy 306 (2022) 117936 2 E. Masero et al. Fig. 1. Schematic of a parabolic-through solar collector plant. to implement. Various studies in this context propose distributed MPC strategies where it is assumed that there are valves at the beginning of each loop and can be controlled individually. For example, a noncooperative distributed MPC is proposed in [19], where each loop controller computes its HTF flow individually; [20] details a logicbased distributed MPC to approach the performance of a centralized controller while reducing the computational load; and [21] studies the manipulation of loop valves to homogenize the temperature of the solar field. There are several challenges regarding the application of distributed control strategies in commercial plants. Firstly, valves are not currently employed to control their plants but for balancing hydraulically the flow of loops in the field. Secondly, the fixed cooperation structure between loop controllers cannot deal with the changing operation conditions, e.g., high solar irradiance differences between loops caused by clouds covering partially large-scale solar fields. A more sophisticated approach is coalitional model predictive control [22,23], which arranges local agents in time-varying cooperating clusters to attain a trade-off between the cooperation effort and performance. This approach can be seen as a dynamically changing system partitioning that deals with the varying conditions of the system operation [24,25]. From this viewpoint, partitions can be computed in a fully distributed manner or using a hierarchical architecture. An example of the former is presented in [26], which uses the Shapley value to find the best partition. Regarding the latter, [27] proposed a two-layer architecture in which local agents compute their control inputs using MPC in the lower layer, and, in the upper one, a supervisor checks and re-balances the inputs as needed to reject disturbances. Potential applications of this approach include microgrids [28], chemical plants [29], robot fleets [30], traffic networks [31], supply chains [32], telecommunication networks [33], wind farms [34], irrigation canals [22], water networks [35], automated vehicles [36], and also parabolic-trough solar collector plants [37]. In this article, we propose a hierarchical market-based coalitional MPC to maximize the thermal energy collected by distributed solar parabolic-trough plants. The upper control layer periodically executes an NMPC controller combined with a market mechanism to select the partition that maximizes the overall utility. In the lower layer, the resulting coalitions work independently to solve their NMPC subproblems, which maximize the thermal energy production. Another objective of the proposed coalitional approach is to simplify the search of the field partition making feasible to control large-scale solar collector plants. Therefore, the contributions of this paper are: •A market-based scheme for selecting the field partition online. We propose to let local controllers offer and demand HTF in a market to mitigate the effects derived from clouds, dusty collectors, and maintenance issues. Once fluid transactions are settled, the agents involved in the trades form coalitions and operate as single entities. This method is inspired by the electricity intraday market [38,39], where prices are set dynamically, and auction-based strategies as those proposed in [40,41]. •Improvements of thermal energy efficiency by up to 12% by implementing the proposed coalitional MPC strategy in two simulated solar collector fields of 10 and 100 loops modeled after the ACUREX field at the Plataforma Solar de Almería (PSA) [16,12]. •Feasible real-time implementation in large-scale fields applying the proposed method due to its superior scalability. Since real plants use a sampling time of 60 s, it would be feasible to control fields of up to 300 loops approximately with our current implementation. •The payback period on software investments for controlling valves in a decentralized way will be between one and six months, depending on the plant’s energy production. Index of contents. Section 2introduces the nonlinear model and constraints of distributed parabolic-trough solar collector fields. Section 3 explains and illustrates the market mechanism to determine the cooperation partition. Section 4describes the control objective and the proposed algorithm. Details of the case studies and results from the coalitional control regarding other strategies are illustrated in Section 5. Conclusions are summarized in Section 6. Finally, a brief design guide for the main parameters of the controller is given in Appendix. 2. Parabolic-trough solar collector field model This section presents the model of a parabolic-trough solar collector field as the one depicted in Fig. 1. It is composed of concave mirrors that concentrate the solar direct normal irradiance (DNI) on a tube located in its focal line. The HTF that circulates through the tube absorbs the thermal energy and is sent to the power conversion system (PCS), where superheated steam is produced to drive turbines and generators. Likewise, the HTF can also be sent to a thermal energy storage (TES) system in some plants for later use. 2.1. Nonlinear system dynamics The parabolic-through solar collector field can be disaggregated into a set of parallel loops  ≜ {1,…, 𝑁loop}, which are split into  ≜ {1,…, 𝑁seg}segments of length 𝛥𝑙, as represented in the field of Fig. 1. The system dynamics can be described by a distributed parameter model of the temperature composed of partial differential equations (PDEs) that depend on time (𝑡) and space (𝑙), describing the energy balance [9]: 𝜌m𝐶m𝐴m𝜕𝑇 m 𝜕𝑡 (𝑡, 𝑙) = 𝜂col(𝑙)𝐺(𝑙)𝐼(𝑡, 𝑙) − 𝜋𝐷m𝐻l(𝑇m(𝑡, 𝑙) − 𝑇a) −𝜋𝐷f𝐻t(𝑡, 𝑙)(𝑇m(𝑡, 𝑙) − 𝑇f(𝑡, 𝑙)), (1) 𝜌f𝐶f𝐴f𝜕𝑇 f 𝜕𝑡 (𝑡, 𝑙) + 𝜌f𝐶f𝑞(𝑡)𝜕𝑇 f 𝜕𝑡 (𝑡, 𝑙) = 𝜋𝐷f𝐻t(𝑡, 𝑙)(𝑇m(𝑡, 𝑙) − 𝑇f(𝑡, 𝑙)),(2) where superscripts m,f, and arepresent the metal tube, the HTF, and the ambient, respectively. The model parameters and variables are summarized in Table 1. The PDEs can be approximated using the forward Euler’s method as: 𝜕𝑇 𝜕𝑡 (𝑡, 𝑙) ≈ 𝑇(𝑡+𝛥𝑡, 𝑙) − 𝑇(𝑡, 𝑙) 𝛥𝑡 ,𝜕𝑇 𝜕𝑙 (𝑡, 𝑙) ≈ 𝑇(𝑡, 𝑙 +𝛥𝑙) − 𝑇(𝑡, 𝑙) 𝛥𝑙 ,(3) so that Eqs. (1) and (2) can be solved by a two-step iterative process with finite differences as proposed by [3]: Applied Energy 306 (2022) 117936 3 E. Masero et al. Table 1 Summary of model parameters and variables. Symbol Description Unit 𝜂col Efficiency of collectors – 𝐺Aperture of collectors m 𝐼Direct solar irradiance W∕m2 𝐷mOutside diameter of the pipe m 𝐷fInside diameter of the pipe m 𝐻tCoef. of transmission metal-fluid W∕(m2◦C) 𝐻lCoef. of thermal losses W∕(m2◦C) 𝑞Flow rate l/s 𝛥𝑡 Discretization time of the model s 𝛥𝑙 Length of segments m 𝜌m, 𝜌fDensities kg∕m3 𝐶m, 𝐶fSpecific heat capacities J∕(kg ◦C) 𝐴m, 𝐴fCross-sectional areas m2 𝑇m, 𝑇 f, 𝑇 aTemperatures ◦C •The metal tube and fluid temperature are computed at the first stage, considering the fluid to be steady-state: 𝑇m 𝑖,𝑗 (𝑘) = 𝑇m 𝑖,𝑗 (𝑘− 1) + 𝛥𝑡 𝜌m𝐶m𝐴m(𝜂col 𝑖,𝑗 𝐺𝑖,𝑗 𝐼𝑖,𝑗 (𝑘) −𝜋𝐷m𝐻l 𝑖,𝑗 (𝑘− 1)(𝑇m 𝑖,𝑗 (𝑘− 1) − 𝑇a(𝑘)) −𝜋𝐷f𝐻t 𝑖,𝑗 (𝑘− 1)(𝑇m 𝑖,𝑗 (𝑘− 1) − 𝑇1f 𝑖,𝑗 (𝑘− 1))), 𝑇f 𝑖,𝑗 (𝑘) = 𝑇1f 𝑖,𝑗 (𝑘− 1) + 𝜋𝐷f𝐻t 𝑖,𝑗 (𝑘− 1)𝛥𝑡 𝜌f 𝑖,𝑗 (𝑘− 1)𝐶f 𝑖,𝑗 (𝑘− 1)𝐴f(𝑇m 𝑖,𝑗 (𝑘− 1) − 𝑇1f 𝑖,𝑗 (𝑘− 1)), (4) where segment 𝑖∈and loop 𝑗∈represent the spatial dependence (𝑙), and 𝑘denotes the discrete-time instant 𝑡=𝑘𝛥𝑡. •At the second stage, the fluid temperature is corrected using the net energy transported by the fluid: 𝑇1f 𝑖,𝑗 (𝑘) = 𝑇f 𝑖,𝑗 (𝑘) − 𝑞𝑗(𝑘)𝛥𝑡 𝛥𝑙 𝐴f(𝑇f 𝑖,𝑗 (𝑘) − 𝑇f 𝑖−1,𝑗 (𝑘)).(5) 2.2. Inlet and outlet temperature The inlet temperature of the initial segments of each loop 𝑗∈is considered to be equal to the inlet temperature of the field, i.e., 𝑇f 1,𝑗 (𝑘) = 𝑇in(𝑘)for all 𝑘. Note that the inlet temperature 𝑇in(𝑘), the corrected solar DNI 𝐼𝑖,𝑗 (𝑘), and the ambient temperature 𝑇a(𝑘)are considered as disturbances that can be measured or estimated at instant 𝑘from the control viewpoint. Additionally, the outlet temperature of the solar field is computed considering the temperature of the final segments of the loops as 𝑇out (𝑘) = ∑𝑗∈𝑇f 𝑁seg,𝑗 𝑞𝑗(𝑘) 𝑞T(𝑘),(6) where 𝑞T(𝑘) = ∑𝑗∈𝑞𝑗(𝑘)is the total flow rate. 2.3. Operational constraints The solar collector field presents several operational constraints that need to be considered. The first one is the flow rate of each loop 𝑗∈, which is bounded by 𝑞min 𝑗≤𝑞𝑗(𝑘)≤𝑞max 𝑗,(7) where the minimum flow rate 𝑞min 𝑗>0is based on the Reynolds number to guarantee a turbulent flow and, therefore, achieve homogeneous HTF temperature along the pipe. The maximum value 𝑞max 𝑗depends on the maximum allowed pressure drop. The total flow 𝑞T(𝑘)is also upper bounded by the maximum total flow rate 𝑞max T. Another constraint is considered in each loop 𝑗∈to keep the overall outlet temperature 𝑇out within its operational limits: 𝑇f,min ≤𝑇f 𝑁seg,𝑗 (𝑘)≤𝑇f,max.(8) where 𝑇f,min and 𝑇f,max are, respectively, the minimum and the maximum temperature of the fluid. Since the outlet HTF temperature of each loop depends on the HTF flow, the constraints on the maximum operating temperatures may impose some limitations on the HTF flow in each loop. These limitations also rely upon the DNI received by the loop and its reflectivity, which depends on the maintenance received because recently cleaned loops will have much higher reflectivity. 3. HTF flow market Our objective is to maximize the resulting thermal power while minimizing control efforts and temperature constraint violations despite disturbances such as the changes in the solar irradiance due to the clouds and the reflectivity differences between the loops. To this end, the HTF flow of each loop is locally controlled by an agent that can cooperate with other agents to reduce its local cost and contribute to the global objective. However, satisfying the coupling constraints in a non-centralized fashion requires setting local constraints that split shared resources as the HTF. Consequently, some loops might improve their performance by increasing their flow over their limits, while others may reach the optimum without activating their constraints. Therefore, we proposed to organize the agents’ supply and demand of HTF in a market following the idea of the electricity intraday market and auction-based mechanisms [39,40]. 3.1. Field partitioning The solar collector field can be considered a cooperative network described by an undirected graph 𝐺≜(,)with a set of loops and a set of cooperating links ⊆×. Each loop 𝑗∈is locally managed by a controller with local information, and each enabled link 𝑒∈ establishes a bidirectional information flow between the corresponding loop controllers. Definition 1 (Field Partition).A partition (𝑘)of the graph is formed by a set of connected, non-overlapping, non-empty clusters of loops (the so-called coalitions): (𝑘)≜{1,…,𝑁},(9) satisfying ⋃∈(𝑘)=. Note that the size of a coalition can vary from a single loop to the entire solar collector field. Given partition (𝑘)at instant 𝑘, the thermal power of each coalition is defined as 𝑊(𝑘) = 𝑊out (𝑘) − 𝑊in (𝑘) = ∑ 𝑗∈ 𝑊out 𝑗(𝑘) − ∑ 𝑗∈ 𝑊in 𝑗(𝑘),(10) where the output and input thermal powers of each loop 𝑗∈are, respectively, 𝑊out 𝑗(𝑘) = 𝜌f 𝑁seg,𝑗 (𝑇f 𝑁seg,𝑗 (𝑘))𝐶f 𝑁seg,𝑗 (𝑇f 𝑁seg,𝑗 (𝑘))𝑞𝑗(𝑘)𝑇f 𝑁seg,𝑗 (𝑘), 𝑊in 𝑗(𝑘) = 𝜌f 1,𝑗 (𝑇f 1,𝑗 (𝑘))𝐶f 1,𝑗 (𝑇f 1,𝑗 (𝑘))𝑞𝑗(𝑘)𝑇f 1,𝑗 (𝑘). The coalition flow rate 𝑞(𝑘)is used as a control input, assuming that the hydraulic dynamics are considerably faster than the thermal dynamics. Moreover, the outlet temperature of any coalition is 𝑇out (𝑘) = ∑𝑗∈𝑞𝑗(𝑘)𝑇f 𝑁seg,𝑗 𝑞(𝑘),(11) where the coalition flow rate 𝑞(𝑘) = ∑𝑗∈𝑞𝑗(𝑘)by aggregating the HTF flow rates of all loops grouped in . Moreover, it is lower and upper delimited by 𝑞min =∑𝑗∈𝑞min 𝑗and 𝑞max =∑𝑗∈𝑞max 𝑗, respectively. Applied Energy 306 (2022) 117936 4 E. Masero et al. 3.2. Control objective The control goal of any coalition ∈(𝑘)is to compute the flow rate sequence in the control horizon 𝑁u 𝑄(𝑘)≜[𝑞(𝑘), 𝑞(𝑘+ 1),…, 𝑞(𝑘+𝑁u− 1)], to minimize the cost function: 𝐽(⋅) = 𝑁p ∑ 𝑛=1 (−𝑊(𝑘+𝑛) + 𝛼𝛹(𝑇f 𝑁seg,𝑗 (𝑘+𝑛))) +𝛽 𝑁p−1 ∑ 𝑛=0 ∑ 𝑗∈(𝑞𝑗(𝑘+𝑛) − 𝑞𝑗(𝑘+𝑛− 1))2, (12) where 𝑁pis the prediction horizon,1(𝑘+𝑛)represents the predicted value 𝑛steps ahead of 𝑘,𝛼, 𝛽∈R>0are tuning parameters, and 𝛹(𝑇f 𝑁seg,𝑗 (𝑘+𝑛)) =∑ 𝑗∈ max (𝑇f 𝑁seg,𝑗 (𝑘+𝑛) − 𝑇f,max 𝑇f,max , 𝑇f,min −𝑇f 𝑁seg,𝑗 (𝑘+𝑛) 𝑇f,max ,0)2 . (13) By minimizing 𝐽, we maximize the thermal power delivered in the prediction horizon while minimizing the control efforts and temperature violations. Notably, the second and the third terms of (12) are soft constraints weighted by the tuning parameters 𝛼and 𝛽. The second term, which corresponds to (13), penalizes the maximum quadratic violation of the outlet temperature constraints of coalition . The third term of (12) penalizes the variations in the inputs to reduce the flow rate oscillation and the wear of valves during daily operation. Remark 1. The control problem of each coalition can be solved by a loop that works as a leader, e.g., the smallest numbered loop, or distributed among the loops in the coalition [42]. 3.3. Market mechanism The set of loops can be divided into two disjoint subsets: 𝑠and 𝑑of supply and demand loops, respectively. To gain performance, •any supply loop 𝑠∈𝑠can decrease 𝑠𝑢𝑝𝑠∈R>0its assigned HTF flow, •any demand loop 𝑑∈𝑑can increase 𝑑𝑒𝑚𝑑∈R>0its assigned HTF flow. The total demand and supply of HTF flow are computed as 𝑠𝑢𝑝T(𝑘) = ∑ ∈(𝑘)∑ 𝑠∈∩𝑠 𝑠𝑢𝑝𝑠(𝑘), 𝑑𝑒𝑚T(𝑘) = ∑ ∈(𝑘)∑ 𝑑∈∩𝑑 𝑑𝑒𝑚𝑑(𝑘), where both are constrained by 0as the lower limit, and 𝑞max T−∑𝑗∈𝑞min 𝑗 as the upper limit. The challenge is to find the best partition (𝑘)of the field to deal with the disturbances at each time instant 𝑘. We propose an auction-based mechanism inspired by the electricity intraday market that allows agents to offer and demand HTF flow. In particular, the HTF transactions are swapped immediately between agents based on their utility gain to HTF flow changes, with the utility function defined as: 𝑈(⋅)=−𝐽(⋅).(14) The market mechanism goes as follows: 1See Appendix for further details. Fig. 2. Example of five demand and five supply loops in an HTF flow market, which are represented in solid red and dashed blue lines, respectively. 1. Each coalition initially contains one agent = {𝑗}, with its utility gain computed using (14) when its HTF flow is increased/decreased a quantum 𝛥𝑞. Those with higher utility gains when their flow is increased contain demand agents 𝑑; otherwise, they contain supply agents 𝑠. 2. The HTF transactions take place in coalitions of two agents following an auction-based mechanism in which the supplier with the highest utility gain shares its surplus of HTF with the highest bidder. Afterward, the supplier with the second-highest utility gain is grouped with the second-highest bidder, and so forth. Therefore, new coalitions are formed by merging entities starting from top-winning pairs to least-losing pairs. 3. Let 𝑈𝑑be the utility of the demand agent and 𝑈𝑠denote that of the supplier. The resulting entity has a utility 𝑈computed by (14), providing surpluses (𝑈−𝑈𝑑)and (𝑈−𝑈𝑠)to the demand and supply agent, respectively. If any surplus is below a given threshold 𝜑∈R0+, we consider that it is not worth forming a coalition so as to reduce the computational burden. The resulting coalitions can also be classified again into supply and demand entities repeating this procedure 𝑀∈N+times to find a better partition (𝑘) = 𝑀, subject to a maximum coalition size constraint ||max ∈N+. Note that the bigger a coalition is, the better its performance becomes but at the expense of increasing the computational effort. For that reason, 𝑀and ||max need to be properly set to find a trade-off between performance improvement and computational burden. Remark 2. The set 𝑑(or 𝑠) may be empty if all agents want to increase (or decrease) their HTF flows. In this situation, agents can also be sorted attending to their utility gains. Coalitions of a feasible size can be formed by grouping agents with low and high utility gains if surpluses greater than 𝜑can be obtained. Remark 3. There are no physical HTF exchanges between loops. That is, the HTF transactions are more related with changes in the local restrictions of the problems than with a physical transfer of HTF. The flow is a shared resource and the market mechanism improves the efficiency of the distribution. Example 1 (HTF Flow Market).Let us consider a solar collector field of 𝑁loop = 10 loops with |𝑑|= 5 demand loops and |𝑠|= 5 supply loops. Applied Energy 306 (2022) 117936 5 E. Masero et al. In Fig. 2, the demand curve is sorted in descending order according to the utility gains of those loops that are willing to increase their HTF. Thus, these utility gains can be interpreted as the maximum amount of utility that these entities are willing to trade to get an additional quantum of HTF flow. Likewise, the supply curve is arranged in increasing order with regard to the utility loss. The value depicted in the figure corresponds to the compensation in utility units that the suppliers require to reduce their local HTF flow. Therefore, it can be interpreted as the price measured in utility units that is demanded to provide the corresponding quantum of HTF flow. A measure of each trade is represented by the purple dotted lines of Fig. 2, which is computed as the arithmetic mean between the utility gains demanded by HTF flow suppliers and the utility offered by demand agents. Note that this could be seen as the average price of the transaction. As shown in Fig. 2, the demand loop 𝑑1offers 280 utility units for a quantum of HTF, whereas supply loop 𝑠1demands, at least, 20 utility units to share its quantum. Therefore, the HTF flow quantum is traded generating an overall a surplus of 260 utility units, which represents an average agent surplus of 130 utility gains. Finally, the resulting partition is 1= {1,2,3,4,5} = {{𝑑1, 𝑠1},{𝑑2, 𝑠2},{𝑑3, 𝑠3},{𝑑4, 𝑠4},{𝑑5, 𝑠5}}. 4. Coalitional algorithm The flow rate of the field can be non-uniformly distributed due to reflectivity differences between loops, giving rise to a loss of performance. To overcome this imbalance, we allow agents to offer and demand HTF flow in a spot market. Therefore, it is proposed a heuristic, hierarchical coalitional NMPC based on a market mechanism for clustering loops in such a way that the resulting energy is maximized. Details of the proposed algorithm are specified hereunder. 4.1. Upper control layer This layer executes Algorithm 1 to decide the new partition every 𝑇up ∈N+time instants, according to the available information (system state, disturbances). Recall that the procedure to obtain the partition can be repeated 𝑀iterations, generating larger coalitions with increased performance and computational burden. Algorithm 1 (Upper control layer). Output: (𝑘) if mod (𝑘, 𝑇up)=0then: 0: Set 𝑀,𝜑, and ||max. 1: Estimate the solar DNI in the prediction horizon (e.g., following [43,44], and measure/estimate temperature 𝑇f 𝑖,𝑗 (𝑘) ∀𝑖∈ , and flow rate 𝑞𝑗(𝑘−1) for any 𝑗∈(e.g. based on 𝑞T(𝑘−1) and the valve position via a hydraulic model). 2: Set initial coalitions as singletons containing one agent, i.e., 𝑗= {𝑗}, and the initial partitioning 0=⋃𝑗∈𝑗. for 𝑚= 1 to 𝑀iteration do: 3: Evaluate the utility gain of each coalition using (14) and classify it as a supply or demand entity. 4: Sort demand and supply coalitions in descending order regarding their utility gains. 5: Merge one supply and one demand coalitions if the resulting surpluses exceed the threshold 𝜑and the resulting size is below ||max. 6: Partition 𝑚is composed of the resulting coalitions. end for 7: After 𝑀iterations, the final partition (𝑘) = 𝑀is sent to the lower control layer. end if 4.2. Lower control layer After the selected partition is received from the upper control layer, coalitions compute their optimal control actions and implement them. The steps in the lower control layer are detailed in Algorithm 2. Algorithm 2 (Lower control layer). At each time instant 𝑘: 1: Each coalition ∈(𝑘)receives its forecasted solar DNI in the prediction horizon 𝑁p, and measures/estimates 𝑇f 𝑖,𝑗 (𝑘),∀𝑖∈of any loop 𝑗∈. 2: For each ∈(𝑘), compute the optimal flow rate of each loop belonging to coalition : 𝑄∗ (𝑘) = arg min 𝑄(𝑘)𝐽(⋅),(15) subject to the dynamics (4)–(5), the constraints (7)–(8), and 𝑞max =𝑞max T|| 𝑁loop , considering that the demand loops of coalition can benefit from the extra flow rate of supply loops. 3: Each coalition ∈(𝑘)applies the first element of its optimal input sequence 𝑄∗ (𝑘), and the rest is discarded. 4: Set 𝑘←𝑘+ 1. 5. Simulation benchmarks Two solar parabolic-trough fields are described here: (i) the ACUREX field, which has ten loops, and (ii) its ten-time extension, i.e., a 100-loop field. Both fields are simulated to analyze the advantages of the proposed coalitional algorithm in terms of scalability and performance for the following methods: 1. No-valves Track denotes a PI controller, which manages the whole field (without valves) to hold the outlet temperature in the set-point 𝑇out ref = 280◦C. 2. No-valves represents a centralized MPC method that considers that loops cannot be controlled independently, and thus, the same flow rate is applied to all of them. Its cost function is (12), being =. 3. Cen refers to a centralized MPC where the controller manages the whole field =. Its cost function is (12). 4. Dec denotes a decentralized MPC with cost function (12) where each loop 𝑗∈is managed by a local controller with no information about the rest of loops, i.e., = {𝑗}for all ∈ (𝑘). Note that the maximum flow in each loop becomes 𝑞max 𝑗= 𝑞max T∕𝑁loop in this case. 5. Coal. Pair indicates the loop-pair clustering algorithm proposed in [37], which promotes coalitions of two loops to balance the field, but the upper control layer is here established to be executed every 𝑇up = 5 time instants. 6. Coal. Market represents the proposed coalitional MPC method. The partition is updated every 𝑇up = 5 time instants by the upper layer after 𝑀= 2 iterations, with parameters 𝛥𝑞 = 0.5l/s and ||max = 3 selected due to their good trade-off between performance and computational burden in the simulations ( Appendix details the procedure performed to tune these parameters). Applied Energy 306 (2022) 117936 6 E. Masero et al. Table 2 Model parameters and constraints of the collector field. Symbol Value Units 𝜌m7800 kg∕m3 𝐶m550 J∕(kg ◦C) 𝐴f7.55 ⋅10−4 m2 𝐷m0.031 m 𝐷f0.0254 m 𝑞min 𝑗0.2l/s 𝑞max 𝑗1.5l/s 𝑞max T6.5l/s 𝑇a25 ◦C 𝑇f,min 220 ◦C 𝑇f,max 300 ◦C These control methods are compared in terms of the average energy produced (  𝐸) for the simulation time 𝑁sim. Another key feature is the maximum computational time, which is computed as the maximum time needed per coalition to solve problem (15) plus the time needed to compute the partition of the field divided among the number of coalitions, i.e., 𝜏= max (𝜏max (𝑘) + 𝜏(𝑘) |(𝑘)|),∀𝑘∈𝑁sim (16) where 𝜏max (𝑘)is the maximum time to solve the optimization problem of all coalitions at time instant 𝑘,𝜏(𝑘)is the time required to compute the partition, and |(𝑘)|indicates the number of coalitions in partition (𝑘). 5.1. Description of solar collector fields The distributed collector field of ten loops considered is ACUREX, which was located at the Plataforma Solar de Almería (PSA) [16]. This field was composed of East–West aligned single axis parabolic-trough collectors. There were 𝑁loop = 10 parallel collector loops of length 𝑑loop = 174 m, which were discretized into 𝑁seg = 174 segments of length 𝛥𝑙 = 1 m. Each loop 𝑗∈ {1,…, 𝑁loop}had 12 collectors connected in series, being the active part (that receiving solar irradiance) 144 m long, and the passive part 30 m. The HTF in ACUREX was the Therminol 55 thermal oil, whose density and its specific heat capacity were, respectively, 𝜌f 𝑖,𝑗 (𝑇f 𝑖,𝑗 (𝑘))= 903 − 0.672𝑇f 𝑖,𝑗 (𝑘), 𝐶f 𝑖,𝑗 (𝑇f 𝑖,𝑗 (𝑘))= 1820 + 3.478𝑇f 𝑖,𝑗 (𝑘). The coefficients of metal-fluid transmission and thermal losses were, respectively, defined as 𝐻t 𝑖,𝑗 (𝑞𝑖(𝑘), 𝑇 f 𝑖,𝑗 (𝑘))=𝑞0.8 𝑖(𝑘)𝛷(𝑇f 𝑖,𝑗 (𝑘)), 𝐻l 𝑖,𝑗 (𝑇f 𝑖,𝑗 (𝑘), 𝑇 a(𝑘))= 0.00249(𝑇f 𝑖,𝑗 (𝑘) − 𝑇a(𝑘))− 0.06133, where 𝛷(𝑇f 𝑖,𝑗 (𝑘))is defined as follows, in accordance with [3]: 𝛷(𝑇f 𝑖,𝑗 (𝑘)) = 2.17 ⋅106− 5.01 ⋅104𝑇f 𝑖,𝑗 (𝑘)+4.53 ⋅102𝑇f 𝑖,𝑗 (𝑘)2 − 1.64𝑇f 𝑖,𝑗 (𝑘)3+ 2.1⋅10−3𝑇f 𝑖,𝑗 (𝑘)4. A summary of the main model parameters and constraints is shown in Table 2. In this case, the ambient temperature 𝑇a(𝑘)is considered as a constant for simplicity. The outlet temperature of the steam generator is considered to have a temperature drop of 90◦C with respect to the inlet temperature and a time constant of 600 s. That is, 𝑇in(𝑠)  𝑇out (𝑠)=1 600𝑠+ 1 ,(17) where  𝑇out (𝑠) = 𝑇out (𝑠) − 90. Applying the z-transform to (17) with a discretization step 𝛥𝑡 = 0.5s, the resulting discrete-time equation of inlet temperature is 𝑇in(𝑘)=0.999167 𝑇in(𝑘− 1) + 8.33 ⋅10−4(𝑇out (𝑘− 1) − 90).(18) The second collector field considered is composed of 𝑁loop = 100 loops. Since it is an extension of ACUREX, the previous equations are held. The constraints and model parameters are equal to these of Table 2 except for the maximum total flow rate, which becomes ten times larger, i.e., 𝑞max T= 65 l/s. 5.2. Design of the parameters and defocus mechanism All simulations have been performed in Matlab®version R2020a with a PC Intel®Core™i7-8700 CPU at 3.20 GHz and 16 GB RAM. Algorithms 1 and 2 have been implemented using the solver fmincon from the Optimization Toolbox with the option: interior-point. The following simplifications are also considered: (i) the coefficient of thermal losses 𝐻lis pre-computed for a given set of flow rates and temperatures, and (ii) the sampling time of the controller 𝛥𝑡ctr is selected as a multiple of the discretization time of the model (i.e., 𝛥𝑡ctr =𝛾𝛥𝑡 with 𝛾∈ N+). For all MPC strategies, the prediction and control horizons are, respectively, 𝑁p= 12 and 𝑁u= 10 samples. These values have been selected by trial and error to reach the maximum performance with reasonable computational effort. To keep the outlet temperature within its bounds, tuning parameters for the 10-loop and 100-loop fields have been selected as 𝛼𝑗= 5⋅10−4 and 𝛼𝑗= 8⋅10−3, respectively. Furthermore, a simple defocus mechanism is employed to prevent HTF degradation from exceeding the maximum temperature 𝑇f,max. To this end, mirror defocusing is applied in the last six collectors of each loop as needed. Hence, the efficiency of the last six collectors 𝜂𝑐𝑜𝑙, which depends on the defocusing angle 𝜃as described in [45], have been made dependent on loop temperature as follows: Defocusing algorithm. 𝜂col 0=𝜂col if 𝑇f 𝑁seg,𝑗 >290◦Cthen 𝜃= 1.6◦,𝜂col ←0.75 𝜂col 0 end if if 𝑇f 𝑁seg,𝑗 >292◦Cthen 𝜃= 2.25◦,𝜂col ←0.5𝜂col 0 end if if 𝑇f 𝑁seg,𝑗 >294◦Cthen 𝜃= 2.75◦,𝜂col ←0.25 𝜂col 0 end if if 𝑇f 𝑁seg,𝑗 >295◦Cthen 𝜃= 5◦,𝜂col ←0𝜂col 0 end if where 𝜂col 0corresponds to the collector efficiency when it is defocused. The tuning parameter 𝛽𝑗= 3 has been chosen to weigh the valve control effort for both fields, considering that a small and large value will lead to more and less flow rate oscillation in all the valves during daily operation, respectively. In the control model, 𝛥𝑙 = 6 m and 𝛥𝑡 = 3 s are considered to speed up the computational time. Thus, the sampling time of the controllers is 𝛥𝑡ctr =𝛾𝛥𝑡 = 60 s by setting 𝛾= 20. Initially, it is considered 𝑇in, 𝑇 f 1,𝑗 , 𝑇 1f 1,𝑗 = 155 ◦Cand 𝑇m 1,𝑗 = 165 ◦C,∀𝑗∈. 5.3. Results of the 10-loop field The simulation length is 𝑁sim = 120 min, and the two-hour effective solar DNI (𝜂col 𝑖,𝑗 𝐺𝑖,𝑗 𝐼𝑖,𝑗 (𝑘)) profile used as disturbances is shown in Fig. 3. As observed, it is considered that 20% of the ACUREX field has 480 W/m (loops #5 and #10), another 20% has 900 W/m (loops #3 and #4) and the rest has 800 W/m. Additionally, a moving cloud entering from 𝑡= 63 min to 𝑡= 76 min is considered (see Fig. 3 (b)–(d)). As observed, the effective solar DNI at 𝑡= 70 min decreases significantly. This moving Applied Energy 306 (2022) 117936 7 E. Masero et al. Fig. 3. Collectors view of the 10-loop field. Fig. 4. Total flow trajectories in the 10-loop collector field. Fig. 5. HTF flow trajectories of the loops specified in the legend for the 10-loop collector field. Applied Energy 306 (2022) 117936 8 E. Masero et al. Fig. 6. Formation of coalitions with the Coal. Market method in the ACUREX field. Fig. 7. Collectors view of the 100-loop field. Fig. 8. Total flow trajectories in the 100-loop collector field. Table 3 Numerical results obtained by applying different controllers in a two-hour simulation of the 10-loop solar field. Control method 𝜏[s]  𝐸[kWh] Improvement No-valves track 0.43 2093.6– No-valves 28.34 2308.6 10.27% Dec 2.5 2335.2 11.54% Coal. Pair 9.167 2347.8 12.14% Coal. Market 19.17 2349.4 12.22% Cen 275.28 2355.6 12.51% cloud represents small cumulus of clouds that can affect solar collector fields on partially cloudy days. Table 3 displays the mean energy obtained with all control methods in the ACUREX field. The resulting energy is delimited by the No-valves Track method, which provides the lowest energy  𝐸= 2093.6kWh, and the Cen method, which produces the highest energy  𝐸= 2355.6kWh. Conversely, the high computational time 𝜏= 275.28 s becomes centralized control in an unfeasible strategy to be implemented in real-time. Note that the sampling time of the controller is usually around 𝛥𝑡ctr = 60 s in these types of solar fields. Therefore, the No-valves, Dec, Coal. Pair, and Coal. Market schemes (with 10.27%, 11.54%, 12.14%, and 12.22% of mean energy improvement regarding No-valves Track, respectively) are more suitable to manage this solar collector field. Applied Energy 306 (2022) 117936 9 E. Masero et al. Fig. 9. Formation of coalitions with the Coal. Market method for the two-hour simulation. Fig. 4 shows the overall trajectories of the outlet temperature, the thermal power, and flow rate in the ACUREX field.2As observed in Fig. 4 (a), the outlet temperature of the field holds within limits 𝑇f,max = 300◦C and 𝑇f,max = 220◦C for all methods considered thanks to the controller and the defocus mechanism triggered as needed. Observe that the No-valves Track method holds the outlet temperature in the set-point 𝑇f ref = 280◦C. The total HTF flow also decreases when the moving cloud is covering part of the field (see Fig. 4 (c)). As an example, it is shown the outlet temperature and the local flow rate of a clean loop #2 and a dusty loop #10 in Fig. 5. The flow rate is decreased because of the drop in DNI between 𝑡= 63 min and 𝑡= 76 min. Providing the outlet temperature of a loop reaches 290◦C, the defocus mechanism is triggered, e.g., as seen in loop #2 with the No-valves method at 𝑡= 70 min. The coalitions formed by applying the Coal. Market scheme to the ACUREX field are displayed in Fig. 6. The 𝑦-axis represents the loops, the 𝑥-axis indicates the length of the simulation, and the color bar shows the ten different colors, i.e., it represents up to ten coalitions in total. For example, the partition of the system is = {{1},{2},{6},{7}, {4,8,10},{3,5,9}} at time instant 𝑡= 1 min, and = {{1},{4},{6},{7}, {8},{2,5},{3,9,10}} at 𝑡= 120 min. After evaluating the maximum coalition sizes as shown in Appendix, the limit was set to ||max = 3because it provides the best trade-off between performance and computational load. 5.4. Results of the 100-loop field For the 100-loop field, we also simulate two hours using the effective solar DNI profile represented in Fig. 7. As shown, it is considered that 20% of loops have 450 W/m, 16% have 900 W/m, and the rest have 750 W/m. Additionally, a moving cloud is considered passing over the field from 𝑡= 65 min to 76 min (see Fig. 7 (b)–(d)). Fig. 8 shows the overall trajectories of the outlet temperature, the thermal power, and the flow rate obtained in the 100-loop field by implementing the considered controllers. Similar to the previous case, the outlet temperature increases, and the flow rate decreases in the time interval that the moving cloud covers the field to maximize the resulting thermal power. Note that the centralized MPC has not been implemented because of the substantial computational (approximately 𝜏= 410 min). Therefore, the centralized approach is not suitable to handle this field in real-time. The resulting partitions with Coal. Market are depicted in Fig. 9, where the color bar represents up to 100 colors for the total number of coalitions. As an example, there are 60 coalitions at time instant 𝑡= 1 min, 100 coalitions at 𝑡= 66 min, and 43 coalitions at 𝑡= 101 min, with the maximum coalition size ||max = 3 established previously. 2The reader is referred to the electronic version of the manuscript to interpret the color legends of most figures. Table 4 Numerical results by applying the control techniques in the two-hour simulation of the 100-loop solar field. Control method 𝜏[s]  𝐸[kWh] Improvement No-valves track 2.66 20672 – No-valves 228.39 21276 2.92% Dec 3.11 22688 9.75% Coal. Pair 13.86 22778 10.19% Coal. Market 35.19 22826 10.42% The numerical results obtained by implementing the mentioned control strategies in the large-scale field are collected in Table 4. The mean energy obtained on applying the No-valves Track method is outperformed 2.92% by No-valves, 9.75% by Dec, 10.19% by Coal. Pair, and 10.42% by Coal. Market. In this regard, the proposed Coal. Market presents the best improvement of mean energy (10.42%) and can be computed in a reasonable time (𝜏= 35.19 s). Hence, the proposed coalitional MPC scale suitably for controlling this distributed solar collector field. 5.5. Discussion: Applicability to a real scenario The market-based scheme for selecting the solar field partition is a novel method inspired by the electricity intraday market. It allows loop controllers to offer and demand HTF flow aware of the utility gain of the transactions, which can be done almost instantaneously because we only increase/decrease the HTF limits of each loop within the coalition. Nowadays, there is no prototype of this mechanism in practice, but as we have shown it would be profitable and feasible to be implemented the method in a real scenario. Nevertheless, the applicability of the proposed method may present the following limitations in practice: (i) Flow meters are supposed to be installed in each loop. However, most commercial plants do not have flow meters in all loops, but only in various field sectors. In this case, the flow rate of each loop can be estimated from the total flow, which is often measured, and the valve position of all loops with the aid of a hydraulic model. (ii) Valves are already installed in many solar parabolic-trough plants and can be manipulated from the SCADA system. However, this manipulation is only made for balancing the flow of the solar field hydraulically. To the best of our knowledge, some solar power companies are now considering the use of loop valves for controlling their plants and not just for hydraulic balance. (iii) Scalability in large-scale fields. The proposed clustering approach is scalable from the 10-loop to the 100-loop solar collector field, as shown in the results section. In a larger field, the limitation of scalability is given by the computational time, which is composed of: (i) the time that the coalition requires to solve the MPC problem, which is related to the number of loops per coalition, and (ii) the time needed to compute the partition of the field. We have elaborated a simulation study with fields of increasing number of loops, as displayed in Fig. 10. Considering that real plants use a sampling time of 60 s, it would be feasible to control fields of up to 300 loops approximately with our current implementation. Implementations in larger fields would require more powerful computers and more efficient programming languages than Matlab, which focuses on mathematical prototyping rather than on the speed of the solution. In this regard, there are implementations of MPC in C and C++, e.g., QP-OASES [46] and C/GMRES [47], which could be employed. Therefore, our assessment must be considered as a lower bound of the size of the plants that could be calculated with the proposed method.