Lagrangian relaxation based heuristics for a chance-constrained optimization model of a hybrid solar-battery storage system

We develop a stochastic optimization model for scheduling a hybrid solar-battery storage system. Solar power in excess of the promise can be used to charge the battery, while power short of the promise is met by discharging the battery. We ensure reliable operations by using a joint chance constraint. Models with a few hundred scenarios are relatively tractable; for larger models, we demonstrate how a Lagrangian relaxation scheme provides improved results. To further accelerate the Lagrangian scheme, we embed the progressive hedging algorithm within the subgradient iterations of the Lagrangian relaxation. We investigate several enhancements of the progressive hedging algorithm, and find bundling of scenarios results in the best bounds. Finally, we provide a generalization for how our analysis extends to a microgrid with multiple batteries and photovoltaic generators.


Motivation
A photovoltaic power station (PPS) is a system consisting of solar panels and inverters to convert light into electricity. Solar power is uncertain, and, of course, not available for all hours. Hence, coupling a sufficiently large storage device with a PPS can significantly increase its economic value by shifting energy to times when it is in higher demand, creating a hybrid satisfying the above inequality to be no less than 1 − ε, where ε is a predetermined threshold. Thus, if we "fail" in even a single i, we fail to satisfy the JCC. In this sense, the JCC represents a very stringent satisfaction criterion. Computationally, optimization models with JCCs are hard to solve [37], and several tailored decomposition approaches [34,36], as well as heuristics [6,60], are available for their solution. Theoretically, such optimization models are known to be N P-hard [37]; further, the feasible region is generally non-convex [47], with a few exceptions such as [23].
Nonetheless, an optimization model with JCCs is an attractive modeling choice for ensuring reliability of a set of operations. Chance-constrained models find frequent applications, e.g., power systems [28,62], traffic flow [13,59], portfolio selection [1,32], etc. There exists a modeling choice with regards to the timing of observing the uncertainty. In the first case, the entire set of future scenarios is known after making the first-stage decision, x. In the second case, the uncertainty is revealed periodically and subsequent decisions need to be made in multiple stages. Such models are known as two-stage and multi-stage stochastic programs, respectively; see, e.g., [47]. Multi-stage models are computationally harder than two-stage models. Since optimizing over JCCs already poses a significant challenge, in this tradeoff between fidelity and computational effort, we use a two-stage model instead of a multi-stage model. In our model, we make a first-stage decision of 24 hourly day-ahead power promises. Then, we observe the previously uncertain solar power output and make all the second-stage decisions of charging/discharging the battery. We use a highly reliable regime of operations (ε ≤ 0.05). Our motivation comes from the fact that failure to satisfy the promise could lead to curtailment of future contracts with the supplier, and thus ensuring reliability is valued more than increasing profit. We describe our model in greater detail in Sect. 2. Power discharged from battery during hour t in scenario ω [kW] x ω t Energy stored in battery at hour t in scenario ω [kWh] w ω t Binary variable; 1 if battery is charging at hour t, or 0 if discharging

Optimization model
Boundary condition: x ω 1 = x 0 , ∀ω ∈ Ω. The objective function in (1a) aims to maximize profit; i.e., revenue from the promise minus the expected cost of operating (charging or discharging) the battery. A more sophisticated model could include piecewise linear or quadratic costs; see, e.g., [12,63]. The JCC in Eq. (1b) requires that the joint probability of meeting the promised energy to sell, by the available solar power (via s ω ) and discharging the battery (via q ω ), for the entire time horizon meets a threshold; excess solar energy can be used to charge the battery (via p ω ). If the solar power exceeds the need, we might need to curtail it. However, in this article we do not penalize this curtailment, and hence Eq. (1b) is expressed as an inequality. Here, ε is a positive quantity just greater than zero, such as 0.05. Constraint (1c) relates the energy stored by the battery in a subsequent hour with the previous hour; here, the quantity ηΔp ω t − 1 η Δq ω t is the energy charged or discharged by the battery during hour t. Constraints (1d), (1e) and (1g) ensure that both p ω t and q ω t cannot be simultaneously positive. Further, since η < 1, a fraction of energy is lost while both charging and discharging; i.e., we consume more than 100% while charging and supply less than 100% while discharging. Constraint (1f) ensures the minimum and maximum amount of energy stored in the battery, and thus Eq. (1d)-(1h) restrict the amount the battery can be charged or discharged. Lower bounds on x are typically included for electric batteries to maintain long and healthy lifetimes; see, e.g., [40]. Constraints (1g)-(1h) ensure the binary and non-negativity restrictions on the relevant decision variables.
One potential criticism of optimization models with JCCs, such as model (1), is that violations can be large for the 100·ε percent of failed scenarios in Eq. (1b). Although, Eq. (1b) restricts the number of violations to be a few, there is nothing in model (1) that restricts the magnitude of these violations. This well-known drawback is an artifact of chance-constrained programming in general, and is also reflected in our model. Several methods exist to bound the violations. We direct interested readers to [50] and [64] for two remedies using risk measures and envelope constraints, respectively. However, models that only penalize energy violations (or, do not employ chance constraints) can yield solutions that violate the promise with large frequency when it is economically advantageous. Our aim is to decide whether it is feasible to couple a solar plant with an appropriately sized battery unit, and whether this coupling would enable the hybrid unit to participate in the day-ahead market. Now, consider the following optimization model in which we drop the binary variable

Proposition 1
In an optimal solution to model (2), both p ω t and q ω t cannot be positive together for any t ∈ T , ω ∈ Ω.
Proof Assume an optimal solution of model (2) is given by Assume that p ω t and q ω t are positive, for some t ∈ T , ω ∈ Ω. We show that this solution cannot be optimal by proving a better solution (one that gives a larger objective function value) always exists. First, assume that η 2 p ω t > q ω t , for this t and ω pair. Then, consider a new solution We observe that this solution is feasible to model (2); further, it results in a strictly larger objective function than that by the optimal solution. This is a contradiction.
Next, assume that η 2 p ω t ≤ q ω t , for this t and ω pair. Again, consider a new solution This solution is feasible to model (2), and results in a strictly larger objective function than that by the optimal solution. This is a contradiction.
Thus, in any optimal solution to model (2), for any t and ω at most one of p ω t or q ω t is positive.
Proposition 1 offers the advantage that there are no binary w = [w t ] ω t=1,2,...,|T | variables in the second-stage. Although, a reformulation of constraint (1b) would still entail secondstage binaries (see, Eq. (3)), the reformulation in model (2) could save computational effort by reducing |T ||Ω| binary variables. There is some evidence, however, that having additional binary variables can actually speed up the computations by offering mixed-integer program (MIP) solvers additional variables to branch and cut on; see, e.g. [41,56]. However, we choose the reformulation without the binary w variables as it lends itself better to the heuristics we describe in the later sections of this article. We reformulate the probabilistic constraint in Eq. (2b) using a big-M approach as follows: Here, · rounds its argument down to the nearest integer. Henceforth, we use model (2) with Eq. (2b) reformulated as Eq. (3) as our proposed model. Below, we provide a sufficiently large value for the big-M of Eq. (3).

Proposition 2
The following inequality is a valid inequality for model (2): denotes the l th largest solar power value at time t.
Proof We first seek an upper bound on the quantity q ω t − p ω t . If q ω t = 0, a valid upper bound on this quantity is 0. If p ω t = 0, a valid upper bound on this quantity is min{Q, η(X − X )}; this is obtained from Eqs. (2c), (2e) and (2f). Thus, from Proposition 1, q ω t − p ω t ≤ min{Q, η(X − X )}, ∀t ∈ T , ω ∈ Ω. Now, the proof replicates that of Proposition 5.2.1 of [55]. By constraint (3b) effectively at most N ε of the z ω variables can take a value of 1; i.e., for every t at most N ε of the scenarios can be removed. Thus, the proposed inequality is valid.
We denote the upper bound, obtained from Proposition 2, on the y variables as y max , ∀t ∈ T . And, in our computational experiments of Sect. 6 we use the tightest big-M given from Corollary 1, namely M ω t = y max t − s ω t , ∀t ∈ T , ω ∈ Ω. Despite this, as we show later in this article, model (2) is computationally challenging for relatively larger values of ε (such as 0.05) or, for a large number of scenarios. Model (2) is a classic two-stage joint chance-constrained stochastic program with recourse. The second stage has integer restrictions (the variable z is binary) as well, which makes it even more challenging [2]. To this end, we investigate a Lagrangian relaxation scheme as well as a heuristic, in addition to the traditional approach with the big-M constraints.

A Lagrangian dual procedure
First, we note that constraint (3b) links the scenarios together. This motivates a straightforward Lagrangian relaxation of model (2). Lagrangian relaxation schemes have been studied before in various settings; see, e.g., individual chance constraints for unit commitment [42], relaxing non-anticipativity constraints [3], and coupling scenario decomposition with Lagrangian decomposition [60].
iter ← iter+1, update time to the cumulative wall-clock time. 10: end while for model (2) for the step-size update rule in Step 4. In principle, this lower bound can be computed using any feasible solution; setting N ε of the z variables to 1, and the rest to 0, provides a feasible solution. However, since the termination criteria of the Lagrangian relaxation procedure depends on this lower bound, good feasible solution values (those giving a larger value of the objective function) assist in convergence.
Lower Bounding Heuristic: To this end, we obtain a feasible solution using the following simple procedure described in [3]. We first solve model (2) separately for each ω with the corresponding z ω ← 0; i.e., a total of |Ω| problems. We rank the corresponding objective function values from the lowest to the highest and choose the first N ε (i.e., the worst performing scenarios) of these to set z ω ← 1. Then, with N ε of the z values fixed to 1, we solve model (2), to obtain a valid lower bound (LB) to the problem.
For an initial upper bound to the Lagrangian relaxation procedure, we use the linear programming (LP) relaxation of model (2). In the following proposition, we show that model (4) is never infeasible.

A progressive hedging heuristic for the Lagrangian dual
The progressive hedging (PH) algorithm can be used to solve stochastic programs that have easy-to-solve individual scenario problems [48]. While PH cannot be directly applied to model (2) due to the interlinking of the scenarios via constraint (3b), the relaxed model (4) lends itself perfectly to the decomposition structure required for PH. In this section, we leverage this idea to develop computationally cheap heuristics for the Lagrangian dual. First, we convert our maximization problems to a minimization by reversing the signs on the objective coefficients (R, C c , C d ), since most of the existing literature on PH considers minimization [16,60,61]. Then, model (5) is the PH sub-problem; we summarize the PH heuristic in Algorithm 2.
Solve model (5), ∀ω ∈ Ω, obtain optimal y ω t . 5: Here, u ω t is the scenario-dependent weight, y t is the scenario-averaged value for y t , and P H ω is the optimal value of each sub-problem. Solving the LR dual approximately via the PH algorithm offers several advantages and disadvantages compared to the naive solution method. We discuss some of these in Sect. 6.5.
In the presence of discrete decision variables (as is the case with our model), the PH algorithm can only be used as a heuristic as its convergence is not guaranteed [35,61]. Since we are interested in solving the Lagrangian dual, to obtain an upper bound to the true problem given by model (2), any upper bound to the Lagrangian dual is also an upper bound to the true problem. Consider Fig. 1; here, LB is the lower bound obtained via the heuristic procedure described in Sect. 3, z * is the true optimal value of model (2), and LD is the optimal value of model (4). Since convergence of the PH algorithm is not guaranteed, if we solve the Lagrangian dual (LD) via the PH algorithm, the PH optimal objective function value could be larger or smaller than LD.
If the PH algorithm results in a bound larger than LD at any subgradient iteration of Algorithm 1, then the bound is still a valid upper bound to z * . However, if the PH algorithm results in a bound smaller than LD at any subgradient iteration of Algorithm 1, then the bound could also be smaller than z * (or, even smaller than LB). To prevent this case, we seek valid (i.e., guaranteed to be larger than LD) upper bounds to LD via the PH algorithm (which may or may not be optimal to PH). Such bounds can be obtained by solving model (5) without the quadratic term in the objective using the procedure described in [16]. This bound is represented in Fig. 1 by P H U B , and is computed in Algorithm 2 in Steps 1 and 6.

Data sources
For the computational experiments in this article, we use η = 0.9 [29], and values of R t from the Electric Reliability Council of Texas (ERCOT) for the year 2012 [15,56]. Table 1 presents the hourly marginal selling prices. Further, we use X = 960kWh, with a 2V/1000Ah  (2). An arrow from a source node to a sink node indicates the optimal value of the source node is smaller than the sink node. Here, z * is the true optimal value, LB is the lower bound from a heuristic, LD is the optimal value of the Lagrangian dual, P H U B is an upper bound for the Lagrangian dual obtained via PH, and z * and z * are the upper and lower bounds for z * if solved sub-optimally. An optimality gap of β between an upper bound U and a lower bound L indicates β = U −L U . The optimal PH value is not shown and could be larger or smaller than z * rating, from the system with 480 lead-acid batteries described in [68]. To compute the cost of charging or discharging the battery, we again use the formula given by [68]: 390α Q I , where α is an effective weighting factor, I is the initial investment cost to purchase the battery, and 390Q is an approximation for the total Ah throughput of the battery [27]. For the system in [68], a X kWh sized battery is equivalent to 2000 960 X Ah rating. We assume α = 0.5, and use a lead-acid battery cost of $200/kWh from [52]. An approximately similiar battery cost is available from [46]. Then, C c = C d = $0.0256/kWh. Further, we assume X = 0.2X , P = Q = 0.5X , and the battery is 50% charged initially (i.e., x 0 = 0.5X ).
To generate solar power scenarios, we take hourly year-long historical solar power output from a site described in [19]. We use ARMA( p, q) models for the fourteen hours with sunlight, [06:00-20:00). For each hour, we verify the stationarity of the time series and test a number of ARMA( p, q) models to find the best. For details on the forecasting method we used, see [57]. Then, using Monte Carlo sampling from the best ARMA models for each hour, we create a set of hourly scenarios. Figure 2 plots a sampling of 300 scenarios, as well as the hourly medians and the ten percent quantile. All tests in this article are carried out with GAMS 25.1.2 using CPLEX 12.8 [49] on an Intel Core i7 2.8 GHz processor with 16 GB of memory. We attempt to solve the original model (2) to a MIP gap of 0.1%, with a maximum time limit of 40 min. For the rest of the experiments in this article, we also consider an optimality gap of 0.1%; i.e., if the relative gap between the best upper and lower bounds is less than or equal to 0.1% we consider the problem solved. We use multiple sets of scenarios-the smallest instance with 150 scenarios and the largest instance with 1200 scenarios-with two reliability regimes having ε = 0.01, 0.05.

Analysis of naive solution method
In Table 2 we present our first results obtained by solving model (2) naively. The "z * " column denotes the 24-h profit obtained in dollars; if the problem is not solved to the optimal tolerance, we present the best feasible solution. The ε = 0.01 instance is relatively easier (smaller computation time) than the ε = 0.05 instance due to the fewer number of combinatorics involved, and possibly due to a smaller value of the big-M via Corollary 1; the problem can be solved to the optimal tolerance for all instances except |Ω| = 1200. The gaps increase consistently when |Ω| increases; beyond 600 scenarios the ε = 0.05 case cannot be solved to the optimal tolerance.

Analysis: solution profile
Next, we provide a profile of an optimal solution. We choose the ε = 0.05, |Ω| = 300 instance; i.e., the third row of Table 2. To this end, we first solve the instance to a MIP gap of 0%, with an indefinite time limit. The optimal profit for the 24 hour horizon is $283.79; i.e., the value reported in Table 2 is indeed optimal to a MIP gap of 0%. The selling price, R t , exceeds the operational cost of the battery, C c = C d , for the nine hours [13:00-22:00), see Table 1. Solar power is non-zero, for at least one scenario, only for the fourteen hours [6:00-19:00), see Fig. 2. We provide the characteristics of the optimal solution in Table 3.
No power is promised during the seven night hours from midnight to the hour 6:00AM (i.e., until 6:59AM) and the three evening hours from 20:00 to the hour 22:00 (i.e., until 22:59PM). This is because the cost of operating the battery is higher than the reward for the seven night hours, and further there is no solar power available in these hours. In the three evening hours, although the cost of operation is (slightly) less than the reward, the model still chooses to not make a promise. This is because there is no solar power, and the model prefers to completely drain the battery to its lower limit in the last hour. Most of the hours meet the energy promise largely through the solar power, which is dictated by the economics of model (2). However, there is a good synergy between the charging and discharging times of the battery; i.e., (i) the promise made does not rely entirely on the solar power, and (ii) for the same hour the model chooses to charge or discharge based on the observed solar power scenario. An extreme example is the hour 17:00-17:59. In this hour, the battery is discharged (on average) to over 80% of its discharge limit to increase profit-solar power is decreasing beyond this hour, and this hour has the largest marginal profit (R t − C c ).
Also, we note that excess solar power in the afternoon hours is used to charge the battery. Hence, we see that the model is successfully using the battery to shift solar power to times when energy commands a higher price. Finally, as we mentioned before, at the last hour of 23:00-23:59, the battery is completely discharged to its permissible limit. To counter this end-effect, the model can be run in a rolling-time horizon or a boundary condition could be imposed; see, e.g., [30].

Analysis: storage size
Next, we analyze the impact of the size of the battery storage on the optimal profit. Similar to Sect. 6.2.1, we choose the ε = 0.05, |Ω| = 300 instance. We recompute the optimal profit from model (2) by varying the maximum storage capacity of the battery, X , in increments of 10%. Since we set x 0 = 0.5X and X = 0.2X in our computational experiments, the lowest feasible limit of X is 0.5 X corresponding to a 50% decrease. Figure 3 presents the paretooptimal curve for the percentage changes in the optimal profit by varying the maximum storage capacity of the battery. The red dot presents the base case; i.e., a profit of $ 283.79 for a battery size of X = 960 kWh. There is no change in the optimal profit beyond a size of 1344 kWh; this suggests that for the given parameters of the model a battery of size X = 1344 kWh is sufficiently large.

Analysis: quality of the solution
In this section, we perform an out-of-sample validation to analyze the performance of the optimal power promise. Since model (2) is solved using a given set of scenarios by a socalled sample average approximation (SAA) of the "true" model, this solution is naturally sub-optimal or infeasible to a new set of scenarios [4]. Thus, we seek to answer the following question: to what minimum reliability level can we guarantee a given power promise profile to safeguard us against new scenarios? Intuitively, we expect this reliability level to be lower than the chosen 1 − ε level. To this end, we first independently sample ten batches of 300 scenarios each. We fix the promised power output to the optimal power output of the base case instance with |Ω| = 300 instance. And, then we find the minimum number of violations required to ensure feasible operations (i.e., constraints (2c)-(2g)) in each of the ten batches. We construct a 95% confidence interval (CI) around the number of violations using a t-distribution. For the SAA solved with a 95% reliability regime, the 95% CI of the true reliability is (85.1%, 89.5%). Analogously for the 99% reliability regime, the 95% CI is (92.9%, 95.4%).

Analysis of lower bound heuristic
In this section, we briefly analyze the performance of the lower bounding heuristic mentioned in Sect. 3. Table 4 presents our results. The "LB" values are the lower bounds to model (2), see also Fig. 1. We solve all the scenario sub-problems using the Gather-Update-Solve-Scatter extension of GAMS [10]. This scheme is very efficient; 1200 scenario problems plus the fixed scenario problem can be solved in about 30 seconds. To compute the "Gap", we use the best known upper bound-this is available from either the naive solution (Table 2),  Table 4 indicates the gap is smaller than the optimality tolerance. Thus, six of the twelve problems have a gap under 1%, while two are solved to optimality.

Analysis of the Lagrangian relaxation scheme
In this section, we attempt to solve the Lagrangian relaxation of model (2) given by model (4) using Algorithm 1. We use a maximum of 10 iterations (i.e., M = 10). We use the same maximum time limit (time) and optimality tolerance (ψ) of 2400 seconds and 0.1%, respectively, as for the results in Table 2. Theoretically, this gap might not be attained even with a larger time limit as (i) different schemes for obtaining a lower bound give different corresponding gaps, and (ii) even the tightest upper bound from the Lagrangian procedure might not be attained as our problem includes binary variables. However, as we show later in this section our empirical results demonstrate that in some instances an optimality gap smaller than the naive solution is achieved. We use the LB bound from Sect. 6.3 in Step 4 of Algorithm 1. Further, we "warm-start" the z variables for the first LR iteration, using the optimal z values from the LB heuristic. Also, at each subgradient iteration, we warm-start the y variables with their optimal value from the previous subgradient iteration. In CPLEX, this can be done by setting the mipstart parameter to 1. Table 5 presents our results. All 10 subgradient iterations complete in the required time limit, except for |Ω| = 900, 1200. In our experiments, if ten LR iterations perform within the time limit, the optimal Lagrangian dual (denoted by "LR") is typically achieved around the sixth iteration. The four instances (ε = 0.05, |Ω| = 600, 900, 1200 and ε = 0.01, |Ω| = 1200) that are not solved to optimality in Table 3 all have smaller gaps in Table 5 than  Table 2, suggesting the Lagrangian relaxation is better suited for larger problems. Also note, that for these four instances, only a few subgradient iterations are completed. This again lends credence to the observation that good Lagrangian bounds are obtained fairly early in the iterations.
The optimal solution provided by the Lagrangian relaxation is infeasible to model (2). Thus, the algorithm is useful only for providing good bounds to model (2), and not for obtain- Table 5 Computational results for model (2) using Algorithm 1. "LB" and "LR" refer to lower bound and Lagrangian relaxation, respectively. The "Gap" column is the relative gap between the LB and LR. We use a maximum time limit of 2400 seconds (cumulative sum of LP, LB and LR) and a gap of 0.1%. For details, see Sect. 6 ing high quality feasible solutions. At best, the optimal y variables from the last iteration of Algorithm 1 could be used to warm-start model (2). This aspect of Lagrangian decomposition techniques is a well-known drawback; researchers have addressed this issue in a number of ways, see, e.g., [66,71]. A similar issue is encountered in Step 9 of Algorithm 2 as well; model (4) could be infeasible. However, if y t ≤ min{Q, η(X − X )} + s , ∀t ∈ T , then it follows from Proposition 3 that model (4) is assuredly feasible.

Analysis of the progressive hedging heuristic
The performance of the Lagrangian relaxation algorithm in Sect. 6.4 suggests little room for improving the optimality gaps. In this section, we pursue a different approach, specifically aimed at achieving fast bounds, albeit of a lower quality. To this end, we also analyze the advantages and disadvantages of using the PH heuristic presented in Sect. 4 and Algorithm 2. As we describe in Sect. 4, PH offers theoretical difficulties, however they are often marred by the algorithm's practical performance. This makes PH a powerful heuristic for practical problems; similar trends have been reported earlier [51,61].

Choice of
Algorithm 2 is known to be sensitive to the choice of the parameter ρ [16,61]; no consensus seems to exist in the literature for deciding optimal values of ρ. We compare the sensitivity of our results using three values of ρ t , in addition to the ρ t = 1 results reported above. In the first two runs, we use constant ρ t values of 0.1 and 10, and in the third run we use a coefficient specific ρ value described in [61]. This value is ρ t = R t max{1, ω∈Ω (y t −y ω in Step 7 faster. This is similar to the trend observed in [16]. However, ρ t = 0.1 results in weaker upper bounds than ρ t = 1. The variable specific ρ closes the gap fastest, however the improvements compared to ρ t = 10 are minimal. Further, none of the ρ values completely close the gap; and, ρ t = 1 provides the best upper bounds.

Linearization
The quadratic proximal term in the PH objective can be linearized using piece-wise linear segments. While, there are many sophisticated techniques for this linearization [22,58], we attempt a basic first-order Taylor expansion. Dropping the t subscript, with L − 1 segments, we have: y 2 ≈ a 2 + 2a(y − a), where a = (l−1)y max L−1 , l = 1, 2, . . . , L. However, in our experiments, we observe that to achieve solutions comparable to that of the original mixedinteger quadratic constrained program (MIQCP), a large value for L -close to 100 -is needed. Then, the time required to solve the larger MIPs is no significantly different from that of the MIQCP. The reason for this is as follows.
There is only a single binary variable for each subproblem in Step 6 of Algorithm 2; thus, each MIP can be solved in at most two LP iterations (by branching on z ω = 0 and z ω = 1). When L increases, the time required for an LP solution increases due to the increased number of constraints in the LP. Poor values of y reported by PH cascade to poor values for z in Step 9 of the algorithm, resulting in poor values for γ , and ultimately poor performance for the Lagrangian dual. Two possible suggestions for the interested reader, that we do not investigate, are as follows. First, use a dynamic linearization as described in Veliz et al. [58] that gradually increments the number of piecewise segments. Second, use a subgradient step-size update rule that does not depend on γ .

Convergence of algorithm 2: Bundling
The number of PH iterations required for each subgradient iteration, iter, can be large to force PH to converge. Although each iteration involves only a single binary variable, this can make the PH algorithm computationally inefficient. In our computational experiments, PH does not converge using the criteria in Step 7 of Algorithm 2, even with iter = 100. As the number of scenarios grows, the time to solve the MIQCPs in Step 6 grows as well; solving MIQCPs can be significantly more challenging than solving MIPs. However, nearly always the best upper bound, P H U B , is obtained in the first iteration itself. To this end, we begin our discussion by reporting results obtained by setting iter = 1; i.e., without even entering the while loop of Algorithm 2. Table 6 reports these results. The PH upper bounds are on average 10% (and, at most 11.7%) away from the optimal Lagrangian values (the fourth column of Table 5). For instances that are solved optimally in Table 2, the P H U B values are on average 10.9% (and, at most 11.7%) away from the true optimal solutions. We obtain these very quickly, see column "Time" in Table 6; the 1200 scenario instances are solved in under half a minute.
However, if Step 4 of Algorithm 1 is replaced with Algorithm 2 without the while loop, the quality of the Lagrangian relaxation scheme deteriorates significantly due to the cascading effect we describe in Sect. 6.5.2. Again, there is a tradeoff between computational effort required (a large value for iter) versus the quality of the solution (a value close to zero for γ ). To this end, we investigate the following enhancement in the PH algorithm.
Instead of solving individual scenario sub-problems, we bundle together several scenarios and solve this "mini-extensive" form problem as our new sub-problem. Bundling scenarios  Table 7 Computational results for model (2) using Algorithm 2 with Algorithm 1 using bundling. Here, "size" indicates the size of a bundle, "LR" indicates the Lagrangian relaxation bound obtained using the PH algorithm, "iterations" indicates the number of completed LR iterations, and "Gap" is the relative gap between the LB column of could offer the advantage of speeding convergence of PH, but with increased computational effort per PH iteration. For details on implementation of bundling in progressive hedging, see, e.g. [16]. Further, the bounds obtained in Step 6 of Algorithm 2 also extend to the case with bundling [16]. Finally, we emphasize again that since the PH algorithm is a heuristic for solving model (4), the bounds obtained would be worse than those obtained using a naive solution method (Table 5); however, we are interested in the tradeoff between speed and quality. Table 7 compares the results of our experiments using two different bundles. We use Algorithm 1 with Step 3 substituted by Algorithm 2. Further, to see the tradeoff between the computational time and the value of the optimal solution, we set a maximum time limit of 900 seconds. When the completion of a LR iteration would have resulted in the total time Table 8 Comparison of different upper bounds for model (2). "LP" denotes the linear programming relaxation, "LR" denotes the Lagrangian relaxation of model (4), and "LR-PH" denotes the Lagrangian relaxation using progressive hedging algorithm with bundling. The "Gap" columns denote the optimality gaps from max{z * , L B}; see, Fig. 1 50. In all of the 12 instances, the 150 size bundles lead to better (lower) Lagrangian bounds than the 100 size bundles. However, there is no consistent trend in the time required for completion. For example, the 450 scenario instance with ε = 0.05 took lesser time for the 150 bundle case than the 100 bundle case; but, the effect is the opposite for the 450 scenario instance with ε = 0.01.

Summary of bounds
From the preceding discussions, we have two valid upper bounds for model (2) available from plain LR and LR with PH. A third bound is obtained by solving the LP relaxation of model (2). We summarize these bounds in Table 8. The "Gap" here is the optimality gap of the bound from max{z * , LB}. The LP has a gap of 5.9% on average. From the LR column, we note that four of the twelve instances are solved to optimality-but this optimality is not verifiable (see Table 5). Imposing a reduced time limit and solving the LR with the bundled PH algorithm increases this optimality gap by at most 1.5% (with bundles of 150 scenarios). This again indicates that the PH algorithm achieves good bounds in fairly early iterations.

Background
In this section, we present a proof-of-concept for a generalization of the models and analysis presented in the preceding sections of this article. To this end, we consider a microgrid with several photovoltaic units and two storage units. To build our test case, we consider a single-phase simplification of the IEEE 13-bus feeder [25,53]. We adopt a network flow approximation to model the power flow; this is equivalent to the commonly used "DC approximation" of power flow as the network considered is acyclic or radial, as is common for distribution systems. While this is a loose approximation [11], our aim is to demonstrate how the techniques described above can be applied in a different context. We modify the IEEE 13-bus feeder system as follows. We place PV units at these five buses: 652, 680, 675, 611, 634. We assume the solar power available from the PV units at these buses are, respectively, 1, 2, 3, 4, 5 times 1 5 s ω t , where s ω t is the solar power data used in the rest of this article. We also place two storage devices: one at bus 632 and another at bus 671. The interface to the larger grid is at bus 650, which we denote at b gr , and we buy and sell energy from this bus. Figure 4 gives a topological representation of the modified IEEE 13-bus feeder, with the placements of PV and battery storage devices. Further, we reduce the three-phase system to a single phase and set the line capacities, F l , to eight times the values in the IEEE 13-bus feeder system. The system is infeasible for lower values of the thermal line limit, and indeed for the computational instances we present in Sect. 7.3 this thermal limit is almost reached. Further, since this system does not provide load values indexed by time, we scale the load by computing load factors from [24, Table 3], to achieve a temporal resolution. Finally, we set X s to ten times the value in the rest of this article for each storage device s.

Conclusion
We developed a stochastic optimization model to schedule a standalone hybrid system consisting of battery storage with solar power. We used a chance-constrained formulation to ensure highly reliable operations under the uncertainty of solar power. The model makes charging/discharging decisions dictated by a mix of economics and the available solar power. For a few hundred scenarios, the model is relatively tractable benefiting from a strong LP relaxation. For larger number of scenarios, we demonstrate the applicability of a simple Lagrangian relaxation procedure. The procedure is remarkably effective, and tightens the optimality gap when a naive solution method does not. Finally, we present a heuristic, based on the progressive hedging algorithm, that we couple with the Lagrangian relaxation. We find significant benefit when scenarios are bundled and solved as one in the progressive hedging iterations; however, when solved plainly the PH heuristic does not achieve good quality solutions. We conclude with an extension to a larger system with multiple storage devices and batteries connected to form a microgrid. Future work could examine interconnected updates of the subgradient algorithm with the progressive hedging algorithm. The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request. regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.