1. Introduction

• Developed architecture and algorithms for predictive analysis to identify power grid nodes at heightened risk even before wildfire events unfold. Developed algorithm utilizes environmental parameters, historical wildfire occurrences, vegetation types, and voltage data for predictive analysis.
• Developed solution approach and regional-specific risk analysis to wildfires using the Principal Component Analysis (PCA), isolating the most influential determinants of node vulnerability. Also, develop algorithm utilizes Moderate Resolution Imaging Spectroradiometer (MODIS) derived vegetation metrics, Hierarchical Density-Based Spatial Clustering of Applications with Noise (HDBSCAN), Voronoi-HDBSCAN and enhanced proximity analysis using overlay of electric grid and wildfire coordinates. Furthermore, by undertaking a comparative analysis across five distinct regions, our research elucidates region-specific risk profiles, paving the way for tailored future mitigation strategies.

2. Historical Inference from Wildfire Data

Figure 1. Classification of various US Regions using Voronoi diagrams and HDBSCAN.

Figure 1. Classification of various US Regions using Voronoi diagrams and HDBSCAN.

3. Wildfire Risk Factors

Figure 3. Real-time wildfires being displayed in a geographical mapping format.

Figure 3. Real-time wildfires being displayed in a geographical mapping format.

3.1. Enhanced Proximity Analysis between Wildfire Incidents and Electrical Grid

While there are multiple methods to determine the distance between two geographical points, the Haversine formula stands out for its accuracy, especially for significant spans on a spherical body like Earth. This formula, rooted in trigonometric principles, precisely calculates the great-circle distance, which is the shortest distance between any two points on the surface of a sphere.

For two geocoordinates $P_1({\lambda }_1,{\varphi }_1)$ and $P_2({\lambda }_2,{\varphi }_2)$:

$$\begin{array}{cc}\hfill \Delta \varphi & ={\varphi }_2-{\varphi }_1\hfill \\ \hfill \Delta \lambda & ={\lambda }_2-{\lambda }_1\hfill \\ \hfill a& ={sin}^2\left(\frac{\Delta \varphi }{2}\right)+cos\left({\varphi }_1\right)\times cos\left({\varphi }_2\right)\times {sin}^2\left(\frac{\Delta \lambda }{2}\right)\hfill \\ \hfill c& =2\times \mathrm{asin}\left(\sqrt{a}\right)\hfill \\ \hfill d& =R\times c\hfill \end{array}$$

where R is Earth’s radius, and d is the computed distance.

Figure 4. Combined overlay of both electric power grid, renewable generation and real-time wildfires.

Figure 4. Combined overlay of both electric power grid, renewable generation and real-time wildfires.

In the project’s schema, the electrical grid is represented by a graph $G(V,E)$, with V and E being sets of nodes (buses) and edges (transmission lines), respectively. The vulnerability or risk factor of a node, $r_v$, is influenced by various parameters, with the proximity to active wildfire regions, calculated using the Haversine formula, being a principal factor. This is denoted as $d(v,W)$ for a node v and wildfire location W. Specifically:

$$r_v=f(d(v,W),\mathcal{P})$$

Here, f is a risk determination function, and $\mathcal{P}$ represents other parameters influencing risk, such as meteorological conditions, grid load, and infrastructure age.

The application protocol can be articulated as follows:

• Extract real-time geocoordinates of active wildfire incidents.
• For each node $v\in V$, utilize the Haversine formula to compute $d(v,W)$.
• Integrate $d(v,W)$ into the risk function f to update the risk factor $r_v$ for each node.
• Prioritize nodes based on increasing risk values, thereby aiding in real-time grid management decisions.

A critical aspect of risk evaluation in the context of proximity analysis is the normalization of derived distances. While the Haversine formula provides an accurate measure of the great-circle distance between two geographical points, these absolute values may span a vast range. For the purpose of a uniform and comparative risk assessment, it is essential to map these distances to a normalized scale, typically between 0 and 1.

One of the most common methods for normalization is the Min-Max scaling. Given a distance $d(v,W)$ and assuming $d_{\max }$ and $d_{\min }$ are the maximum and minimum distances observed across all nodes, the normalized distance, $d_{\mathrm{norm}}(v,W)$, can be computed as:

$$d_{\mathrm{norm}}(v,W)=\frac{d(v,W)-d_{\min }}{d_{\max }-d_{\min }}$$

(2)

However, considering the context of risk, proximity to a wildfire poses a higher threat. Thus, it might be more intuitive to invert this normalized value, where a value closer to 1 indicates closer proximity and, therefore, higher risk:

$$d_{\mathrm{norm}}^{\prime }(v,W)=1-d_{\mathrm{norm}}(v,W)$$

(3)

Incorporating this normalized distance measure into the risk determination function ensures that the proximity influence is consistent across all nodes, regardless of their absolute geographical separations. Moreover, this allows for a clear comparison among nodes, where nodes with values closer to 1 are in more immediate danger and might require urgent interventions.

The Haversine formula provides a rigorous mathematical foundation to the proximity-based risk analysis of the electrical grid relative to wildfire threats. This, when combined with other risk parameters, offers a holistic and highly accurate risk assessment model essential for the safeguarding and efficient management of contemporary electrical grids.

3.2. Historical Wildfire Frequency as a Risk Factor

The spatial distribution of wildfires across different regions bears testament to not only environmental conditions but also human activities, forest management practices, and various socio-economic factors. However, for the sake of electrical grid robustness, it’s paramount to convert these spatial patterns into quantifiable metrics that can be used as risk indicators. Historical wildfire frequency emerges as a pivotal metric in this scenario.

From the HDBSCAN clustering methodology, clusters of historical wildfires are identified across the regions. Each cluster’s density provides an immediate measure of wildfire frequency for that specific region. Let’s denote the number of wildfires in a given cluster $C_i$ as $W\left(C_i\right)$.

To determine the historical wildfire frequency $F_r$ for a region r, the following equation can be employed:

$$F_r=\frac{{\sum }_{i\in Cluster{s}_r}W\left(C_i\right)}{Area\left(H_r\right)}$$

(4)

Where $Cluster{s}_r$ are the identified clusters in the region r and $Area\left(H_r\right)$ is the total area of region r.

To incorporate the influence zones identified through Voronoi polygons, a weighted frequency can be used. This takes into account not just the number of wildfire events but also their spatial influence:

$$W{F}_r=\frac{{\sum }_{i\in Cluster{s}_r}W\left(C_i\right)\times Area(V\cap H_r)}{Area\left(H_r\right)}$$

(5)

Where $Area(V\cap H_r)$ represents the union of Voronoi regions associated with wildfire events in the region r.

Figure 5. Flowchart describing the process used to integrate the historical wildfire factor into the risk model.

Figure 5. Flowchart describing the process used to integrate the historical wildfire factor into the risk model.

This weighted wildfire frequency, $W{F}_r$, provides a more nuanced understanding of the historical wildfire frequency by integrating spatial influence areas.

Given the vast differences in regional sizes and historical data availability, it’s essential to normalize these frequency values. Again Min-Max scaling is used where the minimum is 0 and the maximum is 1:

$$N{F}_r=\frac{W{F}_r-min\left(WF\right)}{max\left(WF\right)-min\left(WF\right)}$$

(6)

Where $N{F}_r$ is the normalized frequency for region r, and $min\left(WF\right)$ and $max\left(WF\right)$ are the minimum and maximum weighted frequencies among all regions respectively. This normalized metric ensures comparability across regions irrespective of their size and can be directly integrated as a factor in the risk model.

The historical wildfire frequency, especially when weighted by spatial influence, provides an invaluable insight into the regions more likely to experience future events. It is grounded on the notion that past events indicate areas of inherent vulnerability, whether due to local environmental conditions or human factors. By incorporating this frequency into the risk model, the electrical grid’s representation and analysis become more in line with real-world challenges, making it a vital component in assessing potential future vulnerabilities.

3.3. Voltage Analysis in Electrical Grid Nodes and Transmission Lines

Voltage, in power systems, is not merely an electrical parameter but a crucial indicator of the system’s health, stability, and operational efficiency. The magnitude and phase of voltage across nodes (buses) and transmission lines can shed light on a plethora of system attributes, from load dynamics to reactive power compensation.

Figure 6. Buses and transmission lines represented as nodes and edges of a graph.

Figure 6. Buses and transmission lines represented as nodes and edges of a graph.

The granular voltage information for each node and transmission line within our framework is acquired from a comprehensive dataset, meticulously curated from multiple sensors and telemetry equipment distributed across the grid. These sensors, often connected to sophisticated Supervisory Control and Data Acquisition (SCADA) systems, provide near-real-time measurements.

Let’s first denote the voltage at any node i as $V_i$, which can be expressed in polar form as:

$$V_i=\left|V_i\right|\angle {\theta }_i$$

In power systems, the power flow equation that relates voltage and power is:

$$P_i+j{Q}_i=\sum _{k=1}^n|V_i\left|\right|V_k|Y_{ik}\angle ({\theta }_{ik}-{\theta }_i+{\theta }_k)$$

For transmission lines, the voltage drop can be represented using the line impedance $Z=R+jX$:

$$\Delta V=I_{line}\times Z$$

Voltage data for each node and transmission line is dynamically linked within our advanced grid management framework.

For any node i, the voltage profile over time t can be formulated as:

$$V_i\left(t\right)=V_{base}+{\int }_0^t\Delta V_i\left(\tau \right)d\tau$$

The operational state of a node or transmission line is influenced by its voltage magnitude. The evaluation of the eigenvalues of the system’s Jacobian matrix is the following:

$$J=\left[\begin{array}{cc}\partial P/\partial \theta & \partial P/\partial \left|V\right|\\ \partial Q/\partial \theta & \partial Q/\partial \left|V\right|\end{array}\right]$$

Voltage deviations, both sag (when demand is higher than generation) and swell (when generation is higher than demand) [17], and can indicate potential grid vulnerabilities. The normalization of voltage magnitudes ensures that the operational state of each node can be compared, offering a unified measure of vulnerability.

Given a node’s voltage magnitude $|V_i|$ and taking into account the permissible voltage range $[V_{\min },V_{\max }]$ specified by grid standards, the normalized voltage magnitude $V_{\mathrm{norm},\mathrm{i}}$ can be computed using the Min-Max normalization.

However, from a risk perspective, significant deviations from the nominal voltage value (either too high or too low) are more concerning. Thus, it might be valuable to use a modified normalization scheme that accentuates deviations:

$$V_{\mathrm{norm},\mathrm{i}}^{\prime }=1-\left|V_{\mathrm{norm},\mathrm{i}}-0.5\right|\times 2$$

(7)

Here, $V_{\mathrm{norm},\mathrm{i}}^{\prime }$ inverts the normalized voltage such that values approaching 0 indicate nominal operation, and as the voltage deviates from nominal (either due to sag or swell), $V_{\mathrm{norm},\mathrm{i}}^{\prime }$ increases, approaching 1. This renders nodes with a $V_{\mathrm{norm},\mathrm{i}}^{\prime }$ value closer to 1 as more vulnerable, warranting monitoring or corrective action.

This normalization approach ensures that voltage magnitudes, irrespective of their absolute value, contribute consistently to the risk assessment metric. By centering the scale on nominal operational values, and expanding outward to encompass extreme vulnerabilities, grid operators and analysts can prioritize nodes based on their deviation from desired operational standards.

It’s worth noting that while this normalization provides a framework for assessing vulnerability, a comprehensive risk assessment might necessitate further refinements considering other voltage-related parameters like voltage stability margins, phase imbalances, and harmonic distortions. By integrating real-time voltage data and coupling it with mathematical models, our framework offers insights into the operational health of the grid.

4. Modeling Risk of Wildfires to Power Grids

4.1. Data Representation and Principal Component Analysis (PCA)

Given a dataset X with n power grid nodes (rows) and 4 risk factors (columns): distance from the nearest real-time wildfire, vegetation, voltage, and historical wildfire frequency, we aim to understand the significance of each factor in the context of wildfire risk. The matrix representation of the dataset is:

$$X=\left[\begin{array}{cccc}{d}_1& v_1& v{o}_1& h_1\\ d_2& v_2& v{o}_2& h_2\\ ⋮& ⋮& ⋮& ⋮\\ d_n& v_n& v{o}_n& h_n\end{array}\right]$$

where:

• d represents the distance from the nearest real-time wildfire.
• v represents vegetation.
• $vo$ represents voltage.
• h represents historical wildfire frequency.

Each column (risk factor) of X is mean-centered:

$$X_{ij}=X_{ij}-\frac{1}{n}\sum _{i=1}^n{X}_{ij}$$

The covariance matrix is a crucial component in PCA as it captures the pairwise covariances between the different features in the dataset. The covariance between two features indicates how much the features vary together. A positive covariance indicates that as one feature increases, the other also tends to increase, while a negative covariance indicates that as one feature increases, the other tends to decrease.

Given our centered data matrix X of size $n\times m$ (where n is the number of data points and m is the number of features), the covariance matrix C of size $m\times m$ is computed as:

$$C=\frac{1}{n-1}{X}^{T}X$$

Each element $C_{ij}$ of the covariance matrix represents the covariance between the $i^{th}$ feature and the $j^{th}$ feature and is given by:

$$C_{ij}=\frac{1}{n-1}\sum _{k=1}^n(x_{ki}-{\overline{x}}_i)(x_{kj}-{\overline{x}}_j)$$

where $x_{ki}$ is the value of the $i^{th}$ feature for the $k^{th}$ data point, and ${\overline{x}}_i$ is the mean of the $i^{th}$ feature.

The diagonal elements of the covariance matrix, $C_{ii}$, represent the variance of the $i^{th}$ feature. The variance measures the spread or dispersion of a feature around its mean. The covariance matrix provides insights into the relationships between features. Eigen decomposition of the covariance matrix is used in PCA to determine the principal components, which are the directions of maximum variance in the data.

Eigen decomposition is a fundamental operation in linear algebra that involves decomposing a matrix into its constituent eigenvalues and eigenvectors. In the context of PCA, Eigen decomposition of the covariance matrix C reveals the principal components of the data. Given the covariance matrix C, the eigenvalues $\lambda$ and the corresponding eigenvectors v satisfy the equation:

$$Cv=\lambda v$$

The Eigenvector v represents a direction in the feature space, while the corresponding Eigenvalue $\lambda$ indicates the variance of the data along that direction. In other words, the magnitude of the Eigenvalue signifies the importance or the amount of variance captured by its corresponding Eigenvector.

The steps involved in the Eigen decomposition process are as follows:

• The first step is to compute the eigenvalues. The eigenvalues of C are the solutions to the characteristic equation:

  $$\det (C-\lambda I)=0$$

  where I is the identity matrix of the same size as C.
• Next, we have to compute the eigenvectors. For each eigenvalue $\lambda$, the corresponding eigenvector v is found by solving the linear system:

  $$(C-\lambda I)v=0$$
• Once all eigenvalues and eigenvectors are computed, they are ordered in decreasing order of the eigenvalues. The eigenvector corresponding to the largest eigenvalue represents the direction of maximum variance in the data, known as the first principal component. Subsequent eigenvectors represent orthogonal directions of decreasing variance.
• In PCA, it’s common to select the top k eigenvectors (principal components) that capture the most variance in the data. This allows for a reduction in dimensionality while retaining most of the data’s original variance.

The Eigen decomposition of the covariance matrix C provides a basis transformation where the new axes (principal components) are the directions of maximum variance in the data. This transformation is crucial for dimensionality reduction and feature extraction in PCA.

Given the eigenvalues ${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_m$ and their corresponding eigenvectors $v_1,v_2,\cdots ,v_m$, we first sort the eigenvalues in descending order:

$${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_m$$

The sorted eigenvalues have corresponding eigenvectors, which we denote as:

$$v_{\left(1\right)},v_{\left(2\right)},\cdots ,v_{\left(m\right)}$$

To reduce the dimensionality from m dimensions to k dimensions (where $k<m$), we select the first k eigenvectors:

$$F=[v_{\left(1\right)},v_{\left(2\right)},\cdots ,v_{\left(k\right)}]$$

This matrix F is our feature vector, and it will be used to transform the original data matrix X into a reduced dimensionality matrix Y

Given our original data matrix X of size $n\times m$ (where n is the number of data points and m is the number of features) and our feature vector F of size $m\times k$ (where k is the number of selected principal components), we can project the data onto the lower-dimensional space by multiplying X with F.

The transformed data Y is obtained by:

$$Y=XF$$

Each row of Y represents a data point in the original dataset but now transformed into the new lower-dimensional space spanned by the principal components. The columns of Y represent the coordinates of the data points in this new space.

Mathematically, the $i^{th}$ row of Y, denoted as $y_i$, is given by:

$$y_i=x_{i}F$$

where $x_i$ is the $i^{th}$ row of X.

This projection essentially captures the most significant patterns in the data while discarding the less important variations. The principal components in F act as the new axes, and the data is represented concerning these axes, ensuring that the variance (or information) is maximized in this reduced space.

Eigenvalues in PCA represent the variance captured by each principal component. Their magnitude provides insights into the significance of each component.

• Eigenvalue Interpretation: Each eigenvalue ${\lambda }_i$ indicates the variance explained by its corresponding eigenvector. A larger ${\lambda }_i$ denotes greater significance.
• Total Variance: Given by:

  $$\mathrm{Total}\mathrm{Variance}=\sum _{i=1}^m{\lambda }_i$$

  where m is the number of eigenvalues.
• Proportion of Variance: For the $i^{th}$ component:

  $${\mathrm{Proportion}}_i=\frac{{\lambda }_i}{\mathrm{Total}\mathrm{Variance}}$$
• Weight Derivation: The proportion of variance explained by a component represents the weight of the corresponding risk factor. For instance, if a component explains 50% variance, its weight is 0.5.
• Ranking Risk Factors: Risk factors can be ranked by ordering the eigenvalues in descending order.

In the context of wildfire risk, these weights help prioritize interventions based on the significance of each risk factor.

For each risk factor j:

$${\mathrm{Variance}}_j=\frac{1}{n}\sum _{i=1}^n{(X_{ij}-{\overline{X}}_j)}^2$$

$$\mathrm{Std}.{\mathrm{Deviation}}_j=\sqrt{{\mathrm{Variance}}_j}$$

For risk factors j and k:

$${\mathrm{Covariance}}_{jk}=\frac{1}{n}\sum _{i=1}^n(X_{ij}-{\overline{X}}_j)(X_{ik}-{\overline{X}}_k)$$

Each eigenvector v is normalized:

$$v=\frac{v}{\parallel v\parallel }$$

The proportion of variance p captured by the $k^{th}$ principal component is:

$$p_k=\frac{{\lambda }_k}{{\sum }_{i=1}^m{\lambda }_i}$$

The cumulative variance captured by the first k principal components is:

$$P_k=\sum _{i=1}^k{p}_i$$

Table 1. Risk Factors and their Computed Weights.

Table 1. Risk Factors and their Computed Weights.

Risk Factor	Weight
Distance from nearest real-time wildfire	0.4497
Historical wildfire frequency	0.2986
Vegetation Information	0.1984
Voltage	0.0533

5. Results: Computation of Risk Factor Analysis & Discussion

5.6. Node Ranking

The nodes ranked from lowest to highest risk are:

• Node 1 (Western Region - Californian coast): 0.7919
• Node 5 (Southwestern Region - Arizona’s desert): 0.6409
• Node 3 (Northeastern Region - Forests of Vermont): 0.4907
• Node 4 (Southeastern Region - Florida’s wetlands): 0.4462
• Node 2 (Midwestern Region - Plains of Kansas): 0.3259

Table 2. Risk values associated with nodes from various U.S. regions.

Table 2. Risk values associated with nodes from various U.S. regions.

Node	Distance from Wildfire	Historical Wildfire Frequency	Vegetation Info	Voltage
Node 1(California)	0.823764	0.783512	0.764915	0.671243
Node 2(Kansas)	0.250358	0.330914	0.418276	0.592117
Node 3(Vermont)	0.394821	0.452317	0.730256	0.623489
Node 4(Florida)	0.352715	0.418944	0.654727	0.610596
Node 5(Arizona)	0.724938	0.691423	0.374526	0.641278

6. Conclusions

References

1. Sathaye, J.; Dale, L.; Larsen, P.; Fitts, G.; Koy, K.; Lewis, S.; others. Estimating risk to California energy infrastructure from projected climate change. Technical report, 2011.
2. Sathaye, J.A.; Dale, L.L.; Larsen, P.H.; Fitts, G.A.; Koy, K.; Lewis, S.M.; de Lucena, A.F.P. Rising Temps, Tides, and Wildfires: Assessing the Risk to California’s Energy Infrastructure from Projected Climate Change. IEEE Power and Energy Magazine 2013, 11, 32–45. [Google Scholar] [CrossRef]
3. Frame, D.; Rosier, S.; Noy, I.; Harrington, L.; Carey-Smith, T.; Sparrow, S.; others. Climate change attribution and the economic costs of extreme weather events: a study on damages from extreme rainfall and drought. Climatic Change 2020, 162, 781–97. [Google Scholar] [CrossRef]
4. Wang, X.; Bocchini, P. Predicting wildfire ignition induced by dynamic conductor swaying under strong winds. Sci Rep 2023, 13, 3998. [Google Scholar] [CrossRef] [PubMed]
5. Pacific Gas and Electric Company. 2020 Wildfire Mitigation Plan Report 2020.
6. San Diego Gas & Electric Company. Wildfire Mitigation Plan 2020.
7. Southern California Edison Company. Wildfire Mitigation Plan Update 2022.
8. California Public Utilities Commission. Deenergization (PSPS) 2020.
9. Baker, D. Underground power lines don’t cause wildfires. But they’re really expensive. San Francisco Chronicle, 2017.
10. Nazaripouya, H. Power Grid Resilience under Wildfire: A Review on Challenges and Solutions. 2020 IEEE Power & Energy Society General Meeting (PESGM); , 2020; pp. 1–5. [CrossRef]
11. Dale, L.; Carnall, M.; Wei, M.; Fitts, G.; McDonald, S. Assessing the Impact of Wildfires on the California Electricity Grid. Technical Report CCCA4-CEC-2018-002, California Energy Commission, 2018.
12. Rhodes, N.; Ntaimo, L.; Roald, L. Balancing Wildfire Risk and Power Outages Through Optimized Power Shut-Offs. IEEE Transactions on Power Systems 2021, 36, 3118–3128. [Google Scholar] [CrossRef]
13. Zanin Bertoletti, A.; Campos do Prado, J. Transmission System Expansion Planning under Wildfire Risk. 2022 North American Power Symposium (NAPS), 2022, pp. 1–6. [CrossRef]
14. Astudillo, A.; Cui, B.; Zamzam, A.S. Managing Power Systems-Induced Wildfire Risks Using Optimal Scheduled Shutoffs. 2022 IEEE Power & Energy Society General Meeting (PESGM), 2022, pp. 1–5. [CrossRef]
15. Nazemi, M.; Dehghanian, P. Powering Through Wildfires: An Integrated Solution for Enhanced Safety and Resilience in Power Grids. IEEE Transactions on Industry Applications 2022, 58, 4192–4202. [Google Scholar] [CrossRef]
16. Paryati.; Salahddine, K.; Salah-ddine, K. THE IMPLEMENTATION OF KRUSKAL’S ALGORITHM FOR MINIMUM SPANNING TREE IN A GRAPH. 2021.
17. Agüera-Pérez, A.; Palomares-Salas, J.C.; González de la Rosa, J.J.; Sierra-Fernández, J.M.; Ayora-Sedeño, D.; Moreno-Muñoz, A. Characterization of electrical sags and swells using higher-order statistical estimators. Measurement 2011, 44, 1453–1460. [Google Scholar] [CrossRef]
18. Andrews, P.L. BEHAVE: Fire Behavior Prediction and Fuel Modeling System - BURN Subsystem, Part 1. Technical Report INT-194, United States Department of Agriculture - Forest Service, Intermountain Research Station, Ogden, UT, 1986.