1.2pt

# Multi-period Power System Risk Minimization under Wildfire Disruptions

Hanbin Yang, Noah Rhodes, Haoxiang Yang, Line Roald, Lewis Ntaimo Hanbin Yang
and Haoxiang Yang are supported by the National Natural Science Foundation of
China, under the grant numbers 72201232 and 72231008. Noah Rhodes and Line
Roald are supported by the U.S. National Science Foundation ASCENT program
under award 2132904 and by the Department of Energy, Office of Science, Office
of Advanced Scientific Computing Research, Applied Mathematics program under
Contract Number DE-AC02-06CH11347.

###### Abstract

Natural wildfire becomes increasingly frequent as climate change evolves,
posing a growing threat to power systems, while grid failures simultaneously
fuel the most destructive wildfires. Preemptive de-energization of grid
equipment is effective in mitigating grid-induced wildfires but may cause
significant power outages during natural wildfires. This paper proposes a
novel two-stage stochastic program for planning preemptive de-energization and
solves it via an enhanced Lagrangian cut decomposition algorithm. We model
wildfire events as stochastic disruptions with random magnitude and timing.
The stochastic program maximizes the electricity delivered while proactively
de-energizing components over multiple time periods to reduce wildfire risks.
We use a cellular automaton process to sample grid failure and wildfire
scenarios driven by realistic risk and environmental factors. We test our
method on an augmented version of the RTS-GLMC test case in Southern
California and compare it with four benchmark cases, including deterministic,
wait-and-see, and robust optimization formulations as well as a comparison
with prior wildfire risk optimization. Our method reduces wildfire damage
costs and load-shedding losses, and our nominal plan is robust against
uncertainty perturbation.

###### Index Terms:

optimal power flow, stochastic mixed-integer programming, wildfire risk, de-
energization, Lagrangian cut, decomposition algorithm.

##  I Introduction

Wildfire-related disasters have become more frequent and severe in recent
years across the United States [1, 2]. A major aspect of the wildfire research
focuses on wildfire mitigation plans in power systems [3]. There is a mutual
impact between wildfire ignition and power system operations: on the one hand,
many catastrophic fires with significant loss of life and property were
sparked by power infrastructures; for example, the Camp Fire in 2018, sparked
by a transmission line, caused $84$ lives deaths and an estimated $\$9.3$
billion in residential property damage alone [4]; on the other hand,
regardless of the source, wildfires lead to damages in power system
components, which change the topology and resource availability within a power
system [5].

Utility companies have considered approaches in different time scales to
reduce the risk of wildfire ignitions [3]. Long-term approaches can last for
years, such as equipment inspection, hardening, and vegetation management [6],
may require enormous costs and qualified personnel and should be complemented
with day-to-day operations such as public safety power shut-offs (PSPS) [7] to
eliminate risk. In this paper, we focus on optimizing the operational de-
energization decisions over a short time horizon similar to the setup in [8].

When optimizing the de-energization decisions, many previous works mainly
focus on wildfire risks caused by power system components but overlook the
simultaneous natural/man-made wildfires’ damage on power systems [9, 10, 11].
In addition, previous literature and industry practices take a deterministic
approach because it is computationally challenging to solve optimization under
uncertainty models with binary variables for de-energization, such as a
wildfire hazard prevention system based on data mining techniques [12], risk-
based optimization models for de-energization [8, 13], a data-driven wildfire
decision-making framework for power system resilience [14], and industry
practices using thresholds to indicate when to shut off [3, 15]. However,
incorporating wildfire spread dynamics can significantly improve the
performance of de-energization. Recently, more models have started to
incorporate stochastic aspects of wildfire risks, such as a dynamic
programming model to optimize a PSPS [16], a two-stage robust optimization
model for capacity expansion considering PSPS [17], a two-stage stochastic
integer program solved by progressive hedging [18], and rolling horizon
optimization models for joint de-energization and restoration [19, 20]. That
said, most of the literature still carries a relatively simple decision
structure and cannot fully capture the temporal dynamics of the wildfire; for
example, a commonly used uncertainty set in robust optimization, such as a
budgeted uncertainty set [21], can only capture the number of damaged
components, but is not able to accurately characterize the temporal effect of
wildfire progress within a multi-period horizon. Therefore, it requires new
studies to examine the interaction between power system operations with high-
fidelity and dynamic wildfire progress [22].

To holistically analyze the two-way impact between the wildfire ignitions and
the power system operations, we classify wildfires into two categories,
exogenous and endogenous fires: exogenous fires are external to the power
system, while power system faults trigger endogenous fires. Extreme weather
conditions such as drought [23] can lead to exogenous wildfires, which may
cause damages in power system components; e.g., dysfunctional shunts [24],
utility-scale converters [25], and line conductors [26]. Endogenous fires are
ignited when downed transmission lines come into contact with surrounding
vegetation and can damage nearby power system components [8]. A framework to
assess endogenous wildfire risk is proposed in [27], and simulation models
based on such risk measures can provide early warnings [28]. In this paper, we
explicitly model both exogenous and endogenous wildfires’ impact via binary
parameters indicating whether either wildfire damages each power system
component in every period.

Optimizing a sequential decision process that models shut-off and power flow
decisions requires solving a multistage stochastic program. This can be
challenging because the computational complexity increases exponentially with
the planning horizon length [29], and there may be integer variables involved
[8]. We develop a two-stage stochastic mixed-integer programming (SMIP) model
that captures the interaction between power system de-energization and both
exogenous and endogenous wildfires over multiple periods. We minimize our
wildfire risk measure, the expected total cost that incorporates the wildfire
damage cost and the load-shedding cost. In addition, our two-stage SMIP model
incorporates the spatial and temporal uncertainty of wildfires and represents
real-world situations with disruption scenarios, which can handle complex and
dynamic situations better and produce more resilient operations than
deterministic models. We use a cellular automaton process to simulate wildfire
spread and generate scenarios that describe the magnitude and timing of
wildfire disruptions. Since wildfire occurs infrequently, we assume at most
one wildfire occurrence during the planning horizon to reduce the multi-period
problem into two stages [30]. This two-stage scenario-based setting fits the
low-probability high-impact property of wildfire ignition and can capture the
dynamics of wildfire progress with a high fidelity. We design a proactive
scheme that executes an effective de-energization plan before the wildfire and
optimizes power flow operations of the remaining available components
afterward.

We propose a decomposition algorithm for the two-stage SMIP with non-convex
and non-smooth value functions. Our algorithm uses linear cutting planes,
named _Lagrangian cuts_ , to approximate the value functions [31]. The vanilla
version of Lagrangian cuts is tight theoretically but can be inefficient
computationally. To overcome this limitation, we introduce an enhanced version
of the Lagrangian cut, named square-minimization cut, which improves the
quality of the approximations and speeds up the algorithm. Our method shows
high computational accuracy and reliability. The two-stage SMIP model supports
the decision-making process under uncertainty, reduces the likelihood of a
disaster, and maintains as much load delivery as possible after a disruption.

The contributions of our paper are fourfold:

  1. 1.

We propose a cellular automaton model that simulates the onset and spread of
exogenous and endogenous wildfires based on environmental data. This enables
us to assess wildfire risk through scenarios used in the SMIP model.

  2. 2.

We formulate a two-stage SMIP model that manages wildfire risk in power system
operations. Unlike the existing literature, our model captures the random
nature of wildfire ignitions using disruption scenarios with random onset
times to keep the model size manageable.

  3. 3.

We develop a cutting plane algorithm that solves the SMIP model with a finite
convergence guarantee. We improve its computational performance significantly
with a new square-minimization cut.

  4. 4.

We show our two-stage SMIP model yields significantly more resilient solutions
than its deterministic counterparts and more economical solutions than the
robust optimization model. We also analyze the solution under different
parametric settings to gain operational insights.

The rest of the paper is organized as follows: Section II presents the
wildfire simulation model and the SMIP formulation; Section III derives the
cutting plane algorithm and its computational enhancement; Section IV covers
benchmark cases and numerical analyses; and Section V concludes the paper with
a brief discussion of future work.

##  II Two-Stage SMIP Model

We consider multi-period dispatch and de-energization operations under
wildfire disruption over a short-term time horizon ($24$ hours). The power
network is represented by a graph $({\cal B},{\cal L})$, where ${\cal B}$ is
the set of buses and ${\cal L}$ is the set of lines. We use $\mathcal{D}_{i}$,
$\mathcal{G}_{i}$, and ${\cal L}_{i}$ to represent the subset of loads,
generators, and transmission lines connected to bus $i$, respectively. For
each period $t\in\mathcal{T}$, we formulate DC power flow constraints to
decide the active power generation $P^{G}_{gt}$ of generator $g\in{\cal G}$,
the power flow $P^{L}_{ijt}$ on transmission line $(i,j)\in{\cal L}$ from bus
$i$ to bus $j$, and the phase angle $\theta_{it}$ of bus $i\in{\cal B}$.

In our model, we define wildfire risk as the combined cost of load-shedding
and wildfire damage under uncertain wildfire disruption over the specified
time horizon. In model (1), a de-energized generator will have zero generation
capacity, and a de-energized line will be considered open. If a bus is shut
off, all generators and lines connected to it will also be de-energized. Our
model uses binary decision variables $z_{ct}$ to represent whether a component
$c\in{\cal C}$ is de-energized at time period $t$. We assume that a de-
energized component will remain off until the end of the time horizon due to
safety considerations [8]. Model (1) obtains a first-stage plan, i.e., a
nominal plan, that should be implemented until a disruption occurs or the time
horizon ends, whichever occurs first. If a fire disruption is observed at
period $\tau^{\omega}$, given the current shut-off state
$z_{\cdot\tau^{\omega}-1}$, the model enters the second stage and incurs the
second-stage value function $f^{\omega}$, where we assume that the components
at revealed ignition locations will be shut off.

###  II-A Wildfire Scenario Modeling

On top of the constructed electric power network, we model wildfire
uncertainty as a stochastic disruption with random occurrence timing,
location, and spread. We consider a short-term operation and thus assume at
most one disruption will occur over the time horizon. We form a portfolio of
disruption scenarios indexed by $\omega\in\Omega$, within which we model i)
endogenous fires, which are caused by faults of operating electric components,
and ii) exogenous fires, which are caused by random natural wildfire
ignitions.

We consider the smallest rectangular area that can cover the entire power
network and divide it into a two-dimensional grid of a specified resolution,
in which the set of cells is denoted by $\mathcal{K}$. We then construct a
mapping between the grid cells and the buses and transmission lines based on
their geographic locations. We let $k_{i}\in\mathcal{K}$ denote the cell where
the bus $i$ is located. A transmission line, $l=(i,j)$, can span multiple
cells, which is denoted by $\mathcal{K}_{ij}$. Furthermore, a cell
$k\in\mathcal{K}$ can contain multiple transmission lines, and we denote the
set of transmission lines passing cell $k$ by ${\cal L}^{k}$. Finally, we let
$\mathcal{K}_{{\cal B}}$ and $\mathcal{K}_{{\cal L}}$ be the collection of
cells containing at least a bus and at least a transmission line,
respectively.

We employ the cellular automaton method, with propagation rules listed in
[32], to simulate the spread of exogenous wildfires. Overall, this setup leads
to a wildfire disruption in approximately $95\%$ of scenarios. Each cell is
characterized by four states which evolve in discrete time. The four states
are as follows: (i) cell has no forest fuel; (ii) cell has forest fuel that
has not ignited; (iii) cell has ignited forest fuel; (iv) cell has fully
burning forest fuel. Endogenous disruptions are caused by power system
component failures. The status of each component is normal, at fault, or
ignited. We define the onset time of a wildfire disruption $\tau^{\omega}$ as
the earlier time between an exogenous wildfire occurrence and any electric
component’s fault.

##### Exogenous Wildfires

We assume that there are no exogenous wildfires in the grid at the start of
the time horizon. For each cell $k\in\mathcal{K}$ in period $t\in{\cal T}$, it
receives information from the environment to update its state and interacts
with surrounding cells based on a set of defined rules:

  1. 1.

if its state in period $t-1$ is the unburned fuel state (ii), it has a
probability $p_{tk}$, influenced by environmental factors, to be ignited and
its state will transfer from (ii) to the ignited fuel state (iii);

  2. 2.

if its state in period $t-1$ is (iii), then the state will transfer to the
fully burning state (iv);

  3. 3.

if it is a fully burning cell, the state will stay at (iv) and can spread
flames to its $8$ adjacent cells indexed by $k^{\prime}$, igniting them with
probabilities $q_{tk^{\prime}}$ if they are in status (ii), where parameters
$q_{tk^{\prime}}$ are subject to environmental factors like vegetation, ground
elevation, and wind.

We record the exogenous wildfire damage after the end of time horizon ${\cal
T}$, letting $v_{c}^{\omega}=1$ if a wildfire (state (iii) or (iv)) ever
occurs on $k_{i}$ if $c$ is bus $i$ or any generator connected to bus $i$, or
on any cell in $\mathcal{K}_{ij}$, if $c$ is transmission line $(i,j)$.

##### Endogenous Wildfires

We simulate endogenous wildfires using a two-step procedure:

  1. 1.

First, we simulate electrical component faults. To initialize the fault
simulation process, we assume there is no fault component in the power system
at the start of the time horizon, similar to the exogenous wildfire
simulation. For every period $t\in{\cal T}$, a fault may occur at a component
$c\in{\cal C}$ following a Bernoulli distribution with probability $p_{tc}$.
Once a component faults, we record the state by letting $u_{c}^{\omega}=1$.

  2. 2.

Our second step runs an independent cellular automaton process for each fault
component $c$ from the fault period to the end of the time horizon,
initialized with the cell(s) that contains the component in state (iii). After
the cell evolves to state (iv) in the next period, it will spread the
simulated wildfire to the adjacent cell $k^{\prime}$ with a probability
$q_{t,k^{\prime}}$ at time period $t$ and further propagate. The components
damaged by this spread in the time horizon form a set $I_{c}^{\omega}$.

Fig. 1 provides an illustration of a set of components affected by an
endogenous fire: a fault occurs at an energized bus $i$, igniting an
endogenous wildfire and impacting two transmission lines, a bus, and a
generator (marked by dashed lines). Notice that we assume a fault directly
leads to a wildfire because we aim to simulate how an endogenous fire would
spread if the component $c$ were not de-energized. In addition, we build
independent simulation processes, each starting with only the cell(s)
containing the fault component $c$, in order to clearly trace the spread of a
wildfire to its origin and eliminate the interference of endogenous wildfires
started elsewhere.

Bus $j_{1}$Bus $i$Bus $j_{3}$$g_{1}$$g_{2}$$g_{3}$$I_{i}^{\omega}=\\{i,\
g_{2},\ (i,j_{1}),\ (i,j_{3}),\ j_{3}\\}$. fire spread direction Figure 1: A
random set $I_{i}^{\omega}$ illustration. A fault occurs at bus $i$
(energized), which causes an endogenous wildfire to spread and affect two
transmission lines $(i,j_{1})$ and $(i,j_{2})$, bus $j_{3}$, and generator
$g_{2}$, which are marked by blue dashed lines.

In summary, for a specific scenario $\omega\in\Omega$, we use a cellular
automaton procedure to generate the wildfire disruption uncertainty
$\xi^{\omega}=\\{\tau^{\omega},v_{\cdot}^{\omega},u_{\cdot}^{\omega},I_{\cdot}^{\omega}\\}$,
which includes the disruption time $\tau^{\omega}$, the binary parameters
$v_{c}^{\omega}$ indicating the exogenous fire locations, the binary
parameters $u_{c}^{\omega}$ indicating the fault locations, and the set
$I_{c}^{\omega}$ characterizing components damaged by the potential endogenous
fire at component $c$. We detail the parameter setup:

#### Mapping onto Grid

We convert the location of each bus $i\in{\cal B}$ from latitude and longitude
[33] to Universal Transverse Mercator(UTM) coordinates [34]. The smallest
rectangle that contains the entire power network has a length of $180,000$ and
a width of $138,000$ UTM unit. We divide such a rectangular area into square
cells with an edge length of $1,000$ UTM units. Such edge length ensures a
high resolution of our simulation model and is in the same order of magnitude
as the average distance that a wildfire can spread within a time period.

#### Probability of Exogenous Wildfire Ignition $p_{tk}$

We assume that only wildfires that occur in cells with at least one power
system component will threaten the electric power system. As stated in Section
II-A, we model the ignition of exogenous wildfires at time $t$ by a
probability $p_{tk}$, which describes the likelihood of a cell
$k\in\mathcal{K}$ transitioning from the state (ii), un-burned fuel, to state
(iii), ignited the fuel.To estimate $p_{tk}$, we use the Wildland Fire
Potential Index (WFPI) [35] as a measure of the relative probability of
ignition. The probability of an exogenous wildfire occurring in cell $k$ in
which the set of all transmission lines passing through cell $k$ is ${\cal
L}^{k}$ is

$p_{tk}=\dfrac{\sum_{l\in{\cal
L}^{k}}\texttt{WFPI}_{l}}{\sum_{l^{\prime}\in{\cal
L}}{\texttt{WFPI}_{l^{\prime}}}}.$

#### Probability of Wildfire Spread $q_{tk^{\prime}}$

Suppose there is a fully burning cell. Each of the adjacent cells
$k^{\prime}\in\mathcal{K}$ will be ignited with probability $q_{tk^{\prime}}$.
This probability $q_{tk^{\prime}}$ depends on: 1) a reference probability that
cell $k^{\prime}$ will be burning under no wind and flat terrain,
$q_{0k^{\prime}}$; 2) the vegetation type, $q^{\text{veg}}_{k^{\prime}}$, and
vegetation density, $q^{\text{den}}_{k^{\prime}}$; 3) geographical
information, $q^{s}_{k^{\prime}}$; and (iv) wind speed and direction,
$q^{w}_{tk^{\prime}}$, as follows:

$q_{tk^{\prime}}=q_{0k^{\prime}}\left(1+q^{\text{veg}}_{k^{\prime}}\right)\left(1+q^{\text{den}}_{k^{\prime}}\right)q^{s}_{k^{\prime}}q^{w}_{tk^{\prime}}.$

We refer to Ref. [32] and [36] for detailed parameter values.

#### Probability of Fault $p_{tc}$

Transmission lines are more prone to failures than other electrical components
due to environmental factors such as lightning and wind [37]. Therefore, in
our implementation, we focus on the fault at transmission lines and set
$p_{tc}=0,\forall c\in{\cal B}\cup{\cal G}$. This setup can be easily extended
if we want to consider the risk of fault at buses and generators. In any time
period, we assume a constant fault probability for each line $l\in{\cal L}$ as
$p_{tl}=1-e^{-\bar{\lambda}}$, where $\bar{\lambda}$ is the estimated hourly
failure rate in the Poisson regression model of [37].

###  II-B Two-Stage SMIP Model

When a fault occurs at an energized component, the associated endogenous
wildfire may incur high damage costs. To reduce the endogenous risk of
igniting a wildfire, the power system operator can de-energize, i.e., shut
off, power equipment. However, suppose multiple components are de-energized to
prevent endogenous fires. In that case, the power supply capacity may be
greatly reduced, resulting in the inability to meet crucial load demand and
causing serious secondary disasters. To improve the reliability of the power
system and mitigate wildfire damage, operators should de-energize some
potentially dangerous electrical equipment under high-risk conditions to
ensure enough power supply while significantly reducing risk.

To achieve a balanced risk trade-off between endogenous fire and load-
shedding, we build a two-stage SMIP model (1) which optimizes shut-off
decisions under wildfire disruption uncertainty over the given time horizon.
In model (1), a de-energized generator will have zero generation capacity, and
a de-energized line will be considered open. If a bus is shut off, all
generators and lines connected to it will also be de-energized. Our model uses
binary decision variables $z_{ct}$ to represent whether a component $c\in{\cal
C}$ is de-energized at time period $t$. We assume that a de-energized
component will remain off until the end of the time horizon due to safety
considerations [8]. Model (1) obtains a first-stage plan, i.e., a nominal
plan, that should be implemented until a disruption occurs or the time horizon
ends, whichever occurs first. If a fire disruption is observed at period
$\tau^{\omega}$, given the current shut-off state $z_{\cdot\tau^{\omega}-1}$,
the model enters the second stage and incurs the second-stage value function
$f^{\omega}$, where we assume that the components at revealed ignition
locations will be shut off.

$\displaystyle{Z^{*}=\min\
\sum_{\omega\in\Omega}p^{\omega}\left[\sum_{t=1}^{\tau^{\omega}-1}\sum_{d\in{\cal
D}}w_{d}(1-x_{dt})+f^{\omega}(z_{\cdot\tau^{\omega}-1},\xi^{\omega})\right]}$
$\displaystyle\mbox{s.t.}\quad\forall t\in{\cal T}:$ $\displaystyle
P^{L}_{ijt}\leq-
b_{ij}\left(\theta_{it}-\theta_{jt}+\bar{\theta}(1-z_{ijt})\right)\quad\quad\
\forall(i,j)\in{\cal L}$ (1a) $\displaystyle P^{L}_{ijt}\geq-
b_{ij}\left(\theta_{it}-\theta_{jt}+\underaccent{\bar}{\theta}(1-z_{ijt})\right)\quad\quad\
\forall(i,j)\in{\cal L}$ (1b) $\displaystyle-W_{ij}z_{ijt}\leq P^{L}_{ijt}\leq
W_{ij}z_{ijt}\quad\quad\quad\quad\quad\quad\forall(i,j)\in{\cal L}$ (1c)
$\displaystyle{\sum_{g\in{\cal G}_{i}}P^{G}_{gt}+\sum_{l\in{\cal
L}_{i}}P^{L}_{lt}=\sum_{d\in{\cal
D}_{i}}D_{dt}x_{dt}\qquad\quad\quad\quad\forall i\in{\cal B}}$ (1d)
$\displaystyle\underaccent{\bar}{P}^{G}_{g}z_{gt}\leq
P^{G}_{gt}\leq\bar{P}^{G}_{g}z_{gt}\quad\quad\quad\quad\quad\quad\quad\quad\quad\
\forall g\in{\cal G}$ (1e) $\displaystyle z_{it}\geq
x_{dt}\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\;\forall
i\in{\cal B},d\in{\cal D}_{i}$ (1f) $\displaystyle z_{it}\geq
z_{gt}\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\ \forall
i\in{\cal B},g\in{\cal G}_{i}$ (1g) $\displaystyle z_{it}\geq
z_{ijt}\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\forall
i\in{\cal B},(i,j)\in{\cal L}$ (1h) $\displaystyle z_{it}\geq
z_{jit}\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\forall
i\in{\cal B},(j,i)\in{\cal L}$ (1i) $\displaystyle z_{ct}\geq
z_{c\min\\{t+1,T\\}}\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\
\forall c\in{\cal C}$ (1j) $\displaystyle
z_{ct}\in\\{0,1\\}\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\
\forall c\in{\cal C}.$ (1k)

The objective function considers the expected cost, which consists of the
load-shedding cost before the disruption and the total cost after the
disruption. Constraints (1a)-(1e) correspond to the DC power flow model and
constraints (1f)-(1i) model the component interactions, equivalent to
constraint (7)-(9) in [8]. Constraints (1j) describe the temporal logic of
components’ status: once a component is shut off, it stays off, and the
corresponding $z$ variable is $0$. The scenario tree of this two-stage SMIP
model is shown as Fig. 2, in which each node represents the decisions at a
time period, and each branch extending to the right represents a disruption
scenario. The black box of the left diagonal branch represents the nominal
plan obtained if no disruption occurs. Suppose we have a scenario $\omega$
with $\tau^{\omega}=2$. The red box on the right describes the decisions
associated with such value function $f^{\omega}$, and the red node implies the
disruption time $\tau^{\omega}=2$. The branches in the blue box are the
scenarios that are disrupted in period $4$ and share the same historical
information (the nodes in blue).

Period$1$$2$$3$$4$$5$$..$$..$$\overset{\tau^{\omega}}{\rightsquigarrow}$
Figure 2: A scenario-tree illustration of the wildfire disruption problem with
$T=5$ and $\tau^{\omega}=2$. The branch in the red box represents the
disruptive scenarios where the disruption occurs in period $2$. The branches
in the blue boxes represent the scenarios sharing the same historical
information (the nodes in blue).

We define the expected total cost after the disruption
$\sum_{\omega\in\Omega}p^{\omega}f^{\omega}(\cdot,\xi^{\omega})$ as our
wildfire risk measure. For each scenario, the second-stage value function
$f^{\omega}$ takes in the shut-off state variables
$\hat{z}_{\cdot\tau^{\omega}-1}$ and the disruption uncertainty $\xi^{\omega}$
as its input. It focuses on the operations after $\tau^{\omega}$ and can be
evaluated by solving the following mixed-integer program:

$\displaystyle{f^{\omega}(\hat{z}_{\cdot\tau^{\omega}-1},\xi^{\omega})=}$
$\displaystyle\min\quad\sum_{t=\tau^{\omega}}^{T}\sum_{d\in{\cal
D}}w_{d}(1-x^{\omega}_{dt})+\sum_{c\in{\cal C}}r_{c}{\eta^{\omega}_{c}}$
$\displaystyle\mbox{s.t.}\quad\forall t\in\\{\tau^{\omega},\dots,T\\}:$
$\displaystyle P^{L,\omega}_{ijt}\leq-
b_{ij}\left(\theta^{\omega}_{it}-\theta^{\omega}_{jt}+\bar{\theta}(1-y^{\omega}_{ij})\right)\qquad\forall(i,j)\in{\cal
L}$ (2a) $\displaystyle P^{L,\omega}_{ijt}\geq-
b_{ij}\left(\theta^{\omega}_{it}-\theta^{\omega}_{jt}+\underaccent{\bar}{\theta}(1-y^{\omega}_{ij})\right)\qquad\forall(i,j)\in{\cal
L}$ (2b) $\displaystyle-W_{ij}y^{\omega}_{ij}\leq P^{L,\omega}_{ijt}\leq
W_{ij}y^{\omega}_{ij}\quad\quad\quad\quad\qquad\forall(i,j)\in{\cal L}$ (2c)
$\displaystyle\sum_{g\in{\cal G}_{i}}P^{G,\omega}_{gt}+\sum_{l\in{\cal
L}_{i}}P^{L,\omega}_{lt}=\sum_{d\in{\cal D}_{i}}D_{dt}x^{\omega}_{dt}\
\quad\qquad\forall i\in{\cal B}$ (2d)
$\displaystyle\underaccent{\bar}{P}^{G,\omega}_{g}y^{\omega}_{g}\leq
P^{G,\omega}_{gt}\leq\bar{P}^{G,\omega}_{g}y^{\omega}_{g}\quad\quad\quad\quad\quad\quad\quad\
\ \forall g\in{\cal G}$ (2e) $\displaystyle y^{\omega}_{i}\geq
x^{\omega}_{dt}\
\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\qquad\forall i\in{\cal
B},d\in{\cal D}_{i}$ (2f) $\displaystyle y^{\omega}_{i}\geq y^{\omega}_{g}\
\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\qquad\ \forall i\in{\cal
B},g\in{\cal G}_{i}$ (2g) $\displaystyle y^{\omega}_{i}\geq y^{\omega}_{ij}\ \
\ \quad\quad\quad\quad\quad\quad\quad\quad\qquad\forall i\in{\cal
B},(i,j)\in{\cal L}$ (2h) $\displaystyle y^{\omega}_{i}\geq y^{\omega}_{ji}\ \
\ \quad\quad\quad\quad\quad\quad\quad\quad\qquad\forall i\in{\cal
B},(j,i)\in{\cal L}$ (2i) $\displaystyle y^{\omega}_{c}\leq z^{\omega}_{c}\
\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\qquad\forall
c\in{\cal C}$ (2j) $\displaystyle y^{\omega}_{c}\leq 1-{\eta^{\omega}_{c}}\ \
\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\qquad\forall c\in{\cal
C}$ (2k) $\displaystyle{\eta^{\omega}_{c}}\geq v_{c}^{\omega}\
\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\qquad\forall
c\in{\cal C}$ (2l) $\displaystyle{\eta^{\omega}_{k}}\geq
u_{c}^{\omega}z^{\omega}_{c}\
\quad\quad\quad\quad\quad\quad\quad\quad\quad\qquad\forall c\in{\cal C},k\in
I^{\omega}_{c}$ (2m) $\displaystyle z^{\omega}_{c}=\hat{z}_{c\tau^{\omega}-1}\
\qquad\quad\qquad\qquad\qquad\qquad\qquad\ \forall c\in{\cal C}$ (2n)
$\displaystyle y^{\omega}_{c},{\eta^{\omega}_{c}},z^{\omega}_{c}\in\\{0,1\\}\
\ \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\ \forall c\in{\cal C}.$
(2o)

The objective function includes the load-shedding cost after the disruption
and the wildfire damage cost. The wildfire damage cost for a component
$c\in{\cal C}$, denoted by $r_{c}$, consists of the replacement cost of the
electric components and the financial loss to the nearby communities. In model
(2), the energization status will stay the same in the remaining time horizon,
as we assume that all wildfire damages reveal at period $\tau^{\omega}$ and no
recovery decisions take place afterward. Constraints (2a)-(2e) model an
equivalent form of the DC power flow constraints as in model (1). Constraints
(2f)-(2i) indicate the functioning state of components, similar to their
first-stage counterparts (1f)-(1i). We create a local copy of the first-stage
shut-off decisions, $z_{c}^{\omega}$, via the nonanticipativity constraint
(2n), which ensures that the shut-off state of each component at the time of
wildfire disruption is identical for all scenarios that share the same
history. The binding relationship is illustrated by Fig. 2 and discussed in
detail in [30]. With $z_{c}^{\omega}$ indicating whether component $c$ has
been shut off, constraint (2j) makes sure that the shut-off components stay
down through the second stage. Constraint (2k) states that damaged components
no longer function, where we model the exogenous fire damage by constraint
(2l) and the endogenous fire damage by constraint (2m). Notice that an
endogenous fire started at component $c$ requires both a fault occurrence
$u_{c}^{\omega}=1$ and the component not being shut off $z_{c}^{\omega}=1$ and
spreads to components $k\in I_{c}^{\omega}$.

##  III Decomposition Algorithm Based on Lagrangian Cut

The value function $f^{\omega}$ is a non-convex function without an analytical
form, as we need to solve a MIP subproblem (2) to evaluate it. A common way to
solve model (1) is to replace $f^{\omega}$ with closed-form approximation,
such as linear cutting planes. If the state variables are binary, we can
generate Lagrangian cuts to equivalently reformulate $f^{\omega}$ as a convex
piecewise linear function for each $\omega\in\Omega$ [31]. Such Lagrangian
cuts can be applied to model (1) as well because our state variables $z$ are
binary. In this section, we first show the derivation of the Lagrangian cut
and the convergence properties of the cutting plane method. In addition, we
propose an improved Lagrangian cut, _square-minimization cut_ (SMC), and
illustrate how it can accelerate the cutting plane algorithm.

At $\ell$-th iteration, we solve a Lagrangian dual problem (3) to obtain the
cut intercept $v^{\omega,\ell}$ as its optimal value $v^{*}$ and gradient
$\lambda^{\omega,\ell}$ as the optimal solution $\lambda^{*}$. It is worth
noting that model (3) is a convex program and we can efficiently solve it
using convex programming algorithms, such as the bundle method with
stabilization techniques [38]:

$v^{*}=\max_{\lambda}R^{\omega}(\hat{z},\lambda),$ (3)

where the Lagrangian relaxation problem $R^{\omega}(\hat{z},\lambda)$ is
obtained by relaxing the nonanticipativity constraint (2n):

$\displaystyle R^{\omega}(\hat{z},\lambda)=\min\quad$
$\displaystyle\sum_{t=\tau^{\omega}}^{T}\sum_{d\in\mathcal{D}}w_{d}(1-x_{dt}^{\omega})+$
$\displaystyle\sum_{c\in\mathcal{C}}[r_{c}{\eta^{\omega}_{c}}+{\lambda_{c}(\hat{z}_{c,\tau^{\omega}-1}-z_{c}^{\omega})}]$
(4) s.t. $\displaystyle\text{Constraints
}\eqref{eqn:pfconslt}-\eqref{eqn:ignitiont}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0},\eqref{eqn:binaryrestriction}}.$

Since we can always find a feasible solution to subproblem (2) by setting all
variables in model (2) to zero, we do not need to generate feasibility cuts to
characterize the domain of $f^{\omega}$ and can focus on generating optimality
cuts to approximate $f^{\omega}$ [29]. For a specific $\omega$, we denote the
lower approximation of $f^{\omega}$ by $V^{\omega}$, and write the $\ell$-th
cut as follows:

$V^{\omega}\geq
v^{\omega,\ell}+(\lambda^{\omega,\ell})^{\top}(z_{\tau^{\omega}-1}-\hat{z}^{\ell}_{\tau^{\omega}-1}).$
(5)

With construction of $L-1$ Lagrangian cuts (5), we can substitute the function
$f^{\omega}$ by to obtain the lower approximation for model (1) as $(M_{L})$.
Since the cuts form a lower approximation of $f^{\omega}$, the optimal value
$Z^{*}_{L}$ is a lower bound of $Z^{*}$.

$\displaystyle(M_{L})\quad Z_{L}^{*}=\min\quad$
$\displaystyle\sum_{\omega\in\Omega}{p^{\omega}}\left[\sum^{\tau^{\omega}-1}_{t=1}w_{d}(1-x_{dt})+V^{\omega}\right]$
(6) s.t. $\displaystyle\text{Constraints
}\eqref{eqn:pfconsl1}-\eqref{eqn:componenttime1}$ $\displaystyle
V^{\omega}\geq(\lambda^{\omega,\ell})^{\top}(z_{\cdot\tau^{\omega}-1}-\hat{z}^{\ell}_{\cdot\tau^{\omega}-1})+$
$\displaystyle\quad\quad
v^{\omega,\ell},\quad\forall\,\omega\in\Omega,\ell=1,\dots,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>
Once we obtain the optimal solution $(\hat{x},\hat{z},\hat{V},\hat{P})$ to
$(M_{L})$, it is a feasible solution to model (1) and we can calculate an
upper bound by evaluating
$f^{\omega}(\hat{z}_{\cdot\tau^{\omega}-1},\xi^{\omega})$ for each
$\omega\in\Omega$. We present the detailed steps of the decomposition
algorithm in Algorithm 1 below, which iteratively updates the bounds. In the
end, we terminate Algorithm 1 once the relative gap is within a predefined
tolerance threshold $\epsilon\geq 0$. We can show that Algorithm 1 converges
to the optimal value in finite steps.

1Initialization cut iteration number
$\ell={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1}$,
lower bound $LB=0$, upper bound $UB=\infty$ and $\epsilon\geq 0$;

2 while  _$\frac{UB-LB}{UB}\geq\epsilon$_ do

/* Forward Steps */

3 Solve the master problem $(M_{\ell})$ and obtain the first-stage shut-off
solution $\hat{z}^{\ell}$, optimal value $Z_{\ell}^{*}$, and the
approximations of value function $\hat{V}^{\omega,\ell}$ ;

4 Update $LB=Z_{\ell}^{*}$;

5 for  _$\omega\in\Omega$_ do

6 Solve model (2) with $\hat{z}^{\ell}$ and obtain
$f^{\omega}(\hat{z}^{\ell}_{\cdot\tau^{\omega}},\xi^{\omega})$;

7

8 end for

9 Let
$\bar{Z}=Z_{\ell}^{*}+\sum_{\omega\in\Omega}p^{\omega}(f^{\omega}(\hat{z}^{\ell}_{\cdot\tau^{\omega}},\xi^{\omega})-\hat{V}^{\omega,\ell}$)
;

10 if  _$\bar{Z} <UB$_ then

11 Update $UB=\bar{Z}$ and incumbent solution as $z^{*}=\hat{z}^{\ell}$;

12

13 end if

/* Backward Steps */

14 for  _$\omega\in\Omega$_ do

15 Solve the Lagrangian dual problem (3), and obtain the optimal solution
$\lambda^{\omega,\ell}$ and the optimal value
$v^{\omega,\ell}=R^{\omega}(\hat{z}^{\ell},\lambda^{\omega,\ell})$;

16 And add cut (5) to $(M_{\ell})$;

17

18 end for

19

20 Let $\ell=\ell+1$;

21

22 end while

Output the $\epsilon$-optimal value $UB$ and solution $z^{*}$.

Algorithm 1 Decomposition algorithm based on Lagrangian cuts to solve model
(1)

######  Theorem 1 (Convergence).

When $\epsilon=0$, Algorithm 1 terminates in a finite number of iterations and
outputs an optimal solution to model (1) after finitely many iterations $L$.

###### Proof.

According to Theorem $3$ in Ref. [31], for a given incumbent first-stage
solution $\hat{z}$ and scenario $\omega$, the Lagrangian cut generated at
$\hat{z}$ is proven to be tight, meaning that its value at $\hat{z}$ is equal
to the value function
$f^{\omega}(\hat{z}_{\cdot,\tau^{\omega}-1},\xi^{\omega})$ at this particular
point. Since all state variables are binary, there is a finite number of
distinct first-stage solutions $z$, allowing us to fully describe the value
function using a finite number of Lagrangian cuts. Consequently, there exists
a finite iteration number $k_{1}$, after which no new cuts are generated. For
any iteration $k>k_{1}$, with the resulting first-stage state variable
$\hat{z}^{k}$ and approximate second-stage value function values
$\\{V^{\omega,k}\\}_{\omega\in\Omega}$, two cases may arise: 1) for each
scenario $\omega$,
$V^{\omega,k}=f^{\omega}(\hat{z}^{k}_{\cdot,\tau^{\omega}-1},\xi^{\omega})$;
in this scenario, the resulting upper bound $\bar{Z}$ equals the optimal value
of the master problem $Z^{*}_{k}$, leading the algorithm to terminate with a
global optimum; 2) there exists a scenario $\omega$ such that
$V^{\omega}<f^{\omega}(\hat{z}^{k}_{\cdot,\tau^{\omega}-1},\xi^{\omega})$; in
this case, a Lagrangian cut is generated to enhance the approximation of the
second-stage value function. However, this contradicts the result that no new
cuts are generated after iteration $k_{1}$. Consequently, the algorithm
converges to a global optimum in a finite number of steps. ∎

Although the Lagrangian cut (5) is tight and valid at the incumbent binary
solution $\hat{z}$ and results in finite convergence, the cut may be steep and
not provide a good lower approximation at other solutions (shown as the black
line in Fig. 3). Therefore, We propose a rotated cut, the SMC, which is
illustrated by the red line in Fig. 3. Instead of solving model (3) to obtain
the cut coefficients $\lambda$ and $v$, we solve an alternative model:

$\displaystyle\min_{\lambda}\quad$ $\displaystyle\lambda^{\top}\lambda$ (7a)
s.t. $\displaystyle
R^{\omega}(\hat{z},\lambda)\geq(1-\delta)f^{\omega}(\hat{z}_{\cdot\tau^{\omega}-1},\xi^{\omega}),$
(7b)

and we let $\lambda^{\omega,\ell}$ equal to its optimal solution $\lambda^{*}$
and $v^{\omega,\ell}=R^{\omega}(\hat{z},\lambda^{*})$. Fig. 3 shows that the
difference between a steep cut and a flat cut lies in the angle they make with
the horizontal plane. We prefer a flat cut with a smaller
$\lambda^{\top}\lambda$, as $\lambda$ represents the linear cut’s
coefficients. While we minimize $\lambda^{\top}\lambda$ and rotate the cut, it
is possible that the cut becomes non-tight at the incumbent solution
$\hat{z}$. Therefore, we set up constraint (7b) to force the cut value to be
within a $\delta$ neighborhood of $f^{\omega}$ at $\hat{z}$, which can be
considered an “anchor point.” Model (7) is always feasible as the Lagrangian
cut coefficients obtained from solving model (3) can serve as a feasible
solution. As the function $R^{\omega}$ is a concave function of $\lambda$
given $\hat{z}$, constraint (7b) characterizes a convex set and we can again
use convex programming solution methods to solve model (7).

$0.2$$0.4$$0.6$$0.8$$1$$1.5$$2$$2.5$$3$$\hat{z}$$z$$f^{\omega}$non-tight
SMCtight & steep Figure 3: Illustration of different types of cuts: the black
cut and the red cut are tight at $\hat{z}=0$, but the blue cut is not since
its value at $\hat{z}$ is smaller than the function value; the red cut has
better performance than the black cut as it also characterizes the function
value exactly at another $\hat{z}=1$.

##  IV Numerical Results

In this section, we first describe the numerical experiment setup, with which
we evaluate the efficiency of our decomposition method and the solution
quality. We then compare our two-stage SMIP model (1) with four benchmark
cases: a deterministic model that ignores wildfire disruption, a wait-and-see
model with perfect wildfire information, a model based on an aggregate
wildfire risk metric, and a robust optimization model that seeks to minimize
the total cost of the worst-case scenario. Finally, we examine the robustness
of our nominal solution against potential inaccuracy in wildfire disruption
probability.

###  IV-A Experiment Setup

We use the RTS-GMLC $73$-bus case [33] and run experiments on a $24$-hour
horizon with each time period accounting for an hour ($T=24$). We assume a
constant electricity demand except for the peak hours, $9$ am to $12$ pm and
$3$ pm to $7$ pm, during which the demand is $1.2$ times the non-peak value.

The economic loss caused by fire and sudden power outages depends on the fire
intensity and the amount of load in the affected area. To quantify the
economic impact of wildfire damage and load loss, we assign each electrical
component a numerical rating based on its importance and impact on the area,
following a similar approach to the Value Response Index in Texas Wildfire
Risk Rating System [39]. We assign each electrical component a scaled cost
based on value data from [40], [41], [42], and [43] to create comparable cost
values for equipment damage and shed load. We obtain the following cost
parameters:

  1. 1.

Load priorities $w_{d}$ range from $50$ to $1000$;

  2. 2.

The damage costs $r_{c}$ of wind turbines, thermal and nuclear power plants
are $50$, $1000$, and $2500$, respectively;

  3. 3.

The damage cost $r_{c}$ of each bus is $50$;

  4. 4.

The damage cost $r_{c}$ of a transmission line is $0.285\ell$, where $\ell$ is
the length of the transmission line.

Fig. 4 color codes the load and generator priority levels based on their load-
shedding and damage costs, respectively.

Figure 4: Illustration of the RTS-GMLC system.

All optimization models were implemented using JuMP [44] in Julia v1.6 [45]
and solved by Gurobi v9.5.1 [46] on a computer with a 10-core M1 Pro CPU and
$32$GB memory. The network plots are generated using PowerPlot.jl, which
depends on PowerModels.jl package [47]. Scenario simulation is constructed by
an agent-based model package, Agent.jl [48]. For Algorithm 1 and SMC, we set
$\epsilon=1\%$ and $\delta=10^{-4}$.

###  IV-B Computational Performance and Solution Analysis

We examine how the algorithm performance and solution quality change with
sample size. Our two-stage SMIP model, utilizing a set of scenarios $\Omega$,
can be considered sample average approximation (SAA) and serve as a lower
bound estimator of the optimal value for the true stochastic process of
wildfire progress. As the number of scenarios increases, the approximation
precision improves, but the required computational effort also increases [29].
Our objective is to determine a sample size that yields a high-quality
solution within an acceptable computational budget. To achieve this, we first
obtain twenty first-stage SAA solutions
$\hat{X}=\\{\hat{x},\hat{z},\hat{\theta},\hat{P}\\}$ and optimal values with
sample sizes of $20,50,100,200$, and $500$.

Next, we generate $n=5,000$ wildfire scenarios $\tilde{\Omega}$ using the
procedure outlined in Section II-A for the out-of-sample test. We evaluate the
expected total cost of each SAA solution using those $5,000$ scenarios as
follows:

$\displaystyle g_{n}(\hat{X})=$ $\displaystyle\min\
\sum_{\omega\in\tilde{\Omega}}\frac{1}{n}\left[\sum_{t=1}^{\tau^{\omega}-1}\sum_{d\in{\cal
D}}w_{d}(1-\hat{x}_{dt})+f^{\omega}(\hat{z}_{\cdot\tau^{\omega}-1},\xi^{\omega})\right].$
(8)

Such cost can serve as an upper-bound estimator for the true wildfire
progress. We take the mean of those twenty lower- and upper-bound estimators
and calculate their $95\%$ confidence intervals for each sample size, shown in
Fig. 5. As the figure illustrates, the gap between the lower- and upper-bound
estimators narrows, and both estimators become less variable with increasing
sample size. We consider this gap acceptable for a sample size of 500.

As the sample size increases, solving the SMIP model (1) becomes increasingly
challenging. Table I shows the effectiveness of Algorithm 1 with SMC. With a
sample size of $500$, this algorithm is capable of solving the model within a
reasonable time frame ($18,000$ seconds). In contrast, state-of-the-art
solvers like Gurobi are unable to handle such a problem. Table I highlights
the necessity of computational improvement through the use of SMC, as the
algorithm with the original Lagrangian cuts (LC) fails to reach the optimal
solution within the given time constraint.

We simulate scenarios and solve the two-stage SMIP offline in advance of a
specific day of operations, which provides a nominal plan. During the day, we
execute the nominal plan until a wildfire disruption occurs. Given the
wildfire disruption realization, we can solve the second-stage problem to
obtain recourse decisions, a process that typically takes less than a second.
Therefore, the run-time of Algorithm 1 is sufficient for the practical
implementation of our model.

Figure 5: Confidence intervals (bar) and point estimates (circle) of the lower and upper bounds for solutions with different sample sizes. The shaded area indicates the range of extreme values (maximum and minimum). TABLE I: Run-time (sec.) of Algorithm 1 utilizing SMC and LC with a tolerance of $1.0\%$, compared with a run-time of solving model (1) by Gurobi with MIPGap $=1\%$. Sample size | Algorithm 1 Time | Gurobi Time
---|---|---
SMC (Sec.) | LC (Sec.) | (Sec.)
$20$ | $712$ | $3287$ | $84$
$50$ | $1248$ | $8741$ | $376$
$100$ | $2604$ | $17911$ | $1510$
$200$ | $5323$ | $>18000$ | $9186$
$500$ | $12776$ | $>18000$ | $>18000$

The nominal plan obtained by model (1) balances the trade-off between
minimizing load-shedding cost and reducing wildfire risk. For example,
transmission line $72$ has a high risk of starting an endogenous wildfire, as
it ranks $10$-th among all transmission lines regarding the number of
scenarios where an endogenous wildfire occurs. However, in a $24$ time-period
horizon, the nominal plan does not de-energize it early because it connects to
a bus with multiple generators that supply a large amount of load; see Fig.
6(a). To illustrate this trade-off, we compare the outcomes of de-energizing
and keeping line $72$ energized. If we intentionally de-energize line $72$
early at $t=12$, it causes a significant loss of generation and load-shed; see
Fig. 6(b). The load-shedding cost outweighs the expected damage cost caused by
line $72$, so the nominal plan prefers to keep it energized. However, after
period $22$, the nominal plan de-energizes it (Fig. 6(c)) because the expected
damage cost from its fault exceeds the load-shedding cost from less
generation. We observe a similar “end-of-horizon” effect when we extend the
time horizon length to $T=36$ and $T=48$.

(a) The nominal load plan for the period $12$ and its associated load-shedding
percentage illustrates how the nominal plan optimizes the load distribution
and minimizes the risk exposure under uncertainty.

(b) The additional load shed in period $12$ when transmission line $72$ is de-
energized intentionally, compared to when the nominal plan keeps it on.

(c) The nominal load plan for period $22$, illustrates how the nominal plan
optimizes the load distribution and minimizes the risk exposure under
uncertainty.

Figure 6: Visualization of trade-off in the nominal plan. The red line in the
bottom-right corner represents transmission line $72$.

###  IV-C Benchmark Comparison

We choose the best optimal solution obtained from twenty batches with a sample
size of $500$ in Section IV-B as the solution to our two-stage SMIP model (1),
denoted as $X^{*}=\\{x^{*},z^{*},\theta^{*},P^{*}\\}$. We then compare this
optimal solution, $X^{*}$, to the following three benchmark models:

  1. 1.

A deterministic model where we maximize load satisfaction without considering
possible disruption:

$\displaystyle(M^{det})\;\;\min\quad$
$\displaystyle\sum_{t=1}^{T}\sum_{d\in{\cal D}}w_{d}(1-x_{dt})$ s.t.
$\displaystyle\text{Constraints
}\eqref{eqn:pfconsl1}-\eqref{eqn:componenttime1}.$

Since $(M^{det})$ does not consider the wildfire risk, the optimal solution
will return a nominal plan
$X^{det}=\\{x^{det},z^{det},\theta^{det},P^{det}\\}$ with no shut-off.

  2. 2.

A wait-and-see model where we can make a scenario-specific plan assuming
$\xi^{\omega}$ is known before we make decisions. For each $\omega\in\Omega$,
we solve

$\displaystyle(M^{ws,\omega})\qquad$$\displaystyle
g^{ws,\omega}=\qquad\qquad\qquad\qquad$ $\displaystyle\min\quad$
$\displaystyle\sum_{t=1}^{\tau^{\omega}-1}\sum_{d\in{\cal
D}}w_{d}(1-x_{dt})+f^{\omega}(z_{\cdot\tau^{\omega}-1},\xi^{\omega})$ s.t.
$\displaystyle\text{Constraints
}\eqref{eqn:pfconsl1}-\eqref{eqn:componenttime1},$

and obtain a scenario-specific solution
$X^{ws,\omega}=\\{x^{ws,\omega},z^{ws,\omega},\theta^{ws,\omega},P^{ws,\omega}\\}$.
An expected wait-and-see cost can then be computed as
$g_{n}^{*}=\sum_{\omega\in\Omega}{\frac{1}{n}}g^{ws,\omega}$.

  3. 3.

A wildfire-risk-based model proposed by [19] to simultaneously reduce wildfire
risk and load-shedding:

$\displaystyle(M^{rb})\;\;\min\quad$
$\displaystyle\sum_{t=1}^{T}\left[\dfrac{\alpha
R_{Fire,t}}{R_{Tot}}-(1-\alpha)\dfrac{\sum_{d\in{\cal
D}}x_{dt}w_{d}D_{dt}}{D_{Tot}}\right]$ s.t. $\displaystyle\mbox{Constraints
}\eqref{eqn:pfconsl1}-\eqref{eqn:componenttime1}$ $\displaystyle
R_{Fire,t}=\sum_{c\in{\cal C}}R_{ct}z_{ct}\quad\quad\forall t\in{\cal T},$

The component-wise risk measure $R_{ct}$ is the mean number of components
affected by an endogenous fire started at component $c$ across all scenarios.
Fig. 7 illustrates the value of $R_{ct}$ for $t=13$. The components that face
significant wildfire risk are marked in the figure. By controlling the number
of energized components, we can compute the total wildfire risk $R_{Fire,t}$
for the period $t$. We use the total load $D_{Tot}=\sum_{d\in{\cal
D},t\in{\cal T}}D_{dt}$ and total wildfire risk $R_{Tot}=\sum_{c\in{\cal
C},t\in{\cal T}}R_{ct}$ as the normalization factors such that we can combine
the wildfire risk and the load-shedding in the objective function. The
parameter $\alpha\in[0,1]$ adjusts the trade-off between serving loads and
reducing wildfire risk. Solving model $(M^{rb})$ returns a nominal plan
$X^{rb,\alpha}=\\{x^{rb,\alpha},z^{rb,\alpha},\theta^{rb,\alpha},P^{rb,\alpha}\\}$.

Figure 7: Illustration of wildfire risk $R_{ct}$ at $t=13$.

  4. 4.

A scenario-based robust optimization model, in which we minimize the total
cost of the worst-case scenario:

$\displaystyle(M^{ro})\ \min$ $\displaystyle\ \max_{\omega\in\Omega}\
\sum_{t=1}^{\tau^{\omega}-1}\sum_{d\in{\cal
D}}w_{d}(1-x_{dt})+f^{\omega}(z_{\cdot\tau^{\omega}-1},\xi^{\omega})$ s.t.
$\displaystyle\
\eqref{eqn:pfconsl1}-~{}\eqref{eqn:binaryrestriction1}\qquad\qquad\qquad\qquad\qquad\forall
t\in{\cal T}.$

Solving model $(M^{ro})$ returns a nominal plan
$X^{ro}=\\{x^{ro},z^{ro},\theta^{ro},P^{ro}\\}$.

Each model is solved, then evaluated with the same $n=5,000$ out-of-sample
scenarios $\tilde{\Omega}$ used to evaluate the two-stage SMIP model (1). For
the deterministic model and the wait-and-see model, we obtain two following
ratios to show their relative difference compared to model (1):

$\displaystyle\dfrac{g_{n}(X^{det})-g_{n}(X^{*})}{g_{n}(X^{*})}=0.639\mbox{
and }\dfrac{g_{n}(X^{*})-g_{n}^{*}}{g_{n}(X^{*})}=0.416.$

The first ratio emphasizes that the two-stage SMIP can significantly improve
the operations compared to the deterministic approach agnostic to potential
wildfires. The second ratio quantifies the value of perfect information. We
can save $41.6\%$ of the cost if we can forecast the wildfire perfectly before
we make the shut-off decisions. This large ratio reflects the high variance in
wildfire scenarios and the difficulty for a single nominal shut-off plan to
perform optimally in all scenarios. While it is not entirely possible to
obtain the perfect information, it is helpful if we can obtain auxiliary data,
conditioned on which we can identify a better representation of the scenarios.
For example, we may exclude the scenarios with wildfires occurring in areas
with forecast precipitation.

Fig. 8 illustrates such effect by selecting a subset of twenty scenarios: we
can observe that the cost of the SMIP solution is significantly lower than the
cost of the deterministic solution in all but one of the selected scenarios
and can approach the best possible cost in two scenarios. If the auxiliary
data happen to show those two scenarios are likely to happen, our SMIP
solution will perform well. On the contrary, if the auxiliary data tell the
other way, we may want to re-solve the two-stage SMIP model 1 with a refined
set of scenarios.

Figure 8: Scenario-specific costs for model (1), the deterministic model
$(M^{det})$, and the wait-and-see model $(M^{ws,\omega})$. For each scenario,
the green star and orange circle represent the costs incurred by the nominal
plans generated by $(M^{det})$ and model (1). The blue star marks the wait-
and-see cost.

For the wildfire-risk-based model $(M^{rb})$, we use the set of scenarios
$\Omega^{*}$, with which we obtain $X^{*}$, to calculate the wildfire risk
value $R_{Fire,t}$ for all $t\in{\cal T}$. We obtain a solution
${X}^{rb,\alpha}$ for $\alpha\in\\{0.0,0.1,\dots,0.9\\}$. We list the out-of-
sample mean load-shedding cost and wildfire damage cost under non-disruptive
scenarios (no wildfire ignitions) and disruptive scenarios in Table II,
evaluated with $\tilde{\Omega}$.

If the wildfire risk is ignored ($\alpha=0.0$), it is equivalent to the
deterministic solution $X^{det}$, which leads to a large out-of-sample cost
$g_{n}(X^{det})$. By increasing $\alpha$, we can first reduce the total cost
via the preventive shut-offs of potentially risky components, resulting in
less damage cost. However, as $\alpha$ increases further, the nominal plan
becomes too conservative, and many components are de-energized such that the
increasing load-shedding cost outweighs the reduced damage costs. This
illustrates a major challenge in the wildfire-risk-based model $(M^{rb})$ as
there is little guidance to select a proper $\alpha$.

We then compare our SMIP model (1) with the wildfire-risk-based model
$(M^{rb})$, via the ratio of relative improvements (RRI)
$(g_{n}({X}^{rb,\alpha})-g_{n}(X^{*}))/g_{n}(X^{*})$ in Table II. The results
show that the nominal plan obtained by solving model (1) performs
significantly better, and the RRI values for $X^{rb,\alpha}$ exceed $45\%$ for
all $\alpha$ selected. This indicates that, without considering wildfire
scenarios and potential adaptive recourse actions, static non-robust planning
is far from optimum, regardless of the level of conservatism. The reason is
that $(M^{rb})$ only considers minimizing the total wildfire risk and
satisfying the power demand, which leads to the nominal plans that de-energize
any component with $R_{ct}>0$ for some $t\in{\cal T}$ as early as possible,
only if this action does not affect the future supply of electricity. However,
this property may compromise the resilience of the power grid by de-energizing
an excessive number of components. If a disruption occurs, there might be very
high load-shedding costs.

TABLE II: Load-shedding costs and damage costs of nominal plans $X^{rb,\alpha}$ under non-disruptive and disruptive scenarios.  | Non-Disr. | Disruptive | Total Cost
---|---|---|---
$\alpha$ | Load shed | Load shed | Damage | $g_{n}(\cdot)$ | RRI
$0.0/X^{det}$ | $0.0$ | $3247.0$ | $2814.0$ | $5936.2$ | $63.9\%$
$0.1$ | $148.6$ | $3843.5$ | $1767.5$ | $5495.4$ | $51.7\%$
$0.2$ | $162.9$ | $3850.7$ | $1739.8$ | $5480.8$ | $51.4\%$
$0.3$ | $221.4$ | $3891.0$ | $1713.9$ | $5489.4$ | $51.6\%$
$0.4$ | $171.6$ | $3786.3$ | $1783.4$ | $5454.9$ | $50.6\%$
$0.5$ | $251.1$ | $3712.8$ | $1773.9$ | $5373.7$ | $48.4\%$
$0.6$ | $106.7$ | $3686.8$ | $1817.0$ | $5390.5$ | $48.9\%$
$0.7$ | $469.3$ | $3721.4$ | $1754.1$ | $5425.9$ | $50.0\%$
$0.8$ | $883.6$ | $3949.0$ | $1700.1$ | $5532.7$ | $52.8\%$
$0.9$ | $1505.8$ | $4285.3$ | $1341.1$ | $5510.5$ | $52.2\%$
$X^{*}$ | $1678.2$ | $2355.2$ | $1333.6$ | $3612.7$ | $0.0\%$

We use a scenario-appending algorithm similar to that in Ref. [49]. The
algorithm starts with solving a relaxed model considering only a prespecified
subset of scenarios $\hat{\Omega}\subset\Omega^{*}$ and obtains a solution
$\hat{X}$. An oracle is run to find the worst-case scenario
$\hat{\omega}\in\Omega^{*}$ given $\hat{X}$. We append such a scenario to the
set $\hat{\Omega}$ and repeat the above process until the obtained
$\hat{\omega}$ is already in $\hat{\Omega}$. The experiment demonstrates that
the subset $\hat{\Omega}$ includes $23$ scenarios, much smaller than
$|\Omega^{*}|=500$, allowing us to run the algorithm in $109$ seconds.
Although the robust optimization algorithm runs faster, the robust solution
$X^{ro}$ presents worse out-of-sample test results using scenario set
$\tilde{\Omega}$ in TABLE III, compared with those of the SMIP solution
$X^{*}$. The robust nominal plan $X^{ro}$ de-energizes more components early
in the time horizon to mitigate endogenous wildfire risks, which is
illustrated in Fig. 9. Consequently, this conservative strategy results in a
substantial amount of load-shedding in both nominal and disruptive scenarios,
leading to a total cost of up to four times that of $X^{*}$. The worst-case
in-sample cost among scenarios in $\Omega^{*}$ for $X^{ro}$ is $32,374$
compared to $34,016$ for $X^{*}$, which is reasonable as $X^{ro}$ optimizes
the worst-case scenario. However, the out-of-sample worst-case costs, as
demonstrated in the last column of TABLE III, indicate that $X^{ro}$ is
sensitive to the scenario selection and performs poorly when facing new
samples. In contrast, our stochastic programming solution $X^{*}$ shows
robustness as the out-of-sample worst-case cost is even smaller than the in-
sample one. Moreover, the performance of $X^{ro}$ significantly lags behind
that of $X^{*}$ across other scenarios. Although the robust model requires
less time to generate a feasible nominal plan, the numerical results show that
the robust nominal plan is conservative because it relies on limited scenario
information and prioritizes minimizing worst-case costs.

TABLE III: Load-shedding costs and damage costs of nominal plans $X^{ro}$ under non-disruptive and disruptive scenarios.  | Non-Disr. | Disruptive | Cost
---|---|---|---
| Load shed | Load shed | Damage | $g_{n}(\cdot)$ | Worst-case
$X^{ro}$ | $30130.2$ | $18067.3$ | $2506.3$ | $20770.5$ | $56848.3$
$X^{*}$ | $1678.2$ | $2355.2$ | $1333.6$ | $3612.7$ | $28871.1$
Figure 9: Illustration of the power system under the nominal plan $X^{ro}$ at
time $t=16$.

###  IV-D Sensitivity Analysis of Disruption Probability

We simulate the set of wildfire scenarios with a probabilistic model described
in Section II-A, which then serves as the uncertainty model in the SMIP. It is
necessary to examine the robustness of the optimal solution to the SMIP
against the potential inaccuracy in such an uncertainty model. Therefore, we
perform a sensitivity analysis in this section to check the out-of-sample
performance of SMIP solutions obtained with different uncertainty assumptions.

In Section II, we construct a scenario set $\Omega^{*}$ in the SMIP model (1),
each scenario with an equal probability of occurrence, $1/|\Omega^{*}|$. Such
a scenario set can be partitioned as
$\Omega^{*}=\Omega^{*}_{D}\cup\Omega^{*}_{N}$ , where $\Omega^{*}_{D}$
$(\Omega^{*}_{N})$ represents the set of scenarios with(out) a disruption.
When there is no disruption, the nominal plan is carried out over the entire
time horizon. Therefore, the partition is equivalent to
$\Omega^{*}=\Omega^{*}_{D}\cup\\{{\omega_{0}}\\}$ with probabilities
$p^{\omega_{0}}=|\Omega^{*}_{N}|/|\Omega^{*}|$ and
$p^{\omega}=(1-p^{\omega_{0}})/|\Omega^{*}_{D}|$ for
$\omega\in\Omega^{*}_{D}$, where ${\omega_{0}}$ is a scenario without
disruption.

For the sensitivity analysis, we add a perturbation $\Delta p$ to
$p^{\omega_{0}}$ and adjust $p^{\omega}$ for $\omega\in\Omega^{*}_{D}$
accordingly. We solve the SMIP for each $\Delta p$ and test its solution’s
out-of-sample performance over the scenario set $\tilde{\Omega}$, displayed in
Table IV. The out-of-sample no-disruption probability
$|\tilde{\Omega}_{N}|/|\tilde{\Omega}|$ is approximately $0.05$, similar to
that of $\Omega^{*}$.

We observe from Table IV that the optimal SMIP solution can yield a decent
out-of-sample performance even when the input perturbation is quadruple the
original no-disruption probability ($\Delta p=0.2$). This indicates that the
SMIP optimal solution is robust against input inaccuracy. When $\Delta p$
keeps increasing, the SMIP values a nominal plan that shuts off fewer
components and thus minimizes the load-shedding cost under scenarios without
wildfire disruption. It leads to higher out-of-sample load-shedding and damage
costs simultaneously as insufficient shut-offs are executed to prevent
endogenous wildfires. The sensitivity analysis results show that the SMIP
optimal solution is reliable as long as the input uncertainty model is not too
liberal when estimating the wildfire risk.

TABLE IV: Load-shedding costs and damage costs of nominal plans $X^{\Delta p}$ under non-disruptive and disruptive scenarios.  | Non-Disr. | Disruptive | Total Cost
---|---|---|---
$\Delta p$ | Load shed | Load shed | Damage | $g_{n}(\cdot)$ | RRI
$-0.05$ | $1805.1$ | $2553.0$ | $1242.7$ | $3782.8$ | $4.5\%$
$-0.015$ | $1741.5$ | $2513.7$ | $1353.6$ | $3787.3$ | $4.6\%$
$0.0/X^{*}$ | $1678.2$ | $2355.2$ | $1333.6$ | $3612.7$ | $0.0\%$
$0.015$ | $1660.0$ | $2343.9$ | $1341.7$ | $3609.7$ | $-0.1\%$
$0.05$ | $1504.9$ | $2360.2$ | $1325.4$ | $3609.6$ | $-0.1\%$
$0.1$ | $1250.7$ | $2344.8$ | $1317.1$ | $3586.5$ | $-0.7\%$
$0.2$ | $828.4$ | $2196.5$ | $1578.8$ | $3697.5$ | $2.3\%$
$0.3$ | $207.3$ | $2181.1$ | $1752.3$ | $3882.5$ | $7.3\%$
$0.5$ | $185.1$ | $2343.5$ | $1845.6$ | $4109.1$ | $13.5\%$
$0.7$ | $109.7$ | $2410.8$ | $1851.0$ | $4153.3$ | $15.0\%$
$0.9$ | $99.0$ | $2476.8$ | $1878.1$ | $4276.7$ | $18.4\%$

###  IV-E Interaction between Endogenous and Exogenous Wildfires

In this section, we investigate how endogenous and exogenous wildfires
interact and affect the shut-off plans. Intuitively speaking, considering
endogenous wildfires tends to encourage more shut-offs to prevent faults
ignited by power system components, while the presence of exogenous wildfires
introduces additional load-shedding costs which can be mitigated by keeping
more components energized. To quantify this interaction, we compare three
nominal plans: 1) $X^{exo}$ obtained by solving model (1) with $500$ scenarios
that only contain exogenous wildfires; 2) $X^{end}$ obtained by solving model
(1) with $500$ scenarios that only contain endogenous wildfires; and 3)
$X^{*}$ obtained in Section IV-C with $500$ scenarios that contain a mixture
of exogenous and endogenous wildfires. Fig. 10 illustrates the power system
status, highlighting the de-energized components and load-shed in period
$t=16$. The nominal plan $X^{exo}$, shown in the left plot, exhibits the
fewest number of de-energized components. In contrast, the right plot shows
the nominal plan $X^{end}$ leads to the highest number of shut-offs. The
nominal plan $X^{*}$, considering both exogenous and endogenous wildfires,
represents a trade-off between the two. We perform three additional out-of-
sample tests to evaluate those nominal plans with $5,000$ scenarios of: 1)
exogenous wildfires only; 2) endogenous wildfires only; and 3) a mixture of
endogenous and exogenous wildfires, as shown in Fig. 10. The nominal plan
$X^{exo}$ exhibits inferior performance in both load-shedding and damage costs
when endogenous wildfires are possible. The nominal plan $X^{end}$ becomes
conservative and shuts off too many components, resulting in consistently high
load-shedding costs. While the nominal plan $X^{*}$ achieves a balance in
between. Our finding confirms the intuition of how exogenous and endogenous
wildfires interact and highlights the importance of modeling both types of
wildfires in power system operations.

(a) $15$ transmission lines are de-energized in $X^{exo}$ in period $16$.

(b) $26$ transmission lines are de-energized in $X^{*}$ in period $16$.

(c) $31$ transmission lines are de-energized in $X^{end}$ in period $16$.

Figure 10: Visualization of interaction between endogenous and exogenous wildfire risk. TABLE V: The performance of the nominal plans generated from different scenarios under different test sets.  | Exogenous Wildfire Scenario Set | Mixture Wildfire Scenario Set | Endogenous Wildfire Scenario Set
---|---|---|---
| Non-Disr. | Disruptive | Total | Non-Disr. | Disruptive | Total | Non-Disr. | Disruptive | Total
Type | Load shed | Load shed | Damage | $g_{n}(\cdot)$ | Load shed | Load shed | Damage | $g_{n}(\cdot)$ | Load shed | Load shed | Damage | $g_{n}(\cdot)$
$X^{exo}$ | $0.0$ | $707.1$ | $636.5$ | $1289.8$ | $0.0$ | $3231.8$ | $2961.0$ | $5960.2$ | $0.0$ | $3380.4$ | $2979.8$ | $6106.0$
$X^{*\quad}$ | $1678.1$ | $1892.8$ | $636.5$ | $2495.2$ | $1678.2$ | $2355.2$ | $1333.6$ | $3612.7$ | $1678.2$ | $3107.1$ | $1814.0$ | $4792.8$
$X^{end}$ | $1998.5$ | $2249.9$ | $636.5$ | $2850.9$ | $1998.5$ | $2748.6$ | $1240.3$ | $3909.1$ | $1998.5$ | $3147.0$ | $1514.0$ | $4554.4$

##  V Conclusions & Future Work

This work contributes to the literature on stochastic optimization for power
system operation under wildfire risk. We propose a novel two-stage stochastic
program that captures the uncertainty of wildfire disruptions and temporal
dynamics. We also develop an efficient decomposition algorithm that exploits
the binary structure of state variables and generates valid cuts for solving
large-scale instances. Our numerical experiments demonstrate the effectiveness
and robustness of our approach in terms of solution quality and computational
speed. They illustrate the benefits of the stochastic programming model by
comparing it with deterministic counterparts and the scenario-based robust
optimization model.

For future work, we plan to extend our research in several directions. One
direction is to investigate a theory-driven approach for generating effective
and efficient cuts in our decomposition algorithm. Another direction is to
expand the current modeling setting: 1) to relax the assumption of at most one
disruption and extend our model to a multistage problem; 2) to include the
unit commitment decisions, which would require an algorithmic development to
handle the continuous state variables to model the ramping constraints; 3) to
formulate a realistic multi-period uncertainty set/ambiguity set based on
high-fidelity wildfire simulators and model the wildfire disruption in
robust/distributionally robust optimization. A third direction is to address
the challenge of ensuring feasibility with AC power flow when using a linear
approximation of power flow in our model. Finally, we aim to develop a rolling
horizon framework that can handle the wildfire forecast updates and
dynamically re-optimize the de-energization schedule.

## References

  * [1] J. W. Muhs, M. Parvania, and M. Shahidehpour, “Wildfire risk mitigation: A paradigm shift in power systems planning and operation,” _IEEE Open Access Journal of Power and Energy_ , vol. 7, pp. 366–375, 2020.
  * [2] T. P. Holmes, R. J. Huggett Jr, and A. L. Westerling, “Statistical analysis of large wildfires,” _The Economics of Forest Disturbances_ , vol. 79, pp. 59–77, 2008.
  * [3] D. A. Z. Vazquez, F. Qiu, N. Fan, and K. Sharp, “Wildfire mitigation plans in power systems: A literature review,” _IEEE Transactions on Power Systems_ , vol. 37, no. 5, pp. 3540–3551, 2022.
  * [4] “2019 wildfire risk report,” CoreLogic, Tech. Rep., 2019. [Online]. Available: https://www.corelogic.com/wp-content/uploads/sites/4/downloadable-docs/wildfire-report_0919-01-screen.pdf
  * [5] M. Choobineh, B. Ansari, and S. Mohagheghi, “Vulnerability assessment of the power grid against progressing wildfires,” _Fire Safety Journal_ , vol. 73, pp. 20–28, 2015.
  * [6] S. Jazebi, F. de León, and A. Nelson, “Review of wildfire management techniques—Part I: Causes, prevention, detection, suppression, and data analytics,” _IEEE Transactions on Power Delivery_ , vol. 35, no. 1, pp. 430–439, 2020.
  * [7] “Utility public safety power shutoff plans (de-energization).” California Public Utilities Commission (CPUC), 2022. [Online]. Available: https://www.cpuc.ca.gov/psps/
  * [8] N. Rhodes, L. Ntaimo, and L. Roald, “Balancing wildfire risk and power outages through optimized power shut-offs,” _IEEE Transactions on Power Systems_ , vol. 36, no. 4, pp. 3118–3128, 2021.
  * [9] S. Taylor and L. A. Roald, “A framework for risk assessment and optimal line upgrade selection to mitigate wildfire risk,” _Electric Power Systems Research_ , vol. 213, p. 108592, 2022.
  * [10] J. W. Muhs, M. Parvania, and M. Shahidehpour, “Wildfire risk mitigation: A paradigm shift in power systems planning and operation,” _IEEE Open Access Journal of Power and Energy_ , vol. 7, pp. 366–375, 2020.
  * [11] M. Davoudi, B. Efaw, M. Avendaño-Mora, J. L. Lauletta, and G. B. Huffman, “Reclosing of distribution systems for wildfire prevention,” _IEEE Transactions on Power Delivery_ , vol. 36, no. 4, pp. 2298–2307, 2021.
  * [12] T. Zhou, B. Li, C. Wu, Y. Tan, L. Mao, and W. Wu, “Studies on big data mining techniques in wildfire prevention for power system,” in _2019 IEEE 3rd Conference on Energy Internet and Energy System Integration (EI2)_ , 2019, pp. 866–871.
  * [13] R. Bayani, M. Waseem, S. D. Manshadi, and H. Davani, “Quantifying the risk of wildfire ignition by power lines under extreme weather conditions,” _IEEE Systems Journal_ , vol. 17, no. 1, pp. 1024–1034, 2022.
  * [14] W. Hong, B. Wang, M. Yao, D. Callaway, L. Dale, and C. Huang, “Data-driven power system optimal decision making strategy under wildfire events,” Lawrence Livermore National Lab, Tech. Rep. LLNL-CONF-831390, 2022.
  * [15] C. Huang, Q. Hu, L. Sang, D. D. Lucas, R. Wong, B. Wang, W. Hong, M. Yao, and V. Donde, “A review of public safety power shutoffs (psps) for wildfire mitigation: Policies, practices, models and data sources,” _IEEE Transactions on Energy Markets, Policy and Regulation_ , 2023.
  * [16] A. Lesage-Landry, F. Pellerin, J. A. Taylor, and D. S. Callaway, “Optimally scheduling public safety power shutoffs,” _arXiv preprint arXiv:2203.02861_ , 2022.
  * [17] R. Bayani and S. D. Manshadi, “Resilient expansion planning of electricity grid under prolonged wildfire risk,” _IEEE Transactions on Smart Grid_ , 2023\.
  * [18] Y. Zhou, K. Sundar, D. Deka, and H. Zhu, “Optimal power system topology control under uncertain wildfire risk,” _arXiv preprint arXiv:2303.07558_ , 2023.
  * [19] N. Rhodes and L. Roald, “Co-optimization of power line de-energization and restoration under high wildfire ignition risk,” _arXiv preprint arXiv:2204.02507_ , 2022.
  * [20] A. Kody, A. West, and D. K. Molzahn, “Sharing the load: Considering fairness in de-energization scheduling to mitigate wildfire ignition risk using rolling optimization,” in _2022 IEEE 61st Conference on Decision and Control (CDC)_. IEEE, 2022, pp. 5705–5712.
  * [21] G. Y. Ke, “Managing reliable emergency logistics for hazardous materials: A two-stage robust optimization approach,” _Computers & Operations Research_, vol. 138, p. 105557, 2022.
  * [22] S. Jazebi, F. de León, and A. Nelson, “Review of wildfire management techniques—Part II: Urgent call for investment in research and development of preventative solutions,” _IEEE Transactions on Power Delivery_ , vol. 35, no. 1, pp. 440–450, 2020.
  * [23] S. Li and T. Banerjee, “Spatial and temporal pattern of wildfires in California from 2000 to 2019,” _Scientific Reports_ , vol. 11, no. 1:8779, 2021.
  * [24] M. Nagpal, R. P. Barone, T. G. Martinich, Z. Jiao, S.-H. Manuel, and S. Merriman, “Wildfire trips de-energized line shunt reactor,” _IEEE Transactions on Power Delivery_ , vol. 34, no. 2, pp. 760–768, 2018.
  * [25] H. N. V. Pico and B. B. Johnson, “Transient stability assessment of multi-machine multi-converter power systems,” _IEEE Transactions on Power Systems_ , vol. 34, no. 5, pp. 3504–3514, 2019.
  * [26] J. W. Mitchell, “Power line failures and catastrophic wildfires under extreme weather conditions,” _Engineering Failure Analysis_ , vol. 35, pp. 726–735, 2013.
  * [27] J. H. Scott, M. P. Thompson, and D. E. Calkin, “A wildfire risk assessment framework for land and resource management,” U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station, Tech. Rep. RMRS-GTR-315, 2013\.
  * [28] S. Dian, P. Cheng, Q. Ye, J. Wu, R. Luo, C. Wang, D. Hui, N. Zhou, D. Zou, Q. Yu, and X. Gong, “Integrating wildfires propagation prediction into early warning of electrical transmission line outages,” _IEEE Access_ , vol. 7, pp. 27 586–27 603, 2019.
  * [29] J. R. Birge and F. Louveaux, _Introduction to Stochastic Programming_. New York, NY, USA: Springer-Verlag, 1997.
  * [30] H. Yang and H. Nagarajan, “Optimal power flow in distribution networks under N-1 disruptions: A multistage stochastic programming approach,” _INFORMS Journal on Computing_ , vol. 34, no. 2, pp. 690–709, 2022.
  * [31] J. Zou, S. Ahmed, and X. A. Sun, “Stochastic dual dynamic integer programming,” _Mathematical Programming_ , vol. 175, no. 1-2, pp. 461–502, 2019.
  * [32] J. Freire and C. Dacamara, “Using cellular automata to simulate wildfire propagation and to assist in fire prevention and fighting,” _Natural Hazards and Earth System Sciences Discussions_ , pp. 1–17, 2018.
  * [33] C. Barrows, A. Bloom, A. Ehlen, J. Ikäheimo, J. Jorgenson, D. Krishnamurthy, J. Lau, B. McBennett, M. O’Connell, E. Preston, A. Staid, G. Stephen, and J.-P. Watson, “The IEEE reliability test system: A proposed 2019 update,” _IEEE Transactions on Power Systems_ , vol. 35, no. 1, pp. 119–127, 2020\.
  * [34] MapTools, “Why use UTM coordinates,” 2023. [Online]. Available: https://www.maptools.com/tutorials/utm/why_use_utm
  * [35] “Wildland fire potential index (WFPI).” [Online]. Available: https://www.usgs.gov/fire-danger-forecast/wildland-fire-potential-index-wfpi
  * [36] A. Alexandridis, D. Vakalis, C. Siettos, and G. Bafas, “A cellular automata model for forest fire spread prediction: The case of the wildfire that swept through Spetses island in 1990,” _Applied Mathematics and Computation_ , vol. 204, no. 1, pp. 191–201, 2008.
  * [37] S. Yang, W. Zhou, S. Zhu, L. Wang, L. Ye, X. Xia, and H. Li, “Failure probability estimation of overhead transmission lines considering the spatial and temporal variation in severe weather,” _Journal of Modern Power Systems and Clean Energy_ , vol. 7, no. 1, pp. 131–138, 2019.
  * [38] C. Lemaréchal, A. Nemirovskii, and Y. Nesterov, “New variants of bundle methods,” _Mathematical programming_ , vol. 69, no. 1, pp. 111–147, 1995\.
  * [39] “Texas A&M forest service: Texas wildfire risk assessment portal.” [Online]. Available: https://texaswildfirerisk.com/
  * [40] World Nuclear Association, “Economics of nuclear power,” 2022. [Online]. Available: https://world-nuclear.org/information-library/economic-aspects/economics-of-nuclear-power.aspx
  * [41] D. Blewett, “Wind turbine cost: How Much? Are They Worth It In 2022?” 2022. [Online]. Available: https://weatherguardwind.com/how-much-does-wind-turbine-cost-worth-it/
  * [42] World Nuclear Association, “Financing nuclear energy,” 2020. [Online]. Available: https://world-nuclear.org/information-library/economic-aspects/financing-nuclear-energy.aspx
  * [43] “2021 transmission cost report.” Australian Energy Market Operator, Tech. Rep., 2021. [Online]. Available: https://aemo.com.au/-/media/files/major-publications/isp/2021/transmission-cost-report.pdf?la=en
  * [44] I. Dunning, J. Huchette, and M. Lubin, “JuMP: A modeling language for mathematical optimization,” _SIAM Review_ , vol. 59, no. 2, pp. 295–320, 2017.
  * [45] S. K. Jeff Bezanson, Alan Edelman and V. Shah, “Julia: A fresh approach to numerical computing,” _SIAM Review_ , 2014.
  * [46] Gurobi Optimization, Inc., _Gurobi optimizer reference manual._ , 2014. [Online]. Available: https://www.gurobi.com
  * [47] C. Coffrin, R. Bent, K. Sundar, Y. Ng, and M. Lubin, “Powermodels.jl: An open-source framework for exploring power flow formulations,” in _2018 Power Systems Computation Conference (PSCC)_ , June 2018, pp. 1–8.
  * [48] G. Datseris, A. R. Vahdati, and T. C. DuBois, “Agents.jl: a performant and feature-full agent-based modeling software of minimal code complexity,” _Simulation_ , p. 003754972110688, 2022.
  * [49] H. Yang, D. P. Morton, C. Bandi, and K. Dvijotham, “Robust optimization for electricity generation,” _INFORMS Journal on Computing_ , vol. 33, no. 1, pp. 336–351, 2021.