# Conditional Value at Risk-Sensitive Solar Hosting Capacity Analysis in
Distribution Networks

Avinash N. Madavan Nathan Dahlin Subhonmesh Bose Lang Tong A. N. Madavan, N.
Dahlin, and S. Bose are with the Department of Electrical and Computer
Engineering at the University of Illinois at Urbana-Champaign, Urbana IL, USA
61801. L. Tong is with the School of Electrical and Computer Engineering at
Cornell University, Ithaca NY, USA 14680. This work was partially supported by
the NSF CAREER grant 2048065. Emails: {madavan2, dahlin<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Solar hosting capacity analysis (HCA) assesses the ability of a distribution
network to host distributed solar generation without seriously violating
distribution network constraints. In this paper, we consider risk-sensitive
HCA that limits the risk of network constraint violations with a collection of
scenarios of solar irradiance and nodal power demands, where risk is modeled
via the conditional value at risk (CVaR) measure. First, we consider the
question of maximizing aggregate installed solar capacities, subject to risk
constraints and solve it as a second-order cone program (SOCP) with a standard
conic relaxation of the feasible set with power flow equations. Second, we
design an incremental algorithm to decide whether a configuration of solar
installations has acceptable risk of constraint violations, modeled via CVaR.
The algorithm circumvents explicit risk computation by incrementally
constructing inner and outer polyhedral approximations of the set of
acceptable solar installation configurations from prior such tests conducted.
Our numerical examples study the impact of risk parameters, the number of
scenarios and the scalability of our framework.

###### Index Terms:

Hosting capacity analysis, distributed solar, conditional value at risk

##  I Introduction

According to a recent study in [1], usable rooftop space in the United States
can produce enough solar energy to meet the nation’s current electricity
demands. In 2022, new utility-scale solar is projected to account for almost
half of all planned capacity expansions, per the Energy Information
Administration. Although small-scale distributed rooftop solar (with
capacities less than 1MW) only produces about 1% of total electricity
generation today, this sector is rapidly growing, driven in part by ambitious
targets set forth by renewable portfolio standards (e.g., in [2]). California
has exceeded 10GW of installed capacity of distributed solar in 2022.
Distributed solar accounts for more than $60$% of solar projects in New York,
Hawaii, Maryland, New Jersey and Massachusetts. These numbers are expected to
grow, aided by legislation such as Title 24 in California and bill S2165 in
Massachusetts that require new building constructions to install rooftop
solar. Given this rapid growth, we ask: _is the grid ready to host distributed
solar at scale_?

Solar production curtailments, exceeding 15% of total potential output, have
often been necessary in California to avoid violation of system constraints,
including voltage and harmonic ranges, and protection system failures, etc.
See [3, 4] for details. Today, these interventions impact utility-scale solar
installations. However, similar questions arise for distributed solar in
low/medium-voltage distribution networks. Specifically, how much installed
solar capacity can the distribution network host, without jeopardizing safe
operational limits for the network? Such quantification studies are
categorized as _hosting capacity analyses_ (HCA). This paper aims to provide
efficient mathematical frameworks for the same.

Traditional approaches to HCA have relied on historic peak load levels or
point estimates of solar generation and network load profiles, e.g., see [5,
6]. Distributed solar is highly intermittent (especially during the spring
season), and network load characteristics are expected to change significantly
with the increased adoption of new technologies such as home batteries and
electric vehicles. Consequently, deterministic HCA, conducted with a few
representative scenarios, will fail to capture the variety of operating
conditions the network might face. See [5, 7] for a discussion. Instead, one
must account for the _statistics_ of operating conditions, both observed in
historical data and those obtained with reasonable projections of load
characteristics into the future. Several recent works have indeed proposed to
consider said statistics in HCA through stochastic programming; see [8] for a
survey.

The central challenge of stochastic optimization-based HCA is _scalability_.
Such techniques must consider large numbers of possible grid scenarios, drawn
from sampled or synthesized data. Using parametric families of distributions
to describe such data is often difficult, especially if synthetic data is
obtained from generative models. HCA methods that use scenario-based
optimization typically entail the solution of a large number of power flow
problems to compute network-relevant variables such as voltages and power
flows in each scenario. Then, the statistics of possible network constraint
violations with these voltages/power flows are processed, possibly through a
risk measure, e.g., the probability of constraint violations in [8]. See
Figure 1 for an illustration of the work-flow of stochastic HCA methods. The
key innovation of our work lies in the use of the _conditional value at risk_
(CVaR) measure to process the statistics of constraint violation, as we
explain in Section II. While the use of this risk measure is popular in
finance (e.g., see [9, 10]), its favorable mathematical properties are
popularizing its use in engineering domains, including power systems. As will
become evident, these properties allow efficient algorithm design to tackle
the scalability challenge of HCA.

The specific contributions of this paper are as follows. First in Section III,
we propose a CVaR-sensitive hosting capacity maximization problem that seeks
to identify the maximum aggregate capacity of solar installations in a
distribution network, and describe a scenario-based approach to solve it in
Section IV via _convex optimization_. This convexity follows from our use of
CVaR, as well as a standard second-order cone programming-based relaxation of
the DistFlow power flow model. The vast majority of risk measures introduce
nonconvexity. Moreover, such measures require explicit evaluation of the risk
of constraint violations across scenarios in order to assess the risk of a
particular installation configuration. As a result, one must consider solar
installations one at a time, assess the risk level for each, and update the
configuration suitably to optimize it. In contrast, CVaR preserves the
convexity of our power flow model and allows for formulation of an
optimization problem that directly optimizes over realizable installations
according to risk level, defined through scenarios, and quality as captured by
the chosen objective function. Convexity allows us to leverage mature off-the-
shelf convex optimization tools and avoids the pitfalls of constraint
linearization, mixed integer programming, or the design of complex, parameter-
sensitive algorithms to explore the space of capacity configurations that
require careful parameter-tuning, e.g., in [11, 12, 13]. Our numerical
experiments demonstrate the impact of risk parameters and the number of
scenarios on the solution quality and runtime.

Figure 1: The work-flow of risk-sensitive solar HCA.

Second, we propose a method in Section V to quickly determine whether a
candidate solar installation configuration has acceptable risk of network
constraint violations, as evaluated via CVaR. Our proposed method is
incremental, i.e., we utilize the results of prior tests to expedite
determination as to whether a new candidate capacity configuration is
acceptable. The properties of the CVaR measure dictate that the resulting set
of acceptable configurations is convex. Thus, we construct and refine convex
inner and outer approximations to this set as more configurations are tested.
These approximations often permit us to quickly certify whether a new
candidate configuration is acceptable, without having to explicitly test for
it using all scenarios. Our numerical examples illustrate that acceptability
certification time decreases rapidly, as the knowledge builds over time from
prior tests. We end the paper in Section VI with concluding remarks.

##  II CVaR-Sensitive HCA Problems

Solar HCA assesses whether a distribution network can support distributed
solar generation from installed capacity configurations without seriously
violating distribution network constraints. In this section, we present
abstract formulations of the CVaR-sensitive solar HCA problems that we aim to
study in the rest of the paper, utilizing Figure 1.

Consider a set of solar irradiance and load scenarios $\Omega$, i.e., $\Omega$
captures all uncertainty that the network faces from solar productions and
consumer power demands. For a scenario $\omega\in\Omega$, its solar irradiance
${\bm{\alpha}}(\omega)$ translates into solar power generation across the
network via a configuration of installed solar capacities ${\bm{\psi}}$. For
example, each entry of ${\bm{\psi}}$ might specify the surface area of rooftop
photovoltaic panels at the network buses. Accompanying the conditions
affecting solar generation, scenarios also include power demands at each bus.
Thus, scenarios in $\Omega$ and solar capacity configuration ${\bm{\psi}}$
yield a set of corresponding net nodal power injections ${\cal
S}_{\textrm{inj}}$, identified as solar generation less power demands at each
bus. These injection scenarios then serve as inputs to a power flow model,
which identifies a set ${\cal G}$ of voltages and power flows throughout the
network that are consistent with the real/reactive power injections in ${\cal
S}_{\textrm{inj}}$.

For deterministic hosting capacity analysis, $\Omega$ is a singleton. In such
an analysis, an installed capacity configuration ${\bm{\psi}}$ is deemed
acceptable, if there exist voltages and power flows over the distribution
network, given the power flow model, that satisfy a collection of network
constraints, e.g., thermal limits, voltage limits, etc. For the risk-sensitive
counterpart, we do not certify the acceptability of ${\bm{\psi}}$, based on a
single scenario ${\bm{\omega}}\in\Omega$, but do so _across_ $\Omega$. Thus in
the bottom portion of Figure 1, we consider the _statistics_ of the voltages
and power flows arising from the scenarios in $\Omega$ under an installation
configuration ${\bm{\psi}}$. We utilize a risk-measure $\rho$ to process the
statistics of network constraint violations, and call ${\bm{\psi}}$ feasible,
if the “risk” of network constraint violation is low enough. In this
framework, acceptability is determined in a statistical sense, and we consider
two questions:

  * •

$\mathscr{A}({\bm{\psi}})$: Decide whether an installation configuration
${\bm{\psi}}$ is acceptable.

  * •

$\mathscr{O}$: Optimize a suitable metric of installed capacities
${\bm{\psi}}$ over acceptable configurations.

Figure 2: Mean and the conditional value at risk (CVaR) at level $\delta$ for
a random variable $\zeta$.

In this work, we use the _conditional value at risk_ (CVaR) measure for
$\rho$. For a scalar random variable $\zeta$, we have

${\rm
CVaR}{}_{\delta}[\zeta]:=\underset{t\in\mathbb{R}}{\text{minimum}}\left\\{t+\frac{1}{1-\delta}\mathds{E}[[\zeta-t]_{+}]\right\\},$
(1)

given a parameter $\delta\in[0,1)$. Here, $\mathds{E}[\cdot]$ computes the
expectation with respect to the probability distribution over $\zeta$. We use
the notation $[z]_{+}=\max\\{0,z\\}$ for any scalar $z$. Driving
$\delta\downarrow 0$, $\text{{\rm CVaR}{}}_{\delta}[\zeta]$ approaches the
mean value of $\zeta$. For $\delta\uparrow 1$, $\text{{\rm
CVaR}{}}_{\delta}[\zeta]$ becomes the highest value (essential supremum) that
$\zeta$ can take. For a random variable $\zeta$ that admits a smooth
cumulative distribution function $F_{\zeta}$, as in Figure 2, CVaR gives the
expectation over the $(1-\delta)$-tail of the distribution, i.e.,

${\rm CVaR}{}_{\delta}[\zeta]:=\mathbb{E}[\zeta\,|\,\zeta\geq
F_{\zeta}^{-1}(\delta)],$ (2)

The definition in (1) applies to random variables with arbitrary
distributions. See [14] for details. Thus, given an installation configuration
${\bm{\psi}}$, $\mathscr{A}({\bm{\psi}})$ seeks to decide, if

$\displaystyle{\rm CVaR}_{\delta}[{g}_{i}({\bm{\psi}};\Omega)]\leq 0,\quad
i=1,\ldots,I,$ (3)

where $g_{i}$ encodes the network constraints on voltages and power flows that
support the power injections, given solar capacities ${\bm{\psi}}$ and a
scenario $\omega\in\Omega$. Out of all acceptable configurations,
$\mathscr{O}$ identifies those which optimize a particular objective:

$\displaystyle\text{maximize}\ \text{aggregate capacity }f({\bm{\psi}}),\
\text{subject to }\eqref{eq:P1.abstract}.$ (4)

One can show that, as Figure 2 suggests,

$\displaystyle{\rm CVaR}{}_{\delta}[g_{i}({\bm{\psi}};\Omega)]\leq
0\implies\mathbb{P}\\{g_{i}({\bm{\psi}};\Omega)>0\\}\leq 1-\delta.$ (5)

In other words, the CVaR-sensitive constraint in (3) automatically limits the
violation probability of network constraints. Additionally, it also tightly
controls the _extent_ of such violations. CVaR has emerged as a popular tool
for risk management in finance and engineering, e.g., see [10, 9]. See [15,
16] for its recent use in power system operation. As will become evident,
several properties of ${\rm CVaR}$ make it particularly suitable to our
problem setting.

In Section III, we define the distribution network model and present concrete
formulations of $\mathscr{A}({\bm{\psi}})$ and $\mathscr{O}$ as optimization
problems. These problems are nonconvex, owing to Kirchhoff’s laws. Then, we
_convexify_ them to $\mathscr{A}^{\textrm{C}}({\bm{\psi}})$ and
$\mathscr{O}^{\textrm{C}}$, by considering a standard second-order cone
programming-based relaxation of the feasible set that is defined by the power
flow equations. In Section IV, we develop a scenario-based counterpart
$\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$ of $\mathscr{O}^{\textrm{C}}$ by
exploiting the properties of ${\rm CVaR}$. In Section V, we propose an
incremental algorithm to solve the scenario-based variant
$\mathscr{A}^{\textrm{C}}_{\textrm{sample}}({\bm{\psi}})$ of
$\mathscr{A}^{\textrm{C}}({\bm{\psi}})$ to decide acceptability of
${\bm{\psi}}$.

##  III Formulating $\mathscr{A},\mathscr{O}$ and their Convex Relaxations
$\mathscr{A}^{\textrm{C}},\mathscr{O}^{\textrm{C}}$

We begin by encoding the single-phase equivalent of a three-phase radial
distribution network as a directed graph ${\mathfrak{G}}({\cal V},\vec{{\cal
E}})$. Let ${\cal V}=\\{1,\ldots,n\\}$ describe the distribution buses, where
the first bus is the substation. We say $j\to j^{\prime}\in\vec{{\cal E}}$, if
bus $j$ is connected to bus $j^{\prime}$ via a distribution line. Assign the
directions of these edges arbitrarily. Given that the radial network
${\mathfrak{G}}$ has $n-1$ edges, the edge-to-node adjacency matrix
${\bm{B}}\in\mathbb{R}^{(n-1)\times n}$ is described by

$\displaystyle{\bm{B}}_{ej}=\begin{cases}+1,&\text{if }e=j\to
j^{\prime}\in\vec{{\cal E}}\text{ for some }j^{\prime}\in{\cal V},\\\
-1,&\text{if }e=j^{\prime}\to j\in\vec{{\cal E}}\text{ for some
}j^{\prime}\in{\cal V},\\\ 0,&\text{otherwise}.\end{cases}$ (6)

Consider a probability space $(\Omega,{\cal B},\mathbb{P})$, where $\Omega$
encodes all possible scenarios of solar irradiance and real/reactive power
demand profiles, ${\cal B}$ describes a suitable $\sigma$-algebra over
$\Omega$ and $\mathbb{P}$ describes the probability measure over $\Omega$.
Precisely, define

$\displaystyle{\bm{\alpha}}(\omega)\in\mathbb{R}^{n-1},\
{\bm{p}}_{\textrm{D}}(\omega)\in\mathbb{R}^{n-1},\
{\bm{q}}_{\textrm{D}}(\omega)\in\mathbb{R}^{n-1}$ (7)

as the profiles of solar irradiance, real power demand and reactive power
demand across buses $\overline{{{\cal V}}}:=\\{2,\ldots,n\\}$ in scenario
$\omega\in\Omega$, i.e.,
${\bm{\alpha}},{\bm{p}}_{\textrm{D}},{\bm{q}}_{\textrm{D}}$ are ${\cal
B}$-measurable maps (random variables). If
${\bm{\psi}}\in\mathbb{R}^{n-1}_{+}$ describes the installed capacities of PV
panels across buses in $\overline{{{\cal V}}}$, then the real and reactive
power generation from solar panels across the same buses in scenario $\omega$
is given by

$\displaystyle{\bm{p}}_{\textrm{G}}(\omega)={\bm{\alpha}}(\omega)\odot{\bm{\psi}},\quad{\bm{q}}_{\textrm{G}}(\omega)={\bm{\eta}}_{\textrm{G}}\odot{\bm{\alpha}}(\omega)\odot{\bm{\psi}},$
(8)

i.e., ${\bm{p}}_{\textrm{G}},{\bm{q}}_{\textrm{G}}$ are ${\cal B}$-measurable
for each ${\bm{\psi}}$. We use the notation ${\bm{z}}_{1}\odot{\bm{z}}_{2}$ to
denote the Hadamard (element-wise) product of two arbitrary matrices/vectors
${\bm{z}}_{1}$ and ${\bm{z}}_{2}$. Here,
${\bm{\eta}}_{\textrm{G}}\in\mathbb{R}^{n-1}$ describes a vector of power
factors for solar power injection. We remark that today’s inverters often
operate at unit power factors, for which ${\bm{\eta}}_{\textrm{G}}=0$. The
generation profile assumes a direct scaling of the irradiance with the
installed capacity–a premise that holds when inverters connecting solar panels
track the maximum power point [17]. Then, the net real and reactive power
injections into buses $\overline{{{\cal V}}}$ are, respectively,

$\displaystyle{\bm{p}}_{\textrm{inj}}(\omega)={\bm{p}}_{\textrm{G}}(\omega)-{\bm{p}}_{\textrm{D}}(\omega),\
{\bm{q}}_{\textrm{inj}}(\omega)={\bm{q}}_{\textrm{G}}(\omega)-{\bm{q}}_{\textrm{D}}(\omega),$
(9)

where ${\bm{p}}_{\textrm{inj}},{\bm{q}}_{\textrm{inj}}$ are ${\cal
B}$-measurable.

Derived from Kirchhoff’s laws, the power flow equations identify voltages and
real/reactive power flows over distribution lines that can sustain these power
injections across buses in $\overline{{{\cal V}}}$ in scenario $\omega$. To
describe these equations, we introduce additional notation. For the $n-1$
(directed) distribution lines in $\vec{{\cal E}}$, let
${\bm{R}}\in\mathbb{R}^{n-1}$ and ${\bm{X}}\in\mathbb{R}^{n-1}$ collect their
resistances and reactances. Define ${\bm{P}}(\omega)\in\mathbb{R}^{n-1}$ and
${\bm{Q}}(\omega)\in\mathbb{R}^{n-1}$ as the _sending-end_ real and reactive
power flows in scenario $\omega$ over said lines. Also, let
${\bm{L}}(\omega)\in\mathbb{R}^{n-1}$ describe the squared current magnitude
over these lines. For the $n$ buses in the distribution network, let
${\bm{W}}(\omega)\in\mathbb{R}^{n}$ collect the _squared_ voltage magnitudes
in scenario $\omega$.

The following matrix operations will prove useful in presenting the power flow
equations. For any matrix/vector ${\bm{z}}$, let ${\bm{z}}_{+}$ compute the
positive part of each entry of ${\bm{z}}$. Similarly, define
${\bm{z}}_{-}:=[-{\bm{z}}]_{+}$. Where unambiguous, we use
${\bm{z}}^{2}:={\bm{z}}\odot{\bm{z}}$. Finally, for a vector ${\bm{z}}$ with
more than one row, let ${\bm{\pi}}[{\bm{z}}]$ return the same vector with its
first row removed.

With this notation, we write the power flow equations as

$\displaystyle{\bm{p}}_{\textrm{inj}}(\omega)={\bm{\pi}}\left[{\bm{B}}_{+}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}{\bm{P}}(\omega)-{\bm{B}}_{-}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}\left({\bm{P}}(\omega)-{\bm{R}}\odot{\bm{L}}(\omega)\right)\right],$
(10a)
$\displaystyle{\bm{q}}_{\textrm{inj}}(\omega)={\bm{\pi}}\left[{\bm{B}}_{+}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}{\bm{Q}}(\omega)-{\bm{B}}_{-}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}\left({\bm{Q}}(\omega)-{\bm{X}}\odot{\bm{L}}(\omega)\right)\right],$
(10b)
$\displaystyle{\bm{B}}{\bm{W}}(\omega)=2\left[{\bm{R}}\odot{\bm{P}}(\omega)+{\bm{X}}\odot{\bm{Q}}(\omega)\right]$
$\displaystyle\qquad\qquad-({\bm{R}}^{2}+{\bm{X}}^{2})\odot{\bm{L}}(\omega),$
(10c)
$\displaystyle\left[{\bm{B}}_{+}{\bm{W}}(\omega)\right]\odot{\bm{L}}(\omega)=[{\bm{P}}(\omega)]^{2}+[{\bm{Q}}(\omega)]^{2},$
(10d)

This form of the power flow equations for radial distribution networks has
been derived in [18]. The first two relations encode power balance at each bus
in $\overline{{{\cal V}}}$. The last two equations describe Kirchhoff’s
voltage law, expressed in terms of squared line currents and real/reactive
power flows.

For a given installed capacity configuration ${\bm{\psi}}$, a scenario
$\omega$ defines the real and reactive power injections
${\bm{p}}_{\textrm{inj}}(\omega),{\bm{q}}_{\textrm{inj}}(\omega)$ via (7)–(9).
We say that ${\bm{\Delta}}(\omega)$ given by

$\displaystyle{\bm{\Delta}}(\omega):=\left({\bm{P}}(\omega),{\bm{Q}}(\omega),{\bm{L}}(\omega),{\bm{W}}(\omega)\right),$
(11)

_explains_ the power injections, when they collectively satisfy (10). To
streamline notation, define the set

$\displaystyle{\cal
G}\left({\bm{p}}_{\textrm{inj}}(\omega),{\bm{q}}_{\textrm{inj}}(\omega)\right):=$
$\displaystyle\left\\{{\bm{\Delta}}(\omega)\text{ explains
}{\bm{p}}_{\textrm{inj}}(\omega),{\bm{q}}_{\textrm{inj}}(\omega)\right\\}.$
(12)

Similarly, we say that a ${\cal B}$-measurable ${\bm{\Delta}}$ _explains_ an
installed capacity configuration ${\bm{\psi}}$, if it explains the induced
injection profiles $\mathbb{P}$-almost surely, i.e.,

$\displaystyle{\bm{\Delta}}\in{\cal
G}\left({\bm{\alpha}}\odot{\bm{\psi}}-{\bm{p}}_{\textrm{D}},{\bm{\eta}}_{\textrm{G}}\odot{\bm{\alpha}}\odot{\bm{\psi}}-{\bm{q}}_{\textrm{D}}\right)\
\mathbb{P}-\text{a.s.}$ (13)

For a given ${\bm{\psi}}$, a scenario $\omega$ respects the network
constraints, if it is explained by ${\bm{\Delta}}(\omega)$ that satisfies

$\displaystyle\underline{{{\bm{W}}}}\leq{\bm{W}}(\omega)\leq\overline{{{\bm{W}}}},\quad[{\bm{P}}(\omega)]^{2}+[{\bm{Q}}(\omega)]^{2}\leq\overline{{{\bm{S}}}}^{2}.$
(14a)

The first among these relations define the set of acceptable nodal voltages
(typically within 5% of nominal voltage levels). The second relation imposes
thermal limitations on line flows. The inequalities are defined element-wise,
where $\underline{{{\bm{W}}}},\overline{{{\bm{W}}}},\overline{{{\bm{S}}}}$ are
constant vectors of appropriate dimensions. In the risk-sensitive hosting
capacity problem, we do not impose such a constraint almost surely, but
constrain the risk associated with the statistics of constraint violation. To
mathematically describe such a constraint, consider some ${\cal B}$-measurable
${\bm{\Delta}}$ that explains the installed capacity configuration
${\bm{\psi}}$. Then, ${\bm{W}}-\overline{{{\bm{W}}}}$ is a random variable,
for which a positive value indicates violation of network constraints. We
impose the constraint

$\displaystyle{\rm CVaR}_{\nu}\left[{\bm{W}}-\overline{{{\bm{W}}}}\right]\leq
0,$ (15)

where ${\rm CVaR}$ for a random variable is defined in (1). If ${\bm{W}}$ has
a smooth cumulative distribution function, the above constraint imposes that
the probability of violating the upper bound on voltage magnitudes is bounded
above by $1-\nu$ and the expectation of ${\bm{W}}-\overline{{{\bm{W}}}}$ over
its $\nu$-tail is bounded above by zero. Thus, such a constraint not only
limits the probability of violations but seeks to control the statistics of
violations. ${\rm CVaR}$ is translation-invariant, implying that (15) holds if
and only if ${\rm CVaR}_{\nu}\left[{\bm{W}}\right]\leq\overline{{{\bm{W}}}}$.

Repeating the above exercise for all network constraints in (14), we define
the following risk-sensitive test for acceptability $\mathscr{A}({\bm{\psi}})$
of the configuration ${\bm{\psi}}\in\mathbb{R}^{n-1}_{+}$. We say that
${\bm{\psi}}$ is acceptable, if there exists ${\cal B}$-measurable
${\bm{\Delta}}$, that satisfies

$\displaystyle{\bm{\Delta}}=({\bm{P}},{\bm{Q}},{\bm{L}},{\bm{W}})\text{
explains }{\bm{\psi}},\ \mathbb{P}-\text{a.s.},$ (16a) $\displaystyle{\rm
CVaR}_{\nu}\left[{\bm{W}}\right]\leq\overline{{{\bm{W}}}},$ (16b)
$\displaystyle{\rm CVaR}_{\nu}[-{\bm{W}}]\leq-\underline{{{\bm{W}}}},$ (16c)
$\displaystyle{\rm
CVaR}_{\gamma}[{\bm{P}}^{2}+{\bm{Q}}^{2}]\leq\overline{{{\bm{S}}}}^{2},$ (16d)

The formulation can be easily extended to consider different risk levels $\nu$
and $\gamma$ for different buses and distribution lines, to indicate the
presence of sensitive equipment. Building on the description of
$\mathscr{A}({\bm{\psi}})$, we next define the risk-sensitive hosting capacity
maximization problem $\mathscr{O}$ that seeks to maximize some function of
installed capacities ${\bm{\psi}}$ over a set of acceptable subsets of
$[0,\overline{{{\bm{\psi}}}}]$. Here, $\overline{{{\bm{\psi}}}}$ encodes the
maximum surface area available for solar panel installations across buses in
$\overline{{{\cal V}}}$.

$\displaystyle\underset{{\bm{\psi}}\in[0,\overline{{{\bm{\psi}}}}],{\bm{\Delta}}}{\text{maximize}}\
\mathds{1}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}{\bm{\psi}},\
\ \text{subject to}\ \
\eqref{eq:risk.hcp.Delta.explains},\eqref{eq:risk.hcp.constraints.voltage.up},\eqref{eq:risk.hcp.constraints.voltage.down},\eqref{eq:risk.hcp.constraints.flow.up}.$
(17)

over ${\bm{\psi}}\in\mathbb{R}^{n-1}_{+}$ and ${\cal B}$-measurable
${\bm{\Delta}}$. The objective function aims to maximize the total installed
capacity of PV panels across the distribution network, where $\mathds{1}$
stands for a vector of all ones of appropriate size.

The feasible set of ${\bm{\Delta}}$ for $\mathscr{A}({\bm{\psi}})$, and hence
for $\mathscr{O}$, is nonconvex, because of (16a). Specifically, this
nonconvexity arises due to the nature of the power flow equations in (10). For
a given power injection
${\bm{p}}_{\textrm{inj}}(\omega),{\bm{q}}_{\textrm{inj}}(\omega)$, the set
${\cal G}$ defined in (12) is generally nonconvex, owing to the quadratic
equality constraint in (10d). It is common practice to _convexify_ the power
flow model, either via linearization or convex relaxation. For example, the
hosting capacity analysis in [19] takes the linearization route. We adopt a
second-order cone programming (SOCP)-based relaxation of the set, similar to
that in [8], and consider the alternate to (10d) for each $\omega$ in

$\displaystyle\left[{\bm{B}}_{+}{\bm{W}}(\omega)\right]\odot{\bm{L}}(\omega)\geq[{\bm{P}}(\omega)]^{2}+[{\bm{Q}}(\omega)]^{2}.$
(18)

Then, define the _convex_ set

$\displaystyle\begin{aligned} {\cal
G}^{\textrm{C}}({\bm{p}}_{\textrm{inj}}(\omega),{\bm{q}}_{\textrm{inj}}(\omega))&:=\left\\{{\bm{\Delta}}(\omega)|({\bm{p}}_{\textrm{inj}}(\omega),{\bm{q}}_{\textrm{inj}}(\omega),{\bm{\Delta}}(\omega))\right.\\\
&\left.\quad\text{ satisfies
}\eqref{eq:hcp.flow.real},\eqref{eq:hcp.flow.reactive},\eqref{eq:hcp.flow.voltage},\eqref{eq:hcp.flow.loss.convex}\right\\}.\end{aligned}$
(19)

The inequality in (18) for each $\omega$ can be written as a collection of
second-order cone constraints using the relationship

$\displaystyle\begin{aligned} &z_{1}z_{2}\geq Z_{1}^{2}+Z_{2}^{2}\\\
&\iff\left\|\begin{pmatrix}2Z_{1},2Z_{2},z_{1}-z_{2}\end{pmatrix}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}\right\|_{2}\leq
z_{1}+z_{2}\\\
&\iff\begin{pmatrix}2Z_{1},2Z_{2},z_{1}-z_{2},z_{1}+z_{2}\end{pmatrix}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}\in{\cal
K}_{\textrm{SO}}\end{aligned}$ (20)

for any scalars $Z_{1},Z_{2},z_{1},z_{2}$, where $\left\|\cdot\right\|_{2}$
computes the 2-norm of a vector and ${\cal K}_{\textrm{SO}}$, defined as

$\displaystyle{\cal
K}_{\textrm{SO}}=\\{({\bm{y}}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}},\upsilon)^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}|\left\|{\bm{y}}\right\|_{2}\leq\upsilon\\},$
(21)

is the convex second-order cone. The set of ${\bm{\Delta}}$’s that belong to
the convex set ${\cal G}^{\textrm{C}}$, $\mathbb{P}$-almost surely, then
becomes convex. This convex relaxation of the set of feasible injections
defined via Kirchhoff’s laws has been widely studied in [20, 21, 22].

The set of ${\bm{\Delta}}(\omega)$ that satisfies the network constraints in
(14) is convex. Thanks to the coherence of the ${\rm CVaR}$ measure (see
[23]), this implies that the set of ${\bm{\Delta}}$ satisfying (16b), (16c),
(16d) is a convex set. Using the notation introduced above, we define the test
for acceptability $\mathscr{A}^{\textrm{C}}({\bm{\psi}})$ with a “relaxed”
power flow model as follows. We say ${\bm{\psi}}$ is acceptable, if there
exists ${\cal B}$-measurable ${\bm{\Delta}}$, that satisfies (16b), (16c),
(16d) and

$\displaystyle{\bm{\Delta}}\in{\cal
G}^{\textrm{C}}\left({\bm{\alpha}}\odot{\bm{\psi}}-{\bm{p}}_{\textrm{D}},{\bm{\eta}}_{\textrm{G}}\odot{\bm{\alpha}}\odot{\bm{\psi}}-{\bm{q}}_{\textrm{D}}\right)\
\mathbb{P}-\text{a.s.}$ (22)

The corresponding convexified hosting capacity maximization problem
$\mathscr{O}^{\textrm{C}}$ over ${\bm{\psi}}\in\mathbb{R}^{n-1}_{+}$ and
${\cal B}$-measurable ${\bm{\Delta}}$ becomes

$\displaystyle\underset{{\bm{\psi}}\in[0,\overline{{{\bm{\psi}}}}],{\bm{\Delta}}}{\text{maximize}}\
\mathds{1}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}{\bm{\psi}},\
\ \text{subject to}\ \
\eqref{eq:risk.hcp.Delta.explains.conv},\eqref{eq:risk.hcp.constraints.voltage.up},\eqref{eq:risk.hcp.constraints.voltage.down},\eqref{eq:risk.hcp.constraints.flow.up}.$
(23)

##  IV A Scenario-Based Approach to Solve $\mathscr{O}^{\textrm{C}}$

The set of ${\bm{\Delta}}$’s satisfying
$\mathscr{A}^{\textrm{C}}({\bm{\psi}})$ and the set of
${\bm{\psi}},{\bm{\Delta}}$ satisfying $\mathscr{O}^{\textrm{C}}$ are convex.
However, finding/optimizing ${\cal B}$-measurable ${\bm{\Delta}}$ can be
challenging. If $\Omega$ is compactly supported, ${\bm{\Delta}}$’s are
infinite-dimensional functions defined over a compact set. To circumvent the
difficulty of optimizing over function spaces, we take a _sampling_ approach
to solve $\mathscr{O}^{\textrm{C}}$.

Consider $K$ independent and identically distributed (i.i.d.) samples
$\omega^{1},\ldots,\omega^{K}$ from $\Omega$, drawn according to $\mathbb{P}$.
With $K$ samples of solar irradiance, real and reactive power demands across
buses in $\overline{{{\cal V}}}$, denoted by the vector tuples
$\\{({\bm{\alpha}}^{k},{\bm{p}}_{\textrm{D}}^{k},{\bm{q}}_{\textrm{D}}^{k})\\}_{k=1}^{K}$,
consider

$\displaystyle{\bm{\Delta}}(\omega^{k})={\bm{\Delta}}^{k}=({\bm{P}}^{k},{\bm{Q}}^{k},{\bm{L}}^{k},{\bm{W}}^{k}).$
(24)

The sampled variant of (22) becomes

$\displaystyle{\bm{\Delta}}^{k}\in{\cal
G}^{\textrm{C}}\left({\bm{\alpha}}^{k}\odot{\bm{\psi}}-{\bm{p}}_{\textrm{D}}^{k},{\bm{\eta}}_{\textrm{G}}\odot{\bm{\alpha}}^{k}\odot{\bm{\psi}}-{\bm{q}}_{\textrm{D}}^{k}\right),$
(25)

for $k=1,\ldots,K$. Leveraging the variational form of ${\rm CVaR}$ from (1),
we write the scenario-based counterpart of (16b) as

$\displaystyle\begin{aligned} &\min_{\overline{{{\bm{w}}}}}\
\left\\{\overline{{{\bm{w}}}}+\frac{1}{1-\nu}\frac{1}{K}\sum_{k=1}^{K}\left[{\bm{W}}^{k}-\overline{{{\bm{w}}}}\right]_{+}\right\\}\leq\overline{{{\bm{W}}}}\\\
&\iff\begin{cases}\overline{{{\bm{w}}}}+\frac{1}{1-\nu}\frac{1}{K}\sum_{k=1}^{K}\overline{{{\bm{t}}}}^{k}_{w}\leq\overline{{{\bm{W}}}},\\\
\overline{{{\bm{t}}}}^{k}_{w}\geq{\bm{W}}^{k}-\overline{{{\bm{w}}}},\
\overline{{{\bm{t}}}}^{k}_{w}\geq 0,\end{cases}\end{aligned}$ (26)

for some
$\overline{{{\bm{w}}}}\in\mathbb{R}^{n},\overline{{{\bm{t}}}}^{k}_{w}\in\mathbb{R}^{n}$,
$k=1,\ldots,K$. Thus, one can drop the sampled ${\rm CVaR}$ constraint, and
instead include the latter pair of deterministic constraints listed above, and
optimize over $\overline{{{\bm{w}}}},\overline{{{\bm{t}}}}^{k}_{w}$,
$k=1,\ldots,K$. We repeat this procedure for (16c) and (16d) and present the
sampled counterpart $\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$ of
$\mathscr{O}^{\textrm{C}}$.

$\displaystyle{\text{maximize}}\quad\mathds{1}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}{\bm{\psi}},$
subject to $\displaystyle\ 0\leq{\bm{\psi}}\leq\overline{{{\bm{\psi}}}},$
(27a) $\displaystyle\
{\bm{\alpha}}^{k}\odot{\bm{\psi}}-{\bm{p}}_{\textrm{D}}^{k}={\bm{\pi}}\left[{\bm{B}}_{+}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}{\bm{P}}^{k}-{\bm{B}}_{-}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}\left({\bm{P}}^{k}-{\bm{R}}\odot{\bm{L}}^{k}\right)\right],$
(27b) $\displaystyle\
{\bm{\eta}}_{\textrm{G}}\odot{\bm{\alpha}}^{k}\odot{\bm{\psi}}-{\bm{q}}_{\textrm{D}}^{k}={\bm{\pi}}\left[{\bm{B}}_{+}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}{\bm{Q}}^{k}-{\bm{B}}_{-}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}\left({\bm{Q}}^{k}-{\bm{X}}\odot{\bm{L}}^{k}\right)\right],$
(27c) $\displaystyle\
{\bm{B}}{\bm{W}}^{k}=2\left[{\bm{R}}\odot{\bm{P}}^{k}+{\bm{X}}\odot{\bm{Q}}^{k}\right]-({\bm{R}}^{2}+{\bm{X}}^{2})\odot{\bm{L}}^{k},$
(27d) $\displaystyle\
\left[{\bm{B}}_{+}{\bm{W}}^{k}\right]\odot{\bm{L}}^{k}\geq[{\bm{P}}^{k}]^{2}+[{\bm{Q}}^{k}]^{2},$
(27e) $\displaystyle\
\overline{{{\bm{w}}}}+\frac{1}{1-\nu}\frac{1}{K}\sum_{k=1}^{K}\overline{{{\bm{t}}}}^{k}_{w}\leq\overline{{{\bm{W}}}},\
\overline{{{\bm{t}}}}^{k}_{w}\geq{\bm{W}}^{k}-\overline{{{\bm{w}}}},$ (27f)
$\displaystyle\
\underline{{{\bm{w}}}}+\frac{1}{1-\nu}\frac{1}{K}\sum_{k=1}^{K}\underline{{{\bm{t}}}}^{k}_{w}\leq-\underline{{{\bm{W}}}},\
\underline{{{\bm{t}}}}^{k}_{w}\geq-{\bm{W}}^{k}-\underline{{{\bm{w}}}},$ (27g)
$\displaystyle\
{\bm{s}}+\frac{1}{1-\gamma}\frac{1}{K}\sum_{k=1}^{K}{\bm{t}}^{k}_{s}\leq\overline{{{\bm{S}}}}^{2},$
(27h) $\displaystyle\
{\bm{t}}^{k}_{s}\geq[{\bm{P}}^{k}]^{2}+[{\bm{Q}}^{k}]^{2}-{\bm{s}},$ (27i)
$\displaystyle\
\overline{{{\bm{t}}}}^{k}_{w}\geq{\bm{0}},\,\underline{{{\bm{t}}}}^{k}_{w}\geq{\bm{0}},\,{\bm{t}}^{k}_{s}\geq{\bm{0}},\
\text{for }k=1,\dots,K$ (27j)

over the variables
${\bm{\psi}},{\bm{P}}^{k},{\bm{Q}}^{k},{\bm{W}}^{k},{\bm{L}}^{k},\overline{{{\bm{w}}}},\underline{{{\bm{w}}}},{\bm{s}},\overline{{{\bm{t}}}}^{k}_{w},\underline{{{\bm{t}}}}^{k}_{w},{\bm{t}}^{k}_{s}$.
All constraints in $\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$ are linear,
except those in (27e) and (27i). One can utilize (20) to cast them as second-
order cone constraints. As a result,
$\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$ becomes an SOCP–a _convex_
optimization problem. Such SOCP’s can be solved using off-the-shelf solvers
that employ interior-point algorithms.

The CVaR-based formulation of $\mathscr{O}^{\textrm{C}}$ and its scenario-
based variant in $\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$ provide a
structured framework to identify optimal installation configurations. The
statistics of network constraint violation across scenarios impact the choice
of the optimal installation profile. While our approach shares many parallels
with the BayesOpt framework with probabilistic constraints in [8], there are
important differences. To explain the difference, consider
$\mathscr{O}^{\textrm{C}}$ with probabilistic constraint enforcement, i.e.,
replace ${\rm CVaR}_{\nu}[{\bm{W}}]\leq\overline{{{\bm{W}}}}$ with
$\mathbb{P}\\{{\bm{W}}-\overline{{{\bm{W}}}}>0\\}\leq 1-\nu$. With scenarios,
the probabilistic constraint amounts to imposing

$\frac{1}{K}\sum_{k=1}^{K}\mathbb{I}_{\\{{\bm{W}}^{k}-\overline{{{\bm{W}}}}>0\\}}\leq
1-\nu,$ (28)

where $\mathbb{I}$ is the indicator function that evaluates to unity when its
argument is true, and becomes zero otherwise. Tackling such constraints with
indicator functions within our framework requires the solution of a nonconvex
mixed-integer SOCP, which is typically more computationally challenging than
solving the SOCP problem $\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$ with
CVaR-based constraints. Note that a solver for
$\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$ will typically update all
optimization variables together towards an optimal solution. In contrast, the
probability of constraint violation in the left-hand-side of (28) must be
computed, each time the candidate installation profile ${\bm{\psi}}$ changes
in search of an optimal solution. More precisely, as ${\bm{\psi}}$ is updated
to ${\bm{\psi}}^{\prime}$, the already-computed ${\bm{\Delta}}^{k}$ for
scenario $\omega^{k}$ with installation profile ${\bm{\psi}}$ is discarded,
and must be recomputed from scratch with an installation profile
${\bm{\psi}}^{\prime}\neq{\bm{\psi}}$. Lacking a convex structure for the set
of ${\bm{\Delta}}^{k}$’s that satisfy (28), it becomes difficult to recycle
the computation of ${\bm{\Delta}}^{k}$ _across_ candidate installation
profiles. Our formulation does not explicitly evaluate CVaR over scenarios for
each ${\bm{\psi}}$. Rather, we leverage the variational characterization of
CVaR in (1) to formulate $\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$ over
${\bm{\psi}}$, ${\bm{\Delta}}^{k}$’s and auxiliary variables, where a solver
can exploit the convexity of the problem to iteratively update these variables
together.

The lack of convexity in the left-hand-side of (28) makes it challenging to
find natural directions to update ${\bm{\psi}}$. The Bayesian optimization
framework in [8] tackles this challenge through a combination of two steps.
First, the constraint in (28) is included through a penalty in the objective
function to construct an unconstrained optimization problem. Then, the
resulting objective is minimized via kernel smoothing of the objective
function. Running such algorithms requires careful choice of the penalty
factor and the kernel. Our CVaR-based formulation being convex, avoids such
manual tuning.

Our formulation of $\mathscr{O}^{\textrm{C}}$ and
$\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$ utilizes the SOCP-based
relaxation of the set of feasible injections, characterized by Kirchhoff’s
laws. The setup is flexible in that one can utilize any “convexified” power
flow model to obtain a convex programming formulation of the CVaR-sensitive
hosting capacity problem. This approach offers a sharp contrast to those that
rely heavily on the choice of a particular power flow model. For example, the
sample reduction technique in [19] crucially depends on a linear power flow
model, where constraint violations are processed through a quadratic penalty
and optimized, permitting the use of multi-parametric quadratic programming
theory to assess the feasibility of an installation configuration
${\bm{\psi}}$. Such a method does not naturally generalize to SOCP-based power
flow models.

###  IV-A Numerical Experiments with
$\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$

We solved $\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$ on two networks: a
3-bus network for illustrative purposes and the Southern California Edison’s
56-bus network from [24], studied previously in [8], to demonstrate
scalability. Simulations on the 3-bus network were performed on an i5-4690k
machine with 16GB RAM, while the 56-bus simulations were performed on an AWS
m4.10xlarge instance. All results were generated using MOSEK version 9.3. Load
data and DER energy profiles were taken at 10 minute intervals over a 365-day
period (52.56K intervals) using data from Southern California Edison. For all
experiments, DER installation power factors were set to $0.97$ ($\eta=0.251$),
and the feeder voltage was fixed at $w_{1}=1$ p.u. Data and code are available
on GitHub at https://github.com/amadavan/RSHC.

####  IV-A1 On a three-bus network

In Figure 3a, we examine the effect of the risk parameters ($\gamma=\nu$) on
both the solution and the run-time. The optimal aggregate capacity
$\mathds{1}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}{\bm{\psi}}^{\star}$
decreases with risk aversion. This is anticipated as a higher risk parameter
leads to a tighter risk-sensitive constraint, shrinking the feasible set. In
turn, this leads to a decreased maximum total installed capacity. We believe
that increasing run-times are a consequence of the tightening constraints.

(a)  (b)

Figure 3: Experiments on the 3-bus network, showing LABEL:sub@fig:3bus.time
the effect of the CVaR parameter $\gamma$ on the maximum capacity and runtime,
and LABEL:sub@fig:3bus.ssa the deviation in solutions with increasing samples
with $\gamma=\nu=0.8$.

The problem dimension scales linearly with the number of scenarios $K$. Each
scenario introduces $9n-5$ constraints and $7n-4$ variables, where $n$ is the
number of buses. In Figure 3b, we evaluate the effect of sub-sampling on the
resulting solution, plotting the mean optimal solution with deviations over 20
runs per sample size. Fewer samples result in greater variance in the
solution, but reduced run-times.

####  IV-A2 On a 56-bus distribution network

Figure 4 shows the effect of the risk parameter $\gamma$ (keeping $\nu$
constant) on both run-times and optimal objective values attained with
$K=$52.56K scenarios. As in the 3-bus example, the objective value decreases
with $\nu$. The algorithm run-times are 1-2 hours for each value of risk
parameter considered, and increase significantly for $\gamma$ close to $1$.

(a)  (b)

Figure 4: Experiments on the 56-bus network with $\nu=0.9$, showing
LABEL:sub@fig:56bus.time the effect of the CVaR parameter $\gamma$ on the
maximum capacity and runtime, and LABEL:sub@fig:56bus.ssa the deviation in
solutions with increasing samples with $\gamma=0.8$.

As previously noted, for a random variable $\zeta$, the constraint ${\rm
CVaR}{}_{\delta}[{\zeta}]\leq 0$ is a convex inner approximation of the
constraint $\mathbb{P}[\zeta\geq 0]\leq 1-\delta$. In other words, a
configuration ${\bm{\psi}}$ for which the CVaR of voltage variations and power
flows satisfies the constraints in $\mathscr{O}^{\textrm{C}}$, they must
satisfy the corresponding probabilistic constraints. To illustrate the same
empirically with the optimal solution ${\bm{\psi}}^{\star}$, consider the
solution obtained with $\nu=0.9$ and $\gamma=0.8$. For this solution, we
evaluate the proportion of the 52.56K ${\bm{\Delta}}^{k}$’s for which each of
the 110 upper and lower voltage limits and 55 power flow limits are violated.
Figure 5 shows the cumulative distribution histogram of the proportions of
voltage violations across the buses and power flow violations across the
lines. The support of voltage violation proportions lie fully to the left of
$1-\nu$ and the power flow violation proportions to the left of $1-\gamma$.
That is, our CVaR-based constraint enforcement automatically guarantees
probabilistic constraint enforcement.

Figure 5: Cumulative frequency histogram of proportions of 52,560 scenarios
over which upper/lower voltage constraints across buses and power flow
constraints across distribution lines for the 56-bus network are violated with
$\nu=0.9$, $\gamma=0.8$.

While our reported solutions used MOSEK, we have compared our results with an
implementation in Gurobi. For high risk-aversion parameters, e.g.,
$\nu=\gamma=0.9$, Gurobi took almost twice as long to produce an optimal
solution, compared to MOSEK, even with small number of scenarios. However, in
performing our simulations, MOSEK tended to yield a stable optimal solution
($\approx$ the same result as Gurobi) the fastest, when the objective function
was scaled by a large parameter $\sim 10^{9}$. Run-times reduced by as much as
by a factor of two with such scaling. The Gurobi implementation remained
immune to objective function scaling.

Compared to the BayesOpt algorithm in [8], our implementation is _centralized_
, and naturally requires more memory than a solution architecture that solves
scenario-wise power flow problems. A side-by-side comparison with earlier
works remains difficult, given that optimization with probabilistic constraint
enforcement has no optimality guarantees and its outcome depends strongly on
implementation details.

One might surmise that our formulation of
$\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$ can be parallelized using an
algorithmic architecture such as generalized Benders’ decomposition in [25].
Indeed, one can devise such a method that processes scenarios $\omega^{k}$’s
in parallel. However, our implementation of generalized Benders’ decomposition
for $\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$ led to increasingly larger
master problems, resulting from feasibility cuts from individual scenarios,
that proved computationally prohibitive. We abandoned that route, given that
the centralized implementation ran in reasonable time. A more careful
implementation with lazy constraint generation protocols is left for future
endeavors.

##  V An Incremental Algorithm to Solve
$\mathscr{A}^{\textrm{C}}({\bm{\psi}})$ with Scenarios

Setting aside the issue of optimality, we now study a scenario-based approach
to solve $\mathscr{A}^{\textrm{C}}({\bm{\psi}})$ and certify whether the risk
of network constraint violations associated with certain configuration of
solar installation capacities ${\bm{\psi}}\in\mathbb{R}^{n-1}_{+}$ is
acceptable. Recall that ${\bm{\psi}}$ is deemed acceptable, if there exists a
${\cal B}$-measurbable ${\bm{\Delta}}$ that satisfies (16b), (16c), (16d) and
(22). Similar to the design of $\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$,
one can certify the acceptability of ${\bm{\psi}}$ with a list of $K$ i.i.d.
samples $\omega^{1},\ldots,\omega^{K}$ from $\Omega$ via the convex
optimization problem
$\mathscr{A}^{\textrm{C}}_{\textrm{sample}}({\bm{\psi}})$, described by

$\displaystyle\mathscr{A}^{\textrm{C},\star}_{\textrm{sample}}({\bm{\psi}}):={\text{maximize}}\
0,\ \text{subject to }\ \eqref{eq:rshcp.net.1}-\eqref{eq:rshcp.nonneg},$ (29)

over ${\bm{P}}^{k},{\bm{Q}}^{k},{\bm{W}}^{k},{\bm{L}}^{k}$,
$\overline{{{\bm{w}}}},\underline{{{\bm{w}}}},{\bm{s}},\overline{{{\bm{t}}}}^{k}_{w},\underline{{{\bm{t}}}}^{k}_{w},{\bm{t}}^{k}_{s}$
for $k=1,\ldots,K$. If (29) admits a feasible point, then ${\cal
A}^{\star}({\bm{\psi}})=0$ and ${\bm{\psi}}$ is acceptable. Otherwise, ${\cal
A}^{\star}({\bm{\psi}})=-\infty$ and ${\bm{\psi}}$ is unacceptable.

Figure 6: Inner (green) and outer (white) sets portray
${\bm{\Psi}}_{\textrm{in}}$ and ${\bm{\Psi}}_{\textrm{out}}$, respectively,
that approximate the set of acceptable installations ${\bm{\Psi}}$ (gray). The
outer frame encodes the set $[0,\overline{{{\bm{\psi}}}}]$. The sets
$\mathbb{A}_{\psi}$ and $\mathbb{U}_{\psi}$ comprise the points with green and
red borders, respectively.

.

One can always solve (29) to certify whether ${\bm{\psi}}$ is acceptable. This
problem has a similar time-complexity as
$\mathscr{O}^{\textrm{C}}_{\textrm{sample}}$. In the sequel, we develop an
_incremental algorithm_ to certify acceptability. That is, we first attempt to
leverage the knowledge of a list of acceptable and unacceptable configurations
to certify whether a new candidate configuration is acceptable. Our design is
such that this check is often much faster than solving (29). If unsuccessful,
we then solve (29) and use the new configuration to update the knowledge about
acceptable and unacceptable configurations. As will be evident, such a method
drastically reduces the instances of ${\bm{\psi}}$’s for which (29) must be
solved to certify acceptability with more tests. Thus, the more we test, the
less we require to solve (29). Such a testing paradigm is useful, especially
to test acceptability of installation configurations that are close to each
other, e.g., those that lie on likely solar adoption paths.

Let ${\bm{\Psi}}$ denote the set of all acceptable configurations of
physically realizable installed solar capacities, i.e.,
${\bm{\Psi}}\subseteq[0,\overline{{{\bm{\psi}}}}]$. Thanks to the convexity of
CVaR, ${\bm{\Psi}}$ is a convex set, shaded in gray in Figure 6. We maintain
polyhedral inner and outer approximations ${\bm{\Psi}}_{\textrm{in}}$ and
${\bm{\Psi}}_{\textrm{out}}$, respectively, of ${\bm{\Psi}}$ satisfying

$\displaystyle{\bm{\Psi}}_{\textrm{in}}\subseteq{\bm{\Psi}}\subseteq{\bm{\Psi}}_{\textrm{out}}\subseteq[0,\overline{{{\bm{\psi}}}}].$
(30)

These sets are visualized in Figure 6. We construct
${\bm{\Psi}}_{\textrm{in}}$ and ${\bm{\Psi}}_{\textrm{out}}$ incrementally
from prior acceptability tests. For a new candidate ${\bm{\psi}}$, if
${\bm{\psi}}\in{\bm{\Psi}}_{\textrm{in}}$, then ${\bm{\psi}}$ is acceptable.
On the other hand, if ${\bm{\psi}}\notin{\bm{\Psi}}_{\textrm{out}}$, then
${\bm{\psi}}$ is unacceptable. If neither of these tests certify the
acceptability of ${\bm{\psi}}$, then we run (29). The result of the
optimization program either grows ${\bm{\Psi}}_{\textrm{in}}$ or shrinks
${\bm{\Psi}}_{\textrm{out}}$, thus expanding the space of ${\bm{\psi}}$’s,
where acceptability can be certified without running (29). In the rest of this
section, we outline the incremental construction of
${\bm{\Psi}}_{\textrm{in}},{\bm{\Psi}}_{\textrm{out}}$ and how we test
${\bm{\psi}}\in{\bm{\Psi}}_{\textrm{in}}$ or
${\bm{\psi}}\notin{\bm{\Psi}}_{\textrm{out}}$. The method for testing
acceptability of ${\bm{\psi}}$ is summarized in Algorithm 1.

With a list
$\mathbb{A}_{\psi}:=\\{{\bm{\psi}}^{a}[1],\ldots,{\bm{\psi}}^{a}[M]\\}$ of
acceptable configurations, we construct ${\bm{\Psi}}_{\textrm{in}}$ as their
_convex hull_ ,

$\displaystyle\begin{gathered}{\bm{\Psi}}_{\textrm{in}}=\textrm{conv}(\mathbb{A}_{\psi})=\left\\{\sum_{i=1}^{M}\beta_{i}{\bm{\psi}}^{a}[i]:{\bm{\beta}}\in\mathbb{R}^{M}_{+},\mathds{1}^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}{\bm{\beta}}=1\right\\}.\end{gathered}$
(32)

Testing if ${\bm{\psi}}\in{\bm{\Psi}}_{\textrm{in}}$ can be solved as a linear
program over ${\bm{\beta}}\in\mathbb{R}^{M}$. If ${\bm{\psi}}$ is certified
acceptable from running (29), then
$\mathbb{A}_{\psi}\leftarrow\mathbb{A}_{\psi}\cup{\bm{\psi}}$ and define again
${\bm{\Psi}}_{\textrm{in}}=\textrm{conv}(\mathbb{A}_{\psi})$. We grow
${\bm{\Psi}}_{\textrm{in}}$ incrementally, starting from the null set.

Next, we describe how we construct ${\bm{\Psi}}_{\textrm{out}}$, starting with
$[0,\overline{{{\bm{\psi}}}}]$. In the event that (29) finds ${\bm{\psi}}$
unacceptable, we construct a _feasibility cut_ , borrowing the idea from
generalized Benders’ decomposition in [25]. To explain cut generation, define
${\bm{x}}$ as the vectorized concatenation of all variables in (29) and write
(29) as

$\displaystyle\mathscr{A}^{\textrm{C},\star}_{\textrm{sample}}({\bm{\psi}})=\,$
$\displaystyle\underset{{\bm{x}}}{\text{maximize}}$ $\displaystyle 0,$ subject
to
$\displaystyle{\bm{C}}_{\psi}{\bm{\psi}}={\bm{D}}_{\psi}{\bm{x}}+{\bm{E}}_{\psi},$
(33a) $\displaystyle{\bm{D}}{\bm{x}}\leq{\bm{E}},$ (33b)
$\displaystyle\begin{pmatrix}{\bm{F}}_{i}{\bm{x}}+{\bm{{\bm{G}}}}_{i}\\\
{\bm{f}}^{{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}}_{i}{\bm{x}}+g_{i}\end{pmatrix}\in{\cal
K}_{\textrm{SO}},$ (33c) $\displaystyle i=1,\ldots,I.$

The equalities (27b) and (27c) are written compactly as (33a). The equality in
(27d) can be written as two inequalities. These, together with (27f)–(27h) and
(27j) become examples of (33b). The nonlinear constraints in (27e) and (27i)
are written as $I$ SO cone constraints (33c) using (20). To lighten notation,
assume $I=1$ and consider the following Lagrangian function

$\displaystyle\begin{aligned} &{\cal
L}({\bm{x}},{\bm{\lambda}},{\bm{\mu}},{\bm{\mu}}_{1},\mu_{2})\\\
&={\bm{\lambda}}^{{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}}({\bm{C}}_{\psi}{\bm{\psi}}-{\bm{D}}_{\psi}{\bm{x}}-{\bm{E}}_{\psi})+{\bm{\mu}}^{{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}}({\bm{D}}{\bm{x}}-{\bm{E}})\\\
&\qquad-{\bm{\mu}}_{1}^{{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}}({\bm{F}}{\bm{x}}+{\bm{G}})-\mu_{2}({\bm{f}}^{{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}}{\bm{x}}+g)\end{aligned}$
(34)

for Lagrange multipliers ${\bm{\lambda}},{\bm{\mu}},{\bm{\mu}}_{1},\mu_{2}$ of
compatible dimensions. Then, (33) can be written as

$-\mathscr{A}^{\textrm{C},\star}_{\textrm{sample}}({\bm{\psi}})=\underset{{\bm{x}}}{\text{min}}\
\underset{\begin{subarray}{c}{\bm{\lambda}},{\bm{\mu}}\geq 0,\\\
({\bm{\mu}}_{1},\mu_{2})\in{\cal K}_{\textrm{SO}}\end{subarray}}{\text{max}}\
{\cal L}({\bm{x}},{\bm{\lambda}},{\bm{\mu}},{\bm{\mu}}_{1},\mu_{2}).$ (35)

The above representation leverages the fact that ${\cal K}_{\textrm{SO}}$ is a
self-dual convex cone (see [26] for details). Weak duality implies

$-\mathscr{A}^{\textrm{C},\star}_{\textrm{sample}}({\bm{\psi}})\geq\underset{\begin{subarray}{c}{\bm{\lambda}},{\bm{\mu}}\geq
0,\\\ ({\bm{\mu}}_{1},\mu_{2})\in{\cal
K}_{\textrm{SO}}\end{subarray}}{\text{max}}\ \
\underset{{\bm{x}}}{\text{min}}\ \ {\cal
L}({\bm{x}},{\bm{\lambda}},{\bm{\mu}},{\bm{\mu}}_{1},\mu_{2}),$ (36)

where the right hand side of the above inequality equals the optimal cost of
the dual problem. Solving the inner minimization over ${\bm{x}}$ in the above
relation, yields

$\displaystyle-\mathscr{A}^{\textrm{C},\star}_{\textrm{sample}}({\bm{\psi}})\geq$
$\displaystyle\underset{{\bm{\lambda}},{\bm{\mu}},{\bm{\mu}}_{1},\mu_{2}}{\text{maximum}}$
$\displaystyle({\bm{C}}_{\psi}{\bm{\psi}}-{\bm{E}}_{\psi})^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}{\bm{\lambda}}$
(37)
$\displaystyle\quad-{\bm{E}}^{{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}}{\bm{\mu}}-{\bm{G}}^{{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}}{\bm{\mu}}_{1}-\mu_{2}g,$
subject to
$\displaystyle{\bm{D}}^{{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}}{\bm{\mu}}={\bm{D}}_{\psi}^{{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}}{\bm{\lambda}}+{\bm{F}}^{{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}}{\bm{\mu}}_{1}+{\bm{f}}\mu_{2},$
$\displaystyle{\bm{\mu}}\geq
0,\,({\bm{\mu}}_{1}^{{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}},\mu_{2})^{\mathchoice{\raisebox{0.75346pt}{$\displaystyle\intercal$}}{\raisebox{0.75346pt}{$\textstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptstyle\intercal$}}{\raisebox{0.75346pt}{$\scriptscriptstyle\intercal$}}}\in{\cal
K}_{\textrm{SO}}.$

The dual program for (33) is an SOCP, for which the origin is always feasible.
The feasible set for this dual problem is a convex cone.111${\cal S}$ is a
convex cone, if $s_{1},s_{2}\in{\cal S},\kappa_{1},\kappa_{2}\geq 0$ implies
$\kappa_{1}s_{1}+\kappa_{2}s_{2}\in{\cal S}$. Call the objective function of
this problem as ${\cal
D}({\bm{\lambda}},{\bm{\mu}},{\bm{\mu}}_{1},\mu_{2};{\bm{\psi}})$. Then, for
any feasible
${\bm{\lambda}}^{\circ},{\bm{\mu}}^{\circ},{\bm{\mu}}_{1}^{\circ},\mu_{2}^{\circ}$
yielding a positive dual objective, we have

$\displaystyle\begin{aligned} &{\cal
D}({\bm{\lambda}}^{\circ},{\bm{\mu}}^{\circ},{\bm{\mu}}_{1}^{\circ},\mu_{2}^{\circ};{\bm{\psi}})>0\\\
&\implies\lim_{\kappa\to\infty}{\cal
D}(\kappa{\bm{\lambda}}^{\circ},\kappa{\bm{\mu}}^{\circ},\kappa{\bm{\mu}}_{1}^{\circ},\kappa\mu_{2}^{\circ};{\bm{\psi}})=\infty\\\
&\implies\mathscr{A}^{\textrm{C},\star}_{\textrm{sample}}({\bm{\psi}})=-\infty.\end{aligned}$
(38)

In other words, existence of a dual feasible point with positive objective
implies the existence of a _dual improving ray_ , along which the dual
objective grows unbounded. In turn, existence of such a ray with a given
${\bm{\psi}}$ certifies that (33) is infeasible and ${\bm{\psi}}$ is
unacceptable. Assuming that (33) is strongly infeasible (see [27]) for all
unacceptable ${\bm{\psi}}$, each such ${\bm{\psi}}$ must admit a dual
improving ray. Interior point solvers for (29) typically produce a dual
improving ray, when infeasible. See [28] for the mechanics of producing such a
ray. Consider a list of unacceptable configurations
$\mathbb{U}_{\psi}:=\left\\{{\bm{\psi}}^{u}[1],\ldots,{\bm{\psi}}^{u}[M^{\prime}]\right\\}$.
For an unacceptable ${\bm{\psi}}^{u}[i]\in\mathbb{U}_{\psi}$, let the dual
improving ray be given by
${\bm{\lambda}}^{\circ}[i],{\bm{\mu}}^{\circ}[i],{\bm{\mu}}_{1}^{\circ}[i],\mu_{2}^{\circ}[i]$.
Then, define

$\displaystyle\begin{aligned}
{\bm{\Psi}}_{\textrm{out}}&=\left\\{{\bm{\psi}}:{\cal
D}({\bm{\lambda}}^{\circ}[i],{\bm{\mu}}^{\circ}[i],{\bm{\mu}}_{1}^{\circ}[i],\mu_{2}^{\circ}[i];{\bm{\psi}})\leq
0,\right.\\\
&\qquad\qquad\left.i=1,\ldots,M^{\prime},{\bm{\psi}}\in[0,\overline{{{\bm{\psi}}}}]\right\\}.\end{aligned}$
(39)

If ${\bm{\psi}}$ violates at least one among the $M$ affine constraints, then
${\bm{\psi}}\notin{\bm{\Psi}}_{\textrm{out}}$. When (29) returns a new
unacceptable ${\bm{\psi}}$, then
$\mathbb{U}_{\psi}\leftarrow\mathbb{U}_{\psi}\cup{\bm{\psi}}$, and the cut
defined via its dual improving ray is added to the description of
${\bm{\Psi}}_{\textrm{out}}$ via (39). Addition of this cut shrinks
${\bm{\Psi}}_{\textrm{out}}$.

Algorithm 1 crucially relies on the convex nature of the set ${\bm{\Psi}}$–a
property that exploits the properties of the CVaR measure. With a
probabilistic constraint enforcement as in [8], this set may _not_ be convex,
rendering the construction of polyhedral inner/outer approximations of
${\bm{\Psi}}$ challenging.

Algorithm 1 Incremental Scenario-Based Test
$\mathscr{A}^{\textrm{C}}_{\textrm{sample}}$.

1:Input:  ${\bm{\psi}}$,
$(\mathbb{A}_{\psi},{\bm{\Psi}}_{\textrm{in}}),(\mathbb{U}_{\psi},{\bm{\Psi}}_{\textrm{out}})$.

2:if ${\bm{\psi}}\notin{\bm{\Psi}}_{\textrm{out}}$ then return unacceptable.

3:else if ${\bm{\psi}}\in{\bm{\Psi}}_{\textrm{in}}$ then return acceptable.

4:else

5: Solve (29).

6: if infeasible and dual improving ray
$({\bm{\lambda}}^{\circ},{\bm{\mu}}^{\circ},{\bm{\mu}}_{1}^{\circ},\mu_{2}^{\circ})$
exists then

7: $\mathbb{U}_{\psi}\leftarrow\mathbb{U}_{\psi}\cup{\bm{\psi}}$, update
${\bm{\Psi}}_{\textrm{out}}$ via (39).

8: return unacceptable.

9: else

10: $\mathbb{A}_{\psi}\leftarrow\mathbb{A}_{\psi}\cup{\bm{\psi}}$, update
${\bm{\Psi}}_{\textrm{in}}$ via (32).

11: return acceptable.

12: end if

13:end if

Figure 7: Visualization of ${\bm{\Psi}}_{\textrm{in}}$ (green) and
$[0,\overline{{{\bm{\psi}}}}]\setminus{\bm{\Psi}}_{\textrm{out}}$ (red)
regions after testing for acceptability of 10, 50, 100 and 1000 configurations
on the 3-bus network example with $\nu=\gamma=0.6$ on the top and
$\nu=\gamma=0.8$ on the bottom. Hollow points required solving (29), while
solid points could be certified using
${\bm{\Psi}}_{\textrm{in}}/{\bm{\Psi}}_{\textrm{out}}$.

###  V-A Numerical Experiments with
$\mathscr{A}^{\textrm{C}}_{\textrm{sample}}({\bm{\psi}})$

####  V-A1 On the three-bus network

We tested acceptability of uniformly randomly generated installation
configurations $(\psi_{2},\psi_{3})\in[0,4]^{2}$ p.u. under risk parameters
$\nu=\gamma=0.8$. At the beginning of the testing process,
${\bm{\Psi}}_{\textrm{in}}$ and ${\bm{\Psi}}_{\textrm{out}}$ do not contain
relevant information to judge acceptability of candidate configurations, and
so (29) gets solved for these configurations. However, as more configurations
are tested, ${\bm{\Psi}}_{\textrm{in}}$ grows and ${\bm{\Psi}}_{\textrm{out}}$
shrinks, as Figure 7 reveals. As a result, more configurations get quickly
certified as acceptable/unacceptable using ${\bm{\Psi}}_{\textrm{in}}$ and
${\bm{\Psi}}_{\textrm{out}}$. As Table I shows, the run-time differences
between solving (29) and testing via ${\bm{\Psi}}_{\textrm{in}}$ and
${\bm{\Psi}}_{\textrm{out}}$ are substantial.

TABLE I: Summary of results from acceptability tests over the 3-bus distribution network example. Test | Mean Runtime (s) | Frequency
---|---|---
${\bm{\psi}}\in{\bm{\Psi}}_{\textrm{in}}$ test | $9.678\times 10^{-4}$ | 451
${\bm{\psi}}\notin{\bm{\Psi}}_{\textrm{out}}$ test | $1.195\times 10^{-6}$ | 489
Solving (29), when ${\bm{\psi}}\in{\bm{\Psi}}_{\textrm{in}}$ | $46.288$ | 47
Solving (29), when ${\bm{\psi}}\notin{\bm{\Psi}}_{\textrm{out}}$ | $133.107$ | 13

Only a small fraction (6%) of the 1000 tested configurations, per Table I,
needed the solution of (29) to certify whether they are acceptable. The sets
${\bm{\Psi}}_{\textrm{in}}$ and ${\bm{\Psi}}_{\textrm{out}}$ therefore
effectively capture the knowledge of past tests, highlighting the efficacy of
our incremental approach to solve
$\mathscr{A}^{\textrm{C}}_{\textrm{sample}}({\bm{\psi}})$.

Figure 7 portrays how the sets ${\bm{\Psi}}_{\textrm{in}}$ and
${\bm{\Psi}}_{\textrm{out}}$ change with tests. As one might expect, the set
of acceptable configurations with lower risk parameters ($\nu=\gamma=0.6$) in
the top half of Figure 7 is larger than that with higher risk parameters
($\nu=\gamma=0.8$) in the bottom half of Figure 7. These figures with
increasing number of points were generated with the same set of test
configurations, and hence, illustrate how they differently contribute to the
construction of the sets ${\bm{\Psi}}_{\textrm{in}}$ and
${\bm{\Psi}}_{\textrm{out}}$ with different risk parameters.

In practice, installation configurations chosen for certification will likely
not be generated uniformly randomly. Rather, these configurations will lie
along projected solar adoption paths, starting from an initial acceptable
configuration. These capacities being nearby, we expect
${\bm{\Psi}}_{\textrm{in}}$ to be effective in quickly declaring acceptability
for most points. Again, as the adoption path leaves ${\bm{\Psi}}$, a few tests
near the boundary should construct cutting planes that later points on the
trajectory can be quickly ruled out using ${\bm{\Psi}}_{\textrm{out}}$.

####  V-A2 On the 56-bus distribution network

We ran our tests for acceptable configurations with 15K scenarios over
capacity configurations that are uniformly distributed within $[-0.2,0.2]^{5}$
p.u. box around ${\bm{\psi}}^{\star}$ obtained with $\nu=0.9$ and
$\gamma=0.8$. We first ran 32 tests corresponding to the corners of this
5-dimensional hypercube. The corner configurations of the hypercube provide
valuable information to the construction of ${\bm{\Psi}}_{\textrm{in}}$ and
${\bm{\Psi}}_{\textrm{out}}$. After these 32 corners, we considered an
additional $118$ data points sampled uniformly over the aforementioned
hypercube. Table II records the results of our numerical experiments. Checking
acceptability via ${\bm{\Psi}}_{\textrm{in}}$ or ${\bm{\Psi}}_{\textrm{out}}$
takes less than a second, while running (29) takes substantially longer than
for the 3-bus network example, owing to the increased network size. Even with
a small number (118) of tests, majority of the points are certified to be
acceptable or unacceptable using ${\bm{\Psi}}_{\textrm{in}}$ and
${\bm{\Psi}}_{\textrm{out}}$, without requiring explicit solution of (29). The
run-time differences and the frequency of certification using
${\bm{\Psi}}_{\textrm{in}}$ and ${\bm{\Psi}}_{\textrm{out}}$ illustrate the
value of our incremental testing procedure for acceptability of installation
configurations for large distribution networks. The gains in run-times can
only increase with larger number of scenarios, network sizes, and number of
tests conducted.

TABLE II: Summary of results from acceptability tests over the 56-bus distribution network example. Test | Mean Runtime (s) | Frequency
---|---|---
${\bm{\psi}}\in{\bm{\Psi}}_{\textrm{in}}$ test | $1.473\times 10^{-3}$ | 6
${\bm{\psi}}\notin{\bm{\Psi}}_{\textrm{out}}$ test | $1.419\times 10^{-6}$ | 79
Solving (29), when ${\bm{\psi}}\in{\bm{\Psi}}_{\textrm{in}}$ | $978.764$ | 49
Solving (29), when ${\bm{\psi}}\notin{\bm{\Psi}}_{\textrm{out}}$ | $1297.369$ | 16

##  VI Conclusions

Having already experienced several years of sustained growth, the continued
adoption of distributed and utility scale solar generation is in many
jurisdictions not only expected, but now mandated by legislation. The primary
question facing distribution grid operators and planners has therefore become
not if, but how and when solar will be integrated at scale. HCA assesses the
limits of a network for safely accommodating solar generation capacity. Given
the highly intermittent nature of solar generation, together with changing
consumption patterns, moving forward, notions of safety must account for the
statistics of both present and forecasted operating conditions. In order to
address the scalability issues arising from consideration of large numbers of
scenarios, in this work we proposed an HCA formulation based on the use of
CVaR as a network constraint risk measure. CVaR provides guarantees as to the
frequency and extent of constraint violations, and moreover possesses
desirable mathematical properties yielding a tractable approach to HCA.
Further leveraging these properties, we provided an expedited method for
assessing the acceptability of a given installation configuration, using prior
knowledge of installations certified as either acceptable or unacceptable.

There are a number of interesting research directions we aim to pursue. First,
while our optimization approach requires no specialized algorithm design or
tuning, the resulting problems are still large and thus, require considerable
memory resources. We want to develop decomposition techniques that alleviate
such difficulties, without sacrificing quality of the solutions obtained.
Second, we aim to explore the potential of scenario reduction, e.g., via
clustering, to help relax the computational burden. Finally, we want to model
and evaluate the effect of active control of distributed energy resources in
CVaR-sensitive HCA with possibly unbalanced three-phase distribution networks.

## References

  * [1] S. Joshi, S. Mittal, P. Holloway, P. R. Shukla, B. Ó Gallachóir, and J. Glynn, “High resolution global spatiotemporal assessment of rooftop solar photovoltaics potential for renewable electricity generation,” _Nature communications_ , vol. 12, no. 1, pp. 1–15, 2021.
  * [2] J. Spector, “Hawaii is ahead of schedule for renewable power adoption,” February 2021, [Online; posted 22-February-2022]. [Online]. Available: https://www.greentechmedia.com/articles/read/hawaii-is-ahead-of-schedule-for-renewable-power-adoption
  * [3] U. E. I. Administration, “California’s curtailments of solar electricity generation continue to increase,” https://www.eia.gov/todayinenergy/detail.php?id=49276, accessed: 2022-02-22.
  * [4] B. Azibek, A. Abukhan, H. K. Nunna, B. Mukatov, S. Kamalasadan, and S. Doolla, “Hosting capacity enhancement in low voltage distribution networks: challenges and solutions,” in _2020 IEEE International Conference on Power Electronics, Smart Grid and Renewable Energy (PESGRE2020)_. IEEE, 2020, pp. 1–6.
  * [5] “Small generator interconnection procedures (sgip),” FERC, Tech. Rep., 2020.
  * [6] N. Baldenko and S. Behzadirafi, “Determination of photovoltaic hosting capacity on radial electric distribution feeders,” in _2016 IEEE International Conference on Power System Technology (POWERCON)_. IEEE, 2016, pp. 1–4.
  * [7] F. Capitanescu, L. F. Ochoa, H. Margossian, and N. D. Hatziargyriou, “Assessing the potential of network reconfiguration to improve distributed generation hosting capacity in active distribution systems,” _IEEE Transactions on Power Systems_ , vol. 30, no. 1, pp. 346–356, 2014.
  * [8] X. Geng, L. Tong, A. Bhattacharya, B. Mallick, and L. Xie, “Probabilistic hosting capacity analysis via bayesian optimization,” _arXiv preprint arXiv:2011.05193_ , 2020.
  * [9] J. Kisiala, “Conditional value-at-risk: Theory and applications,” _arXiv preprint arXiv:1511.00140_ , 2015.
  * [10] A. E. Lim, J. G. Shanthikumar, and G.-Y. Vahn, “Conditional value-at-risk in portfolio optimization: Coherent but fragile,” _Operations Research Letters_ , vol. 39, no. 3, pp. 163–171, 2011.
  * [11] F. Ding and B. Mather, “On distributed pv hosting capacity estimation, sensitivity study, and improvement,” _IEEE Transactions on Sustainable Energy_ , vol. 8, no. 3, pp. 1010–1020, 2016.
  * [12] M. S. S. Abad, J. Ma, D. Zhang, A. S. Ahmadyar, and H. Marzooghi, “Probabilistic assessment of hosting capacity in radial distribution systems,” _IEEE Transactions on Sustainable Energy_ , vol. 9, no. 4, pp. 1935–1947, 2018.
  * [13] S. F. Santos, D. Z. Fitiwi, M. Shafie-Khah, A. W. Bizuayehu, C. M. Cabrita, and J. P. Catalão, “New multistage and stochastic mathematical model for maximizing res hosting capacity—part i: Problem formulation,” _IEEE Transactions on Sustainable Energy_ , vol. 8, no. 1, pp. 304–319, 2016.
  * [14] R. T. Rockafellar and S. Uryasev, “Conditional value-at-risk for general loss distributions,” _Journal of banking & finance_, vol. 26, no. 7, pp. 1443–1471, 2002.
  * [15] S. Bruno, S. Ahmed, A. Shapiro, and A. Street, “Risk neutral and risk averse approaches to multistage renewable investment planning under uncertainty,” _European Journal of Operational Research_ , vol. 250, no. 3, pp. 979–989, 2016.
  * [16] A. N. Madavan, S. Bose, Y. Guo, and L. Tong, “Risk-sensitive security-constrained economic dispatch via critical region exploration,” in _2019 IEEE Power & Energy Society General Meeting (PESGM)_. IEEE, 2019, pp. 1–5.
  * [17] A. H. Arab, B. Taghezouit, K. Abdeladim, S. Semaoui, A. Razagui, A. Gherbi, S. Boulahchiche, and I. H. Mahammed, “Maximum power output performance modeling of solar photovoltaic modules,” _Energy Reports_ , vol. 6, pp. 680–686, 2020.
  * [18] M. E. Baran and F. F. Wu, “Optimal capacitor placement on radial distribution systems,” _IEEE Transactions on power Delivery_ , vol. 4, no. 1, pp. 725–734, 1989.
  * [19] S. Taheri, M. Jalali, V. Kekatos, and L. Tong, “Fast probabilistic hosting capacity analysis for active distribution systems,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 3, pp. 2000–2012, 2020.
  * [20] S. Bose, S. H. Low, T. Teeraratkul, and B. Hassibi, “Equivalent relaxations of optimal power flow,” _IEEE Transactions on Automatic Control_ , vol. 60, no. 3, pp. 729–742, 2014.
  * [21] M. Farivar and S. H. Low, “Branch flow model: Relaxations and convexification—part i,” _IEEE Transactions on Power Systems_ , vol. 28, no. 3, pp. 2554–2564, 2013.
  * [22] R. A. Jabr, “Radial distribution load flow using conic programming,” _IEEE transactions on power systems_ , vol. 21, no. 3, pp. 1458–1459, 2006\.
  * [23] F. Delbaen and S. Biagini, _Coherent risk measures_. Springer, 2000.
  * [24] L. Gan, N. Li, U. Topcu, and S. H. Low, “Exact convex relaxation of optimal power flow in radial networks,” _IEEE Transactions on Automatic Control_ , vol. 60, no. 1, pp. 72–87, 2014.
  * [25] A. M. Geoffrion, “Generalized benders decomposition,” _Journal of optimization theory and applications_ , vol. 10, no. 4, pp. 237–260, 1972.
  * [26] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [27] B. F. Lourenço, M. Muramatsu, and T. Tsuchiya, “Weak infeasibility in second order cone programming,” _Optimization Letters_ , vol. 10, no. 8, pp. 1743–1755, 2016.
  * [28] B. O’donoghue, E. Chu, N. Parikh, and S. Boyd, “Conic optimization via operator splitting and homogeneous self-dual embedding,” _Journal of Optimization Theory and Applications_ , vol. 169, no. 3, pp. 1042–1068, 2016.