New methods and results in the optimisation of solar power tower plants

Programa de doctorado Matemáticas PhD Dissertation New Methods and Results in the Optimisation of Solar Power Tower Plants Author Thomas Ian Ashley Supervisors Prof. Dr. Emilio Carrizosa Priego Prof. Dr. Enrique Fernández-Cara some text to complete some text to complete some text to complete some text to complete some text to complete some text to complete some text to complete some text to complete some text to complete A todos los que se han implicado en mi formación, tanto profesional como personal, en los últimos tres años. Espero devolver el mismo apoyo que me habéis demostrado. Gracias. Agradecimientos A Thesis is a journey of self discovery that is completed by oneself, but also profoundly aected by the people with whom one is surrounded. For me, these people have been especially important, and here I will attempt to express my gratitude to a select few. En primer lugar, me gustaría agradecer a mis directores Emilio y Enrique, quienes me han guiado tanto durante los últimos tres años. En vez de contratar a un estudiante conocido de la universidad, tomaron el riesgo de contratar a un inglés desconocido, y por eso siempre estaré agradecido. Me han enseñado no sólo a investigar, sino también la cultura española; hasta me he convertido en un nativo de Los Palacios. Quiero agradecer a los trabajadores de la secretaría del IMUS, quienes me han ayudado mucho con los trámites complejos del doctorado, recibiéndome siempre con una sonrisa y sin quejarse nunca por mis dicultades con el español. Llegando a Sevilla con una lengua y cultura nuevas era todo un reto. Tuve la inmensa suerte de ser bienvenido en un grupo de gente increíble en el IMUS, que ha aceptado las diferencias culturales y mi español mediocre, y todo ello me ha hecho sentir como en casa. La Familia IMUS es el logro más importante de mi tiempo en Sevilla, y aunque sois demasiados para mencionaros individualmente a todos, espero que todos sepáis el gran impacto que habéis generado en mi vida; Asun, Víctor, Tana, Nati, Mónica, Sandra, Alba, Moisés, FJ, Vanesa, Marina, Elena, Carlos, Edu, Jose, Marithania, etc. Gracias por los cafés, congresos, risas, estas, bodas (pronto tres) y momentos especiales. Aunque estoy cerrando esta etapa en el IMUS, yo sé que aun quedan muchas más historias memorables con vosotros, y con los próximos miembros de la Familia IMUS (bebé Rodrigo). Por supuesto, quiero agradecer a Las Pueblecitas, Reme y Cristina, por siempre haber estado cuando las he necesitado, aunque hacer un resumen aquí no les haría justicia. Las dos me han apoyado al máximo durante estos tres años, tanto en el trabajo como enseñándome a ser de Posadas. Gracias a Reme por haber aguantado mis preguntas tontas constantes, y por ser la voz de la razón. Especialmente me gustaría dar las gracias a Cristina, ya que ella fue quien me recogió en el parking de la universidad mi primer día, y ha sido una de mis mejores amigas hasta mi último día en el IMUS. IV V Eres la luz de La Familia y una fuente de alegría e imaginación para todos. El IMUS no sería igual sin ti. I must also thank my family, who have raised me to be the person I am today, and have always strived to provide me with the best opportunities to succeed. Although we now live in separate countries and are unable to meet as much as we would like, their continued love and support has helped me immensely over the years of this work. También tengo que agredecer a la familia de Marta, mi nueva familia, que me ha ayudado mucho a acostumbrarme a vivir en España y por aceptarme sin dudas. To Marta I must give the greatest thanks, and beg the most forgiveness. It is thanks to you that I came to Sevilla, and thanks to you that I complete this Thesis. You have supported me through the bad moments, and been there to enjoy all of the good ones. We have grown together, learned together, and ultimately married in this latest etapa of our lives. It is cliche, but true; this part of my life would not have been possible without you, and for that I thank you from the bottom of my heart. Este trabajo ha sido nanciado por:  Optimal Heliostat Fields for Solar Tower Power Plants (PCIN-2015-108). Ministerio de Economía y Competitividad, Spain.  Mathematical Optimization for Data Visualization and Decision Making (MTM2015- 65915-R). Ministerio de Economía y Competitividad, Spain.  Optimización Global. Nuevos Algoritmos y Aplicaciones (P11-FQM-7603). Consejería de Innovación, Ciencia y Empresas. Junta de Andalucía.  Optimización (FQM-329). Junta de Andalucía.  Ecuaciones Diferenciales, Simulación Numérica y Desarrollo del Software (FQM - 131) Junta de Andalucía  COST Action Mathematics for industry network MI-NET (TD1409). European Union.  Research and Innovation Sta Exchange Network of European Data Scientists - NeEDS (H2020-EU RISE 822214). European Union. VI Resumen La tecnología aplicada a la generación de energía renovable ha experimentado un gran avance en las últimas décadas, propiciando un mayor interés en la literatura. Plantas conocidas como Planta Solar de Torre Central (SPT, del inglés Solar Power Tower) son un tipo de tecnología termosolar de concentración (CSP, del inglés Concentrated Solar Power) que continúan en pleno desarrollo por todo el mundo, y consisten en subsistemas que están abiertos a su optimización. Esta Tesis se enmarca en el desarrollo de nuevos métodos y resultados para la optimización de plantas SPT, con un interés particular en la optimización de las operaciones. El Capítulo 1 ofrece información relevante sobre el sector energético actual y continúa con la descripción del diseño y la modelización de una planta SPT. En este capítulo también se describe la teoría óptica que determina la transferencia de radiación incidente en el sistema, y se presentan sus ecuaciones más relevantes. En el Capítulo 2, las operaciones de limpieza del campo de heliostatos están optimizadas por un horario de duración ja, usando programación lineal entera binaria (BILP, del inglés Binary Integer Linear Programming). El problema dimensional se aborda con un algoritmo de agrupamiento, antes de encontrar una solución inicial para el problema de asignación. Por último, se presenta una búsqueda local novedosa mediante técnicas heurísticas que mide el atractivo de una ruta a través del uso de un procedimiento de optimización secuencial por pares que minimiza una medida ponderada de benecio, mientras penaliza la pérdida de energía total. En los Capítulos 3-5 se investiga la estrategia de enfoque adoptada para el campo de heliostatos cuando consideramos un perl de la distribución del ujo deseado, además de incluir restricciones operacionales. En el Capítulo 3, se desarrolla un modelo BILP, donde se seleccionan unos puntos de enfoque predeterminados en el receptor. La función objetiva es lineal y se restringe por equivalencias lineales que se relacionan por una distribución suavizada (para proteger los componentes del receptor de cargas de ujo anormales) a través del uso de una penalización. En el Capítulo 4, se extiende este modelo considerando las variables continuas en el receptor en vez de jadas en puntos de enfoque predeterminados. El resultado es un problema de optimización con una función objetiva no-lineal, no-convexa y con restricciones no-lineales. En este caso, VIII Nomenclature (u, v) Coordinates on receiver surface (ui, vi) Test point i on receiver surface αclu Clustering distance weighting constant αls Local search distance weighting constant αsun Solar angle ¯ H Set of primary heliostats βclu Clustering energy weighting constant βls Local search energy weighting constant βsun Solar angle ˙ p(t) Velocity of aim points at time t  Armijo's Rule parameter γ Step size γ0 Initial step size γk Step size k κ Heliostat eciency L Lagrangian E Euclidean space F Euclidean space P Hilbert space p(t) Collection of aim points at time t XVI Nomenclature XVII Vp Target velocity vector P Projection function µ Penalisation constant Ω Receiver Surface σ1 Velocity constraint threshold σ2 Radiation constraint threshold τ Radiation distribution limit ˜ K Heliostat eciency ~p Receiver normal vector ~v Solar vector ~w Heliostat-Receiver vector ~wh Heliostat-Receiver vector for heliostat h ξ Receiver inclination ζ Degradation function A Weighting parameter C Set of periods in schedule c Convergence limit C∗ Maximum energy constant C∗ Minimum energy constant ce Heliostat eective area cr Heliostat reectance Etar Target distribution Etar u,v Target distribution at point (u, v) f(h, p) Total radiation on receiver for heliostat h aiming at p fatm Atmospheric attenuation fcos Cosine eciency XVIII Nomenclature fref Heliostat reectivity fsp Spillage eciency Fu,v(h, p) Gaussian distribution of heliostat h aiming at p G Balanced objective function g Objective function H Total number of heliostats in eld h Heliostat h I Number of test points i Receiver test point J Combined objective function K Average heliostat eciency M1 Velocity constraint operator M2 Radiation constraint operator P Number of heliostat clusters p Aiming point p ph Aiming point for heliostat h S Size of heliostat clusters s Slack variable Sh Cardinality of cluster with primary heliostat h si Period i of schedule Th Tower height Tol Simulation tolerance Vp Velocity constraint constant Acronyms BILP Binary Integer Linear Programming. CSP Concentrated Solar Power. DNI Direct Normal Irradiance. HTF Heat Transfer Fluid. NRES Non-Renewable Energy Sources. PV Photovoltaic. SPT Solar Power Tower. TES Thermal Energy Storage. UN United Nations. VRP Vehicle Routing Problem. XX Chapter 1 Introduction 2 Introduction 1.1 Energy Sources Energy demand across the globe continues to rise with population growth [Jones and Warner, 2016], and in tandem with our ever-increasing dependence on technology [Sorrell, 2015]. To meet this demand, various natural resources are exploited, including:  Crude Oil  Coal  Natural Gas  Nuclear  Biomass  Solar  Wind  Hydroelectric The current dominant resources used in the production of electricity are fossil fuels (Crude oil, coal, natural gas) [Martins et al., 2018]. In 2016, the authors in Jones and Warner [2016] noted that 91% of worldwide electricity was produced by non-renewable energy sources (NRES) and that 87% was from fossil fuel. Research suggests that whilst there is a signicant amount of crude oil remaining, more than two-thirds of production may need to be replaced by 2030 [Sorrell et al., 2010]. There is much discussion on the theoretical date when the world will reach Peak Oil [Chapman, 2014], that is when the crude oil easiest to extract has been depleted and costs begin to rise, and these uncertainties indicate the need to progress alternative resource adoption. However, the improvement of alternative fossil fuel extraction techniques (for example shale gas extraction [Soeder, 2018]) promises to extend potential supply for many years to come [Helm, 2016]. The environmental cost of fossil fuel use, for example land degradation and emissions (e.g. carbon dioxide and sulphur dioxide), has come to the forefront in the public domain in recent years and cast a negative light on their continued use. Specically, the gradual warming of our planet and its consequent eects on the environment has become a global issue, where it is projected that surface warming increases nearly linearly with carbon emissions [Williams et al., 2017]. 4 1.6. Contributions of this Thesis 11 The cleanliness of a heliostat decreases over time due to particulate deposition, and can be measured from historical data for a particular region and used to create an expected degradation function. In this work we assume a linear degradation function δ=δ(t) , but note that any other function can also be implemented into the model. The optimisation model applied in this work considered a 3 stage heuristic approach to maximise the eciency of heliostats across a xed-length schedule: rst, a clustering optimisation to reduce problem dimensionality is performed; then, the optimisation of the cleaning schedule to produce an initial solution; nally, this initial solution is improved upon through the use of a pair-wise local search heuristic, which takes into account route attractiveness. Heliostat clustering This can be modelled as a p -median problem, where heliostats can be viewed as clients and p potential plants have to be selected, the primary heliostats, which are the central heliostats of the groups. The formulation corresponds to a binary integer linear programming (BILP) problem in which the overall dissimilarity between each heliostat and its associated primary heliostat is to be minimised. The dissimilarity between two heliostats h and h0 is measured as the dierence of the generated potential energies Eh and Eh0 , combined with the physical distance Dist(h, h0) between heliostats h and h0 in a weighted objective function λhh0=αDist(h, h0)+β|Eh− Eh0| where α , β > 0 . In this way, a user can dene relative importance in the clustering phase between energy generation similarity, or proximity of heliostats in clusters. Once the optimal clustering of heliostats has been found, the next stage of the algorithm looks to optimise the cleaning schedule. Schedule optimisation The scheduling optimisation problem is considered with a BILP technique, where the objective is the maximisation of total eciency of the heliostat eld across the schedule. The objective function comprises a summation of integer variables for the allocation of clusters to periods, multiplied by the energy reected by each cluster in the corresponding period, where the degradation function is taken into account. In order to consider the cumulative degradation of eciency across periods, binary linear constraints are included, and are explained in detail in Chapter 2. In addition, linear constraints impose limits on the number of clusters that may be cleaned per day, as well as assuring that each heliostat is cleaned at least once across the schedule. The resulting optimal cleaning schedule for the heliostat eld considers the radiation potential of each heliostat, its degradation of eciency over time and its grouping. With 12 Chapter 1. Introduction this information, the cleaning vehicle operator has a daily schedule of heliostats to clean, where the specic route taken between selected groups must also be considered. Local Search When applied to large heliostat elds, the presented approach may result in optimal cleaning schedules where physically distant groups are allocated to the same day, which may be unappealing to the SPT plant operator. Therefore, a route improvement heuristic is provided, whereby the optimal schedule found in the initial solution is altered whilst considering route attractiveness [Rossit et al., 2016]. The algorithm developed for this heuristic is a rolling optimisation procedure with a penalised integer linear objective function. The objective function minimises the physical distance between the groups in a single period, by permitting swaps from adjacent periods, with a penalty function applied that minimises the loss of eciency caused by swapping groups from the optimal schedule previously found. The solution is the optimal allocation of groups between the two adjacent periods considered, where the weighted energy/distance objective function has been minimised. The optimisation procedure is then repeated over subsequent pairs of periods, until no overall increase in the objective function is found. This local search heuristic looks to improve the schedule in terms of route attractiveness, but does so at a cost to the overall eciency of the cleaning schedule. Such optimisation is scalable to any heliostat eld shape and size, and can be tailored according to specic local constraints (for example impassible terrain) set by the operator. As a nal step, a greedy algorithm is implemented to nd the shortest path for the route taken by the cleaning vehicle for each period of the schedule. The proposed optimisation algorithm is demonstrated on the SPT plant described in Section 1.6.6, producing a cleaning schedule for a 16 day period and highlighting model functionality. 1.6.2 A rst approach to the optimisation of aiming strategies In Chapter 3, the work based on Ashley et al. [2017] is presented, which considers the optimisation of aiming strategies for the heliostat elds in SPT plants. The aiming strategy utilised is known to be important in SPT plant operation [Relloso and García, 2015], and has been a topic of interest in the literature. In Berenguel et al. [2018], the authors approach heliostat aiming with a two-layered optimisation procedure, utilising a genetic algorithm. Further examples can be seen in Besarati et al. [2014] and Astol et al. [2017a]. Such methods often lack a full mathematical description of the problem, and utilise 1.6. Contributions of this Thesis 13 heuristics for simplication. We look to thoroughly model the SPT plant system as an optimisation problem which can be extended to any size and shape plant, and nd the optimal aiming strategy under local weather conditions. Optimisation formulation In this work, a BILP model is presented to optimise the aiming strategy of heliostats in the eld to a set of xed aiming points on the receiver surface. The objective function is formed from the optical model given in Section 1.7.4 that describes the reection of incident solar radiation by the heliostats onto the receiver surface, and calculates the total energy depending upon the aiming strategy chosen. Accurate optical models are essential for simulation of SPT plants and have been analysed extensively, for example Igel and Hughes [1979]; Collado et al. [1986]. The overall objective of an SPT plant operator is to maximise output (and therefore prot), where considerations on storage and current market value of electricity are taken into consideration. Therefore, the aiming point for heliostat h must be chosen so that it maximises radiation reaching the receiver, calculated by the aforementioned optical model. The summation of these individual radiation contributions for all h gives the total radiation reaching the receiver for the given aiming strategy. However, maximisation of the thermal energy reaching the receiver surface does not directly correlate to maximising eciency. Depending on the type of receiver used, higher eciencies are obtained when a certain distribution of energy is held across the receiver [Yu et al., 2014]. This distribution is dynamic, in that it changes with time and conditions of the receiver, and therefore must be taken into consideration within an optimisation procedure. Hence, a set of linear constraints are added to the optimisation model that homogenise the distribution of energy across the receiver surface, whilst the objective seeks the maximal amount of energy. A binary integer linear program is developed with a penalised objective function that assigns heliostats to aiming points at a specic time point, and the Gurobi solver package [Gurobi Optimization Solver, 2019] is used to nd the optimal solution. Finally, an illustrative example is presented for multiple time points across a day utilising the SPT plant data presented in Section 1.6.6, where behaviour of results in relation to solar angles is discussed. 1.6.3 Continuous aiming strategies (I): the stationary case Chapter 4, based upon the article Ashley et al. [2018], concerns the reformulation of the work presented in Chapter 3 to consider the aiming points as continuous variables 14 Chapter 1. Introduction across the feasible region Ω . In this work, the objective again considers both the maximisation of energy reaching the receiver, and the minimisation of deviation from a specic target distribution, where the variables are now continuous across the receiver surface Ω . No pre-dened set of aiming points is used. Instead, each heliostats' aiming point is represented as an individual continuous variable in Ω . In this way, the solution found is more representative of the real-world case, but comes at the cost of higher problem complexity and dimensionality. Optimisation model Due to the use of continuous variables, a BILP technique can no longer be applied, and instead continuous optimisation techniques must be considered. In the case of optimisation of aiming strategies, our objective function is large dimensional, non-linear, non-convex and is subject to convex constraints on the aiming points ph∈Ω for all h in the set H of heliostats in the eld. The objective function, given in Equation (1.6.1) and explained in Chapter 4, is dierentiable across Ω , which permits the use of gradient based optimisation techniques. AX h∈H f(h, ph)−(1 −A)Z"X h∈H Fu,v(h, ph)−Etar u,v #2 dΩ. (1.6.1) The relative weighting parameter A is added to both parts of the objective function, so that an analysis into the relative importance of each may be conducted and the Pareto Front approximated. To solve the optimisation problem, a gradient ascent algorithm is developed that maximises the objective function by computing its gradient. To accomplish this, rapid numerical calculations of the partial derivatives of a non-linear and non-convex function must be performed, for each heliostat in the eld. Considering the aiming points ph to be continuous across Ω introduces integrals across the domain, which for numerical purposes are replaced by nite dimensional approximations in space. Computational time is tackled in the algorithm by the development of a nonstandard step size selection routine, where each variable of the optimisation problem (aiming point) receives independent treatment. The algorithm is written in Python and rapidly calculates the gradient for all heliostats in the eld, and subsequently optimises the objective function. Due to the multi-modal nature of the considered function, it is necessary to incorporate techniques ensuring that the global solution is found, rather than a local maximum. In this work, a multi-start procedure from randomised starting positions is applied. The analysis shows 1.6. Contributions of this Thesis 15 that, for this function, with a relatively low number of multi-starts, a suitable solution is found. To demonstrate the functionality of this algorithm, a numerical illustration is given with realistic SPT plant data dened in Section 1.6.6, where the weighting parameter is varied across multiple runs to generate an approximation to the Pareto Front . Finally, conclusions are drawn and further work is discussed. 1.6.4 Continuous aiming strategies (II): the dynamic case In Chapter 5, the work in Chapter 4 is further extended to a dynamic formulation, where the optimal aiming strategy for a SPT plant is calculated across a time period. In dynamic optimisation, the objective function must be integrated over continuous time, where the unknown is a time-dependent function that indicates the points in Ω aimed at by the heliostats and a discretisation must be implemented to split the problem into manageable subsections of discrete time. The work presented in this chapter, based upon the paper Ashley et al. [2019b], considers a theoretical investigation into various aspects of this dynamic problem. Speci- cally, the existence of a solution is proven, and the optimality conditions of the problem are presented. When applied to the SPT plant framework, the existence of a solution and optimality conditions are again conrmed. In the case of SPT plants, the change to a dynamic problem allows the user to better simulate the physical behaviour of the system. One of the main concerns is the ever-changing meteorological conditions over the heliostat eld, where incident radiation changes due to solar angle shift cause the optimal aiming strategy to change. If local weather conditions are known, and the future radiation gures for the eld can be estimated (from forecasting or historical data analysis), then an optimisation procedure can be applied to optimise future aiming strategies. Physical limitations on the rotational speed of the heliostats restrict the possible movement speed of aiming points that can be allowed. Without limitation, large oscillations could be found in aiming point position over time, which would negatively aect the simulation. Therefore, a constraint is added to the optimisation problem that limits the movement speed of aiming points to a prescribed constant. On the other hand, thermal uctuations across the receiver surface can provoke permanent damage if left unchecked [Relloso and García, 2015]. The optimisation problem looks to minimise the dierence from a desired target distribution over time, where rapid changes may be encountered due to local weather conditions and radiation input. Therefore, it is necessary to include a constraint on the change in radiation on the receiver surface. This is achieved by the discretisation of the receiver into a large quantity of test points, where the change in radiation over time at any point may not be above a certain constant. In this work, two algorithms are presented; in the rst one, we incorporate the 16 Chapter 1. Introduction dynamic constraints on the speed of rotation of the aim points and thermal uctuations as penalty terms in the objective function; in the second one, we apply an Augmented Lagrangian technique. In the penalisation algorithm, a similar approach to that applied in Chapter 4 is utilised. The partial derivatives of the penalised objective function are computed analytically over a time period, and a gradient ascent algorithm is applied to nd the optimal solution. Once the optimal solution to the initial time step is found, it is used as an initial guess for the second time step and the constraints are updated with the new radiation information. As in Chapter 4, the objective function considered here is non-linear, non-convex and formed of continuous variables across the feasible domain Ω . The extension to the dynamic case expands the problem dimensionality, and a second numerical approximation in the time variable is performed. In the Augmented Lagrangian algorithm, the constraints are considered together with the objective function in a duality-penalty method. As before, a gradient ascent approach is used, where the gradient is calculated analytically. Again, a discretised time period is considered, and a recurring optimisation procedure is implemented to calculate the global optima across all time steps. Both algorithms are implemented in Python and are demonstrated with an illustrative example utilising realistic SPT plant data from Section 1.6.6, with radiation information for the duration of a day. Finally, conclusions on the algorithms presented are drawn, and intentions of future work are discussed. 1.6.5 Inclement weather eects on optimal aiming strategies Chapter 6 investigates in more depth the aiming point strategy of a heliostat eld, considering the eects of inclement weather. Transient clouds not only have a profound impact on the eciency of an SPT plant, but also pose a threat to the integrity of the receiver components, as rapid thermal uctuations can cause permanent damage [Martínez-Chico et al., 2011]. A cloud may be modelled according to the location of its shadow on the heliostat eld, with information on its approximate shape and size. The loss of eciency associated with cloud cover on a heliostat is proportional to the density of the cloud, and is dependent on weather conditions. When considering the aiming strategy of an SPT plant across a day, the varying solar radiation input is taken into account, and the strategy adjusted accordingly. If information on local cloud cover is available (see [Lopez-Martinez and Rubio, 2002; Martínez-Chico et al., 2011] for related research), from observation or forecasting, then the optimisation of aiming strategies can be updated to include eciency losses due to cloud cover. 1.6. Contributions of this Thesis 17 In this work, a BILP problem was developed that optimised the aiming strategy for a heliostat eld at a given time point, when cloud information was known. With this model, the aiming strategy was optimised in order to approach a target distribution of energy on the receiver, whilst also preventing large thermal gradients between test points. Due to the uncertainty of cloud cover, a stochastic linear integer programming technique was applied, where the cloud characteristics (location, size and density) were considered to be stochastic in nature. A variability percentage was allowed for these characteristics over a set of scenarios, and the optimisation problem updated to account for the uncertainty. The resulting optimal aiming strategy then describes the best solution when considering cloud characteristic uncertainty. This methodology is illustrated with an example utilising realistic SPT plant data from Section 1.6.6, and assumed cloud characteristic parameters. 1.6.6 An illustrative example The work presented in Chapters 2-6 of this Thesis develops optimisation models for various aspects related to SPT plant operation and design. In each chapter, the functionality and performance of the designed algorithms are presented using an illustrative example, where realistic data for a particular SPT plant is used. In this section we summarise the SPT plant chosen (PS10 in Seville, Spain [Abengoa, 2019]) and the imposed assumptions. This was the rst operational commercial SPT plant and therefore has an abundance of data available, and is frequently used in the literature. The general specications for the PS10 SPT plant are given in Table 1.1, and the layout of the 624 heliostats in the eld is shown in Figure 1.1. In Figure 1.2 the reected radiation reaching the receiver for each heliostat, averaged across the design point day, is given. This averaged value is utilised in Chapter 2 in the optimisation of cleaning schedules, where the contribution across time of a heliostat must be considered. This averaged assumption is a sensible simplication for the purposes of optimisation model design, but would benet from further investigation at a later date. The incident radiation on the heliostat eld is assumed equal for all heliostats, which is a common assumption in the literature due to almost negligible changes across short distances. The radiation prole for the design point day, noted in Table 1.1, is given in Figure 1.3. It is worth noting that Figure 1.2 shows a distinct dierence in the contribution to the system across heliostats, which demonstrates the need for optimisation of heliostat management. The receiver of the PS10 SPT has been modelled as a circular at plate in this Thesis, as a convenient simplication for the model. However, it should be noted that 18 Chapter 1. Introduction all models developed in this work are directly applicable to any other receiver shape or size. Attribute Value Note General Cost e35 m Capacity 11 MW Solar hours per day 9−15 320 days per year Location 37◦2603200N 06◦1501500W Latitude/Longitude Design point day 21st March To calculate incident radiation Direct normal irradiance (DNI) 823.9W/m2 Approximated at midday Tower/Receiver Tower Height 115 m Receiver Height 6 m Approximated Receiver Width 6 m Approximated Receiver Type Cavity Receiver Inclination 12.5◦C Approximated HTF Water Storage Yes 1 hour Heliostat Number of Heliostats 624 Heliostat size 120m2 Table 1.1: PS10 data Figure 1.1: PS10 Heliostat Locations Figure 1.2: Average reected radiation per heliostat 1.7. Technical Background 19 Figure 1.3: Design point radiation prole 1.7 Technical Background In this section, we present a technological background of SPT plant subsystems and relevant radiation modelling information. As introduced previously, an SPT plant is formed of a eld of reective surfaces, known as heliostats, which concentrate incident solar radiation onto a centrally located receiver, mounted atop a tower. The thermal energy generated by this process is transported by a HTF to a generator, whereby electricity is produced. As demonstrated in Section 1.5, there are various subsystems of an SPT plant design that can be improved upon by mathematical optimisation, and what follows is a brief summary of the properties and equations that are important in an optimisation model of these subsystems. 1.7.1 Tower In an SPT plant, the tower is usually the central focus of the plant, which supports the receiver at a designated height Th . The gathered heat from the receiver is transported by the HTF via a system of tubes in the tower to the next stage of the cycle (either directly to the generator or to storage). The height of the tower, Th , inuences the angle to which the heliostats must orientate themselves in order to maintain solar focus on the receiver, and also inuences the optimal design of the heliostat eld. 20 Chapter 1. Introduction Tower design has been researched in terms of optimisation [Carrizosa et al., 2015a], but for the purposes of this Thesis, has been assumed prescribed. 1.7.2 Receiver The receiver constitutes one of the most important elements of an SPT plant, and its design inuences the type of heliostat eld that will be employed. For the purposes of the work in this Thesis, a at-plate circular receiver has been assumed. However, as previously mentioned, it should be noted that this was chosen for convenience, and the work presented is directly extendible to other receiver shapes and sizes. For a at plate receiver (and by extension any shape or size receiver) important properties include its dimensions, which directly inuence the amount of radiation captured and spillage incurred. The feasible region for the optimisation problems developed in this Thesis, Ω , is dened as the continuous region contained within the circle r2=u2+v2 centered at the origin of the (u, v) plane on the receiver surface. Here, r is the given value of the radius of the circular receiver, where r=1 2RecW and RecW is the receiver width. The inclination of the receiver, ξ , aects the distribution of radiation on the receiver surface. This is assumed constant in this work, however it is noted that in a coupled optimisation of aiming strategies and receiver design, this parameter will be of importance. The materials and internal design of the receiver inuence the transport of thermal energy from the surface to the HTF in the interior. Aside from maximising durability of the receiver (high thermal ux causes rapid deterioration [Relloso and García, 2015; Sánchez-González et al., 2015]), this also has implications on the required optimal radiation distribution. 1.7.3 Heliostat Field The shape of the heliostat eld largely depends upon the local topography, and the type of receiver employed. However, in almost all cases in existing SPT plants, a patternbased approach has been used to form the layout of the heliostats. Examples of this approach include radially staggered, sunower and spiral [Mutuberria et al., 2015]. Pattern-based elds are the easiest to develop and construct, but are not optimal in the sense of maximising eciency. Non-pattern-based heliostat elds have been shown to improve overall eciency of an SPT plant [Carrizosa et al., 2015a], and on a caseby-case basis could be used in the design of future plants. For the purposes of this Thesis, we assume the concentric circle pattern utilised in the PS10 SPT plant, as described in Section 1.6.6. In a South facing eld (with a North facing receiver) the heliostats are dened with 1.7. Technical Background 27 where f5(φi, x, y) = (1 + f2 4(x, y)) + cos φi(1 −f2 4(x, y)) 2||~w||2. (1.7.16) see [Collado et al., 1986; Collado, 2008] for more information. The total radiation captured by the receiver for heliostat h∈H at a specic time can then be calculated by the product of eciency losses and radiation distribution: f(h) = Rad ·fsp(h)·fatm(h)·fref (h)·fsb(h)·fcos(h), (1.7.17) where Rad is the level of incident solar radiation on the heliostat eld. Chapter 2 Heliostat eld cleaning scheduling for Solar Power Tower plants: A heuristic approach 28 Heliostat eld cleaning scheduling for Solar Power Tower plants: A heuristic approach 2.1 Introduction Accessibility is a driving factor in the use of standard structured heliostat eld designs, such as the radially staggered pattern [Sánchez and Romero, 2006]. Structured designs allow easy access to all heliostats for maintenance, and also partially pre-determine the way in which the heliostat eld will be cleaned by their structure of rows of heliostats. On the other hand, utilising a pattern-free eld design is shown to increase the eciency of the heliostat eld [Carrizosa et al., 2015a; Cruz et al., 2018b], but complicates the accessibility of heliostats within the eld for both maintenance and cleaning. Heliostats are required to be cleaned regularly, as accumulation of dust and foreign debris will lower the reectivity of the mirror and therefore lower the eciency of the solar plant [Roth and Pettit, 1980; Sarver et al., 2013]. Various strategies have been developed for cleaning heliostat elds [Fernández- García et al., 2013]. However, the most widely implemented method is the use of a vehicle with a cleaning arm, which cleans heliostats with a mixture of water and brushing. Cleaning all heliostats daily is impractical due to the number within the eld and water scarcity is common in regions with high solar radiation. The frequency of cleaning for each heliostat is partially determined by its physical location within the eld, as the energy generated by the heliostat is strongly dependent on its location, see Figure 1 in Carrizosa et al. [2017], and heliostats that provide more energy are of more importance to be kept clean. Moreover, the proximity between heliostats and to structures will cause shielding from wind, which can aect dust deposition [Singh et al., 2015]. The vehicle used for the cleaning activities in a CSP plant has a limited water carrying capacity, and will need to return to the water depot once empty. This limited water capacity, and length of time taken to clean each individual heliostat, will determine a maximum number of heliostats able to be cleaned in a certain cleaning period. 30 2.1. Introduction 31 Therefore in order to clean an entire eld of heliostats, a cleaning schedule is desired, where the objective is to maximise the overall amount of incident energy reected onto the receiver, under constraints on the number of heliostats cleaned per time period. The goal of this chapter is to optimise the cleaning schedule for CSP plants with any size and shape heliostat eld. To do this, we will consider the optical eciency of the heliostats in the eld, and also the subsequent routing problems for each period in the schedule. Whilst this study is applicable to heliostat elds in any general solar power plant, the case of most interest is in an SPT plant. In this type of solar plant, the heliostat eld layout generally has a more complex geometry and larger size, causing the optimal scheduling of cleaning activities to be complicated. The global cleaning strategy combines an allocation problem, whereby heliostats are allocated to particular cleaning periods, with a routing problem for each period. These problems are not independent, as the route taken by a cleaning vehicle in a particular period will aect the heliostats allocated across the cleaning schedule. The main innovative contribution of this work, summarised in the pseudocode in Figure 2.1, is precisely the way this complex scheduling-routing problem is addressed: we rst perform a clustering analysis to divide the eld into homogeneous groups, then we determine the optimal cleaning schedule, and nally a local search is performed to improve this sequential solution. Figure 2.1: Clustering Pseudocode The structure of this chapter is as follows. In Section 2.2, degradation in optical 32 Chapter 2. Heliostat eld cleaning scheduling for Solar Power Tower plants: A heuristic approach eciency due to soiling is investigated; a grouping optimisation problem is presented, followed by an optimisation procedure to nd the cleaning schedule. A subsequent local search heuristic is then applied to improve the obtained solution. Section 2.3 illustrates the presented approach with real SPT data given in Section 1.6.6. Finally, conclusions are given in Section 2.4. 2.2 Modelling In this section, the eect of heliostat soiling on eciency of energy transfer is detailed and the optimisation problem described. The incident solar radiation reected by each heliostat h at any given time instant t , labelled f(h, t) , is calculated and stored in a preprocessing step. The method used to model the movement of incident radiation within the system is detailed in Section 1.7. 2.2.1 Eciency Degradation If heliostat soiling is not considered, the energy reected by heliostat h at time t is f(h, t) and thus the energy generated in a period of length T , [0, T] is RT 0f(h, t)dt . In order to consider heliostat soiling, a degradation function ζh(t) for heliostat h at time t must be dened, where the optical eciency and therefore the energy generated by a heliostat will be reduced over time, unless cleaning operations are performed. The degradation function chosen directly models the eect of soiling on the optical eciency of a given heliostat and is dependent on local weather conditions. In this study, the degradation function is assumed to be linear in time. Due to both routing and environmental costs, not all heliostats should be cleaned every day. If the analysis is performed for a single heliostat within the eld, which is to be cleaned r times during a cleaning schedule, one can determine the optimal periods to perform cleaning operations in order to maximise optical eciency of the heliostat. Ignoring routing issues, the solution is simple: the cleaning operations should be done periodically. Indeed, if we denote the eciency of a xed heliostat by κs , where s is the allocated cleaning period, we can dene the total eciency across time as: Ks:= Zs 0 κtdt, (2.2.1) where K0≡κx is decreasing and K is therefore a concave function. We can nd the total eciency of a heliostat, across a cleaning schedule with multiple cleaning instants, by calculating: ˜ Ks1, ..., sr=Zs1 0 κtdt +Zs2 0 κtdt +... +Zsr 0 κtdt, (2.2.2) 2.2. Modelling 33 Our aim is to maximise the eciency of the heliostats across the cleaning schedule: Maximise ˜ Ks1, ..., sr, (2.2.3) subject to: s1, ..., sr≥0, (2.2.4) r X i=1 si=T. (2.2.5) Since E is a concave function 1 r˜ Ks1, ..., sr=1 rPr i=1 K(si)≤K(T r) which implies ˜ Ks1, ..., sr≤rK(T r) = ˜ KT r,T r, ..., T r . Therefore, the best strategy is to perform periodic cleaning operations. While the previous analysis shows that the optimal cleaning strategy for one single heliostat is given by a periodic schedule, the problem to be addressed involves many heliostats. In the next section, we formulate the clustering procedure as an optimisation problem for the allocation of heliostats to groups. 2.2.2 Clustering heliostats Addressing the scheduling problem and subsequent routing problems may be unmanageable, since the considered heliostat eld may contain thousands of units. For this reason it is worthwhile applying a clustering strategy [Grotschel and Wakabayashi, 1989; Hansen and Jaumard, 1997] before the optimisation of a cleaning schedule, reducing problem size and complexity. We look to cluster the heliostats in the eld into distinct groups. Groups should not be geographically disperse, since this would increase routing costs. Moreover, since the cleaning schedule will be identical for heliostats in the same group, such heliostats should generate a similar amount of energy. For this analysis, the energy f(h, t) generated by a heliostat h is averaged along one day and labelled Eh , as shown in Figure 2.2 for the PS10 SPT plant in Sanlúcar la Mayor, Seville [Abengoa, 2019]. A dissimilarity function is then introduced: for any pair of heliostats, h, h0 , let λhh0 denote the dissimilarity between h and h0 , given by λhh0=αcluDist(h, h0) + βclu|Eh− Eh0| , where αclu, βclu >0 are given constants and Dist(h, h0) is the physical distance between the heliostats h and h0 . The choice of constants αclu and βclu determines the importance of distance and energy in the clustering optimisation. These values are chosen according to the plant being modelled and the interest of the user in maximising energy or minimising distance. 34 Chapter 2. Heliostat eld cleaning scheduling for Solar Power Tower plants: A heuristic approach Figure 2.2: Average reected radiation per heliostat This clustering can be expressed as a p -median problem, as in Avella et al. [2007]; Daskin and Maass [2015]; Mladenovi¢ et al. [2007]. Heliostats can be viewed as clients and p potential plants have to be selected, the primary heliostats, which are the central heliostats of the groups. This is modelled as a binary linear problem in which the overall dissimilarity between each heliostat and its associated primary heliostat is to be minimised. More precisely, let us denote by H the set of heliostats in the SPT eld. For any heliostat h∈H , let us set yh=   1, if h is a primary heliostat 0, otherwise. (2.2.6) For heliostats h and h0∈H , let us set ahh0=   1, if h0 is allocated to group with primary heliostat h 0, otherwise. (2.2.7) We then constrain the optimisation problem by requiring that each heliostat h may only be allocated to one group, with primary heliostat h : X h∈H ahh0= 1 ∀h0∈H. (2.2.8) We also limit the number of clusters to a constant value P : X h∈H yh=P. (2.2.9) Finally, we set a limit, S , on the number of heliostats allocated to each group. 2.2. Modelling 35 X h0∈H ahh0≤Syh∀h∈H. (2.2.10) Finally, the objective function can be written in the form X h0∈HX h∈H λhh0ahh0, (2.2.11) which is to be minimised. 2.2.3 Schedule Optimisation Once heliostats have been clustered into groups, as detailed in Section 2.2.2, we allocate those groups to cleaning periods in order to optimise the overall cleaning schedule. The objective of the optimisation procedure is to maximise the total energy generated over the schedule duration, whilst considering the allocation of heliostat groups to cleaning periods and the subsequent degradation of eciency of each heliostat. Using the reected energy from each heliostat in the eld (which is calculated in a preprocessing step) and the assumed degradation function, we maximise the total energy reected by all heliostats across the cleaning period. In what follows we show how to address the scheduling problem by using Mathematical Optimisation. We denote by C the set of periods and ¯ H the set of primary heliostats, obtained by solving the p-median problem described in Section 2.2.2. Moreover, we dene the binary variable xhc , where: xhc =   1, if group with primary heliostat h is cleaned in period c 0, otherwise. (2.2.12) for any h∈¯ H and c= 1, ..., C . In order to determine the loss of eciency for a specic period, it is necessary to also dene the binary variable zhcr , which accounts for whether a heliostat group has been cleaned in the previous periods with relation to the current period. Thus, for any h∈¯ H , c∈C and r= 0, ..., c , we set: zhcr =   1, if group with primary heliostat h was cleaned r periods before period c 0, otherwise. (2.2.13) For instance, if we are considering whether or not to clean heliostat group 3 during a schedule of length 2, we would include the variables z310 , z311 , z320 , z321 and z322 . 36 Chapter 2. Heliostat eld cleaning scheduling for Solar Power Tower plants: A heuristic approach By considering the energy generated by a heliostat with full optical eciency, subtracting the amount of energy lost per period it is not cleaned, Ehcr , and summing for all heliostats over all periods, one can calculate the overall energy generated by the eld. Using the binary variables zhcr , our objective is then X h∈¯ HX c∈CX r Ehcrzhcr, (2.2.14) to be maximised. We constrain this objective function by requiring that each heliostat group must be cleaned at least once across all periods: X c∈C xhc ≥1∀h∈¯ H. (2.2.15) This constraint is applied, as it is assumed that all groups of heliostats in the eld provide signicant energy to the system, and must all be maintained in the long term. In the case where a heliostat eld has not been optimised, it may be of interest to not enforce this constraint, and investigate whether or not some heliostats are never cleaned and therefore may not be of signicant value to the overall system. We also assume that each cleaning route may clean at most τ heliostats: X h∈¯ H Sh·xhc ≤τ∀c∈C, (2.2.16) where Sh denotes the cardinality of the group with primary heliostat h . This constraint models the water carrying capacity of the cleaning vehicle, as well as the time taken to complete a cleaning period, and it may be changed according to each particular SPT plant specication. We also introduce the following constraints: zhcc ≤1−xhi ∀i= 1, ..., c, ∀c= 1, ..., C, ∀h∈¯ H, (2.2.17) where zhcc will be zero if heliostat group with primary heliostat h is cleaned in any period; zhcr ≤xh(c−r)∀r= 0, ..., c −1,∀c= 1, ..., C, ∀h∈¯ H, (2.2.18) where zhcr will be zero if heliostat group with primary heliostat h is not cleaned r periods before period c ; zhcr ≤1−xh(c−r+1) ∀r= 1, ..., c −1,∀c= 1, ..., C, ∀h∈¯ H, (2.2.19) where zhcr will be zero if heliostat group with primary heliostat h is cleaned in the 2.4. Conclusions 43 Figure 2.8: Optimised Cleaning Schedule Figure 2.9: Calculated Route Period 2 44 Chapter 2. Heliostat eld cleaning scheduling for Solar Power Tower plants: A heuristic approach typical cleaning technology and assumptions on the eciency degradation of heliostats due to soiling mechanisms. The presented method provides an increase of nearly 5% in total energy of the schedule, when compared to an assumed cleaning schedule which followed the rows of the heliostats in the plant. Whilst an increase in total energy was found, the real benet for the method would be found in the design phase of an SPT plant, where the location of heliostats and their access roads could take into account the optimal cleaning schedule. A p -median type linear integer program was developed to perform a clustering analysis, reducing large problem spaces by nding the optimal grouping strategy of the heliostats, considering a weighted objective function of physical distance and energy prole between heliostats in the eld. The heliostat grouping is then used to obtain an initial cleaning schedule of the heliostat groups over a period of time, in order to maximise the energy reaching the receiver. In our numerical illustration, this program was presented for a schedule of 16 periods, with a route limit of 52 heliostats per period, where the results show an optimal energy prole, but a sub-optimal routing solution for each period. The initial solution is then rened with a local search heuristic, which pairwise swaps groups of heliostats between consecutive days in the schedule. Each move amounts to solving a linear integer problem, where operational costs are reduced in the routing problems, at a cost of total energy gained during the schedule. The optimisation procedure developed in this work can be utilised by current SPT plant operators in order to optimise their cleaning routes. It can also be utilised in the planning phase of a new SPT plant, where pattern-free heliostat elds will require more complex cleaning routes, which need to be optimised to ensure maximum energy generation. The application of Operations Research techniques to SPT plant design and operation has many possibilities, where this work could be extended to include; time dependency if cleaning operations are conducted during the day, routing problems with depots due to water carrying capacity of cleaning vehicles, stochastic processes for weather events which have varying eects on the soiling of heliostats, and large scale problems with multiple towers and vehicles. The length of the schedule optimised in Section 2.3 was 16 periods, which was chosen in order to present a study of an interesting length. This choice will aect the result of the optimisation procedure, due to variable local weather conditions, and further research of interest is to optimise the schedule length chosen, as well as the schedule itself. The example route shown in Figure 2.9 does not follow the rows formed during the construction of the SPT plant, and assumes that the cleaning vehicle can navigate between heliostats. In the case of the PS10 SPT plant, the application of this route 2.4. Conclusions 45 is possible due to spacing between heliostats and the ground topology. However, some SPT plants have obstructions, both natural and man-made, within their elds. These obstructions will limit the paths that the cleaning vehicle can take, and must be taken into account when optimising the cleaning schedule. This can be implemented in the method developed in this study with additional constraints on the feasible region. This investigation could be used to inuence future design of SPT plants, in order to maximise their energy generation and reduce overall costs, by aiding in the design of cleaning schedules whilst adapting the layout of the heliostat eld using techniques as in Carrizosa et al. [2015a]. The results presented use a novel swapping approach for pairs of periods in the schedule. However this is directly extendible for any number of periods at a time. Another extension to this work could be the investigation into the eects of optimal cleaning schedules, where we look to maximise optical eciency of the heliostats, against optimal aiming strategies as developed in Ashley et al. [2017]. In the presented work, the resulting cleaning schedule consists of a set of heliostats to be cleaned during each period. The routing problem for each period has been solved utilising a greedy algorithm, however an improvement to this study could consider a dierent approach. The degradation of heliostat eciency over time with respect to soiling can have large impacts on overall SPT plant performance [Sarver et al., 2013]. Many factors inuence the exact degradation curve, including geographical location of the SPT plant and proximity of heliostats to one another. The eect of realistic degradation information on heliostat eciency when combined with an optimisation procedure would allow plant operators to make the best decisions when approaching cleaning operations. For the purposes of this study, a linear degradation curve has been assumed, although a stochastic case with weather events was also investigated and should be implemented in future work. Finally, as technology advances, more potential cleaning techniques become available, such as robots. Autonomous robots are being developed to clean heliostat elds, and the optimisation of their operations could further extend the model considered in this work. Chapter 3 A rst approach to the optimisation of aiming strategies 46 A rst approach to the optimisation of aiming strategies 3.1 Introduction The aiming strategy commonly used in research into the optimisation of SPTs assumes that all heliostats in the eld aim at the centre of the receiver, see for example Carrizosa et al. [2015a], and a summary of various recent SPT optimisation techniques by Barberena et al. [2016]. This assumption allows for easier computation of the ux distribution across the receiver surface and reduces complexity of the adjustment of the heliostats. Using a central aim point for all heliostats leads to a large heat ux at the centre of the receiver and large ux gradients towards the edge of the receiver, which can cause strong heat loads and potentially lead to damage over time and costly repairs [Lopez- Martinez and Rubio, 2002; Relloso and García, 2015; Salomé et al., 2013; Sánchez- González et al., 2015; Yu et al., 2014]. An uneven ux distribution across the receiver surface also lowers the eciency of the energy transfer to the thermal uid within the receiver [Fend et al., 2004; Yu et al., 2014]. Therefore, maintaining an even distribution will increase the eciency and allow for greater energy production. Some research has been conducted where more complex aiming strategies are considered for dierent receiver types [Astol et al., 2016; Belhomme et al., 2013], as well as closed-loop feedback mechanisms to provoke changes in aiming strategy [Berenguel et al., 2004; Kribus et al., 2004]. Applications of alternate optimisation algorithms for the aiming strategy have also been exhibited [Salomé et al., 2013; Sánchez-González et al., 2015; Yu et al., 2014; Cruz et al., 2016] and a summary of optimisation techniques collected [Baños et al., 2011]. The distribution of reected energy from the heliostats onto the receiver is assumed to be a Gaussian, as detailed in Section 1.7. The distribution can be written in the form: f1(x, y)exp−f2(u, v, x, y) 2f2 3(x, y), (3.1.1) 48 3.2. Changing Aim Point 49 In order to reduce damage to receiver components and to optimise the energy reaching the receiver, multiple aim points across the receiver surface in a grid can be considered, as depicted in Figure 3.1: a set of possible points within the receiver is given, and the solver must choose, for each heliostat and time instant, the most appropriate aiming point. -3 -2 -1 0 1 2 3 Receiver Width (m) -3 -2 -1 0 1 2 3 Receiver Height (m) Figure 3.1: Aim Point Grid The eld of heliostats would then be split, so that a proportion of them aim at each point. The result of this would be to smooth the ux gradient across the receiver surface and to reduce the peak heat ux at the centre. Another eect of this strategy would be to increase the amount of energy that is being lost due to spillage, where aiming towards the edge of the receiver causes some of the energy to miss completely. The goal of the study in this chapter is to design an aiming strategy that minimises the ux gradient across the surface of the receiver, whilst minimising the spillage and maintaining a minimum amount of energy. 3.2 Changing Aim Point The spillage eciency fsp for a particular heliostat is calculated as the integral of the distribution given in (3.2.1) across the receiver surface, as detailed in Section 1.7.4. fsp(x, y) = f1(x, y)Z2π 0Zr 0 exp−˜ f2(ρ, φ, x, y) 2f2 3(x, y)ρ dρ dφ , (3.2.1) where r is the radius of the circular receiver and ˜ f2(ρ, φ, x, y)≡f2(ρcos φ, ρ sin φ, x, y) . In order to model the distribution of energy across the receiver surface when we consider an aiming point oset from the centre, the bounds of integration in (3.1.1) 50 Chapter 3. A rst approach to the optimisation of aiming strategies need to be altered, as shown below. Consider the Gaussian distribution of the reected energy as centered at the origin, and that the receiver is oset by the aim point oset. fsp(t, x, y) = f1(t, x, y)Z2π 0Zˆr(φ) 0 exp −˜ f2(ρ, φ, x, y) 2f2 3(t, x, y)ρ dρ dφ , (3.2.2) where ˆr(φ) is given by the solution of the equation R2= ˆr2+r2 0−2ˆrr0cos(θ−φ) We shall assume that the heliostat is always pointing within the receiver boundaries, if it is to be included in calculations (i.e. unless it is turned away for safety or other reasons and is then not producing output). The spillage eciency in (3.2.2) then becomes: fsp(t, x, y)≈2π 10 f1(t, x, y)f2 3(t, x, y) 10 X i=1 1 f5(φi, x, y)1−exp−ˆr(φi)2f5(φi, x, y) 2f2 3(t, x, y), (3.2.3) where, f5(φi, x, y) = (1 + f2 4(x, y)) + cos φ(1 −f2 4(x, y)) 2||~w||2. (3.2.4) Multiplying the spillage eciency, fsp(t, x, y) , by the solar radiation at a specic time point and by other losses inherent in the system (as detailed in Section 1.7) gives us the total energy that reaches the receiver at time t . This value is for one heliostat, aiming at a specic aimpoint. 3.3 Optimisation In this section our goal is to develop an optimisation model for the selection of aiming points for the heliostats in a xed eld, in order to maximise the energy produced under some homogeneity constraints on radiation. There are various optimisation techniques that could be applied to solve this problem, including heuristic methods such as Ant Colony or Genetic Algorithm techniques, however the BILP technique was implemented in this case. BILP techniques lead to the optimal solution within a nite time period, whereas this is not guaranteed with purely heuristic methods and is therefore not the best choice for this problem. Applying the BILP technique with a heuristic time limit allows a solution to be obtained quickly, whilst the gap between the obtained solution and the optimal solution will be between a dened upper bound. Changing the aim point of a heliostat from the centre of the receiver to another point aects the amount of energy reaching the receiver by changing the slant range 3.3. Optimisation 51 (that is, the modulus of ~w ), the spillage eciency, the cosine eciency and the light distribution across the surface. Let Ap be the set of aiming points on the receiver surface, where Ap∈R2 and let H be the set of heliostats aiming at the points in Ap . Let us dene an optimisation procedure for any xed time instant, t . For h∈H , a∈Ap , set zha to the boolean variable dened as: zha =   1, if heliostat h is allocated to aiming point a 0, otherwise . (3.3.1) The reected radiation pattern Fb(h, a) is the radiation point value at aiming point b received from heliostat h aiming at aiming point a . The total integrated radiation, f(h, a) , is the total radiation received across the receiver from heliostat h aiming at aiming point a . We therefore look to maximise the total incident energy on the receiver: Maximise X a,h f(h, a)zha. (3.3.2) We constrain this objective function by requiring that no heliostat may be looking at more than one aiming point on the receiver, but may be stowed in case of high winds or potential damage to the receiver or mirror. This gives the constraint: X h∈H zha ≤1∀a∈Ap. (3.3.3) We also constrain the received energy at the aiming points: C∗≤X h∈H a∈Ap Fb(h, a)zha ≤C∗∀b∈Ap, (3.3.4) where C∗ is a xed maximum energy and C∗ is a xed minimum energy. These constraints prevent the receiver being subject to excessive temperatures (which could cause permanent damage) and also ensure that a minimum amount of energy is being collected at each aiming point. In order to approximate a uniform distribution of energy across the receiver, we will look to also constrain the range of energy received between any two aiming points by imposing max a(X h,a Fi(h, a)zha)−min a(X h,a Fj(h, a)zha)≤τ , (3.3.5) 52 Chapter 3. A rst approach to the optimisation of aiming strategies where τ is a given constant. This can also be written in the form max a(X h∈H a∈Ap Fi(h, a)zha)≤τ+ min a(X h∈H a∈Ap Fj(h, a)zha), (3.3.6) which is equivalent to the following set of linear constraints: X h∈H a∈Ap Fi(h, a)zha ≤τ+X h∈H a∈Ap Fj(h, a)zha ∀i, j ∈Ap, with i 6=j (3.3.7) The optimisation problem to be solved at each time instant t can then be summarised as follows: Maximise X h∈H a∈Ap f(h, a)zha Subject to: X h zha ≤1∀a∈Ap, C∗≤X h∈H a∈Ap Fi(h, a)zha ≤C∗∀i∈Ap, X h∈H a∈Ap Fi(h, a)zha ≤τ+X h∈H a∈Ap Fj(h, a)zha ∀i, j ∈Apwith i 6=j zha ∈ { 0,1 } ∀h∈H, ∀a∈Ap The optimal aiming strategy for an SPT plant is dependent on the time of day, as well as the day of the year. For an optimal aiming strategy to be achieved, it must be optimised at a rate that will capture the changing radiation pattern over time. In terms of the optimisation problem, both the objective function and the constraints will change as functions of time, caused by the variable incident radiation and physical constraints on the SPT plant. A rapid change in incident radiation at one point on the receiver surface, for instance caused by passing clouds, could potentially cause damage, indicating that frequent updates to the aiming strategy are needed. As there is no guarantee that the optimal aiming pattern will remain optimal, or even feasible, over time, the optimisation procedure must be repeated frequently, using knowledge of local weather to constrain the problem in real-time. Continuous aiming strategies (I): the stationary case 4.1 Introduction In this chapter, as in Chapter 3, we develop an optimisation model for the aiming strategy of an SPT plant. In Chapter 3, a discrete approximation of the receiver surface Ω into a set of potential aiming points is used, thereby reducing problem dimensionality and allowing a BILP technique to be applied. A more accurate approach is to allow the aiming point variables to take any value and be continuously dierentiable across the receiver surface. This is equivalent to dening a separate aiming point for each heliostat in the eld, within the set created by the boundaries of the receiver dimensions. This allows the simulation of SPT plant operations to better reect real-life conditions, where aiming points for heliostats will not be limited to pre-set locations. This will be the viewpoint (and the main contribution) in this chapter. Whilst the rst objective is the maximisation of total radiation reaching the receiver surface, operational limitations of the SPT plant must be taken into consideration. As previously noted, these include inhomogeneous heating of the receiver surface, where large thermal uxes can cause non-optimal energy generation, or even permanent damage to the receiver components. Hence, the objective in this chapter is to maximise the radiation captured by the SPT plant receiver, whilst taking into account the deviance from a desired radiation distribution across the receiver using continuous optimisation models. The two objectives will be combined into one single non-convex, non-linear criterion via additive weighting. By varying the weights, an approximation to the corresponding Pareto Front will be obtained. Some related heuristic methods can be found in Wagner and Wendelin [2018]. The structure of this chapter is as follows: Section 4.2 details the problem formulation; Section 4.3 describes how to numerically formulate the problem; Section 4.4 provides an illustrative example of the construction of the Pareto frontier of the biobjective problem, and Section 4.5 presents concluding remarks. 60 4.2. Model 61 4.2 Model The rst component of the criterion in our aiming strategy optimisation model is the radiation generated by the heliostats on the receiver surface. As in Ashley et al. [2017], the radiation passing through the system is modelled using a Gaussian distribution on the receiver, a non-empty bounded open convex set Ω⊂R2 . We will denote by |Ω| (resp. ¯ Ω ) the measure (resp. the closure) of Ω . In this work, any heliostat h∈H will be required to aim at point ph∈¯ Ω . The value of this distribution evaluated at a cartesian point (u, v)∈Ω is denoted by Fu,v(h, p) , where the total radiation at any point (u, v) is the sum of contributions from all heliostats. The total radiation captured over the whole receiver surface associated to heliostat h aiming at point ph∈¯ Ω can be written as f(h, ph) = RFu,v(h, ph)dΩ . Therefore the total radiation corresponding to all heliostats in H can be expressed in the form: X h∈H f(h, ph). (4.2.1) The second criterion in this problem considers the dierence between the radiation reaching the receiver and a desired target distribution Etar u,v . This distribution will in practise be decided by the SPT plant operators, depending on weather conditions and the thermal status of the receiver. This second objective can be expressed as the integral of the square of the calculated radiation Fu,v(h, ph) minus the target distribution Etar u,v : Z"X h∈H Fu,v(h, ph)−Etar u,v #2 dΩ. (4.2.2) Other criteria, such as the overall radiation excess with respect to a target distribution, may be considered instead. The reason to choose (4.2.2) in this formulation is that, as seen below, it leads to a continuously dierentiable objective function, which allows us to consider gradient based algorithms for the problem resolution. Combining (4.2.1) and (4.2.2), we arrive at a suitable objective function, where we consider a parameter A∈(0,1) that controls the relative importance between maximising energy and minimising deviation from the target distribution: AX h∈H f(h, ph)−(1 −A)Z"X h∈H Fu,v(h, ph)−Etar u,v #2 dΩ. (4.2.3) which is to be maximised. Optimising the objective function in (4.2.3) over a range of values of A produces an approximation to the Pareto frontier of the bi-objective problem of simultaneous optimisation of energy generated and deviation with respect to the target distribution. The following section devises a numerical method for the solution of the continuous 62 Chapter 4. Continuous aiming strategies (I): the stationary case optimisation model proposed. 4.3 Numerical Methods The continuous optimisation model proposed in Section 4.2 involves calculating integrals over the receiver surface Ω . For numerical purposes, these integrals must be replaced by summations over a nite set of equally spaced test points (ui, vi)∈Ω with i= 1, ..., I . The size of I directly aects the precision to which the problem in Equation (4.2.3) is approximated and, also, aects the numerical complexity of the problem solution. The deviation from the desired ux distribution Etar u,v is then calculated and furnished only for the test points (ui, vi) . It is important to note that we `discretise' the receiver Ω only at the numerical integration level (remember that ph can take any value in Ω ). The resulting objective function is the following: g(P):=AX h∈H fsp(h, ph)−(1 −A)|Ω| I I X i=1 "X h∈H Fui,vi(h, ph)−Etar ui,vi#2 , (4.3.1) where P= (ph)h∈H , |Ω| denotes the measure of Ω and |I| is the cardinality of I . The coecient |Ω| I stems from the numerical approximation of the integral in (4.2.3). Thus, we want to solve a non-linear non-convex optimisation problem with continuous variables of large dimension (twice the number of heliostats in the eld), subject to the convex constraints ph∈¯ Ω . Function g in Equation (4.3.1) can then be maximised using a gradient ascent algorithm with projection, see Le Floch et al. [2015]; Ranganathan et al. [2011]. We start with an initial solution P0 and then, at each iteration of the algorithm, update the components of P in the direction of steepest ascent of the objective function g by step size γ . The nal step in each iteration of the algorithm utilises a projection method to correct any values of P to ensure heliostats aim at Ω . The selection of the stepsize γ taken at each update to the gradient ascent algorithm is an important factor in the convergence and much research has focused on this choice, for example Liu and Liu [2018]. If the stepsize is too large, the algorithm may diverge, and if it is too small, it will take too long to converge. A method to nd the optimal stepsize at each iteration can be found using Armijo's Rule [McCormick, 1977], where a constant value ∈(0,1) is used to iteratively reduce the stepsize until an improvement on the objective function is no longer given. This is tested at the k -th iteration of the algorithm against the (k−1) -th iteration as follows: g(Pk)> g(Pk−1) (4.3.2) 4.3. Numerical Methods 63 Traditionally, the value of γ is xed for all elements of the system. However, in this work the algorithm is adapted to allow dierent γ for each heliostat h : γk,h =γk−1,h ·. (4.3.3) This permits heliostats to take the greatest stepsize independently of each other, thereby potentially increasing the speed of the algorithm. However, it is also important to note that this procedure may in fact increase running times, so careful selection of  is required. Therefore, the gradient ascent algorithm reads: ˜ Pk+1 =Pk+γk,h∇g(Pk),Pk+1 =P(˜ Pk+1), (4.3.4) where, for each P , P(P) denotes the component-wise projection of P onto Ω . The considered objective function and its gradient are complex. It is therefore important to customise related eective numerical techniques leading to reasonable computational eort. The function being modelled in this chapter is highly multi-modal and, consequently, convergence to the global optimum is signicantly dependent on the starting P0 given to the algorithm. We therefore apply a multistart procedure where the algorithm is run multiple times, utilising a dierent random allocation of starting heliostat aiming points. Applying this method for each value of A , it is possible to approximate the Pareto Front of the model. It will be accepted that the aforementioned algorithm converges to a local optimum when the Euclidean norm of the gradient of the energy function is below a selected precision value c . The selection of this value determines how close to the local optimum the algorithm must nish, whilst also heavily inuencing the computation time required. Summarising, the convergence test used in the algorithm can be written in the form k ∇g(P)k≤ c, where ∇g(P) = AX h∈H∇fsp(h, ph) −2(1 −A)|Ω| I I X i=1 "(X h∈H Fui,vi(h, ph)−Etar ui,vi)X h∈H∇Fui,vi(h, ph)#. (4.3.5) Note that this approach can serve to consider many prescribed ux distributions Fu,v (not necessarily Gaussian). 64 Chapter 4. Continuous aiming strategies (I): the stationary case In the following section we apply the numerical discretisation of the continuous optimisation model developed in this chapter to a real SPT plant and we present an illustrative example. 4.4 Illustrative Example The developed optimisation procedure is illustrated using the PS10 SPT plant in Sanlúcar la Mayor, Seville [Abengoa, 2019], as detailed in Section 1.6.6. This SPT plant has a eld of 624 heliostats in a South facing eld, arranged as shown in Figure 4.1. The receiver Ω is a disk and the grid points (ui, vi) have been xed equally spaced and of equal number in both dimensions. Figure 4.1: PS10 Heliostat Locations The algorithm developed in Section 4.3 is implemented in Python on a standard specication desktop computer, for a chosen time point of midday. The Armijo's Rule parameter  is set to a value of 0.8, as the suggested value in McCormick [1977], and the initial value for γ is set to 0.01. A lower limit of γ equal to 10−8 is set, in order to prevent unnecessary computations occurring in the algorithm. As mentioned in Section 4.3, since the problem is highly multi-modal, a multistart procedure has been implemented, where the best solution across 30 runs is selected for each variation of the parameter A . Figure 4.2 shows the peak of the objective function for one value of A across 100 multistart runs, where the multi-modality of the problem can be clearly seen. The variance in solutions found over the 100 multistart runs is also aected by the stop criterion and step size used in each application of the algorithm, meaning that whilst only a few solutions may be found, the level of convergence may not be the same. This can be seen in Figure 4.2 by the number of peaks and troughs at similar, but 4.4. Illustrative Example 65 slightly dierent, values. These possibly represent the same solution, reached at dierent levels of convergence. Therefore, when utilising an adequate selection of step size and convergence test, it is found that a multistart operation with 30 runs suces in practise to nd a solution. The convergence of the gradient to zero and the objective function to the solution can be seen in Figures 4.3-4.4 for one particular simulation. Figure 4.2: Multistart analysis Figure 4.3: Gradient convergence Figure 4.4: Objective function convergence The parameter A in the objective function in Equation (4.3.1) has been tested between 0 and 1 in steps of 0.01. A value of 0 indicates that the deviation from the target distribution objective is the most important contribution. Contrarily, a value of 1 indicates that the total radiation objective is the most relevant. The target distribution has been assumed constant across the receiver surface. It is important to note that the chosen target distribution must be tailored to the SPT plants current conditions, and a constant distribution has been implemented here as an example. The values of both objectives for each value of A are shown in Figure 4.5, where 66 Chapter 4. Continuous aiming strategies (I): the stationary case the Pareto Front has also been marked. The set of Pareto equilibria identies those solutions which cannot be improved in terms of both objectives, and therefore give the best optima in terms of minimising target distribution deviation and maximising total radiation capture. Figure 4.5: Objective values with Pareto Front The results shown in Figure 4.5 indicate that the choice of A can produce highly diering results in the objective function, and this can generate quite dierent aiming strategies trading o energy maximisation and minimisation of deviation from the target radiation distribution. This is illustrated in Figures 4.6a-4.7b. Figure 4.6a shows the resultant aiming strategy when the value of A is set to 0.9, causing the captured radiation to be signicantly more important than adhering to the target distribution. The aiming strategy will capture more radiation, as shown in Figure 4.6b, but will fail to yield a homogeneous ux distribution, due to the centrally focused heliostats. Figure 4.7a shows the resultant aiming strategy when A is set to 0.3, which creates a much more homogeneous ux distribution, as shown in Figure 4.7b. These Figures demonstrate the importance of the value of making a judicious choice of A , since improvement in both objectives is highly dependent on its choice. The results illustrated in this section are for one particular time point, and two examples of the weighting variable A . Each run of this simulation takes less than 10 seconds, which is then multiplied by the number of multistarts performed. With advance knowledge of local weather conditions, the optimal aiming strategy can be calculated without conditions on time. However, to account for rapid changes in weather, it may be useful to re-calculate the optimal aiming strategy in short time scales. Therefore, the rapid computation of this algorithm is advantageous, and also lends itself to applications in SPT plants with larger heliostat eld sizes. 4.5. Conclusions 67 (a) Optimised aiming strategy (b) Optimised ux distribution with maximum ux 2.8e6 Figure 4.6: Optimised result with A= 0.9 (a) Optimised aiming strategy (b) Optimised ux distribution with maximum ux 2e6 Figure 4.7: Optimised result with A= 0.3 The next section draws conclusions from the method and numerical illustration developed in this work, and also discusses possible extensions and current research. 4.5 Conclusions In this work, a bi-objective optimisation model has been implemented to nd the optimal aiming strategy for an SPT plant of any size or shape, and a numerical illustration for a real SPT plant is presented. For the PS10 SPT plant with a eld of 624 heliostats, the optimal aiming strategy has been found using the objective function given in (4.3.1), that must be viewed as 68 Chapter 4. Continuous aiming strategies (I): the stationary case a numerical approximation to (4.2.3). In this illustrative example, we have solved the optimisation problem with a multistart procedure for multiple values of the parameter A in (0,1) . Utilising the optimisation model and numerical method developed in this work, it is possible for the operators of an SPT plant to identify the optimal aiming strategy, considering current weather and plant requirements. Depending on the choice of the weighting parameter A in the objective function, it is possible to seek a desired balance between maximising overall radiation captured and minimising deviation from a desired distribution across the receiver. The method developed in this work signicantly improves upon the method presented in Chapter 3, as here the problem space is not constrained by pre-set aiming points, and allows for realistic representation of SPT plant operations, where it is possible to monitor relative performance of each objective. The numerical method developed for this work made use of Armijo's rule to iterate the step-size in the algorithm, with independent step-sizes for all components in the set space. This modication can cause longer run times for the simulation. However, a careful application of the algorithm in this case allows it to improve performance and converge to the optimal solution. This performance increase over other methods is of critical importance when simulations in an SPT plant will be re-run during the course of a day for changing local weather conditions. For the PS10 SPT plant, the optimal aiming strategy can be easily found using the method outlined in this work. For larger SPT plants with more heliostats in the eld, the same arguments and techniques can be applied. However, if the number of heliostats is large, and the associated computational cost becomes untenable, an amendment can be made to reduce the problem dimensionality. A possible modication could be the use of a clustering algorithm, such as in Ashley et al. [2017]; Carrizosa et al. [2013]. In such research, the heliostats in the eld are clustered using an optimisation procedure which takes into account potential radiation generation as well as physical location and, then, the same aiming strategy is chosen for all heliostats within the same cluster. Another extension to this work that could assist in reducing computational cost is the application of stochastic techniques [Fonseca et al., 2017; Schmidt et al., 2013; Wang, 2017] where, at each step, the gradient is calculated for just a random sample of heliostats and then used to update the general population. Such methods reduce the calculation time required for each iteration of the algorithm at the cost of utilising several (or many) potentially erroneous components of the gradients. 5.3. Some Iterative Algorithms 75 The set Pad is closed in P , thanks to (5.1.5) and (5.1.7). Hence, ˆ p∈Pad . On the other hand, p7→ J(p) is continuous in C0([0, T]; F) , thanks to (5.1.2). Therefore, one has J(ˆ p) = max p∈Pad J(p). This ends the proof.  Remark 1. In view of the possible non-convexity of J and M , uniqueness is out of scope in general.  Recall that P0 is a non-empty closed convex set of H1(0, T;F) . In our second result, we present suitable necessary optimality conditions that must be satised by the solutions to (5.1.1). For simplicity, we will assume from now on that E=R2 and we will denote by M1 and M2 (resp. σ1 and σ2 ) the components of M (resp. σ ). Theorem 2. Let ˆ p be a solution to (5.1.1) . Assume that the constraints associated to M are qualied at ˆp , that is:  ∃q1,q2∈H1(0, T;F) such that hM0 i(ˆ p),qii<0 for i= 1,2 Then, there exist λ1, λ2≥0 such that the triplet (ˆ p, λ1, λ2) satises    hJ0(ˆ p),p−ˆ pi−λ1hM0 1(ˆ p),p−ˆ pi−λ2hM0 2(ˆ p),p−ˆ pi ≤ 0∀p∈P0, λi(Mi(ˆ p)−σi)=0, i = 1,2. (5.2.1) For the proof, it suces to apply directly the Karush-Kuhn-Tucker principle; see for instance Theorem 9-2.3 in Nocedal et al. [1999]. 5.3 Some Iterative Algorithms 5.3.1 Penalisation In this section, we introduce an iterative method for the solution of the dynamic optimisation problem (5.1.1) based on penalisation techniques. The advantage of this approach is that it reduces the task to the solution of another optimisation problem whose constraints are very easy to handle. The drawback is that a (small) parameter must be introduced and this can have a signicant (undesired) inuence in the results. Thus, let us x µ > 0 and let us set Jµ(p) := J(p)−1 2µ|(M(p)−σ)+|2, (5.3.1) Here, z+= (M(p)−σ)+ stands for the positive part of z (understood in the component-wise sense) and, as before, |·| denotes the Euclidean norm in E . 76 Chapter 5. Continuous aiming strategies (II): the dynamic case Then, the aforementioned approximation of (5.1.1) is the following:    Maximise Jµ(p) Subject to p∈P0. (5.3.2) It is reasonable to expect that, for small µ > 0 , any solution to the penalized problem (5.3.2) solves approximately (5.1.1). In the remainder of this section, we will be concerned with the numerical solution of (5.3.2). To this purpose, we will rst introduce a time discretisation and then an iterative algorithm of the gradient ascent kind. Thus, let us begin by replacing H1(0, T;F) by a nite dimensional subspace. The easiest and most natural way is to introduce a large integer N , set τ:= T/N , consider a uniform partition {t0= 0 < t1< t2< ... < tN=T} with tn=nτ for all n and then work in the corresponding space of functions pN: [0, T]7→ F which are continuous and piecewise linear. This space will be denoted by XN . We will also consider the closed convex set P0,N :=XN∩P0 and the orthogonal projector P0,N :XN7→ P0,N . Observe that P0,N (pN)(t) is, for each t , the projection of each component of pN(t) onto Ω . The N -th approximated problem is then    Maximise Jµ(pN) Subject to pN∈P0,N . (5.3.3) In order to solve this problem, we can apply a gradient ascent algorithm with variable step size and projection. Accordingly, in the n -th iterate we compute the function pn+1 N , with pn+1 N=P0,N (˜ pn+1 N),˜ pn+1 N=pn N+γn∇Jµ(pn N). (5.3.4) Here, γn is a conveniently chosen positive number and ∇Jµ(pn N) denotes the gradient of the objective function in (5.3.3), that is, ∇Jµ(pN) = ∇J(pN)−1 µ(M(pN)−σ)+·∇M(pN). (5.3.5) Therefore, the gradient ascent algorithm requires, at each step, the calculation of the derivative of the objective function and a projection of each component of pN(t) for each nodal time t=tj . Obviously, the complexity of this computation depends on the properties of the particular function G=G(t, p) and the mapping M:H1(0, T;F)7→ E and the sizes of dim F and N . The described penalisation algorithm, denoted Algorithm 1 in this chapter, is outlined in the pseudocode shown below. 5.3. Some Iterative Algorithms 77 Algorithm 1 Penalisation µ← Set penalisation parameter (small) oat γ0← Set initial step size oat Tol ← Simulation tolerance oat p0← Initial variable set oat for N = 1 to T do while ObjDi > Tol do while StepDi ≤0 do J(p)← Calculate objective function oat M(p)← Calculate constraint set oat ∇Jµ(pN) = ∇J(pN)−1 µ(M(pN)−σ)+·∇M(pN). ˜ pn+1 N=pn N+γn∇Jµ(pn N). pn+1 N=P0,N (˜ pn+1 N) J(˜ pn+1 N)← Calculate updated objective function oat StepDi =|Jµ(pn+1 N)−Jµ(pn N)| γn+1 =γn· Next ; end while ObjDi =|Jµ(pn+1 N)−Jµ(pn+1 N−1)| Next ; end while end for end ; 5.3.2 Augmented Lagrangian The optimisation problem (5.1.1) can also be solved using the information furnished by Theorem 2. This is indicated in this section, again after time discretisation, through the so called Augmented Lagrangian techniques. Recall that E=R2 and M1 and M2 denote the components of M . Thus, let us rst introduce the so called Augmented Lagrangian Lµ:H1(0, T;F)× E7→ R , with Lµ(p;λ) := J(p)− 2 X i=1 ψ(Mi(p)−σi, λi;µ), (5.3.6) where we have introduced ψ(z, β;µ) :=      z·β+1 2µ|z|2 if z+µβ ≥0 −µ 2|β|2 otherwise, (5.3.7) 78 Chapter 5. Continuous aiming strategies (II): the dynamic case and, again, µ > 0 . Then, it can be proved that the optimisation problem (5.1.1) is equivalent to the following      Minimise sup p∈P0Lµ(p;λ) Subject to: λ∈E, λ ≥0. (5.3.8) An explanation of this equivalence can be given arguing as follows. The inequality constraint M(p)≤σ is equivalent to the equality M(p)+s=σ , with the so called slack variable s belonging to E and s≥0 . In view of the Karush-Kuhn-Tucker optimality condition (5.2.1), it makes sense to consider the Modied Lagrangian ˆ L(p, s;λ) := J(p)−λ·(M(p) + s−σ) (5.3.9) and its penalised version ˆ Lµ(p, s;λ) := J(p)−λ·(M(p) + s−σ)−1 2µ|M(p) + s−σ|2. (5.3.10) Then, recalling (5.3.6)-(5.3.7), it is easy to check that sup s∈E, s≥0 ˆ Lµ(p, s;λ) = Lµ(p;λ), whence we see that (5.3.8) is an appropriate reformulation of (5.1.1). As in the previous section, in practice, in order to solve (5.3.8) we must provide a nite dimensional approximation. With the notation used in Section 5.3.1, a suitable choice is the following:      Minimise sup P0,N Lµ(pN;λ) Subject to: λ∈E, λ ≥0. (5.3.11) This problem can be solved with a duality-penalty algorithm that, at the n -th step, furnishes the multiplier λn+1 according to the following:  Compute a solution pn N to the problem    Maximise Lµ(pN;λn) Subject to pN∈P0,N (5.3.12)  Then, take λn+1 = (λn+1 µ(Mi(pn N)−σ))+. (5.3.12 0 ) This will be denoted Algorithm 2 in this work, and is given in the pseudocode given 5.3. Some Iterative Algorithms 79 below. As in Section 5.3.1, the solution of (5.3.12) can be obtained through a variable step gradient ascent method with projection. Thus, at the k -th step, we compute pn,k+1 N , with pn,k+1 N=P0,N (˜ pn,k+1 N),˜ pn,k+1 N=pn,k N+γk∇Lµ(pn,k N;λn). (5.3.13) Algorithm 2 Augmented Lagrangian λ← Lagrangian multiplier oat σ← Constraint constant oat µ← Penalty constant oat Tol ← Simulation tolerance oat p← Variable set oat J ← Objective function oat M ← Constraint set oat for t = 1 to T do while ObjDi > Tol do while StepDi ≤0 do J(p)← Calculate objective function oat M(p)← Calculate constraint set oat ∇Jµ(pN) = ∇J(pN)−1 µ(M(pN)−σ)+·∇M(pN). Maximise ˆ Lµ(p, s;λ) Fix s Maximise Lµ(pN;λn) ˜ pn,k+1 N=pn,k N+γk∇Lµ(pn,k N;λn) pn,k+1 N=P0,N (˜ pn,k+1 N) Jµ(˜ pn+1 N)← Calculate updated objective function oat StepDi =|Jµ(pn+1 N)−Jµ(pn N)| γk+1 =γk· Next ; end while ObjDi =|ˆ Lµ(pn,k+1 N;λn)−ˆ Lµ(pn,k N;λn)| λn+1 = (λn+1 µ(Mi(pn N)−σ))+ Next ; end while end for end ; In the following section, we detail the model for an SPT plant and verify that all assumptions made in the formulation of the general problem in Section 5.1 remain valid. 80 Chapter 5. Continuous aiming strategies (II): the dynamic case 5.4 The Model for an SPT Plant As considered in Chapters 3-4, the aiming point of each heliostat in the eld on the receiver surface directly aects the energy generation. In particular, it has been shown that, although maximising incident radiation is important, it is also benecial to maintain a desired ux distribution to aid thermal transfer [Relloso and García, 2015; Yu et al., 2014]. In this chapter we extend the models previously discussed to the dynamic case, considering time dependent variables and constraints. These constraints, as introduced in Sections 5.1-5.2, concern physical aspects of the SPT. One of them must be viewed as a limitation on the rotational speed of the heliostat and is obviously justied by operational reasons. On the other hand, the radiation at any point of the receiver surface should not vary drastically over short time periods, in order to prevent ash heating and this is considered in the second constraint. The radiation passing through the system is modelled using a Gaussian distribution, as presented in Section 1.7. We assume that there are H heliostats in the eld and we denote by ph(t) the point aimed at by the h -th heliostat at time t . We have ph(t)∈¯ ΩR for all h and t ; ΩR⊂R2 is bounded, open and convex with 0∈ΩR and we take Ω = ¯ ΩR . Accordingly, the Euclidean space F will have dimension 2H : F=R2H. (5.4.1) The usual Euclidean norm in F will be denoted by k·k . At time t , the radiation measured at a point (u, v)∈Ω and furnished by the h -th heliostat is given by Fu,v(h, ph(t), t) . For each h , the real-valued function (u, v, ph, t)7→ Fu,v(h, ph, t) is smooth. The associated total radiation captured over the receiver surface Ω for a given heliostat h is thus given by F0(h, ph(t), t) := ZZΩ Fu,v(h, ph(t), t)du dv. Therefore, the total radiation supplied by all heliostats at time t can be written in the form f(p(t), t) := H X h=1 F0(h, ph(t), t) = ZZΩH X h=1 Fu,v(h, ph(t), t)du dv, (5.4.2) where we have used the notation p:= (p1, . . . , pH) . In order to limit the motion of heliostats, so that they do not move faster than their velocity limits, we introduce the velocity ˙ph= ˙ph(t) of each heliostat h , the velocity vector ˙ p=˙ p(t) and a target velocity vector Vp∈F with Vp≥0 and we impose the 5.4. The Model for an SPT Plant 81 collective velocities to be approximately below Vp in the sense that ZT 0k(˙ p(t)−Vp)+k2dt ≤σ1 for some σ1≥0 . The change in time in radiation at each point of the receiver must also be limited and kept below a xed value σ2 . This will be witten in the form ZT 0 ZZΩ H X h=1 d dtFu,v(h, ph(t), t) 2du dv!dt ≤σ2. Our objective is to maximise the radiation captured by the receiver, whilst maintaining a target energy prole which maximises absorption. Thus, we deal with the constrained optimisation problem (5.1.1), where              G(t, p) := AZZΩH X h=1 Fu,v(h, ph, t)du dv −(1 −A)ZZΩ H X h=1 Fu,v(h, ph, t)−Etar u,v (t) 2du dv (5.4.3) ( A∈[0,1] is a weighting parameter), M1(p):=ZT 0k(˙ p(t)−Vp)+k2dt, (5.4.4) M2(p):=ZT 0ZZΩ H X h=1 d dtFu,v(h, ph(t), t) 2du dvdt, (5.4.5) the σi are prescribed non-negative constants and Etar u,v =Etar u,v (t) is the desired target distribution at time t . For convenience, we will assume that σ2> M2(0) = ZT 0ZZΩ H X h=1 d dtFu,v(h, 0, t) 2du dvdt. (5.4.6) This is sucient to guarantee that the zero function in H1(0, T;F) belongs to Pad and, consequently, Pad is non-empty. In the following sections, we will check that the theoretical analysis in Section 5.2 and the proposed iterative algorithms in Section 5.3 are valid in this framework. In particular, we consider the existence of optimal aiming strategies, their characterisation through optimality conditions and nite dimensional approximations and penalty and duality-penalty solution methods are described. 82 Chapter 5. Continuous aiming strategies (II): the dynamic case 5.4.1 Existence and Optimality The following result holds: Proposition 1. Let F , G and M be given by (5.4.1) and (5.4.3)  (5.4.5) and let P and Pad be given by (5.1.3)  (5.1.4) . Then, the assumptions (5.1.2) and (5.1.5)  (5.1.8) are satised. Proof: It is clear that (5.1.2), (5.1.5) and (5.1.6) are satised. For convenience, let us introduce ah(t, ph;u, v):=∂ ∂phFu,v(h, ph, t) , bh(t, p;u, v):=PH h=1 Fu,v(h, ph, t)−Etar u,v (t) and c(t, p,q;u, v):=PH h=1 d dt ∂ ∂phFu,v(h, ph, t)qh . Then one has ∂G ∂ph (t, p) = AZZΩ ah(t, ph;u, v)du dv −2(1 −A)ZZΩ bh(t, p;u, v)ah(t, ph;u, v)du dv, for all h and (t, p) and hM0 1(p),qi= 2 ZT 0 (˙ p(t)−Vp)+·˙ q(t)dt (5.4.7) and hM0 2(p),qi= 2 ZT 0hZZΩH X h=1 d dtFu,v(h, ph(t), t)c(t, p(t),q(t); u, v)du dvidt (5.4.8) for all p,q∈H1(0, T;F) . Let us now check that (5.1.7) holds. Assume that pn→p weakly in H1(0, T;F) . Since this implies the uniform convergence of the pn in [0, T] , we have lim n→+∞M2(pn) = M2(p) (5.4.9) On the other hand, since the function p7→ M1(p) is convex and continuous in H1(0, T;F) , we also have lim inf n→+∞M2(pn)≥M2(p). (5.4.10) From (5.4.9) and (5.4.10), we get (5.1.7). Finally, it is clear that (5.1.8) is satised. Indeed, if p∈P0 and M(p)≤σ , we have in particular that kp(t)k ≤ C∀t∈[0, T] (5.4.11) 5.4. The Model for an SPT Plant 83 and ZT 0||˙ p(t)||2dt ≤2ZT 0||(˙ p(t)2−Vp)+||2dt +C≤C, (5.4.12) whence p belongs to a bounded set in H1(0, T;F) . This ends the proof.  From Proposition 1, we get the following consequences:  Theorem 1 can be applied and there exists at least one optimal aiming strategy for the modelled SPT plant.  Theorem 2 can also be applied and, assuming that at an optimal ˆ p the constraints associated to the Mi are qualied, we get the necessary conditions (5.2.1) for some multipliers λ1, λ2≥0 . 5.4.2 Finite Dimensional Approximation and Iterative Algorithms In practice, as in Sections 5.3.1 and 5.3.2, we must approximate the innite dimensional problem (5.1.1) corresponding to (5.4.3)-(5.4.5) and replace H1(0, T;F) by XN and P0 by P0,N . The resulting tasks are thus to solve (5.3.3) and/or (5.3.11). Note that, in both cases, the computations of ∇J(pN) , ∇M1(pN) and ∇M2(pN) are needed and this requires integrals on Ω , with respect to (u, v) of several functions. For this reason, it is convenient to x a set of test points (ui, vi) , 1≤i≤I and replace these integrals by appropriate nite sums. The total radiation reaching the receiver at time t is therefore approximated by: ZT 0ZZΩH X h=1 Fu,v(h, ph(t), t)du dv dt ≈ N X n=1 I X i=1 H X h=1 Fui,vi(h, ph(tn), tn) (5.4.13) Similarly, the deviation from the target distribution is approximated by: ZT 0ZZΩ b(t, p(t); u, v)2du dv dt ≈ N X n=1 I X i=1 b(tn,p(tn); ui, vi)2. (5.4.14) Finally, the constraint mappings M1 and M2 given by (5.4.4)-(5.4.5) are approximated as follows: M1(p)=ZT 0k(˙ p(t)−Vp)+k2dt ≈MN 1(pN):= N X n=1 (1 2|p(tn)−p(tn−1)|−Vp)2 +, (5.4.15) 84 Chapter 5. Continuous aiming strategies (II): the dynamic case M2(p)=ZT 0 ZZΩH X h=1 d dtFu,v(h, pn h(t), t)2du dv!dt ≈MN 2(pN):= N X n=1 I X i=1  H X h=1 Fui,vi(h, ph(tn), tn)−Fui,vi(h, ph(tn−1), tn−1) 2. (5.4.16) Penalisation Algorithm After these approximations, we see that, in the case of the SPT model, the objective function in (5.3.3) is given by: JN µ(pN) = A N X n=1 I X i=1 H X h=1 Fui,vi(h, ph(tn), tn)−(1 −A) N X n=1 I X i=1 b(tn,p(tn); ui, vi)2 −1 2µh N X n=1 (1 2|p(tn)−p(tn−1)|−Vp)2 +−σ12 + +N X n=1 I X i=1  H X h=1 Fui,vi(h, ph(tn), tn)−Fui,vi(h, ph(tn−1), tn−1) 2−σ22 +i, (5.4.17) It is not dicult to compute from this identity the partial derivatives of JN µ and, accordingly, its gradient ∇Jµ . This makes it possible to apply the gradient ascent algorithm in the present context. The choice of step size parameters γk used for the gradient ascent algorithm have a large impact on the speed of convergence. In the SPT plant context, the objective function is non-linear, non-convex and of large dimensions, indicating that step size choice will also factor heavily into likelihood of convergence of the algorithm. As in Chapter 4, we apply a variant to the Armijo's Rule , whereby the parameter γk is dierent for each component, that is, for each heliostat h∈H . Consequently, the γk h are computed according to the rules γk h=γk−1 h·h, h = 1, . . . , H. (5.4.18) The gradient ascent method with projection is then: pn+1 N=P0,N (˜ pn+1 N),˜ pn+1 N=pn N+ Γk∇Jµ(pn N), (5.4.19) where Γk= diag (γk 1, . . . , γk N). The gradient ascent method applied to the penalisation algorithm, as described in the pseudocode given in Algorithm 1, is applied at each step of the discretised time period. The computed optimal result at a given time is then used as the initial guess 5.5. A Numerical Experiment 91 This is repeated until no improvement on the solution can be found. The resultant solution is then used as the initial solution to the next time point. The computed aiming strategies using the augmented lagrangian algorithm are then given in Figures 5.17-5.41. Figure 5.16: t=0 aiming distribution Figure 5.17: t=1, λ0 aiming distribution Figure 5.18: t=1, λ1 aiming distribution Figure 5.19: t=1, λ2 aiming distribution Figure 5.20: t=1, λ3 aiming distribution 92 Chapter 5. Continuous aiming strategies (II): the dynamic case Figure 5.21: t=2, λ0 aiming distribution Figure 5.22: t=2, λ1 aiming distribution Figure 5.23: t=2, λ2 aiming distribution Figure 5.24: t=3, λ0 aiming distribution Figure 5.25: t=3, λ1 aiming distribution 5.5. A Numerical Experiment 93 Figure 5.26: t=3, λ2 aiming distribution Figure 5.27: t=4, λ0 aiming distribution Figure 5.28: t=4, λ1 aiming distribution Figure 5.29: t=5, λ0 aiming distribution Figure 5.30: t=5, λ1 aiming distribution 94 Chapter 5. Continuous aiming strategies (II): the dynamic case Figure 5.31: t=5, λ2 aiming distribution Figure 5.32: t=6, λ0 aiming distribution Figure 5.33: t=6, λ1 aiming distribution Figure 5.34: t=7, λ0 aiming distribution Figure 5.35: t=7, λ1 aiming distribution 5.5. A Numerical Experiment 95 Figure 5.36: t=8, λ0 aiming distribution Figure 5.37: t=8, λ1 aiming distribution Figure 5.38: t=9, λ0 aiming distribution Figure 5.39: t=9, λ1 aiming distribution Figure 5.40: t=10, λ0 aiming distribution Figure 5.41: t=10, λ1 aiming distribution The maximum ux on the receiver surface and the incident radiation level on the heliostat eld are given in Figure 5.42. 96 Chapter 5. Continuous aiming strategies (II): the dynamic case Figure 5.42: Maximum ux on receiver (black) and incident solar radiation (red) The results furnished by the Augmented Lagrangian algorithm agree with those provided by the penalisation algorithm. 5.5.3 Numerical Considerations The PS10 plant eld contains 624 heliostats, however there exist other larger SPT plants (with potentially thousands of heliostats) which greatly increases the dimensionality of the problem. The aiming points are described on the receiver surface in polar coordinates, necessitating 2H variables for the gradient ascent technique. Within each step of the algorithm, the gradients must be calculated, and aiming points updated, before the radiation calculations are performed. For larger values of H , the application of the algorithms presented in this chapter may become infeasible, due to computation times. For the numerical approximations presented in this work, the heliostat eld size will have an eect on performance. 5.6 Conclusions In this chapter, a general dynamic optimisation problem has been investigated. Theoretical properties of the problem have been discussed, including the demonstration of solution existence and optimality. Furthermore, a real-world example for an SPT plant has been considered, where a mathematical model that describes transfer of energy within the system is developed. The theoretical properties of the general optimisation problem are shown applicable to the real world problem, and dynamic constraints that describe SPT plant limitations are given. The goal of this work was the dynamic optimisation of aiming strategies in an SPT plant, where the objective function was non-linear, non-convex and of large dimension. Two algorithms have been considered to nd the optimal solution, and a numerical illustration is given using real data. 5.6. Conclusions 97 The numerical illustration nds the optimal aiming strategy for an SPT plant over the period of one day, considering the change in incident radiation from the Sun as input. The physical limitations of the plant have been introduced as dynamic constraints, in terms of heliostat rotation speed and ux homogeneity on the receiver surface over time. The two algorithms presented (Penalisation and Augmented Lagrangian) provide similar results in the numerical experiment, in similar computational times. The Augmented Lagrangian algorithm is a modication of a penalty technique, which should provide better numerical stability in some cases. With larger size problems, and adequately chosen parameters, this model could increase the reliability of the algorithm. The dynamic optimisation problem considered in this work must be viewed as an improvement of other research concerning the optimisation of aiming strategies in SPT plants, for example Astol et al. [2017b]; Wang [2017], due to the inclusion of dynamic constraints. Instead of optimising the aiming strategy at certain xed times, the method presented in this work looks to optimise across a time period. This approach can arrive to a solution that more closely reects the true optimum when considering problems in dynamic systems. The methods presented in this chapter can be adapted to all types of SPT plants, and even other forms of concentrating solar power technology. The inclusion of more heliostats, or multiple receivers, has been carried out in several real-world plants and the work presented here can be directly extended to consider these cases. The use of a modied Armijos' Rule for the step size in the gradient ascent method can lead to faster convergence. However, with a poor parameter selection, this may actually reduce the convergence speed. Further adaptation to the algorithms presented here should consider carefully the eects of such techniques, in conjunction with a highly multi-modal objective function. Considering the problem dimensionality when increasing the number of variables, specically the number of heliostats in our example, one method that could be used to extend the presented approach could be to use a clustering algorithm. Thus, the heliostats could be clustered into groups, considering a weighted objective of dierence in location and energy, as detailed in Chapter 2. As remarked in Section 5.5.3, as the dimensionality of the problem increases, the performance of the proposed algorithms deteriorates. An extension to this work could consider a stochastic gradient method, which would reduce the required calculations at each step of the algorithm, thereby allowing larger problems to be tackled with the proposed method. Finally, the integration of dynamic aiming strategy optimisation with a three dimensional thermal transfer model of the receiver is considered as the next step in this work. Thermal transfer from the incident radiation on the receiver surface, through to the HTF in the interior, is a three dimensional dynamic process that is highly depen- 98 Chapter 5. Continuous aiming strategies (II): the dynamic case dent on both solar input and current temperature distribution in the receiver materials. Varying factors such as the aiming strategy and the ow rate of the HTF will aect this process, and should be included in the dynamic coupled model. Chapter 6 Inclement Weather Eects on Optimal Aiming Strategies in Solar Power Tower plants 100 6.4. Conclusions 107 The optimal aiming strategy for the heliostat eld changes with a cloud implemented, but does not exhibit a specic pattern based upon the cloud location. This is expected, as the allocation of aiming points is primarily aected by the cosine angle of the incident radiation, and less by the magnitude of radiation. Therefore the results show that to maintain the constraints of homogeneity on the receiver surface, the overall aiming strategy must change, not only the heliostats covered by the cloud. 6.4 Conclusions The eect of uncertainty in cloud location, size and density on the optimal aiming strategy for an SPT plant heliostat eld has been investigated. The method has been demonstrated for a time point at the PS10 SPT plant and shows how the eciency of the plant can be maximised when there is local inclement weather. A stochastic linear integer programming technique has been applied with short time limits, which provide near-optimal solutions in near real-time. It is expected that increasing computational power and providing the solver with an initial heuristic solution could further increase the speed of the program. The uncertainty of the cloud parameters has been modelled in this chapter using a uniform distribution, considering 5 possible scenarios. This method is directly extendible to include more scenario possibilities and a dierent probability distribution for their occurrence. These factors can be determined by considering historical weather data for the region of interest, as well as the inclusion of weather predicting technology such as satellite data and cloud detecting cameras. The optimisation method used in this research intends to nd a near-optimal solution, that is solved within short timescales. It was found that solutions suciently close to the optimal solution can be found with a simulation time of 30s, allowing this method to be implemented in real time with local weather predictions for an SPT plant. Depending on the location of the SPT plant, other types of inclement weather may be typical, such as rain, snow and sand storms. The method implemented in this article may be extended to account for such weather conditions, and demonstrate their eect on the eciency of an SPT plant. The next step in this work is the inclusion of stochastic inclement weather in the dynamic aiming strategy optimisation presented in Chapter 5. In this case, not only will the objective of maximising energy (or reaching a target distribution) be reached when considering inclement weather, but additional dynamic constraints can be added to prevent ash heating on the receiver due to cloud passage over the heliostat eld. List of Figures 1.1 PS10 Heliostat Locations . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.2 Average reected radiation per heliostat . . . . . . . . . . . . . . . . . . 18 1.3 Design point radiation prole . . . . . . . . . . . . . . . . . . . . . . . . 19 1.4 Standard heliostat design [Abengoa, 2019] . . . . . . . . . . . . . . . . . 21 1.5 DirectionalVectors.............................. 23 1.6 Cosineangles ................................. 23 1.7 Linear degradation function . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.1 Clustering Pseudocode . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.2 Average reected radiation per heliostat . . . . . . . . . . . . . . . . . . 34 2.3 PS10 Heliostat eld layout . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.4 Average reected radiation . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.5 Optimised Heliostat Grouping (groups 1-4) . . . . . . . . . . . . . . . . 41 2.6 Optimal Cleaning Schedule Allocation . . . . . . . . . . . . . . . . . . . 41 2.7 Optimal Cleaning Schedule Energy Prole . . . . . . . . . . . . . . . . . 42 2.8 Optimised Cleaning Schedule . . . . . . . . . . . . . . . . . . . . . . . . 43 2.9 Calculated Route Period 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.1 AimPointGrid................................ 49 3.2 PS10Layout ................................. 53 3.3 AimingPoints................................. 53 3.4 6amAllocation ................................ 55 3.5 6amDistribution............................... 55 3.6 12pmAllocation ............................... 55 3.7 12pmDistribution .............................. 55 3.8 6pmAllocation................................ 55 3.9 6pmDistribution............................... 55 3.10Cloudexample ................................ 56 3.1112pmAllocation ............................... 56 3.12 12pm Allocation with cloud . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.13 12pm Allocation (Greedy eld) . . . . . . . . . . . . . . . . . . . . . . 56 108 List of Figures 109 3.14 12pm Distribution (Greedy eld) . . . . . . . . . . . . . . . . . . . . . . 56 4.1 PS10 Heliostat Locations . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.2 Multistartanalysis .............................. 65 4.3 Gradientconvergence............................. 65 4.4 Objective function convergence . . . . . . . . . . . . . . . . . . . . . . . 65 4.5 Objective values with Pareto Front ..................... 66 4.6 Optimised result with A= 0.9 ....................... 67 4.7 Optimised result with A= 0.3 ....................... 67 5.1 HeliostatLayout ............................... 86 5.2 Reected radiation for midday . . . . . . . . . . . . . . . . . . . . . . . . 86 5.3 Initial aiming point distribution . . . . . . . . . . . . . . . . . . . . . . . 87 5.4 t=1 aiming distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.5 t=2 aiming distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.6 t=3 aiming distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.7 t=4 aiming distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.8 t=5 aiming distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.9 t=6 aiming distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.10 t=7 aiming distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.11 t=8 aiming distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.12 t=9 aiming distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.13 t=10 aiming distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.14 Incident radiation on SPT plant eld . . . . . . . . . . . . . . . . . . . . 90 5.15 Maximum ux on receiver over time . . . . . . . . . . . . . . . . . . . . 90 5.16 t=0 aiming distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.17 t=1, λ0 aimingdistribution ......................... 91 5.18 t=1, λ1 aimingdistribution ......................... 91 5.19 t=1, λ2 aimingdistribution ......................... 91 5.20 t=1, λ3 aimingdistribution ......................... 91 5.21 t=2, λ0 aimingdistribution ......................... 92 5.22 t=2, λ1 aimingdistribution ......................... 92 5.23 t=2, λ2 aimingdistribution ......................... 92 5.24 t=3, λ0 aimingdistribution ......................... 92 5.25 t=3, λ1 aimingdistribution ......................... 92 5.26 t=3, λ2 aimingdistribution ......................... 93 5.27 t=4, λ0 aimingdistribution ......................... 93 5.28 t=4, λ1 aimingdistribution ......................... 93 5.29 t=5, λ0 aimingdistribution ......................... 93 5.30 t=5, λ1 aimingdistribution ......................... 93 110 List of Figures 5.31 t=5, λ2 aimingdistribution ......................... 94 5.32 t=6, λ0 aimingdistribution ......................... 94 5.33 t=6, λ1 aimingdistribution ......................... 94 5.34 t=7, λ0 aimingdistribution ......................... 94 5.35 t=7, λ1 aimingdistribution ......................... 94 5.36 t=8, λ0 aimingdistribution ......................... 95 5.37 t=8, λ1 aimingdistribution ......................... 95 5.38 t=9, λ0 aimingdistribution ......................... 95 5.39 t=9, λ1 aimingdistribution ......................... 95 5.40 t=10, λ0 aimingdistribution......................... 95 5.41 t=10, λ1 aimingdistribution......................... 95 5.42 Maximum ux on receiver (black) and incident solar radiation (red) . . . 96 6.1 Aiming point layout on receiver . . . . . . . . . . . . . . . . . . . . . . . 105 6.2 PS10 Heliostat eld layout . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.3 Example of cloud implementation . . . . . . . . . . . . . . . . . . . . . . 106 6.4 Optimal aiming strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 6.5 Energy distribution on receiver surface . . . . . . . . . . . . . . . . . . . 106 6.6 Optimal aiming strategy with cloud uncertainty . . . . . . . . . . . . . . 106 6.7 Energy distribution on receiver surface with cloud uncertainty . . . . . . 106 List of Tables 1.1 PS10data................................... 18 2.1 Cleaning optimisation parameter values . . . . . . . . . . . . . . . . . . 40 5.1 Penalisation algorithm parameter values . . . . . . . . . . . . . . . . . . 87 5.2 Augmented Lagrangian algorithm parameter values . . . . . . . . . . . . 90 Bibliography Glenn A. Jones and Kevin J. Warner. The 21st century population-energy-climate nexus. Energy Policy, 93:206212, 2016. ISSN 03014215. doi: 10.1016/j.enpol.2016. 02.044. URL http://dx.doi.org/10.1016/j.enpol.2016.02.044 . Steve Sorrell. Reducing energy demand: A review of issues, challenges and approaches. Renewable and Sustainable Energy Reviews, 47:7482, 2015. ISSN 18790690. doi: 10. 1016/j.rser.2015.03.002. URL http://dx.doi.org/10.1016/j.rser.2015.03.002 . Florinda Martins, Carlos Felgueiras, and Miroslava Smitková. Fossil fuel energy consumption in European countries. Energy Procedia, 153:107111, 2018. ISSN 18766102. doi: 10.1016/j.egypro.2018.10.050. URL https://doi.org/10.1016/j. egypro.2018.10.050 . Steve Sorrell, Jamie Speirs, Roger Bentley, Adam Brandt, and Richard Miller. Global oil depletion: A review of the evidence. Energy Policy, 38(9):52905295, 2010. ISSN 03014215. doi: 10.1016/j.enpol.2010.04.046. URL http://dx.doi.org/10.1016/j. enpol.2010.04.046 . Ian Chapman. The end of Peak Oil? Why this topic is still relevant despite recent denials. Energy Policy, 64:93101, 2014. ISSN 03014215. doi: 10.1016/j.enpol.2013. 05.010. URL http://dx.doi.org/10.1016/j.enpol.2013.05.010 . Daniel J. Soeder. The successful development of gas and oil resources from shales in North America. Journal of Petroleum Science and Engineering, 163(December 2017): 399420, 2018. ISSN 09204105. doi: 10.1016/j.petrol.2017.12.084. URL https: //doi.org/10.1016/j.petrol.2017.12.084 . Dieter Helm. The future of fossil fuels-is it the end? Oxford Review of Economic Policy, 32(2):191205, 2016. ISSN 14602121. doi: 10.1093/oxrep/grw015. Richard G. Williams, Vassil Roussenov, Philip Goodwin, Laure Resplandy, and Laurent Bopp. Sensitivity of global warming to carbon emissions: Eects of heat and carbon uptake in a suite of earth system models. Journal of Climate, 30(23):93439363, 2017. ISSN 08948755. doi: 10.1175/JCLI-D-16-0468.1. 113 114 Bibliography David C. Mowery, Richard R. Nelson, and Ben R. Martin. Technology policy and global warming: Why new policy models are needed (or why putting new wine in old bottles won't work). Research Policy, 39(8):10111023, 2010. ISSN 00487333. doi: 10.1016/j.respol.2010.05.008. Payam Nejat, Fatemeh Jomehzadeh, Mohammad Mahdi Taheri, Mohammad Gohari, and Muhd Zaimi Muhd. A global review of energy consumption, CO2 emissions and policy in the residential sector (with an overview of the top ten CO2 emitting countries). Renewable and Sustainable Energy Reviews, 43:843862, 2015. ISSN 13640321. doi: 10.1016/j.rser.2014.11.066. URL http://dx.doi.org/10.1016/j. rser.2014.11.066 . Peter Christo. The promissory note: COP 21 and the Paris Climate Agreement. Environmental Politics, 25(5):765787, 2016. ISSN 17438934. doi: 10.1080/09644016. 2016.1191818. URL http://dx.doi.org/10.1080/09644016.2016.1191818 . Kevin J. Warner and Glenn A. Jones. A population-induced renewable energy timeline in nine world regions. Energy Policy, 101(August 2016):6576, 2017. ISSN 03014215. doi: 10.1016/j.enpol.2016.11.031. URL http://dx.doi.org/10.1016/j. enpol.2016.11.031 . Patrick Moriarty and Damon Honnery. Can renewable energy power the future? Energy Policy, 93:37, 2016. ISSN 03014215. doi: 10.1016/j.enpol.2016.02.051. URL http: //dx.doi.org/10.1016/j.enpol.2016.02.051 . Ehsanul Kabir, Pawan Kumar, Sandeep Kumar, Adedeji A. Adelodun, and Ki Hyun Kim. Solar energy: Potential and future prospects. Renewable and Sustainable Energy Reviews, 82(August 2017):894900, 2018. ISSN 18790690. doi: 10.1016/j. rser.2017.09.094. Ted Trainer. Some problems in storing renewable energy. Energy Policy, 110(July): 386393, 2017. ISSN 03014215. doi: 10.1016/j.enpol.2017.07.061. URL https:// doi.org/10.1016/j.enpol.2017.07.061 . Eric Viardot. The role of cooperatives in overcoming the barriers to adoption of renewable energy. Energy Policy, 63:756764, 2013. ISSN 03014215. doi: 10.1016/j.enpol. 2013.08.034. URL http://dx.doi.org/10.1016/j.enpol.2013.08.034 . Ken Willis, Riccardo Scarpa, Rose Gilroy, and Neveen Hamza. Renewable energy adoption in an ageing population: Heterogeneity in preferences for micro-generation technology adoption. Energy Policy, 39(10):60216029, 2011. ISSN 03014215. doi: 10.1016/j.enpol.2011.06.066. URL http://dx.doi.org/10.1016/j.enpol.2011.06. 066 . Bibliography 115 Staan Jacobsson and Anna Johnson. The diusion of renewable energy technology: An analytical framework and key issues for research. Energy Policy, 28(9):625640, 2000. ISSN 03014215. doi: 10.1016/S0301-4215(00)00041-0. John K. Kaldellis and D. Zarakis. The wind energy (r)evolution: A short review of a long history. Renewable Energy, 36(7):18871901, 2011. ISSN 09601481. doi: 10.1016/j.renene.2011.01.002. Bent Sørensen. A history of renewable energy technology. Energy Policy, 19(1):812, 1991. ISSN 03014215. doi: 10.1016/0301-4215(91)90072-V. Christiane Zar, Alexander E. Lumsdon, Jürgen Berlekamp, Laura Tydecks, and Klement Tockner. A global boom in hydropower dam construction. Aquatic Sciences, 77(1):161170, 2014. ISSN 14209055. doi: 10.1007/s00027-014-0377-0. Akhtar Hussain, Syed Muhammad Arif, and Muhammad Aslam. Emerging renewable and sustainable energy technologies: State of the art. Renewable and Sustainable Energy Reviews, 71(December 2016):1228, 2017. ISSN 18790690. doi: 10.1016/j. rser.2016.12.033. URL http://dx.doi.org/10.1016/j.rser.2016.12.033 . Yogesh Kumar, Jordan Ringenberg, Soma Shekara Depuru, Vijay K. Devabhaktuni, Jin Woo Lee, Efstratios Nikolaidis, Brett Andersen, and Abdollah Afjeh. Wind energy: Trends and enabling technologies. Renewable and Sustainable Energy Reviews, 53:209224, 2016. ISSN 18790690. doi: 10.1016/j.rser.2015.07.200. Bhubaneswari Parida, S. Iniyan, and Ranko Goic. A review of solar photovoltaic technologies. Renewable and Sustainable Energy Reviews, 15(3):16251636, 2011. ISSN 13640321. doi: 10.1016/j.rser.2010.11.032. URL http://dx.doi.org/10.1016/j. rser.2010.11.032 . Christian Breyer, Dmitrii Bogdanov, Arman Aghahosseini, Ashish Gulagi, Michael Child, Ayobami Solomon Oyewo, Javier Farfan, Kristina Sadovskaia, and Pasi Vainikka. Solar photovoltaics demand for the global energy transition in the power sector. Progress in Photovoltaics: Research and Applications, 26(8):505523, 2018. ISSN 1099159X. doi: 10.1002/pip.2950. Ze Dong Cheng, Xue Ru Zhao, Ya Ling He, and Yu Qiu. A novel optical optimization model for linear Fresnel reector concentrators. Renewable Energy, 129:486499, 2018. ISSN 18790682. doi: 10.1016/j.renene.2018.06.019. URL https://doi.org/ 10.1016/j.renene.2018.06.019 . Evangelos Bellos, Emmanouil Mathioulakis, Elias Papanicolaou, and Vassilis Belessiotis. Experimental investigation of the daily performance of an integrated linear Fresnel 116 Bibliography reector system. Solar Energy, 167(April):220230, 2018. ISSN 0038092X. doi: 10. 1016/j.solener.2018.04.019. URL https://doi.org/10.1016/j.solener.2018.04. 019 . Md Tasbirul Islam, Nazmul Huda, A. B. Abdullah, and R. Saidur. A comprehensive review of state-of-the-art concentrating solar power (CSP) technologies: Current status and research trends. Renewable and Sustainable Energy Reviews, 91 (November 2017):9871018, 2018. ISSN 18790690. doi: 10.1016/j.rser.2018.04.097. URL https://doi.org/10.1016/j.rser.2018.04.097 . Craig S. Turchi, Judith Vidal, and Matthew Bauer. Molten salt power towers operating at 600−650◦ C: Salt selection and cost benets. Solar Energy, 164(November 2017): 3846, 2018. ISSN 0038092X. doi: 10.1016/j.solener.2018.01.063. URL https:// doi.org/10.1016/j.solener.2018.01.063 . Alexander Bonk, Salvatore Sau, Nerea Uranga, Marta Hernaiz, and Thomas Bauer. Advanced heat transfer uids for direct molten salt line-focusing CSP plants. Progress in Energy and Combustion Science, 67:6987, 2018. ISSN 03601285. doi: 10.1016/j. pecs.2018.02.002. URL https://doi.org/10.1016/j.pecs.2018.02.002 . Alejandro Calderón, Anabel Palacios, Camila Barreneche, Mercè Segarra, Cristina Prieto, Alfonso Rodriguez-Sanchez, and A. Inés Fernández. High temperature systems using solid particles as TES and HTF material: A review. Applied Energy, 213(June 2017):100111, 2018. ISSN 03062619. doi: 10.1016/j.apenergy.2017.12.107. URL https://doi.org/10.1016/j.apenergy.2017.12.107 . Abdelrahman El-Leathy, Sheldon Jeter, Hany Al-Ansary, Syed Noman Danish, Rageh Saeed, Said Abdel-Khalik, Matthew Golob, Eldwin Djajadiwinata, and Zeyad Al- Suhaibani. Thermal performance evaluation of lining materials used in thermal energy storage for a falling particle receiver based CSP system. Solar Energy, 178 (December 2018):268277, 2019. ISSN 0038092X. doi: 10.1016/j.solener.2018.12.047. URL https://doi.org/10.1016/j.solener.2018.12.047 . Simone Polimeni, Marco Binotti, Luca Moretti, and Giampaolo Manzolini. Comparison of sodium and KCl-MgCl2 as heat transfer uids in CSP solar tower with sCO2 power cycles. Solar Energy, 162(September 2017):510524, 2018. ISSN 0038092X. doi: 10.1016/j.solener.2018.01.046. M. R. Rodríguez-Sánchez, A. Sánchez-González, and D. Santana. Feasibility study of a new concept of solar external receiver: Variable velocity receiver. Applied Thermal Engineering, 128:335344, 2018. ISSN 13594311. doi: 10.1016/j.applthermaleng.2017. 08.173. Bibliography 123 Maaike Hoogeboom, Maria Battarra, Güne³ Erdogan, and Daniele Vigo. Erratum- Exact Algorithms for the Clustered Vehicle Routing Problem. Operations Research, 64(2):456457, 2016. ISSN 0030-364X. doi: 10.1287/opre.2016.1491. URL http: //pubsonline.informs.org/doi/10.1287/opre.2016.1491 . Caroline Prodhon and Christian Prins. A survey of recent research on location-routing problems. European Journal of Operational Research, 238(1):117, 2014. ISSN 03772217. doi: 10.1016/j.ejor.2014.01.005. URL http://dx.doi.org/10.1016/j. ejor.2014.01.005 . Thomas Fend, Bernhard Hoschmidt, Robert Pitz-Paal, Oliver Reutter, and Peter Rietbrock. Porous materials as open volumetric solar receivers: Experimental determination of thermophysical and heat transfer properties. Energy, 29(5-6):823833, 2004. ISSN 03605442. doi: 10.1016/S0360-5442(03)00188-9. Marco Astol, Marco Binotti, Simone Mazzola, Luca Zanellato, and Giampaolo Manzolini. Heliostat aiming point optimization for external tower receiver. Solar Energy, pages 116, 2016. ISSN 0038-092X. doi: 10.1016/j.solener.2016.03.042. URL http://dx.doi.org/10.1016/j.solener.2016.03.042 . Boris Belhomme, Robert Pitz-Paal, and Peter Schwarzbözl. Optimization of Heliostat Aim Point Selection for Central Receiver Systems Based on the Ant Colony Optimization Metaheuristic. Journal of Solar Energy Engineering, 136(February 2014):011005, 2013. ISSN 0199-6231. doi: 10.1115/1. 4024738. URL http://solarenergyengineering.asmedigitalcollection.asme. org/article.aspx?doi=10.1115/1.4024738 . M. Berenguel, F. R. Rubio, A. Valverde, P. J. Lara, M. R. Arahal, E. F. Camacho, and M. López. An articial vision-based control system for automatic heliostat positioning oset correction in a central receiver solar power plant. Solar Energy, 76(5):563575, 2004. ISSN 0038092X. doi: 10.1016/j.solener.2003.12.006. Abraham Kribus, Irina Vishnevetsky, Amnon Yogev, and Tatiana Rubinov. Closed loop control of heliostats. Energy, 29(5-6):905913, 2004. ISSN 03605442. doi: 10.1016/S0360-5442(03)00195-6. N. C. Cruz, J. L. Redondo, J. D. Álvarez, M. Berenguel, and P. M. Ortigosa. A parallel Teaching-Learning-Based Optimization procedure for automatic heliostat aiming. The Journal of Supercomputing, 2016. ISSN 0920- 8542. doi: 10.1007/s11227-016-1914-5. URL http://link.springer.com/10.1007/ s11227-016-1914-5 . 124 Bibliography R. Baños, F. Manzano-Agugliaro, F. G. Montoya, C. Gil, A. Alcayde, and J. Gómez. Optimization methods applied to renewable and sustainable energy: A review. Renewable and Sustainable Energy Reviews, 15(4):17531766, 2011. ISSN 13640321. doi: 10.1016/j.rser.2010.12.008. J. Alonso, F. J. Batlles, G. López, and A. Ternero. Sky camera imagery processing based on a sky classication using radiometric data. Energy, 68(May):599608, 2014. ISSN 03605442. doi: 10.1016/j.energy.2014.02.035. URL http://dx.doi.org/10. 1016/j.energy.2014.02.035 . J. Alonso-Montesinos and F. J. Batlles. Solar radiation forecasting in the short- and medium-term under all sky conditions. Energy, 83:387393, 2015. ISSN 03605442. doi: 10.1016/j.energy.2015.02.036. URL http://dx.doi.org/10.1016/j.energy.2015. 02.036 . E. Carrizosa, C. Domínguez-Bravo, E. Fernández-Cara, and M. Quero. Optimization of multiple receivers solar power tower systems. Energy, 90:20852093, 2015b. ISSN 03605442. doi: 10.1016/j.energy.2015.08.005. Mark Schmitz, Peter Schwarzbözl, Reiner Buck, and Robert Pitz-Paal. Assessment of the potential improvement due to multiple apertures in central receiver systems with secondary concentrators. Solar Energy, 80(1):111120, 2006. ISSN 0038092X. doi: 10.1016/j.solener.2005.02.012. Caroline Le Floch, Francois Belletti, Samveg Saxena, Alexandre M. Bayen, and Scott Moura. Distributed optimal charging of electric vehicles for demand response and load shaping. Proceedings of the IEEE Conference on Decision and Control, 54rd IEEE(Cdc):65706576, 2015. ISSN 07431546. doi: 10.1109/CDC.2015.7403254. Raghuram Ranganathan, Wasfy Mikhael, Nasser Kutkut, and Issa Batarseh. Adaptive sun tracking algorithm for incident energy maximization and eciency improvement of PV panels. Renewable Energy, 36(10):26232626, 2011. ISSN 09601481. doi: 10.1016/j.renene.2010.06.011. URL http://dx.doi.org/10.1016/j.renene.2010. 06.011 . Zexian Liu and Hongwei Liu. An ecient gradient method with approximate optimal stepsize for large-scale unconstrained optimization. Numerical Algorithms, 78(1): 2139, 2018. ISSN 15729265. doi: 10.1007/s11075-017-0365-2. Garth P. McCormick. A Modication of Armijo ' S Step-Size Rule for Negative Curvature. Mathematical Programming, 13:111115, 1977. Emilio Carrizosa, Nenad Mladenovi¢, and Raca Todosijevi¢. Variable neighborhood search for minimum sum-of-squares clustering on networks. European Journal of Bibliography 125 Operational Research, 230(2):356363, 2013. ISSN 03772217. doi: 10.1016/j.ejor. 2013.04.027. Rahul Rahul Mark Fonseca, Bailian Chen, Jan Dirk Jansen, and Albert Reynolds. A Stochastic Simplex Approximate Gradient (StoSAG) for optimization under uncertainty. International Journal for Numerical Methods in Engineering, 109(13):1756 1776, 2017. ISSN 10970207. doi: 10.1002/nme.5342. Mark Schmidt, Nicolas Le Roux, and Francis Bach. Minimizing Finite Sums with the Stochastic Average Gradient. pages 83112, 2013. ISSN 14364646. doi: 10.1007/ s10107-016-1030-6. URL http://arxiv.org/abs/1309.2388 . Honggang Wang. Multi-objective retrospective optimization using stochastic zigzag search. European Journal of Operational Research, 263(3):946960, 2017. ISSN 03772217. doi: 10.1016/j.ejor.2017.06.039. URL http://dx.doi.org/10.1016/j. ejor.2017.06.039 . Verena Schmid. Solving the dynamic ambulance relocation and dispatching problem using approximate dynamic programming. European Journal of Operational Research, 219(3):611621, 2012. ISSN 03772217. doi: 10.1016/j.ejor.2011.10.043. URL http://dx.doi.org/10.1016/j.ejor.2011.10.043 . Victor Pillac, Michel Gendreau, Christelle Guéret, and Andrés L. Medaglia. A review of dynamic vehicle routing problems. European Journal of Operational Research, 225(1):111, 2013. ISSN 03772217. doi: 10.1016/j.ejor.2012.08.015. URL http: //dx.doi.org/10.1016/j.ejor.2012.08.015 . Ivar Ekeland and Roger Temam. Convex analysis and variational problems. Elsevier, 1976. ISBN 9780080875224. Jorge Nocedal, Stephen J Wright, and Stephen M Robinson. Numerical Optimization. Springer, 1 edition, 1999. ISBN 0387987932. doi: 10.1007/b98874. Marco Astol, Marco Binotti, Simone Mazzola, Luca Zanellato, and Giampaolo Manzolini. Heliostat aiming point optimization for external tower receiver. Solar Energy, 157:11141129, 2017b. ISSN 0038092X. doi: 10.1016/j.solener.2016.03.042. URL https://doi.org/10.1016/j.solener.2016.03.042 . H. Breitkreuz, M. Schroedter-Homscheidt, and T. Holzer-Popp. A case study to prepare for the utilization of aerosol forecasts in solar energy industries. Solar Energy, 81(11): 13771385, 2007. ISSN 0038092X. doi: 10.1016/j.solener.2007.01.009. Bruno Dürr. Automatic cloud amount detection by surface longwave downward radiation measurements. Journal of Geophysical Research, 109(D5):D05201, 2004. ISSN 126 Bibliography 0148-0227. doi: 10.1029/2003JD004182. URL http://doi.wiley.com/10.1029/ 2003JD004182 . Germain Augsburger and Daniel Favrat. Modelling of the receiver transient ux distribution due to cloud passages on a solar tower thermal power plant. Solar Energy, 87 (1):4252, 2013. ISSN 0038092X. doi: 10.1016/j.solener.2012.10.010. R. Tapakis and A. G. Charalambides. Equipment and methodologies for cloud detection and classication: A review. Solar Energy, 95:392430, 2013. ISSN 0038092X. doi: 10.1016/j.solener.2012.11.015. URL http://dx.doi.org/10.1016/j.solener.2012. 11.015 .