AbstractEnvironmental constraints in hydropower systems serve to ensure sustainable use of water resources. Through accurate treatment in hydropower scheduling, one seeks to respect such constraints in the planning phase while optimizing the utilization of hydropower. However, many environmental constraints introduce state-dependencies and even nonconvexities to the scheduling problem, making them challenging to represent in stochastic hydropower scheduling models. This paper describes how the state-dependent maximum discharge constraint, which is widely enforced in the Norwegian hydropower system, can be embedded within the stochastic dual dynamic programming (SDDP) algorithm for hydropower scheduling without compromising computational time. In this work, a combination of constraint relaxation and time-dependent auxiliary lower reservoir volume bounds is applied, and the modeling is verified through computational experiments on two different systems. The results demonstrate that the addition of an auxiliary lower bound on reservoir volume has significant potential for improved system operation, and that a bound based on the minimum accumulated inflow in the constraint period is the most efficient.IntroductionHydropower is a dominant generation technology in the Nordic power system, accounting for 143  TWh/year or 96% of the total power production in Norway in 2017 (Statistisk Sentralbyrå). In the future, the Nordic power system will have tighter connections with Europe and an increasing proportion of intermittent renewable generation from, for example, wind and solar power. Rapid and unpredictable fluctuations in intermittent generation will offer new possibilities for controllable generation, such as regulated hydropower, to be able to respond to these fluctuations.In symphony with the ongoing power market changes, the physical and environmental requirements associated with hydropower operation are changing, through proposed revisions of hydropower concessions and the implementation of EU Water Framework Directive, for example. The directive strives to ensure sustainable use of water resources, balancing the multiple uses such as water hydropower, irrigation, water supply, flood control, and recreation (European Commission). Consequently, hydropower producers need to both adjust their operational schedules according to the new price patterns seen in the market and at the same time relate to new operational constraints. In this context, the hydropower producers need scheduling models that represent physical and environmental constraints in a precise and consistent manner.Operational planning (or scheduling) models have been widely used by Nordic hydropower producers for several decades. Methods and models for hydropower scheduling have been developed along different methodological tracks, as documented in the reviews in (Labadie 2004; de Queiroz 2016; Macian-Sorribes and Pulido-Velazquez 2020; Pérez-Díaz et al. 2021; Giuliani et al. 2021) and benchmarked in Côté and Leconte (2015). The early works on scheduling of hydropower reservoirs used the principles of stochastic dynamic programming (SDP) (see Stage and Larsson 1961). SDP decomposes the multistage planning problem into a sequence of single-stage subproblems that can be solved by backward induction (Tejada-Guibert et al. 1995; Kim and Palmer 1997; Nandalal and Bogardi 2007). SDP algorithms are still widely applied and have mild requirements to system representation (Giuliani et al. 2021), allowing representation nonconvexities. Variants of the SDP model, such as sampling stochastic dynamic programming (SSDP) (Kelman et al. 1990; Côté and Arsenault 2019) and robust stochastic dynamic programming (RSDP) (Kim et al. 2021) has extended the applicability of this type of algorithm. However, because state variables need to be discretized in the SDP algorithm, it suffers from the “curse of dimensionality,” leading to computationally intractable problems when considering systems with a large number of state variables. The stochastic dual dynamic programming (SDDP) introduced in Pereira and Pinto (1991) allows solving the scheduling problem without discretizing the state variables, and is therefore computationally tractable for systems with multiple reservoirs. The SDDP algorithm can be classified as a sampling-based variant of multistage Benders decomposition method (Birge and Louveaux 2011), requiring a convex problem formulation, and is normally solved by using linear programming (LP). Extensions of the SDDP algorithm have been frequently reported in recent literature (see Rebennack 2016; Poorsepahy-Samian et al. 2016; Helseth et al. 2016; Street et al. 2017; Zou et al. 2018). The SDDP algorithm is currently the most widely used method for medium-term hydropower scheduling in the Nordic market (Helseth et al. 2018). In the Nordic context, medium-term scheduling refers to the process of computing strategies for local or regional hydropower resources providing a long-term price forecast. In particular, the extended algorithm to incorporate uncertainty in power price presented in Gjelsvik et al. (1999) has become popular and is implemented in the software ProdRisk (SINTEF 2022).Environmental constraints in water resources systems come in many flavors, depending physical and legislative conditions. The technical literature tend to emphasize on ecologically acceptable flows, in terms of magnitude and rate of change (Pérez-Díaz and Wilhelmi 2010; Chakrabarti et al. 2011; Olivares et al. 2015; Guisandez et al. 2016). The impact of minimum flows and maximum ramping rates certainly limits the flexibility of the hydro system, but from a modeling point of view these constraints fit well into algorithms relying on a convex model formulation, such as SDDP. Other environmental constraints involve state-dependencies, which are not easily treated in the SDDP algorithm, as reported in Pereira-Bonvallet et al. (2016) and Helseth et al. (2020), and recently reviewed in Schäffer et al. (2020). The main complicating factor is the nonconvexities associated with such constraints. Ideally, the assessment of nonconvex environmental constraints within hydropower scheduling models should be performed using methodologies that allow explicit modeling of the nonconvexities. A recent reformulation of the SDDP algorithm, known as stochastic dual dynamic integer programming (SDDiP), has proven convergence also in the nonconvex case. SDDiP was first presented in Zou et al. (2018), and has later been applied to hydropower scheduling in Hjelmeland et al. (2018), and was further tested to deal with state-dependent discharge constraints in Helseth et al. (2020). Although the numerical results were promising, the computation times of the SDDiP method are prohibitive, and it is therefore fair to claim that the algorithm is not generally suited for operational use in its current form. Recent research on the application of SDP to nonconvex environmental constraints is presented in Schäffer et al. (2022).This work concerns the widely applied SDDP algorithm, and investigates different approaches for representing a particular type of environmental constraint that comes in many flavors in hydropower systems, namely the state-dependent maximum discharge. In particular, this work emphasizes on cases where the maximum discharge capacity depends on the reservoir level. There are numerous examples of this type of constraint in the Norwegian and Swedish hydropower systems, and two such cases from Norway are presented in the “Computational Experiments” section. This type of constraint is commonly motivated by the reservoir’s abilities to serve recreational purposes (tourism, lodging, fishing, canoeing, etc.) in the summer season.This paper elaborates on how the formulation of this constraint can be approximated in the SDDP algorithm in a transparent manner and without compromising computation time. A problem relaxation supported by a time-dependent auxiliary lower reservoir volume bound is proposed, and the performance of this approach is demonstrated in two separate case studies resembling Norwegian hydropower systems. When properly defined, the auxiliary reservoir bounds serve to tighten the problem relaxation, providing improved SDDP strategies. Results in terms of water values, simulated reservoir trajectories, and expected profits are presented to validate our approach. This work contributes with new knowledge on how to model this type of constraint within SDDP and quantify its potential impact on scheduling strategies and simulation results.Model DescriptionThis study assumes a risk-neutral, price-taking hydropower producer aiming at maximizing the expected profit from operating a hydropower system while respecting all physical and legislative constraints. This is the typical market context for producers in the liberalized Nordic market. The decision problem can be formulated as a multistage stochastic optimization problem, where the expectation is typically taken over the stochastic variables’ inflow and market price. The combined SDP/SDDP algorithm presented in Gjelsvik et al. (1999, 2010) is currently the most widely applied algorithm for solving this type of scheduling problem in the Nordic market. It combines the classical SDDP algorithm described in Pereira and Pinto (1991) with an outer layer based on SDP to incorporate the market price uncertainty.The model presentation emphasizes on the SDDP part of the algorithm, assuming a stochastic inflow and a deterministic price sequence. By treating the price as deterministic, the authors attempt to clarify the tight relationship between the studied constraint and the price pattern, without loss of generality. That is, the SDP-part of the algorithm can be easily added to the model formulation later on. It is assumed that the probability distributions of the inflow can be discretized, and that problem can be decomposed into weekly decision stages. Consequently, the realizations of the inflow for an entire week are known at the beginning of that week. The use of weekly decision stages is the standard approach when computing water values for producers in the Nordic market (Gjelsvik et al. 2010; Helseth et al. 2018), and provides a compromise between computational complexity and accuracy. There are no further time steps within each decision stage t in the model applied in this work.For all stages in the planning horizon 1⋯T a vector Xt is defined, comprising all decision variables for that stage. Associated with Xt there is a price vector λt. The overall objective is then to find an operating strategy to obtain (1) maxE{∑t=1Tλt⊺Xt+Φ(ST)}where E denotes the expectation operator; and Φ(ST) the end of horizon valuation of state variables in ST.Inflow ModelInflows to the hydropower reservoirs are considered as uncertain in Eq. (1) and are represented by a stochastic model explained in the following. A set of inflow time series Is,s∈S are available as a starting point. Each time series can be associated with one or multiple reservoirs. These values are first normalized to remove seasonal variation and resemble a weakly stationary process (2) where Ist is the vector of observed inflows for scenario s and stage t; μt the vector of average weekly inflows; and σt the vector of standard deviations. Next, a vector autoregressive inflow model of first order (VAR-1) is fitted to the normalized inflow according to Eq. (3) (3) In Eq. (3), the inflow zt is explained by the product of the correlation matrix ϕ and the inflow from the previous stage zt−1 and residuals ϵt. To make sure that the inflow model always generates nonnegative inflows to the mathematical model subsequently described, a three-parameter lognormal distribution is fitted to the residuals as outlined in Maceira and Bezerra (1997) and de Matos and Finardi (2012). For the interested reader, further discussions on inflow modeling within the SDDP algorithms can be found in Mbeutcha et al. (2021), and in particular approaches guaranteeing nonnegative inflows in Poorsepahy-Samian et al. (2016) and Raso et al. (2017).Weekly Stage ProblemIn the further model formulation, a watercourse is referred to as a set of connected hydropower modules h∈H. Each module h can have a set of upstream connected modules Ωh and a downstream module receiving its discharged water. A power station can be associated with h where the conversion from discharge to power is described by a piecewise linear PQ-curve, where water discharge through the station is modeled using one variable qhmtD per discharge segment m∈Mh. These segments will be used in decreasing order according to their discharge efficiency ηhm, assuming that ηhm decreases with m, as illustrated with three segments in Fig. 1. A set of decision stages t with no further time steps within each decision stage is considered. The weekly stage problem for stage t can be formulated as in Eqs. (4a)–(4j). Symbols related to the inflow model formulation were previously described (4a) αt=max(λt∑h∈Hpht−∑h∈HCvht−+αt+1)(4b) vht+∑m∈MhqhmtD+qhtS−∑u∈Ωh∑m∈MuqumtD−Fsh[σt⊺ϕ]szs,t−1=vh,t−1+Fsh(σstϵst+μst): πhvs=s(h),∀  h∈H(4c) pht−∑m∈MhηhmqhmtD=0∀  h∈H(4d) qhmtD−γhtQ¯hmD≤0∀  h∈H,m∈Mh(4e) vht+vht−−γhtV˜ht≥0∀  h∈H(4f) V_ht≤vht+vht−≤V¯ht(4g) αt+1−∑h∈Hπhtcvvht−∑s∈S[ϕπtcz]szs,t−1≤βtc+(πtcz)⊺ϵt∀  c∈C(4h) zs,t−1=zs,t−1*: πsz∀  s∈S(4i) pht,vht,vht−,qhmtD,qhtS≥0(4j) where λt is the power price; pht is the power generation; vht− is a slack variable for the lower reservoir bound with a marginal cost C; αt+1 is the future expected profit; vht is the reservoir volume; qhmtD is the discharge at PQ-curve segment m; qhtS is the spillage; Ωh is the set of modules upstream h; Fsh is the scaling factor converting inflow from series s to module h; ηhm is the discharge efficiency at the PQ-curve segment m; Q¯hmD is the maximum discharge limit for segment m; γht is a binary variable indicating if the discharge waterway is open (1) or not (0); V_ht and V¯ht are the lower and upper reservoir volume boundaries; V˜ht is the reservoir threshold; and πhtcv, πtcz, and βc are coefficients for Benders cut c. The notation [σt⊺ϕ]s is used to indicate the s element of the vector resulting from the multiplication within the brackets [ ]. Dual values πhv and πsz are obtained from Eqs. (4b) and (4h) and further applied in the SDDP algorithm as subsequently described.The objective in Eq. (4a) is to maximize profit from selling power, subtracted the cost of violating the reservoir boundary in Eq. (4f). Temporal discounting for time preference is not considered in this work. The variables (left-hand side) are separated from the parameters (right-hand side) in each constraint and the matrix notation is adapted to highlight the presence of the inflow state variable. With this setup, it is easier to see how constraints are being updated before re-solving (“warm starting”) the stage problem. Reservoir balances are stated in Eq. (4b) where the normalized inflow state zt−1 is transformed back to physical inflow, according to Eqs. (2) and (3). Production functions are described in Eq. (4c). Eqs. (4d) and (4e) control the environmental constraint, saying that discharge is not allowed if reservoir volume is lower than a predefined threshold volume V˜ht. If vht reaches the threshold V˜ht, the binary variable γht can be “switched on” to allow discharge qhmtD according to the maximum discharge boundary Q¯hmD per discharge segment m in Eq. (4d), as illustrated in Fig. 1.Benders cuts are given in Eq. (4g), constraining the future expected profit. Although the copy constraints in Eq. (4h) are not an active part of the optimization problem, it is included here for ease of computing the dual value πsz=(∂αt)/(∂zs,t−1). All variables are nonnegative except αt. The slack variables vht− are penalized with a cost C ensuring that they are only used if unavoidable, as will be further discussed in this paper.The presence of the binary variable γht makes the problem in Eqs. (4a)–(4j) a mixed integer linear programming (MIP) problem. If γht is relaxed so that so that γht∈[0,1], the problem of Eqs. (4a)–(4j) becomes an LP problem. This is a convex relaxation of the original formulation, facilitating the possibility of the SDDP algorithm to converge at an optimal policy for the relaxed problem formulation.Solution ProcedureBy relaxing γht the overall optimization problem becomes convex and can be solved by SDDP, as discussed and illustrated in this section. The overall problem is decomposed into LP stage problems defined by Eqs. (4a)–(4j), for each week and with given values of inflows. The stage problems are solved repeatedly in the SDDP algorithm and the reservoir balances in Eq. (4b) and Benders cuts in Eq. (4g) are updated according to changes in state variables. An SDDP iteration comprises a forward simulation and a backward recursion as subsequently described.Forward simulation: First, a sequence of inflow scenarios n=1,…,N covering the period of analysis from t=1,…,T are computed by randomly sampling residuals from the fitted three-parameter lognormal distribution and by use of Eq. (3). Subsequently, the weekly stage problem of Eqs. (4a)–(4j) is solved for each stage t along the simulated scenario i, and results are collected and state variables are updated for the next stage. The simulated state at the end of the stage is used as the initial state for the next stage. The forward simulation provides an updated set of state trajectories, as illustrated by following the black (thicker) lines and black dots in Fig. 2 forward in time for each scenario sample i. The forward simulation is used to obtain an upper bound J+ in Eq. (5) representing the first stage profit plus the future expected profit seen from the first stage, and the lower bound J− in Eq. (6) representing the expected simulated profit (5) J+=∑h∈H(λ1ph1−Cvh1−)+α1(6) J−=1I∑t=1T∑i=1I∑h∈H(λtpiht−Cviht−)Backward recursion: For each state trajectory obtained in the previous forward iteration, one starts from the state at the end of week T−1, and for each of the K sampled inflow realizations, often referred to as backward openings and illustrated by the white dots in Fig. 2, one computes the optimal operation for week T. From the solution of each of the K stage problems of Eqs. (4a)–(4j), contributions to cut coefficients are gathered as the dual values πhv from Eq. (4b) and πsz from Eq. (4h), and the right-hand side βct=αt*−∑h∈Hπhvvh,t−1*−∑s∈Sπszzs,t−1*. These cut contributions are later averaged over all inflow samples to form a Benders cut of type Eq. (4g). Thus, inflow realizations in a week T for scenario sample i1 are used to construct one cut for stage T−1, illustrated by the red horizontal line in Fig. 2. The residual vector ϵt=[ϵt1,…,ϵtK] should be state independent, and thus, ϵt is uncorrelated from stage to stage. Therefore, the cut constructed for scenario i1 and for stage T−1 will be valid for all the other scenarios at stage T−1, as illustrated by the red arrows and dashed horizontal lines in Fig. 2. This is often referred to as cut sharing (Infanger and Morton 1996) and is of crucial importance for the convergence characteristics of the SDDP algorithm. One then repeats the procedure for week T−1, and so on, to obtain an updated operating strategy.Convergence can be declared when the upper bound is inside the confidence interval of the expected simulated profit or when the upper bound stabilizes (Shapiro 2011). If this condition fails, new sample points are added and another backward recursion is performed to refine the approximation at each time step. More details on the theory and applications behind the SDDP method and its convergence can be found in Tilmant and Kelman (2007) and Homem-de-Mello et al. (2011).RelaxationThis section presents a discussion on how the environmental constraint controlled by Eqs. (4d), (4e), and (4j) can best be approximated to meet the convexity requirement of the SDDP algorithm without compromising the computational performance. In the “Enhanced Relaxation” section, it is elaborated on that the use of an auxiliary lower reservoir volume bound can be included to improve the approximation.Standard RelaxationThe relationships between the upper boundary (γQ¯D) applied to qD in Eq. (4d) and variables γ and v are shown in Figs. 3(a and b), respectively. From this figure, the discontinuity introduced in Eq. (4e) becomes apparent.The dots in Figs. 3(a and b) indicate positions where γ takes the value 0 or 1. If one restricts the solution space to binary values of γ, the upper boundary on qD follow the solid drawn line in the two figures. That is, γQ¯D equals Q¯D when γ=1 and 0 when γ=0, and γ can only be 1 when v≥V˜. If one relaxes γ to a continuous variable, the solution space is relaxed as shown with the dashed black lines in Figs. 3(a and b). The upper boundary on discharge gradually increases with an increasing γ [Fig. 3(a)], and a nonzero γ is allowed whenever v>0.Enhanced RelaxationIn many practical applications of such constraints, the threshold V˜ corresponds to a relatively high reservoir level that forms a target level in the summer season, and consequently the relaxation discussed so far is not a tight approximation of the original formulation. It may seem intuitive and tempting to tighten the relaxed constraint Eq. (4e) by linearizing from the initial volume vh,t−1 as shown in Eq. (7) (7) vht+vht−−γht(V˜ht−vh,t−1)≥vh,t−1∀  h∈HBecause vh,t−1 is a state variable, this would leave us with an expression comprising the product of the two variables γht and vh,t−1, which would violate the convexity requirement of the SDDP algorithm. The bilinear term (γhtvh,t−1) in Eq. (7) can, for example, be approximated using McCormick envelopes (Cerisola et al. 2012). However, it can be shown that such approximation of Eq. (7) leaves us with the challenge of finding a tight and valid time-dependent lower bound estimate on vh,t−1. This lower bound estimate is referred to as an auxiliary and time-dependent lower reservoir volume boundary denoted as V_h,t−1* in the following. Reformulating Eq. (7) gives (8) vht+vht−−γht(V˜ht−V_h,t−1*)≥V_h,t−1*∀  h∈HIf V_h,t−1* for the reservoir volume is defined, the linearization could start from these limits [and not vh,t−1=0, as in Eq. (4e)] to provide a tighter relaxation. This is illustrated with the stapled red line in Fig. 3(b), and implies changing Eqs. (4e)–(8). Note that with γht=0 in Eq. (8), V_h,t−1* becomes a lower bound on reservoir volume. The challenge is then to define the parameters V_h,t−1* to be used with Eq. (8) that allows the tightest possible relaxation of Eqs. (4a)–(4j) without compromising the SDDP convexity requirement.Due to the cut sharing property of the SDDP algorithm, cuts should be valid for a wide range of values of the state variables vh,t−1. Typically, this range becomes narrower with an increasing number of SDDP iterations. However, to facilitate cut sharing, the parameters V_h,t−1* cannot be adjusted dynamically. Moreover, the auxiliary lower bound should be low enough to avoid constraining system operation more than necessary, because this would conflict with the purpose of the environmental constraint. Because inflow is uncertain, situations where the slack variable vht− needs to be activated at marginal cost C may occur, directly impacting the strategies (cut parameters). Thus, any V_h,t−1*>0 should be carefully evaluated because it potentially impacts the water values according to the penalty C, providing an overly strong incentive to store water. From this discussion the authors recommend that the auxiliary and time-dependent lower reservoir volume bounds V_h,t−1* are defined prior to running the SDDP model and that the procedure for defining these parameters is transparent.The authors suggest basing the parameter V_h,t−1* on the minimum accumulated inflow, starting from the first week of the environmental constraint. Fig. 4 shows the development of min(V_h,t−1*,V˜t) in the constraint periods (weeks 18–35 and 18–48) for the two systems described in “Computational Experiments,” where V_h,t−1* is represented by the minimum (red curve) and average (gray curve) accumulated inflow. The inflow time series from which the minimum and average values are found comprise large (10,000) sets of scenarios sampled from the inflow model previously described prior to running the actual SDDP model. If one lets V_h,t−1* follow the minimum accumulated inflow (colored red in Fig. 4), the volume trajectories considered both in the forward simulation and backward recursion should be able to meet this auxiliary lower bound in the vast majority of encountered states. Thus, the replacement of Eq. (4e) with Eq. (8) will improve the SDDP strategies. If V_h,t−1* follows a higher trajectory, such as the average accumulated inflow (colored gray in Fig. 4), the problem formulation will be further constrained. This may provide better approximation of the environmental constraint, but may also lead to overly conservative strategies, as discussed further in “Computational Experiments.”Finally, it should be noted that a tighter relaxation will be possible for other types of scheduling models that do not rely on the cut sharing property, such as scenario-based approaches (Helseth et al. 2018) or multistage Benders decomposition (Diniz et al. 2018).Computational ExperimentsComputational Setup and OrganizationA computer model was established implementing the aforementioned SDDP algorithm. The model was used to obtain and verify strategies (cuts) obtained using the different approximations previously outlined for two Norwegian cascaded hydropower systems. The environmental regulation differs in the two systems, exemplifying flavors of the regulation regime.A scheduling horizon of 3 years was applied with weekly decision stages. For each of the watercourses an initial reservoir volume equal to 50% of the maximum reservoir volume was assumed. Inflow models were fitted for both systems using the corresponding inflow time series comprising 62 historical years. A total of 100 inflow scenarios were re-sampled in each forward iteration, and 12 backward openings were used at each stage in the backward iterations of the algorithm.Before applying the aforementioned fitted statistical model in the SDDP model, it was compared against the historical observations. Results from this comparison for the inflow series used in System 1 are shown in Fig. 5. The black lines are the 0 and 100 percentiles (dashed) and the mean (solid drawn) obtained when sampling 10,000 time sequences of one year from the statistical model. The blue lines are the 0 and 100 percentiles (dashed) and the mean (solid drawn) obtained from the historical data. The statistical model avoids negative values which is an important point in the experiments. Moreover, it provides the correct mean value and a good fit for the 0 percentile. Extremal inflow values are not captured that well. In summary, we find that the stochastic inflow model serves its purpose but has potential for improvements with respect to capturing the extremal (wet) values.The model was run with two deterministic price scenarios according to the weekly average NordPool system prices for the years 2018 and 2020 (Nord Pool), shown in Fig. 6, each of them repeated over the scheduling horizon. Possible correlations between inflow and market prices are neglected in the computational experiments. The model parameter C was tuned according to the maximum price level to avoid the use of vht− while producing.For each price scenario, four cases were run: 1.Reference (REF) using the formulation Eqs. (4a)–(4j) setting γht=1 in Eq. (4d) and γht=0 in Eq. (4e) (i.e., ignoring the environmental constraint).2.Standard Relaxation (SLIN) using the formulation Eqs. (4a)–(4j) with γht∈[0,1].3.Enhanced Relaxation Minimum (ELIN-MIN) using the formulation Eqs. (4a)–(4j) with γht∈[0,1], and replacing Eq. (4e) with Eq. (8). V_ht* is set equal to the minimum accumulated inflow, obtained from a set of 10,000 sampled scenarios from the SDDP inflow model.4.Enhanced Relaxation Average (ELIN-MED) with similar mathematical modeling as ELIN-MIN, but with V_ht* set equal to the average accumulated inflow, obtained from the same sampled scenarios as in ELIN-MIN.The sequence of model runs to prepare results for each of the two considered systems is subsequently explained. For both price scenarios and for all four cases, the SDDP model was run to obtain strategies [cuts of type Eq. (4g)] for each of them. Recall that the SDDP model is an optimization model based on LP where the nonconvex environmental constraints are relaxed. Then, for both price scenarios and for each of the four cases, system operation was simulated by using the respective cuts from the SDDP runs. In the simulations a fixed set of 4,000 out of sample inflow scenarios were applied. For each inflow scenario and for each week in the scheduling horizon, the system simulation solved the MIP problem described by Eqs. (4a)–(4j). Thus, the four different SDDP strategies were tested using a common and presampled set of scenarios and a common simulation model.The SDDP and simulation models were implemented in Julia, using the JuMP package (Dunning et al. 2017) and CPLEX 12.10 solver for solution of the optimization problems. All tests were carried out on an Intel Core i7-9850H processor with maximum frequency of 4.60 GHz and 64 GB RAM. Parallel processing was not applied.The stabilization of the upper bound was used as the convergence criterion for the SDDP model. As discussed in Shapiro (2011), the stabilization of the lower (for minimization problem) bound indicates that further runs of the algorithm do not significantly improve the constructed policy, and may serve as a practical convergence criterion.System 1–Description and ResultsThe first system description resembles the Bergsdalen watercourse located in the western part of Norway. An illustration of the topology and specification of some technical characteristics is provided in Fig. 7. For each reservoir shown in the figure the average annual inflow and the storage capacity are stated, both in mm3. Each power plant is identified with its installed capacity in MW. An environmental constraint is associated with reservoir 2, stating that no water should be discharged from the reservoir in between weeks 18 to 35 if the reservoir volume is lower than 145  mm3.Strategy from SDDP RunsFig. 8 shows the convergence process for both price scenarios and for each of the four cases, comparing the upper bounds on expected profit as defined in Eq. (5) for the first 80 iterations. The total computation time was approximately 90 minutes for all three cases with the environmental constraint, while case REF was 10% faster. Consequently, the different treatments of Eq. (8) do not imply notable differences in computational complexity. As expected, the converged upper bounds are lowest for the most constrained cases. Recall that case REF does not respect the environmental constraint and thus provides a too optimistic expectation of the future profit. The reduced expected profit seen in Fig. 8 when gradually tightening the approximation of the environmental constraint can be seen as the cost of a stricter constraint.Next, the strategies obtained for reservoir 2 from the four cases were studied by inspecting their cuts from the SDDP runs. The expected marginal value of water (water value) at a given stage (week) and state (reservoir volume) can be found as the coefficients (πv) of the binding cut for that stage and state. Note that water values provide important information for the producers regarding the valuation of stored water applied in the shorter-term operational decision making. Figs. 9 and 10 show the water values for reservoir 2 for weeks 15 and 30, respectively. These values are plotted as a function of the volume in reservoir 2, while fixing all other reservoir volumes to their corresponding expected values. The considered reservoir intervals in the two figures reflect the bands of reservoir states visited in the final SDDP forward simulations for the week in question (these bands are shown in Fig. 11 for the case REF). Water values differ most for 2018 prices, being generally higher and more sensitive to reservoir volume in the most constrained cases for both week 15 and 30. As a general note, the differences in water values for this system indicate that the improved representation of the environmental constraint will be more important to consider if one expects the price level to be relatively high in the constraint period.SimulationIn this section, the results obtained when simulating operation following the different strategies are discussed. The simulated average reservoir volume trajectories for reservoir 2 are shown in Figs. 11(a and b). The 0 and 100 percentiles for the REF case are included in the figure to indicate the variation. As expected, one observes in Fig. 11(a) that the tightening of constraint Eq. (8) leads to a higher simulated average reservoir trajectory, reflecting the importance of being able to produce at high summer prices. With the low summer prices for 2020, similar but far less pronounced differences can be seen.The expected increases in annual profit compared to the REF case are shown in Table 1. Relative to the annual profit in the range of EUR 30–50 million, the increases reported in Table 1 are modest. However, one should keep in mind that the improvements can be harvested simply by improving the representation of Eq. (8) and that producers will often hunt for such marginal improvements in a competitive market.Table 1. Increase in expected annual profit for System 1 compared to case REF (in 103  €)Table 1. Increase in expected annual profit for System 1 compared to case REF (in 103  €)CasePrice 2018Price 2020SLIN2059ELIN-MIN25211ELIN-MED3018In general, increased precision in constraint representation leads to a higher expected profit. With 2018 prices, the enhanced relaxation based on the highest auxiliary lower reservoir volume bound (ELIN-MED) provides the highest expected profit. For 2020 prices, the ELIN-MIN strategy is slightly more profitable. Because prices are lower in the summer season, reservoir levels according to the 2020 price scenario are lower as well, due to the reduced incentive to generate power in this period. As a consequence, the lower reservoir bound applied in ELIN-MED constrains the SDDP strategy more than what is considered optimal, and thus leading to a lower simulated profit.System 2–Description and ResultsThe second system resembles the main string of the Rjukan watercourse located in a mountainous area in Southern Norway. Note that the full system comprises two strings and several additional reservoirs. However, to ease the presentation this study emphasizes on the string comprising the reservoir with an environmental constraint. The topology and technical data for this system are provided in Fig. 12. The environmental constraint is associated with the upper reservoir, but is more complicated than in System 1. In the period between weeks 18 and 48, no discharge should take place if the reservoir volume is lower than 734  mm3. However, according to an interpretation of the concession, the reservoir volume available at the beginning of the constraint period could be discharged in the constraint period. The water available in week 18 will therefore be used differently and have a different value than water arriving to the reservoir at a later stage. The modeling of this constraint was facilitated by introducing a “virtual reservoir” (1b) in parallel with the “physical reservoir” (1), as illustrated in the box in Fig. 12. Operation of the physical and virtual reservoirs are set to obey the following rules: •All inflow to the upper reservoir enters the physical reservoir.•A separate volume constraint is introduced to ensure that the sum volume v˜1=v1+v1b does not exceed the maximum boundary of the physical reservoir.•The maximum reservoir boundary of the physical reservoir is set to zero in week 18 to ensure that all water is transferred to the virtual reservoir when reaching that week.•The maximum reservoir boundary of the virtual reservoir is set to zero at the end of week 48 to ensure all water is transferred back to the physical reservoir by then.•Water transfer between the physical and the virtual reservoir is only made possible in week 18 (from 1 to 1b) and week 48 (from 1b to 1).•A copy of the power station is associated with the virtual reservoir, and additional constraints are imposed to ensure the sum production does not exceed the physical limitations of the original station.•Discharge from the physical reservoir 1 is constrained according to Eqs. (4e) or (8) in the constraint period, referring to the sum volume v˜1=v1+v1b.•Discharge from the virtual reservoir 1b is not constrained by the environmental constraint.Strategy from SDDP RunsFig. 13 shows the convergence process for all cases and for both price scenarios, comparing the upper bounds for all cases for the first 80 iterations. The total computation time was approximately 130 minutes for all three runs with the environmental constraint, while case REF was 12% faster. As for System 1, the different treatments of Eq. (8) do not imply notable differences in computational complexity, and the converged upper bounds are lowest for the most constrained cases. The curves for ELIN-MED are substantially lower than the others. These large differences indicate that the model is not always able to meet the auxiliary lower reservoir boundary in the ELIN-MED case, imposing high penalties that are reflected in the cuts.The use of the slack variable vht− recorded from the backward recursion per SDDP iteration is illustrated in Fig. 14 for cases ELIN-MIN and ELIN-MED considering 2018 prices. The numbers should be related to the total number of LP problems solved in the backward recursion in a single iteration (156×100×7=109.200). As expected, the ELIN-MED case imposes more frequent use of the slack variable in the backward recursion, explaining the differences in upper bounds shown in Fig. 13.Figs. 15 and 16 show how the water values for reservoir 1 differ between cases REF and ELIN-MIN for both the physical (1) and the virtual (1b) reservoir in the beginning of the constraint period (week 20). Each figure is accompanied with a probability density plot of the two reservoir volumes according to the simulated values in week 20. The water is valued significantly higher for both reservoirs in case ELIN-MIN, demonstrating the impact of the environmental constraint. Water in the auxiliary reservoir can be operated more freely in the constraint period for case ELIN-MIN and is therefore valued higher than water in the physical reservoir in both price scenarios. Note that the REF case values water in both reservoirs similarly because the environmental constraint is ignored, giving the model no incentive to prefer one of the reservoirs over the other.SimulationThe simulated average reservoir volume trajectories for reservoir 1 (the sum of 1 and 1b) are shown in Fig. 17. The four cases differ significantly considering 2018 prices, with ELIN-MED strategies leading to the highest reservoir volumes. With 2020 prices, ELIN-MED still suggests high reservoir volumes, while differences between the three other cases are muted.The expected increases in annual profit compared to the REF case are shown in Table 2. The enhanced relaxation based on the minimum auxiliary lower volume boundary (ELIN-MIN) provides the highest expected profit in both price scenarios, improving the increased profit by EUR 458,000 and EUR 165,000 for price scenarios 2018 and 2020, respectively. The ELIN-MED cases provided overly conservative strategies, due to the frequent use of penalties in the backward recursion illustrated in Fig. 14, leading to lost opportunities for power generation at favorable prices compared to ELIN-MIN. Moreover, by further inspecting the results, it was found that the ELIN-MED cases led to higher spillage than ELIN-MIN.Table 2. Increase in expected annual profit compared to case REF (in 103  €) for System 2Table 2. Increase in expected annual profit compared to case REF (in 103  €) for System 2CasePrice 2018Price 2020SLIN2,107463ELIN-MIN2,565628ELIN-MED2,329534ConclusionsThis work concerns the treatment of volume-dependent maximum discharge constraints in hydropower scheduling models based on SDDP. This state-dependent environmental constraint introduces a pronounced nonconvexity in the scheduling problem formulation, challenging the use of the SDDP algorithm. A combination of constraint relaxation and auxiliary and time-dependent lower reservoir volume boundaries is proposed to deal with this constraint without sacrificing the advantageous computational performance of the SDDP algorithm. If carefully defined, the auxiliary volume boundaries serve to tighten the relaxation of the environmental constraint, leading to a more accurate model representation.The proposed modeling approach is applied to two hydropower systems resembling Norwegian watercourses, demonstrating the impact on water values and simulation results as well as the expected profit improvements, considering uncertainty in inflow and deterministic prices. The water values, which inform the producer when and where to discharge water for power generation, are found to be sensitive to the environmental constraint.The computational experiments demonstrate the importance of defining the auxiliary lower reservoir boundaries high enough to provide a tight constraint relaxation, and low enough to prevent exaggerated use of slack variables for the environmental constraint. Thus, the definition of these boundaries should serve to naturally constrain the solution space while always being a relaxation of the original (MIP) formulation. For the general case, the authors recommend using the minimum accumulated inflow over the constraint period as a time-dependent auxiliary lower reservoir volume bound.The authors find that the use of the proposed approach increases the expected annual profits with EUR 47,000 for System 1 and EUR 458,000 for System 2 compared to the standard relaxation approach in a scenario with high prices in the constraint period. Although the type of constraint considered is typically active in the low-load summer period, the ongoing transition to a low-carbon power system may lead to more frequent occurrences of high prices throughout the year, underlining the importance of considering high constraint-period prices. Moreover, the results demonstrate that the accuracy provided by using auxiliary lower reservoir bounds based on the average accumulated inflow is system dependent, and that this approach may lead to overly conservative strategies. In particular, results from System 1 show that a volume bound based on the average accumulated inflow may in some cases further improve results, while this was not found for System 2.Further validation of the presented approach is needed to provide more robust estimates of result improvements. Further work may involve validation using an SDDP model with stochastic modeling of both price and inflow, and the application to different case studies. A detailed comparison between the proposed SDDP-based modeling against other algorithms (such as SDP or SDDiP) for the presented systems could provide important knowledge regarding the importance of accurate treatment of the nonconvex environmental constraint.Data Availability StatementAll data and code that support the findings of this study are available from the corresponding author upon reasonable request.AcknowledgmentsThis work was funded by The Research Council of Norway through Project No. 257588.References Birge, J. R., and F. Louveaux. 2011. Introduction to stochastic programming. 2nd ed. Berlin: Springer. Cerisola, S., J. M. Latorre, and A. Ramos. 2012. “Stochastic dual dynamic programming applied to nonconvex hydrothermal models.” Eur. J. Oper. Res. 218 (3): 687–697. Chakrabarti, B. B., N. Newham, D. Goodwin, and C. Edwards. 2011. “Wind-hydro firming with environmental constraints in New Zealand.” In Proc., of IEEE Power and Energy Society General Meeting. Piscataway, NJ: IEEE Power & Energy Society. Côté, P., and R. Arsenault. 2019. “Efficient implementation of sampling stochastic dynamic programming algorithm for multireservoir management in the hydropower sector.” J. Water Resour. Plann. Manage. 145 (4): 05019005. Côté, P., and R. Leconte. 2015. “Comparison of stochastic optimization algorithms for hydropower reservoir operation with ensemble streamflow prediction.” J. Water Resour. Plann. Manage. 142 (2): 04015046. de Matos, V. L., and E. C. Finardi. 2012. “A computational study of a stochastic optimization model for long term hydrothermal scheduling.” Int. J. Electr. Power Energy Syst. 43 (1): 1443–1452. Diniz, A. L., F. D. S. Costa, M. E. Maceira, T. N. dos Santos, L. C. B. dos Santos, and R. N. Cabral. 2018. “Short/mid-term hydrothermal dispatch and spot pricing for large-scale systems—The case of Brazil.” In Proc., 20th Power System Computation Conf. Piscataway, NJ: IEEE. Dunning, I., J. Huchette, and M. Lubin. 2017. “JuMP: A modeling language for mathematical optimization.” SIAM Rev. 59 (2): 295–320. Giuliani, M., J. R. Lamontagne, P. M. Reed, and A. Castelletti. 2021. “A state-of-the-art review of optimal reservoir control for managing conflicting demands in a changing world.” Water Resour. Res. 57 (12): 1–26. Gjelsvik, A., M. M. Belsnes, and A. Haugstad. 1999. “An algorithm for stochastic medium-term hydrothermal scheduling under spot price uncertainty.” In Proc., 13th Power System Computation Conf. ETH Zürich, Switzerland: Fachgruppe Energieübertragungssysteme. Gjelsvik, A., B. Mo, and A. Haugstad. 2010. “Long- and medium-term operations planning and stochastic modelling in hydro-dominated power systems based on stochastic dual dynamic programming.” In Handbook of power systems, 33–55. Berlin: Springer. Guisandez, I., J. I. Perez-Diaz, and J. R. Wilhelmi. 2016. “The influence of environmental constraints on the water value.” Energies 9 (6): 1–21. Helseth, A., M. Fodstad, and B. Mo. 2016. “Optimal medium-term hydropower scheduling considering energy and reserve capacity markets.” IEEE Trans. Sustainable Energy 7 (3): 934–942. Helseth, A., B. Mo, and H. O. Hågenvik. 2020. “Nonconvex environmental constraints in hydropower scheduling.” In Proc., Int. Conf. on Probabilistic Methods Applied to Power Systems. Liege, Belgium. Helseth, A., B. Mo, A. L. Henden, and G. Warland. 2018. “Detailed long-term hydro-thermal scheduling for expansion planning in the Nordic power system.” IET Gener. Transm. Distrib. 12 (2): 441–447. Hjelmeland, M. N., J. Zou, A. Helseth, and S. Ahmed. 2018. “Nonconvex medium-term hydropower scheduling by stochastic dual dynamic integer programming.” IEEE Trans. Sustainable Energy 10 (1): 481–490. Homem-de-Mello, T., V. L. de Matos, and E. C. Finardi. 2011. “Sampling strategies and stopping criteria for stochastic dual dynamic programming: A case study in long-term hydrothermal scheduling.” Energy Syst. 2 (14): 1–31. Infanger, G., and D. P. Morton. 1996. “Cut sharing for multistage stochastic linear programs with interstage dependency.” Math. Program. 75 (2): 241–256. Kelman, J., J. R. Stedinger, L. A. Cooper, E. Hsu, and S. Q. Yuan. 1990. “Sampling stochastic dynamic programming applied to reservoir operation.” Water Resour. Res. 26 (3): 447–454. Maceira, M. E. P., and C. V. Bezerra. 1997. “Stochastic streamflow model for hydroelectric systems.” In Proc., 5th Conf. on Probabilistic Methods Applied to Power Systems. Vancouver, BC, Canada: Tom Gutwin and B. C. Hydro. Mbeutcha, Y., M. Gendreau, and G. Emiel. 2021. “Benefit of PARMA modeling for long-term hydroelectric scheduling using stochastic dual dynamic programming.” J. Water Resour. Plann. Manage. 147 (3): 05021002. Nandalal, K. D. W., and J. J. Bogardi. 2007. Dynamic programming based operation of reservoirs, applicability and limits. Cambridge, MA: Cambridge University Press. Olivares, M. A., J. Haas, R. Palma-Behnke, and C. Benavides. 2015. “A framework to identify Pareto-efficient subdaily environmental flow constraints on hydropower reservoirs using a grid-wide power dispatch model.” Water Resour. Res. 51 (1): 3664–3680. Pereira, M. V. F., and L. M. V. G. Pinto. 1991. “Multi-stage stochastic optimization applied to energy planning.” Math. Program. 52 (10): 359–375. Pereira-Bonvallet, E., S. Püschel-Løvengreen, M. Matus, and R. Moreno. 2016. “Optimizing hydrothermal scheduling with non-convex irrigation constraints.” Energy Procedia 87 (12): 132–140. Pérez-Díaz, J. I., M. Belsnes, and A. L. Diniz. 2021. “Optimization of hydropower operation.” Comprehensive Renewable Energy 6: 84–104. Pérez-Díaz, J. I., and J. R. Wilhelmi. 2010. “Assessment of the economic impact of environmental constraints on short-term hydropower plant operation.” Energy Policy 38 (12): 7960–7970. Raso, L., P. O. Malaterre, and J. C. Bader. 2017. “Effective streamflow process modeling for optimal reservoir operation using stochastic dual dynamic programming.” J. Water Resour. Plann. Manage. 143 (4): 11. Rebennack, S. 2016. “Combining sampling-based and scenario-based nested Benders decomposition methods: Application to stochastic dual dynamic programming.” Math. Program. 156 (1): 343–389. Schäffer, L. E., A. Adeva-Bustos, T. H. Bakken, A. Helseth, and M. Korpås. 2020. “Modelling of environmental constraints for hydropower optimization problems—A review.” In Proc., Int. Conf. on the European Energy Market (EEM). Piscataway, NJ: IEEE. Schäffer, L. E., A. Helseth, and M. Korpås. 2022. “A stochastic dynamic programming model for hydropower scheduling with state-dependent maximum discharge constraints.” Renewable Energy 194 (July): 571–581. Stage, S., and Y. Larsson. 1961. “Incremental cost of water power.” Trans. Am. Inst. Electr. Eng. 80 (3): 361–364. Street, A., A. Brigatto, and D. M. Valladão. 2017. “Co-optimization of energy and ancillary services for hydrothermal operation planning under a general security criterion.” IEEE Trans. Power Syst. 32 (6): 4914–4923. Tejada-Guibert, J. A., S. A. Johnson, and J. R. Stedinger. 1995. “The value of hydrologic information in stochastic dynamic programming models of a multireservoir system.” Water Resour. Res. 31 (10): 2571–2579. Tilmant, A., and R. Kelman. 2007. “A stochastic approach to analyze trade-offs and risks associated with large-scale water resources systems.” Water Resour. Res. 43 (6): 7.

Source link

Leave a Reply

Your email address will not be published. Required fields are marked *