HydroPowerModels.jl: A Julia/JuMP Package for Hydrothermal Economic Dispatch Optimization

HydroPowerModels.jl is a Julia package for solving multistage, steady-state, hydro-dominated, power network optimization problems with stochastic dual dynamic programming (SDDP). Our state-of-the-art open source tool is flexible enough for practitioners in the electrical sector to test new ideas in an efficient way. This tool was made possible by the Julia language and the surrounding ecosystem of packages. We use JuMP, a package for mathematical programming modeling; PowerModels.jl, a JuMP-extension for power network optimization; and SDDP.jl, another extension that implements the SDDP algorithm.


Introduction
The hydrothermal dispatch problem is important for the planning and operation of hydro-dominated electrical systems such as the Brazilian national grid. The objective is to coordinate generation, energy distribution, and hydro-storage management in order to minimize the cost of operation. In this context, the hydrothermal dispatch problem is a medium-term planning problem where uncertainties related to the hydrology (i.e., rainfall and other inflows) of the hydroelectric plants and consumer-demand have fundamental importance. Generation takes two main forms: i) hydrogeneration with a low marginal cost, and ii) thermal-generation, with a high marginal cost. However, in hydro-dominated systems such as Brazil, there is often insufficient capacity of thermal generation to meet demand. Thus, the system-operator faces a trade-off between using water for cheap generation in the present, against conserving water for future periods of drought. Because of this trade-off, the hydrothermal dispatch problem is often modeled as a multistage stochastic problem [16,14,17]. Solving multistage stochastic programs, however, is a challenging numerical problem. The solution of the problem is intractable and common approximations such as stochastic dynamic programming [2] [3] can suffer from high dimensionality (frequently referred as the curse of dimensionality). One method that partially overcomes the curse of dimensionality is the stochastic dual dynamic programming algorithm of [16]. A number of programs implementing SDDP are in-use around the world, ranging from unpublished implementations in academic institutions to professional software such as the product developed by PSR, a software and consulting company, also called SDDP 1 [18]. However, until recently, there was no fast, reliable, and open-source implementation of the SDDP algorithm. Without such a tool, researchers and practitioners have not had a common ground for the discussion and analysis of different hydrothermal dispatch formulations and their solutions. The objective of this work is to build an open source tool, called Hy-droPowerModels.jl, that can be this common ground. HydroPow-erModels.jl can be used to assess the impact of modeling choices during the planning of a hydrothermal power system. These choices include the usage of different network formulations, the consideration of different risk measures, and the planning horizons for uncertain future costs. Addressing these issues provides the research community and the energy industry with a powerful tool for the efficient design of hydrothermal power systems. To develop a tool that can be used by both researchers and industry professionals, we take advantage of the Julia language [4] and two main packages: PowerModels.jl [7], which implements power flow models for electrical dispatch, SDDP.jl [8], which implements the stochastic dual dynamic programming algorithm. Both Power-Models.jl and SDDP.jl handle their respective optimization models through JuMP.jl [9], a Julia package for mathematical optimization. JuMP.jl makes it simple to write complex optimization problems and solve them with numerous solvers. HydroPowerModels.jl takes advantage from the fact that PowerModels.jl and SDDP.jl were not only developed in Julia, but also deeply rely on JuMP.jl to build and solve mathematical optimization problems. SDDP.jl is used to specify the hydro storage dynamics and stochastics of inflows, renewables and loads. PowerModels.jl is used to provide multiple network dispatch formulations as a starting point for the HydroPow-erModels.jl formulation that couples the electrical constraints with hydro constraints and uncertainty. The next sections are organized as follows: 1) A brief understanding of multistage stochastic problems and decisions under uncertainty; 2) A presentation of the modeling formulation and problem specification; 3) Different simplifications of the model; 4) An explanation of the solution method (SDDP); 5) A comparison of existing alternative solution implementations and listing of the proposed package functionalities; 6) A case study to clarify usage.

Multistage Stochastic Optimization
Stochastic Programming (SP) [20] is a branch of optimization under uncertainty, where the realization of some random variables (ω) influence the conditions of the problem and consequently the optimal decisions. Uncertainty relates to the probability distributions of parameters and may be incorporated to the problem in various manners. One class of simple stochastic programs is known as a two-stage problem with recourse. It can be formulated as follows: The objective of the first stage is composed of an immediate cost term c T x and a cost-to-go function Q(x). Q is a function of some decision variables of the first stage x that fix the state of the second stage (its feasible region). These variables are called state variables. First stage decisions are made under uncertainty, while the decision variables of the second stage u are chosen after the realization of the variable ω (under a deterministic scope) and are so called recourse variables. The function ρ is known as a risk measure [1] and is commonly assumed to be the expectation operator E.
The goal of this problem is to find an optimal stage decision x and an optimal second-stage decision u for each realization of ω conditioned on x. Collectively, this set of decisions is known as a policy.
The two-stage problem discussed above naturally extends to a multistage problem via recursion. A multistage stochastic program with T stages can be formulated as follows: Assuming the problem is linear, we have: In this setting, the uncertain data ω 1 , ..., ω T is revealed gradually over time. The sequence ω t ∈ R d t of data vectors is viewed as a stochastic process, i.e., as a sequence of random variables with a specified probability distribution. Just like in the two stage problem, there are state variables at every stage that partially impact the objective through the risk measure of the subsequent stages optimal values. The policy optimized by this problem is a mapping from the realized stochastic process to the decisions for each stage. Using dynamic programming [2], the nested formulation of a multistage stochastic program may be represented by the following Bellman recursion for each stage: For the purpose of simplicity, we assume that Q T +1 (·, ·) = 0. In these equations, the optimal value at stage t depends on previous decision x t−1 and the realization of the data process ω t . Finally, the optimal value of the first stage problem gives the optimal value of the corresponding multistage problem.

The Hydrothermal Dispatch Problem
The hydrothermal dispatch problem is a multistage stochastic optimization problem that comprises an optimal power flow (OPF) problem and the hydro storage management for multiple periods and scenarios. Introduced by [6], the OPF problem extends the economic dispatch problem, where the goal is to plan the operation of an electrical power system by determining the contribution of available energy sources in supplying demands, to include constraints representing the power flow equations, resulting in a more realistic model of the dispatch. Since its introduction, the OPF problem has received additional constraints to better represent power systems. The adopted version of the PowerModels.jl package, is a AC Optimal Power Flow (AC-OPF) problem. The AC Power Flow constraints implement voltage bounds, generation bounds, nodal conservation of power, power flow on lines and thermal limit of the lines (power flow limit) for an alternating current system. Different from original OPF formulations, support was also added for multiple load and shunt components on each bus together with a line charging that supports conductance and asymmetrical values. Besides the modifications from PowerModels.jl, deficit variables to work as slack to the attendance of each load have frequently been used by the electrical sector and will be incorporated in our package. This variables allow a more detailed evaluation of the cost of not attending some demands. The optimal dispatch tries to find the most economic use of system resources and is the objective function of an OPF. However, in a hydrothermal power system, where water is a main resource and its inflow is uncertain, risk averse planning is an important task to ensure the lowest cost operation. The necessary hydro storage management to attend demand during periods of scarcity adds an extra layer of decisions and constraints. Reservoir operation is a problem composed by the balance of incoming and outgoing flow of water as well as its usage to generate power. Since many storage facilities are linked through rivers, this problem may be simply viewed as a directed flow graph where each node in the graph represents a storage reservoir. A constraint binds the power generation of generators using water as fuel to the out-flow of the respective storage. Operation is also restricted by the limits of reservoir volume and outflow. The solution of this problem returns an operation policy representing the optimal generation and hydro management possible given the horizon studied and the scenarios considered. As it is common for multistage problems, we will define the underlying sub-problem (i.e., Q) that unifies the OPF problem with the hydro-thermal dispatch problem.

The Mathematical Model
Following the notation chosen by the PowerModels.jl package, the sets and parameters used to define the sub-problem are listed with the addition of those created for the hydro storage management and the data for the state of the sub-problem:

Sets.
N -buses R -reference buses E, E R -branches, forward and reverse orientation G, G i -generators and generators at bus i ∈ N L, L i -loads and loads at bus i ∈ N S, S i -shunts and shunts at bus i ∈ N H, H G -reservoirs and reservoirs with power generation Operators.
-real part of a complex number ∠ -angle of the polar representation of a complex number (·) * -complex conjugate | · | -absolute value States.
For simplicity the stage index t of the sub-problem is omitted from every data and variable with exception of the incoming reservoir volume ν h,t−1 (in order to differentiate it from the outgoing reservoir volume ν h ). Notice that, as usual, alternating components (i.e. voltage, power generation, power demand and power flow) are modeled using complex numbers to fully represent the process. Bellow, the decision variables of our problem: Variables. A complete mathematical model of the sub-problem is as follows, subject to: The objective of the sub-problem (1a) is to minimize the costs of real power generation, cost of deficit from the real energy supply, cost of spillage (in order to try and avoid degenerate solutions) and the cost-to-go function Q. Constraint (1b) fixes reference buses complex voltage angles to zero, as the remaining angles will be defined accordingly. Constraint (1j) bounds the complex power generation, representing the physical limitation of generators and fuel source availability. The magnitude of the complex voltage is bounded in constraint (1k) by restricting the absolute square of its value. The upper limit alone defines a circular feasible region for each voltage, while the lower limit reshapes the region as a ring, bringing a non-convexity to the problem. The Branch complex power flow is formulated in (1d) and (1e), dependent on the voltage at each end and implementing elements of line charging and the effects of transformers. The power flow is bounded in (1f) through its absolute value. These power limits of the lines represent thermal limits and stability limits.
Constraint (1c) implements Kirchhoff's Current law (KCL), which refers to power preservation at each node, balancing generation, demand, flow and shunt. Although, deficit variables have been added in order to guarantee feasibility in case of lack of power availability. Angle difference between buses are bounded in (1g). The reason for the limits is to approximate the transient stability constraints of power flowing in branches. These restrictions refer to the synchronism among machines at each end of a line. The limits depend on the equipment installed and the system configuration. An important variable in an economic dispatch problem is the marginal cost of energy at each bus, which, in optimality, is determined by the dual value of (1c). This value is also referred to as a shadow price, local marginal price (LMP), or nodal price. Regardless, this value represents the cost of an extra unit of energy in a bus. The conservation of water equation is implemented in (1h), where the water stored at a reservoir should equal the water previously stored plus the incoming flow (precipitation and water from upstream reservoirs), minus the portion used to generate energy and the one spilled away. The binding of the hydro real power generation and water used is done in (1i), which depends on the efficiency of the generator modeled through a production factor. Constraints (1l) and (1m) bound respectively the volume of water stored and the amount of water used in generation. These limits are defined by the capacity limit of storage facility and the equipment installed.
Note that problem (1) is a nonlinear, non-convex optimization problem. This will have important implications for our solution approach.

Network Formulations
The AC Optimal Power Flow problem, defined in (1), captures the nonlinear and complex nature of the power flow. However, solving a nonlinear problem (NLP) is hard. Guaranteeing global optimality often involves an extensive number of comparisons and, thus, it is not done. Local optimums are used instead, providing inconsistent solutions. Along with the scarcity of efficient NLP solvers and the numerical issues created by large problem instances, the depth of research using AC-OPF is reduced. Besides, the requirement of convexity frequently needed for applications limit its usage. Hence, many approximations and relaxations have been developed for the AC-OPF [15]. In general these are simplifications and are easier to solve, but they ignore some parts of the problem. As a result, it is imperative to understand and weigh the advantages and compromises of each formulation when choosing one to use in a model. The Linear DC approximation is a linear formulation which partially represents power voltage. This formulation makes some assumptions for linearization purposes: voltage magnitudes are sufficiently near nominal value (one), angle differences are close to zero, and there are practically no power losses. In all effects, a purely active power system linear model. Still it is important to notice that, while this model seeks to approximate the feasible region of the AC power flow, it may not include the entire feasible region, including the global optimum. Instead of approximating parts of the power system, it is possible to relax some of the nonlinear power flow constraints. The resulting relaxations, when solved to optimality, provide lower bounds to the original problem because their feasible sets include all the solutions of the original problem. Different convex relaxations were proposed for the optimal power flow problem, here we highlight some of the most important ones. Convex relaxations are specially useful because their solution can be proved to be optimal. The Copper Plate is the simplest linear relaxation and models the system by a centralized energy pool, relaxing transmission lines limits and the Kirchhoff's laws. Easily implementable and solvable, this is the most simplified linear formulation, neglecting the entire grid of the problem. However, because of the simplifications of this model it produces only one shadow price, i.e. the marginal price of energy in each node of the network are equal. This may send a distorted signal to the agents responsible for the system, and may return the most unfeasible solution across all formulations. Another linear relaxation, the Transportation model, or tube model, extends Copper Plate by adding line limits. Locational restrictions are better represented and the value of transmission lines are clearer. Although, by completely ignoring the power voltage, network design becomes less relevant and incorrect system analysis is possible.
The SOC relaxation is a non-linear convex relaxation that is tighter than the linear versions, i.e., its feasible region is strictly contained inside them. This formulation relaxes the non-convex constraints of the problem, composed by the bilinear product of the voltage variables, by neglecting the phases of the voltages and saving only their branch wise difference and magnitudes. The resulting problem may be specified as a second order cone formulation, hence its name.
The SDP relaxation deals with the non-convex constraints of the problem by using the fact that they define a positive semi-definite matrix with rank 1. The relaxation comes from removing the rank 1 restriction. The feasible region of this formulation is contained within the feasible region of the SOC relaxation, providing a better bound to the original problem. The QC relaxation exploits the polar form of the non-convex constraints and uses known convex envelopes to relax each non-convex term present. These envelopes retain stricter links between the voltage variables, producing a tighter relaxation than the SOC formulation. These, and other relaxations, provide a good alternative to solve the original problem. Besides bounding the optimal value of the original problem, they have sufficiently good solutions for real world applications. In fact, by being convex, these relaxations are suitable to be used in many solution methods for multistage problems, for instance SDDP for which convergence is guaranteed if subproblem are convex. The PowerModels.jl package, a framework for steady-state power network optimization, is able to construct these and other mathematical programming formulations of the OPF problem. Since PowerModels.jl uses JuMP.jl to construct these formulations, they can be passed to a variety of solvers. This allows the user to easily choose an approximation or relaxation, solve it, and then discuss and compare the impacts of using different relaxations and approximations in the planning and operation of the economic dispatch problem.

Solution method
As we have seen, the hydrothermal dispatch is a complicated problem with different network formulations. The OPF problem, is only a part of a sub-problem composing a multistage stochastic program. As discussed previously, solving a multistage stochastic program has its own difficulties, and requires specific and efficient algorithms.
In a multistage stochastic program, we are faced with a cost-to-go function: ω t+1 )] also depends on a cost-to-go ρ t+2 [Q t+2 (x t+1 , ω t+2 )], and the evaluation of those functions can be expensive. The crucial step that facilitates the solution of these problems is to construct approximations of the cost-to-go functions, recursively, going backward in time. Thus, the optimal value of the first stage problem approximates the optimal value of the corresponding multistage problem.
For the construction of this approximation a widely used method is dynamic programming, which evaluates the function in a range of discrete values of the state variable for further interpolation. However, this method becomes intractable with the growth of the state dimension (commonly referred as the curse of dimensionality of dynamic programming). A solution for this was proposed by [16] with a method called stochastic dual dynamic programming (SDDP). The methodology, simply posed, approximates the cost-to-go function by the maximum of a set of linear hyper-planes called cuts. SDDP is based on an interactive construction of the cost-to-go function approximations. The procedure may be divided in subsequential forward and backward passes, where the first chooses points in which the second will update the current approximation of the functions. For a detailed discussion of the SDDP algorithm, see [16,8].
Other solution algorithms have been proposed to solve multistage problems, such as progressive hedging and nested Benders. However, for most large instances, none has proven to be viable and efficient alternatives. SDDP has been extensively used to solve and plan hydrothermal dispatch operations since its original publication.
On the other hand, this method is a complex and difficult to implement algorithm. Moreover, to make it flexible enough for different applications while not compromising performance is a main issue. Hence, the historical scarcity of reliable and open software versions has limited development and discussions of its contributions. Yet, new and efficient implementations have risen with the advance of open-source languages as the Julia language. One such implementation is SDDP.jl, a Julia/JuMP package for solving large multistage convex stochastic optimization problems using SDDP. It provides a practical and efficient way to find the solution to the hydrothermal dispatch problem. With this package it was possible to define the dynamics and random variables of our problem and solve the instantiated model with the SDDP algorithm.

HydroPowerModels.jl
HydroPowerModels.jl uses the PowerModels.jl and SDDP.jl packages to implement and solve different hydrothermal dispatch formulations. It provides an interface to easily solve and simulate hydrothermal dispatch models and allows the creation of a collection of hydrothermal problems described in input files for the package (following the PowerModels.jl standard), thereby helping the discussion of methodology and the resulting policies for specific case instances.
In contrast to the previously available official software for hydrothermal dispatch models, the proposed package is part of an academic open-source effort. This helps to promote the continuous improvement of models and solution algorithms for the research community.
Other academic implementations of SDDP have being developed and may be applied to the hydrothermal dispatch problem. However, the advances of the Julia Language and JuMP.jl are a much more adequate framework than those of MATLAB [5] or Python [19,12]. Moreover, a free and open-source tool can be of great help for the research community alternatively to commercial solvers as [13] and the renowned version from PSR Inc. Additional implementations of SDDP are also available in Julia [10] [11], but SDDP.jl [8] has proven an easy to learn, efficient version of SDDP that is flexible enough for the purposes of the HydroPowerModels.jl package. HydroPowerModels.jl is composed of different and useful functionalities, from compact case sharing to dispatch solution results visualizations. A work-flow of a simple usage of the package helps to give a basic overview: -Load case data from input files describing: Power network data; Reservoir facilities details and water network data; Inflow scenarios. A code example is presented in the next subsections to help clarify the usage of the package. Although, for a more extensive tutorial of the package, a detailed documentation is made available in https://andrewrosemberg.github.io/ HydroPowerModels.jl/latest/.

Usage
The usage of HydroPowerModels.jl follow the paradigms of the Julia language and the structure of the dependent packages. In order to access the available functionalities, first import Hy-droPowerModels.jl and an adequate solver: § ¤  Figure 1 shows the installed power available in the network (grouped by bus) using a logarithmic scale. The red nodes represent the thermal generators, the blue represent the hydro genera-tors. For comparison purposes, orange nodes have been added that are equivalent to average real power demands. As we can see from the plot of the grid 1, this appears to be a well balanced case, with similar installed hydro and thermal power capacity and with a reasonable average demand. In addition, it is a well distributed network, without any evident critical sections susceptible to impacting power flow problems. Those facts are indications of a significant hydro-generation optimal dispatch without large complications. For this study, a 52 stage horizon planning and simulation have been executed using the following case parameters:

Results
The simulate command returns a detailed dictionary of the execution. In order to plot those results returned by the simulate function, you may choose from a variety of methods, including the function plot_aggregate_results(). This function receives the dictionary results and generates the most common aggregated variable plots, which best summarize simulations of a hydrothermal dispatch: § ¤ p l o t _ a g g r e g a t e d _ r e s u l t s ( r e s u l t s ) ¦ ¥ Figures 2 to 9 show the output from the above command. As mentioned, the plots are of aggregated quantities, but the methods used to aggregate were chosen in order to help analysis. For example: the final nodal price is an average of nodal prices weighted by the contribution of local loads to the total demand; reservoir volume was grouped weighted by the amount of energy that could be produced by the stored water (as was the inflow of water). As expected the optimal dispatch of the simulations uses more hydro-generators, however it needs thermal-generators to met all demands without deficit. On this hydro-dominated system, the uncertain inflow is a driving factor of optimal dispatch. As we can see in Figure 9, the inflow has a strong seasonality component, resulting the significant seasonality trait observable in the variables of the policy simulations. 2-3, 5-8. Similar studies are possible for any case and formulation chosen, helping to analyze existing realistic cases and assess impacts of future system changes. For example, we reran the study using a SOC relaxation to compare results. Most quantities analyzed, as those presented in the aggregated results, are similar, but operation costs differ as it is observable in figures 10,11. This is expected, since the SOC relaxation is a more realistic representation of reality. Moreover, it is important to point out this extra cost would be even bigger if we used the policy of the transportation model to operate the SOC model in the simulations. In other words, there is a hidden cost when planning a dispatch associated with the fact we are not evaluating costs on the actual system we want to operate. Therefore, measuring the impacts of possible simplifications is a needed step in discussing hydrothermal economic dispatch. Hy-droPowerModels.jl intends to provide a common ground for discussions and analysis and a easy to use tool for research and applications.