Power Plant Operation Optimization: Economic Dispatch of Combined Cycles with Machine Learning

1| Power Plant Operation Optimization Economic dispatch of combined cycle power plants STEFANO ROSSO SUPERVISOR: MONIKA TOPEL CAPRILES EXAMINER: BJÖRN LAUMERT, EDUARDO PRIETO |2 Abstract As electricity production from renewable sources increases, higher flexibility is required by fossil fuel generation to cope with the inherent fluctuations of solar and wind power. This results in shorter operating cycles and steeper ramps for the turbines, and more uncertainty for the operators. This thesis work applies mathematical optimization and statistical learning to improve the economic dispatch of a combined cycle power plant composed by two separate blocks of two gas turbines and one steam turbine. The goal is to minimize the input fuel to the gas turbines while respecting a series of constraints related to the demand the plant faces, power generation limits etc. This is achieved through the creation of a mathematical model of the plant that regulates how the plant can operate. The model is then optimized to reduce fuel consumption at a minimum. Machine learning techniques have been applied to sensor data from the plant itself to realistically simulate the behavior of the turbines. Input-Output curves have been obtained for power and exhaust heat generation of all the turbines using ordinary least squares on monthly data with a ten minutes sampling rate. The model is cross-validated and proven statistically valid. The optimization problem is formulated through generalized disjunctive programming in the form of a mixed-integer linear problem (MILP) and solved using a branch-and-bound algorithm. The output of the model is a one-week dispatch, in fifteen minutes intervals, carried out for two months in total. Lower fuel consumption is achieved using the optimization model, with a weekly reduction of fuel consumed in the range of 2-4%. A sensitivity analysis and a correlation matrix are used to highlights the demand and the maximum available capacity as critical parameters. Results show that the most efficient machines (alternatively, the ones with highest available capacity) should be operated at maximum load while still striving for an efficient utilization of the exhaust gas. 3| Contents Sammanfattning .............................................. Error! Bookmark not defined. Abstract .................................................................................................................. 2 1. Introduction ..................................................................................................... 5 1.1 Objectives ....................................................................................................... 6 1.2 Scope and limitations ............................................................................... 7 1.3 Thesis Structure ........................................................................................ 8 2. Economic dispatch in literature ...................................................................... 9 3. Optimization ................................................................................................... 13 3.1 Mathematical Optimization ......................................................................... 13 3.2 Generalized Disjunctive Programming .................................................... 17 3.3 Branch-and-Bound ................................................................................... 19 3.4 Indicator Variables ................................................................................... 21 3.5 Linearization ............................................................................................ 23 4. Statistical learning ......................................................................................... 25 4.1 Linear Regression........................................................................................ 25 4.2 Simple linear regression .......................................................................... 26 4.3 Fitness of the model ................................................................................ 28 4.4 Evaluation metrics and cross-validation ................................................. 31 5. Optimization Model....................................................................................... 34 5.1 Correlation Matrix ....................................................................................... 34 5.2 Additional Parameters ................................................................................ 36 5.3 Optimization Model .................................................................................... 37 6. Results & Discussion ...................................................................................... 41 6.1 Block 1 ........................................................................................................... 41 6.2 Block 2 ......................................................................................................... 44 6.3 Full Plant ..................................................................................................... 46 6.4 Discussion ................................................................................................... 48 |4 7. Conclusions .................................................................................................... 52 References ............................................................................................................ 53 Appendix A – Pyomo Examples ........................................................................... 57 Example 1 - A Linear Programming (LP) model .............................................. 57 Example 2- Production Scheduling, linear formulation .................................. 58 Example 3 - Capacity Expansion Planning ..................................................... 59 5| Introduction 1. Introduction An increasing magnitude of renewable capacity is being installed year after year. This is causing major shifts in the pattern of electricity consumption and generation by source. In the next years, fossil fuel generation is forecasted to be more and more important as a back-up producer, i.e. to compensate for renewables fluctuation and unpredictability, and ensure stability in the grid. The implications are shorter operating cycles, steeper ramps for the turbines, and more uncertainty for the operators. Combined cycle power plants (CCPP) in particular are the most suited to cover this stabilizing role given the high ramping capacity of the gas turbines, which allows them to start-up fast and reach high loads in a matter of minutes, and their relatively low environmental impact. In comparison with coal-fueled power plant, CCPP produce less emissions thanks to a cleaner fuel, while also achieving greater overall efficiency. In an attempt to understand the effects that these changes will have on the life cycle of the turbines, on their operation patterns and how to prepare the operators to face sudden change in demand, Siemens Industrial Turbomachinery (SIT) would like to put at better use the vast amount of data collected from the hundreds of sensors installed in its machines all over the world. Several projects have been set up in the past years to carry investigation about machines lifetime - studying components lifetime and how they are affected by various kinds of failures, - and operations. Operation planning problems are characterized by minimization or maximization of some objective functions (e.g. fuel cost, profit, machines downtime, emissions) given a combination of input parameters and several constraints. Determining the optimal power output under certain conditions is at the core of the economic dispatch (ED) problem. The extension of the ED to multiple timestamps is the unit commitment (UC), which usually provides an hourly overview of the load to be provided by the generators. The classic way of solving these problems is through physical models only, which means that thermodynamic equations are used, with some level of approximation Introduction |6 to allow the optimization process to finish in a reasonable time. While this method generally works, it doesn’t take advantage of the individual characteristics of the single generation units. Case-specific solutions have been usually investigated and, while more general software able to handle different problem sets exists already in the market, a definitive answer is yet to come. This thesis project aims to integrate the physical way of modelling with the sensor data from the plant to reach an improved solution, builds on the previous work carried out in Bahilo Rodriguez (2018) to create an optimization model for gas turbines operation and it’s situated within a bigger and longer strategy with the final objective of expanding the model to the CCPP while providing short-term and long-term profit increase to the customers. The work presented in this thesis will be implemented with the Pyomo optimization package in Python, which presents some advantages over other well-established algebraic modeling languages like AMPL and GAMS. First and foremost, Pyomo is open source, making it a great candidate for scientific studies and thesis projects; being based on and fully integrated within the Python language it is easy to declare and manipulate objects, integrate different modules and create a clean pipeline from the raw data to the final problem. Pyomo interfaces to open-source and commercial solvers like GLPK and CPLEX, providing all the degrees of freedom and customization that these softwares allow for (changing solving algorithm, optimization parameters and so on). 1.1 Objectives The main goal of this thesis work is to analyze fuel consumption patterns from a power plant and recognize opportunity for improvement that could lead to fuel savings while providing the same load to the grid. This process usually goes under the name of descriptive analytics. The power plant being used as a case study comprises two separate combined cycle units, where each unit is composed by two gas turbines (each with its own recovery boiler) and one steam turbine. To reach this goal, the author has developed an optimization model that simulates plant operation under certain constraints. It should be noticed that the model will serve as the basis for other two master thesis projects that include profit optimization as objective function (Ahmed, 2019) and maintenance constraints (Bhatt, 2019). In this respect, the model should be: 7| Introduction • Simple and fast. Being the main building block, it is important for it to not act as a bottleneck, to require too many parameters or to render other applications impractical or difficult to implement. • Accurate. Simplicity should not be a priority over achieving good results, in any case. The model will need to achieve significant improvement on the status quo in order to justify its deployment. • Able to generalize. While the model is designed on a specific case study, the results should be reproducible on a different plant with minimal changes. The output of the model gives the optimal power output for each turbine at a minimum cost, which is a solution to the ED problem. 1.2 Scope and limitations Dispatch optimization is a multi-faceted problem with multiple outcomes on many different levels, to make for an interesting problem that is easy to expand indefinitely. Just as an example, the duration of the optimization period will allow for different objectives such as fuel consumption minimization in the short term, long and short-term profit maximization, availability (up-time) maximization, and asset management. It is then especially important to put some clear boundaries to what will be achieved and what not. Data for a duration of two years was collected. Out of these, a total period of two months was chosen for the investigation. This is due to a general lack of good data in the first eighteen months and a low demand in the last four that often resulted in only one turbine being operated. The analysis would have given the same results for these months. The frequency chosen for the optimization of the dispatch is fifteen minutes, since it has been considered the best trade-off between simplicity (the model doesn’t need too many constraints) and accuracy (the model is not as vague as a one-hour dispatch). This influences what needs to be coded into the model. Regarding the gas turbines, the author of the thesis will: • Analyze sensor data, trying to reconstruct missing values and create new information from it (e.g. deriving the input fuel from the power output and the efficiency). This is known as data wrangling. • Build input-output (IO) curves. These curves map the relationship between two variables like the power output and the efficiency or the power output and the exhaust heat. Introduction |8 Regarding the steam turbines, the author of the thesis will: • Create a start-up model from the sensor data and the start-up curves. Since the start-up of a turbine depends on how long the turbine has been shutdown and can thus happen in a range that varies between few minutes (hot start-up) to hours (cold start-up), this is not a trivial task. • Investigate the relationship between the exhaust heat of the gas turbines and the steam turbines power output. This was complicated by the lack of data from the recovery boilers. As already briefly mentioned in the objectives paragraph, other two thesis projects were developed in parallel to this one and with very clear demarcation lines among them. Any consideration related to increased thermal stress, reduced lifetime, wear of the components etc. was deferred, as well as any economic consideration regarding bidding, fuel price, and costs in general. Finally, environmental considerations regarding maximum allowed emissions, the possibility of introducing a cap and trade scheme or a carbon tax were impossible to formulate due to general lack of information regarding emissions of different kind (e.g. NOx, CO2). While the general goal of this work, reducing fuel consumption, will lead to lower emissions for same equipment facing the same load, the real impact cannot be properly estimated (Rifaat, 1998). 1.3 Thesis Structure • Chapter 2 presents a literature review of the ED problem, showing how different scholars have approached it and what are the most common techniques in use. • Chapter 3 explains more in detail what an optimization problem is and how it can be formulated and solved. Multiple examples will be provided, using the economic dispatch as a reference. • Chapter 4 introduces some concepts of machine learning and statistical methods. This chapter covers only the basis of the subject, in order to show how the sensor data helps the optimization problem. • Chapter 5 is dedicated to the model structure. The chapter deals with all the additional details of building a real optimization model. • Chapter 6 shows the results obtained and various sensitivity analysis. • Chapter 7 lays down some conclusions and basis for future work. 9| Economic dispatch in literature 2. Economic dispatch in literature There is an abundance of literature for ED problems in CCPP targeting the problem from many perspectives. One of the main differentiators is how the researchers approach the realization of a coherent, realistic and solvable model. Compared to the simple cycle of traditional thermal units, the number of possible interactions among the CC components and the variety of constraints at the level of the single turbine and the whole plant create a much harder problem to solve. The need for simplification leads to the use of different models and, according to how the plant is represented, it is possible to group them into four kinds (Bayon et al., 2014): • Aggregated models. The CC is represented as a single thermal unit and different constraints regarding units’ interactions and limits are ignored. This highly simplistic model leaves the determination of the unit commitment to the dispatch operator. • Pseudo-unit models. Pseudo-units are used to represent the gas turbines and their proportional share of steam generation (Lopez et al.,2010) with all the units facing the same constraints. This could be considered a slightly more complex aggregated model. • Configuration or mode models. The CC is represented as mutually exclusive combinations of gas turbines and steam turbines. The transition to one state to the other is regulated by a predetermined state transition diagram. • Physical models. Each of the CC component is represented by different constraints and parameter, with its own ramping capacity, maximum power production etc. These last two methods are considered the most realistic, with the configuration model being adopted by a vast majority of the researchers. One of the main advantages of the configuration model is the apparent simplicity of the idea, as each configuration is close to a pseudo-unit facing its own set of constraints. In comparison, the physical model requires a bigger set of constraints for each unit and a higher level of detail. Optimization |16 This expression (3.9) is non-linear, and its associated problem will be much harder to solve than its linear counterpart. It is not in the scope of these thesis to demonstrate so, but more information can be found in Hochbaum (2007) and Sorenson (1976). If the efficiency is non-linear and the fuel input relies on it, how can we formulate a linear objective function for our problem? The answer to this question will be introduced in the next sections, after talking about linear regression. For the moment, it is sufficient to say that we want to minimize the fuel input for each turbine i at each interval t given a certain power output P (3.10). 𝑚𝑖𝑛 ∑∑𝐹𝑢𝑒𝑙𝐶𝑜𝑠𝑡𝑖,𝑡 𝑁 𝑖=1 𝑇 𝑡=1 (𝑃𝑖,𝑡) 3.10 What we need to define now are some general constraints the model will be subject to. As said before, we want to meet a certain demand, so our power output should match that demand. For each turbine the power output should stay Figure 3.1 – Block Diagram of a combined cycle (Gicquel, 2011) 17| Optimization between a lower and an upper limit, and at each instant there is a maximum power output that can be obtained depending on the ambient condition that will affect various parameters e.g. the compressor inlet temperature (CIT). Moreover, there is a ramping rate that should not be exceeded by the turbine, so that the difference in power output for any two subsequent instants cannot be higher than a certain value. These constrains are represented by equations 3.11 to 3.15: ∑𝑃𝑖,𝑡 𝑁 𝑖=1 =𝐷𝑖 3.11 𝑃𝑖,𝑡≤𝑃𝑖,𝑡≤𝑃𝑖,𝑡     3.12 𝑃𝑖,𝑡     ≤𝛼𝑃𝑚𝑎𝑥,𝑖,𝑡(𝐶𝐼𝑇𝑖,𝑡,𝑅𝐻𝑖,𝑡,𝑝𝑖,𝑡) 3.13 𝑃𝑖,𝑡− 𝑃𝑖,𝑡−1≤𝑅𝑈𝑖 3.14 𝑃𝑖,𝑡−1−𝑃𝑖,𝑡≤𝑅𝑈𝑖 3.15 The meaning of 𝛼 in 3.13 will be discussed in the data quality subsection. Other generally valid constraints are slightly harder to express: the turbines have different start-up times according to the duration of the shut-down, there is a minimum time the turbine must be functioning after the start-up and a minimum time it must be off after the shut-down, if a turbine is off eq. 3.12 doesn’t apply, the steam turbine can only be on if any of the gas turbines is, and many more. To model these constrains we need to introduce indicator variables 𝑦𝑖,𝑡 which are 1 when the turbine is on and 0 when it’s off, time related variables ℎ𝑜𝑛,𝑖,𝑡, ℎ𝑜𝑓𝑓,𝑖,𝑡 that are used to count for how long the turbine is been on/off, and possibly other extra variables. It is easy to see that while power can assume any continuous value within its constraints, having 𝑦=0.3 has no meaning. This is an example of a MILP model. An interesting and simple way to write down MILPs is generalized disjunctive programming (GDP). 3.2 Generalized Disjunctive Programming Disjunctive programming is a major subfield of mathematical optimization where at least one constraint must be satisfied, but not necessarily all the others. Scheduling problems can often be written in a GDP form. It must be noted that this is simply an higher-level formulation of a MILP problem. We first show the formulation of a GDP problem (Raman and Grossmann, 1994), then define its terms: min𝑍=𝑓(𝑥)+ ∑𝑐𝑘 𝐾 𝑘 Optimization |18 𝑠.𝑡.𝑔(𝑥)≤0 ⋁ 𝑌𝑖𝑘 𝑟𝑖𝑘(𝑥)≤0 𝑐𝑘=𝛾𝑖𝑘 𝑖∈𝐷𝑘 𝑘∈𝐾 𝛺(𝑌)=𝑇𝑟𝑢𝑒 𝑥≤𝑥≤𝑥 𝑥∈𝑅𝑛,𝑐𝑘∈𝑅1,𝑌𝑖𝑘∈{𝑇𝑟𝑢𝑒,𝐹𝑎𝑙𝑠𝑒} First, Z is our objective function. The function g represents the set of global constraints. Each k denotes a disjunction composed by 𝐷𝑘terms connected by an OR. In each term is present a Boolean variable 𝑌𝑖𝑘 that determine whether the inequality and equality constraints ( 𝑟𝑖𝑘(𝑥)≤0,𝑐𝑘=𝛾𝑖𝑘) are enforced. Each 𝑌𝑖𝑘 variable is contained in 𝛺(𝑌), a set of logic proposition used to enforce various logic (in most of the applications this is an exclusive OR operator with one of the disjunct necessarily active). For more details on these terms see Sawaya (2006). Example 3.3 - GDP model of a gas turbine Suppose we want to model a gas turbine with two states: on and off. There is no reference to what the objective function of the problem is or how the turbine interacts with other component of the model, so that only the two disjuncts are going to be shown. Off state (𝑌1): 𝑃𝑖,𝑡=0 3.16 𝑦𝑖,𝑡=0 3.17 ℎ𝑜𝑛,𝑖,𝑡=0 3.18 ℎ𝑜𝑓𝑓,𝑖,𝑡= ℎ𝑜𝑓𝑓,𝑖,𝑡−1+1 3.19 𝑢𝑝𝑖𝑦𝑖,𝑡−1−ℎ𝑜𝑛,𝑖,𝑡−1≤0 3.20 On state (𝑌2): 𝑃𝑖,𝑡≤𝑃𝑖,𝑡≤𝑃𝑖,𝑡     3.21 𝑃𝑖,𝑡−𝑃𝑖,𝑡−1≤𝑟𝑎𝑚𝑝𝑖,𝑡 3.22 𝑃𝑖,𝑡−1−𝑃𝑖,𝑡≤𝑟𝑎𝑚𝑝𝑖,𝑡 3.23 𝑦𝑖,𝑡=1 3.24 ℎ𝑜𝑛,𝑖,𝑡=ℎ𝑜𝑛,𝑖,𝑡−1+1 3.25 19| Optimization ℎ𝑜𝑓𝑓,𝑖,𝑡= 0 3.26 𝑑𝑤𝑖−𝑑𝑤𝑖𝑦𝑖,𝑡−1−ℎ𝑜𝑓𝑓,𝑖,𝑡−1≤0 3.27 Disjunction: 𝑌1+𝑌2 = 1 3.28 Equation 3.20 and 3.27 enforce the minimum up-time (𝑢𝑝𝑖) and down-time (𝑑𝑤𝑖) of the machine If the reader wants to know more about how this constraint was obtained, they can jump to the Indicator Variables paragraph. Equation 3.28 forces one and only one of the two sets of constraints to be active. In order to implement this in practice, a big-M method is used (again to be found in the Indicator Variables paragraph). 3.3 Branch-and-Bound We have seen how an integer programming problem is in its essence an LP problem with the additional constraint of some (or all) of the variables being integer. The first step to solve any IP problem is to first solve its associated LP problem (from now on called relaxation or relaxed problem) where the condition of having integer variables is removed. If all the necessary variables end up being integer in the first optimal solution the problem is solved, but this is not usually the case and additional steps are needed. To illustrate how the Branch-and-Bound algorithm works, consider this IP: maximize 18𝑥1+11𝑥2 subject to 𝑥1+𝑥2 ≤7 13𝑥1+7𝑥2≤65 𝑥1,𝑥2 ≥0 𝑥1,𝑥2 𝑖𝑛𝑡𝑒𝑔𝑒𝑟𝑠 The solution to the relaxed LP (which will be noted with 𝑃0) gives 𝑥1=83,𝑥2=13 3, with an objective function ζ = 96.67. ζ here represents our upper limit, since adding constraints cannot increase the value of the objective function. For 𝑥1 to be an integer, it must satisfy either 𝑥1≤2 or 𝑥1≥4. We add these two constraints, separately, to create other two problems (𝑃1 𝑎𝑛𝑑 𝑃2). Optimization |20 The solution for 𝑃1 gives 𝑥1=2,𝑥2=5 and ζ = 91. This is our lower limit ζ and it can only be updated if a new integer solution with a higher objective function value is found. The solution for 𝑃2 gives 𝑥1=3,𝑥2=4.71 and ζ = 95.86. This solution updates our upper limit and creates two new problem ( 𝑃3 𝑎𝑛𝑑 𝑃4). If the new upper limit would have been lower than our lower limit, we could have declared to have found the optimal solution, since 𝑃3 ,𝑃4 and all the other resulting problems would necessarily have a lower optimum value than 𝑃2. Fig. 3.2 shows how we are starting to develop a tree of LPs called enumeration tree (Vanderbei, 1996). 𝑃1 is a leaf of this tree and is noted with a double box, while 𝑃2 is called a node of the tree. Figure 3.2 – The beginning of the enumeration tree Solving 𝑃3 gives us additional branches and the next step could be to investigate the deepest node in the tree 𝑃5(Depth-First Search, DFS) or the earliest found unsolved node 𝑃4(Breadth-First Search, DFS). The algorithm generally follows a DFS logic since it has been proven to be optimal (Edelkamp & Schrodl, 2012). 21| Optimization The solution to 𝑃5 𝑥1=3,𝑥2=3 has an objective function value ζ = 87 which is less than our current lower limit and our best solution is not updated. Fig. X.2 shows the state of the enumeration tree after solving 𝑃6. Eventually, the path following 𝑃7 brings to two integers solution with a lower objective function than ζ, 𝑃8 is infeasible and so 𝑃4. The enumeration tree has been exhausted and 𝑃1 is proven to be our optimal solution. 3.4 Indicator Variables It is common practice to introduce in LPs binary variables that when coupled with other continuous variables are used to indicate certain states or conditions. Let’s take the case in which we want to identify a turbine as on if its power is higher than zero and off otherwise. We can introduce an indicator variable δ in a constraint of the form: 𝑃− Mδ≤0 3.29 Where M is an upper bound for the power 𝑃. Equation 3.29 enforces the condition ‘if 𝑃>0 then δ=1’. It is not possible to completely represent the condition ‘δ=1 if and only if 𝑃>0’ unless we define a threshold ε above which the turbine will be considered to be on. This is usually not a limitation when representing real life problems, since we generally don’t want to distinguish between e.g. having no impurities in a fuel or having three molecules of impurities. Thus, we can write: 𝑃− εδ≥0 3.30 Successfully representing the condition ‘𝑖𝑓 δ=1 then 𝑃≥ε’. 3.29 and 3.30 give us the condition ‘δ=1 if and only if 𝑃≥ε’. We have seen how to use indicator variables to represent states, but they can also be used to represent conditional logic statements. In a GDP formulation of the MILP, each disjunct has an associated binary variable that determines whether to apply the associated constraints or not. Optimization |22 Figure 3.3 – The enumeration tree after solving 𝑷𝟔 23| Optimization Let us consider a generic inequality: ∑𝑎𝑖𝑥𝑖≤𝑏 𝑖 3.31 We want to use an indicator variable δ that will cause the constraint to be enforced when δ=1 and to ignore the constraint when δ=0. We divide this task in two parts: first, we model the condition ‘if δ=1 then the constraint is enforced’. This can be done considering that what we are trying to express is equivalent to ‘∑𝑎𝑖𝑥𝑖−𝑏≤0 𝑖 if (1−δ)=0′, that is: ∑𝑎𝑖𝑥𝑖−𝑏≤𝑀(1−δ) 𝑖 3.32 ∑𝑎𝑖𝑥𝑖+𝑀δ≤𝑀+𝑏 𝑖 3.33 To model the second condition ‘if δ=0 then constraint is not enforced’ we need once again to introduce a tolerance ε, so to have: ∑𝑎𝑖𝑥𝑖− 𝑏>0 𝑖 3.34 ∑𝑎𝑖𝑥𝑖≥ε+b 𝑖 3.35 ∑𝑎𝑖𝑥𝑖−(𝑚−ε)δ≥ε+b 𝑖 3.36 Where m is a lower bound for the expression ∑𝑎𝑖𝑥𝑖−𝑏 𝑖. Equations 3.33 and 3.36 together ensure the expected behavior. 3.5 Linearization We have already introduced, without making it explicit, an example of linearization. Eq. 3.22 and 3.23 are effectively the linear version of the condition we wanted to achieve: |𝑃𝑖,𝑡− 𝑃𝑖,𝑡−1|≤𝑅𝑈𝑖 3.37 In general, it is not possible to linearize any given function, but it is possible to approximate them to piecewise linear function, in other words functions that are linear for a certain interval. Fig. 3.4 shows how log(x) is approximated to straight lines between points (1,5,10,20,35,50). This results in the following equations: X = 1𝜆1+5𝜆2+10𝜆3+20𝜆4+35𝜆5+50𝜆6 3.38 y = 0𝜆1+1.609𝜆2+3.302𝜆3+3.995𝜆4+4.555𝜆5+4.912𝜆6 3.39 𝜆1+𝜆2+𝜆3+𝜆4+𝜆5+𝜆6=1 3.40 Additionally, for these equations to be valid, at most two adjacent variables 𝜆𝑖 can be non-zero. This condition is defined as 𝜆𝑖 being part of a special ordered set of type 2 (SOS2) and while it cannot be expressed in linear terms, it can be modelled using binary variables and integer programming (Beale and Tomlin, 1969). Optimization |24 Figure 3.4 – Piecewise approximation of the logarithmic function We have introduced the basic concepts of linear programming, giving an idea of many different ways to build mathematical optimization models, the limitations of each of them and some ways to overcome them. In particular, it has been discussed how any sort of logical consideration can be introduced into an LP transforming it into an IP (using indicator variables) and how this affect solving time (by having to solve a slightly different version of the same problem multiple times). A way to transform any function in a linear fashion has been provided. 25| Statistical learning 4. Statistical learning This paragraph is dedicated to a machine learning technique called linear regression, what it is and why it is needed in this particular problem. The mathematical foundations are laid down together with some considerations for creating a statistically valid model. 4.1 Linear Regression From our knowledge of mathematical optimization, it would be possible to build a thermodynamical model (TDM) of the plant, which equations have been known for years, taking into account some of the specific characteristics of the turbine models (nominal capacity, ramping capacity and so on) to customize it for different sites. This model would be, in theory, the most accurate possible thanks to its sound theoretical foundations. We are going to argue here that there are two shortcomings that make this model impractical to use and less accurate than expected. The first problem regards solving time. In the TDM each turbine’s component has its inputs and outputs (mass flow, heat, power) and the relationship with each other component is mapped, so that for each load the turbine is facing, the equations are solved again and again to ensure the best behavior that is able to satisfy all the constraints. If we were to analyze only one gas turbine, this problem could be solved in a somehow acceptable time since the demand faced by the turbine would be identical to the plant demand and the problem would be only solved once, but when introducing (at least) another gas turbine and (at least) a steam turbine, the possible combinations become easily too many to be solved in reasonable time even once. From what we know about the branch and bound algorithm, in reality the problem would be solved thousands of times before being able to declare optimality of the solution. The second concern is about accuracy of the model over time and its difficulties in dealing with plant degradation. The TDM relies on a plethora of coefficients usually determined experimentally and directly related to the conditions of the plant. With time, the components start degrading and some of their properties change. A pure thermodynamical approach wouldn’t be suited to model this changes (which are statistical in nature) and it could result not only in nonoptimal dispatches but even in non-feasible ones. The penalties for not being able Statistical learning |32 Table 3 – Single outlier Index Error |Error| Error^2 1 0 0 0 2 0 0 0 3 0 0 0 4 0 0 0 5 0 0 0 6 18 18 225 Table 4 – Results of different magnitudes of error in RMSE and MAE Table 1 Table 2 Table 3 RMSE 3 4.266 7.348 MAE 3 3 3 This can be considered an advantage of the RMSE when in fact larger errors are detrimental to the model performance in a non-linear way (i.e. an error twice as big as another one is more than twice as bad). There is anyway an important consideration to be done here: the MAE is the lower limit of RMSE, while its upper limit is given by the MAE times the square root of the sample size 𝑛 (Chai and Draxler, 2014). The tendency of the RMSE to increase with test sample size makes it hard to compare among different models. Once one (or a series of) evaluation metric is chosen it can be used to improve underfitting or overfitting problems. Generally speaking, underfitting is a condition in which the model is not good enough to produce clear patterns so that using more data, data with higher quality, a more powerful regression technique or introducing more features in the model will likely solve the problem. Overfitting is more interesting, since many algorithms tend to be skewed in this direction. Training a model on a set of data and evaluating the results on the same set has little value, as shown before. The model would be not as good when presented with new data, sometimes drastically worse. Common practice is to divide the set in two partitions called test and validation set. The parameters are learned fitting the training set and the model is scored on the validation set. While an improvement on the previous situation, it is still possible to overfit on the validation set, especially when using more complex techniques than linear 33| Statistical learning regression. This make sense, since we are trying to optimize for it and repeated use of the validation set will indirectly give information to the model. A possible solution is to use thee partitions instead of two. The new partition, called the test set, is never shown to the model until all the parameters are learned, and it is thus possible to have an unbiased estimation of our error. The drawback is that doing so reduces the amount of data we can use in the training phase, decreasing the overall quality of the model. Another possible solution is cross-validation, a procedure that splits the training data into smaller sets and repeatedly evaluate the model fitting on different combinations of these sets. The way the splitting is done determines the exact cross-validation method used. K-Fold validation, for example, splits the set in equal parts with no resampling and uses K-1 folds to train the model and one fold to evaluate prediction errors. Finally, the model is tested on the real test set. Figure 4.5 shows how this process looks like: Figure 4.5 – K-Fold Validation Other cross-validation techniques are leave-one-out (LOO) validation, - which is an extreme case of K-Fold where K = N and the model is trained with the whole training set except one sample, - and Shuffle & Split, which generates random folds of equal length with resampling. Optimization Model |34 5. Optimization Model This section presents how the model is built from the selection of the decision variables and parameters to the structure of the optimization. 5.1 Correlation Matrix As mentioned before, our model makes little use of thermodynamic parameters and constraints. To map the relationships between variables and the objective function, we build linear regressions among them, but a question stands still: how to select these variables? One approach could be to just use them all since the model will automatically recognize the ones that have a higher influence (their associated weights will be higher) and the ones that do not (whose weights will be closer to zero). This is a disingenuous solution for two reasons: the first one is that the model will likely overfit, the second that having to track many parameters will slow down optimization. A better approach is to use a correlation matrix, using Pearson correlation coefficients to determine the degree of correlation among variables (Pearson, 1895). We will not delve into the theory behind this particular technique, it suffices to say that a correlation matrix expresses the correlation among two variables using a number between -1 (perfectly negative correlation) and 1 (perfectly positive correlation), with 0 implying that no correlation exists. Figure 5.1 shows the results of running a correlation analysis on data from one gas turbine (plus the efficiency of the whole cycle). The high correlation between efficiency, active load, heat input and exhaust heat content should not surprise since they all depend on the same hidden variable (turbine inlet temperature). It is not possible to use two of these variables together when performing regression on another variable. If we want to minimize fuel consumption, our dependent variable will be the heat input (‘Heat In’). If we want to maximize the efficiency of the combined cycle, that (‘Eff CC’) will be the variable. 35| Optimization Model Figure 5.1 – Correlation matrix for one gas turbine In both cases, it is advised to select the active load as the independent variable, the reasons being: 1. This is the variable we are trying to optimize for in our original problem (i.e. dispatch optimization); 2. The other parameters are usually calculated ex post and there is no way for the plant operator to decide them in the dispatch; Optimization Model |36 3. It is not immediately clear how the model itself would use the exhaust heat content or the efficiency of the turbine as decision variables. Also, compressor inlet temperature (CIT) and relative humidity (RH) will be used, since they are easy to add to the model (they do not depend on the results in any way) and show a discrete negative correlation to our independent variables. Another parameter that will be obtained through linear regression is the exhaust heat content. The same variables will be used in this case. 5.2 Additional Parameters There are two additional considerations to be added before creating our regression curves. They both concern the maximum output of our turbines. Siemens is using an algorithm to calculate the thermal load percentage of its machines and it is be possible to calculate the maximum available capacity (𝑃) from it at any instant, since 𝐿𝑜𝑎𝑑%= 𝑃𝑖 𝑃𝑖 5.1 Unfortunately, for this thesis’ case study, the thermal load percentage doesn’t always match the electric one due to unexpected customer behavior. To overcome this problem, another evaluation of the maximum available power must be done, using a different technique that relies on calculations made over some specific conditions over eleven variables. It was not possible to collect all the needed values for all the parameters, - and even if it would have, the problem of component degradation (thus, altering the behavior of the machines) would still be present. In other words, we have an algorithm providing correct values for some known conditions and another algorithm providing values that must be corrected to take into account all the relevant factors. These two algorithms should provide the same results, so we are able to calculate the correction factor 𝛼. 𝑃1=𝛼𝑃2 5.2 Similarly, the maximum output of the steam turbine is controlled by the exhaust heat content and the effectiveness of the HRSG. No data is present for this last component, so that we had to consider together the efficiency of the boiler and 37| Optimization Model the steam turbine in a single parameter 𝛽. With the assumption that the steam turbine is always producing at maximum power for a given heat input (a conservative but realistic assumption), we have: 𝑃𝑠𝑡𝑒𝑎𝑚=𝛽𝑄𝑣󰇗 5.3 The value of 𝛽 varies highly in the dataset, which is a problem mainly because affect whether the problem will solve (i.e. whether it will be possible to satisfy the demand). Four different approaches were tested, from most to least realistic: • Using the median of 𝛽. This is the most conservative approach, but it allows the optimization to fail. • Using the median of 𝛽 introducing a variable for the unmet demand. Assumes a fixed relationship between the exhaust energy content and the steam turbine output and allow the system to miss part of the load in order to make the problem always solvable. • Trying to predict 𝛽 introducing a variable for the unmet demand. This is the most complex option and it also relies more heavily on the assumption that the steam turbine output is equivalent to its maximum capacity. • Linear regression on 𝛽𝑄𝑣󰇗 as a unique parameter. This doesn’t work because the model is able to find ways to maximize steam output beyond what is reasonable expected in order to save fuel. The second solution proved to be the best when working with a configuration model, while the unit model never fails even when using simply the median of 𝛽. 5.3 Optimization Model Once the variables have been selected it is possible to use linear regression to obtain the input-output curves for each turbine. Optimization Model |38 Figure 5.2 – Linear regression on inlet fuel Figure 5.3 – Linear regression on exhaust heat content 39| Optimization Model These curves map the power output to our variables for given ambient conditions and are fed to the model. It is possible to see, from the shape of the regression curves, that the square of the power was used together with the aforementioned variables. Figure 5.4 shows the unit model optimization layout, how the turbines are connected and what the inputs and outputs are. The equations governing the units are given in chapter 3. In the next chapter we are going to show the results of the model and how changing some of the parameters (namely, demand and maximum available capacity) can produce different outcomes. Optimization Model |40 Figure 5.4 – Layout of the optimization model 41| Results & Discussion 6. Results & Discussion This section introduces the results of the optimization problem. The economic dispatch in a period of eight weeks is obtained for the two individual blocks. The algorithm is also validated for one-week period considering the whole plant. To further investigate the validity of the algorithm, the demand faced by the plant is artificially increased or decreased. An additional load profile is investigated. The related results are discussed in the last section of this chapter. 6.1 Block 1 In Fig. 6.1 is possible to see how the optimization results compare against the real operations of the first combined cycle block. One week has been chosen as the most representative, similar plots can be obtained selecting any time period. Two things can be immediately noticed: 1. The orange and blue line in the upper graph almost exactly overlap. This means that the customer always operates the two gas turbines at the same load. 2. The steam turbine is operated in both cases at very similar levels. This is an indicator that the algorithm is able to realistically reproduce the interaction between the gas and the steam cycle. Regarding the optimization results, we are concerned with them being sensible, realistic and to not present any immediate and clear drawback. Results & Discussion |48 Figure 6.9 - Comparison of fuel consumption, full plant 6.4 Discussion In order to verify that the results obtained are not overfitting the data (as explained in chapter 3), the trained model is used to optimize different loads. These loads are obtained multiplying the existing load by a factor of 0.6, 0.8, 1.1 and 1.15 (sensitivity analysis), plus an extra load obtained using the demand faced by another plant being investigated in the wider scope of this project (model fitness analysis). While there are no reference results for these additional loads, it is possible to consider whether the results are sensible or not in comparison with the reference state (multiplication factor of 1) The results of the sensitivity analysis are shown in Fig. 6.10. Figure 6.10 – Sensitivity Analysis with multiple multiplication factors 49| Results & Discussion As expected, the best results are obtained when the load approaches the maximum available capacity of the combined cycle (factor of 1.15). This allows for both of the turbines to run at their best efficiencies and minimize fuel consumption in respect to the power output. The low efficiencies obtained by lowering the demand by 0.8 was also expected and can be easily explained by having to run the most efficient turbine at a load lower than the maximum available capacity or having to sacrifice some of the steam turbine load, both options necessarily decrease the efficiency of the plant. This case also presents the wider swings in efficiency, a characteristic due to the load being around the point in which one of the turbines could be shut down. Figure 6.11 – Production patterns for 0.6 and 0.8 multiplication factors An interesting result is that multiplying the demand by 0.6 and 1.1 doesn´t move them much from the reference case. For the 0.6 case this is due to the plant being able to operate without the inefficient turbine most of the time (when this is not true anymore, the efficiency gets closer to the most inefficient case) while for the multiplication factor of 1.1 it seems to simply depend on the load level itself and it may be related to the inaccuracies of the regression (both in fuel consumption and steam production). Results & Discussion |50 The results of figure 6.12 look familiar at this point, with the only difference of the load decreasing so much that one of the turbines shuts-down. Another noticeable difference is that GT01 approaches often its lower limit of 20% load, which enforces the steam turbine to regulate its power output according to the demand. Thus, the efficiency profile becomes really similar to the shape of the steam turbine power output (considering all other factors equal, a better utilization of the exhaust heat naturally increases the efficiency). Figure 6.12 – One week optimum dispatch, alternative load Finally, figure 6.13 represents the variation in efficiency with the variation in maximum capacity. As expected, when increasing the maximum power output of the turbines, the efficiency of the whole cycle increase. The opposite is true when reducing the maximum capacity. Overall, the model behaves as expected when varying the load to different degrees and when presented with a notably different load profile, showing it to be accurate and general enough. The behavior in comparison to the maximum output is the expected one. 51| Results & Discussion Figure 6.13 – Maximum available capacity sensitivity analysis |52 7. Conclusions As a result of the work presented in these pages, is possible to conclude that: • There is significant space for improvement in operation of CCPP in respect to the fuel consumption of the plant; • Using few data sources is possible to create a realistic model that can provide an optimum economic dispatch; • Adding complexity to the model can bring even higher savings at the cost of increased computational power needed. • In particular, being able to model the HRSG and the chiller unit are clear priorities to increase confidence in the model results. • It is very important to determine correctly demand faced by the turbines, as it is one of the main variables in the resulting plant efficiency. Future work can be done starting from this thesis, in particular: • The lack of data for the heat recovery steam generation (HRSG) boilers did not allow to work on combined heat and power (CHP) plant since too many assumptions were needed for the model to be considered realistic. Future work should consider cogeneration plants. • Chillers and additional equipment can be included in the model to improve the results of the optimization. • If possible, analysis should be done on the carbon content of the fuel used, in order to allow for environmental considerations to be made. 53| References References Ahmed, M. E. (2019). Short-term profit operation optimization of combined cycle power plants using Machine Learning (Master thesis), unpublished. Anscombe, F.J. (1973). Graphs in Statistical Analysis. The American Statistician. Bahilo Rodriguez, E. (2018). Unit commitment of gas turbines using Machine learning and MILP programming (Master thesis), unpublished. Bayon, L., Garcia Nieto, P.J., Grau, J.M., and Suarez, P.M. (2014). An economic dispatch algorithm of combined cycle units. International Journal of Computer Mathematics. Beale, E.M.L. and J.A. Tomlin (1969). Special facilities in a general mathematical programming system for non-convex problems using ordered sets of variables, in J. Lawrence (Ed.), Proc. 5th Int. Conf. of Operations Research, Tavistock, London. Bellman, R. (1957). Dynamic Programming, Princeton University Press. Bhatt, D. (2019). Short-term profit operation optimization of combined cycle power plants using Machine Learning (Master thesis), unpublished. Byelogrlic, M.R. (2000). Inclusion of Combined Cycle Plants into Optimal Resource Scheduling. Power Engineering Society Summer Meeting. Carlisle, A., Dozier, G. (2001). An Off-The-Shelf PSO. Proceedings of the Particle Swarm Optimization Workshop. Chai, T., and Draxler, R.R. (2014). Root mean square error (RMSE) or mean absolute error (MAE)? – Arguments against avoiding RMSE in the literature. Geosci. Model Dev. Conejo, A.J., Castillo, E., Garcia-Berstrand, R., and Minguez, R. (2006). Decomposition Techniques in Mathematical Programming: Engineering and Science Applications. Springer. Correa, S.M. (1993). A Review of NOx Formation Under Gas-Turbine Combustion Conditions. Combustion Science and Technology. References |54 Floudas, C.A., and Pardalos, P.M. (2009). Encyclopedia of optimization. Springer. Gaing, Z.L. (2003). Particle swarm optimization to solving the economic dispatch considering the generator constraints. IEEE Transactions on Power Systems. Gicquel, R. (2011). Energy Systems - A New Approach to Engineering Thermodynamics. CRC Press. Edelkamp, S., and Schrodl, S. (2012). Heuristic Search - Theory and Applications. Academic Press. Hochbaum, D.S. (2007). Complexity and algorithms for nonlinear optimization problems. Springer. Hui, H., Yu, C., Gao, F., and Surendran, R. (2011). Combined Cycle Resource Scheduling in ERCOT Nodal Market. IEEE Power and Energy society General Meeting. Kim, B.M., Kim, Y.B., and Oh, C.H. (1997). A study on the convergence of genetic algorithms. Computers & Industrial Engineering. Kirk, D.E. (1970). Optimal Control Theory: An Introduction. Englewood Cliffs, NJ: Prentice-Hall. Konash, O., El-Sharakawi, M. (2009). Economic Dispatch using Particle Swarm Optimization for combined cycle generators. IEEE/PES Power Systems Conference and Exposition. Liu C., Shahidehpour, M., Li, Z., and Fotuhi-Firuzabad, M. (2009). Component and Mode Models for the Short-Term Scheduling of Combined-Cycle Units. IEEE Transactions on Power Systems. Lopez, J.A., Gomez, R.N., and Moya, I.G. (2010). Practical Commitment of Combined Cycle Plants using Dynamic Programming. IEEE Electrical Power & Energy Conference. Michaud, M.G., Westmoreland, P.R., Feitelberg, A.S. (1992). Chemical mechanisms of NOx formation for gas turbine conditions. Symposium (international) on combustion. 55| References Othman, F., Rahmani, R., Yusof, R. (2012). Solving economic dispatch problem using particle swarm optimization by an evolutionary technique for initializing particles. Journal of Theoretical and Applied Information Technology. Pancholi, R.K., Swatup, K.S. (2003). Particle swarm optimization for economic dispatch with line flow and voltage constraints. TENCON 2004. Conference on Convergent Technologies for Asia-Pacific Region. Pearson, K. (20 June 1895). Notes on regression and inheritance in the case of two parents. Proceedings of the Royal Society of London. Raffensperger. J.F. (1999). The marriage of dynamic programming and integer programming. Proceedings of the ORSNZ 24th Annual Conference. Raman, R. and Grossmann, I.E. (1994). Modelling and Computational Techniques for Logic Based Integer Programming, Computers & Chemical Engineering. Rayaprolu, K. (2009). Boilers for power and process. Boca Raton: CRC Press. Rifaat, R.M. (1998) . Economic dispatch of combined cycle cogeneration plants with environmental constraints. IEEE Proceedings of the 1998 International Conference on Energy and Management and Power Delivery. Sawaya, N. (2006). Reformulations, Relaxations and Cutting Planes for Linear Generalized Disjunctive Programming. PhD thesis, Chemical Engineering Department, Carnegie Mellon University. Sorenson, H.W. (1976). An introduction to nonlinear programming—I: Necessary and sufficient conditions. Computer & Electrical engineering. Talbi, E.G. (2009). Metaheuristics: From Design to Implementation. Wiley. Vanderbei, R.J. (1996). Linear Programming: Foundations and Extensions. Springer. Yalcinoz, T.,Altun, H., Uzam, M. (2001). Economic Dispatch Solution Using A Genetic Algorithm Based on Arithmetic Crossover. IEEE Porto Power Tech conference. References |56 Yousif, M., Ai, Q.,Gao, Y.,Watoo, W.A., Jiang, Z. and Hao, R. (2018). Application of Particle Swarm Optimization to a Scheduling Strategy for Microgrids Coupled with Natural Gas Networks. Energies. 57| Appendix A – Pyomo Examples Appendix A – Pyomo Examples Example 1 - A Linear Programming (LP) model from pyomo.environ import * # First we create the model to which everything will be assigned model = ConcreteModel() # Variables declaration model.x = Var([0,1,2,3], within = NonNegativeReals) # Constraint declaration model.c1 = Constraint(expr = 3*model.x[0]+2*model.x[1]+4*model.x[2]+3*model.x[3] <= 96) model.c2 = Constraint(expr = 2*model.x[2]-model.x[0] <= 0) model.c3 = Constraint(expr = 3*model.x[3]+model.x[2]-model.x[1]<= 0) model.c4 = Constraint(expr = model.x[3] <= 5) #Objective function model.obj = Objective(sense = maximize,expr = 120*model.x[0]+85*model.x[1]+200*model.x[2]+160*model.x[3]) # Calling the solver solver = SolverFactory('cplex') solver.solve(model)