A Scalable Stochastic Programming Approach
for the Design of Flexible Systems

For a given sample $\boldsymbol{\theta}^{k}$, problem (1.3) can be written as:

where $U\in\mathbb{R}_{+}$ is a sufficiently large constant. Note that the continuous variable for a particular sample $u^{k}$ is equivalently replaced with $y^{k}U$. The binary variable $y^{k}$ takes a value of zero if the realization $\boldsymbol{\theta}^{k}$ is feasible (corresponding to $\psi(\mathbf{d},\boldsymbol{\theta}^{k})=0$) or takes a value of one if it is infeasible (corresponding to $\psi(\mathbf{d},\boldsymbol{\theta}^{k})>0$).

By combining all realizations $k\in K$, we obtain the aggregated problem:

Since $\mathbf{d}$ is fixed, this problem is fully decoupled in the set $K$ (i.e., solving Problem (2.11) is equivalent to solving Problem (2.10) for each $k\in K$) and the optimal values of the binary variables $y^{k}$ satisfy $SF_{K}(\mathbf{d})=\frac{1}{|K|}\sum_{k\in K}(1-y^{k})$. Consequently, (2.11) can be used to compute the sample average approximation of the $SF$ index and the approximation becomes exact as $|K|\rightarrow\infty$.

The constant $U$ in the objective can be eliminated without affecting the solution. Moreover, we obtain an equivalent problem by maximizing $SF_{K}(\mathbf{d})=\frac{1}{|K|}\sum_{k\in K}(1-y^{k})$ directly. These modifications give the problem:

An important observation is that a continuous relaxation of problem (2.11) is given by:

This provides a relaxation because $u^{k}$ can be modeled as $u^{k}=y^{k}U$ with $0\leq y^{k}\leq 1$ and sufficiently large $U\in\mathbb{R}_{+}$. In summary, the continuous problem (2.13) is a relaxation of the mixed-integer problem (2.11) (or, equivalently, (2.12)). We also note that $u_{k}=\psi(\mathbf{d},\boldsymbol{\theta}^{k})$ is a measure of infeasibility for sample $\boldsymbol{\theta}^{k}$ and thus the continuous problem (2.13) minimizes the mean (expected) infeasibility. In other words, the expected infeasibility is a surrogate measure of flexibility.

We now proceed to find a design $\mathbf{d}$ that minimizes the design cost and that maximizes the stochastic flexibility index. A sample average approximation of this problem is given by:

A Pareto solution of this problem can be obtained by solving an $\epsilon$-constrained problem of the form:

where $\epsilon_{s}$ is the $\epsilon$-parameter. This formulation is precisely a sample average approximation of the joint chance-constrained problem (1). This problem is computationally expensive to solve due to its mixed-integer nature and due to the large number of potential MC samples.

A key observation is that Pareto solutions for the conflict resolution problem can also be obtained by solving the mixed-integer problem:

Note that this is simply the counterpart to (2.15). We denote the solution of this problem as $\mathbf{d}^{*},{y^{k}}^{*}$ and we have that $c(\mathbf{d}^{*})=\epsilon_{f}$ holds at the solution (because the objective is conflicting) and the Pareto pair is given by $(c(\mathbf{d}^{*}),SF_{K})=(\epsilon_{f},SF_{K})$, where $SF_{K}=\frac{1}{|K|}\sum_{k\in K}(1-{y^{k}}^{*})$.

A continuous relaxation of this problem is given by:

The solution of this problem delivers a continuous solution $\bar{y}^{k}$ from which we recover integer values as $\bar{y}^{k}\leftarrow\mathbbm{1}_{\bar{y}^{k}}$ (i.e., we recover $\bar{y}^{k}=0$ if $k$ is feasible or $\bar{y}^{k}=1$ otherwise). With this, we compute an estimate of the SF index as $\bar{SF}_{K}=\frac{1}{|K|}\sum_{k\in K}(1-\mathbbm{1}_{\bar{y}^{k}})$. We solve (2.17) again (by adjusting $\mathbf{d}\in\mathcal{D}$ and $\mathbf{z}^{k}$) by fixing the binary variables to ensure that the selection of active constraints provides a feasible solution. Note that the objective function $\bar{SF}_{K}$ is fixed since the binaries are fixed. The solution of this feasibility problem is given by $\bar{\mathbf{d}},\bar{\mathbf{z}}^{k}$. We highlight that problem (2.17) can also be solved by using the continuous variables $u^{k}$. We will see that this approach delivers Pareto pairs $(c(\bar{\mathbf{d}}),\bar{SF}_{K})$ that closely match the optimal Pareto pairs $(c({\mathbf{d}^{*}}),{SF}_{K})$ of the mixed-integer formulation. We also note an approximation of comparable quality is achieved when the counterpart problem of (2.17) is solved.

A precise theoretical justification for this continuous relaxation providing such a high quality approximation is the subject of ongoing research and is beyond the scope of this work. We hypothesize that the high quality of the solutions is due to the degenerate nature of the SF index. Specifically, different combinations of active constraints give the same or a very similar SF index. This behavior has been recently reported in [3], where the authors note that chance constraints give rise to a wide range of local minima with similar values. This behavior becomes more evident when one moves the SF index into the objective (as opposed to imposing it in the constraints). This is precisely why it is important to think about the design problem as a conflict resolution problem. Recent work has also provided evidence that continuous (Lagrangian) relaxations of chance-constrained problems provide high quality approximations [1].

The three-node distribution network features a centralized supplier configuration. Figure 1 details this network and provides the arc and supplier capacities. The network is subjected to multivariate Gaussian demands $\boldsymbol{\theta}=(r_{1},r_{2},r_{3})\sim\mathcal{N}(\bar{\boldsymbol{\theta}},V_{\boldsymbol{\theta}})$. The mean $\bar{\boldsymbol{\theta}}$ is taken to be $\bar{\boldsymbol{\theta}}=(0.0,60.0,10.0)$, and the covariance matrix $V_{\boldsymbol{\theta}}$ is assumed to be:

We also set the parameter $U=10000$.

The IEEE 14-node network was originally provided in Dabbagchi in [4]. The system data is obtained from MATPOWER. This test case does not provide arc capacities so we enforce $a^{C}=100$ for all the arcs. A schematic of this system is provided in Figure 2. This system is subjected to a total of 11 uncertain parameters (the network demands). The demands are assumed to be $\boldsymbol{\theta}\sim\mathcal{N}(\bar{\boldsymbol{\theta}},V_{\boldsymbol{\theta}})$, where $\bar{\boldsymbol{\theta}}=(87.3,50.0,25.0,28.8,50.0,25.0,\allowbreak 0,0,0,0,0)$ and $V_{\boldsymbol{\theta}}$ is symmetric matrix with $(V_{\boldsymbol{\theta}})_{ii}=1200$ and $(V_{\boldsymbol{\theta}})_{ij}=240,\ \forall i\neq j$. Also, we set $U=10000$.

The 141-node power distribution network corresponds to an urban area in Caracas, Venezuela and was originally developed in [10]. The network data is extracted from MATPOWER, but again no arc capacities are provided (we assume $a^{C}=100$). Figure 3 provides a schematic of the network. This system is subjected to 84 uncertain disturbances corresponding to the demands. The demands are assumed to be $\boldsymbol{\theta}\sim\mathcal{N}(\bar{\boldsymbol{\theta}},V_{\boldsymbol{\theta}})$, where $\bar{\boldsymbol{\theta}}=\bar{\boldsymbol{\theta}}_{fc}$ and $V_{\boldsymbol{\theta}}=100\mathbb{I}$. The point $\boldsymbol{\theta}_{fc}$ is the feasible center and is the instance of $\boldsymbol{\theta}\in\Theta(\mathbf{d})$ that minimizes the feasibility function $\psi(\mathbf{d},\boldsymbol{\theta})$ [15]. We set the upper bounding constant to $U=10000$.

For each network, the cost function $c(\mathbf{d})$ is defined to be linear of the form:

Here the unit cost for each design variable is taken to be $1/\sqrt{n_{d}}$, where $n_{d}=|\mathcal{S}|+|\mathcal{A}|$. Furthermore, a relatively large value of $U$ is selected in the above problems to ensure a sufficient amount of slack is provided for the inequality constraints. This is all done for convenience in conducting the analysis below.

Combining (2.16) with the model equations in (3.18c), we obtain:

We substitute the node balances into the capacity constraints to obtain a system of inequalities that are expressed solely in terms of $\boldsymbol{\theta}$. Thus, the three-node network solutions can be visualized in three dimensions. Figure 5 shows three solutions corresponding to $\epsilon_{c}=\{0,2.75,10.5\}$. Here, we observe how the diagonal plane (on the right side of the figures), which corresponds the supplier capacity constraint, is shifted as the value of $\epsilon_{c}$ is increased to minimize the number of infeasible instances and thus maximize the $SF_{K}$ index.

The same analysis is done using the IEEE-14 network with 2,000 MC samples. This creates a MILP with 50,025 continuous variables, 2,000 binary variables, and 128,001 constraints. The Pareto set is created by varying the value of $\epsilon_{c}$ from 0 to 65 in increments of 2.5. All of the corresponding numerical results are given in Table 4 in the appendix. Figure 6 shows the Pareto set corresponding to these results. Here, the minimum $SF_{K}$ index associated with the base case is 88.10% and the maximum $SF_{K}$ index is 91.50%. This means that, regardless of how much the line and supplier capacities are increased, we can only improve $SF_{K}$ by 3.4% (relative to this MC sample set). This behavior occurs because in certain sampled instances a number of demands are sufficiently negative such that this is a surplus in the system that cannot be rectified since there are no sinks in the network. Such behavior is not readily obvious and thus highlights the utility of this optimal design framework in understanding how to promote the flexibility of a system and in determining to what extent its flexibility can be improved. We also observe that the average solution time for the MILP problems is 40.47s.

Problem (3.21) is also applied to the 141-node power distribution network using 10,000 MC samples and is evaluated for each value of $\epsilon_{c}$ from 0 to 300 in increments of 10. Each MILP contains 1,410,141 continuous variables, 10,000 binary variables, and 4,230,001 constraints. A time limit of 3,600s is also imposed. All 30 results are summarized in Table 5 in the Appendix and Figure 7 shows these results. Here, a much more appreciable range of $SF_{K}$ index values is obtained where the minimum $SF$ index of the system is 35.60% and the maximum $SF$ index is 70.64%. This limited maximum $SF$ index can be attributed to topological limitations in the network. Specifically, the sampled demands are Gaussian random variables and thus can be negative which in certain cases means the network is confronted with excess power production that it cannot shed. Thus, we observe that in this case increasing the network capacities alone cannot provide a $SF$ index near 100% due to topographical limitations. Instead, the system topology would need to be modified to enhance its flexibility (e.g., add a sink to shed excess power production). This highlights how this framework is useful in guiding system design by demonstrating to what degree the system flexibility can be bettered via continuous design variables such as capacity.

A key observation is that only 12 of the 30 problems converged within the time limit. This unfavorable behavior is likely due to the large-U constraints which can lead to weak linear program (LP) relaxations. This scalability limitation clearly demonstrates that Problem (2.16) can become impractical for moderate and large sized systems. We note some advanced methods such as those proposed in [22, 17] can be used to select smaller $U$ values that might strengthen the LP relaxations. However, our proposed continuous relaxation entails a more efficient strategy since it does not require the solution of a MIP. Another observation drawn from the solution times in Table 5 is that Pareto pairs near minimum or maximum design cost require significantly less computational time relative to the intermediate pairs. This trend is also prevalent in the computational results of the IEEE-14 network and the three-node network.

We compare the Pareto pairs of the mixed-integer formulation to those obtained with the continuous formulation:

The $\bar{SF}_{K}$ index is determined after solving this problem via the indicator function $\mathbbm{1}_{y^{k}}$ as described in Section 2.3. The three-node network is again analyzed under the same conditions described in Section 3.2 (i.e., the same samples and $\epsilon_{c}$ values). For the small 3-node network, the problem is an LP with 4,003 variables and 9,001 constraints. Each solution is also verified to be integer feasible by ensuring the formulation is feasible with the values of the $y^{k}$ variables fixed to their indicated values as described in Section 2.3. The numerical results are provided in Tables 2 and 3 in the Appendix and Figure 8 juxtaposes these results with those obtained with the mixed-integer formulation. Here, we note that the Pareto pairs are equivalent and that the mixed-integer $y^{k}$ values exactly match the indicated $y^{k}$ values obtained with the continuous formulation (i.e., the same active constraints are identified). Also, each continuous formulation solution only required 0.0049s on average to converge which equates to a 74% reduction in computational cost.

We also applied Problem (3.22) to the IEEE-14 node power network with the same samples and conditions used in Section 3.2. This is an LP with 52,025 variables and 128,001 constraints. The indicated values of $y^{k}$ are all verified to be integer feasible. These results are detailed in Table 4 in the Appendix and Figure 9 compares the results with those obtained with the mixed-integer formulation. Again, we note that the continuous formulation is able to recover the same Pareto set as the mixed-integer counterpart. Only the last Pareto pair differs (by one $y^{k}$ out of 2,000). All other pairs are identical and identify exactly the same sets of active constraints. Furthermore, the continuous formulation only requires 1.74s on average to solve, which corresponds to a 96% reduction in computational cost relative to the mixed-integer formulation.

Finally, we applied Problem (3.22) to the 141-node power network under the same conditions detailed in Section 3.2. This creates a large-scale LP with 1,420,141 variables and 4,230,001 constraints. The indicated values of $y^{k}$ are again all verified to be integer feasible. These results are provided in Table 5 in the Appendix, and Figure 10 shows the Pareto pairs obtained from both formulations. In this case, the Pareto solutions are very similar but differences are more perceptible but rather small. Specifically, the identified values of $y^{k}$ differ by 0.49% on average relative to the optimal values obtained via the mixed-integer formulation (i.e., 49 differences in the set of $y^{k}$ values out of the 10,000). The continuous Pareto pairs only required 13.41s on average to solve (in contrast to the mixed-integer results for which the majority were unable to converge within the 3,600s time limit).

Our results highlight that the continuous optimal design formulation provides high quality solutions with significant reductions in computational time. As mentioned in Section 2.3, the high quality of the solutions can likely be attributed to the degenerate nature of the SF index.