# Gradient Methods for Scalable Multi-value Electricity Network Expansion
Planning

Anthony Degleris  Abbas El Gamal  Ram Rajagopal

###### Abstract

We consider multi-value expansion planning (MEP), a general bilevel
optimization model in which a planner optimizes arbitrary functions of the
dispatch outcome in the presence of a partially controllable, competitive
electricity market. The MEP problem can be used to jointly plan various grid
assets, such as transmission, generation, and battery storage capacities;
examples include identifying grid investments that minimize emissions in the
absence of a carbon tax, maximizing the profit of a portfolio of renewable
investments and long-term energy contracts, or reducing price inequities
between different grid stakeholders. The MEP problem, however, is in general
nonconvex, making it difficult to solve exactly for large real-world systems.
Therefore, we propose a fast stochastic implicit gradient-based heuristic
method that scales well to large networks with many scenarios. We use a strong
duality reformulation and the McCormick envelope to provide a lower bound on
the performance of our algorithm via convex relaxation. We test the
performance of our method on a large model of the U.S. Western Interconnect
and demonstrate that it scales linearly with network size and number of
scenarios and can be efficiently parallelized on large machines. We find that
for medium-sized 16 hour cases, gradient descent on average finds a 5.3x lower
objective value in 16.5x less time compared to a traditional reformulation-
based approach solved with an interior point method. We conclude with a large
example in which we jointly plan transmission, generation, and storage for a
768 hour case on 100 node system, showing that emissions penalization leads to
additional 40.0% reduction in carbon intensity at an additional cost of
$17.1/MWh.

##  1 Introduction

Expanding the electric power system is essential to future grid
electrification and greenhouse gas curtailment [1, 2, 3, 4]. However, there
are numerous barriers inherent to developing new grid assets, including high
upfront capital costs, long and uncertain development timelines, and extensive
regulatory restrictions [5]. There are also significant challenges associated
with value assessment and alignment. First, neglecting benefits outside of
economic cost and reliability may under-estimate the value of a new
transmission, generation, or storage project. Second, various utilities,
private organizations, and government agencies generally participate in
planning grid expansions, all with competing values and objectives.

To address these challenges, modern expansion planning studies attempt to
assess a range of benefits outside of cost and reliability, such as greenhouse
gas emissions reductions or increased competition among generating assets. For
example, the California Independent System Operator (ISO) has developed the
Transmission Economic Assessment Methodology (TEAM) to holistically evaluate
transmission projects by assessing economic benefits separately for different
market participants, accounting for downstream impacts on market prices and
generator competition, and measuring environmental impacts [6].

Methodologies such as TEAM are effective tools for evaluating new grid
investments in a holistic, market-aware framework. In this paper, we focus on
optimization-based expansion planning models, a complementary approach for
identifying potential grid investments from a large set of candidate
solutions. Specifically, we study multi-value expansion planning (MEP), a
general optimization-based framework for expansion planning problems that
allows the planner to optimize arbitrary functions of the dispatch outcome in
an (partially controllable) competitive electricity market. MEP supports
jointly planning transmission, generation, and storage assets in either a
centralized environment, e.g., a vertically integrated system, or as part of a
decentralized market with multiple agents optimizing distinct objectives.

We formulate multi-value expansion planning as a bilevel optimization problem,
in which the upper-level problem is the planner’s investment decision and the
lower-level problem is market clearing process. We give several examples of
MEP problems, including minimizing carbon emissions in lieu of imposing a
carbon tax or maximizing the profit of a portfolio of renewable investments
and long-term energy contracts. Using strong duality [7, 8], we reformulate
the MEP problem as a single-level optimization problem with bilinear equality
constraints. When the dispatch costs are convex quadratic or linear, the
resulting problem is a nonconvex quadratic program, which is NP-hard in
general. Therefore, inspired by its success in machine learning [9], we
propose applying gradient descent directly to the implicit form of multi-value
expansion planning and show that it converges to local stationary point in
polynomial time. We also use a simple convex relaxation of the MEP problem to
bound the performance of our algorithm and generate good initial candidate
solutions. We show that our approach scales well to large problems with many
scenarios and outperforms applying interior point methods to the reformulated
problem. Finally, we solve MEP problems on a large model of the U.S. Western
Interconnect, demonstrating that gradient descent can find good solutions to
realistically-sized grid planning problems.

##  2 Literature review

In this work, we focus on expansion planning for bulk power systems and
consider resources such as generation, transmission, and storage. Several
comprehensive reviews are available on electricity grid expansion planning,
usually specific to an individual resource; see, for example, [10, 11, 12].
Expansion planning studies classically model the system as vertically
integrated and minimize total system cost. Recently, however, numerous studies
directly model market behavior using hierarchical equilibrium models; see [13]
and references therein. Here, we adopt the latter market-based approach. We
also limit the scope of our study to single investment period expansion
planning with scenario-based uncertainty, although much of our work is
straightforward to extend to more general settings.

#### Bilevel (market-based) expansion planning.

Deregulated electricity markets make expansion planning significantly more
complex because different agents have distinct, competing interests. In this
setting, the resulting expansion planning problem is formulated as bilevel
optimization problem or a Stackelberg game [14, 13], which can be cast as a
mathematical program with equilibrium constraints (MPEC). Most bilevel
expansion planning problems fall into one of two classes. The first are
proactive expansion planning problems which maximize social welfare while
anticipating strategic (profit-maximizing) generator investment and bidding
behavior [14, 15, 16, 17]. Under perfect competition, however, no generation
entity has sufficient market power to affect prices and the Stackelberg
equilibrium can be computed as the result of a (standard) single-level
planning problem [7, 18]. The second class of problems is profit maximization
problems which consider a single or coalition of market participants choosing
expansion decisions to maximize individual profit, as opposed to overall
welfare [19, 20, 21, 22]. Although the same reduction to single-level planning
under perfect competition still applies, it is generally less relevant since
considerations about competition are essential to pricing for an individual
firm. Finally, a recent work [23] considers a system operator expanding
transmission and storage to reduce greenhouse gas emissions in a market
environment without a carbon tax.

One contribution of our work is to provide a unified framework for solving a
large class of bilevel expansion planning problems. In particular, the
implementation and performance of our gradient algorithms generally do not
change significantly across a broad range of bilevel expansion planning
problems, regardless of their specific problem details.

#### Solution algorithms.

Generally, the standard approach for solving a bilevel expansion planning
problem is to reduce it to a nonconvex single-level problem, then apply
standard techniques for solving nonconvex programs. There are two common
approaches for this: to write the equilibrium conditions for the lower level
problem explicitly using the Karush-Kuhn-Tucker (KKT) conditions [24, 25], or
to use strong duality to reformulate the optimality conditions as an
inequality constraint [7, 26, 27]. Both approaches introduce several bilinear
terms and make the resulting problem nonconvex; see [8] for a comparison
between the two. In this work, we use the strong duality-based approach to
derive our convex relaxation.

Our main solution algorithm, inspired by similar work on bilevel optimization
in machine learning [28, 29, 30], is to apply gradient descent directly to the
original implicit form of the bilevel problem. This requires differentiating
through an optimization problem, which can be performed efficiently for a
large class of problems using implicit differentiation [31, 29]. The main
advantages of our approach are that it scales well to high-dimensional and
stochastic problems with many scenarios [32], can be easily warm started if
the problem is perturbed, and is simple to implement.

##  3 Multi-value expansion planning

In this section, we first introduce the dispatch model, which is used to make
real-time operational decisions given long-term planning decisions, and
describe classical expansion planning. We then formulate multi-value expansion
planning and show how classical planning can be viewed as a special case.
Finally, we show how to extend multi-value expansion planning to uncertain
environments.

###  3.1 Dispatch model

Consider an electricity network characterized by long-term physical parameter
vector $\eta\in\mathcal{H}\subset\mathbf{R}^{K}$, which we refer to as the
network parameters or network investments, which may include the values of any
long-term investment decisions. The network parameters include, for example,
transmission line capacities and susceptances, generator capacities, and
battery capacities. We refer to the set $\mathcal{H}\subset\mathbf{R}^{M}$ as
the set of allowable network investments, which we assume to be compact. The
dispatch model is the problem of choosing real-time operational decisions to
minimize the system cost,

$\begin{array}[]{ll}\textrm{minimize}&c(x)\\\\[2.5pt] \textrm{subject
to}&A(\eta)x\leq b(\eta),\end{array}$ (1)

where the variable is $x\in\mathbf{R}^{N}$, which specifies all short-term
operation decisions, such as the power output of each generator. The cost
function $c:\mathbf{R}^{N}\rightarrow\mathbf{R}$ specifies the cost of
operating the system, e.g., the total cost of producing power. The functions
$A:\mathbf{R}^{K}\rightarrow\mathbf{R}^{M\times N}$ and
$b:\mathbf{R}^{K}\rightarrow\mathbf{R}^{M}$ describe the $M$ affine inequality
constraints of the system as a function of the network parameters. These
constraints may include, for example, limits on generation and transmission
capacity. For ease of notation, we only consider problems with inequality
constraints; however, all our technical results should be straightforward to
extend to dispatch models with equality constraints as well. Throughout the
text, we may refer to the operator, the entity that solves the dispatch model
and sets $x$. We will make the following assumptions about the dispatch model.

######  Assumption 1 (Strict feasibility, i.e., Slater’s condition).

For all $\eta\in\mathcal{H}$, there exists $x\in\mathbf{R}^{N}$ such that
$A(\eta)x<b(\eta)$.

######  Assumption 2 (Differentiable strongly convex costs).

The cost function $c$ is strongly convex and twice differentiable with
Lipschitz continuous second derivatives.

#### Dispatch map.

The first two assumptions imply that for all $\eta\in\mathbf{R}^{K}$, problem
(1) is convex and has a unique solution. Let
$x^{*}(\eta):\mathbf{R}^{K}\rightarrow\mathbf{R}^{N}$ denote the optimal
choice of $x$ in (1) given network parameters $\eta$. The function $x^{*}$ is
well-defined, since (1) has a unique solution. We call the function $x^{*}$
the dispatch map, since it maps the network parameters to the solution of the
dispatch model.

#### Dual problem.

The dual problem of (1) is

$\begin{array}[]{ll}\textrm{maximize}&u(\lambda,\eta)\\\\[2.5pt]
\textrm{subject to}&\lambda\geq 0,\end{array}$ (2)

where the variable is $\lambda\in\mathbf{R}^{M}$, and
$u(\lambda,\eta)=-\lambda^{T}b(\eta)-c^{\star}(-A(\eta)^{T}\lambda)$. Here the
function $c^{\star}(y)=-\inf_{x}\left(c(x)-y^{T}x\right)$ is the convex
conjugate of $c$, which is concave and, in many cases, has a simple analytical
expression [33]. The variable $\lambda$ is the dual variable of the constraint
$A(\eta)x\leq b(\eta)$ in (1). We make two additional assumptions that
guarantee dual uniqueness.

######  Assumption 3 (Non-redundant constraint).

Denote $\tilde{a}_{1}(\eta),\ldots,\tilde{a}_{M}(\eta)\in\mathbf{R}^{n}$ the
rows of $A(\eta)$ and let $m_{1}(\eta),\ldots,m_{p}(\eta)\in\\{1,\ldots,M\\}$
be the indices of binding constraints at optimality, i.e., the indices $m$
such that $\tilde{a}_{m}(\eta)^{T}x^{*}(\eta)=b(\eta)$. Then, for all
$\eta\in\mathcal{H}$, the vectors
$\tilde{a}_{m_{1}(\eta)}(\eta),\ldots,\tilde{a}_{m_{p}(\eta)}(\eta)$ are
linearly independent.

######  Assumption 4 (Strict complementarity).

Denote $\mathcal{B}(\eta)=\\{m_{1}(\eta),\ldots,m_{p}(\eta)\\}$ the indices of
the binding constraints at optimality, as before and let $\lambda^{*}$ be a
solution to (2). Then, for all $\eta\in\mathcal{H}$, $\lambda_{m}^{*}>0$ for
$m\in\mathcal{B}(\eta)$.

Assumptions 2, 3, and 4 together imply that the solution to (2) exists
uniquely for all $\eta\in\mathbf{R}^{K}$; see Appendix A.2. We call (2) the
dual dispatch problem and define
$\lambda^{*}(\eta):\mathbf{R}^{K}\rightarrow\mathbf{R}^{N}$, which we call the
dual dispatch map, the optimal choice of $\lambda$ in (2) as a function of
$\eta$. We make two final assumptions to simplify several technical results.

######  Assumption 5 (Bounded domain).

There exist convex, compact sets $\mathcal{X}\subset\mathbf{R}^{N}$ and
$\Lambda\subset\mathbf{R}^{M}$ such that $x^{*}(\eta)\in\mathcal{X}$ and
$\lambda^{*}(\eta)\in\Lambda$ for all $\eta\in\mathcal{H}$.

######  Assumption 6 (Affine coefficients).

The functions $A$ and $b$ are affine (and hence differentiable) in $\eta$.

#### Discussion of assumptions.

Assumption 1 says that for every allowable network investment, the dispatch
model is strictly feasible. If we assume load curtailment is a feasible,
albeit expensive, dispatch outcome, then we can generally guarantee
feasibility by allowing the ‘zero power’ solution, i.e., no power is produced
or consumed by any device in the network. Given feasibility, we can also
guarantee strict feasibility by perturbing the right hand side of the
inequality constraints by a small constant. Assumptions 3, 4, and 5 can also
generally be guaranteed through careful problem construction and small
perturbations. Perhaps most restrictive is the strong convexity requirement in
Assumption 2. However, this too can be guaranteed through problem
perturbation, e.g., adding a small quadratic penalty to each variable.
Moreover, empirical results suggest that this assumption can be relaxed, so
long as the resulting problem has a unique solution for every
$\eta\in\mathcal{H}$. Finally, Assumption 6 is mainly intended to simplify
various technical results and can be relaxed to just twice differentiability
if needed.

###  3.2 Expansion planning

Expansion planning is the problem of choosing both long-term investments and
real-time operation decisions to minimize overall system cost,

$\begin{array}[]{ll}\textrm{minimize}&\gamma^{T}\eta+c(x)\\\\[2.5pt]
\textrm{subject to}&A(\eta)x\leq b(\eta)\\\\[2.5pt]
&\eta\in\mathcal{H},\end{array}$ (3)

where the variables are $x\in\mathbf{R}^{N}$ and $\eta\in\mathbf{R}^{K}$. The
vector $\gamma\in\mathbf{R}^{K}$ gives the cost of the network investment. Two
common types of constraints on network investments include continuous box
constraints, where
$\mathcal{H}^{\textrm{cont}}=\\{\eta:\eta^{\min}\leq\eta\leq\eta^{\max}\\}$
for $\eta^{\min},\eta^{\max}\in\mathbf{R}^{K}$, and binary box constraints,
where
$\mathcal{H}^{\textrm{int}}=\\{\eta^{\min}_{1},\eta^{\max}_{1}\\}\times\cdots\times\\{\eta^{\min}_{k},\eta^{\max}_{k}\\}$.
Throughout the remainder of the paper, we will consider continuous investments
decisions, $\mathcal{H}=\mathcal{H}^{\textrm{cont}}$, but we will refer to the
discrete case when relevant.

#### Bilevel form.

The classical problem (3) can be rewritten as

$\begin{array}[]{ll}\textrm{minimize}&\gamma^{T}\eta+c(x^{*}(\eta))\\\\[2.5pt]
\textrm{subject to}&\eta\in\mathcal{H},\end{array}$ (4)

where the variable is now just $\eta\in\mathbf{R}^{K}$. We refer to (4) as the
bilevel form of the classical expansion problem (and refer to (3) as the joint
form), since it involves optimizing a function of the outcome of another
optimization problem [17]. The bilevel form of classical expansion planning is
useful for deriving decomposition-based algorithms that consist of generating
cuts on $c^{\textrm{opt}}(\eta)=c(x^{*}(\eta))$ [34].

###  3.3 Multi-value expansion planning

Multi-value expansion planning is the problem of choosing long-term
investments to minimize some arbitrary function of the dispatch outcome,

$\begin{array}[]{ll}\textrm{minimize}&\gamma^{T}\eta+h(x^{*}(\eta),\lambda^{*}(\eta))\\\\[2.5pt]
\textrm{subject to}&\eta\in\mathcal{H},\end{array}$ (P1)

where the variable is $\eta\in\mathbf{R}^{K}$. The function
$h:\mathbf{R}^{N}\times\mathbf{R}^{M}\rightarrow\mathbf{R}$ is the planner’s
objective and may differ from the dispatch cost $c$. Multi-value expansion
planning generalizes classical expansion planning, since we can set
$h(x,\lambda)=c(x)$. When $h(x,\lambda)\neq c(x)$, the planner chooses $\eta$
to minimize some function of $x^{*}(\eta)$ and $\lambda^{*}(\eta)$, given that
the dispatch outcomes $x^{*}(\eta)$ and $\lambda^{*}(\eta)$ are still chosen
to minimize cost. We refer to (P1) as the implicit form of multi-value
expansion planning because the objective includes an implicit function. We
also refer to the entity that solves multi-value expansion planning as the
planner, which, in general, need not be the same entity as the operator.
Multi-value expansion planning can be considered as a general class of bilevel
expansion planning problems.

The name multi-value comes from the fact that the planner may consider
multiple values related to the dispatch outcome, not just cost. In a modern
wholesale electricity market, the operator chooses the dispatch outcome
$x^{*}(\eta)$ and prices $\lambda^{*}(\eta)$ to minimize total cost, since it
coincides with the welfare maximizing and market clearing solution [35].
However, the planner may not be interested in minimizing total cost; if the
planner, for instance, is a utility purchasing power on behalf of their
consumers, then they are interested in minimizing the price they pay for
power, i.e., the surplus they capture as a consumer, not the overall cost.
Similarly, other agents investing in network upgrades will have other
differing and possibly competing objectives, e.g., a state policy maker
funding grid investment might wish to minimize system emissions. These agents
cannot change the cost-minimizing nature of the dispatch model, but they can
strategically choose network investments that account for it. In this way, the
implicit form of multi-value expansion planning can be interpreted as a direct
formulation of the underlying Stackelberg game: the planner is the leader and
chooses $\eta$ in the first round, and the operator is the follower and
chooses $x^{*}$ and $\lambda^{*}$ in the second round following a
deterministic strategy known to both agents.

###  3.4 Stochastic multi-value expansion planning

In many settings, there is uncertainty about the state of the system, e.g.,
future load, technology costs, and renewable production, after the investment
decisions have been made. We can explicitly model this uncertainty by solving
a stochastic version of the multi-value expansion planning problem in which
the dispatch model is selected uniformly at random from one of $S$ scenarios.
The dispatch model in scenario $s$ is

$\begin{array}[]{ll}\textrm{minimize}&c_{s}(x)\\\\[2.5pt] \textrm{subject
to}&A_{s}(\eta)x\leq b_{s}(\eta),\end{array}$ (5)

where the variable is $x\in\mathbf{R}^{N}$, and
$c_{s}:\mathbf{R}^{N}\rightarrow\mathbf{R}$,
$A_{s}:\mathbf{R}^{K}\rightarrow\mathbf{R}^{M\times N}$, and
$b:\mathbf{R}^{K}\rightarrow\mathbf{R}^{M}$ define the costs and constraints
of the scenario. Let $x_{s}^{*}(\eta):\mathbf{R}^{K}\rightarrow\mathbf{R}^{N}$
and $\lambda_{s}^{*}(\eta):\mathbf{R}^{K}\rightarrow\mathbf{R}^{M}$ denote the
solution to the $s$th dispatch problem and its dual, respectively. Then
stochastic multi-value expansion planning is the problem of choosing long-term
investments to minimize the expected value of some function of the dispatch
outcome,

$\begin{array}[]{ll}\textrm{minimize}&\gamma^{T}\eta+(1/S)\sum_{s=1}^{S}h(x_{s}^{*}(\eta),\lambda_{s}^{*}(\eta))\\\\[2.5pt]
\textrm{subject to}&\eta\in\mathcal{H},\end{array}$ (6)

where the variable is $\eta\in\mathbf{R}^{K}$.

##  4 Illustrative examples

In this section we give three examples that illustrate multi-value expansion
planning problems on a small 6-bus test network. We first introduce a simple
dispatch model in which generators with linear costs dispatch power across a
transmission network using linearized power flow constraints. Then we consider
optimal transmission investment for three different objectives: cost
minimization across the system, emissions minimization, and profit
maximization for an individual generator. Finally, we give an example of a
dynamic dispatch model with multiple time periods, for which the same planning
problems can be solved.

#### Dispatch model.

Consider an electricity network with $n$ nodes.111Although $n$ is used as an
index for the primal dispatch variables in the previous section, we use it
here as the number of nodes in the network in order to remain consistent with
standard power systems notation. Other symbols, namely $c$ and $K$, are
slightly overloaded but in a fashion consistent with Section 3.  Each node $i$
has electricity demand $p_{i}^{\textrm{desired}}\in\mathbf{R}_{+}$, which can
be curtailed (if needed) at a large cost $\delta\in\mathbf{R}_{+}$. Node $i$
may also have a generator with capacity $g_{i}^{\max}\in\mathbf{R}^{+}$ and
linear cost $c_{i}\in\mathbf{R}_{+}$; nodes without generators have
$g_{i}^{\max}=0$. The nodes are connected by a set of $K$ transmission lines.
Each line $k$ has electrical susceptance $\beta_{k}\in\mathbf{R}_{+}$ and per-
susceptance capacity $f_{k}^{\max}\in\mathbf{R}^{+}$. Each line also has a
source node and a sink node $\ell_{1}(k)\neq\ell_{2}(k)\in\\{1,\ldots,n\\}$,
respectively. The topology of the whole network can be specified by the
incidence matrix $L\in\mathbf{R}^{n\times K}$ with entries

$L_{ik}=\begin{cases}1&i=\ell_{1}(k)\\\ -1&i=\ell_{2}(k)\\\
0&\textrm{otherwise}\end{cases}.$

The transmission line flows $f_{k}\in\mathbf{R}$ must satisfy both capacity
constraints, $-\eta_{k}f_{k}^{\max}\leq f_{k}\leq\eta_{k}f_{k}^{\max}$, and
voltage phase angle constraints, $f=\mathbf{diag}(\eta)L^{T}\theta$, where
$\theta\in\mathbf{R}^{n}$ are the nodal voltage phase angles. The phase angles
are relative to an arbitrary reference node, so we set $\theta_{1}=0$. The
full dispatch model is

$\begin{array}[]{ll}\textrm{minimize}&\delta\mathbf{1}^{T}p+c^{T}g+\epsilon\left(\|p\|_{2}^{2}+\|g\|_{2}+\|f\|_{2}^{2}+\|\theta\|_{2}^{2}\right)\\\\[2.5pt]
\textrm{subject to}&f=\mathbf{diag}(\beta)L^{T}\theta\\\\[2.5pt]
&g-p=Mf\\\\[2.5pt] &\theta_{1}=0\\\\[2.5pt] &0\leq p\leq
p^{\textrm{desired}}\\\\[2.5pt] &0\leq g\leq g^{\max}\\\\[2.5pt]
&-\mathbf{diag}(\beta)f_{j}^{\max}\leq
f\leq\mathbf{diag}(\beta)f^{\max}\end{array}$ (7)

where the variables are $p,g,\theta\in\mathbf{R}^{n}$ and
$f\in\mathbf{R}^{K}$. Above, $\epsilon>0$ is a small constant coefficient of a
quadratic term used to ensure the objective is strongly convex. We consider
transmission expansion planning problems in which the susceptances can be
modified to meet some objective. Each line susceptance has a box constraint
$\beta_{k}^{\min}\leq\beta_{k}\leq\beta_{k}^{\max}$ and investment cost
$\gamma_{k}\in\mathbf{R}_{+}$. Here, the operation variable is
$x=(p,g,\theta,f)$ and the planning variable is $\eta=\beta$. The dispatch
cost function is
$c(x)=\delta\mathbf{1}^{T}p+c^{T}g+\epsilon\left(\|p\|_{2}^{2}+\|g\|_{2}+\|f\|_{2}^{2}+\|\theta\|_{2}^{2}\right)$,
and the constraint functions $A(\eta)$ and $b(\eta)$ include all the linear
and constant coefficients in the constraints, respectively; each coefficient
is an affine function of $\eta$.

To make things concrete, we consider a classical 6-node test network from [36,
37] modified to include generator fuel costs and emissions rates; see Figure
1. The network has 6 nodes and 3 generators; initially, there 6 transmission
lines and one of the generators is disconnected, so not all demand can be
satisfied; see Panel A. The planner may upgrade any of the existing 6
transmission lines or build new transmission lines between any two nodes in
order to optimize their objectives. However, the generator power outputs and
transmission line flows themselves are operational decisions chosen by a
system operator to minimize cost and clear the market in (7).

We consider the following three examples of planning problems on this network:
cost minimization, emissions minimization, and individual-generator profit
maximization.

Figure 1:  Illustrative examples of multi-value expansion planning. (A)
Initial network before transmission expansion. (B) Transmission expansion
(blue lines) that minimizes total cost. (C) Expansion that minimizes emissions
(i.e., emissions weight $w=\infty$). (D) Expansion that maximizes the price at
node 6.

#### Cost minimization.

In the cost minimization problem, the planner solves

$\begin{array}[]{ll}\textrm{minimize}&\gamma^{T}\eta+c(x^{*}(\eta))\\\\[2.5pt]
\textrm{subject to}&\eta\in\mathcal{H},\end{array}$ (8)

where the variable is $\eta\in\mathbf{R}^{K}$. This problem corresponds to
standard expansion planning and can be reformulated as a joint optimization
problem over $p$, $g$, $\theta$, and $\eta$. However, the resulting problem is
nonconvex because the constraint $f=\mathbf{diag}(\beta)L^{T}\theta$ is
bilinear in $\eta$ and $\theta$. The transmission expansion decisions for the
6-bus system are displayed as blue lines in Figure 1, Panel B.

#### Emissions minimization.

Suppose each generator has nonnegative emissions rate
$e_{i}\in\mathbf{R}_{+}$. A planner may wish to invest in transmission in
order to alleviate congestion that limits the use of low emissions generators.
In the emissions-aware expansion planning problem [23], the planner solves

$\begin{array}[]{ll}\textrm{minimize}&\gamma^{T}\eta+c(x^{*}(\eta))+we^{T}g^{*}(\eta)\\\\[2.5pt]
\textrm{subject to}&\eta\in\mathcal{H},\end{array}$ (9)

where the variable is $\eta\in\mathbf{R}^{K}$. (Here, $g^{*}(\eta)$ represents
the optimal choice of $g$ in (7) given $\eta$.) The weight $w>0$ can be
interpreted as the implied cost of emissions. Although emissions are not
explicitly incorporated into the generator costs $c$ (e.g., via a carbon tax),
the planner has chosen a dollar value to associate with mitigated emissions.
Emissions-aware planning is especially useful for policy makers seeking to
reduce emissions in the power grid, and provides a unified approach for
jointly planning economic and policy-driven expansions.

We display results for an emissions-aware expansion planning example in Figure
1, Panel C. Notably, the optimal transmission investment is different from
that of cost minimization in Panel B: instead of connecting to the loads at
nodes 2 and 4, the islanded generator is connected to a high emissions
generator at node 3. To better understand the trade-off between cost and
emissions, we plot their approximate Pareto frontier in Figure 2. The results
suggest strategically selecting transmission investments can lead to a large
decrease in emissions at a slight (10-30%) increase in total cost.

Figure 2:  Trade-off between total cost (investment cost plus dispatch cost)
and emissions for Problem (9) using the 6-node network in Figure 1. (Black
dots) Solutions generated by varying the emissions weight $w\in[0,10]$ and
solving the resulting multi-value expansion planning problem using Algorithm
1. (Green line) Pareto frontier for the cost-emissions trade-off. (Red dots)
The cost-emissions profile of the networks in Figure 1, Panels B and C.

#### Profit maximization.

If a merchant generator is responsible for funding transmission investment,
they are naturally more interested in their own profit than the overall cost
of the system. Let $\nu^{*}(\eta)\in\mathbf{R}^{n}$ be the optimal dual
variable of the constraint $g-d=Mf$ in (7). The variable $\nu_{i}^{*}(\eta)$
gives the locational marginal price at node $i$; after the operator solves the
dispatch problem, they publish the prices $\nu^{*}(\eta)$ to correctly
incentivize each generator to produce the correct amount of power. The profit
maximization problem is

$\begin{array}[]{ll}\textrm{minimize}&\gamma^{T}\eta+\sum_{i\in\mathcal{G}}\nu_{i}^{*}(\eta)\cdot
g_{i}^{*}(\eta)\\\\[2.5pt] \textrm{subject to}&\eta\in\mathcal{H},\end{array}$
(10)

where the variable is $\eta\in\mathbf{R}^{K}$, and
$\mathcal{G}\subset\\{1,\ldots,n\\}$ is the set of generators owned by the
merchant.222 The outer objective $h$ in the profit maximization problem is
bilinear and hence nonconvex. As we will see later, this means the relaxation
derived in Section 6.2 will not be convex either. Fortunately, because the
objective is bilinear, the same techniques can be applied to relax the terms
in the objective as well.  In the 6-bus example, we consider a profit
maximization problem for a merchant that owns generator $\mathcal{G}=\\{6\\}$
and is paying for a transmission upgrade to connect to the system. We show the
resulting expansion decisions in Figure 1, Panel D. Note that profit
maximization leads to a different set of network upgrades than either cost
(Panel B) or emissions minimization (Panel C); the additional transmission
allows for increased export capacity from node six and therefore higher local
prices.

###  4.1 Dynamic model

We can extend the dispatch model in (7) to be dynamic and model multiple time
periods. The model is exactly the same as (7), except that it is solved over
$T$ time periods with battery storage devices. It is given by the optimization
problem,

$\begin{array}[]{ll}\textrm{minimize}&\sum_{t=1}^{T}\left(\delta\mathbf{1}^{T}p_{t}+c^{T}g_{t}+\epsilon\left(\|p_{t}\|_{2}^{2}+\|g_{t}\|_{2}+\|f_{t}\|_{2}^{2}+\|\theta_{t}\|_{2}^{2}+\|s_{t}\|_{2}^{2}+\|v_{t}\|_{2}^{2}\right)\right)\\\\[2.5pt]
\textrm{subject to}&f_{t}=\mathbf{diag}(\beta)L^{T}\theta_{t}\\\\[2.5pt]
&g_{t}-p_{t}-v_{t}=Mf_{t}\\\\[2.5pt] &s_{t}=s_{t-1}+v_{t}\\\\[2.5pt]
&\theta_{1}=0\\\\[2.5pt] &0\leq p_{t}\leq p_{t}^{\textrm{desired}}\\\\[2.5pt]
&0\leq g_{t}\leq g_{t}^{\max}\\\\[2.5pt] &0\leq
s_{t}\leq\mathbf{diag}(\omega)v^{\max}\\\\[2.5pt] &-v^{\max}\leq v\leq
v^{\max}\\\\[2.5pt] &-\mathbf{diag}(\beta)f_{j}^{\max}\leq
f_{t}\leq\mathbf{diag}(\beta)f^{\max}\end{array}$ (11)

where the variables are $p_{t},g_{t},\theta_{t},s_{t},v_{t}\in\mathbf{R}^{n}$
and $f_{t}\in\mathbf{R}^{K}$ for $t=1,\ldots,T$. Note that the desired load at
node $i$ is now the time varying quantity
$(p_{t}^{\textrm{desired}})_{i}\in\mathbf{R}_{+}$. Similarly, the maximum
nodal generation capacity $(g_{t}^{\max})_{i}$ also varies with time. Costs
and transmission capacities are considered fixed across time, but could be
made time varying if desired. The model also includes lossless battery storage
at each node with maximum (dis)charge power $v^{\max}_{i}\in\mathbf{R}_{+}$
and duration $\omega_{i}\in\mathbf{R}_{+}$. The battery dynamics are given by
the constraint $s_{t}=s_{t-1}+v_{t}$, with $s_{0}=0$ as an initial condition.
In this problem, we set the planning variable to be
$\eta=(\beta,g^{\max},v^{\max})$, jointly planning all three device
capacities.

The main difference between the static and dynamic models is the addition of
storage (with expandable capacities). Planning problems for the dynamic model
are almost identical to their static counterparts. For example, the emissions
minimization problem is

$\begin{array}[]{ll}\textrm{minimize}&\gamma^{T}\eta+c(x^{*}(\eta))+w\sum_{t}e^{T}g_{t}^{*}(\eta)\\\\[2.5pt]
\textrm{subject to}&\eta\in\mathcal{H},\end{array}$ (12)

where the variable is $\eta\in\mathcal{H}$. The dynamic model is, of course,
straightforward to extend in the usual ways. For example, we could add other
devices like DC transmission lines or flexible loads. We could also add
additional constraints, such as ramping constraints or lossy batteries. One
benefit of our general planning framework is that it is simple to extend all
our results to these cases by changing $c$, $A$, $b$, and the resulting
Jacobian $\partial z^{*}(\eta)$.

##  5 Scalable gradient algorithm

In this section, we propose a local gradient algorithm that can be applied
directly to the implicit form of multi-value expansion planning in (P1). We
first explain how to differentiate the solution maps $x^{*}(\eta)$ and
$\lambda^{*}(\eta)$ using the implicit function theorem. Then, we derive a
simple gradient-based algorithm that converges to a local optimum in
polynomial time. We generalize this algorithm to the stochastic setting and
give two heuristic improvements that empirically seem to improve the
convergence rate and the solution quality.

###  5.1 Implicit differentiation

The vectors $x\in\mathbf{R}^{N}$ and $\lambda\in\mathbf{R}^{M}$ are solutions
to the primal and dual dispatch models (1) and (2), respectively, if and only
if they satisfy $\kappa(x,\lambda,\eta)=0$, where

$\kappa(x,\lambda,\eta)=\begin{bmatrix}\nabla c(x)+A(\eta)^{T}\lambda\\\
\mathbf{diag}(\lambda)(A(\eta)x-b(\eta))\end{bmatrix}.$ (13)

The operator $\kappa:\mathbf{R}^{N+M+K}\rightarrow\mathbf{R}^{N+M}$ is the
Karush-Kuhn-Tucker (KKT) operator; finding a solution to a convex optimization
problem and its dual is equivalent to solving the system of equations defined
by the problem’s KKT operator. In this way, we can think of the primal and
dual dispatch maps $x^{*}(\eta)$ and $\lambda^{*}(\eta)$ as the solutions to
the nonlinear system of equations $\kappa(x,\lambda,\eta)=0$ parameterized by
$\eta$. This gives a method for differentiating the dispatch maps.

######  Theorem 1 (Implicit Function Theorem, [31], Chapter 1B).

Consider a continuously differentiable function
$\kappa:\mathbf{R}^{D}\times\mathbf{R}^{K}\rightarrow\mathbf{R}^{D}$. Fix
$\eta_{0}\in\mathbf{R}^{K}$ and suppose there exists a unique point
$z_{0}\in\mathbf{R}^{D}$ such that $\kappa(z_{0},\eta_{0})=0$. Then there is a
function $z^{*}(\eta):\Omega\rightarrow\mathbf{R}^{D}$ defined on an open set
$\Omega\subset\mathbf{R}^{K}$ containing $z_{0}$ such that
$\kappa(z^{*}(\eta),\eta)=0$ for all $\eta\in\Omega$. Moreover, if
$\partial_{1}\kappa(z^{*}(\eta),\eta)$ is invertible, then $z^{*}(\eta)$ is
differentiable with Jacobian,

$\partial
z^{*}(\eta)=-\partial_{1}\kappa(z^{*}(\eta),\eta)^{-1}\cdot\partial_{2}\kappa(z^{*}(\eta),\eta).$

In the context of the dispatch model,
$z^{*}(\eta)=(x^{*}(\eta),\lambda^{*}(\eta))$. From here forth, we use $z$
instead of $(x,\lambda)$ whenever possible. The differentiability and
uniqueness assumptions follow directly from Assumptions 1, 2, and 6. The
gradient of the multi-value expansion planning objective
$J(\eta)=\gamma^{T}\eta+h(z^{*}(\eta))$ is then

$\displaystyle\nabla J(\eta)$ $\displaystyle=\gamma+\partial
z^{*}(\eta)^{T}\cdot\nabla h(\eta).$

#### Fast computation.

The Jacobians of the solution maps are only used as part of a matrix-vector
product; therefore, we need not compute $\partial z^{*}$ explicitly and,
hence, need not invert $\partial_{1}\kappa$ directly. Instead, we only need to
solve the $(M+N)$-dimensional linear system $(\partial_{1}\kappa)^{T}y=\nabla
h$, which can be accomplished efficiently either using a sparse direct
factorization or an indirect method [38].

###  5.2 Gradient descent

Assume the investment constraint set $\mathcal{H}$ is the continuous box
$\mathcal{H}^{\textrm{cont}}$. Theorem 1 implies that solving Problem (P1)
amounts to minimizing a differentiable function subject to box constraints. A
standard tool for solving this type of problem is (projected) gradient
descent, which iteratively updates the network investments using gradient
steps projected onto $\mathcal{H}$. Evaluating the gradient $\nabla J(\eta)$
requires solving the dispatch model with parameters $\eta$, computing
Jacobians of the KKT operator, and solving a linear system. We describe the
full approach, called implicit (projected) gradient descent in Algorithm 1.

Despite the fact that $h(z^{*}(\eta))$ is defined implicitly through the
solution to an optimization problem, gradient descent still provably converges
to a local optimum in polynomial time.

given: an initial point $\eta^{(0)}\in\mathcal{H}$ and step size $\alpha>0$
---
while: stopping criterion is not satisfied
| 1\. solve the dispatch model (1) and its dual (2) with parameter
$\eta^{(i)}$
| 2\. compute the gradient $\nabla J(\eta^{(i)})=\gamma+\partial
z^{*}(\eta^{(i)})^{T}\cdot\nabla h(z^{*}(\eta^{(i)}))$
| 3\. set
$\eta^{{(i+1)}}:=\mathbf{proj}_{\mathcal{H}}\left(\eta^{(i)}-\alpha\nabla
J(\eta^{(i)})\right)$

Algorithm 1 Projected gradient descent

######  Theorem 2 (Gradient descent convergence).

Assume $h$ is differentiable with Lipschitz continuous gradient. Then, for
sufficiently small step size, projected gradient descent (Algorithm 1) finds
an $\epsilon$-stationary point,

$\|\eta^{(i)}-\mathbf{proj}_{\mathcal{H}}\left(\eta^{(i)}-\alpha\nabla
J(\eta^{(i)})\right)\|_{2}^{2}\leq\epsilon,$

in $O(1/\epsilon^{2})$ iterations, where the complexity is independent of the
network parameter dimension $K$.

###### Proof.

See Appendix A.3. ∎

Since the dispatch model (1) is convex, the dispatch map $z^{*}(\eta)$ and its
gradient $\partial z^{*}(\eta)$ can be evaluated in time polynomial in the
problem dimension using, for example, an interior point method [33].
Therefore, each iteration of Algorithm 1 has polynomial complexity, and,
hence, the algorithm reaches an $\epsilon$-stationary point in polynomial
time.

Implicit gradient descent has been used for bilevel optimization in a number
of settings [28, 29, 39]. However, as far as we are aware, this method has yet
to be applied to power system expansion planning problems, save a recent work
applying gradient methods to emissions-aware problems [23]. The closest
related methods are sensitivity methods, e.g., [40], which use gradients in a
heuristic manner to solve classical grid expansion planning problems. These
methods, however, use duality to explicitly derive the gradient of operation
cost $c(x^{*}(\eta))$, and therefore cannot be generally applied to generic
bilevel objectives.

#### Interpretation as sensitivity-based investment.

In some cases, gradient descent has a natural interpretation as a sensitivity-
based investment process. In this process, planners observe the current state
of the network parameters $\eta^{(i)}$ and its dispatch outcome
$z^{*}(\eta^{(i)})$. They also observe $\partial z^{*}(\eta^{(i)})$ the
sensitivities of the various grid outcomes with respect to said network
parameters. For example, the sensitivity of operation cost $c(x^{*}(\eta))$
with respect to additional generation capacity at a node in the network is a
function of that node’s locational marginal price, and the sensitivity of
emissions is a function of the node’s locational marginal emission rate [41,
23, 42]. After observing these metrics, the planner makes a decision to invest
in slightly more or less of each network parameter. Then, the new state of the
network is revealed and the process repeats. This process is similar to how
the grid is incrementally upgraded and modified in practice in response to
changing economic and environmental signals.

#### Parallelization.

Algorithm 1 can be applied to stochastic multi-value planning problems (6) as
well. The gradient of (6) is

$\nabla J^{\textrm{stochastic}}(\eta)=\gamma+\frac{1}{S}\sum_{s=1}^{S}\partial
z_{s}^{*}(\eta)^{T}\cdot\nabla h(z_{s}^{*}(\eta)).$

The summands $\partial z_{s}^{*}(\eta)^{T}\cdot\nabla h(z_{s}^{*}(\eta))$ can
be computed in parallel over $S$ processors. If a machine with many cores is
available, this can drastically reduce the complexity of gradient descent
applied to large stochastic problems.

#### Stochastic gradients.

We can further accelerate gradient descent for stochastic multi-value planning
by sampling a subset of the scenarios $\mathcal{B}\subset\\{1,\ldots,S\\}$ and
estimating the gradient stochastically. This method scales well to problems
with hundreds or even thousands of scenarios and has been extensively studied
for solving large-scale nonconvex problems in machine learning [43]. We
describe this approach, stochastic gradient descent, in Algorithm 2.

#### Initialization and stopping criterion.

We use two simple initialization schemes: uniform random sampling from the
feasible set $\mathcal{H}$ and using the solution to the convex relaxation in
Section 6.2. When initializing using random sampling, we find it is often
helpful to sample multiple times—say, five to ten—and choose the sample with
the lowest objective value as an initialization. As a stopping criterion, we
terminate the algorithm when either: (i) $\eta^{(i)}$ is approximately
stationary as defined in Theorem 2, (ii) the objective value is sufficiently
close to the lower bound produced by the relaxation in Section 6.2, or (iii) a
fixed iteration limit is reached.

given: an initial point $\eta^{(0)}\in\mathcal{H}$, step size $\alpha>0$, and
a batch size $0<R<S$
---
while: stopping criterion is not satisfied
| 1\. sample $s_{1},\ldots,s_{R}$ uniformly at random from $\\{1,\ldots,S\\}$
| 2\. solve (5) and its dual with parameter $\eta^{(i)}$ for each scenario
$s_{1},\ldots,s_{R}$
| 3\. compute a stochastic estimate of the gradient
|  $\Delta^{(i)}=\gamma+\frac{1}{R}\sum_{r}\partial
z_{s_{r}}^{*}(\eta^{(i)})^{T}\cdot\nabla h(z_{s_{r}}^{*}(\eta^{(i)}))$
| 4\. set
$\eta^{{(i+1)}}:=\mathbf{proj}_{\mathcal{H}}\left(\eta^{(i)}-\alpha\Delta^{(i)}\right)$

Algorithm 2 Stochastic projected gradient descent

#### Integer constraints.

There is no straightforward extension of gradient descent to the integer-
constrained case. A simple heuristic is to solve the continuous problem, then
use probabilistic rounding [44] to generate a few integer solutions. In many
cases, this method works well and finds high quality solutions, especially
when the number of decisions is large.

##  6 Convex relaxation-based lower bound

We now show how to reformulate multi-value expansion planning as a nonconvex
quadratic program using strong duality and derive a simple convex relaxation
of this problem using the McCormick envelope. This convex relaxation gives a
lower bound on the optimal objective value of (P1) and often yields a good
initial network investment for the gradient algorithm. Finally, we show that
the relaxation is exact when the network parameters are binary constrained and
approximate when the network parameters are continuous and constrained over a
small box.

###  6.1 Strong duality reformulation

To reduce a bilevel optimization problem to an single level problem, we can
formulate the optimality conditions of the lower level problem as explicit
constraints in the upper level problem. Specifically, the upper level problem
can be written as joint optimization problem over $\eta$, $x$, and $\lambda$,
then either the KKT or strong duality conditions of the lower level problem
can be used to enforce lower level optimality. See [8, 45] for an overview and
[46] for a comprehensive treatment. In this work, we will use the strong
duality-based approaches because of their simplicity and favorable performance
in [8].

For a given $\eta\in\mathbf{R}^{K}$, weak duality implies that

$c(x)\geq u(\lambda,\eta),$ (14)

for any $x\in\mathbf{R}^{N}$ and $\lambda\in\mathbf{R}^{M}$ such that
$A(\eta)x\leq b(\eta)$ and $\lambda\geq 0$. Since the dispatch model (1) and
its dual (2) are convex and concave, respectively, strong duality implies that

$c(x^{*}(\eta))=u(\lambda^{*}(\eta),\eta),$

and, conversely, that if $c(x)=u(\lambda,\eta)$ for $x\in\mathbf{R}^{N}$ and
$\lambda\in\mathbf{R}^{M}$ such that $A(\eta)x\leq b(\eta)$ and $\lambda\geq
0$, then $x=x^{*}(\eta)$ and $\lambda=\lambda^{*}(\eta)$, i.e., $(x,\lambda)$
are optimal if and only if $c(x)=u(\lambda,\eta)$. Because (14) holds for any
$x\in\mathbf{R}^{N}$ and $\lambda\in\mathbf{R}^{M}$, the above statement
implies that $c(x)\leq u(\lambda,\eta)$ is equivalent to $(x,\lambda)$ being
optimal. This allows us to reformulate multi-value expansion planning as a
joint optimization problem over $\eta$, $x$, and $\lambda$. The
$\delta$-strong duality formulation of multi-value expansion planning is

$\begin{array}[]{ll}\textrm{minimize}&\gamma^{T}\eta+h(x,\lambda)\\\\[2.5pt]
\textrm{subject to}&c(x)\leq u(\lambda,\eta)+\delta\\\\[2.5pt] &A(\eta)x\leq
b(\eta)\\\\[2.5pt] &\lambda\geq 0\\\\[2.5pt] &\eta\in\mathcal{H},\end{array}$
(P2)

where the variables are $x\in\mathbf{R}^{N}$, $\lambda\in\mathbf{R}^{M}$, and
$\eta\in\mathbf{R}^{K}$. Above, $\delta\geq 0$ is a small constant used to
relax the strong duality constraint. When $\delta=0$, (P2) is equivalent to
the implicit form (P1) of multi-value expansion planning; any feasible point
$(x,\lambda,\eta)$ of (P2) is such that $\eta$ is feasible for (P1),
$x=x^{*}(\eta)$, and $\lambda=\lambda^{*}(\eta)$. Since the objective of both
problems is the same, this implies the two problems are equivalent.

In practice, we set $\delta$ to a small positive constant so that the
constraint $c(x)\leq u(\lambda,\eta)+\delta$ is strictly feasible; this is
done to guarantee that certain technical conditions hold for the relaxed
problem in Section 6.2. Fortunately, in this case, Problem (P2) is still a
close approximation of (P1), in the following sense.

######  Definition ($\epsilon$-approximation).

Consider two optimization problems,

$\min_{(\eta,y)\in\mathcal{C}_{1}}\
f(\eta,y)\qquad\textrm{and}\qquad\min_{(\eta,z)\in\mathcal{C}_{2}}\
g(\eta,z),$

which we refer to as (A) and (B), respectively, where
$\mathcal{C}_{1}\subset\mathbf{R}^{K+N_{1}}$ and
$\mathcal{C}_{2}\subset\mathbf{R}^{K+N_{2}}$. We say (B) is an
$\epsilon$-approximation to (A) over $\eta$ if for all
$(\eta,z)\in\mathcal{C}_{2}$, there exists $y\in\mathcal{Y}$ such that
$(\eta,y)\in\mathcal{C}_{1}$ and

$\displaystyle|f(\eta,y)-g(\eta,z)|\leq\epsilon.$

######  Lemma 3.

Suppose $h$ is $L$-Lipschitz continuous. For all $\epsilon>0$, there exists
some $\delta>0$ such that the $\delta$-strong duality formulation (P2) is an
$\epsilon$-approximation of (P1) over $\eta$, i.e., for all $(\eta,x,\lambda)$
feasible for (P2),

$|h(x,\lambda)-h(x^{*}(\eta),\lambda^{*}(\eta))|\leq\epsilon.$

###### Proof.

See Appendix B.1. ∎

Even when $h$ is convex, Problem (P2) is generally not convex because of the
bilinear terms $A(\eta)x$ and
$u(\lambda,\eta)=-\lambda^{T}b(\eta)-c^{\star}(-A(\eta)^{T}\lambda)$. If $c$
and $h$ are quadratic (or linear), then (P2) is a nonconvex quadratic problem
and is hard to solve in general. However, generic off-the-shelf solvers for
nonconvex quadratic programs—for example, a nonlinear interior point solver
such as Ipopt [47]—can often find good locally optimal solutions in reasonable
amounts of time.

###  6.2 McCormick relaxation

Consider the bilinear function $\phi(\alpha,\beta)=\alpha\beta$, where
$\alpha,\beta\in\mathbf{R}$. The function $\phi$ is nonconvex, but has an
intuitive convex (and concave) envelope when there are box constraints on
$\alpha$ and $\beta$.

######  Definition (Convex envelope).

Consider a function $\phi:\mathbf{R}^{N}\rightarrow\mathbf{R}$ and a convex
set $\mathcal{S}\subset\mathbf{R}^{N}$. Then the convex envelope of $\phi$ is
the largest convex under-estimator of $\phi$ over $\mathcal{S}$,

$\mathbf{conv}_{\mathcal{S}}(\phi)=\sup\\{f:f\textrm{ convex and
}f(x)\leq\phi(x)\textrm{ for all }x\in\mathcal{S}\\}.$

Similarly, the concave envelope of $\phi$ is the smallest concave over-
estimator of $\phi$ over $\mathcal{S}$,

$\mathbf{conc}_{\mathcal{S}}(\phi)=\sup\\{f:f\textrm{ concave and
}f(x)\geq\phi(x)\textrm{ for all }x\in\mathcal{S}\\}.$

Let
$\mathcal{B}=\\{\alpha:\alpha^{\min}\leq\alpha\leq\alpha^{\max}\\}\times\\{\beta:\beta^{\min}\leq\beta\leq\beta^{\max}\\}$
and suppose $(\alpha,\beta)\in\mathcal{B}$. The box constraints imply that the
product $(\alpha-\alpha^{\min})(\beta-\beta^{\min})$ is always nonnegative,
which means that

$\alpha\beta-\alpha\beta^{\min}-\alpha^{\min}\beta+\alpha^{\min}\beta^{\min}\geq
0,$

or, equivalently,

$\phi(\alpha,\beta)\geq\alpha\beta^{\min}+\alpha^{\min}\beta-\alpha^{\min}\beta^{\min}.$
(15)

Equation (15) gives an affine lower bound on $\phi(\alpha,\beta)$. Similarly,
the nonnegative product $(\alpha^{\max}-\alpha)(\beta^{\max}-\beta)$ generates
the lower bound

$\alpha\beta\geq\alpha^{\max}\beta+\alpha\beta^{\max}-\alpha^{\max}\beta^{\max}.$
(16)

Equations (15) and (16) define the convex envelope of $\phi$ over the
intervals $\alpha^{\min}\leq\alpha\leq\alpha^{\max}$ and
$\beta^{\min}\leq\beta\leq\beta^{\max}$. We can use similar logic to derive
the concave envelope of $\phi$ as well.

######  Lemma 4 (McCormick envelope, [48]).

The convex and concave envelopes of $\phi(\alpha,\beta)=\alpha\beta$ on
$\mathcal{B}$ are

$\displaystyle\mathbf{conv}_{\mathcal{B}}(\phi)$
$\displaystyle=\max\left(\alpha\beta^{\min}+\alpha^{\min}\beta-\alpha^{\min}\beta^{\min},\
\alpha^{\max}\beta+\alpha\beta^{\max}-\alpha^{\max}\beta^{\max}\right),$
$\displaystyle\mathbf{conc}_{\mathcal{B}}(\phi)$
$\displaystyle=\min\left(\alpha^{\min}\beta+\alpha\beta^{\max}-\alpha^{\min}\beta^{\max},\
\alpha\beta^{\min}+\alpha^{\max}\beta-\alpha^{\max}\beta^{\min}\right),$

respectively.

######  Remark 1.

For
$(\alpha,\beta)\in\mathcal{B}=[\alpha^{\min},\alpha^{\max}]\times[\beta^{\min},\beta^{\max}]$,
the bilinear constraint $y=\phi(\alpha,\beta)$ has convex relaxation
$\mathbf{conv}_{\mathcal{B}}(\alpha\beta)\leq
y\leq\mathbf{conc}_{\mathcal{B}}(\alpha\beta),$ which, for short, we write as
$y\in\mathbf{env}(\alpha\beta)$.

The McCormick envelope is a well-know relaxation for bilinear terms in
optimization problems; see, for example, [49] for an application to gas
expansion planning and [50, 51] for applications classical power system
expansion planning. Here, we use Lemma 4 to derive a convex relaxation on the
strong duality formulation of multi-value expansion planning when the
planner’s objective $h$ is convex. We first introduce three new variables
$Y_{ij}=A(\eta)_{ij}x_{j}$, $V_{ji}=-A(\eta)_{ij}\lambda_{i}$, and
$\mu_{i}=\lambda_{i}b(\eta)_{i}$ that each model a product of two variables.
We then introduce box constraints on any unbounded variables and, finally,
replace each of these products with upper and lower bounds on their convex and
concave envelopes, respectively. The $\delta$-relaxed formulation of multi-
value expansion planning is

$\begin{array}[]{ll}\textrm{minimize}&\gamma^{T}\eta+h(x,\lambda)\\\\[2.5pt]
\textrm{subject to}&c(x)\leq-\mu-c^{\star}(-V\mathbf{1})+\delta\\\\[2.5pt]
&Y\mathbf{1}\leq b(\eta)\\\\[2.5pt] &x^{\min}\leq x\leq x^{\max}\\\\[2.5pt]
&0\leq\lambda\leq\lambda^{\max}\\\\[2.5pt]
&Y_{ij}\in\mathbf{env}(A_{ij}(\eta)x_{j}),\quad i\in\\{1,\ldots,m\\},\
j\in\\{1,\ldots,n\\}\\\\[2.5pt]
&V_{ji}\in\mathbf{env}(A_{ij}(\eta)\lambda_{i}),\quad i\in\\{1,\ldots,m\\},\
j\in\\{1,\ldots,n\\}\\\\[2.5pt]
&\mu_{i}\in\mathbf{env}(\lambda_{i}b_{i}(\eta)),\quad
i\in\\{1,\ldots,m\\}\\\\[2.5pt] &\eta\in\mathcal{H},\end{array}$ (P3)

where the variables are $\eta\in\mathbf{R}^{K}$, $x\in\mathbf{R}^{N}$,
$\lambda\in\mathbf{R}^{M}$, $Y\in\mathbf{R}^{M\times N}$,
$V\in\mathbf{R}^{N\times M}$, and $\mu\in\mathbf{R}^{M}$. Here we drop the
subscript on the $\mathbf{conv}$ and $\mathbf{conc}$ operators and use the box
constraints implied by other constraints of the problem. Box constraints on
$A_{ij}(\eta)$ and $b_{i}(\eta)$ are implied from the constraint $\eta\in H$
and the linearity of $A$ and $b$. The box constraints on $x$ and $\lambda$
arise from Assumption 5 and are assumed to either arise from problem structure
or, if needed, can be selected as arbitrarily large values. Problem (P3) is
convex, has the same objective as Problem (P2), and relaxes its feasible set:
any feasible point of (P2) is also feasible for (P3).

#### Lower bound and initialization.

Since Problem (P3) is a convex relaxation of (P1), its optimal objective value
is a lower bound on that of the original multi-value problem (P1). Moreover,
its solution $\eta^{\textrm{relax}}\in\mathcal{H}$ is, in many cases, close to
the solution to (P1). Therefore, we can use $\eta^{\textrm{relax}}$ as an
initialization for Algorithm (1). In this context, we can interpret Algorithm
(1) as a local refinement applied to the output of our convex relaxation.

$1$$2$$3$$-2$$0$$2$$\alpha$$\beta$ Figure 3:  McCormick envelope of the
constraint $\phi=\alpha\beta$ over the box
$(\alpha,\beta)\in[1,3]\times[-2,2]$. projected onto the plane
$\\{(\alpha,\beta,\phi):\phi=1.5\\}$. (Black line) The nonconvex feasible set.
(Gray dashed lines) Linear inequalities defined by convex and concave
(McCormick) envelopes of $f(\alpha,\beta)=\alpha\beta$ over the boxed region.
(Gray shaded region) Convexified feasible set.

When the network parameters are binary box constrained
$\eta\in\mathcal{H}^{\textrm{int}}$, the relaxation (P3) is tight.
Intuitively, this is because the convex envelope of a bilinear function meets
the function itself at the endpoints of the interval constraints; see Figure
3. As a result, integer-constrained multi-value planning can be converted into
a convex mixed-integer program and solved with off-the-shelf methods, if
desired.

######  Theorem 5 (Tight integer relaxation).

Consider an integer constrained variable
$\alpha\in\\{\alpha^{\min},\alpha^{\max}\\}$ and a box constrained variable
$\beta\in[\beta^{\min},\beta^{\max}]$ (which may or may not have integer
constraints as well). Then for any $y\in\mathbf{R}$, the equation
$y=\alpha\beta$ holds if and only if $y\in\mathbf{env}(\alpha\beta).$ As a
consequence, if $\mathcal{H}=\mathcal{H}^{\textrm{int}}$, then the feasible
sets of (P2) and its relaxation (P3) are equal, i.e., the relaxation is tight.

###### Proof.

See Appendix B.2. ∎

When the network parameters are continuous, (P3) is not a tight relaxation in
general. However, we can derive an analogous result to the binary case when
the width of the box $\mathcal{H}^{\textrm{cont}}$ is small: as the width of
$H^{\textrm{cont}}$ goes to zero, the Hausdorff distance between the relaxed
set and the original set approaches zero, and so the solution to (P3)
approaches that of (P2).

######  Theorem 6 (Approximate continuous relaxation).

For all $\epsilon>0$, there exists some $\alpha>0$ such that when
$\|\eta^{\min}-\eta^{\max}\|_{\infty},\leq\alpha,$ then (P3) is a
$\epsilon$-approximation of (P2) over $\eta$.

###### Proof.

See Appendix B.3 ∎

Theorem 6 implies that the lower bound from the solution to (P3) approaches
the optimal value of (P2) as $\|\eta^{\min}-\eta^{\max}\|_{\infty}\rightarrow
0$. As a result, (P3) can be used as part of a branch-and-bound scheme to
solve multi-value expansion planning exactly [52], albeit in exponential time
in the worst-case.

#### Reformulation and relaxation for stochastic problems.

Both the strong duality formulation (P2) and its relaxation (P3) are
straightforward to extend to stochastic multi-value expansion planning. In
this case, the strong-duality formulation is a joint optimization problem over
the network parameters $\eta$ and the primal variables $x_{1},\ldots,x_{S}$
and dual variables $\lambda_{1},\ldots,\lambda_{S}$ for each scenario. The
resulting optimization problems have $O(SN)$ variables and $O(SM)$
constraints; this may make them prohibitively expensive to solve, especially
for the nonconvex strong-duality formulation and for problems with integer-
constrained network investments.

##  7 Computational experiments

In this section, we describe our software implementation of the methods
developed in Sections 5 and 6 and analyze their performance on various sized
synthetic models of the U.S. Western Interconnect in four experiments. First,
we examine the convergence properties of Algorithm 1 and show how it can be
used to refine solutions produced by the convex relaxation. Second, we compare
the performance of Algorithm 1 with that of a traditional interior point
method applied directly to the nonconvex strong duality formulation of multi-
value expansion planning. Third, we study how the runtime of Algorithm 1
scales with both the dimension of the problem and the number of scenarios.
Finally, we explore the resulting solutions on a planning problem with a long
time horizon.

In our experiments, we solve multi-value expansion planning problems using the
same dispatch model from Section 4. Unlike Section 4, however, we jointly plan
both transmission and generation capacities, i.e., $\eta=(\beta,g^{\max})$ in
the notation of Section 4. In each experiment, we solve the emissions-aware
expansion planning problems from (9) with an emissions penalty of
$400-equivalent per metric ton $\textrm{CO}_{2}$. In the first three
experiments, Section 7.1 through Section 7.3, we treat a single hour as one
scenario, ignoring dynamic devices like batteries for simplicity, using the
dispatch model in (7). In the large example in Section 7.4, we use the dynamic
dispatch model from (11) with 24 hourly periods, treating a full day as a
single scenario and jointly planning generation, transmission, and storage
capacities.

#### Software implementation.

We implement the methods developed in Sections 5 and 6 in Julia [53]. We build
the dispatch model from Section 4 using `Convex.jl` [54] and solve it using
Mosek [55] to implement the dispatch map $z^{*}(\eta)$. The Jacobian $\partial
z^{*}(\eta)$ is implemented manually using a modular construction that
supports a variety of devices in a (DC) linearized power flow model. To
compute the gradient $\nabla h(z)$ for various planning objectives, we use the
automatic differentiation library `Zygote.jl` [56]. Finally, we model the
strong duality formulation of multi-value expansion planning and its
relaxation using `JuMP.jl` [57]. The nonconvex strong duality formulation is
solved with Ipopt [47], and the relaxation can be solved with Mosek. For
stochastic problems from Section 3.4, we parallelize (multithread) the
evaluation of both the dispatch map and its Jacobian. All our experiments are
run on a single AMD EPYC 7763 processor with 64 cores (128 virtual cores) and
256 GB of RAM at NERSC.

#### Western Interconnect dataset.

We use PyPSA-USA [58], an open source model of the U.S. Western Interconnect,
to load and process network models. The dataset includes a synthetic network
with 4,786 nodes with 8760 hours of data from 2019. The network granularity
can also be adjusted using a cluster algorithm down to as few as 30 nodes,
allowing us to create planning cases with dispatch models of different
dimensions. For example, clustering the network to 100 nodes results in a
problem with 222 transmission lines and 379 generators, whereas clustering to
500 nodes results in 1108 and 1356 generators, respectively. In all of our
experiments, we consider a future grid scenario with 50% load growth relative
to 2019. We reduce generation and transmission capacities by 20% to account
for retirements and outages, leaving all other data the same. The dataset
includes realistic fuel costs, investment costs, hypothetical expansion
capacities, electrical network parameters, and energy storage devices.
However, we emphasize that our grid model is not a finely calibrated dataset
suitable for long-term policy recommendations, but rather a realistic grid
model for testing the methods developed in this paper. Additional details on
the dataset are available in Appendix C.

###  7.1 Convergence analysis

In our Experiment 1, we analyze the convergence of gradient descent for
varying step sizes. We test this on two networks, one with 100 nodes and
another 500 nodes, both for a deterministic planning problem and a stochastic
one with four scenarios. In each example, we use the convex relaxation from
Section 6.2 to generate a lower bound on the optimal objective value and to
initialize Algorithm 1.

In this experiment and Experiment 2, we consider the suboptimality gap,

$\textrm{gap}=\frac{J^{(i)}}{J^{\mathrm{lower}}}-1,$ (17)

where $J^{{(i)}}$ is the objective value of the gradient algorithm and
$J^{\mathrm{lower}}$ is the lower bound given by the solution to Problem (P3).
The lower bound may not be tight, in which case the suboptimality gap at a
global optima will be greater than zero. However, a suboptimality gap of zero
guarantees global optimality, and a suboptimality gap of $\alpha$ guarantees
that objective value of the solution found is no worse than $1+\alpha$ times
the objective value at a global optimum.

We test step sizes from $10^{-3}$ to $10^{-1}$ for each case, displaying our
results in Figure 4. We terminate gradient descent after 2000 iterations,
although it is clear that many of the runs converge more quickly. As expected
by Theorem 2, gradient descent always converges to a stationary point when the
step size is sufficiently small. Larger step sizes may not necessarily
converge (see the 100 node, 4 scenario case), but generally converge faster
than smaller step sizes. Regardless, gradient descent with every tested step
size improves upon the relaxation-based solution. Finally, we note that the
suboptimality gap (17) is less than 10% in the four tested cases.

Figure 4:  Convergence results for gradient descent with varying step sizes.
(Dashed black line) Lower bound produced by the relaxation from Section 6.2.
(Solid black line) Initial choice of the network parameter, produced by the
same relaxation. The width of the line is the time to solve the relaxation.
(Solid colored lines) Loss curves for gradient descent. Vertical dashed lines
mark every 250th iteration.

###  7.2 Comparison to local interior point methods

In Experiment 2, we compare gradient descent to solving the strong duality
formulation directly using a local interior point solver. Specifically, we
apply the Ipopt solver to Problem (P2), comparing its objective value and
runtime to running 3,000 iterations of gradient descent initialized at the
solution to the convex relaxation (P3) for problems with 100 nodes and 1, 4,
and 16 scenarios. For each number of scenarios, we sample 10 different sets of
hourly data to obtain robust estimates of performance. Because the interior
point solver’s runtime can vary significantly from case to case, we limit its
runtime to 10 minutes times the number of scenarios. If Ipopt has not
converged, we use the best objective value solution found so far, with a
relatively high feasibility tolerance of $1.0$.

We show the results of this experiment in Figure 5, normalizing objective
values using the convex relaxation-based lower bound. For 1 scenario problems,
the interior point algorithm outperforms gradient descent both in objective
value and runtime by a small margin. For 4 scenario problems, the interior
point algorithm and gradient descent find similar quality solutions, but
gradient descent takes significantly less time. Finally, for 16 scenario
problems, gradient descent significantly outperforms the interior point
algorithm both in objective value and runtime. Note that the interior point
algorithm is terminated early for every instance of the 16 scenario problem.
This highlights the practical importance of developing planning algorithms
that scale well to many scenarios.

Figure 5:  Performance comparison between a local interior point solver (IP),
gradient descent (GD), and the solution to the relaxation in Section 6.2 (RL).
(Left side) Box plots of final objective values (divided by the lower bound)
for different cases with 1, 4, or 16 scenarios. Box bottom and top lines are
the 25th and 75th percentiles, respectively, and midlines are medians. (Right
side) Box plot of algorithm runtimes for different cases. For visibility, box
bottom and top lines are set to be 10 minutes apart if the IQR is less than 10
minutes.

###  7.3 Computational scaling

One of the primary advantages of using gradient descent to solve planning
problems is that it scales well to large networks with many scenarios. In
Experiment 3, we verify this by empirically studying the computational time
per iteration and, separately, the number of iterations required to converge
to a stationary point.

#### Computational performance.

First, we study the computational time per iteration for cases with network
size varying from 200 to 1000 nodes, solving the resulting planning problems
with a single scenario. Each case is run 3 times for 500 iterations, and the
per-iteration time averaged across the duration of the run. We plot these
results in Figure 6 (left panel). The per-iteration runtime appears to grow
approximately linearly with network size. This roughly agrees with theoretical
intuition: the main computational bottleneck is factorizing the sparse KKT
matrix while solving the dispatch model, which grows near-linearly with
problem dimension for many structured matrices (see, for example, [59]).

Second, we study how the per-iteration computational time scales with number
of scenarios. Since we are able to parallelize the majority of the algorithm —
solving the dispatch model and evaluating gradients for different scenarios —
we expect a near-constant scaling on a machine with many processors. We use a
100-node network and scale the number of scenarios from 1 to 32, again running
each case 3 times. Experiments are run with 64 threads on a 128 core (64
physical cores) cluster machine. Therefore, in principle, runtime should be
constant as a function of the number of scenarios, barring competing
parallelism between linear algebra libraries and the solver. The results are
plotted in Figure 6 (right panel). For convenience, we plot worst (linear) and
best (constant) case theoretical scaling in the same panel. Results
demonstrate that evaluating gradients in parallel is significantly faster than
a serial implementation and achieves a near-constant scaling with the number
of scenarios.

Figure 6:  Per-iteration computational time for various cases. Black dots
denote runtimes of individual cases and blue solid lines are averaged runtime.
(Left) Time per iteration as a function of network size. (Right) Time per
iteration as a function of number of scenarios. Serial performance (orange
dotted line), i.e., linear scaling and ideal parallel performance (green
dotted line), i.e., constant scaling, are plotted for reference. Gradient
descent benefits significantly from parallelization.

#### Convergence rate.

Finally, we study the number of iterations until convergence. Here, we define
convergence as when the best objective value so far decreases by less than
0.1% over 100 iterations. If the algorithm has not converged after 5,000
iterations, we terminate it early. We consider a variety of cases as before;
we vary the number of nodes from 200 to 1000 on single scenario cases, then
vary the number of scenarios from 1 to 16 on 100 node cases. For each network
size, we solve 5 distinct scenarios separately to ensure our results are not
an artifact of a particular time period; for each number of scenarios, we also
solve 10 distinct cases for the same reasons. For each individual case, we
select sweep step sizes from $10^{-3}$ to $10^{-1}$ and choose the step size
with the best loss at convergence. We display these results in Figure 7. In
general, we observe a slight increase in the number of iterations required for
convergence as a function of either network size or number of scenarios. This
may suggest that larger cases have gradients with a larger Lipschitz constant.
In practice, we believe this could be addressed through more sophisticated
problem preconditioning or non-constant step size schedules.

Figure 7:  Iterations until convergence for different cases on the PyPSA-USA
network. Black dots denote individual cases and blue solid lines are averages.
(Left) Iterations until convergence as a function of network size. (Right)
Iterations until convergence as a function of number of scenarios. Gradient
descent appears to converge within a constant number of iterations, as
suggested by the theoretical analysis in Section 5. Faster convergence rates
for large cases may be indicative of a shrinking Lipschitz constant.

###  7.4 Large example

In Experiment 4, we demonstrate the potential utility of our method by solving
a large case study. We now consider dynamic dispatch models with multiple time
periods. Each dispatch problem is solved over $T=24$ hours coupled by battery
storage systems. We then solve two expansion planning problems using $S=32$
scenarios (768 hours total) using the 100-node PyPSA-USA network: an
emissions-aware problem with a penalty of $400 per metric ton CO2 (as in the
previous experiments) and a cost-only problem without a carbon penalty. We
specifically consider the 16 peak load days and the 16 peak renewable days as
our scenarios. In both problems we jointly plan transmission, generation, and
storage.

Because we encounter memory issues when solving the entire problem, we use
stochastic gradient descent as described in Algorithm 2 with a batch size of
$R=4$. Instead of sampling batches randomly, we randomly partition the
scenarios into batches and iterate through each batch serially, as is standard
practice in machine learning. We set the step size to $0.01$ and run the
algorithm for 1000 iterations. In both cases, the algorithm converges after
just a few hundred iterations. The final objective values are within 2% and
33% of the lower bound for the cost-only and emissions-aware cases,
respectively.

We plot the resulting transmission, storage, and generation capacity
expansions for both cases in Figure 8. As expected, the emissions-aware
problem invests in significantly more renewable, transmission, and storage
capacity. In order to understand each solution qualitatively, we also plot the
generation dispatch for the scenario with the highest peak load in Figure 9.
The results show that greater investment in both transmission and storage
enable higher renewable penetration, reducing the need for gas capacity. We
reemphasize that these conclusions are not to be interpreted as actual
forecasts or policy suggestions for the Western U.S., just as a demonstration
of the potential applicability of our methods to large scale case studies.

Quantitatively, the cost of investment is higher in the emissions-aware case:
$\$44.3/\textrm{MWh}$ as opposed to $\$15.9/\textrm{MWh}$ in the cost-only
case. The greater investment in renewable capacity leads to a significant
reduction in operation cost ($\$5.0/\textrm{MWh}$ versus
$\$16.3/\textrm{MWh}$) and carbon intensity ($47.4\
\textrm{kgCO}_{2}/\textrm{MWh}$ versus $186.8\
\textrm{kgCO}_{2}/\textrm{MWh}$). Both cases improve upon the 2019 base case
carbon intensity of $348.3\ \textrm{kgCO}_{2}/\textrm{MWh}$, but the
emissions-aware case achieves an additional $40.0\%$ reduction in carbon
intensity at an additional $\$17.1/\textrm{MWh}$ in total cost.

Figure 8:  Expanded capacity for transmission, storage, and generation using
cost-only (orange) and emissions-aware (green) expansion planning models.
Initial capacities are shown in blue, and non-expandable resources in gray.
Emissions-aware expansion leads to significantly more transmission, storage,
and generation investment, even in a market dispatched solely on cost.  Figure
9:  Generator outputs on the peak summer load day using new capacities from
the cost-only and emissions-aware expansions. (Black line) Total hourly load
across the network. (Shaded solid regions) Generator outputs by resource type.
(Transparent yellow region) Load shift from battery charging or discharging.
Cost-only investment reduces carbon intensity by 46.4% relative to the 2019
baseline, whereas emissions-aware investment leads to a 86.4% reduction
through greater renewable, transmission, and storage investment.

##  8 Conclusion

We consider multi-value expansion planning, a general bilevel optimization
framework for jointly optimizing generation, transmission, and storage assets
in settings where the planner’s objective may differ from the objective used
to dispatch power. We give several examples of multi-value expansion planning,
such as reducing emissions, even in the absence of a carbon tax, or maximizing
renewable profits. In general, this problem is high dimensional, nonconvex,
and may need to be solved across many scenarios. Therefore, we propose
applying gradient descent to the implicit form of multi-value expansion
planning, which converges to a local stationary point in polynomial time. We
demonstrate, both theoretically and empirically, that this method is scalable
to high dimensional problems with many scenarios. Moreover, we show gradient
descent outperforms commonly used strong duality-based reformulation methods
even on modest sized cases.

One strength of our work is that the multi-value formulation is fairly general
and can be extended to a variety of related expansion planning problems. For
example, we can consider planning electricity system expansions jointly with
other resources, e.g., gas or water, with no change to our gradient algorithm
other than the implementation of $x^{*}(\eta)$ and $\partial x^{*}(\eta)$. We
can also consider more advanced variants of expansion planning, such as
problems with multiple investment periods.

## Acknowledgements

This research used resources of the National Energy Research Scientific
Computing Center (NERSC), a U.S. Department of Energy Office of Science User
Facility located at Lawrence Berkeley National Laboratory, operated under
Contract No. DE-AC02-05CH11231 using NERSC award DDR-ERCAP0026889.

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research,
Department of Energy Computational Science Graduate Fellowship under Award
Number DE-SC0021110.

## References

  * [1] Eric D Larson, Chris Greig, Jesse D Jenkins, Erin N Mayfield, Stephen W Pascala, and others.  Net-zero america: Potential pathways, infrastructure, and impacts: Final report.  Princeton University, 2021.
  * [2] Paul L Joskow.  Transmission capacity expansion is needed to decarbonize the electricity sector efficiently.  Joule, 4(1):1–3, January 2020.
  * [3] Patrick R Brown and Audun Botterud.  The value of inter-regional coordination and transmission in decarbonizing the US electricity system.  Joule, 5(1):115–134, January 2021.
  * [4] Tom Brown, David Schlachtberger, Alexander Kies, Stefan Schramm, and Martian Greiner.  Synergies of sector coupling and transmission reinforcement in a cost-optimised, highly renewable european energy system.  Energy, 160:720–739, October 2018.
  * [5] Paul L Joskow.  Facilitating transmission expansion to support efficient decarbonization of the electricity sector.  Economics of Energy & Environmental Policy, 10(2), 2021.
  * [6] Mohamed L Awad, Keith E Casey, Anna S Geevarghese, Jeffrey C Miller, A Farrokh Rahimi, Anjali Y Sheffrin, Mingxia Zhang, Eric Toolson, Glenn Drayton, Benjamin F Hobbs, and F A Wolak.  Using market simulations for economic assessment of transmission upgrades: Application of the california ISO approach.  Restructured Electric Power Systems, pages 241–270, 2010.
  * [7] Lina P Garces, Antonio J Conejo, Raquel Garcia-Bertrand, and Ruben Romero.  A bilevel approach to transmission expansion planning within a market environment.  IEEE Transactions on Power Systems, 24(3):1513–1522, August 2009\.
  * [8] M Hosein Zare, Juan S Borrero, Bo Zeng, and Oleg A Prokopyev.  A note on linearized reformulations for a class of bilevel linear integer problems.  Annals of Operations Research, 272(1-2):99–117, January 2019.
  * [9] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  Nature, 521(7553):436–444, May 2015.
  * [10] Nikolaos E Koltsaklis and Athanasios S Dagoumas.  State-of-the-art generation expansion planning: A review.  Applied energy, 230:563–589, November 2018.
  * [11] Meisam Mahdavi, Carlos Sabillon Antunez, Majid Ajalli, and Ruben Romero.  Transmission expansion planning: Literature review and classification.  IEEE systems journal, 13(3):3129–3140, September 2019.
  * [12] Mohammad Reza Sheibani, Gholam Reza Yousefi, Mohammad Amin Latify, and Sarineh Hacopian Dolatabadi.  Energy storage system expansion planning in power systems: a review.  IET renewable power generation, 12(11):1203–1221, August 2018.
  * [13] Isaac-Camilo Gonzalez-Romero, Sonja Wogrin, and Tomas Gómez.  Review on generation and transmission expansion co‐planning models under a market environment.  IET Generation, Transmission and Distribution, 14(6):931–944, March 2020.
  * [14] Enzo E Sauma and Shmuel S Oren.  Proactive planning and valuation of transmission investments in restructured electricity markets.  Journal of regulatory economics, 30(3):358–387, December 2006.
  * [15] David Pozo, Javier Contreras, and Enzo Sauma.  If you build it, he will come: Anticipative power transmission planning.  Energy economics, 36:135–146, March 2013.
  * [16] David Pozo, Enzo E Sauma, and Javier Contreras.  A three-level static MILP model for generation and transmission expansion planning.  IEEE transactions on power systems : a publication of the Power Engineering Society, 28(1):202–210, February 2013.
  * [17] David Pozo, Enzo Sauma, and Javier Contreras.  Basic theoretical foundations and insights on bilevel models and their applications to power systems.  Annals of operations research, 254(1-2):303–334, July 2017.
  * [18] Francisco D Munoz, Benjamin F Hobbs, and Saamrat Kasina.  Efficient proactive transmission planning to accommodate renewables.  In 2012 IEEE Power and Energy Society General Meeting, pages 1–7. IEEE, July 2012.
  * [19] Luis Baringo and Antonio J Conejo.  Wind power investment within a market environment.  Applied Energy, 88(9):3239–3247, September 2011.
  * [20] Lajos Maurovich-Horvat, Trine K Boomsma, and Afzal S Siddiqui.  Transmission and wind investment in a deregulated electricity industry.  IEEE Transactions on Power Systems, 30(3):1633–1643, May 2015.
  * [21] Tohid Akbari and Mohammad Tavakoli Bina.  Merchant transmission capacity investment: A mathematical program with equilibrium constraint formulation.  Electric Power Components & Systems, 44(1):82–89, January 2016\.
  * [22] Bolun Xu.  Batteries in Electricity Markets: Economic Planning and Operations.  PhD thesis, University of Washington, 2018.
  * [23] Anthony Degleris, Lucas Fuentes Valenzuela, Abbas El Gamal, and Ram Rajagopal.  Emissions-aware electricity network expansion planning via implicit differentiation.  In NeurIPS 2021 Workshop on Tackling Climate Change with Machine Learning. Climate Change AI, 2021.
  * [24] Salvador Pineda, Henrik Bylling, and Juan M Morales.  Efficiently solving linear bilevel programming problems using off-the-shelf optimization software.  Optimization and engineering, 19(1):187–211, March 2018.
  * [25] Gonzalo Constante-Flores, Antonio J Conejo, and S Constante-Flores.  Solving certain complementarity problems in power markets via convex programming.  Top, 30(3):465–491, October 2022.
  * [26] Pablo Garcia-Herreros, Lei Zhang, Pratik Misra, Erdem Arslan, Sanjay Mehta, and Ignacio E Grossmann.  Mixed-integer bilevel optimization for capacity planning with rational markets.  Computers & chemical engineering, 86:33–47, March 2016.
  * [27] Xiaohu Zhang, Di Shi, Zhiwei Wang, Zhe Yu, Xinan Wang, Desong Bian, and Kevin Tomsovic.  Bilevel optimization based transmission expansion planning considering phase shifting transformer.  In 2017 North American Power Symposium (NAPS), pages 1–6, September 2017.
  * [28] Fabian Pedregosa.  Hyperparameter optimization with approximate gradient.  arXiv [stat.ML], pages 737–746, February 2016.
  * [29] Jonathan Lorraine, Paul Vicol, and D Duvenaud.  Optimizing millions of hyperparameters by implicit differentiation.  International Conference on Artificial Intelligence and Statistics, 108:1540–1552, November 2019.
  * [30] Ling Yang, Zhilong Zhang, Yang Song, Shenda Hong, Runsheng Xu, Yue Zhao, Wentao Zhang, Bin Cui, and Ming-Hsuan Yang.  Diffusion models: A comprehensive survey of methods and applications.  arXiv [cs.LG], September 2022.
  * [31] Asen L Dontchev and R Tyrrell Rockafellar.  Implicit functions and solution mappings: A view from variational analysis, volume 616 of Springer Series in Operations Research and Financial Engineering.  Springer, New York, NY, 2 edition, June 2014.
  * [32] Sébastien Bubeck.  Convex optimization: Algorithms and complexity.  Foundations and Trends® in Machine Learning, 8(3-4):231–357, 2015\.
  * [33] Stephen Boyd and Lieven Vandenberghe.  Convex Optimization.  Cambridge University Press, Cambridge, England, March 2004.
  * [34] Mario V F Pereira, Leontina Mvg Pinto, Sergio H F Cunha, and Gerson C Oliveira.  A decomposition approach to automated generation/transmission expansion planning.  IEEE Transactions on Power Apparatus and Systems, PER-5(11):30–31, November 1985.
  * [35] Paul Gribik, William W Hogan, and Susan L Pope.  Market-clearing electricity prices and energy uplift.  Technical report, Harvard Electricity Policy Group, 2007.
  * [36] Len Garver.  Transmission network estimation using linear programming.  IEEE transactions on power apparatus and systems, PAS-89(7):1688–1697, September 1970.
  * [37] Ruben Romero and Alcir Monticelli.  A hierarchical decomposition approach for transmission network expansion planning.  IEEE transactions on power systems : a publication of the Power Engineering Society, 9(1):373–380, February 1994.
  * [38] Lloyd N Trefethen and David Bau, III.  Numerical Linear Algebra.  Society for Industrial and Applied Mathematics, Philadelphia, PA, January 1997.
  * [39] Priya L Donti, Aayushya Agarwal, Neeraj Vijay Bedmutha, Larry Pileggi, and J Zico Kolter.  Adversarially robust learning for security-constrained optimal power flow.  arXiv [math.OC], November 2021.
  * [40] Mario Pereira and Leontina Pinto.  Application of sensitivity analysis of load supplying capability to interactive transmission expansion planning.  IEEE transactions on power apparatus and systems, PAS-104(2):381–389, February 1985.
  * [41] Pablo Ruiz and Aleksandr Rudkevich.  Analysis of marginal carbon intensities in constrained power networks.  Hawaii International Conference on System Sciences, pages 1–9, January 2010.
  * [42] Lucas Fuentes Valenzuela, Anthony Degleris, Abbas El Gamal, Marco Pavone, and Ram Rajagopal.  Dynamic locational marginal emissions via implicit differentiation.  IEEE transactions on power systems : a publication of the Power Engineering Society, 39:1138–1147, February 2023.
  * [43] Léon Bottou.  Large-scale machine learning with stochastic gradient descent.  In Proceedings of COMPSTAT’2010, pages 177–186, Heidelberg, 2010\. Physica-Verlag HD.
  * [44] Prabhakar Raghavan and Clark D Tompson.  Randomized rounding: A technique for provably good algorithms and algorithmic proofs.  Combinatorica. An International Journal on Combinatorics and the Theory of Computing, 7(4):365–374, December 1987.
  * [45] Thomas Kleinert and Martin Schmidt.  Why there is no need to use a big-M in linear bilevel optimization: a computational study of two ready-to-use approaches.  Computational management science, 20(1):1–12, December 2023.
  * [46] Stephan Dempe.  Foundations of bilevel programming.  Nonconvex Optimization and Its Applications. Springer, New York, NY, 2002 edition, May 2002.
  * [47] Andreas Wächter and Lorenz T Biegler.  On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming.  Mathematical programming, 106(1):25–57, March 2006.
  * [48] Faiz A Al-Khayyal and James E Falk.  Jointly constrained biconvex programming.  Mathematics of operations research, 8(2):273–286, May 1983.
  * [49] Conrado Borraz-Sánchez, Russell Bent, Scott Backhaus, Hassan Hijazi, and Pascal Van Hentenryck.  Convex relaxations for gas expansion planning.  INFORMS journal on computing, 28(4):645–656, November 2016.
  * [50] Quentin Ploussard, Luis Olmos, and Andres Ramos.  An operational state aggregation technique for transmission expansion planning based on line benefits.  IEEE transactions on power systems : a publication of the Power Engineering Society, 32(4):2744–2755, July 2017.
  * [51] Shahin Goodarzi, Mohsen Gitizadeh, and Ali Reza Abbasi.  Efficient linear network model for TEP based on piecewise McCormick relaxation.  IET Generation, Transmission and Distribution, 13(23):5404–5412, December 2019.
  * [52] Ralf Lenz and Felipe Serrano.  Tight convex relaxations for the expansion planning problem.  Journal of optimization theory and applications, 194(1):325–352, July 2022.
  * [53] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah.  Julia: A fresh approach to numerical computing.  SIAM Review, 59(1):65–98, January 2017.
  * [54] Madeleine Udell, Karanveer Mohan, David Zeng, Jenny Hong, Steven Diamond, and Stephen Boyd.  Convex optimization in julia.  In 2014 First Workshop for High Performance Technical Computing in Dynamic Languages, pages 18–28. ieeexplore.ieee.org, November 2014.
  * [55] MOSEK ApS.  Mosek.jl, 2022.
  * [56] Michael Innes.  Don’t unroll adjoint: Differentiating SSA-form programs.  CoRR, abs/1810.07951, October 2018.
  * [57] Iain Dunning, Joey Huchette, and Miles Lubin.  JuMP: A modeling language for mathematical optimization.  SIAM Review, 59(2):295–320, January 2017.
  * [58] Kamran Tehranchi, Trevor Barnes, and Ines Azevedo.  PyPSA-USA: An open-source optimisation model of the US power system, 2023.
  * [59] Michael B Cohen, Jonathan Kelner, Rasmus Kyng, John Peebles, Richard Peng, Anup B Rao, and Aaron Sidford.  Solving directed laplacian systems in nearly-linear time through sparse LU factorizations.  In 2018 IEEE 59th Annual Symposium on Foundations of Computer Science (FOCS), pages 898–909. IEEE, October 2018.
  * [60] Xingyu Zhou.  On the fenchel duality between strong convexity and lipschitz continuous gradient.  arXiv [math.OC], March 2018.
  * [61] Boris T Poljak.  Introduction to optimization.  Optimization Software, 1987.
  * [62] G W Stewart.  On the continuity of the generalized inverse.  SIAM journal on applied mathematics, 17(1):33–45, January 1969\.
  * [63] Ernest K Ryu and Stephen Boyd.  Primer on monotone operator methods.  Applied and Computational Mathematics, 15(1):3–43, 2016.
  * [64] Danish Energy Agency and Energinet.  Technology data - energy storage.  Technical report, Danish Energy Agency, 2018.

##  Appendix A Technical results on gradient descent

###  A.1 Definitions

######  Definition (Operator norm).

The operator norm of a matrix $A\in\mathbf{R}^{M\times N}$ is

$\|A\|_{\mathrm{op}}=\sup_{\|x\|_{2}\leq 1}\ \|Ax\|_{2}.$

######  Definition (Lipschitz continuity).

A function $f:\Omega\rightarrow\mathbf{R}^{M}$ is $L$-Lipschitz continuous if
for all $x,y\in\Omega$,

$\|f(x)-f(y)\|_{2}\leq L\|x-y\|_{2}.$

Similarly, a matrix-valued function $A:\Omega\rightarrow\mathbf{R}^{M\times
N}$ is $L$-Lipschitz continuous if for all $x,y\in\Omega$,

$\|A(x)-A(y)\|_{\mathrm{op}}\leq L\|x-y\|_{2}.$

###  A.2 Unique dual map

######  Lemma 7.

The solution to the dual dispatch model (2) is unique.

###### Proof.

By Assumption 2, the conjugate function $c^{\star}$ is both differentiable
with $(1/\epsilon)$-Lipschitz continuous gradient and $(1/L)$-strongly convex
[60]. Recall the objective of (2) is,

$u(\lambda)=-\lambda^{T}b(\eta)-c^{\star}(-A(\eta)^{T}\lambda).$

Now consider two optimal solutions
$\lambda_{1}^{*},\lambda_{2}^{*}\in\mathbf{R}^{M}$. Since $-c^{\star}$ is
strongly concave, it must be that

$A(\eta)^{T}\lambda_{1}^{*}=A(\eta)^{T}\lambda_{2}^{*}.$ (18)

By complementary slackness, $(\lambda_{i})^{*}_{m}>0$ only if the constraint
$\tilde{a}(\eta)^{T}x^{*}(\eta)\leq b(\eta)$ is binding. Assumption 3 tells us
that these rows of $A(\eta)$ are linearly independent, so the sub-matrix of
$A(\eta)^{T}$ corresponding to the non-zero entries of $\lambda_{1}^{*}$ and
$\lambda_{2}^{*}$ has linearly independent columns. Therefore, (18) implies
that $\lambda_{1}^{*}=\lambda_{2}^{*}$, i.e., the solution is unique. ∎

###  A.3 Gradient descent convergence

Theorem 2 is a restatement of the following result.

######  Theorem 8 (Gradient descent convergence, full).

Let $h$ have Lipschitz continuous gradient. Then $J(\eta)$ has Lipschitz
continuous gradient with Lipschitz constant $L>0$ and a lower bound $J^{*}\leq
J(\eta)$ for all $\eta\in\mathcal{H}$. Now suppose Algorithm 1 is run with
$\alpha\leq 1/L$. Then there exists an iteration

$i\leq\left\lceil\frac{2(J(\eta^{1})-J^{*})}{L\epsilon^{2}}\right\rceil$

such that $\eta^{i}$ is $\epsilon$-stationary,

$\|\eta^{i}-\mathbf{proj}_{\mathcal{H}}\left(\eta^{i}-\alpha\nabla
J(\eta^{i})\right)\|_{2}\leq\epsilon.$

At a high level, our approach is to show that the objective function
$J(\eta)=\gamma^{T}\eta+h(z^{*}(\eta))$ has Lipschitz continuous gradients,
then apply a standard analysis for the convergence of projected gradient
descent for Lipschitz functions, e.g., as in [28] or, for a more comprehensive
treatment, [61]. Since we assume the Lipschitz continuity of $\nabla h(z)$,
the main technical challenge is showing that $\partial z^{*}(\eta)$ is
Lipschitz continuous as well (in the metric induced by the operator norm).
Fortunately, this follows from our initial smoothness and regularity
assumptions in Section 3.1.

######  Lemma 9.

Let $D=N+M$ and consider the KKT operator $\kappa$ defined in Equation 13 Then
the matrices $\partial_{1}\kappa(z,\eta)\in\mathbf{R}^{D\times D}$ and
$\partial_{2}\kappa(z,\eta)\in\mathbf{R}^{D\times K}$ are Lipschitz
continuous.

###### Proof.

Recall that,

$\kappa(z,\eta)=\begin{bmatrix}\nabla c(x)+A(\eta)^{T}\lambda\\\
\mathbf{diag}(\lambda)(A(\eta)x-b(\eta))\end{bmatrix},$

where $(x,\lambda)=z$. Therefore, the partial Jacobians are,

$\partial_{1}\kappa(z,\eta)=\begin{bmatrix}\nabla^{2}c(x)&A(\eta)^{T}\\\
\mathbf{diag}(\lambda)A(\eta)&\mathbf{diag}(A(\eta)x-b(\eta))\end{bmatrix},\quad\partial_{2}\kappa(z,\eta)=\begin{bmatrix}\sum_{n}e_{n}\lambda^{T}\partial
a_{n}(\eta)\\\ \mathbf{diag}(\lambda)\left(-\partial b(\eta)+\sum_{n}\partial
a_{n}(\eta)x_{n}\right)\end{bmatrix},$

where $a_{n}(\eta)\in\mathbf{R}^{M}$ is the $n$th column of $A(\eta)$ and
$e_{n}\in\mathbf{R}^{N}$ is the $n$th standard unit vector. By Assumptions 2
and 6, each block of the two matrices is Lipschitz continuous. It follows that
the matrices themselves are Lipschitz continuous (via the triangle
inequality). ∎

######  Lemma 10.

For any $\eta\in\mathcal{H}$, the matrix
$\partial_{1}\kappa(z^{*}(\eta),\eta)$ is invertible.

###### Proof.

Throughout the proof, we will drop the argument of functions of $\eta$ for
readability. Strong convexity implies that the Hessian $H=\nabla^{2}c(x)$ is
invertible. Therefore, by the Schur’s formula [33], $\partial_{1}\kappa$ is
invertible if and only if

$S=\mathbf{diag}(Ax-b)-\mathbf{diag}(\lambda)AH^{-1}A^{T}$

is invertible. To show this, we will show that $Sv=0$ implies $v=0$, i.e., $S$
has linearly independent columns. So suppose $Sv=0$. Denote
$\mathcal{B}\subset\\{1,\ldots,M\\}$ the set of indices of binding
constraints, i.e., the indices $m\in\mathcal{B}$ such that
$\tilde{a}_{m}^{T}x-b_{m}=0$.

First, note that for $m\not\in\mathcal{B}$, complementary slackness implies
that $\lambda_{m}=0$. Consequently,
$(\mathbf{diag}(\lambda)AH^{-1}A^{T}v)_{m}=0$. Therefore, since $Sv=0$, it
must be that $(\mathbf{diag}(Ax-b)v)_{m}=(\tilde{a}_{m}^{T}x-b_{m})v_{m}=0$ as
well. But since $m\not\in\mathcal{B}$, we know $\tilde{a}_{m}^{T}x-b_{m}<0$,
so it must be that $v_{m}=0$ instead. As a result, we know that $v_{m}=0$ for
all $m\not\in\mathcal{B}$.

Now consider $m\in\mathcal{B}$. By definition, $\tilde{a}_{m}^{T}x-b_{m}=0$,
so $(\mathbf{diag}(Ax-b)v)_{m}=0$. Therefore, it must be that

$(\mathbf{diag}(\lambda)AH^{-1}A^{T}v)_{m}=\lambda_{m}(AH^{-1}A^{T}v)_{m}=0$

as well. By strict complementarity (Assumption 4), $\lambda_{m}>0$, so it must
be that $(AH^{-1}A^{T}v)_{m}=0$. Now denote
$A_{\mathcal{B}}\in\mathbf{R}^{|\mathcal{B}|\times N}$ the submatrix of $A$
that includes each of the rows in $\mathcal{B}$ and similarly denote
$v_{\mathcal{B}}\in\mathbf{R}^{|\mathcal{B}|}$ the entries of $v$ in
$\mathcal{B}$. Then we can rewrite $(AH^{-1}A^{T}v)_{m}=0$ for
$m\in\mathcal{B}$ as

$A_{\mathcal{B}}H^{-1}A_{\mathcal{B}}^{T}v_{\mathcal{B}}=0.$

Of course, this implies
$v_{\mathcal{B}}A_{\mathcal{B}}H^{-1}A_{\mathcal{B}}^{T}v_{\mathcal{B}}=0$ as
well. Since $H$ is positive-definite, it has invertible square root $H^{1/2}$,
so we can rewrite the above as,

$\|H^{-1/2}A_{\mathcal{B}}^{T}v_{\mathcal{B}}\|_{2}^{2}=0,$

which then implies $H^{-1/2}A_{\mathcal{B}}^{T}v_{\mathcal{B}}=0$. Multiplying
both sides by $H^{1/2}$, we conclude $A_{\mathcal{B}}^{T}v_{\mathcal{B}}=0$.
By Assumption 3, this can only hold if $v_{\mathcal{B}}=0$.

Since $v_{m}=0$ for $m\not\in\mathcal{B}$ and $m\in\mathcal{B}$, it must be
that $v=0$. We conclude that $S$ has linearly independent columns and is
therefore invertible. As a result, we conclude
$\partial_{1}\kappa(z^{*}(\eta),\eta)$ is invertible as well. ∎

######  Lemma 11.

There exists $\epsilon>0$ such that for all $\eta\in\mathcal{H}$, the Jacobian
$\partial_{1}\kappa(z^{*}(\eta),\eta)$ has minimum singular value

$\rho(\eta)=\sigma_{\min}(\partial_{1}\kappa(z^{*}(\eta),\eta))\geq\epsilon.$

###### Proof.

We start by noting that Lemma 10 implies that the Implicit Function Theorem
(Theorem 1) holds for every $\eta\in\mathcal{H}$. So $z^{*}(\eta)$ is
differentiable and, hence, continuous, for all $\eta\in\mathcal{H}$. As a
consequence of this and Lemma 9, the Jacobian
$\partial_{1}\kappa(z^{*}(\eta),\eta)$ is a continuous function of $\eta$ as
well.

Since the minimum singular value is a continuous function (by Weyl’s
inequality), the function $\rho$ is continuous. Since $\mathcal{H}$ is also
compact, the Extremal Value Theorem implies that there exists
$\hat{\eta}\in\mathcal{H}$ such that

$\rho(\hat{\eta})=\inf_{\eta\in\mathcal{H}}\rho(\eta).$

Of course, $\hat{\eta}\in\mathcal{H}$, so $\rho(\hat{\eta})=\epsilon>0$
(otherwise, $\partial_{1}\kappa(z^{*}(\eta),\eta)$ is not invertible). So
$\inf_{\eta\in\mathcal{H}}\rho(\eta)=\epsilon$, which proves the claim. ∎

######  Corollary 12.

The function $M(\eta)=(\partial_{1}\kappa(z^{*}(\eta),\eta))^{-1}$ is
Lipschitz continuous and bounded.

###### Proof.

Boundedness and, hence, compactness, follows directly from Lemma 11, since

$\displaystyle\left\|(\partial_{1}\kappa(z^{*}(\eta),\eta))^{-1}\right\|=\frac{1}{\sigma_{\max}\left((\partial_{1}\kappa(z^{*}(\eta),\eta))^{-1}\right)}=\frac{1}{\sigma_{\min}\left(\partial_{1}\kappa(z^{*}(\eta),\eta)\right)}\leq
1/\epsilon.$

Continuity follows from the fact that the inverse map $A\rightarrow A^{-1}$ is
continuous [62]. So $M(\eta)$ is a continuous function on a compact set, and
is hence Lipschitz continuous as well. ∎

######  Lemma 13.

The Jacobian $\partial z^{*}(\eta)$ and the gradient of the implicit planning
objective $\nabla J(\eta)$ are Lipschitz continuous functions over
$\eta\in\mathcal{H}$.

###### Proof.

Both functions are compositions, sums, and products of bounded, Lipschitz
continuous functions on compact sets. Therefore, they are Lipschitz
continuous.

For clarity, we will write this out explicitly for the Jacobian of the
solution map $\partial z^{*}(\eta)$. Choose $\eta_{1},\eta_{2}\in\mathcal{H}$.
For ease of notation, let $\partial_{i}\kappa_{j}$ denote
$\partial_{i}\kappa(z^{*}(\eta_{j}),\eta_{j})$ for $i,j\in\\{1,2\\}$. Then

$\displaystyle\|\partial z^{*}(\eta_{1})-\partial z^{*}(\eta_{2})\|$
$\displaystyle=\left\|\left(\partial_{1}\kappa_{1}\right)^{-1}\partial_{2}\kappa_{1}-\left(\partial_{1}\kappa_{2}\right)^{-1}\partial_{2}\kappa_{2}\right\|$
$\displaystyle\leq\left\|\left(\left(\partial_{1}\kappa_{1}\right)^{-1}-\left(\partial_{1}\kappa_{2}\right)^{-1}\right)\partial_{2}\kappa_{1}\right\|+\left\|\left(\partial_{1}\kappa_{2}\right)^{-1}\left(\partial_{2}\kappa_{1}-\partial_{2}\kappa_{2}\right)\right\|$
$\displaystyle\leq\left\|\left(\partial_{1}\kappa_{1}\right)^{-1}-\left(\partial_{1}\kappa_{2}\right)^{-1}\right\|\left\|\partial_{2}\kappa_{1}\right\|+\left\|\left(\partial_{1}\kappa_{2}\right)^{-1}\right\|\left\|\partial_{2}\kappa_{1}-\partial_{2}\kappa_{2}\right\|.$

Of course, both $\partial_{2}K$ and $(\partial_{1}K)^{-1}$ are both Lipschitz
continuous and bounded by Lemma 9 and Corollary 12, respectively. So let $L$
and $M$ be a Lipschitz constant and an upper bound for both matrix-valued
functions. Then the expression above is bounded by

$\displaystyle\leq
M\left\|\left(\partial_{1}\kappa_{1}\right)^{-1}-\left(\partial_{1}\kappa_{2}\right)^{-1}\right\|+M\left\|\partial_{2}\kappa_{1}-\partial_{2}\kappa_{2}\right\|$
$\displaystyle\leq 2ML\|\eta_{1}-\eta_{2}\|,$

which shows $\partial z^{*}(\eta)$ is Lipschitz continuous with Lipschitz
constant $2ML$. A nearly identical proof can be used to show the gradient
$\nabla J(\eta)=\gamma+\partial z^{*}(\eta)^{T}\cdot\nabla h(z^{*}(\eta))$ is
Lipschitz continuous as well. ∎

We are now ready to prove the main result.

###### Proof of Theorem 8..

Our proof is essentially the same of that in [28], except that we ignore error
in the evaluation of the gradient. Let $L$ be the Lipschitz constant of
$\nabla J(\eta)$ and let $J^{*}=\inf_{\eta\in\mathcal{H}}J(\eta)$ (which is
finite by continuity and compactness). To show convergence, we will
recursively apply the standard quadratic upper bound for functions with
Lipschitz continuous gradients,

$J(v)\leq J(u)+\nabla J(u)^{T}(v-u)+(L/2)\|u-v\|_{2}^{2},$

where $u,v\in\mathbf{R}^{K}$, which can be found in [61], for example. We also
use the fact that projections onto convex sets are nonexpansive [63],

$\|\mathbf{proj}_{\mathcal{H}}(u)-\mathbf{proj}_{\mathcal{H}}(v)\|_{2}^{2}\leq(u-v)^{T}(\mathbf{proj}_{\mathcal{H}}(u)-\mathbf{proj}_{\mathcal{H}}(v)),$

for $u,v\in\mathbf{R}^{K}$. In particular, for $u=\eta^{i}-\alpha\nabla
J(\eta^{i})$ and $v=\eta^{i}$, we obtain

$\|\eta^{i+1}-\eta^{i}\|_{2}^{2}\leq-\alpha\nabla
J(\eta^{i})^{T}(\eta^{i+1}-\eta^{i}),$ (19)

Then we obtain the following inequality,

$\displaystyle J(\eta^{i+1})$ $\displaystyle\leq J(\eta^{i})+\nabla
J(\eta^{i})^{T}(\eta^{i+1}-\eta^{i})+(L/2)\|\eta^{i+1}-\eta^{i}\|_{2}^{2}$
$\displaystyle\leq
J(\eta^{i})+\left(\frac{L}{2}-\frac{1}{\alpha}\right)\|\eta^{i+1}-\eta^{i}\|_{2}^{2}$
by (19) $\displaystyle\leq J(\eta^{i})-(L/2)\|\eta^{i+1}-\eta^{i}\|_{2}^{2}$
$\displaystyle\textrm{since $\alpha\leq 1/L$}.$

Applying this bound recursively, we obtain

$J(\eta^{i+1})\leq
J(\eta^{1})-\frac{L}{2}\sum_{j=1}^{i}\|\eta^{j+1}-\eta^{j}\|_{2}^{2}.$

Since $J$ is a continuous function being minimized over a compact set, it has
a finite lower bound $J^{*}\leq J(\eta)$ for all $\eta\in\mathcal{H}$. This
gives a finite upper bound on the series of square norms,

$\sum_{j=1}^{i}\|\eta^{j+1}-\eta^{j}\|_{2}^{2}\leq\frac{2(J(\eta^{1})-J^{*})}{L}.$

In particular, this implies at least one of the $j$s in the summand must be
such that $\|\eta^{j+1}-\eta^{j}\|_{2}^{2}\leq 2(J(\eta^{1})-J^{*})/Li$, which
proves the theorem. ∎

##  Appendix B Technical results on convex relaxation

###  B.1 Approximate strong duality formulation

We restate Lemma 3 here for convenience.

######  Lemma.

Suppose $h$ is $L$-Lipschitz continuous. For all $\epsilon>0$, there exists
some $\delta>0$ such that the $\delta$-strong duality formulation (P2) is an
$\epsilon$-approximation of (P1) over $\eta$, i.e., for all $(\eta,x,\lambda)$
feasible for (P2),

$|h(x,\lambda)-h(x^{*}(\eta),\lambda^{*}(\eta))|\leq\epsilon.$

To prove Lemma 3, we first show that the feasibility for the $\delta$-strong
duality problem implies that $(x,\lambda)$ is close to
$(x^{*}(\eta),\lambda^{*}(\eta))$. Then, it is straightforward to show main
result.

######  Lemma 14.

For all $\epsilon>0$, there exists some $\delta>0$ such that for all
$\eta\in\mathcal{H}$,

$\displaystyle c(x)+b(\eta)^{T}\lambda+c^{*}(-A(\eta)^{T}\lambda)\leq\delta$
(20) $\displaystyle A(\eta)x\leq b(\eta)$ $\displaystyle\lambda\geq 0$

implies

$\max\left(\|x-x^{*}(\eta)\|_{2},\|\lambda-\lambda^{*}(\eta)\|_{2}\right)\leq\epsilon$

###### Proof.

Fix $\epsilon>0$. Now suppose for sake of contradiction that for all
$\delta>0$, there exists
$(\eta,x,\lambda)\in\mathbf{R}^{N}\times\mathbf{R}^{M}$ such that the
constraints in (20) are satisfied but $\|x-x^{*}(\eta)\|_{2}\geq\epsilon$. We
can therefore build sequence of points
$(\eta_{1},x_{1},\lambda_{1}),(\eta_{2},x_{2},\lambda_{2}),\ldots$ where
$x_{i}$ is bounded away from $x_{i}^{*}(\eta_{i})$ and $A(\eta_{i})x_{i}\leq
b(\eta_{i})$, $\lambda_{i}\geq 0$, and
$c(x_{i})+b(\eta_{i})^{T}\lambda_{i}+c^{*}(-A(\eta_{i})^{T}\lambda_{i})\leq
1/i$. By continuity of $c$ and $c^{\star}$ (which follow from Lipschitz
continuous gradients and strong convexity of $c$, respectively) and the
compactness of the domain of $(x,\lambda)$, there exists a subsequence of this
sequence that has a limiting point $(\hat{\eta},\hat{x},\hat{\lambda})$ such
that

$c(\hat{x})+b(\hat{\eta})^{T}\hat{\lambda}+c^{*}(-A(\hat{\eta})^{T}\hat{\lambda})=0.$

Similarly, by the continuity of $x^{*}(\eta)$ (see Lemma 10), we know
$\|\hat{x}-x^{*}(\hat{\eta})\|\geq\epsilon$. This implies that
$x^{*}(\hat{\eta})$ is not the unique solution to the dispatch model with
parameter $\hat{\eta}$, which is a contradiction. Using an identical line of
reasoning, we can show the analogous result for
$\|\lambda-\lambda^{*}(\eta)\|_{2}$, proving the claim. ∎

###### Proof of Lemma 3.

First, lets show that the bound in Lemma 14 is uniform. Indeed

Fix $\epsilon>0$. By Lemma 14, we can pick $\delta$ such that feasibility in
(P2) implies

$\max\left(\|x-x^{*}(\eta)\|_{2},\|\lambda-\lambda^{*}(\eta)\|_{2}\right)\leq\frac{\epsilon}{2L}$

Now consider some $(\eta,x,\lambda)$ feasible for the $\delta$-strong duality
formulation in (P2). By Lipschitz continuity,

$|h(x,\lambda)-h(x^{*}(\eta),\lambda^{*}(\eta))|\leq
L\left(\|x-x^{*}(\eta)\|_{2}+\|\lambda-\lambda^{*}(\eta)\|_{2}\right)\leq\epsilon,$

which proves the claim. ∎

###  B.2 Tight integer relaxation

We restate the theorem here for convenience.

######  Theorem 15 (Tight integer relaxation).

Consider an integer constrained variable
$\alpha\in\\{\alpha^{\min},\alpha^{\max}\\}$ and a box constrained variable
$\beta\in[\beta^{\min},\beta^{\max}]$ (which may or may not have integer
constraints as well). Then for any $y\in\mathbf{R}$, the equation
$y=\alpha\beta$ holds if and only if

$y\in\mathbf{env}(\alpha\beta).$

As a consequence, if $\mathcal{H}=\mathcal{H}^{\textrm{int}}$, then the
feasible sets of (P2) and its relaxation (P3) are equal, i.e., the relaxation
is tight.

###### Proof.

This result may be found in Theorem 3 of [48], but we prove it here for ease
of access. By the definition of the convex and concave envelope, we already
know that $y=\alpha\beta$ implies $y\in\mathbf{env}(\alpha\beta)$. So it
suffices to show that $y\in\mathbf{env}(\alpha\beta)$ implies $y=\alpha\beta$.
Recall that $y\in\mathbf{env}(\alpha\beta)$ if and only if,

$\displaystyle y$
$\displaystyle\geq\alpha\beta^{\min}+\alpha^{\min}\beta-\alpha^{\min}\beta^{\min},$
$\displaystyle y$
$\displaystyle\geq\alpha^{\max}\beta+\alpha\beta^{\max}-\alpha^{\max}\beta^{\max},$
$\displaystyle y$
$\displaystyle\leq\alpha^{\min}\beta+\alpha\beta^{\max}-\alpha^{\min}\beta^{\max},$
$\displaystyle y$
$\displaystyle\leq\alpha\beta^{\min}+\alpha^{\max}\beta-\alpha^{\max}\beta^{\min}.$

Now first consider the case when $\alpha=\alpha^{\min}$. Then the first and
third equations imply $y\geq\alpha\beta$ and $y\leq\alpha\beta$, respectively.
So $y=\alpha\beta$. If $\alpha=\alpha^{\max}$, the the second and fourth
equations imply $y\geq\alpha\beta$ and $y\leq\alpha\beta$, respectively, so
$y=\alpha\beta$ in this case as well. Since $\alpha$ can only be either
$\alpha^{\min}$ or $\alpha^{\max}$ by assumption, it holds that
$y=\alpha\beta$ if and only if $y\in\mathbf{env}(\alpha\beta)$. ∎

###  B.3 Approximate continuous relaxation

We restate Theorem 6 here for convenience.

######  Theorem 16 (Approximate continuous relaxation).

For all $\epsilon>0$, there exists some $\alpha>0$ such that when

$\|\eta^{\min}-\eta^{\max}\|_{\infty},\leq\alpha,$

then (P3) is a $\epsilon$-approximation of (P2) over $\eta$.

To prove this, we first show the relaxed feasible set approximates the
original feasible set in the limiting sense.

######  Lemma 17.

Fix $z^{\min},z^{\max},\eta^{\min},\eta^{\max}\in\mathbf{R}$. Then for all
$z\in[z^{\min},z^{\max}]$ and $\eta\in[\eta^{\min},\eta^{\max}]$, it holds
that

$\sup_{y\in\mathbf{env}(z\eta)}|y-z\eta|\leq|z^{\max}-z^{\min}||\eta^{max}-\eta^{\min}|.$
(21)

###### Proof.

Fix $z\in[z^{\min},z^{\max}]$ and $\eta\in[\eta^{\min},\eta^{\max}]$. First
note that the function $\psi(y)=|y-z\eta|$ is convex. Since the maximum of a
convex function over a convex set is attained at an extreme point of the set,
it suffices to consider the boundary of $\mathbf{env}(z\eta)$. At one of the
lower convex envelopes we have,

$\displaystyle|\eta
z^{\min}+\eta^{\min}z-\eta^{\min}z^{\min}-z\eta|=|(z-z^{\min})(\eta-\eta^{\min})|$
$\displaystyle\leq|z^{\max}-z^{\min}||\eta^{\max}-\eta^{\min}|.$

The same line of reasoning can be used for the three other inequalities that
define $\mathbf{env}(z\eta)$, proving the claim. ∎

By Assumption (5), the primal and dual dispatch solutions are assumed to lie
in a compact domain. Thus, we can assume they lie in some finite (albeit
possibly large) box. This implies that when
$\|\eta^{\max}-\eta^{\min}\|_{\infty}$ is sufficiently small, the set
$\mathbf{env}(z\eta)$ approaches the point $z\eta$.

Next, we state and prove a general result that we will later show includes
Theorem 6 as a special case.

######  Lemma 18.

Let $f:\mathbf{R}^{D}\times\mathbf{R}^{D}\rightarrow\mathbf{R}$ and
$g:\mathbf{R}^{D}\times\mathbf{R}^{D}\rightarrow\mathbf{R}$ be convex,
$L$-Lipschitz continuous functions and let
$z^{\min},z^{\max},\eta^{\min},\eta^{\max}\in\mathbf{R}^{D}$. Assume that
there exists $\alpha>0$ such that for all $\eta\in[\eta^{\min},\eta^{\max}]$,
there exists $z\in[z^{\min},z^{\max}]$ such that

$g(z,\eta,\mathbf{diag}(\eta)z)<-\alpha.$ (22)

Now consider the optimization problem

$\begin{array}[]{ll}\mathrm{minimize}&f(z,\eta)\\\\[2.5pt] \mathrm{subject\
to}&g(z,\eta,y)\leq 0\\\\[2.5pt] &y=\mathbf{diag}(\eta)z\\\\[2.5pt]
&\eta^{\min}\leq\eta\leq\eta^{\max}\\\\[2.5pt] &z^{\min}\leq z\leq
z^{\max},\end{array}$ (O)

where the variables are $z,\eta,y\in\mathbf{R}^{D}$. Also consider the relaxed
problem

$\begin{array}[]{lll}\mathrm{minimize}&f(z,\eta)\\\\[2.5pt] \mathrm{subject\
to}&g(z,\eta,y)\leq 0\\\\[2.5pt]
&y_{i}\in\mathbf{env}(\eta_{i}z_{i}),&i=1,\ldots,D\\\\[2.5pt]
&\eta^{\min}\leq\eta\leq\eta^{\max}\\\\[2.5pt] &z^{\min}\leq z\leq
z^{\max},\end{array}$ (R)

where the variables are also $z,\eta,y\in\mathbf{R}^{D}$. Then for all
$z^{\min},z^{\max}\in\mathbf{R}^{D}$ and $\epsilon>0$, there exists $\delta>0$
such that,

$\|\eta^{\max}-\eta^{\min}\|_{\infty}\leq\delta$

implies that (R) is an $\epsilon$-approximation of (O) over $\eta$, i.e., for
all $(z,\eta,y)$ feasible for (R), there exists $z^{\prime}\in\mathbf{R}^{D}$
such that $(z^{\prime},\eta,\mathbf{diag}(\eta)z^{\prime})$ is feasible for
(O) and

$|f(z,\eta)-f(z^{\prime},\eta)|\leq\epsilon.$

###### Proof.

Without loss of generality, assume $L\geq 1$, $\|z^{\max}-z^{\min}\|_{2}\geq
1$, and $\epsilon<1$. Fix $z^{\min},z^{\max}\in\mathbf{R}^{D}$ and
$\epsilon>0$ and suppose

$\|\eta^{\max}-\eta^{\min}\|_{\infty}\leq\delta=\frac{\alpha\epsilon}{2L^{2}\sqrt{D}\|z^{\max}-z^{\min}\|_{\infty}\|z^{\max}-z^{\min}\|_{2}}.$

Now consider some point $v=(z,\eta,y)$ feasible for (R). Our proof will
proceeds in three steps.

#### Step 1.

By (22), there exists $\tilde{z}\in[z^{\min},z^{\max}]$ such that
$g(\tilde{z},\eta,\mathbf{diag}(\eta)\tilde{z})\leq-\alpha<0$. Let
$\tilde{v}=(\tilde{z},\eta,\mathbf{diag}(\eta)\tilde{z})$. The point
$\tilde{v}$ is feasible for (O), so it must be feasible for (R) as well. Now
choose $\hat{v}=(1-\theta)v+\theta\tilde{v}$. where

$\theta=\frac{\epsilon}{L\|z^{\max}-z^{\min}\|_{2}},$

and note that $\hat{\eta}=\eta$. Since (R) is convex, $\hat{v}$ is feasible
for (R) and

$g(\hat{v})\leq(1-\theta)g(v)+\theta g(\tilde{v})\leq-\theta\alpha.$

#### Step 2.

Let $v^{*}=(\hat{z},\eta,\mathbf{diag}(\eta)\hat{z})$. We claim $v^{*}$ is
feasible for (O). Of course, this point satisfies the second, third, and
fourth constraints of (O) by construction. But since
$\|\eta^{\max}-\eta^{\min}\|_{\infty}\leq\delta$, then by Lemma 17 it must be
that

$\|\hat{v}-v^{*}\|_{2}=\|\hat{y}-\mathbf{diag}(\eta)\hat{z}\|_{2}\leq\sqrt{D}\|z^{\max}-z^{\min}\|_{\infty}\delta.$

So, since $g$ is $L$-Lipschitz, we have

$\displaystyle g(v^{*})$ $\displaystyle\leq g(\hat{v})+L\|\hat{v}-v^{*}\|_{2}$
$\displaystyle\leq-\theta\alpha+L\sqrt{D}\|z^{\max}-z^{\min}\|_{\infty}\delta$
$\displaystyle\leq
L\sqrt{D}\|z^{\max}-z^{\min}\|_{\infty}\delta-\frac{\alpha\epsilon}{L\|z^{\max}-z^{\min}\|_{2}}$
$\displaystyle\leq-\frac{\alpha\epsilon}{2L\|z^{\max}-z^{\min}\|_{2}}\leq 0,$

which means the second constraint is satisfied, and hence $v^{*}$ is feasible
for (O).

#### Step 3.

Lastly, we claim that $|f(z,\eta)-f(z^{*},\eta^{*})|\leq\epsilon$. Recall that
$\eta^{*}=\eta$ and $z^{*}=\hat{z}=(1-\theta)z+\theta\tilde{z}$. Since $f$ is
$L$-Lipschitz, we have

$\displaystyle|f(z,\eta)-f(z^{*},\eta^{*})|$ $\displaystyle\leq
L\|z-\hat{z}\|_{2}$ $\displaystyle=L\theta\|z-\tilde{z}\|_{2}$
$\displaystyle\leq L\theta\|z^{\max}-z^{\min}\|_{2}\leq\epsilon,$

which proves the claim. ∎

Finally, we can prove Theorem 6 by showing it satisfies the conditions of
Lemma (18).

###### Proof of Theorem 6.

This is simply a special case of Lemma 18 with variable $z=(x,\lambda)$,
objective $f(z,\eta)=\gamma^{T}\eta+h(x,\lambda)$, and constraint

$g(z,\eta,y)=\begin{bmatrix}c(x)-\mu-c^{\star}(-V\mathbf{1})+\delta\\\
Y\mathbf{1}-b(\eta)\\\ \end{bmatrix},$

(since each of $Y$, $V$, and $\mu$ are linear functions of
$y=\mathbf{diag}(\eta)z$). The only other step is to show that that there
exists $\alpha>0$ such that $g(z,\eta,y)<-\alpha$. Of course, we can simply
set $\alpha<\delta$ to show this for the first term in $g$. For the second
term, this follows directly from Assumption 1, compactness assumptions on
$\eta$ and $x$ (Assumption 5). ∎

##  Appendix C Experiment Details

We use the PyPSA-USA [58] dataset for all our experiments. The software is
freely available on Github at https://github.com/pypsa/pypsa-usa. We
specifically use commit `065f864` from December 8th, 2023. Using that version
of the package should exactly reproduce the data as in the paper. We generate
data for the U.S. Western Interconnect for the full year of 2019 and use the
default settings to cluster the network to the desired size. The dataset
includes fuel costs, capital, renewable generator maximum expansion
capacities, a realistic electrical network with susceptances and capacities,
battery locations and limits, and other data.

We scale all loads by 150% in order to simulate future load growth. Then, we
decrease generation and transmission capacities by 20% to model retirements
and outages. Existing transmission lines and batteries are allowed to be
expanded by a factor of 10; expansion costs are given by the dataset itself.
Load is curtailable with a curtailment cost of $500 per MWh. Most solutions
have no load curtailment. Renewable expansion capacities are set to be 0.01%
of the PyPSA provided capacity limit. This is because PyPSA uses technical
potentials, i.e., the physical capacity of the land, ignoring zoning and
existing construction. All generators are given a maximum potential capacity
of 50 GW for numerical stability reasons. In practice, this limit is never
encountered. Finally, due to a bug in the PyPSA-USA code, we set battery
capacity expansion costs manually. The battery costs, from the Danish Energy
Agency [64], are $151,940/MWh of battery capacity, annualized over 25 years,
and $171,200/MW of maximum power output, annualized over 10 years. The
discount rate is set to 7.0% and an annual fixed operations and maintenance
cost of $52,845/MW is also added. Batteries with less than 1 hour of duration
are ignored, since these are likely used for ancillary services in practice.
When solving cases for $T$ hours, all capital costs are scaled to $T/8760$
times their annualized cost in order to properly scale capital and operation
costs.