A novel two-model based approach for optimal scheduling in CSP plants

Available online at www.sciencedirect.com

ScienceDirect Solar Energy 126 (2016) 73–92 www.elsevier.com/locate/solener

A novel two-model based approach for optimal scheduling in CSP plants Manuel Jesu´s Vasallo ⇑, Jose´ Manuel Bravo Department of Electronic Engineering, Computer Systems, and Automatics, University of Huelva, Carretera Palos de La Frontera s/n, 21819 Palos de La Frontera, Huelva, Spain Received 29 May 2015; received in revised form 18 September 2015; accepted 22 December 2015 Available online 13 January 2016 Communicated by: Associate Editor Aguilar Jose´ Gonza´lez-Aguilar

Abstract Thermal energy storage (TES) systems allow concentrated solar power (CSP) producers to participate in a day-ahead market. Then, the optimal power scheduling problem can be posed, whose objective is the maximisation of profits derived from electricity sales. Most papers in literature use a mixed-integer linear programming (MILP) approach to solve this type of problems. This paper proposes a novel approach based on the use of two models: a detailed model and a MILP model. This approach combines MILP capabilities and the accuracy of a detailed model. The proposed approach is applied to a 50 MW parabolic-trough-collector based CSP plant with molten-salt-based TES. A detailed model available in literature and validated against operating plant data is used, but some improvements are included for its use in optimal scheduling problems. Moreover, the MILP model was developed to adjust as much as possible to the features of the detailed model. The improvements regarding other scheduling strategies for a specific example are shown. Ó 2016 Elsevier Ltd. All rights reserved.

Keywords: Concentrating solar power plant; Thermal energy storage; Electricity market; Optimized operation strategy; Parabolic trough collector; Mixed-integer linear programming

Variables associated with a time period refer generally to average values. 1. Introduction Amongst renewable energy sources, concentrating solar power (CSP) is a promising technology that is already being exploited in some countries such Spain and USA, where subsidy policies promote its development. There are mainly four commercially-available CSP technologies: (1) linear Fresnel reflector, (2) paraboloid dish, (3) solar ⇑ Corresponding author. Tel.: +34 959217376; fax: +34 959217348.

E-mail addresses: [email protected] [email protected] (J.M. Bravo). http://dx.doi.org/10.1016/j.solener.2015.12.041 0038-092X/Ó 2016 Elsevier Ltd. All rights reserved.

(M.J.

Vasallo),

power tower, and (4) parabolic trough collector (PTC). PTC-based plants, using either synthetic or organic oil as the heat transfer fluid (HTF), are the most commerciallyattractive option (Purohit et al., 2013), as well as the most-widely installed CSP technology (Zhang et al., 2013). Interest in CSP technology is mainly based on the semidispatchable nature of CSP plants with thermal energy storage (TES) and/or fossil-fuel backup systems. These systems allow reducing real-time net power variability in case of insufficient solar energy, extending the whole production period and displacing production periods towards highprice periods. Specifically, the advantages of CSP with TES have been experimentally proven in Dinter and Gonzalez (2014), where a similar operability to that of a conventional fossil-fired power plant is achieved.

74

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

Nomenclature approximated average value for the fraction of thermal power in the PB inlet used to increase thermal state in the PB (kW) during a Dto -long period with startup (kW) a1 ; b1 parameters of linear approximation of function turbine-generated gross electric power against thermal power entering the PB (dimensionless and kW) M high enough numeric value DP HTF stu max maximum average value of the fraction of thermal power in PB inlet used to generate electricity during an hourly period when a startup is performed (kW) pc percentage to establish a minimum point when the turbine starts T min u ; T min d minimum up and down times for the turbine (in time steps of the MILP model) L; N 0 initial times during which the turbine must be online and offline (in time steps of the MILP model)

P start Indices, superscripts and time steps t; k; j time indices for the SF submodel, the detailed model and the MILP model t0 time since the beginning of the active turbine operating phase MILP ;  superscripts for the MILP model and the three-step procedure outputs Dt; Dtp ; Dto time steps of the solar field submodel, the detailed model and the MILP model Parameters N loops number of loops in the SF tsunset ; tphase1 end sunset time and time of the end of phase 1 in the SF state machine tstartup duration of the turbine startup N st number of simulation steps (Dtp ) within the startup time T antifreeze ON ; T antifreeze OFF hysteresis values for the freeze-protection system (°C) T in loop des ; T out loop des design HTF temperatures for loop inlet and outlet (°C) T in disch des ; T out disch des design HTF temperatures for TES inlet and outlet during discharge (°C) DT margin temperature margin parameter in the SF state machine (°C) P abs loop thr threshold parameter for loop-absorbed useful thermal power in the SF state machine (kW) gexch TES efficiency of the HTF-TES heat exchangers gexch efficiency of the HTF-water/steam heat exchangers ggross2net ratio for considering parasitic electric consumption P HTFmin ; P HTFmax minimum and maximum values for thermal power entering the PB (kW) P cmin ; P dmin minimum values for charging and discharging thermal power on the HTF side due to minimum value in salt pump flow rate (kW) P cmax exch ; P dmax exch maximum values for charging and discharging thermal power on the HTF side due to heat-exchanger limitations (kW) P dmax onlyTES maximum value for discharging thermal power on the HTF side in the TES-only mode due to the input limitations of the turbine in this mode (kW) Emin ; Emax minimum and maximum values for energy in the TES Parameters only for the MILP model N number of steps in optimisation time horizon C DP e cost of the gradient of turbine-generated gross power (Eur/kW h)

up

Variables DNIðt=jÞ direct normal irradiance either in instant t or in period j (kW/m2) T a ðtÞ ambient temperature in instant t (°C) T 1 ðtÞ to T 4 ðtÞ; T pipes ðtÞ; T in ðtÞ average temperatures of solar collector assemblies and insulated pipes section, and the SF inlet temperature in instant t (°C) qmloop ðt=kÞ HTF mass flow rate in a loop either in instant t or in period k (kg/s) qTES2PB ðkÞ HTF mass flow rate from the TES to the PB during discharge in period k (kg/s) p1 ðkÞ; p2 ðkÞ fraction of total HTF mass flow rate to the PB in mixed mode from the SF and from the TES in period k P abs loop ðtÞ loop-absorbed useful thermal power in instant t (kW) P SFmax ðt=k=jÞ maximum thermal power available from the SF on the HTF side either in instant t or in periods k or j (kW) P SF ðk=jÞ thermal power actually transferred from the SF on the HTF side in periods k or j (kW) P SF def ðk=jÞ defocused thermal power in the SF on the HTF side in periods k or j (kW) P HTF in ðt0 =kÞ thermal power entering the PB either in instant t0 or in period k (kW) 0 P HTF ðt =k=jÞ fraction of thermal power in PB inlet used to generate electricity either in instant t0 or in periods k or j (kW) P startup ðt0 =k=jÞ fraction of thermal power in PB inlet used to increase thermal state in the PB during startup either in instant t0 or in periods k or j (kW)

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

in ðjÞ

steady-state thermal power in the PB inlet necessary to meet setpoint in period j Estartup setpoint-dependent energy required to increase the PB thermal state during the startup process P e ðk=jÞ; P enet ðk=jÞ turbine-generated gross and net electric power in periods k or j (kW) P el ðjÞ value of P e ðjÞ calculated by linear approximation only valid when the turbine is running (kW) P e SP ðk=jÞ setpoint for the turbine-generated gross electric power in periods k or j (kW) ePe ðkÞ electricity generation error in relation to setpoint P e SP ðkÞ in period k (kW) DP e abs ðjÞ absolute value of the gradient of turbinegenerated gross power in period j (kW) P c ðk=jÞ; P d ðk=jÞ charging and discharging thermal power on the HTF side in periods k or j (kW) P cmax ðkÞ; P dmax ðkÞ maximum values for charging and discharging thermal power on the HTF side in period k (kW) P cmax tank ðkÞ; P dmax tank ðkÞ maximum values for charging and discharging thermal power on the HTF side due to the energy level in the TES in period k (kW) P dmax mixed ðk=jÞ maximum value for discharging thermal power in mixed mode on the HTF side in periods k or j (kW) P dmax mixed l ðjÞ previously calculated value of P dmax mixed ðjÞ only valid when the TES is discharging (kW) Eðk=jÞ TES energy level at the beginning of periods k or j (kW h) pðjÞ market price prediction in period j (Eur) P HTF

Binary variables in detailed model Phase6 equal to 1 if the SF is ready to transfer energy Back2phase1 equal to 1 if the transition between phases 7 to 1 has already been made on the current day (SF state machine) ready to turb startup equal to 1 if enough energy is available for turbine startup stop to turb startup equal to 1 if a shortage of energy happens during turbine startup stop to turb equal to 1 if an order to stop the turbine has been given SFstupready equal to 1 if the startup can be executed from the SF

TES systems can be classified as active or passive types depending on the storage medium: it is liquid in the active type and solid in the passive one. The liquid in the active type can also be used as the HTF. These systems are known as direct-active systems. On the contrary, the expression indirect-active systems is used when both fluids are different

75

TESstupready equal to 1 if the startup can be executed from the TES mixedstupready equal to 1 if the startup can be executed in mixed mode Binary variables in MILP model aðjÞ equal to 1 if the turbine is working during period j bðjÞ equal to 1 if the turbine has started during period j kðjÞ equal to 1 if the turbine has stopped during period j cðjÞ equal to 1 if the TES is charging during period j dðjÞ equal to 1 if the TES is discharging during period j Functions gmixed ð:Þ; gonlySF ð:Þ; gonlyTES ð:Þ the turbine energy efficiency curves in mixed, solar-only and TESonly modes H ðT Þ HTF enthalpy function versus temperature (kJ/kg) Metrics UCENE unmet committed net electric energy per day rUCENE relative unmet committed net electric energy compared to the net electric energy that is effectively served PC penalty cost per day OP obtained profit per day rPI relative profit increase regarding the strategy d-MILPs rPEE relative profit estimation error Acronyms CSP DNI HCE HTF MILP PB PTC SCA SF TES

concentrated solar power direct normal irradiance heat collection element heat transfer fluid mixed-integer linear programming power block parabolic trough collector solar collector assembly solar field thermal energy storage

and an additional heat exchanger is needed. Nowadays the most commercially mature technology for TES is an indirect two-tank molten salt storage system (Kuravi et al., 2013). Simulation is an useful tool to design and analyse engineering systems. Intensive effort has been made on modelling

76

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

and simulation of CSP plants (Reddy et al., 2013). A complete list of software tools to simulate the varied components in CSP systems is presented in Ho (2008). One of the most widely-used tools is Solar Advisor Model, SAM (SAM, 2015), developed by the National Renewable Energy Laboratory in USA. This tool can simulate the hourly output of photovoltaic, concentrating photovoltaic and CSP plants to provide annual performance. The tool also includes a financial model to facilitate decision making processes. As mentioned above, TES allows CSP systems to be a semi-dispatchable power source, as TES’s limited storage capacity prevents it from being totally dispatchable. In spite of this limitation, CSP systems may participate in electricity markets. An important property of TES is their capability to displace production from low-price to highprice periods. Therefore, the optimal generation scheduling problem can be approached. In a deregulated market, the objective of electricity producers is to maximise their total profits from energy sales. Moreover, in electricity markets, power plant owners must offer in advance a daily production schedule. Consequently, electricity price and weather forecast must be taken into account in the optimisation problem. Literature on optimal CSP plant operation is still limited. One of the first works was Sioshansi and Denholm (2010). This work studies the profits obtained by a CSP plant with TES when participating in electricity market in various American regions. Two models are used to solve the optimisation problem: the SAM tool, which provides the thermal power produced in the solar field, and an optimisation model based on mixed-integer linear programming (MILP). The optimisation problem for CSP plant operation in the Spanish market are formulated in Wittmann et al. (2011), but no details about models used are given. An optimisation approach to maximise profits of CSP with TES and gas backup systems in Spain is presented in Usaola (2012), where the MILP approach is used again. Another example of a MILP approach applied to the Spanish market can be found in Kost et al. (2013), where the model in Sioshansi and Denholm (2010) was extended to include a gas backup system and an electric heater for grid energy storage. Different approaches to MILP can be found in Lizarraga-Garcia et al. (2013) and Powell et al. (2014), where Nonlinear Programming is used, and in Channon and Eames (2014), where Dynamic Programming is used. One of the disadvantages of CSP plants is the predictability of its electricity production, which is limited by the forecasting accuracy of direct normal irradiance (DNI). As a result, operation in electric markets generates risks of penalties for deviation from the scheduled generation. These costs are directly related to the size of deviation. Therefore, the importance of DNI forecast accuracy is clear (Law et al., 2014). A study on the mitigation of this impact on deviation penalties offered by a state-of-the-art DNI forecast tool for the Spanish market can be found in Kraas et al. (2013).

Uncertainty in CSP production forecast – together with the limited forecasting accuracy of market prices – justifies the use of robust and stochastic approaches to solve the optimal scheduling problem (Dominguez et al., 2012; Pousinho et al., 2015). Both works solve optimal operation by a MILP context. Despite the foregoing, it is difficult to find today’s commercial dispatchable CSP plants operating with an optimal strategy such as the ones mentioned above (Alliance C, 2014). That is, optimising storage use according to market price forecast is not a common practice. According to literature review, MILP is probably the most widely-used method to model optimal scheduling problems in CSP plants. In general, MILP is a powerful method for optimal scheduling problems in power systems (Simoglou et al., 2012). Its ability to formulate unit commitment, startup, shutdown and ramp-rate constraints, as well as piecewise linear functions, explains its use in this type of problems. Moreover, some efficient MILP solvers can currently find global optimal solutions in short computation times. On the other hand, accurate CSP-plant performance must be described by a hybrid dynamical model. This type of model poses a hierarchical structure characterised by continuous-variable dynamics and logical decision making at the lowest and highest level, respectively. The highest level describes the evolution of operating phases in the plant. An example of this hybrid nature can be found in Garcı´a et al. (2011), where the performance of a detailed model of a 50 MW PTC-based CSP plant with moltensalt-based TES is presented. The model was validated using real operating data. The conclusions drawn from this work include the need to consider the actual plant’s operation philosophy in the model. This philosophy can be easily implemented at the logical decision-making level. Nevertheless, the MILP approaches for CSP plants found in literature use relatively simple optimisation models and do not take into account the studied plant’s actual operation philosophy. Moreover, time step is usually an hour, when all variables are considered constant. A mechanism to consider a more detailed model for a specific CSP plant in the aforementioned sense could provide more accurate results in scheduling problems. This paper proposes a novel MILP-based approach to find optimal scheduling for a CSP plant that maximises profits from energy sales in the electricity market. The novelty is based on using two models in the problem: (1) an enough-detailed CSP plant model and (2) a MILP model. This way, MILP capabilities and a detailed model’s accuracy are combined. The detailed model must be aimed at capturing the significant aspects that sufficiently describe the studied plant’s performance. On the other hand, the MILP model must be formulated to adjust as much as possible to the features of the detailed model, but preserving enough simplicity to avoid increasing computational burden in a prohibitive way. A three-step procedure is proposed to use both models: (1) generation of the thermal

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

production in solar field by means of the detailed model, (2) solution of the MILP model, and (3) refinement of the previous solution by a new simulation run on the detailed model. The context of this work is deterministic (i.e., market price and DNI predictions are assumed to be perfect). The effect of inaccurate forecast is currently under research. The proposed approach is applied to a 50 MW PTCbased CSP plant with molten-salt-based TES. The detailed model of this plant, available in Garcı´a et al. (2011), was improved to complete certain aspects which were not described and to provide flexibility for its use in optimal scheduling problems. The original model is not valid for this purpose. In summary, the contribution of this paper is threefold: 1. A detailed dynamic model of a PTC-based CSP plant with TES taken from literature was improved to be used in optimal scheduling problems. 2. A novel two-model based approach is proposed for optimal scheduling problems of CSP plants in the electricity market. This approach is based on: (a) the formulation of an enough-detailed model of the studied plant, (b) the development of a MILP model adjusted as much as possible to the detailed model, and (c) the application of a three-step procedure. 3. This new approach is applied to a specific case of a 50 MW PTC-based CSP plant. PTC-based CSP plants with TES are described in general in Section 2. The proposed approach is explained in Section 3, and then applied to a specific example in Section 4. The original, detailed and MILP models are described in Sections 4.1, 4.2, and 4.3, respectively. Certain aspects on data transfer between the models are addressed in Section 4.4. Results, comparison with other methods and discussion are put forward in Section 4.5. Finally, conclusions are drawn in Section 5.

Fig. 1. Simplified diagram of a PTC-based CSP plant with TES.

77

2. Description of a PTC-based CSP plant with TES The approach proposed in this work is generally applicable to CSP plants with TES. Specifically, it shall be applied to a 50 MW PTC-based CSP plant with moltensalt-based TES. This type of plant shall be described next. Fig. 1 shows a simplified diagram of these plants, which comprise three main blocks: the Solar Field (SF), the Power Block (PB), and the TES system. The liquids circulating in each of these blocks are the HTF; the HTF, water and steam; and the HTF and molten salt, respectively. A heat exchanger system allows exchanging energy between the HTF and the molten salt in a bidirectional way, and another exchanger system allows transferring energy from the SF and/or TES to the PB. The SF is composed of a number of PTC loops connected in parallel to each other. Several solar collector assemblies (SCA) in series form a loop. Each SCA consists in turn of several solar collector elements (SCE), which are formed by a reflective parabolic mirror and receiver tubes installed at the focal line of the parabolic surface. The receiver tubes, also called heat collection elements (HCE), contain the circulating HTF. Insulated pipes distribute the HTF to the loops. The collectors follow the sun thanks to a single-axis tracking mechanism. This way, a large portion of solar power is concentrated onto the HCEs and transferred to the HTF by increasing its temperature. The PB is composed of a heat exchanger train, a steam turbine coupled to an electricity generator, a condenser, a cooling system and other auxiliary elements. The TES system is an indirect two-tank molten salt system. Molten salt circulates from the cold (hot) tank to the hot (cold) tank in energy charging (discharging) modes. As mentioned above, a bidirectional heat exchanger system allows energy transfer between the HTF and the molten salt. Ultimately, a fossil fuel backup system can be used to prevent HTF from freezing, to keep the temperature of molten salt, or to assist production in case of insufficient solar energy. Several power operating modes can be run in CSP plants (e.g., SF ? PB, SF ? TES, SF ? PB + TES, SF + TES ? PB, TES ? PB), and others when a fossil fuel backup system is taken into account. Moreover, plant is operated according to a sequence of phases. For example, nocturnal recirculation, freezing protection, HTF warm-up period, and sunlight period are typical operating phases in the SF. Different HTF temperature conditions guide the transition between each of these phases, in which a value for HTF mass flow rate is fixed. All these details set the SF’s operation philosophy. In the sunlight period, the main SF pumps are regulated to control HTF temperature in loops’ output. A value around 390 °C is the typical setpoint that guarantees HTF protection against high temperatures. High solar radiation or relatively low PB and/or TES power demand lead to increased risk of high HTF temperature despite HTF flow rate having reached its maximum value. In this situation, a number of collectors must be deliberately defocused. In

78

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

relation to TES operation, charging (discharging) power to (from) the TES is set by HTF flow rate to (from) the TES. Furthermore, a control system regulates salt flow rate to keep the values of salt temperature coming into both tanks close to design values. 3. Description of the proposed approach for optimal scheduling problems in CSP plants From the point of view of an energy producer in the electricity market, the optimal scheduling problem usually has a time horizon of one day. The production schedule for a specific day must be given the previous day. Inputs for the optimal scheduling problem are the weather and market price forecast. The objective is to maximise total profits from energy sales. This work assumes the day-ahead market and the producer’s price-taking property as common (i.e., its production schedules have no feedback on market prices). Moreover, the context of this work is deterministic (i.e., market price and DNI predictions are assumed to be perfect). The effect of inaccurate forecast is currently under research. As seen in Section 2, CSP plant performance has a hybrid nature. That is, it shows both continuous and discrete dynamic behaviour. Whilst continuous variables represent physical quantities such as temperature, mass flow, and power, discrete variables are associated with operating phases. The discrete dynamic behaviour is related to the plant’s operation philosophy. The need to consider actual plant’s operation philosophy in the model to attain accurate simulation results is pointed out in Garcı´a et al. (2011). Nevertheless, the MILP approaches for CSP plants found in literature do not take into account the studied plant’s actual operation philosophy. Moreover, they use relatively simple MILP models (e.g., with simple equipment performance constraints). In literature, SF thermal production is commonly calculated by a detailed model, e.g., the SAM tool (Sioshansi and Denholm, 2010; Madaeni et al., 2012; Kost et al., 2013). This information is used later as an input to the MILP model. This is the only function of the SAM tool. The SAM tool model is very detailed and is parameterized to adapt to a specific plant, but it is relatively far from being an ad hoc model for the studied plant. For example, it is difficult to include the plant’s operation philosophy in the SAM tool. In addition, time steps in SAM and MILP models are of one hour. Shorter time steps can be proposed to improve accuracy. The approach proposed in this paper is based on MILP to solve the optimisation problem. Thus, MILP capabilities are exploited: (1) ability to formulate unit commitment, startup, shutdown and ramp-rate constraints, as well as piecewise linear functions; (2) capability to find global optimal solutions; and (3) relatively short computation times. Moreover, a detailed model of the plant is also used to assist in finding the final solution. The detailed model must be aimed at capturing the significant aspects that describe sufficiently the studied plant’s performance. Specifically, a

hybrid dynamic model, in which the plant’s operation philosophy is represented, is proposed. Another suggestion is using shorter time steps (below 1 h). Furthermore, the MILP model must be formulated to adjust as much as possible to the features of the detailed model, but preserving enough simplicity to avoid increasing computational burden in a prohibitive way. Therefore, the equipment performance constraints used in the detailed model are suggested to be added to the MILP model. In this sense, the MILP models found in literature are very simple in relation to constraints. Note that the detailed model is a simulation model whilst the MILP model is an optimisation model. Once both models are available, a three-step procedure must be performed (see Fig. 6, where the proposed approach was applied to a specific example). These three steps are described next: 1. Calculation of the thermal production in the SF by means of the detailed model. This way, some aspects of the SF are not required in the MILP model, simplifying it. 2. Obtaining a previous solution using the MILP model and information from step 1. As a result, a previous production schedule for CSP plant is generated. 3. Refinement of the previous solution. Taking the production schedule of step 2 as the setpoint for production, simulation is run again by the detailed model. The resulting production is the final solution, which is offered to the market as the production schedule. This last step applies the whole complexity of the detailed model to the MILP solution, thus making it closer to reality. An approximation is assumed in this three-step procedure (as in all works using the MILP approach found in literature): decisions about production do not affect the thermal evolution of HTF in the SF (information from step 1). Since partial defocus is a decision variable in optimal scheduling problems, its application has impact on the HTF thermal state by reducing loop-absorbed thermal power. This approximation is based on the fact that the aim of partial defocus is also maintaining loop-outlet HTF temperature around its setpoint (see Sections 4.2.1 and 4.2.2 for further information). Obviously, the solution found by the proposed approach is suboptimal because the MILP model is an approximation of the detailed model. The proximity to the optimal solution is related to the complexity of the MILP model. The ability of the MILP approach to formulate piecewise linear functions would allow a good fit of the MILP model to the detailed model. Nevertheless, computational time may increase significantly due to the higher number of binary variables required. In this sense, the last step of the threestep procedure applies the whole complexity of the detailed model to the MILP solution and, then, allows maintaining a certain level of simplicity in the MILP model, necessary for

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

computation reasons. Note that, in literature, MILP approaches use the solution of relatively simple MILP models as the offered production schedule. This paper aims to gain accuracy compared with these approaches. Another alternative to solve the optimisation problem is using only the detailed model. As the detailed model is a nonlinear hybrid dynamic model, the application of this method involves the use of a mixed-integer non-linear programming approach. Nevertheless, this method does not guarantee finding the global optimal solution. Moreover, computational time could increase meaningfully. In fact, mixed-integer programming problems are classified as NP-complete, which means that in the worst-case scenario, the computational time increases exponentially alongside the problem size (Bemporad, 2003). Therefore, it is advisable to formulate a relatively small MILP model. 4. Aplication of the proposed approach to a 50 MW PTCbased CSP plant with molten-salt-based TES In this section, the general approach proposed in this work is applied to a 50 MW PTC-based CSP plant with molten-salt-based TES. A fossil fuel backup system is also considered, but its only function is to prevent HTF from freezing. The detailed model needed to apply the approach is an improved version of the model available in Garcı´a et al. (2011). A brief description of this original model is given in Section 4.1. This model cannot be used in scheduling problems, and the necessary improvements are also indicated in Section 4.1. The detailed model is described in Section 4.2. The MILP model is explained in Section 4.3, and Section 4.4 details certain aspects about data transfer between models. Finally, results and discussion are put forward in Section 4.5.

79

approximation of the HTF circuit in five sections. Therefore, the main time step of the model is set to, at least, 10 min, but the time step in the SF submodel is fixed to, at most, 10 s. To conclude this brief description of the features of the original model, it must be mentioned that a state machine governs the evolution of the plant’s operation phases. This state machine is obviously based on the plant’s operation philosophy. Some improvements have been made in the original model for its application in an optimal scheduling problem. These improvements are: 1. Possibility of changing the setpoint for the electricity production instead of the full-power setpoint of the original version. 2. Use of two state machines instead of a single one. In the original version, a single state machine characterises the evolution of operating phases. Its configuration only allows performing the following power operating modes: SF ? PB, SF ? PB + TES, SF + TES ? PB, and TES ? PB (nigh-time period). In the new model, one state machine describes the evolution of operating phases in the SF, and another one is in charge of turbine phases, thus allowing some decoupling. This way, the following power operating modes are added: SF ? TES and TES ? PB (at any time). Moreover, turbine startup can be performed from the SF and/or TES instead of startup only from the SF. 3. Some details not described in Garcı´a et al. (2011) were included into the SF state machine (see Section 4.2.3). 4. The turbine startup model is reformulated completing certain aspects not described in Garcı´a et al. (2011) (see Section 4.2.6.2).

4.2. Detailed model 4.1. Original model and necessary improvements The detailed model in Garcı´a et al. (2011) and the improvements necessary for its application in an optimal scheduling problem are described briefly in this subsection. The original model, which was validated by its authors using real operating data, predicts the electricity output by taking into account the detailed component level analysis in the SF. The description of the other components is performed by energy flow balances. The dominant continuous dynamics are those related to the following processes: (1) the thermal evolution of the HTF in the SF, (2) turbine startup, and (3) thermal energy storage in tanks. Other relationships amongst variables are described by steadystate equations. The model is discretised using time steps of, at least, 10 min. The SF submodel is based on the discretisation of a representative collector loop into four sections in longitudinal direction. Another section is added to represent the insulated pipes. Time steps to discretise the differential equations describing the thermal evolution in the SF require values equal or below 10 s to validate the

In this Subsection, the detailed model of a 50 MW PTCbased CSP plant with molten-salt-based TES is described. Since the detailed model is an improved version of the model in Garcı´a et al. (2011), an in-depth description of the whole model is not given. Fig. 2 shows a block diagram of the detailed model. Each box represents a calculation block. Shadow boxes are associated with dynamic equations, whilst the others boxes are associated with steady-state equations. The dynamic equations are discretised with a 10 min time step (Dtp ), except for the SF thermal model, which has a 10 s time step (Dt). The block Power balance consists of steadystate equations but also includes the ramping behaviour of turbine startup. The arrows entering and leaving each box mark information flow between the blocks. For the sake of clarity, only the main flows are indicated. The nomenclature used in models to refer variables is: (1) xðtÞ when time step is 10 s, (2) xðkÞ when time step is 10 min, and (3) xðjÞ when time step is 1 h (in the MILP model).

80

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

Fig. 2. Block diagram of the detailed model. Time step in blocks with underlined text is Dtp . Time step in other blocks is Dt.

The types of input data considered are indicated in Fig. 2. As meteorological data, DNI and ambient temperature are taken into account. The setpoint for the electricity production is also included as input data. The simulation algorithm implementing the detailed model was developed only to allow hourly changes of the setpoint. This way, frequent operation changes in the turbine are avoided. A brief description of the dynamic blocks and the block Power balance are given next: 1. The SF state machine sets the active operating phase in the SF, which lays down the configuration of the HTF circuit and the HTF flow rate, both in the SF. The three possible configurations of the HTF circuit are explained in Section 4.2.2. The time step of this block is Dt. 2. The turbine state machine sets the active operating phase in the turbine. Its time step is Dtp . 3. The SF thermal model calculates the evolution of HTF temperatures in the SF. Its time step is Dt. 4. The block Energy in TES determines the energy level in molten salt tanks according to the previous state and the current charging/discharging thermal power. Its time step is Dtp . 5. The block Power balance takes into account the setpoint for the electricity production, the maximum power available from the SF, the active turbine operating phase, the TES state, the HTF flow rate from the SF, the equipment’s technical limitations and performance curves, and other information. On a 10 min basis, average values for powers are calculated according to the energy conservation law and the decisions of the plant operator and control systems.

Further information of blocks are given in the following Subsections. Values for parameters of the detailed model can be found in Table 1. 4.2.1. Loop-absorbed useful thermal power and SF insulated piping losses From DNI and other input data, power absorbed by a collector loop is calculated. A maximum value around 1.8 MW is applied. This limit is imposed in real plants by defocusing a number of collectors. This way, risk of high HTF temperature due to very high radiation, even in full operation, is avoided. HCE losses in a loop are obtained from loop temperatures and ambient temperature. Then, loop-absorbed useful thermal power is calculated. SF insulated piping losses are obtained from the temperature of insulated pipes and ambient temperature. All the equations used can be found in Garcı´a et al. (2011). 4.2.2. SF thermal model, HTF circuit configuration and maximum thermal power available from the SF A representative loop and a section of insulated pipes make the HTF circuit in the SF (see Fig. 3). The loop is formed by four SCAs. Five nodes are considered in the lumped parameter thermal equations of the HTF circuit (four SCAs and one section of insulated pipes). The HTF circuit is closed in two different ways, depending on the type of active operating phase. If HTF is recirculated in the active operating phase, the HTF flow bypasses the PB (see schema 1 in Fig. 3). On the contrary, HTF temperature at the PB output is fixed at design loop inlet temperature (see schema 2 in Fig. 3). Finally, when the freeze-protection system for HTF is active, three nodes

Table 1 Values for parameters of the detailed model. T in loop des T out loop des T antifreeze OFF T antifreeze ON gexch TES N loops

296 °C 390 °C 100 °C 60 °C 0.95 156

T in disch des T out disch des DT margin P abs loop thr gexch tstartup

290.5 °C 360 °C 80 °C 20 kW 0.95 20 min

P HTFmin P cmin P cmax exch P dmax onlyTES ggross2net Dtp

20 MW 20.24 MW 100 MW 119 MW 0.95 10 min

P HTFmax P dmin P dmax exch Emin Emax Dt

140 MW 15 MW 124 MW 0 1010 MW h 10 s

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

81

that indicates that the operating phase SF ready is active (see Section 4.2.3). Notice that defocus of some collectors is only considered in all these calculations in case of very high radiation. However, decisions concerning production could also involve partial defocus. Therefore, the calculated thermal evolution in the SF is an approximation, and the power available from the SF is described as maximum. The approximation is based on the fact that the aim of partial defocus is also to keep loop outlet HTF temperature around its setpoint.

Fig. 3. Configuration schemas for the HTF circuit in SF.

replace the node of insulated pipes and the intermediate node receives a heat input of 10 MW from the heater (see schema 3 in Fig. 3). Three configuration schemas for the HTF circuit summarise the aforementioned discussion (see Fig. 3), when T 1 –T 4 are the SCA average temperatures and T pipes is the average temperature in the section of insulated pipes. These five temperatures (seven temperatures figures when the freeze-protection system is active) are the state variables of the thermal model. Loop-absorbed useful thermal power, insulated piping losses and HTF mass flow rate are also used in the thermal model. The equations of the thermal model can be found in Garcı´a et al. (2011). Last, the maximum thermal power available from the SF is calculated using Eq. (1). P SFmax ðtÞ ¼ N loops qmloop ðtÞðH ðT pipes ðtÞÞ  H ðT in ðtÞÞÞ  Phase6

ð1Þ

where H ðT pipes ðtÞÞ  H ðT in ðtÞÞ is the enthalpy change between SF outlet and inlet and Phase6 is a binary variable

4.2.3. SF state machine Fig. 4 shows the diagram of the SF state machine, where t is the current time in hours of the current day (e.g., when 25 h have passed in simulation time, t ¼ 1h) and Back2phase1 is a binary variable that indicates that transition between phases 7 and 1 has already been made on the current day. This last binary variable allows returning to phase 1 when t > tsunset regardless of the previous state. Each phase sets a configuration for the HTF circuit in the SF and a value for the HTF mass flow rate in the SF (see Fig. 4). In phases 5 and 6, HTF mass flow rate changes dynamically looking for the setpoint of the loop outlet temperature (function f T ð:Þ expresses this control). The value for the setpoint is set to the design loop outlet temperature. Transition conditions are also shown in Fig. 4. In Phase 6, the SF is ready to transfer thermal energy to the PB and/or TES. Phase 7 is a one-iteration state before returning to phases 2 or 4. 4.2.4. Turbine state machine Fig. 5 shows the diagram of the turbine state machine, where ready to turb startup is a binary variable that indicates that enough energy is available for turbine startup, stop to turb startup is a binary variable that indicates that shortage of energy happens during turbine startup, and stop to turb is a binary variable that indicates that an order to stop the turbine has been given. Variables

Fig. 4. Diagram of the SF state machine.

82

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

Fig. 5. Diagram of the turbine state machine.

stop to turb startup and stop to turb are calculated in block Power balance. Parameter tstartup is set to 20 min (assuming a hot start (Garcı´a et al., 2011)). Variable ready to turb startup needs further explanation, which is given in Section 4.2.6.2. The active turbine operating phase is used in block Power balance to select the calculation method for power flows. 4.2.5. Energy in TES The energy level in molten salt tanks is obtained using Eq. (2): Eðk þ 1Þ ¼ EðkÞ þ PcðkÞDtp gexch  PdðkÞDtp =gexch

TES

TES

ð2Þ

where EðkÞ is the TES energy level at the beginning of simulation step k. 4.2.6. Power balance Block Power balance calculates average values for power flows on a 10 min basis according to the energy conservation law, the decisions of the plant operator and control systems, and the equipment’s technical limits. The setpoint for the electricity production, the maximum power available from the SF, the active turbine operating phase, the TES state and the HTF flow rate from the SF are input data for this calculation (see Fig. 2). This block applies different methods depending on the active turbine operating phase: a first method when Turbine shutdown or Turbine running is active, and a second method when Turbine startup is active. 4.2.6.1. Turbine running or shutdown. In case of turbine running or shutdown, the set of steady-state equations and inequalities (3)–(24) are applied to obtain all power flows in step k. Most of these expressions are obtained from Garcı´a et al. (2011), but the possibility of changing the setpoint for the electricity production, instead of the full operation in the original model, is also included. All thermal power variables are average values on the HTF side. Known data include the setpoint for turbine

production (P e SP ðkÞ), the average value of maximum power available from the SF (P SFmax ðkÞ), and the TES state. As the degrees of freedom are over zero, criteria must be posed to solve the power balance model. Specifically, two criteria are used with different priority. The high-priority criterion is aimed at minimising production error (ePe ðkÞ). Once this objective is met, the low-priority criterion, which consisted in minimising the defocused thermal power in the SF (P SF def ðkÞ), is applied. Therefore, the decisions of the plant operator and control systems are supposed to be focused on reaching the setpoint for electricity production and making the most of captured solar energy, but respecting equipment limitations. The power balance model is solved by an IF–THEN algorithm, checking all options according to the input data. A first set of steady-state equations is shown next: Pe

SP ðkÞ

¼ P e ðkÞ þ ePe ðkÞ

ð3Þ

P e ðkÞ ¼ gmixed ðgexch P HTF ðkÞÞgexch P HTF ðkÞ

ð4Þ

P enet ðkÞ ¼ ggross2net P e ðkÞ

ð5Þ

gmixed ðx; kÞ ¼ p1 ðkÞgonlySF ðxÞ þ p2 ðkÞgonlyTES ðxÞ

ð6Þ

gonlyTES ðxÞ ¼ gonlySF ðxÞ  0:6=100

ð7Þ

p1 ðkÞ ¼

N loops qmloop ðkÞ N loops qmloop ðkÞ þ qTES2PB ðkÞ

ð8Þ

p2 ðkÞ ¼ 1  p1 ðkÞ

ð9Þ

gonlySF ðxÞ ¼ 0:397  0:243e qTES2PB ðkÞ ¼

H ðT out

disch

x=28:23

P d ðkÞ des Þ  H ðT in

P SFmax ðkÞ ¼ P SF ðkÞ þ P SF

def ðkÞ

ð10Þ disch des Þ

ð11Þ ð12Þ

P SF ðkÞ þ P d ðkÞ ¼ P HTF ðkÞ þ P c ðkÞ

ð13Þ

P d ðkÞP c ðkÞ ¼ 0

ð14Þ

For the sake of simplicity, Eq. (5) assumes a proportional relation between P enet ðkÞ and P e ðkÞ. Parameter ggross2net is set to 0.95, according to the rated point in Garcı´a et al. (2011) in the absence of more information on parasitic electric consumption in that paper (anyway, parasitic losses models can be found in Biencinto et al. (2014) and Wagner and Gilman (2011)). Variable x in Eq. (10) is the thermal power entering turbine in MW. Eq. (7) introduces a correction to turbine efficiency curves in the TES-only mode since HTF temperature in PB inlet is lower in this mode than in the solar-only mode. Eq. (6) calculates gmixed ð:Þ weighting gonlySF ð:Þ and gonlyTES ð:Þ by HTF mass flow rates. Eq. (13) is the power balance on HTF side. Finally, Eq. (14) does not allow simultaneous charge and discharge. The following set of constraints completes the power balance model: ePe ðkÞ P 0 P SF

def ðkÞ

P0

ð15Þ ð16Þ

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92





P HTFmin 6 P HTF ðkÞ 6 P HTFmax or

ð17Þ

P HTF ðkÞ ¼ 0 P cmin 6 P c ðkÞ 6 P cmax ðkÞ or

ð18Þ

P c ðkÞ ¼ 0 P cmax ðkÞ ¼ minðP cmax

exch ; P cmax tank ðkÞÞ

P cmax tank ðkÞ ¼ ðEmax  EðkÞÞ=ðDtp gexch
P dmin 6 P d ðkÞ 6 P dmax ðkÞ or

ð19Þ TES Þ

P d ðkÞ ¼ 0 P dmax ðkÞ ¼ minðP dmax P dmax P dmax

tank ðkÞ

exch ; P dmax tank ðkÞ; P dmax mixed ðkÞÞ

¼ gexch TES ðEðkÞ  Emin Þ=Dtp   P SF ðkÞ ðkÞ ¼ 1  P dmax onlyTES mixed P HTFmax

ð20Þ

 P HTF ðt0 Þ, which is the fraction of P HTF in ðt0 Þ used to generate electricity.  P startup ðt0 Þ, which is the fraction of P HTF in ðt0 Þ used to increase thermal state in the PB.  P HTF in ðjÞ, which is the steady-state thermal power in the PB inlet necessary to meet setpoint P e SP ðjÞ during hour j (when startup is performed). where t0 is the time since the beginning of the startup. The previous variables satisfy the following equations:

ð21Þ ð22Þ ð23Þ ð24Þ

Inequality (15) prevents turbine-generated gross electric power from exceeding its setpoint. Inequality (16) prevents the defocused thermal power in the SF from becoming negative. Inequalities (17), (18), and (21) are non-convex constraints and force null values when minimum and maximum values cannot be met. Equality (24) set a value for P dmax mixed ðkÞ linear with P SF ðkÞ. Parameters Emin and Emax are set to 0 and 1010 MW h, respectively, which is equivalent to around 7.5 h of rated turbine operation. Next, some final details are discussed. Partial SF defocus arises when SF thermal power is high enough to be transferred to the PB and TES. Moreover, defocus is useful to allow discharging power to reach its minimum value when necessary. In another vein, case P HTF ðkÞ ¼ 0 can be accomplished in two situations: (1) P e SP ðkÞ ¼ 0, or (2) P HTFmin is not reached. Then, binary variable stop to turb is set to one (see Section 4.2.4). Finally, another binary variable is added (but not shown in expressions for the sake of simplicity) to select the option of preventing charging when the turbine is stopped and the production setpoint is over zero. 4.2.6.2. Turbine startup. A model to describe turbine startup is developed in this work. This model is based on the startup model in Garcı´a et al. (2011), but certain aspects not described in this paper have been included. During the startup process, PB thermal state increases and steady-state equations are not enough to model this process. Then, a ramping behaviour of the turbine startup is modelled. Moreover, charging is assumed not to be permitted during startup and the production setpoint only allows hourly changes, as mentioned above. First, the continuoustime version of the model is described, then the model is discretised. Next, the algorithm used to guide the startup process is shown. Finally, the method to determine if enough energy is available for turbine startup is explained. Next, three continuous-time power variables and one hourly power variable are defined:  P HTF in ðt0 Þ, which is the thermal power entering the PB on the HTF side.

83

Pe

SP ðjÞ

P HTF

in ðt

¼ gonly 0

SF ðgexch P HTF in ðjÞÞgexch P HTF in ðjÞ

Þ ¼ P HTF

in ðjÞt

P HTF ðt0 Þ ¼ f ðt0 ÞP HTF P HTF

in ðt

0

0

=tstartup

in ðt

0

ð26Þ

Þ

ð27Þ

0

0

Þ ¼ P HTF ðt Þ þ P startup ðt Þ

0

ð25Þ

0

ð28Þ 0

P e ðt Þ ¼ gmixed ðgexch P HTF ðt ÞÞgexch P HTF ðt Þ Z tstartup Estartup ¼ P startup ðt0 Þdt0

ð29Þ ð30Þ

0

Notice that, for the sake of simplicity, P HTF in ðjÞ is calculated in Eq. (25) assuming the solar-only mode. Eq. (26) models P HTF in ðt0 Þ as a ramp as in Garcı´a et al. (2011). Eq. (27) applies factor f ðt0 Þ to calculate P HTF ðt0 Þ. This function ranges from 0 at the beginning to 1 at the end of startup, rising faster with time to approximate the behaviour of commercial turbines (Biencinto et al., 2014). Eq. (28) is derived from definitions. In Eq. (29), function gmixed ð:Þ is obtained with Eqs. (6)–(11), and P e ðt0 Þ is the turbinegenerated gross electric power. Finally, Eq. (30) calculates Estartup , which is the energy necessary to increase the PB thermal state during the startup process. The calculation sequence of expressions (25)–(30) in the displayed order allows obtaining all these variables from the setpoint P e SP ðjÞ. The previous model is continuous-time and must be discretised with a step of Dtp . Parameter N st is defined as the number of steps within the startup time: N st ¼ tstartup =Dtp (equal to 2 for the considered values). The discretised startup model is based on average values. Then, it can be expressed as follows, where k ¼ 1; . . . ; N st indicates the step: P HTF

in ðkÞ

ð2k  1ÞP HTF 2N st

¼

P HTF ðkÞ ¼ f avg ðkÞP HTF P HTF

in ðkÞ

in ðjÞ

in ðkÞ

¼ P HTF ðkÞ þ P startup ðkÞ

P e ðkÞ ¼ gmixed ðgexch P HTF ðkÞÞgexch P HTF ðkÞ Estartup ¼

N st X

P startup ðkÞDtp

ð31Þ ð32Þ ð33Þ ð34Þ ð35Þ

k¼1

Eq. (31) is obtained by averaging the ramp function of Eq. (26). Eq. (32) is the discretised version of Eq. (27). The same can be said for Eqs. (33)–(35) in relation to the

84

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

other equations in the continuous-time model. Notice that Eq. (34) is an approximation because function gmixed ð:Þ is not linear. Values f avg ð1Þ and f avg ð2Þ are set to 0.1 and 0.5, respectively, in the absence of information on commercial turbines. The following algorithm is used to execute the startup process, verifying there is enough energy for this process and calculating the required discharging power: for k ¼ 1 to N st do P HTF in max ðkÞ ¼ kP HTF in ðjÞ=N st P surplus ðkÞ ¼ P SFmax ðkÞ  P HTF in max ðkÞ if P surplus ðkÞ P 0 then P d ðkÞ ¼ 0 else if P dmax ðkÞ P P surplus ðkÞ then P d ðkÞ ¼ maxðP HTF in ðkÞ; P dmin Þ else P d ðkÞ ¼ 0 P HTF in ðkÞ ¼ 0 stop to turb startup ¼ 1 break end if end if end for

4.3. MILP model The MILP model is described in this subsection. Previously, the following aspects are assumed:

Notice that conditions in the algorithm are based on the maximum power required in each step, P HTF in max ðkÞ, but the average value P HTF in ðkÞ is used to calculate the discharging power, respecting minimum value P dmin . When P dmax ðkÞ is not enough, P d ðkÞ and P HTF in ðkÞ are set to null, and binary variable stop to turb startup is set to one (see Section 4.2.4). Finally, the method to determine whether enough energy is available for turbine startup is explained. Binary variable ready to turb startup is used to activate startup (see Fig. 5). This variable is obtained by Eq. (36). ready to turb startup ¼ SFstupready OR TESstupready OR mixedstupready ð36Þ

where binary variables on the right-hand side indicate that startup can be executed from the SF, from the TES and in mixed mode, respectively, and can be calculated using Eqs. (37)–(40). SFstupready ¼ Phase6 AND ðP SFmax ðkÞ P P HTF

in ðjÞÞ

ð37Þ

TESstupready ¼ ðEðkÞ  Emin Þ EHTF in ðjÞ AND ðP dmax ðkÞ gexch TES P P HTF in ðjÞÞ

P

ð38Þ

mixedstupready ¼ Phase6 AND ðEðkÞ  Emin Þgexch TES þ P SFmax ðkÞtstartup P EHTF in ðjÞ AND ðP dmax ðkÞ þ P SFmax ðkÞÞ P P HTF in ðjÞ EHTF

in ðjÞ

¼ 0:55P HTF

in ðjÞt startup

where k is the first step in the possible startup process, Phase6 (defined in SubSection 4.2.2) is evaluated at the beginning of step k, and EHTF in ðjÞ is the energy necessary in the PB inlet during the whole turbine startup process (a security factor of 1.1 is also applied). Power and energy requirements are checked during step k with Eqs. (37)–(40). However, turbine startup interruptions may occur in case solar radiation decreases during the following steps. Notice that power requirements are based on thermal power entering the PB at the end of the startup (i.e., P HTF in ðjÞ, which is the maximum value). Therefore, the power conditions required for startup are conservative. Finally, another binary variable is added (but not shown in the expressions for the sake of simplicity) to prevent startup in mixed mode.

ð39Þ ð40Þ

 This work assumes the day-ahead market and the producer’s price-taking property.  Values for hours for the delivery of the production schedule and the provision of meteorological services are based on Wittmann et al. (2011). According to this source, the production schedule for a specific day must be given during the previous day, specifically before 10.00 a.m. in Spain (note, however, that the deadline for the delivery of the production schedule in the Spanish market is 12.00 a.m. at present). Moreover, the required weather forecast is provided by the meteorological services at about 8.00 a.m. Therefore, the time available to schedule production and, then, to compute the proposed approach, is around two hours.  Market price and weather predictions are assumed to be perfect.  The price paid for electricity is supposed to be the market price (i.e., no premium per MWh is considered).  The objective to maximise is total profits from energy sales. For the sake of simplicity, production costs are not taken into account.  An optimisation horizon based on the next one or two trading days is suggested (Wittmann et al., 2011). The strategy of two-day optimisation is generally used to avoid TES emptiness at the end of the first day. In this strategy, the schedule for the first day is given to the market. Nevertheless, the one-day optimisation strategy is used in the specific example of this paper for the sake of simplicity.  The MILP model resolution is chosen to be hourly, despite the fact that lower time steps would involve a gain in accuracy. In fact, a 20 min-long time step was studied in the specific example of this paper, which did not involve a significant increase in computation time. Nevertheless, no gain in accuracy was achieved regarding hourly resolution because the setpoint in the detailed

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

model only allows hourly changes. This design decision for the simulation algorithm of the detailed model was taken to avoid frequent operation changes in the turbine. Further research in this respect is necessary. According to the proposed approach, the MILP model must be formulated to adjust to the features of detailed model as much as possible, but preserving enough simplicity to avoid increasing computational burden prohibitively. Many variables and parameters defined in Section 4.2 are now used again but with different resolution. Instead of Dtp -long resolution, the variables in the MILP model are hourly average values. The set of equations and inequalities that composes the MILP model is shown next: Maximise N X pðjÞP enet ðjÞDto  C DP e DP e abs ðjÞDto ð41Þ J¼ j¼1

subject to for j ¼ 1; . . . ; N P SFmax ðjÞ ¼ P SF ðjÞ þ P SF P SF

def ðjÞ

def ðjÞ

ð42Þ

P0

ð43Þ

P SF ðjÞ þ P d ðjÞ ¼ P HTF ðjÞ þ P c ðjÞ þ P startup ðjÞ

ð44Þ

P startup ðjÞ ¼ bðjÞP start

up

ð45Þ

P el ðjÞ ¼ a1 P HTF ðjÞ þ b1

ð46Þ

P e ðjÞ ¼ aðjÞP el ðjÞ

ð47Þ

P enet ðjÞ ¼ ggross2net P e ðjÞ

ð48Þ

Eðj þ 1Þ ¼ EðjÞ þ PcðjÞDto gexch

TES

 PdðjÞDto =gexch

TES

ð49Þ Emin 6 Eðj þ 1Þ 6 Emax

ð50Þ

aðjÞP HTFmin 6 P HTF ðjÞ 6 aðjÞP HTFmax

ð51Þ

cðjÞP cmin 6 P c ðjÞ 6 cðjÞP cmax

ð52Þ

exch

dðjÞP dmin 6 P d ðjÞ 6 dðjÞP dmax exch   P SF ðjÞ P dmax mixed l ðjÞ ¼ 1  P dmax P HTFmax P dmax

mixed ðjÞ

¼ dðjÞP dmax

P d ðjÞ 6 P dmax

ð53Þ onlyTES

mixed l ðjÞ

mixed ðjÞ

ð54Þ ð55Þ ð56Þ

cðjÞ þ dðjÞ 6 1

ð57Þ

cðjÞ þ bðjÞ 6 1

ð58Þ

bðjÞ  kðjÞ ¼ aðjÞ  aðj  1Þ

ð59Þ

bðjÞ þ kðjÞ 6 1

ð60Þ

P HTF ðjÞ 6 bðjÞDP HTF

stu max

P HTF ðjÞ P bðjÞpcDP HTF

þ ð1  bðjÞÞM

ð61Þ

stu max

ð62Þ

P e ðjÞ  P e ðj  1Þ 6 DP e

abs ðjÞ

ð63Þ

P e ðj  1Þ  P e ðjÞ 6 DP e

abs ðjÞ

ð64Þ

end for SðT min u ; T min d ; L; N 0 Þ

ð65Þ

85

where j indicates the step, j ¼ 0 expresses the instant immediately prior to the first step, N is the number of hours in studied time horizon and Dto is the time step equal to an hour. Parameter P start up (Eq. (45)) is calculated by expression P start up ¼ Estartup ðP e SP Þ=Dto , where Estartup ’s dependence on P e SP is indicated. To calculate this parameter, firstly a value for P e SP must be selected (it is unknown prior to the resolution of the optimisation problem). Secondly, an approximated value for Estartup is obtained by the detailed model. Lastly, the previous expression is used. Constraint (46) is a linear approximation of Eq. (4) where function gmixed ð:Þ is also approximated by ðgonlySF ð:Þ þ gonlyTES ð:ÞÞ=2. Parameters a1 and b1 are set to 0.3962 (dimensionless) and 3115 kW, respectively, by lineal fit. Since expression (46) is not valid when the turbine is not running, constraint (47) is added, where P el ðjÞ is a previous value only valid when the turbine is running. Constraints (49)–(56) represent the same limitations on P HTF ðjÞ; P c ðjÞ and P d ðjÞ as in the detailed model. The same method of constraints (46) and (47) is used in constraints (54) and (55) because the maximum value defined in constraint (54) is only valid when the TES is discharging. Constraints (47) and (55) are non-linear, but can be converted to linear using the method described in Appendix A, where four linear constraints replace each one. Constraints (57) and (58) avoid, respectively, simultaneous charge and discharge and simultaneous charge and turbine startup. Constraints (59) and (60) set the logical relations amongst aðjÞ; bðjÞ, and kðjÞ. Constraint (61) sets an upper bound for P HTF ðjÞ during turbine startup. Parameter M is a numeric value high enough to remove constraint (61) when bðjÞ ¼ 0. Parameter DP HTF stu max is calculated by Eq. (66). This equation calculates the maximum average value of P HTF ðtÞ during an hourly period when a startup (of N st Dtp duration) is performed and turbine is in full operation. Eqs. (31) and (32) were considered to derive Eq. (66). P  N st ð2k1Þ P HTFmax f ðkÞ Dt þ Dt  N Dt p o st p k¼1 avg 2N st DP HTF stu max ¼ Dto ð66Þ Constraint (62) establishes a minimum point when turbine starts setting a minimum value to P HTF ðjÞ equal to pcDP HTF stu max , where pc is a percentage. Variable DP e abs ðjÞ tends to the absolute value of gradient of P e ðjÞ thanks to constraints (63) and (64). Finally, SðT min u ; T min d ; L; N Þ is a set of constraints that sets minimum up time (T min u ) and minimum down time (T min d ) for the turbine. Parameters L and N 0 are the number of initial hours during which the turbine must be online and offline, respectively. Set Sð:Þ is not shown for the sake of simplicity, but is available in Dominguez et al. (2012). The other variables and constraints not described in this section are discussed in Section 4.2. Other constraints such as maximum gradients of power flows and other costs such as startup and shutdown costs are not considered but can be easily added.

86

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92 Table 3 Computation times. MILP model

Detailed model

Three-step procedure

1.5 s

12 s

25.5 s

and P SF def ðjÞ). The MILP model was accurately adjusted to the features of the detailed model, except for the turbine performance curve, its state machine, the temporal resolution and the external calculation of P SFmax ðjÞ. Notice that fine details of the turbine startup model are lost because of the hourly resolution. 4.4. Data transfer amongst models in the three-step procedure Fig. 6 shows data transfer amongst models in the threestep procedure. Variable T a ðtÞ refers to ambient temperature. Initial conditions for state variables are input data for both models. As already mentioned, the state variables are HTF temperatures in the SF (only in the detailed model), energy in the TES, thermal power in the turbine and phases in state machines. The setpoint generator (see Fig. 6) is a block that calculates values for P e SP ðjÞ, for j ¼ 1; . . . ; N from the hourly gross electric power obtained by the MILP model ðjÞ in Fig. 6). Equation P e SP ðjÞ ¼ P MILP ðjÞ is applied, (P MILP e e except in the following cases:

Fig. 6. Diagram of the three-step procedure.

5

6

x 10

5

4

3 PSFmax(j); kW

2

p(j); Euros/kWh MILP

E(j) 1

0

; kWh

MILP

10.Pe(j)

0

5

10

15

20

1. A turbine startup happens in step j. Thus, the ramp ðjÞ lower than behaviour of the turbine makes P MILP e P e SP ðjÞ. Equations similar to (46) and (66) are used to calculate the setpoint. 2. In the TES-only mode, constraint P d ðjÞ 6 P dmax onlyTES (see Section 4.3) imposes an upper bound below the rated value on the setpoint. 3. The minimum setpoint is related to P HTFmin . However, a slightly higher value is used to avoid certain problems due to differences in the turbine efficiency curves of the detailed and MILP models.

; kW

25

Hours

Fig. 7. Solution of the MILP model (instance A).

In conclusion, parameters, initial conditions for state variables, and variables P SFmax ðjÞ and pðjÞ; for j ¼ 1; . . . ; N are input data. At most, there are two independent decision continuous variables in each step j (e.g., P enet ðjÞ

After the simulation with the detailed model, the hourly average gross electric power (P e ðjÞ) is calculated. This profile (minus the parasitic electric consumption) represents the final solution, which is offered to the market as the production schedule. 4.5. Results and discussion This subsection exposes the results and discussion on the application of the proposed approach to a specific example.

Table 2 Some input data for example. Day tsunset tphase1 end %P start up Initial phases:

21/03/13 19 h 4.75 h 100% 1 and 1

T 10 T 20 T 30 T 40 T pipes0

145 °C 140 °C 135 °C 131 °C 290 °C

E0 (case 1) E0 (case 2) P HTF 0 T min u T min d

Emax =2 0 0 2h 2h

L N0 pc C DP e HTF

0 0 0.5 0.001 Euros/kW h Therminol VP-1

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

87

Table 4 Results of comparison tests when E0 ¼ Emax =2. Strategy

UCENE (MW h)

rUCENE (%)

PC (Eur)

OP (Eur)

rPI (%)

rPEE (%)

d-MILP-d d-MILP d-sMILP-d d-sMILP

0 10.34 0 28.42

0 1.36% 0 3.83%

0 155 0 426

30,906 30,751 30,386 29,959

3.16% 2.64% 1.43% 0

0 1.96% 0 6%

Table 5 Results of comparison tests when E0 ¼ 0. Strategy

UCENE (MW h)

rUCENE (%)

PC (Eur)

OP (Eur)

rPI (%)

rPEE (%)

d-MILP-d d-MILP d-sMILP-d d-sMILP

0 17.94 0 33.59

0 3.15% 0 6.09%

0 269 0 504

22,206 21,937 21,705 21,202

4.74% 3.47% 2.37% 0

0 5.29% 0 10.25%

4

x 10

5.5 5

10

4.5

9

4

Phase in SF Phase in turbine qmloop(t); Kg/s

8

10.DNI(t); kW/m

kW

3.5

2

7

3 6 2.5

MILP

Pe (j) Pe_SP(j)

2

5

*

Pe (j)

1.5

4 3

1

2

0.5 0

1 0

5

10

Hours

15

20

25 0

Fig. 8. Refinement of the MILP solution (instance A).

0

5

10

Hours

15

20

25

Fig. 10. Evolution of phases in SF and turbine.

400 350 300

ºC

250

T (t) 1

T4(t)

200

Tpipes(t) T (t) in

150

T

−AT

out_loop_des

margin

T

100

in_loop_des

Tantifreeze_ON 50 0

T

antifreeze_OFF

0

5

10

Hours

15

Fig. 9. Thermal evolution in SF.

20

25

Moreover, comparison tests with other scheduling strategies have been carried out. The selected day of study is 21-03-2013 and the selected site is Granada (Spain), the same as in the example in Garcı´a et al. (2011). Weather and electricity price forecasts are assumed to be perfect and obtained on a hourly basis. The hourly DNI forecast (DNIðjÞ) is generated from actual data from a weather station relatively near the site. The DNI profile applied to the detail model, DNIðtÞ, is generated by interpolation of DNIðjÞ, but preserving average values. The hourly T a forecast is created from TMY2 data (tmy2s, 2015), and linearly interpolated to be used in the detailed model. The electricity price forecast is obtained from data of the Iberian market operator (OMIE, Omie, 2015). Price and DNI profiles are shown in Figs. 7 and 10, respectively. Action on HTF mass flow rate and collector defocus are the means used by control systems to regulate loop output temperature. The control system based on HTF mass flow

88

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92 5

2.5

4

x 10

x 10 10

PSFmax(k) P

MILP

(k)

P

(k) , P

HTF_in

2

Pc

9

SF

(j)

*

Pc(j)

(k)

HTF

8

Pd(k) Pc(k)

7

1.5

kW

kW

6

1

5 4 3 2

0.5

1 0

0 0

5

10

Hours

15

20

25

0

5

10

15

Hours

20

25

Fig. 14. Evolution of charging thermal power in both models (instance A).

Fig. 11. Power flows on the HTF side (instance A).

4

12

x 10

P

4

(j)

12

e_SP

10

x 10

PMILP (j) d

P (k) e

P

(t)

P

(k)

HTF_in

8

P* (j)

10

d

HTF_in

P

6

startup

8

(k)

kW

kW

PHTF(k)

4

4

2

0 5.8

2 5.9

6

6.1

6.2

Hours

6.3

6.4

6.5

6.6

0

Fig. 12. Details about the startup process (instance A).

900

E(k)

800

E E

700

Emin

5

10

Hours

15

20

25

MILP

(j)

max

rate must be known to run the detailed model, but it has not been defined yet. A feed-forward control structure is used, as in Garcı´a et al. (2011). In another vein, the oneday optimisation strategy is used for the sake of simplicity. Comparison tests of the proposed approach with other three scheduling strategies have been carried out in order to study the advantages of the proposed approach. The scheduling strategies studied in these tests are the following:

600 500 400 300 200 100 0

0

Fig. 15. Evolution of discharging thermal power in both models (instance A).

1000

MWh

6

0

5

10

Hours

15

20

25

Fig. 13. Evolution of energy in TES in both models (instance A).

1. The proposed approach, with the detailed and MILP models of Sections 4.2 and 4.3, respectively (strategy d-MILP-d).

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

3. The actual production schedule, when the strategy dMILP (d-sMILP) is applied, is the schedule obtained by the strategy d-MILP-d (d-sMILP-d). Consequently, these actual production schedules are represented by variable P e ðjÞ (see Section 4.4). 4. Penalty costs only arise in strategies d-MILP and dsMILP.

4

5.5

x 10

5

MILP

Pe 4.5

(j)

*

Pe (j)

4 3.5

kW

89

3

According to the above, the following daily metrics are defined for each strategy:

2.5 2 1.5 1 0.5 0

0

5

10

Hours

15

20

25

Fig. 16. Refinement of the MILP solution (instance B).

2. A strategy based on the detailed and MILP models of Sections 4.2 and 4.3 but without the third step of the three-step procedure (strategy d-MILP). 3. The proposed approach based on the detailed model of Section 4.2 and a simpler MILP model than that of the Section 4.3 (strategy d-sMILP-d). The level of detail of this simpler MILP model is similar to that found in the literature (Sioshansi and Denholm, 2010; Usaola, 2012; Dominguez et al., 2012; Pousinho et al., 2015). 4. A strategy based on the detailed model of Section 4.2 and the simple MILP model, and without the third step of the three-step procedure (strategy d-sMILP).

1. UCENE: unmet committed net electric energy. 2. rUCENE: relative unmet committed net electric energy compared to the net electric energy that is effectively served. 3. PC: penalty cost. 4. OP: obtained profit. 5. rPI: relative profit increase regarding the strategy d-MILPs (representative of the strategies in the literastrategyÞprofitðd-sMILPÞ ture). That is, rPI ¼ profitðevaluated profitðd-sMILPÞ 6. rPEE: relative profit estimation error.

The simple MILP model of the strategies d-sMILP-d and d-sMILP is the same as the MILP model of Section 4.3, except for the elimination of nonzero lower bounds in constraints (52) and (53), and the elimination of constraints (54)–(56). These constraints are not usual in literature. Therefore, the use of a more complex MILP model and the three-step procedure are the elements of the proposed approach that have been evaluated in theses comparison tests. For the definition of metrics used for the evaluation of strategies, the following aspects are considered:

Note that profit estimation error arises because of penalty costs and excessively optimistic profit estimations. Moreover, the penalty cost is calculated assuming 15 Euros per MW h of deviation (an upper bound of penalty cost in the study of Kraas et al. (2013)). Two cases are studied in the comparison tests: (1) E0 ¼ Emax =2 and (2) E0 ¼ 0. Table 2 shows all other input data for the example. Parameter P start up is calculated at rated point and, therefore, under the most demanding circumstances. The initial phases are Nocturnal recirculation I and Turbine shutdown. Parameter C DP e is set to a low, near-zero level to penalise changes in P e ðjÞ, but with little influence on costs. The MILP model was generated in YALMIP (Lo¨fberg, 2004), a free toolbox for MATLAB, and solved by SCIP (v3.0.2) (Achterberg, 2009), a free solver for mixed integer problems. The simulation algorithm of the detailed model was developed in MATLAB R2011a and executed with an Intel Core i5-4440 processor. Table 3 shows the approximated required computation times. Next, Sections 4.5.1 and 4.5.2 expose the comparison results and the analysis of solutions, respectively.

1. Weather forecast and the detailed model are supposed to be perfect. Therefore, production schedules obtained by strategies d-MILP-d and d-sMILP-d can be applied faithfully. 2. Strategies d-MILP and d-sMILP are less restrictive than d-MILP-d and d-sMILP-d, respectively. Consequently, their production schedules (P MILP ðjÞ, see Section 4.4) e and their estimated profits are too optimistic. Therefore, these production schedules cannot be applied truly.

4.5.1. Comparison results with other scheduling strategies Tables 4 and 5 show the comparison results amongst strategies for the two studied cases. As expected, the strategy d-MILP-d is the best for both cases. It can be seen that both the improvement in the accuracy of the MILP model and the three-step procedure yield better results for both cases. The improvements are higher in the E0 ¼ 0 case. This fact is explained in Section 4.5.2. The greatest improvement of the strategy d-MILP-d is achieved compared to the

90

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

strategy d-sMILP and for the E0 ¼ 0 case, where metric rPI is 4.74% against 0 and metric rPEE is 0 (weather forecast and the detailed model are supposed to be perfect, as mentioned above) against 10.25%. Note that the approach developed in this paper can achieve even better results than those exposed in this example in the following situations: 1. A comparison with a scheduling strategy that does not take into account the actual plant’s operation philosophy. Note that the four strategies use the detailed model developed in this paper, which includes the plant’s operation philosophy. 2. The use of a more precise detailed model than the detailed model used in this paper. This way, this greater complexity could prevent a good adjustment of the MILP model to the detailed model (contrary to the good fit achieved in the studied example as seen in Section 4.5.2). This difference between models would give more importance to the three-step procedure and, then, the gain in accuracy would be higher.

4.5.2. Analysis of solutions Figs. 7–15 show some results from the application of strategies d-MILP-d and d-MILP to the E0 ¼ Emax =2 case (instance A), and Fig. 16 compares production schedules of strategies d-sMILP-d and d-sMILP in the E0 ¼ 0 case (instance B). The average values shown in figures are located at the centre of its time interval. Fig. 7 shows the solution of the MILP model. Thanks to the TES, the production is displaced towards high-price periods. The two startups are undertaken in TES-only mode to bring forward the generation. During low-price periods, electricity production is kept at low levels or interrupted from 2.00 h to 6.00 h. It can be seen that in the TESonly mode, upper bound P dmax onlyTES limits production to a point below the rated point. Moreover, the TES is empty at the end of the optimisation horizon, and, therefore, the TES energy is used in full. Fig. 8 shows the refinement of the MILP solution by the last simulation run on the detailed model (the final solution is also shown in hourly average values). It can be seen that the correction is small and focused on startup periods. Notice that the time resolutions of the constraints of the MILP and detailed models are hourly and 10 min, respectively. Thus, the constraints of the detailed model are more stringent, whilst MILP constraints are applied to hourly average values. This difference is more evident in time periods with large variations, such as startup processes. Corrections of the setpoint in the TES-only mode (see Section 4.4) are not shown in Fig. 8 but performed internally in the detailed model. The evolution of HTF temperatures in the SF and operating phases can be observed in Figs. 9 and 10, respectively. Variable qmloop ðtÞ is also shown in Fig. 10. The HTF does

not circulate in the SF until around 4.45 a.m., and HTF temperatures cool down accordingly. Temperatures in loops (T 1 ðtÞ to T 4 ðtÞ) fall much faster than the HTF temperature in the insulated pipes (T pipes ðtÞ) since the receiver tubes are not insulated. At around 4.45 a.m., the HTF starts to circulate through the SF bypassing the PB and at a rate of 1 kg/s per loop. As a consequence, temperatures tend to be equal. When the sun rises at around 7.0– 8.0 a.m., the mass flow rate increases to 2.5 kg/s per loop and the HTF warms thanks to solar radiation. At around 8.50 a.m., the loop outlet temperature (T 4 ðtÞ) starts to be regulated by the temperate controller, which maintains T 4 ðtÞ around its setpoint of 390 °C in case of enough radiation. At around 18 p.m., there is not enough radiation to follow the setpoint and the temperatures fall. At around 19.5 p.m., T pipes ðtÞ is too low and the HTF circulation through the SF stops, cooling down the HTF temperatures. Power flows on the HTF side is shown in Fig. 11, which displays the following sequence of power operating modes: TES ? PB, STOP, TES ? PB, SF + TES ? PB, SF ? PB + TES, SF ? PB, SF + TES ? PB, TES ? PB, STOP. Power limitations during discharge in mixed and TES-only modes and three time periods with partial defocus can also be seen in Fig. 11. The first two periods arise because the upper bound for charging thermal power is reached, whilst the lower bound is reached in the third period. Details on the startup process are shown in Fig. 12 (see Section 4.2.6.2), where variable P HTF inðkÞ is only defined during startup, otherwise it is zero. It can be seen in Fig. 12 that turbine-generated electric gross power (P e ðkÞ) ranges from 0 at the beginning to its setpoint at the end of startup, rising faster with time to approximate the behaviour of commercial turbines. Differences between some MILP output variables (strategy d-MILP) and final output variables (strategy d-MILP-d) are portrayed in Figs. 13–15. The MILP solution is quite close to the final solution as shown in the figures. The good fit of the MILP model to the detailed model can explain this fact. Finally, a comparison between production schedules of strategies d-sMILP-d and d-sMILP in the case E0 ¼ 0 is shown in Fig. 16 (instance B). This is the instance where the difference between production schedules is higher. There are two reasons for this higher difference: 1. During the startup period, details about power flows and about SF operating phases are lost in the MILP model. In instance B, the startup occurs in solar-only mode (because E0 ¼ 0) and the MILP model brings forward the startup regarding the actual situation due to the aforementioned lack of constraint details. This aspect is less critical in TES-only startups, as they are independent from SF operating phases and only the details on power flows are lost.

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

2. Constraints about upper bounds for discharging thermal power in TES-only and mixed modes are removed in the simple MILP model. Then. only the upper bound due to the heat-exchanger limitations (higher than previous bounds) is applied.

5. Conclusion A novel approach is proposed for optimal scheduling problems in CSP plants. The approach is based on using two models: a detailed model and a MILP model. This way, MILP capabilities and the accuracy of a detailed model are combined. It is suggested that the detailed model include the actual plant’s operation philosophy and has a time step below one hour. Moreover, the MILP model must be formulated to adjust as much as possible to the features of the detailed model (e.g., to the equipment’s performance constraints). A three-step procedure is proposed to solve the problem. The detailed model is used in the first and third steps, whilst the MILP model is used to solve the optimisation problem itself in the second step. The first step calculates the thermal production in the SF. The last step applies the whole complexity of the detailed model to the MILP solution, thus bringing it closer to reality and allowing certain level of simplicity in the MILP model, necessary for computation reasons. The proposed approach is applied to a 50 MW PTCbased CSP plant with molten-salt-based TES. The detailed model used is an improved version of a model available in literature. This model was validated by its authors using real operating data, and has been improved in this paper to be used in optimal scheduling problems. The improvements regarding other scheduling strategies for a specific example are shown. The context of this work is deterministic (i.e., market price and DNI predictions are assumed to be perfect). The effect of inaccurate forecast is currently under research. Other future lines include: (1) improvement of the detailed model used and the proposed three-step procedure by including the effect of collector defocus in SF thermal evolution and (2) application of the proposed approach to an example with a fossil fuel backup system that assists production. Appendix A. Converting nonlinear constraints to linear Non-linear constraint (A.1) is equivalent to the four linear constraints (A.2)–(A.5), where zðjÞ and xðjÞ are real variables. dðjÞ is a binary variable, and parameters M and N are the upper and lower bounds of xðjÞ, respectively. zðjÞ ¼ dðjÞxðjÞ

ðA:1Þ

zðjÞ 6 MdðjÞ

ðA:2Þ

zðjÞ P mdðjÞ

ðA:3Þ

91

zðjÞ 6 xðjÞ  mð1  dðjÞÞ

ðA:4Þ

zðjÞ P xðjÞ  Mð1  dðjÞÞ

ðA:5Þ

References Achterberg, T., 2009. Scip: solving constraint integer programs. Math. Program. Comput. 1 (1), 1–41, . Alliance C, 2014. The Economic and Reliability Benefits of CSP with Thermal Energy Storage: Literature Review and Research Needs. Technical Report. Bemporad, A., 2003. Modeling, Control, and Reachability Analysis of Discrete Time Hybrid Systems. University of Siena, DISC School, Dept. Inf. Eng., Italy (Lecture Notes of the DISC Course). Biencinto, M., Bayn, R., Rojas, E., Gonzlez, L., 2014. Simulation and assessment of operation strategies for solar thermal power plants with a thermocline storage tank. Solar Energy 103 (0), 456–472. http://dx. doi.org/10.1016/j.solener.2014.02.037, . Channon, S., Eames, P., 2014. The cost of balancing a parabolic trough concentrated solar power plant in the spanish electricity spot markets. Solar Energy 110 (0), 83–95. http://dx.doi.org/10.1016/j.solener.2014.08.036, . Dinter, F., Gonzalez, D.M., 2014. Operability, reliability and economic benefits of {CSP} with thermal energy storage: first year of operation of {ANDASOL} 3. Energy Proc. 49 (0), 2472–2481. http://dx.doi.org/ 10.1016/j.egypro.2014.03.262, (proceedings of the SolarPACES 2013 International Conference). Dominguez, R., Baringo, L., Conejo, A., 2012. Optimal offering strategy for a concentrating solar power plant. Appl. Energy 98 (0), 316–325. http://dx.doi.org/10.1016/j.apenergy.2012.03.043, . ´ lvarez, J.L., Blanco, D., 2011. Performance model for Garcı´a, I.L., A parabolic trough solar thermal power plants with thermal storage: comparison to operating plant data. Solar Energy 85 (10), 2443–2460. http://dx.doi.org/10.1016/j.solener.2011.07.002, . Ho, C., 2008. Software and Codes for Analysis of Concentrating Solar Power Technologies, Report SAND2008-8053. SANDIA National Laboratories, USA. Kost, C., Flath, C.M., Mst, D., 2013. Concentrating solar power plant investment and operation decisions under different price and support mechanisms. Energy Policy 61 (0), 238–248. http://dx.doi.org/10.1016/ j.enpol.2013.05.040, . Kraas, B., Schroedter-Homscheidt, M., Madlener, R., 2013. Economic merits of a state-of-the-art concentrating solar power forecasting system for participation in the spanish electricity market. Solar Energy 93 (0), 244–255. http://dx.doi.org/10.1016/j.solener.2013.04.012, . Kuravi, S., Trahan, J., Goswami, D.Y., Rahman, M.M., Stefanakos, E. K., 2013. Thermal energy storage technologies and systems for concentrating solar power plants. Prog. Energy Combust. Sci. 39 (4), 285–319. http://dx.doi.org/10.1016/j.pecs.2013.02.001, . Law, E.W., Prasad, A.A., Kay, M., Taylor, R.A., 2014. Direct normal irradiance forecasting and its application to concentrated solar thermal output forecasting a review. Solar Energy 108 (0), 287–307. http://dx. doi.org/10.1016/j.solener.2014.07.008, . Lizarraga-Garcia, E., Ghobeity, A., Totten, M., Mitsos, A., 2013. Optimal operation of a solar-thermal power plant with energy storage and electricity buy-back from grid. Energy 51 (0), 61–70. http://dx.doi.org/

92

M.J. Vasallo, J.M. Bravo / Solar Energy 126 (2016) 73–92

10.1016/j.energy.2013.01.024, . Lo¨fberg, J., 2004. Yalmip: a toolbox for modeling and optimization in MATLAB. In: Proceedings of the CACSD Conference. Taipei, Taiwan . Madaeni, S., Sioshansi, R., Denholm, P., 2012. How thermal energy storage enhances the economic viability of concentrating solar power. Proc. IEEE 100 (2), 335–347. Omie, 2015 . Powell, K.M., Hedengren, J.D., Edgar, T.F., 2014. Dynamic optimization of a hybrid solar thermal and fossil fuel system. Solar Energy 108 (0), 210–218. http://dx.doi.org/10.1016/j.solener.2014.07.004, . Purohit, I., Purohit, P., Shekhar, S., 2013. Evaluating the potential of concentrating solar power generation in northwestern india. Energy Policy 62 (0), 157–175. http://dx.doi.org/10.1016/j.enpol.2013.06.069, . Reddy, V.S., Kaushik, S., Ranjan, K., Tyagi, S., 2013. State-of-the-art of solar thermal power plants a review. Renew. Sustain. Energy Rev. 27 (0), 258–273. http://dx.doi.org/10.1016/j.rser.2013.06.037, . The SAM website, 2015 (last access: 06.04.15).

Simoglou, C.K., Biskas, P.N., Bakirtzis, A.G., 2012. Optimal selfscheduling of a dominant power company in electricity markets. Int. J. Electr. Power Energy Syst. 43 (1), 640–649. http://dx.doi.org/ 10.1016/j.ijepes.2012.05.038, . Sioshansi, R., Denholm, P., 2010. The value of concentrating solar power and thermal energy storage. IEEE Trans. Sustain. Energy 1 (3), 173– 183. http://dx.doi.org/10.1109/TSTE.2010.2052078. User’s manual for tmy2s, 2015 (last access: 22.04.15). Usaola, J., 2012. Operation of concentrating solar power plants with storage in spot electricity markets. Renew. Power Gener., IET 6 (1), 59–66. http://dx.doi.org/10.1049/iet-rpg.2011.0178. Wagner, M.J., Gilman, P., 2011. Technical Manual for the SAM Physical Trough Model, Technical Report. National Renewable Energy Laboratoy of USA. Wittmann, M., Eck, M., Pitz-Paal, R., Mller-Steinhagen, H., 2011. Methodology for optimized operation strategies of solar thermal power plants with integrated heat storage. Solar Energy 85 (4), 653– 659. http://dx.doi.org/10.1016/j.solener.2010.11.024, (solarPACES 2009). Zhang, H., Baeyens, J., Degrve, J., Cacres, G., 2013. Concentrated solar power plants: review and design methodology. Renew. Sustain. Energy Rev. 22 (0), 466–481. http://dx.doi.org/10.1016/j.rser.2013.01.032, .