Clustering-based model predictive control of solar parabolic trough plants

Renewable Energy 216 (2023) 118978 Available online 10 July 2023 0960-1481/© 2023 The Authors. 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 Renewable Energy journal homepage: www.elsevier.com/locate/renene Clustering-based model predictive control of solar parabolic trough plants✩ Paula Chanfreuta,b,∗, José M. Maestrea, Antonio Gallegoa, Anuradha M. Annaswamyc, Eduardo F. Camachoa aDepartment of Systems and Automation Engineering, University of Seville, Spain bDepartment of Mechanical Engineering, Eindhoven University of Technology, The Netherlands cDepartment of Mechanical Engineering, Massachusetts Institute of Technology, USA ARTICLE INFO Keywords: Model predictive control Solar thermal power plants Parabolic trough collectors Control by clustering Coalitional control Hierarchical control ABSTRACT This paper presents a clustering-based model predictive controller for optimizing the heat transfer fluid (HTF) flow rates circulating through every loop in solar parabolic trough plants. In particular, we present a hierarchical approach consisting of two layers: a bottom layer, composed of a set of model predictive control (MPC) agents; and a top layer, which dynamically partitions the set of loops into clusters. Likewise, the top layer allocates a certain share of the total available HTF to each cluster, which is then distributed among the loops by the bottom layer in response to the varying conditions of the solar field, e.g., to deal with passing clouds. The dynamic clustering of the system reduces the number of variables to be coordinated in comparison with centralized MPC, thereby speeding up the computations. Moreover, the loops efficiencies and the heat losses coefficients, which influence the loops control model, are also estimated at the bottom layer. Numerical results on a 10-loop and an 80-loop plant are provided. 1. Introduction Renewable energy generation plays a central role in the path towards a more sustainable and resilient energy system. The benefits of renewable energies go beyond lowering greenhouse gas emissions and also include the diversification of the energy supply and the reduction of the dependency on fossil fuels. In 2021, the share of renewables in the global electricity generation reached about a 28% and is expected to rise until 38% by 2027 according to the International Energy Agency (IEA) [1]. In particular, wind and solar resources reached a 10% of the global electricity mix as reported by the Renewable Energy Policy Network for the 21st Century (REN21) [2]. The radiant energy per unit area that the Sun provides at the top of the Earth’s atmosphere is roughly 1361 W/m2[3]. In this regard, according to the World Energy Council (WEC) [4], the annual solar radiation reaching the Earth’s surface is about 7500 times the world’s annual energy consumption, which illustrates its potential. In recent years, the advances and investment on photovoltaic technology have brought a huge growth of solar energy generation [1], yet there are more 6GW of installed capacity of concentrated solar thermal (CST) power technology [2]. Unlike photovoltaics, which directly convert the sunlight into electricity, CST technologies use mirrors to focus the ✩This work is supported by the European Research Council Advanced Grant OCONTSOLAR under Grant SI-1838/24/2018, and by the Spanish MCIN/AEI/10.13039/501100011033 Project C3PO-R2D2 under Grant PID2020-119476RB-I00. ∗Correspondence to: Department Mechanical Engineering, Gemini Building number 15, Groene Loper, 5612 AE Eindhoven, The Netherlands. E-mail addresses: [email protected],[email protected] (P. Chanfreut), [email protected] (J.M. Maestre), [email protected] (A. Gallego), [email protected] (A.M. Annaswamy), [email protected] (E.F. Camacho). solar rays on a receiver to heat a working fluid, commonly referred to as heat transfer fluid (HTF) [5]. There are four main types of CST systems: parabolic troughs, linear Fresnel reflectors, power towers, and parabolic dishes [6,7]. While the first two use a linear receiver tube [8], the latter two concentrate the solar rays onto a point. The thermal energy gained by the HTF is then converted into mechanical energy using a turbine to produce electricity. Moreover, CST plants can incorporate thermal storage systems, which provides significant advantages to deal with intermittent and variable nature of the solar radiation [9–12]. This article is focused on solar parabolic trough plants [13–15], which represents the most widespread CST technology deployed up to date, with Spain and the United States holding the largest share of the installed capacity [2,16]. The solar field in this type of plants consists of a set of arrays of parabolic collectors, i.e., loops, through which circulates the HTF, which is typically a synthetic oil or a molten salt [17]. To control the HTF temperature around a given setpoint, its flow rate must be dynamically manipulated taking into account the plant dynamics and possible fluctuations of the received solar irradiance. Model predictive control (MPC) has provided positive results to address this problem as reviewed by Camacho et al. [18] (Chapter 5) https://doi.org/10.1016/j.renene.2023.118978 Received 23 February 2023; Received in revised form 12 June 2023; Accepted 29 June 2023 Renewable Energy 216 (2023) 118978 2 P. Chanfreut et al. and illustrated, for example, by Andrade et al. [19], Navas et al. [20], López-Bautista et al. [21] and Gholaminejad and Khaki-Sedigh [22]. MPC is an advanced control method based on iterative computations of the sequence of control inputs that optimizes the system performance along a certain finite time horizon [23,24]. Its inherent optimal nature, and ability to handle complex systems with constraints, are some of the advantages that have made it increasingly popular over the last decades [25,26]. More recently, distributed implementations have allowed scaling the benefits of MPC to large-scale systems, where a centralized control is not viable or practical for real time control [27]. The essence of distributed MPC (DMPC) is to tackle the overall control problem by using multiple local MPC agents that exchange information to coordinate their decisions. Unlike decentralized MPC, DMPC allows for cooperation among agents to improve performance [28], but this comes at the cost of increased communication and computational resources, which may not be desirable or efficient. Motivated by this, a number of recent works have explored distributed approaches with partial and dynamic inter-agent cooperation as surveyed by Chanfreut et al. [29]. The underlying idea to boost the scalability of distributed architectures is often to use time-varying communication topologies, e.g., [30,31]. One of the strategies proposed in this framework is the so-called coalitional control [32], an approach that is characterized by the dynamic formation of clusters of cooperative agents, referred to as coalitions. In this context, the project Optimal Control of Solar Energy Systems (OCONTSOLAR), which has been funded by the European Research Council [33], develops coalitional MPC algorithms to optimize the HTF flow rates in solar parabolic trough plants. Conversely to the classic operation of these plants, the flow in different loops are considered to be manipulable in real-time by using valves. The latter is of special interest when the conditions in the solar field are such that the capacity to warm up the HTF is uneven among the loops, e.g., due to cloud shading [34,35]. Indeed, by using the valves, it is possible to obtain a better thermal balance of the HTF temperature in each loop, which helps to minimize defocusing actions and hence energy losses [36]. However, optimizing the flows in every loop involves dealing with a large set of optimization variables that would also scale with the length of the prediction horizon. The use of a centralized MPC lacks of redundancy and imposes a substantial burden on computation and communication, which motivates the study of noncentralized MPC strategies and, in particular, the mentioned coalitional approaches. In addition, for the flows control to be effective, MPC-based schemes require an accurate model of the solar field. The loops dynamics are though characterized by parameters that are uncertain, difficult to measure, and which may also vary in time, such as their optical efficiencies and heat losses coefficients. The main contribution of this article is a hierarchical cluster-based MPC for optimizing the loops flow rates to track reference outlet temperatures. The proposed approach integrates a bottom layer, comprising a set of MPC agents that control the valves of each loop; and a top layer, which acts as a supervisor. The bottom layer uses approaches based on coalitional MPC algorithms and parameter identification using recursive least squares. On the other hand, the supervisor at the top layer partitions the set of loops in a dynamic manner into clusters, and assigns to these clusters a maximum HTF flow considering the plant constraints. To optimize these maximum flows, a clusters-based lumped description of the overall solar field is considered, which provides a simplified representation of the overall system. Our results on a 10-loop and a 80-loop plant show that the proposed approach can closely approximate the centralized performance while reducing both the computation times and the need for centralized computations. The rest of the article is organized as follows. Section 2presents the system dynamics and constraints, and introduces the control objectives. Section 3describes the architecture of the proposed controller, details the elements comprising the bottom and top layers, and provides a pseudo-code with the control algorithm. Finally, our simulation results and concluding remarks are respectively presented in Sections 4and 5. Fig. 1. Scheme of the solar field. The loops of parabolic collectors are equipped with valves that allow controlling the HTF flow rates that pass through them. 2. Problem formulation Consider a solar parabolic trough plant comprising a set of parallel loops = {1,2,…, 𝑁loops}as shown in Fig. 1. In this respect, throughout this paper we will consider that the HTF flow that is pumped to different loops can be controlled in real-time by using valves. 2.1. System dynamics and constraints Considering the variation of the internal energy of the fluid, the dynamics of the output temperature of a given loop 𝑖∈, i.e., 𝑇out 𝑖[◦C], can be described by the following concentrated parameter model1: 𝐶𝑖 𝑑𝑇 out 𝑖 𝑑𝑡 =𝜂𝑖ℐ𝑖−𝑞𝑖𝑃𝑖(𝑇out 𝑖−𝑇in) − h𝑖,(1) where 𝑇in [◦C] is the inlet temperature, 𝑞𝑖[m3/s] represents the HTF that is pumped to loop 𝑖,𝐶𝑖[J/◦C] is the thermal capacity of the loop, 𝑃𝑖 [J/(m3◦C)] is related to its geometrical and thermal properties, and h𝑖 [W] is a function weighting the heat losses of loop 𝑖. Also, ℐ𝑖=𝑆loop𝐼𝑖, where 𝑆loop [m2] denotes the loops’ reflective surface and 𝐼𝑖[W/m2] is the direct normal irradiance (DNI) received by loop 𝑖. In addition, 𝜂𝑖is the loop efficiency, which considers the optical and geometric efficiency of the collectors in 𝑖. Finally, some of the variables abovementioned are defined as a function of the temperature as follows: 𝜌𝑖= 903 − 0.672𝑇m 𝑖, 𝑐𝑖= 1820 + 3.478𝑇m 𝑖, 𝑃𝑖=𝜌𝑖𝑐𝑖, 𝐶𝑖=𝜌𝑖𝑐𝑖𝐴𝐿, h𝑖=𝑆loop (ℎ2,𝑖(𝑇m 𝑖−𝑇a)2+ℎ1,𝑖(𝑇m 𝑖−𝑇a)), 𝑇 m 𝑖= (𝑇out 𝑖+𝑇in)∕2, (2) where 𝑇a[◦C] is the ambient temperature, 𝐴[m2] is the cross sectional area of the tube, 𝐿[m] is the loops length, and coefficients ℎ1,𝑖 and ℎ2,𝑖 will be estimated as discussed in Section 3.1.2. Note that (2) considers the properties of HTF Therminol 55, which was used in the ACUREX solar plant of the Plataforma Solar de Almería [18, Chapter 4]. Regarding the system constraints, hereafter let 𝑄Tdenote the maximum available flow rate that can be used in the overall plant, 𝑞min the minimum flow rate allowed in the loops, and 𝑞max the corresponding maximum value. Then, it must be guaranteed that ∑ 𝑖∈ 𝑞𝑖≤𝑄T,and 𝑞min ≤𝑞𝑖≤𝑞max,∀𝑖∈.(3) Moreover, the outlet temperature is desired to evolve in range [𝑇min, 𝑇max]for all 𝑖∈, where 𝑇min and 𝑇max denote respectively the minimum and maximum temperatures. 1For the sake of clarity, the time index is omitted in this section. Renewable Energy 216 (2023) 118978 3 P. Chanfreut et al. 2.1.1. Cluster-based model As will be discussed in Section 3.2, the set of loops will be dynamically partitioned into a set of 𝑁cl non-overlapping clusters, i.e., {1,2,…,𝑁cl }, such that 𝑁cl ⋃ 𝑗=1 𝑗=,and 𝑗∩𝑙= ∅ ∀𝑗, 𝑙 ∈ [1, 𝑁cl], 𝑗 ≠𝑙. (4) Similarly to (1), the dynamics of any cluster 𝑗can be modeled as: 𝐶𝑗 𝑑𝑇 out 𝑗 𝑑𝑡 =𝜂𝑗ℐ𝑗−𝑞𝑗𝑃𝑗(𝑇out 𝑗−𝑇in) − h𝑗,(5) where 𝑇out 𝑗is the outlet temperature of cluster 𝑗,𝑞𝑗denotes the total HTF pumped to the loops in the cluster, and 𝐶𝑗,𝑃𝑗and h𝑗are defined analogously to (2). In addition, 𝜂𝑗ℐ𝑗=∑𝑖∈𝑗𝜂𝑖ℐ𝑖represents the effective irradiance in cluster 𝑗. 2.2. Control objectives In this paper, each loop 𝑖∈will be treated as a subsystem controlled by an MPC agent that can manipulate its inlet valve, thus controlling flow 𝑞𝑖. In turn, the set of MPC agents can communicate to coordinate their decisions to some extent. In this framework, the main control goals and premises considered can be summarized as follows: •The overall system should track piecewise constant references on the loops outlet temperature. In this respect, the total available HTF can be unevenly distributed among the set of loops to exploit the time- and space-varying conditions of the solar field. •Centralized and fully coordinated control of all valves should be avoided due to the high computation and communication burden. In this regard, the system performance should approximate that of the centralized controller but reducing the need for centralized computations. •Parameter 𝜂𝑖and the heat losses coefficients, i.e., ℎ1,𝑖 and ℎ2,𝑖, must be estimated for all loops 𝑖∈. The efficiency 𝜂𝑖is assumed to be slowly time-varying, and ℎ1,𝑖 and ℎ2,𝑖 can be considered constants. Regarding the outlet temperature reference, hereafter 𝑇ref , our controller will be designed to be integrated with a setpoint optimization layer as described in Camacho and Gallego [37]. That is, we assume that a different control layer will dynamically provide us with the optimal 𝑇ref to maximize the net electricity production or satisfy a given scheduled electrical production. Remark 1. Parameter 𝜂𝑖considers the geometric and optical efficiency of the collectors in loop 𝑖. The geometric efficiency depends on the position of the collectors with respect to the radiation beam vector, whereas the optical efficiency takes into account aspects as the mirrors reflectivity, the interception factor, and the absorptance of the tube [38, 39]. This efficiency can vary over time due to several reasons, e.g., dust accumulation on the mirrors or deterioration of the selective absorbing coating. Conversely to the heat losses coefficients, which change due to the tubes degradation and vacuum losses on a yearly time scale, changes in the loops efficiencies can be observed over the days. 3. Proposed controller architecture Considering the objectives in Section 2.2, this article proposes the hierarchical approach that is illustrated in Fig. 2. As can be seen, it consists of a top and a bottom layer comprising the following elements: (i) A set of local agents associated to the loops, which optimize the loops flow 𝑞𝑖, and identify the unknown parameters 𝜂𝑖,ℎ1,𝑖, and ℎ2,𝑖, for all 𝑖∈. Fig. 2. Hierarchical architecture of the proposed approach. At the bottom layer, a set of agents control the HTF flow rates in each loop, whereas the top layer adjusts their level of coordination and constraints considering (3). (ii) A supervisor, which partitions dynamically the set of loops into clusters, and updates the flow constraints using the corresponding cluster-based description of the field. In particular, given a partition of the set of loops, say {1,2,…,𝑁cl }, the agents associated with each cluster 𝑗will share data to coordinate their decisions on flow rates 𝑞𝑖for all 𝑖∈𝑗. As introduced before, these groups of cooperative agents will be also denote as coalitions. That is, coordination is limited to subsets of agents, thus providing a half-way solution between a centralized and complete decentralized controllers. The optimization within each coalition will consider the constraints imposed by the supervisor, which will set a maximum on the HTF that they can use. Additionally, the proposed approach considers different time scales. For simulation purposes, we use a discrete-time version of model (1) for all 𝑖∈with an integration step size of 𝛥𝑡sseconds. In this regard, in what follows we use 𝑘as the discrete time index, i.e., time step 𝑘refers to instant 𝑘𝛥𝑡s. On the other hand, the control step size to update the flow rates is 𝛥𝑡c=𝛿c𝛥𝑡sseconds, where 𝛿c∈N+. For the sake of simplicity, we will consider that the sampling time of the model used by the agents and supervisor is also 𝛥𝑡c. The following subsections provide further details on the function and design of the different elements involved in the presented approach. 3.1. Bottom layer The bottom layer determines dynamically the HTF flow rates that circulate through the loops and identifies the unknown parameters. 3.1.1. HTF flow optimization Considering (1), the dynamics of each loop 𝑖∈will be predicted using the following discrete-time model: 𝑇out 𝑖(𝑛+𝛿c|𝑘) = 𝑇out 𝑖(𝑛|𝑘) + 𝛾𝑖(𝑛|𝑘)(𝜂𝑖(𝑘)ℐ𝑖(𝑘) −𝑞𝑖(𝑛|𝑘)𝑃𝑖(𝑛|𝑘)(𝑇out 𝑖(𝑛|𝑘) − 𝑇in(𝑘))− h𝑖(𝑛|𝑘)),(6) where the notation (𝑛|𝑘)indicates the prediction for time step 𝑛≥ 𝑘made at time 𝑘, and 𝛾𝑖(𝑛|𝑘) = 𝛥𝑡c∕𝐶𝑖(𝑛|𝑘). In addition, 𝑃𝑖(𝑛|𝑘) and 𝐶𝑖(𝑛|𝑘)are computed according to (2) considering the predicted mean temperature, i.e., 𝑇m 𝑖(𝑛|𝑘) = (𝑇out 𝑖(𝑛|𝑘) + 𝑇in(𝑘))∕2, and the current ambient temperature 𝑇a(𝑘). Finally, 𝜂𝑖(𝑘)is the estimated efficiency of loop 𝑖at 𝑘, and  h𝑖(𝑛|𝑘) = 𝑆loop( ℎ2,𝑖(𝑘)(𝑇m 𝑖(𝑛|𝑘) − 𝑇a(𝑘))2+  ℎ1,𝑖(𝑘)(𝑇m 𝑖(𝑛|𝑘) − 𝑇a(𝑘))), with  ℎ1,𝑖(𝑘)and  ℎ2,𝑖(𝑘)being the coefficients estimation available at 𝑘. Consider model (6) and take some time step 𝑘such that 𝑘=𝜅𝛥𝑡c, where 𝜅∈N+, i.e., the flow rates are updated at step 𝑘. Also, let (𝑘) be the agents partition selected at 𝑘, and 𝑗be the 𝑗th cluster in (𝑘). Renewable Energy 216 (2023) 118978 4 P. Chanfreut et al. Finally, let us define the prediction horizon used by the MPC as an interval of 𝛿c𝑁pintegration steps ahead, and consider set = {𝑘, 𝑘 + 𝛿c, 𝑘+2𝛿c,…, 𝑘+𝛿c𝑁p}. Then, at time 𝑘, coalition 𝑗solves the following MPC problem: min [𝑞𝑖]𝑖∈𝑗∑ 𝑖∈𝑗∑ 𝑛∈ (𝑤e𝑒2 𝑖(𝑛+𝛿c|𝑘) + 𝑤q𝛥𝑞2 𝑖(𝑛|𝑘)) s.t. Loops prediction model (6), 𝑇out 𝑖(𝑘|𝑘) = 𝑇out 𝑖(𝑘),(7a) 𝑇min ≤𝑇out 𝑖(𝑛|𝑘)≤𝑇max,(7b) 𝑞min ≤𝑞𝑖(𝑛|𝑘)≤𝑞max,(7c) ∑ 𝑖∈𝑗 𝑞𝑖(𝑛|𝑘)≤𝛼𝑗(𝑛)𝑄T,(7d) ∀𝑖∈𝑗,∀𝑛∈,(7e) where 𝑤eand 𝑤qare positive definite weighting scalars, 𝛼𝑗(⋅)is a scale factor that will be set by the top layer, and 𝑒𝑖(𝑛+𝛿c|𝑘) = 𝑇out 𝑖(𝑛+𝛿c|𝑘) − 𝑇ref (𝑛+𝛿c), 𝛥𝑞𝑖(𝑛+𝛿c|𝑘) = 𝑞𝑖(𝑛+𝛿c|𝑘) − 𝑞𝑖(𝑛|𝑘),∀𝑛∈,(8) Notice that at the initial time instant 𝛥𝑞𝑖(𝑘|𝑘) = 𝑞𝑖(𝑘|𝑘) − 𝑞𝑖(𝑘−𝛿c), where 𝑞𝑖(𝑘−𝛿c)is known and equal to 𝑞𝑖(𝑘− 1) since the control inputs are updated each 𝛿csteps. Remark 2. In problem (7), prediction model (6) considers as simplification that temperatures 𝑇in and 𝑇amaintain their value at 𝑘during the prediction horizon. Similarly, the effective irradiance is defined as the product of the last updated estimated efficiency and the DNI at 𝑘, which is assumed to be known. 3.1.2. Parameters identification As introduced in Remark 1, the loops efficiencies and heat losses coefficients are unknown parameters that can vary in time for multiple reasons, e.g. dust accumulation and vacuum losses. To estimate them, this article considers a local system identifier based on the recursive least-square method with exponential forgetting [40] that is designed as follows. Let us consider continuous-time model (1), which can be rewritten as 𝑑𝑇 out 𝑖 𝑑𝑡 = 𝑇out 𝑖=𝜑𝑖+𝜙ℐ 𝑖𝜂𝑖+𝜙h 𝑖[ℎ2,𝑖, ℎ1,𝑖]T,(9) where 𝜑𝑖= −𝑞𝑖𝑃𝑖(𝑇out 𝑖−𝑇in)∕𝐶𝑖,𝜙ℐ 𝑖=ℐ𝑖∕𝐶𝑖, and 𝜙h 𝑖= −𝑆loop∕𝐶𝑖[(𝑇m 𝑖− 𝑇a)2,(𝑇m 𝑖−𝑇a)]. Also, note that parameters 𝜂𝑖,ℎ1,𝑖 and ℎ2,𝑖 may not be uniquely identifiable and, moreover, they change in different time scales (see Remark 1). To deal with this issue, we consider that ℎ1,𝑖 and ℎ2,𝑖 will only be estimated by using data collected in periods of low irradiance, e.g., during defocusing actions or at the end of the day. In that case, 𝜙ℐ 𝑖≈ 0, and hence 𝜂𝑖does not affect the temperature evolution. Once we have an estimation of ℎ1,𝑖 and ℎ2,𝑖, they will be treated as known constant parameters, and we will only adapt 𝜂𝑖during the periods of higher irradiance. For the estimation of ℎ1,𝑖 and ℎ2,𝑖, consider a filtered version of the variables in (9), i.e.,  𝑇out 𝑖,f=𝜆(𝑇out 𝑖−𝑇out 𝑖,f), 𝜙h 𝑖,f=𝜆(𝜙h 𝑖−𝜙h 𝑖,f)and 𝜑𝑖,f=𝜆(𝜑𝑖−𝜑𝑖,f), where 𝜆≥0is the filter inverse constant, and sub-index ‘f’ holds for filtered [41]. Then, we can transform (9) into: 𝑧𝑖=𝜆(𝑇out 𝑖−𝑇out 𝑖,f) − 𝜑𝑖,f=𝜙h 𝑖,f[ℎ2,𝑖, ℎ1,𝑖]T.(10) Considering this, we will use the following least squares adaptive law with exponential forgetting for adjusting these coefficients: [  ℎ2,𝑖   ℎ1,𝑖]=𝛤𝑖𝜙h 𝑖,f(𝑧𝑖−𝜙h 𝑖,f[ ℎ2,𝑖  ℎ1,𝑖]),[ ℎ2,𝑖(0)  ℎ1,𝑖(0)]= ℎ0,(11) where  ℎ0is the initial estimation, and  𝛤𝑖=𝛽𝛤𝑖−𝛤𝑖𝜙h 𝑖,f(𝜙h 𝑖,f)T𝛤𝑖, 𝛤𝑖(0) = 𝛤0,(12) with 𝛽≥0being the continuous time forgetting factor, and 𝛤0>0. Notice that Eqs. (11) and (12) can be discretized providing a discretetime adaptation law that does not require a discretized measure of derivative  𝑇out 𝑖but only samples of 𝑇out 𝑖, i.e., it depends on variables that can be directly measured in the system. In this regard, the leastsquares method with exponential forgetting has been chosen for the simplicity of its implementation and its convergence guarantees of the parameters error under persistently exciting regression vectors. Finally, for the adaptation of 𝜂𝑖along the days, we will proceed analogously and fix the heat losses coefficients to their last estimated value. Note that then (9) can be rewritten in the form of  𝑇out 𝑖=𝜑𝑖+ 𝜙ℐ 𝑖𝜂𝑖, where 𝜑𝑖is treated as a known variable. 3.2. Top layer The top layer sets periodically a partition into coalitions and distributes the total available flow among them for its optimization at the bottom layer. 3.2.1. Clusters formation The set of agents will be dynamically clustered into coalitions such that the dynamics of the loops controlled by the same coalition are similar. Therefore, the elements in the same cluster are approximately characterized by the same lumped parameters. To this end, we build up the following matrix at all control time steps 𝑘: 𝐏(𝑘) = [𝜂𝑖(𝑘)ℐ𝑖(𝑘), 𝑇 m 𝑖(𝑘), h𝑖(𝑘)]𝑖∈.(13) Note that, given (2), variable 𝑇m 𝑖(𝑘)characterizes parameters 𝐶𝑖(𝑘)and 𝑃𝑖(𝑘)that also appear in the loops model. Using matrix 𝐏(𝑘), we can find a partition of the set of loops in the solar field using clustering methods [42]. Without loss of generality, we consider the well-known 𝐾-means algorithm [43], which is an iterative centroid-based algorithm that provides a partition of a data set into 𝐾 clusters. In particular, we will use Matlab®function kmeans to find a partition of the 𝑁loops entries in 𝐏(𝑘)into 𝑁cl clusters.2 The number of clusters 𝑁cl can be selected following the elbow method, which runs the 𝐾-means clustering algorithm for a range of values for 𝐾, and for each of them computes an average score for the resulting clusters. The range of 𝐾to be tested is defined as [2, 𝑁max cl ], where 𝑁max cl denotes the maximum number of clusters. As for the score, different criteria can be used [44], being the Calinski-Harabasz index [45] the one considered in our simulations, which measures the between-cluster and within-cluster variance so as to find dense and well separated clusters. Remark 3. Setting 𝑁max cl close to 𝑁loops generates smaller coalitions, and hence sparser communication topologies, but will increase the number of decision variables of the supervisor (presented in Section 3.2.2). For this reason, 𝑁max cl should be set so that the number of optimization variables handled by the coalitions and those of the supervisor are significantly fewer than those corresponding to the centralized control problem. Also, smaller 𝑁max cl allows the real-time implementation of the elbow method to find the optimal partition. 2Note that 𝑁cl corresponds to the value of 𝐾used when running the 𝐾-means clustering algorithm. Renewable Energy 216 (2023) 118978 5 P. Chanfreut et al. 3.2.2. Clusters-based flow allocation The supervisor will consider a clusters-based model of the overall solar field to find a candidate distribution of the HTF among the coalitions, which will then be used to set the constraints for the bottom layer. The model used at the top layer describes the solar field as an aggregation of clusters of loops modeled by a discrete-time version of (5). Let 𝑘be the current time instant, (𝑘)the chosen partition, and consider some cluster 𝑗∈(𝑘). Then, the prediction model of 𝑗is given by: 𝑇out 𝑗(𝑛+𝛿c|𝑘) = 𝑇out 𝑗(𝑛|𝑘) + 𝛾𝑗(𝑛|𝑘)(𝜂𝑗(𝑘)ℐ𝑗(𝑘) −𝑞𝑗(𝑛|𝑘)𝑃𝑗(𝑛|𝑘)(𝑇out 𝑗(𝑛|𝑘) − 𝑇in(𝑘))− h𝑗(𝑛|𝑘)), (14) where 𝛾𝑗(𝑛|𝑘) = 𝛥𝑡c∕𝐶𝑗(𝑛|𝑘). An overall model of the field can therefore be defined as the aggregation of (14) for all 𝑗∈(𝑘). Considering the latter, the MPC problem solved by the supervisor is the following: min 𝐪(⋅|𝑘)∑ 𝑗∈(𝑘)∑ 𝑛∈top (𝑤top e,𝑗‖𝑒𝑗(𝑛+1|𝑘)‖2+𝑤top q,𝑗‖𝛥𝑞𝑗(𝑛|𝑘)‖2) s.t. Clusters prediction model (14), 𝑇out 𝑗(𝑘|𝑘) = ∑𝑖∈𝑗𝑞𝑖(𝑘− 1)𝑇out 𝑖(𝑘) ∑𝑖∈𝑗𝑞𝑖(𝑘− 1) ,(15a) 𝑇min ≤𝑇out 𝑗(𝑛|𝑘)≤𝑇max,(15b) |𝑗|𝑞min ≤𝑞𝑗(𝑛|𝑘)≤|𝑗|𝑞max,(15c) ∑ 𝑗∈(𝑘) 𝑞𝑗(𝑛|𝑘)≤𝑄T,(15d) ∀𝑗∈(𝑘),∀𝑛∈top,(15e) where 𝑤top e,𝑗and 𝑤top q,𝑗are positive definite weighting scalars, 𝐪(𝑛|𝑘) = [𝑞𝑗(𝑛|𝑘)]𝑗∈(𝑘)is a vector aggregating the predicted flow rates for all clusters at instant 𝑛, and top = {𝑘, 𝑘 +𝛿c, 𝑘 + 2𝛿c,…, 𝑘 +𝛿c𝑁top p}, i.e., the prediction horizon is set as an interval of 𝛿c𝑁top pintegration steps ahead. Also, |𝑗|is the cardinality of set 𝑗, i.e., the number of loops in the cluster, and 𝑒𝑗(𝑛+𝛿c|𝑘) = 𝑇out 𝑗(𝑛+𝛿c|𝑘) − 𝑇ref (𝑛+𝛿c), 𝛥𝑞𝑗(𝑛+𝛿c|𝑘) = 𝑞𝑗(𝑛+𝛿c|𝑘) − 𝑞𝑗(𝑛|𝑘),∀𝑛∈top.(16) Note that constraint (15d) implies that the sum of all clusters’ allocated flow rates do not exceed the total flow rate that can be used in the plant (see (3)). The minimizer obtained by solving (15) is therefore a sequence of HTF flow rates to be pumped to different clusters of loops during the following time steps, say [𝑞∗ 𝑗(𝑘|𝑘), 𝑞∗ 𝑗(𝑘+𝛿c|𝑘),…, 𝑞∗ 𝑗(𝑘+𝛿c𝑁top p|𝑘)], for all 𝑗∈(𝑘). Considering this, let us define 𝛼𝑗(𝑛) = 𝑄T 𝑞∗ 𝑗(𝑛|𝑘) ∑𝑗∈(𝑘)𝑞∗ 𝑗(𝑛|𝑘), 𝑛 ∈top,(17) as a factor that will be used to scale the flow constraints of the MPC agents at the bottom layer. Notice that, given (17),𝛼𝑗(𝑛)>0and ∑𝑗𝛼𝑗(𝑛)=1, for all 𝑛∈top. Remark 4. If the maximum flow rate 𝑄Tis evenly distributed among all clusters, then each 𝑗would have to satisfy 𝑞𝑗(⋅)≤𝑄T|𝑗|∕𝑁loops so that the overall system meets (3). Instead of using ratio |𝑗|∕𝑁loops, these constraints can be made more or less restrictive by introducing time-varying factors 𝛼𝑗(⋅), i.e., 𝑞𝑗(⋅)≤𝛼𝑗(⋅)𝑄T. The latter allows us to readjust the amount of HTF available in different clusters, e.g., to provide greater amount of flow to those clusters receiving higher DNI. Finally, to limit the extent to which the flow constraints are reduced, the scale factors resulting from (17) can be corrected so that they become greater than a certain threshold. For example, consider that at least a flow rate of 𝜁𝑄T∕𝑁loops must be available for every loop at all time steps, where 𝜁 > 0. Then, instead of using 𝛼𝑗(⋅), we consider in (7) scale factors 𝛼𝑗(⋅) = 𝛼𝑗(⋅) + 𝜀𝑗(⋅), where the value of 𝜀𝑗(⋅)is computed by solving the following problem: min [𝜀𝑗(⋅)]𝑗∈(𝑘)∑ 𝑗∈(𝑘)∑ 𝑛∈top 𝜀𝑗(𝑛)2 s.t. 𝛼𝑗(𝑛) = 𝛼𝑗(𝑛) + 𝜀𝑗(𝑛)≥𝜁|𝑗|∕𝑁loops, ∑ 𝑗∈(𝑘) 𝛼𝑗(𝑛)=1. (18) 3.3. Control algorithm The proposed hierarchical MPC algorithm is summarized in Algorithm 1. The overall solar field is simulated in step 11 of this algorithm as an aggregation of parallel loops modeled as 𝑇out 𝑖(𝑘+ 1) = 𝑇out 𝑖(𝑘) + 𝛾(𝑘)(𝜂𝑖(𝑘)ℐ𝑖(𝑘) −𝑞𝑖(𝑘)𝑃𝑖(𝑘)(𝑇out 𝑖(𝑘) − 𝑇in(𝑘))−h𝑖(𝑘)),(19) where 𝛾(𝑘) = 𝛥𝑡s∕𝐶𝑖(𝑘), and 𝛥𝑡sis the simulation step size (see Section 3). The outlet temperature of the overall field at instant 𝑘is accordingly defined as 𝑇out (𝑘+ 1) = ∑𝑖∈𝑞𝑖(𝑘)𝑇out 𝑖(𝑘+ 1) ∑𝑖∈𝑞𝑖(𝑘).(20) In addition, the inlet temperature is governed by the following first order system: 𝑇in(𝑠) 𝑇out (𝑠) − 80 =1 600𝑠+ 1 (21) where 𝑇in is expressed as a function of the outlet temperature of the steam generator, approximated here as 𝑇out − 80 ◦C. Finally, regarding the flows optimization in step 9 of Algorithm 1, it is sufficient for one of the agents in each coalition to solve (7) and communicate the solution to the rest of members. In this case, the agents in each coalition would internally adopt a star communication topology to share the necessary data with the one performing the optimization. Alternatively, we may consider distributed optimization methods, such as ADMM Algorithm 1 Control algorithm Let c= {0, 𝛿c,2𝛿c, ...}be the set of instants in which the flow rates are optimized and updated. Also, assume that the partition into clusters is revised each 𝛥𝑡cl =𝛿cl𝛥𝑡𝑠seconds, and, similarly, define cl = {0, 𝛿cl,2𝛿cl, ...}. Then, at all time instants 𝑘, proceed as follows: 1: if 𝑘∈cthen 2: Compute matrix 𝐏(𝑘)given in (13). 3: if 𝑘∈cl then 4: Update the clusters as described in Section 3.2.1 and define partition (𝑘). 5: end if 6: if ‖𝐏(𝑘) − 𝐏(𝑘− 1)‖∞> 𝜖 or (𝑘)≠ (𝑘− 1) then 7: The supervisor solves (15), and updates scale factors 𝛼𝑗for all 𝑗∈(𝑘)by using (17) and (18). 8: end if 9: All coalitions 𝑗∈(𝑘)solve problem (7) and find the flow rates to be implemented during interval [𝑘, 𝑘 +𝛿c). 10: end if 11: Update the system state and parameters considering (19),(2), and the current flow rates for each loop 𝑖∈. Renewable Energy 216 (2023) 118978 6 P. Chanfreut et al. or ALADIN [46], where all the agents in each coalition will perform iterative computations until converging to a solution. While we will consider the former approach in our simulations, both are equally valid and lead to similar results. 4. Simulation results In this section, we apply the presented coalitional control algorithm to a model of a solar parabolic trough plant, and provide a comparison with decentralized MPC, centralized MPC, and with the case in which no loops valves are used. Moreover, we consider the following scenarios: (i) Parameters 𝜂𝑖,ℎ1,𝑖, and ℎ2,𝑖 are identified as described in Section 3.1.2 for all loops 𝑖∈. (ii) Efficiency 𝜂𝑖is identified following Section 3.1.2, but the heat losses coefficients are set to nominal values ℎ1= −0.06133 and ℎ2= 0.00249 for all loops 𝑖∈.3 (iii) Both the efficiency and heat losses coefficients are set to constant nominal values 𝜂= 0.7,ℎ1and ℎ2for all loops 𝑖∈. For a better illustration of the algorithm performance, we first present results using a 10-loop plant, and subsequently on a extended plant with 80 loops. All simulations have been carried out in a 1.8 GHz quad-core Intel®CoreTM i7/8 GB RAM computer using Matlab®and software CasADi [47] for solving non-linear problems (7) and (15).4 Besides the properties introduced in (2) (i.e., 𝜌𝑖and 𝑐𝑖for all 𝑖∈ {1,2,…, 𝑁loops}), other parameters used in the simulations are given in Table 1. Note that the loops have been modeled as those of the ACUREX solar plant of the Plataforma Solar de Almería (see Chapter 4 of Camacho et al. [18]). Likewise, we considered 𝑄T= 9 ⋅10−4𝑁loops [m3/s], where 𝑁loops is equal to 10 or 80 as appropriate. Finally, 𝑤top e,𝑗and 𝑤top q,𝑗are respectively defined as |𝑗|𝑤eand |𝑗|𝑤q, where |𝑗|is the size of the corresponding cluster. Notice that in the centralized MPC case, a single controller optimizes all loops flow rates considering prediction model (6) for all 𝑖= 1,…, 𝑁loops. Likewise, the flow rate constraints in the centralized problem (recall (3)) are given by 𝑞min ≤𝑞𝑖(⋅|𝑘)≤𝑞max and ∑𝑁loops 𝑖=1 𝑞𝑖(⋅|𝑘)≤ 𝑄T. On the other hand, in the decentralized case, each agent 𝑖considers only its loops dynamics and imposes as flow constraints 𝑞min ≤𝑞𝑖(⋅|𝑘)≤ 𝑄T∕𝑁loops to satisfy (3). Therefore, the centralized performance would match that of a grand-coalition (𝑘) = = {1,2,…, 𝑁loops}during the entire simulation, with 𝛼(𝑘)=1, for all 𝑘≥0. Conversely, the decentralized case corresponds to the case of imposing partition (𝑘) = {{1},{2},…,{𝑁loops}} and 𝛼{𝑖}= 1∕𝑁loops for all 𝑘≥0and 𝑖= 1,…, 𝑁loops. Finally, in the no-valves case, the HTF flow rate is evenly distributed among all loops at all time steps. In particular, the control model considers the loops dynamics, but we impose that the HTF flow rate is equal for all loops. 4.1. 10-loop plant The simulations of the 10-loop plant consider a time period of 7 h (from 8:30 am to 3:30 pm) of a day in which the DNI received by each loop and the ambient temperature evolve as shown in Fig. 3. As can be seen, there are temporary drops in the irradiance evolution due to the effect of passing clouds. Fig. 4 shows the evolution of the loops outlet temperature and the HTF flows in the coalitional, decentralized and no-valves cases considering scenario (i) described at the beginning of this section. 3Coefficients ℎ1and ℎ2are those modeling the heat losses of the ACUREX plant. In our simulations, the actual ℎ1,𝑖 and ℎ2,𝑖 belong in range [0.85ℎ1,1.15ℎ1] and [0.85ℎ2,1.15ℎ2], respectively, for all 𝑖∈. 4The constraints on the maximum outlet temperature have been imposed as soft constraints in the presented simulations. Table 1 Parameters used in the simulation. Value Unit Value Unit Value Unit 𝑞min 0.2 ⋅10−3m3/s 𝑆loop 267.4 m2𝛥𝑡s0.5 s 𝑞max 2⋅10−3m3/s 𝑤e1 – 𝛥𝑡c30 s 𝑇min 220 ◦C𝑤q5⋅108–𝛥𝑡cl 150 s 𝑇max 305 ◦C𝑁p5 – 𝜁0.1 – 𝐴5.067 ⋅10−4m2𝑁top p10 –  𝛤01⋅10−4– 𝐿142 m 𝜖60 – 𝜆0.02 – Fig. 3. Evolution of the ambient temperature (top) and of the DNI for each loop (middle). The bottom plot zooms the irradiance graph in a time period affected by clouds. As can be seen, in the coalitional and decentralized cases, the loops outlet temperature follows closely the reference, and the flow rates decrease as the irradiance falls. However, if no valves are used and the temperature of some loop is close to 𝑇max, then the controller tends to set all loops flows close to the maximum possible, thus decreasing their temperature even for those that were already below the reference. The centralized graphs have not been added to Fig. 4 because they were very similar to the coalitional ones. Indeed, the maximum difference between the loops outlet temperature in both cases was 3 ◦C. This is also reflected in Table 2, which provides the cumulative performance costs in each of the simulations (see rows related to scenario (i)). In particular, it provides the value of the following indexes: 𝐽l=∑ 𝑘 𝑁loops ∑ 𝑖=1 (𝑤e𝑒2 𝑖(𝑘) + 𝑤q(𝑞𝑖(𝑘) − 𝑞𝑖(𝑘−1))2),and 𝐽g=∑ 𝑘(𝑤e‖𝑒(𝑘)‖2+ 𝑁loops ∑ 𝑖=1 𝑤q(𝑞𝑖(𝑘) − 𝑞𝑖(𝑘−1))2), (22) where 𝑒𝑖(𝑘) = 𝑇out 𝑖(𝑘) −𝑇ref (𝑘)and 𝑒(𝑘) = 𝑇out (𝑘) − 𝑇ref (𝑘). Note that 𝐽l weights the local error in the outlet temperature of each loop, whereas 𝐽gconsiders the error in the global outlet temperature of the solar field (see (20)). In view of Fig. 4 and Table 2, it should also be noted that the centralized and coalitional controllers introduce further flexibility in the distribution of the HTF among the loops, which allows for a better performance considering the same 𝑄T. While in the decentralized case all loops flows are kept lower or equal 𝑄T∕10, this threshold is occasionally exceeded in the centralized and coalitional cases without compromising the satisfaction of constraints (3). In this respect, Fig. 5 provides the evolution of the total HTF pumped to the solar field, i.e., ∑10 𝑖=1 𝑞𝑖(𝑘), during part of the simulations of scenario (i). Renewable Energy 216 (2023) 118978 7 P. Chanfreut et al. Fig. 4. Evolution of the loops outlet temperature and of the HTF flow rates in the 10-loop plant with different controllers in scenario (i). The dashed black line indicates the reference temperature. In addition, Figs. 6 and 7focus on the tasks performed by the supervisor in the proposed coalitional approach, i.e., the loops clustering and the constraints update. In particular, Fig. 6 illustrates the moving averages over a sliding window of 15 min of the maximum, mean, and minimum number of loops per clusters over time. On average, there were 3.51 clusters in the system with 2.97 loops. These clusters were selected using the 𝐾-means algorithm and considering 𝑁max cl = 4. Note that that the number and members of these coalitions are dynamically updated leading to time-varying cooperation structures at the bottom layer. In addition, Fig. 7 provides a swarm scatter chart with the values of factors 𝛼set by the supervisor for different clusters sizes. Recall that this factor determines flow constraints (7d), and hence the share of 𝑄Tthat the coalitions can use. Also, note that this type of plot illustrates the distribution of data while also showing the individual points. In particular, for a given coalition size ||, the values obtained for 𝛼throughout the entire simulation are jittered in the x-dimension, i.e., moved aside a bit so that any of them overlap. As can be seen, many of the points concentrate around the value of ||∕10 for each ||, which correspond to an even distribution of 𝑄Tamong the 10 loops. Nonetheless, we also obtained factors 𝛼differing from the latter, e.g., for ||= 3, there are some points around 0.2. This indicates that at some time instants the flow constraints of some loops were tightened in favor of others. Finally, Fig. 8 illustrates the estimation of the loops efficiency and thermal losses coefficients also in scenario (i). These results have been obtained by discretizing adaptation law described in Section 3.1.2 with a sampling time equal to 𝛥𝑡s, and using as initial estimations 𝜂0=𝜂and  ℎ0= [ℎ2, ℎ1]T. During part of the simulation, the efficiency of three of the loops was decreased to simulate a deficient sun tracking. Note that 𝜂𝑖, for all 𝑖= 1,…,10, is adapted in real time in the simulations of the centralized, coalitional, decentralized and no-valves controllers. On the other hand, coefficients  ℎ1,𝑖 and  ℎ2,𝑖 are assumed to be computed during a low irradiance period on some of the previous days. In particular, Fig. 8 (top) shows the evolution of 𝜂𝑖in the coalitional case, and the computation of  ℎ1,𝑖 and  ℎ2,𝑖 is illustrated in Fig. 8 (bottom). That is, the values of  ℎ1,𝑖 and  ℎ2,𝑖 at convergence were then used in the loops control model of all presented simulations of scenario (i). Regarding Fig. 5. Evolution of the total HTF flow rate with different controllers during a cloudy period. Renewable Energy 216 (2023) 118978 8 P. Chanfreut et al. Table 2 Cumulative performance costs in the 10-loop plant simulations. Centr. Coalitional Decent. No-valves Centr. Coalitional Decent. No-valves 𝐽l (i) 2.760 ⋅1052.857 ⋅1055.288 ⋅1051.030 ⋅107 𝐽g (i) 2.649 ⋅1042.681 ⋅1043.701 ⋅1041.485 ⋅105 (ii) 2.812 ⋅1053.600 ⋅1055.344 ⋅1051.030 ⋅107(ii) 2.668 ⋅1042.925 ⋅1043.723 ⋅1041.487 ⋅105 (iii) 1.095 ⋅1061.063 ⋅1061.190 ⋅1061.079 ⋅107(iii) 7.661 ⋅1047.355 ⋅1046.804 ⋅1041.967 ⋅105 Fig. 6. Moving averages over a sliding window of 15 min of the maximum, mean and minimum number of loops per cluster in scenario (i). Fig. 7. Swarm chart of the scale factors 𝛼resulting for different coalitions sizes ||. Note that if 𝑄Tis evenly distributed among all loops, 𝛼would be equal to ||∕10. scenarios (ii) and (iii), Table 2 provides the resulting values of performance indicators (22). As can be seen, the model error introduced by overestimating 𝜂𝑖significantly degrades the system performance and even undermines the advantage of the centralized approach against the decentralized one. Scenario (ii) provided a closer but still suboptimal approximation to the results in case (i), thus also reflecting the benefits of the parameters identification. 4.2. 80-loop plant In the simulations of the 80-loop plant, we consider again a cloudy period of 7 h and the same evolution of the reference and ambient temperatures as before. The resulting overall outlet temperature (see (20)) and HTF flow rates with the proposed coalitional controller in scenario (i) are illustrated in Fig. 9. In contrast to Figs. 4(a), 9indicates the maximum and minimum of these variables among the 80 loops for the sake of clarity. In addition, Table 3 compares the performance indicators in (22) with those obtained in the centralized, decentralized and no-valves cases. Both the results in Fig. 9 and Table 3 are consistent with the performance observed in the simulations of the 10-loop plant in scenario (i). In particular, the coalitional controller approximates again the centralized performance and improves the costs obtained with the decentralized controller. Likewise, it was possible to obtain a much better thermal balance among all loops than in the no-valves case as reflected by costs 𝐽l. As for scenarios (ii) and (iii), a significant decrease of performance was again observed when unknown parameter Fig. 8. Real 𝜂𝑖,ℎ1,𝑖 and ℎ2,𝑖 (dashed lines) and their estimated value (solid lines) for all loops 𝑖. 𝜂𝑖is set to nominal value 𝜂for all loops 𝑖. In this respect, the actual efficiencies were such that 1∕(80||)∑|| 𝑘=0 ∑80 𝑖=1(𝜂𝑖(𝑘) − 𝜂)2= 1.922 ⋅10−4, where ||denotes the number of simulation time steps. Fig. 10 illustrates the formation of coalitions at two different time steps in scenario (i). In particular, the plots axis are associated with the estimated effective irradiance of the loops (𝜂𝑖ℐ𝑖) and their thermal losses (  h𝑖), which are the most relevant features that determine the clustering (see (13)). The markers color indicate the partition obtained implementing the 𝐾-means algorithm and considering 𝑁max cl = 12. With this maximum number of clusters, the mean number of agents in the selected coalitions was 9.1. Regarding the estimation of parameters 𝜂𝑖, ℎ1,𝑖 and  ℎ2,𝑖, the results were also similar to those of the 10-loop plant. After a transient period of about 20 min, the maximum error between the estimated and actual efficiencies was 7.07⋅10−4; and the error on the heat losses coefficients was below 5.28⋅10−8. This larger system with 80 loops allows for a better comparison of the different control approaches in terms of computation time. To this end, consider the following index:  𝑇=1 |c|∑ 𝑘∈c max 𝑗∈(𝑘)𝑡𝑗(𝑘),(23) where 𝑡𝑗(⋅)denotes the time to setup and solve (7) by the 𝑗th coalition in the system. Recall that cis the set of instants in which the flows are updated (see Algorithm 1), and that in the centralized and decentralized cases we always have a single grand-coalition and a set of 𝑁loops singletons, respectively. Also, notice that (23) considers the maximum time among all coalitions since they can optimize their actions in parallel. In scenario (i), the resulting value of (23) was 2.3234 s in the centralized case, 0.6597 s in the coalitional case, and 0.0614 s in the decentralized case. Moreover, Fig. 11 illustrates the time to solve problem (7) for different coalition sizes. As expected, choosing larger clusters increases the number of decision variables involved in (7) and hence the computation times. Finally, notice that the centralized approach requires solving (7) for (𝑘) = at all control time steps 𝑘. That is, 𝑁p⋅𝑁loops (i.e., 400) flow variables were centrally optimized at 840 instants of the simulation. On the other hand, the coalitional Renewable Energy 216 (2023) 118978 9 P. Chanfreut et al. Table 3 Cumulative performance costs in the 80-loop plant simulations. Centr. Coalitional Decent. No-valves Centr. Coalitional Decent. No-valves 𝐽l (i) 2.012 ⋅1062.100 ⋅1062.603 ⋅1062.883 ⋅107 𝐽g (i) 3.813 ⋅1043.818 ⋅1043.969 ⋅1045.677 ⋅104 (ii) 2.025 ⋅1062.353 ⋅1062.616 ⋅1062.883 ⋅107(ii) 3.816 ⋅1043.844 ⋅1043.972 ⋅1045.677 ⋅104 (iii) 5.890 ⋅1065.724 ⋅1066.098 ⋅1063.110 ⋅107(iii) 7.236 ⋅1046.992 ⋅1046.910 ⋅1048.452 ⋅104 Fig. 9. Outlet temperature and HTF flow in the 80-loop plant with the coalitional controller in scenario (i). The gray dotted lines indicate accordingly the maximum and minimum outlet temperature and flow rates among the set of loops. Fig. 10. Clusters formation according to the loops dynamics. Fig. 11. Computation time to optimize the flow rates for different coalition sizes ||. approach only considered centralized computations at 168 instants to solve problem (15), which in turn considers a simplified description of the field with 𝑁top p⋅𝑁cl flow variables (i.e., about 90 flow variables on average). 5. Conclusions A hierarchical clustering-based MPC for optimizing the HTF flow rates in large-scale solar parabolic trough plants has been presented. The proposed approach is built on the idea that the HTF circulating through the loops of solar collectors are controlled in real time by a set of MPC agents. In turn, a supervisor partitions dynamically these agents into cooperative clusters (also called coalitions) and assigns to each of them a maximum flow rate considering the overall conditions and constraints of the solar field. Moreover, the efficiency and heat losses coefficients that determine the loops control model are adapted by using measurements of their outlet temperature. Our results show that the proposed clustering-based MPC can closely approximate the centralized performance, while reducing the need for centralized computations and the coordination efforts. In particular, the loss of performance regarding the overall outlet temperature was under 1.21% in both simulations of a 10-loop and an 80-loop plant. The supervisor allows the MPC agents to take into account the varying conditions among the loops through their constraints adjustment. Therefore, it can restrict the flow that can be pumped to loops covered by clouds to provide more HTF to those receiving higher irradiance. On the other hand, the formation of coalitions among the MPC agents allows for reducing both the communication and the number of variables that they must coordinate, thus also reducing the computation demands. Future research will optimize the outlet temperature reference so as to maximize the net electricity production, as well as studying a fully distributed implementation of the proposed controller. In this regard, we will also explore the potential performance benefits that can be achieved by dynamically adjusting the temperature reference of each of the clusters of loops. CRediT authorship contribution statement Paula Chanfreut: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing – original draft, Writing – review & editing, Visualization. José M. Maestre: Conceptualization, Formal analysis, Writing – review & editing, Supervision, Project administration, Funding acquisition. Antonio Gallego: Conceptualization, Methodology, Formal analysis, Validation, Writing – review & editing. Anuradha M. Annaswamy: Conceptualization, Methodology, Formal analysis, Validation, Writing – review & editing. Eduardo F. Camacho: Conceptualization, Methodology, Formal analysis, Validation, Supervision, Project administration, Funding acquisition. Declaration of competing interest The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Eduardo Fernandez Camacho reports financial support was provided by European Research Council. Jose Maria Maestre reports financial support was provided by Spain Ministry of Science and Innovation. References [1] International Energy Agency (IEA), Renewable electricity, 2022, URL: https: //www.iea.org/reports/renewable-electricity. [2] Renewable Energy Policy Network for the 21st Century (REN21), Renewables 2022: Global status report, 2022, URL: https://www.ren21.net/wp-content/ uploads/2019/05/GSR2022_Full_Report.pdf.