A Computational Framework for Quantifying and
Analyzing System Flexibility

Problem (1.7) can use any compact set $T(\delta)$ whose size can be parameterized in terms of a scalar variable $\delta$. The selection of $T(\delta)$ will depend on the application and on the nature of the uncertain parameters $\boldsymbol{\theta}$. In [22], a number of uncertainty sets are provided that can be used in conjunction with Problem (1.7). Table 1 summarizes such sets. We note that the ellipsoidal norm set in Table 1 is equivalent to $T_{ellip}(\delta)$ as defined in (1.9) provided that $A=V_{\boldsymbol{\theta}}^{-1}$. The ellipsoidal norm set can also be used to approximate polytopes [2]. We also note that $T_{box}(\delta)$ is equivalent to $T_{\infty}(\delta)$ when $\Delta\boldsymbol{\theta}_{i}^{-}=\Delta\boldsymbol{\theta}_{i}^{+}=1$. Li et.al. provide an analysis of the $\ell_{1}$, $\ell_{2}$, and $\ell_{\infty}$-norm uncertainty sets along with their mathematical properties in the context of robust optimization [16]. The conditional value at risk (CVaR) norm, originally proposed by Pavlikov and Uryasev in [18], has gained interest in the robust and stochastic optimization communities because of its ability to approximate $\ell_{p}$ norms via linear programming (precisely reproducing the $\ell_{1}$ and $\ell_{\infty}$ norms as extreme cases) [9].

Li et. al. also analyzed sets that result from intersections (compositions) of $\ell_{p}$-norm sets (e.g., $T_{2}(\delta)\cap T_{\infty}(\delta)$) [16]. Inspired by this, we observe that $T(\delta)$ can be represented by any combination of the uncertainty sets in Table 1 or other sets (provided that the resulting combined set is compact). This provides the ability to incorporate physical knowledge on the uncertain parameters and/or to create a wider range of uncertainty set shapes. For instance, the positive (truncated) ellipsoidal set $T_{ellip+}(\delta):=T_{ellip}(\delta)\cap\mathbb{R}_{+}^{n_{\theta}}$, where $\mathbb{R}_{+}^{n_{\theta}}:=\{\boldsymbol{\theta}:\boldsymbol{\theta}\geq 0\}$ can be used to represent parameters that are known to be non-negative (e.g., demands, prices). Similarly, we can define the positive hyperbox set $T_{box+}(\delta):=T_{box}(\delta)\cap\mathbb{R}_{+}^{n_{\theta}}$. Compositions of sets can be naturally accommodated in conjunction with problem (1.7), by imposing the constraints corresponding to each set. For instance, $T_{ellip+}(\delta)$ is represented by adding the constraints $(\boldsymbol{\theta}-\bar{\boldsymbol{\theta}})^{T}V_{\boldsymbol{\theta}}^{-1}(\boldsymbol{\theta}-\bar{\boldsymbol{\theta}})\leq\delta$ and $\boldsymbol{\theta}\geq 0$ to the MIP formulation.

Most uncertainty sets require a nominal point $\bar{\boldsymbol{\theta}}$. The choice of the nominal point $\bar{\boldsymbol{\theta}}$ is critical; for example, a nominal point that is close to the boundary of the feasible region will have limited flexibility. In some applications, specifying $\bar{\boldsymbol{\theta}}$ might be difficult if not enough data is available or if a system is not routinely operated at any given point (e.g., a power grid). Thus, we introduce methods that can be employed to find well-centered nominal points. We note that one would be tempted to let $\bar{\boldsymbol{\theta}}$ be a free variable in Problem (1.7) to find a nominal point but this approach leads to trivial solutions. This is because the objective is minimized when the nominal point $\bar{\boldsymbol{\theta}}$ is placed on the boundary of $\boldsymbol{\theta}$. Consequently, this naive approach cannot be used to compute a nominal point.

We propose an extension of the analytic center (commonly used in interior-point methods) to compute a nominal point. The analytic center $\bar{\boldsymbol{\theta}}_{ac}$ is given by:

		$\displaystyle{\bar{\boldsymbol{\theta}}}_{ac}:=$	$\displaystyle\mathop{\mbox{argmax}}_{\boldsymbol{\theta},\mathbf{z},\mathbf{x},s_{j}}$			$\displaystyle\sum_{j\in J}\log\left(s_{j}\right)$		(2.11)
	s.t.			$\displaystyle f_{j}(\mathbf{z},\mathbf{x},\boldsymbol{\theta})+s_{j}\leq 0,$		$\displaystyle j\in J$
		$\displaystyle h_{i}(\mathbf{z},\mathbf{x},\boldsymbol{\theta})=0,$		$\displaystyle i\in I.$

Here, $s_{j}\in\mathbb{R}$ are slack variables. This problem pushes the constraints to the interior of the feasible region. Any bounded solution of the above problem satisfies $s_{j}>0$ and thus $f_{j}(\mathbf{z},\mathbf{x},\bar{\boldsymbol{\theta}}_{ac})<0$ for all $j\in J$. Consequently, we have that $\psi(\bar{\boldsymbol{\theta}}_{ac})<0$ and thus $\bar{\boldsymbol{\theta}}_{ac}\in\Theta$ (the analytic center is feasible). We also note that the analytic center is the point that maximizes the geometric mean of the constraints (slack variables). By changing the objective to $\sum_{j\in J}s_{j}$ one finds a center point that maximizes the arithmetic mean.

We also propose a center point $\bar{\boldsymbol{\theta}}_{fc}$ (that we call the feasible center) and that we define as:

		$\displaystyle\bar{\boldsymbol{\theta}}_{fc}:=$	$\displaystyle\mathop{\mbox{argmax}}_{\boldsymbol{\theta},\mathbf{z},\mathbf{x},s}$			$\displaystyle s$		(2.12)
	s.t.			$\displaystyle f_{j}(\mathbf{z},\mathbf{x},\boldsymbol{\theta})+s\leq 0,$	$\displaystyle j\in J$
		$\displaystyle h_{i}(\mathbf{z},\mathbf{x},\boldsymbol{\theta})=0,$	$\displaystyle i\in I.$

Here, $s\in\mathbb{R}$ is a slack variable. The feasible center is the point that maximizes the worst-case slack variable (as opposed to the geometric mean used in the analytic center). We also note that $\bar{\boldsymbol{\theta}}_{fc}=\mathop{\mbox{argmin}}_{\boldsymbol{\theta}}\;\psi(\boldsymbol{\theta})=\mathop{\mbox{argmax}}_{\boldsymbol{\theta}}\;-\psi(\boldsymbol{\theta})$ and that $\psi(\bar{\boldsymbol{\theta}}_{fc})\leq 0$. Moreover, (2.12) is equivalent to (1.5) when the parameter $\theta$ is set as a free variable.

The feasible center is an approximation of the center point:

$\displaystyle\bar{\boldsymbol{\theta}}_{d}:=\mathop{\mbox{argmax}}_{\boldsymbol{\theta}\in\Theta}\;d(\boldsymbol{\theta},\partial\Theta)$

(2.13)

where $d(\boldsymbol{\theta},\partial\Theta)=\min_{\boldsymbol{\theta}^{\prime}\in\partial\Theta}\;\|\boldsymbol{\theta}-\boldsymbol{\theta}^{\prime}\|$. In other words, $\bar{\boldsymbol{\theta}}_{d}$ is the point that maximizes the distance to the boundary of feasible set. This point is a natural center point but it is difficult to compute due to the nested max and min operators. We now proceed to show that $|-\psi(\boldsymbol{\theta})|$ provides a lower bound of $d(\boldsymbol{\theta},\partial\Theta)$. Consequently, because the feasible center $\bar{\boldsymbol{\theta}}_{fc}$ maximizes $-\psi(\boldsymbol{\theta})$, we have that this also implicitly maximizes $d(\boldsymbol{\theta},\partial\Theta)$. This property can be established as follows, if the constraints are Lipschitz continuous (as often assumed in the flexibility analysis literature), we have that $\psi(\cdot)$ is Lipschitz as well and thus:

$\displaystyle|\psi(\hat{\boldsymbol{\theta}})-\psi(\boldsymbol{\theta})|\leq L\|\bar{\boldsymbol{\theta}}_{d}-\boldsymbol{\theta}\|$

(2.14)

where $\hat{\boldsymbol{\theta}}=\mathop{\mbox{argmin}}_{\boldsymbol{\theta}^{\prime}\in\partial\Theta}\;\|\boldsymbol{\theta}-\boldsymbol{\theta}^{\prime}\|$ for a particular value of $\boldsymbol{\theta}$. Because $\hat{\boldsymbol{\theta}}\in\partial\Theta$, we have that $\psi(\hat{\boldsymbol{\theta}})=0$ and thus,

	$\displaystyle|\psi(\hat{\boldsymbol{\theta}})-\psi(\boldsymbol{\theta})|$	$\displaystyle=|-\psi(\boldsymbol{\theta})|$
		$\displaystyle\leq L\|\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}\|$
		$\displaystyle=L\min_{\boldsymbol{\theta}^{\prime}\in\partial\Theta}\;\|\boldsymbol{\theta}^{\prime}-\boldsymbol{\theta}\|$
		$\displaystyle=L\,d(\boldsymbol{\theta},\partial\Theta).$		(2.15)

We thus conclude that:

$\displaystyle d(\boldsymbol{\theta},\partial\Theta)\geq\frac{1}{L}|-\psi(\boldsymbol{\theta})|.$

(2.16)

Consequently, maximizing $-\psi(\boldsymbol{\theta})$ implicitly maximizes $d(\boldsymbol{\theta},\partial\Theta)$. Moreover, we note that $d(\boldsymbol{\theta},\partial\Theta)=-\psi(\boldsymbol{\theta})=0$ if $\boldsymbol{\theta}\in\partial\Theta$.

The flexibility index is a valuable metric that can be used to compare system designs. The utility and interpretation of the index for a system will depend on the choice of $T(\delta)$. For instance, the index can be used to determine the confidence level $\alpha^{*}$ if the set $T_{ellip}(\delta)$ is used (when $\boldsymbol{\theta}~{}\sim\mathcal{N}(\bar{\boldsymbol{\theta}},V_{\boldsymbol{\theta}})$). The index might not have a clear interpretation or utility for some choices of $T(\delta)$ but can still be useful to quantify improvements in flexibility from design modifications (retrofits). We illustrate this feature by using the following example.

Consider two system designs that are each subjected to two normal random variables and that have no recourse variables $\bf{z}$ as summarized in Table 2.

$$\bar{\boldsymbol{\theta}}=\begin{bmatrix}4\\ 5\end{bmatrix}\ \ \ \ \ \ V_{\boldsymbol{\theta}}=\begin{bmatrix}2&1\\ 1&3\end{bmatrix}.$$

(2.17)

Figure 2 depicts the two different systems. It is apparent that Design B has a larger feasible region $\Theta$ and thus has a larger $SF$ index. A key observation is that the choice of $T(\delta)$ significantly affects how it measures system flexibility. For instance, for Design B, the ellipsoidal and hyperbox sets identify different limiting constraints, and thus exhibit distinct behavior. In particular, index $F_{box}$ is limited by constraint $f_{2}$ and would not change if constraints $f_{1}$, $f_{3}$, and/or $f_{4}$ were varied to increase the area of $\Theta$, even though the overall system flexibility would be improved. On the other hand, index $F_{ellip}$ is limited by constraint $f_{1}$. This conservative behavior parallels what is observed with the use of uncertainty sets in robust optimization, where such sets are used to optimize against the worst case scenario such that the solution is robust in the face of uncertainty, and the conservativeness of a solution depends on the chosen shape of the uncertainty set [8].

An important observation that we make is that the flexibility index problem also implicitly identifies inequality constraints that limit flexibility. Such information can be valuable in design or operating procedures. Specifically, at a solution of Problem (1.7), the binary variables $y_{j}$ indicate which inequality constraints are active and therefore limit flexibility. Problem (1.7) can thus be resolved by excluding this first set of limiting constraints, so that the next set of limiting constraints can be identified. This step can be repeated to rank system constraints to a desired extent (as long as the solution of Problem (1.7) is bounded). This is done by noticing that each set of limiting constraints will have an associated flexibility index, which indicates how limiting a particular set of constraints is relative to the other sets. We note that this analysis operates under the assumption that the inequality constraints can be eliminated/relaxed such that the system still has physical meaning. However, this assumption is met in many applications since constraints that cannot be removed (e.g., physical laws) typically correspond to equality constraints.

This methodology can also be used to identify and rank system components that limit flexibility in order to guide design improvements or quantify value of specific components, since such identification is useful to understand and identify system vulnerabilities. Thus, this ranking becomes useful for complex systems where it is difficult to determine vulnerable components. For example, this approach can identify which heat-exchangers in a heat-exchanger network could be retrofitted to most increase the system flexibility, and similarly identify heat-exchangers that would not lead to increased flexibility when improved. Limiting constraint information can also be used to quantify the impact of failure of system components (e.g., a production facility).

This proposed approach shares some common ground with optimal design problems where design variables $\mathbf{d}$ (added to the system constraints $f_{j}(\mathbf{z},\mathbf{x},\boldsymbol{\theta})$ and $h_{i}(\mathbf{z},\mathbf{x},\boldsymbol{\theta})$) are selected to minimize cost and satisfy a specified index $F$ [1, 13, 21]. This provides a straightforward approach for selecting a system design that is sufficiently flexible prior to its determination. However, our approach differs in that it can be used to fully characterize how each component of a particular system design impacts system flexibility instead of choosing design variables to improve flexibility to desired extent. This can becomes useful for intelligently retrofitting existing systems and for understanding which components most limit flexibility (i.e., are most vulnerable to uncertainties) that therefore warrant close monitoring.

We illustrate the constraint ranking methodology using Design A (described in Section 2.3). Here, we choose the set $T_{ellip}(\delta)$. Problem (1.7) is solved to identify the limiting constraints, which are then eliminated to solve the problem again to identify a new set of limiting constraints. In this case, the problem is solved four times as only one constraint becomes active at the solution of each problem. The results are summarized in Table 4. We observe that constraint $f_{1}$ limits the flexibility index the most. Specifically, the value of $F_{ellip}$ corresponding to the next limiting constraint $f_{2}$ is 79.3% larger. Diminishing changes in $F_{ellip}$ are obtained for subsequent limiting constraints, showing that constraints $f_{2}$-$f_{4}$ have little impact on the system flexibility relative to $f_{1}$.

Figure 3 depicts the solutions corresponding to the results presented in Table 4. It is apparent that the uncertainty set is larger after each subsequent ranking step. We again make the observation that the limiting constraints are determined by the shape of the uncertainty set. Thus, different uncertainty sets can yield distinct classifications of limiting constraints.

We now consider the IEEE 14-node power network, 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. This base design is labeled Design 1 and a schematic of this system is provided in Figure 9. Design 2 is obtained by removing the arc that connects node 2 to node 5 and by removing the arc that connects node 9 to node 14. Design 3 uses a centralized power supply scheme, this is accomplished by removing all of the suppliers except for the two located on nodes 1 and 2, by increasing the capacity of each remaining supplier from 332 to 432 and from 140 to 340 respectively, and by increasing the capacity of each arc to 200.

This system is subjected to a total of 10 uncertainty disturbances (the network 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}=(87.3,50.0,25.0,28.8,50.0,25.0,\allowbreak 0,0,0,0,0)$ and $V_{\boldsymbol{\theta}}=1200\mathbb{I}$. We select the sets $T_{box+}(\delta)$ and $T_{ellip+}(\delta)$ (the demands are nonnegative). Each element of the hyperbox deviations $\Delta\boldsymbol{\theta}^{-},\Delta\boldsymbol{\theta}^{+}$ is set to $3\sigma_{i}=103.92$ (corresponding to $\bar{\boldsymbol{\theta}}_{i}\pm 3\sigma_{i}$ confidence bounds). Both of these sets are used to determine the flexibility index from Problem (1.7) for each design. The stochastic flexibility index $SF$ is computed using 1,000,000 MC samples.

Table 9 summarizes these design comparison results. Index $SF$ shows that the flexibility of Design 1 is worsened with the removal of the two arcs, whereas the centralized supply modification improves flexibility. An observation is that $F_{ellip+}$ mirrors the behavior $SF$ at a much lower computational cost. Specifically, $F_{ellip+}$ required less than 20 seconds to compute while computing $SF$ using MC sampling required 13,030 seconds (3.6 hours). Index $F_{box+}$ is not able to mirror the same behavior (it does not indicate that Design 2 is less flexible than Design 1 as would be expected). This analysis also shows that the higher capacity of the centralized supply scheme (Design 3) increases system flexibility relative to Design 1, which again is counter-intuitive.

The limiting constraints of Design 1 are identified and ranked using $T_{ellip+}(\delta)$. The ranking results are shown in Table 10. The capacity constraints corresponding to four of the five suppliers along with three arc capacity constraints are identified as the most limiting ones. Again, this indicates that the uncertainty set touches the boundary of the feasible set at multiple points. We see that the next set of constraints only have an index that is 15.6% larger meaning they are likely to also limit the system flexibility to a certain extent. Subsequent sets have significantly increased flexibility indexes (up to 214% larger) and therefore have little effect on flexibility.

Figure 10 highlights the limiting components of Design 1. The limiting constraints corresponding to the active capacity constraints at each rank level are shaded according to the value of $F_{ellip+}$ (the lower the value, the most limiting the component is). We can see that the suppliers along with arcs attached to such suppliers limit flexibility the most. The rest of the components follow non-intuitive patterns, reflecting complex dependencies due to the network topology.

Finally, we consider a 141-node power distribution network which corresponds to an urban area in Caracas, Venezuela and was originally developed in [15]. The network data is extracted from MATPOWER, but again no arc capacities are provided (we assume $a^{C}=100$). Figure 11 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}$. We use $T_{ellip}(\delta)$ in combination with Problem (1.7) to rank all the inequality constraints. This problem was solved 270 times while iteratively removing the active constraints. We found that several solutions have the same value of $F_{ellip}$. Thus, the active constraints corresponding to solutions with equivalent $F_{ellip}$ values are combined and assigned to the same ranks. A total of 44 rank levels were identified, the ranking data corresponding to the first 30 levels is provided in Appendix A.

Figure 12 shows the limiting components corresponding to the 44 rank levels. We observe that the supplier is the most limiting component along with the arcs that form the spanning tree of the network (the path that connects all nodes in the network). The arcs become less relevant as they branch out towards the network boundaries. We thus see that the flexibility framework reveals topological limitations of the network. We also highlight that the flexibility index can be computed for this large system in 20 seconds. In contrast, given the large number of random parameters, computing the $SF$ index using Monte Carlo sampling is impractical.