Graph-Based Modeling and Decomposition of Energy Infrastructures

The IP method implemented in MadNLP.jl finds the solution of (2) by solving a sequence of barrier subproblems:

with a decreasing sequence for parameter $\mu$. The KKT conditions for (3) give the nonlinear equations:

where $A:=\nabla\boldsymbol{c}(\boldsymbol{x})$, $\boldsymbol{X}:=\mathop{\text{diag}}(\boldsymbol{x})$, and $\boldsymbol{Z}:=\mathop{\text{diag}}(\boldsymbol{z})$. A solution of KKT system (4) is obtained by computing primal-dual Newton steps $\boldsymbol{d}^{*}$ from:

where $W:=\nabla^{2}_{\boldsymbol{x}\boldsymbol{x}}\mathcal{L}(\boldsymbol{x},\boldsymbol{\lambda},\boldsymbol{z})$, $\Sigma:=\boldsymbol{X}^{-1}\boldsymbol{Z}$, and $\delta_{w},\delta_{c}>0$ are regularization parameters. The step $\boldsymbol{d}^{*}$ computed from (5) is safeguarded by a line-search filter procedure to induce global convergence (Wächter and Biegler, 2006). Typically, the solution of the linear system (5) is the most computationally intensive step in the IP method. This system is typically solved using direct linear solvers that are based on $LDL^{\top}$ factorizations (e.g., as those implemented in HSL routines (HSL, 2007)). Decomposition strategies based on Schur complements (Chiang et al., 2014) and iterative strategies (Curtis et al., 2012; Rodriguez et al., 2020) have also been proposed. The solution of (5) based on a direct block $LDL^{\top}$ factorization reveals the inertia (the number of positive, zero, negative eigenvalues) of $\boldsymbol{M}$. This inertia information is crucial in determining the acceptability of the computed step and in triggering the regularization of the linear system. However, inertia is not available when using iterative solution algorithms (as that proposed in this work). In MadNLP.jl, we use an inertia-free regularization strategy to determine the acceptability of the step (Chiang and Zavala, 2016). This method performs a simple negative curvature test to trigger regularization.

We propose a solution strategy for (5) based on an RAS scheme. We define some concepts and quantities that help explain our algorithm. Consider a partition $\{V_{k}\}_{k=1}^{K}$ of $V$; we call $V_{k}$ non-overlapping subdomains. This partition can be obtained by applying a graph partitioning scheme to $G$. We then construct a family of overlapping subdomains $\{V_{k}^{\omega}\}_{k=1}^{K}$ (these are constructed by expanding $V_{k}$). The expansion procedure is performed by progressively incorporating adjacent nodes (Shin et al., 2020b) (the size of overlap $\omega$ represents the expansion level). We observe that:

where $\sqcup$ denotes disjoint union. With $\{V_{k}\}_{k=1}^{K}$ and $\{V^{\omega}_{k}\}_{k=1}^{K}$, we define the corresponding index sets in the space of primal-dual variables in $\mathbb{R}^{\boldsymbol{n}+\boldsymbol{m}}$ as follows:

where $U_{i}\subseteq W:=\mathbb{Z}_{[1,\boldsymbol{n}+\boldsymbol{m}]}$ is the index set of $[x_{i};\lambda_{i}]$ in $[\boldsymbol{x};\boldsymbol{\lambda}]$. We now observe that:

We state the RAS scheme for solving (5) as:

Here, $\ell$ is the RAS iteration counter, $\boldsymbol{r}^{(\ell)}:=\boldsymbol{p}-\boldsymbol{M}\boldsymbol{d}^{(\ell)}$ is the residual; $\boldsymbol{M}_{k}:=\boldsymbol{M}[W^{\omega}_{k},W^{\omega}_{k}]$; $R_{k}:=\{e^{\top}_{i}\}^{\top}_{i\in W_{k}}$; $\widetilde{R}_{k}:=\{(\widetilde{e}^{k}_{i})^{\top}\}^{\top}_{i\in W_{k}}$, where $e_{i}$ is the $i$-th standard basis of $\mathbb{R}^{\boldsymbol{n}+\boldsymbol{m}}$, and $\widetilde{e}^{k}_{i}=e_{i}$ if $i\in W_{k}$ and $0$ otherwise.

The RAS scheme (6) involves the following steps. We first obtain the residual at the current step $\ell$; then, for each overlapping subdomains $\{V^{\omega}_{k}\}_{k=1}^{K}$, the associated residual is extracted as $R^{\top}_{k}\boldsymbol{r}^{(\ell)}$. The $k$-th subsystem is then solved by applying $\boldsymbol{M}^{-1}_{k}$. In MadNLP.jl, a factorization of $\boldsymbol{M}_{k}$ is computed with a direct solver and stored, so that the system can be repeatedly solved whenever the new right-hand-side is given. Subsequently, the solution $\boldsymbol{M}^{-1}_{k}R^{\top}_{k}\boldsymbol{r}^{(\ell)}$ for the $k$-th overlapping subdomain is restricted to the non-overlapping subdomain $V_{k}$ and then mapped back to the full-space by applying $\widetilde{R}^{\omega}_{k}$ (the indices associated with $V\setminus V_{k}$ are set to zero in this step). Key defining features of RAS are the concepts of overlap and restriction. The overlap allows the dampening of the adverse effect of the truncated domain, and the restriction procedure discards the part of the solution where the adverse truncation effect is strong.

It has been empirically observed that the convergence of the RAS algorithm improves as the size of overlap $\omega$ increases and as the conditioning of $\boldsymbol{M}$ improves (Cai and Saad, 1996). For positive definite $\boldsymbol{M}$, an exponential relationship between the convergence rate and the size of overlap has been established (Shin et al., 2020b). When $\omega=0$, the RAS scheme reduces to a block-Jacobi scheme (decentralized) while, when $\omega$ is maximal ($\boldsymbol{M}_{k}=\boldsymbol{M}$), the RAS scheme becomes a direct solution method (centralized). In this sense, RAS provides a bridge between fully decentralized and fully centralized schemes (thus providing flexibility). In the Schwarz submodule of ${\tt MadNLP.jl}$, $\omega$ is set automatically based on the relative size of $V_{k}$, and adaptively adjusted whenever a convergence issue occurs. Algorithm (6) uses a simple static iteration (also called a Richardson iteration), but more sophisticated iterative methods (e.g., the generalized mean residual (GMRES) method) can also be used (by treating $P$ as a preconditioner). Both Richardson and GMRES iterators are implemented in MadNLP.jl.

The construction of partitions for the RAS scheme is illustrated in Figure 2; these partitions correspond to the transient gas network example of Figure 1. Here, by partitioning the problem graph $G=(V,E)$ (first subfigure), the node set $V$ is divided into 4 subdomains $V_{1},V_{2},V_{3},V_{4}$ (second subfigure). These subdomains are associated with the primal-dual index sets $W_{1},W_{2},W_{3},W_{4}$ (third subfigure). After applying expansions, the associated blocks $\boldsymbol{M}_{1},\boldsymbol{M}_{2},\boldsymbol{M}_{3},\boldsymbol{M}_{4}$ of $\boldsymbol{M}$ are identified (last subfigure); these blocks are used by the RAS scheme (6).

The abstraction layers within MadNLP.jl are shown in Figure 3. The problem is modeled as an OptiGraph object (the core modeling object of Plasmo.jl). The OptiGraph is interfaced with multiple JuMP.jl models. The JuMP.jl model objects provide a set of local function oracles (objective, objective gradient, constraint, constraint Jacobian, and Lagrangian hessian). The Plasmo.jl-MadNLP.jl interface collects these local function oracles and creates a set of oracles for the full problem, where the local oracles are evaluated in parallel.

The Solver object of ${\tt MadNLP.jl}$ is created from the OptiGraph object of Plasmo.jl. The Solver object of MadNLP.jl uses a line-search filter IP method (Wächter and Biegler, 2006) to solve the problem. The step computation is performed by the linear solver specified by the user. The linear solver can be specified either as a direct solver or as the RAS solver. When the RAS solver is chosen, multiple subproblem solver objects are created by using standard direct solvers (e.g., by using Ma57 of HSL routines). These subproblem solvers are used for factorization and backsolve for $\boldsymbol{M}_{k}$ blocks. The RAS scheme (6) exploits multi-thread parallelism available in Julia. After termination of the IP solution procedure, the primal-dual solutions are sent back to the OptiGraph object and Model objects from JuMP.jl so that the user can query the solution via the interface provided by Plasmo.jl and JuMP.jl. See Figure 3 for a comparison with a conventional implementation.

We consider a transient gas network problem (Sundar and Zlotnik, 2018) of the form:

where $\dot{\rho}_{it}=\frac{{\rho}_{it}-{\rho}_{it-1}}{\Delta t}$, $P^{a}_{ijt}=W_{a}A_{ij}$, and $f^{a}_{ijt}=A_{ij}\varphi^{a}_{ijt}$. Here, $\mathcal{N}$ is the set of junctions; $\mathcal{P}$ is the set of pipelines; $\mathcal{C}$ is the set of compressors; $\mathcal{R}$ is the set of receipts; $\mathcal{D}$ is the set of demands; $\mathcal{R}(i)$ is the set of receipts at junction $i\in\mathcal{N}$; $\mathcal{D}(i)$ is the set of demands at junction $i\in\mathcal{N}$; $\mathcal{T}$ is the time index set; $\rho$ is the densitiy; $\varphi^{a}$ is the average mass flux; $\varphi^{-}$ is the negative mass flux; $\alpha$ is the compression ratio; $s$ is the supply; $d$ is the demand; $\dot{\rho}$ is the time derivative of density; $P^{a}$ is the power consumption of compressor; $f$ is the mass flow; $c$ is the gas price; $\gamma$ is the economic factor; $\lambda,\hat{L},L,D,A,\Delta t$, and $W_{a}$ are physical parameters. To implicitly enforce the periodicity, we let $\rho_{i0}=\rho_{iT}$, where $T$ is the end time index. The gas network under study consists of 2 compressors, 6 junctions (35 junctions after discretization), 4 pipelines (32 pipelines after discretization), 1 receiving points and 5 transfer points (which work either as receipt or delivery). The model is constructed using GasModels.jl (Bent et al., 2020).

We consider a multi-period AC power flow problem with storage (Geth et al., 2020) of the form:

Here, $\mathcal{G}$ is the set of generators; $\mathcal{N}$ is the set of buses; $\mathcal{E}$ is the set of (directed) branches; $\mathcal{E}^{R}$ is the set of branches with inverted directions; $\mathcal{S}$ is the set of storage; $\mathcal{T}$ is the time index set; $v\in\mathbb{C}$ is the voltage; $e\in\mathbb{R}$ is the state of charge; $s\in\mathbb{C}$ is the power flow; $s^{g}\in\mathbb{C}$ is the power generation; $s^{s}\in\mathbb{C}$ is the complex power injected by the storage; $sc\in\mathbb{R}$ is the charging rate; $sd\in\mathbb{R}$ is the discharging rate; $sqc\in\mathbb{R}$ is the reactive power slack; $c^{0},c^{1},c^{2}\in\mathbb{R}$ are the generation costs; $s^{d}\in\mathbb{C}$ is the power demand; $Y$ is the admittance; $T$ is the branch complex transformation parameter; $\eta$ is the charging efficiency; $s^{loss}$ is the storage energy loss; $\Delta t$ is the time interval. Note that (8) can be reformulated as an NLP with real variables by separately treating the real and imaginary part of the variables and equations (a polar formulation is used here). The power network under study is a variant of IEEE 14 bus test system; this comprises $14$ buses, $5$ generators, $1$ storage, $1$ shunt, and $20$ branches. The detailed model is constructed with PowerModels.jl (Coffrin et al., 2018).

We compare the proposed method (MadNLP.jl interfaced with Plasmo.jl and Schwarz), with the conventional method (MadNLP.jl interfaced with serial/parallel direct solvers Ma57 or MKL-Pardiso along with non-graph based algebraic modeling language JuMP.jl). The conventional methods are referred to as JuMP-Ma57 and JuMP-PardisoMKL, and the proposed method is referred to as Plasmo-Schwarz/Ma57. Furthermore, the mix of proposed/conventional approaches (JuMP-Schwarz/Ma57, Plasmo-Ma57, and Plasmo-PardisoMKL) is also tested together. For JuMP-Schwarz/Ma57, the graph partitioning tool METIS was used to partition the primal-dual coupling graph $\widetilde{G}$ directly. A Richardson scheme was used as an iterator for the RAS scheme. The study was performed by solving the gas (7) and power (8) problems while varying the size of the problems (by increasing the length of the prediction horizon). The code was run on a server computer equipped with 2 CPUs of Intel Xeon CPU E5-2698 v4 running 2.20GHz (20 core for each), and 20 threads are used for the computation. Code to reproduce the results can be found in https://github.com/zavalab/JuliaBox/tree/master/AdchemCaseStudy

For both problems, we found that the graph-based approach can significantly accelerate the solution (see Figure 4). In particular, comparing JuMP-Ma57 and Plasmo-Schwarz/Ma57, Plasmo-Schwarz/Ma57 becomes faster than JuMP-Ma57 when the prediction horizon is 3 days or more. Function evaluations are always faster in Plasmo.jl compared to JuMP.jl because the computational savings from function evaluations directly reduce the total solution time (parallelizing the function evaluation itself has no impact on the other part of the algorithm). On the other hand, one can see that the speed-up from parallel linear algebra is only observed when the problem size is sufficiently large (3 days in the gas network and 60 days in the power network). This is because the reduction in the problem size also reduces the overlap size. In our implementation, we set the size of overlap using the relative size of the block (the size of the overlap is reduced if the overall problem size is reduced). As a result, the RAS scheme (6) may become slow, and the number of required factorization/backsolve steps increases. This indicates that the use of RAS is beneficial only when the problem size is sufficiently large. For the gas problems, the acceleration of linear algebra computations was more pronounced. In contrast, for the power problems, the acceleration of function evaluations was more pronounced. This is because the AC power flow formulation has a large number of nonlinear expressions. By comparing the linear solver time for JuMP-Schwarz/Ma57 and Plasmo-Schwarz/Ma57, we see the advantage of using a graph-based modeling language for obtaining the partitions $\{W_{k}\}_{k=1}^{K}$. We recall that for JuMP-Schwarz/Ma57, the Metis graph partitioning routine is directly applied to $\widetilde{G}$ while Plasmo-Schwarz/Ma57 uses the user-provided problem graph $G$. One can observe that, in general, the linear solver time is shorter for PlasmoNLP-Schwarz/Ma57. This indicates that the user-provided graph information can be leveraged for obtaining high-quality partitions.