Vulnerability and Adaption of Power Grids for Global Warming Induced Energy Demands

Global climate change not only causes extreme weather events but also impacts the consumer’s energy demand, posing 5 a concern for the safe and reliable operation of power grids. Furthermore, the high adaption rate of electric vehicles makes 6 energy transactions volatile. Hence, it is crucial to evaluate the power grid’s reliability using risk assessment studies. Here, 7 we solve this open problem to predict the risk of blackout by developing a real-world data-driven probabilistic framework 8 in conjunction with an expert-based power grid methodology. Utilizing the real-world energy demands, climate station’s 9 temperature data, and IPCC regional climate models, the risk of a city in California increases up to 8% in the summer 10 of 2100 . The winter, spring, summer, and fall seasons’ risks also follow an increasing trend from 2000 to 2100 . We also 11 show that a better localization planning of load hot spots like EV charging stations reduces the risk by 1.8%.

(a) Blackout boundary of the triangle shape 3 bus system in Fig. 2a. (b) Blackout boundary of the power grid and consumers' power demand.
(c) Blackout boundary, consumers' power demand before and after influence of global climate change.

52
The exacerbating change in the climate causes a long-term increase in temperature. This increase in temperature results in a 53 raise in the population's power consumption median value and burdens the power system to operate beyond its original capacity. 54 Thus climate change deteriorates the safety and reliability of the power grid. Several risk assessment methods are proposed 55 based on the historical failure rates of physical components (infrastructure) in the power grid. However, a well-designed power 56 grid may have control to ensure a very low risk of infrastructure failure through redundancy, but it may still experience a 57 system blackout due to the abnormal and less controllable nature of energy demand growth of its users. This is due to climate 58 change or extreme events like heat waves. This abnormal energy demand growth can cause the consumers on the grid to 59 consume more power than the maximum power transfer limits of what the power grid can withstand. Hence, it is required 60 to first identify the "original capacity" i.e., maximum power transfer limits (watts) of the power grid which is referred to as 61 blackout/loadability boundary and distance of the grid's operating condition to the blackout boundary as margin to collapse.  Fig. 1b). The power grid collapses when the power grid's demand due to its users 72 on nodes A and B crosses the blackout boundary in Fig. 1a. Hence, the risk (blackout probability) is defined as the ratio of 73 the shaded area in Fig. 1b to the total area of the green ellipse. Formally, the risk or blackout probability is defined as the 74 probability of the power grid's demand due to its users crossing the blackout boundary of the power grid. of the earth's surface [23]. For example, the increase in temperature may cause consumers A and B to consume more power 79 and their new increased demand consumption pattern is represented by the blue ellipse in Fig. 1c. It can be observed that, 80 when compared to the green ellipse (without climate change impact), the blue ellipse (with climate change impact) moves 81 further away from the blackout boundary. Therefore, due to climate change, the area outside the blackout boundary for the blue 82 ellipse increases more than the green ellipse which exacerbates the risk of power grid collapse. Additionally, not only the user's 83 demand (ellipses) is impacted by the GCC but also the blackout boundary is impacted. For example, when the temperature is 84 increased, the transmission line's conductance and susceptance change causing the power grid to weaken. As shown in Fig. 1c Here, using mesh and tree type networks, we show how demand growth due to climate change can impact the vulnerability of 92 power grids and how the localization of loads can help adapt the power grids to climate change by reducing the vulnerability. 93 In later sections, we extend this demonstration to real-world cases using real-world power demand, climate change data, 94 spatiotemporal temperature data, EV data, and larger power grids.  (a), A three bus mesh type network with two consumers located at nodes 2 and 3 whose demands are given by P 2 and P 3 respectively. (b), Margin to collapse value for every possible pair in (P 2 , P 3 ) plane represented using the color bar. Note that the margin to collapse values of the data points on the border are zeros. The margin to collapse values of the data points (in (P 2 , P 3 ) plane) decreases when moving towards the boundary. (c), The operation of power grid is labeled as desired and undesired operating regions by using a threshold (α = 0.5 p.u.) on the margin to collapse values. For example, if the margin to collapse value of a operating condition (i.e., (P 2 ,P 3 ) pair value) is greater than 0.5 p.u. (blue region) then it is desired and otherwise undesired (orange region). Note that when the threshold (α) is zero, the undesired region indicates the power grid collapse events.

95
System information: Specifically, the 3-bus mesh and tree networks presented in Fig. 2a and Fig. 3a are used as examples 96 whose transmission lines are purely resistive i.e., Y 12 = Y 23 = Y 13 = 1p.u.. For both mesh and tree networks, node 1 has 97 the generator and nodes 2, 3 have the users whose power demands are indicated by P 2 , P 3 respectively.

98
Blackout/loadability boundaries of the mesh and tree networks: The maximum power transfer limit (blackout boundary) 99 of a power grid is dependent on its design and network structure. Using the Pareto front technique (Φ from Methods), the 100 blackout/loadability boundary can be charted for mesh and tree networks as shown in Fig. 2b and Fig. 3b respectively. We can 101 observe in both Fig. 2b and Fig. 3b, a heat map is drawn for every possible pair of (P 2 , P 3 ) and the color of each data point 102 from color bar indicates the distance (margin) to system collapse. It can be noted that the pairs of (P 2 , P 3 ) that contribute to 103 system collapse (located on the outer boundary of these heat maps) in Fig. 2b and Fig. 3b have their margin to collapse (color 104 bar) values (Φ values from Methods) to be zero. The operating points (P 2 , P 3 ) beyond the boundaries in Fig. 2b and Fig. 3b 105 result in a blackout as the power networks cannot host any additional real power demand beyond the loadability boundary. 106 Safe and unsafe operating regions of power grids: In practice, power system operators are more cautious about the worst-107 case scenario (blackout). Hence, to be on a safer margin, the operators prefer to know the probability of power grid operating 108 like 90% close to the blackout boundary. In this article, such an analysis is made possible for the first time, by assigning a 109 positive threshold value (in this example 90%) to α (from Methods). For example, in the case of the 3-bus mesh network in 110 Fig. 2c, the operator may want to identify the operation of the power grid as unsafe (orange region) when the consumers' 111 demand (data points in (P 2 , P 3 ) plane) drives the margin to collapse (color bar) values to range between 0.5 p.u. to 0 p.u. 112 (note that a 0 p.u. indicates the blackout of the power grid). In the same manner, from Fig. 2c, the operation of a power grid 113 is labeled as safe (blue region) when the consumer demand drives the margin to collapse values to be greater or equal to 0.5 114  , A three bus tree type network with two consumers located at nodes 2 and 3 whose demands are given by P 2 and P 3 respectively. (b), Margin to collapse value for every possible pair in (P 2 , P 3 ) plane represented using the color bar. Note that the margin to collapse values of the data points on the border are zeros. The margin to collapse values of the data points (in (P 2 , P 3 ) plane) decreases when moving towards the boundary. (c), The operation of power grid is labeled as desired and undesired operating regions by using a threshold (α = 0.3 p.u.) on the margin to collapse values. For example, if the margin to collapse value of a operating condition (i.e., (P 2 ,P 3 ) pair value) is greater than 0.3 p.u. (blue region) then it is desired and otherwise undesired (orange region). Note that when the threshold (α) is zero, the undesired region indicates the power grid collapse events.
p.u. Similarly, in the case of 3-bus tree network, as shown in Fig. 3c, the safety threshold value of α = 0.3 p.u. is selected.

115
Thus, depending on the power demand (load) of the consumers, the power grid may operate in a safe or unsafe condition.

116
Load scenario types: For the 3-bus mesh and tree networks presented in Fig. 2a and Fig. 3a, the consumers are located at 117 nodes 2 and 3. Generally, the consumers' power demand is dependent on the temperature they are exposed to. For demonstration, 118 in this example we directly consider three load scenario types for the 3-bus mesh and tree networks. The three load types 119 are as follows: 1) load type 1 (from Fig. 4a, green ellipse): power consumption of two users at the temperature T 1 and 120 among which node 3 user's power consumption is highly sensitive to temperature, while node 2 user's power consumption is 121 less sensitive to temperature, 2) load type 2 (from Fig. 4b, blue ellipse): power consumption of the same users discussed in 122 load type 1 but at a higher temperature T 2 (T 2 > T 1 ). We can observe from Fig. 4b that as the temperature increases to T 2 , 123 the consumer located at node 3 tends to consume more power while the consumer at node 2 tends to consume less power, 124 and 3) load type 3 (from Fig. 4c, magenta ellipse): In this setup, we consider same power consumption behaviors of the 125 two users from load type 2 (at temperature T 2 ) but the locations (nodes) of the consumers are swapped. For example, the 126 temperature sensitive consumer is placed at node 2 while the less temperature sensitive consumer is placed at node 3. Hence, 127 for load type 3, we observe that the consumer at node 2 consumes more power while consumer at node 3 consumes less 128 power. Thus, by comparing the calculated blackout risks of load types 1, 2, and 3, we can infer that 1) weather sensitivity 129 (load type 1 versus load type 2 i.e., demand growth) and 2) the location of a consumer (load type 2 versus load type 3 i.e., 130 node sensitivity/localization) impacts the risk of the power grid blackout.

131
Power grid vulnerability assessment: Here, we calculate the risk of a power grid operating in an unsafe region (α < 0.5 ellipse, and orange region and beyond) to the total samples drawn from the joint distribution of load type 1 (green ellipse). A 138 similar approach is repeated to calculate the risks for load types 2 and 3 as shown in Fig. 4f and Fig. 4g. Finally, the risks due 139 to load types 1, 2, and 3 are also calculated for the 3-bus tree network as shown in Fig. 4i, Fig. 4j, and Fig. 4k, respectively.

140
The results are presented in Tab  (f) Impact of load type 2 on mesh network vulnerability.
(g) Impact of load type 3 on mesh network vulnerability.
(h) 3-bus tree network (i) Impact of load type 1 on tree network vulnerability.
(j) Impact of load type 2 on tree network vulnerability.
(k) Impact of load type 3 on tree network vulnerability. , Load type 2: power demand of the consumers when temperature is T 2 (T 2 > T 1 ). The consumer at node 3 is more temperature sensitive than the consumer at node 2 and hence demands more power relatively. Joint distribution of (P 2 , P 3 ) ≈ N (µ, Σ), µ = [0.25, 0.25], Σ = [0.009 0.009; 0.009 0.018]. (c), Load type 3: P 2 and P 3 are real power demand at nodes 2 and 3 respectively when temperature is T 2 (T 2 > T 1 ). Unlike load type 2, the consumer at node 2 is more temperature sensitive than the consumer at node 3 and hence demands more power relatively. Joint distribution of (P 2 , P Climate change-induced demand growth impact on power grid vulnerability: Climate change increases the tempera-146 ture of the earth which directly impacts the power demand patterns of the consumers on the power grid. This increase in 147 temperature and therefore, an increase in power demand can be understood by comparing the load type 1 and 2 scenarios. 148 From Tab I, it can be observed that the probability to enter an unsafe operating region increases for load type 2 (at temperature 149 T 2 ) when compared to that of the load type 1 (at temperature T 1 , T 2 > T 1 ) from 62% to 89% (for mesh network) and from 150 93% to 96% (for tree network). This shows that the risk of the power grid is directly proportional to the power consumption 151 sensitivity of the consumers on the grid and climate change-induced demand growth.

152
Adaption of power grids to climate change using localization of load centers: Similarly, from Tab I, we can also ob-153  TABLE I: Given a specific hour of the day, the probability of power grid failure due to three different load types for the 3-bus mesh and tree networks in Fig. 4d and Fig. 4h respectively. Here, we use the approach described in the above section to analyze a modified real-world urban grid when it serves a 160 hypothetical demand whose patterns are summarized from a real-world city. In this section, we provide the information relevant 161 to the quality and nature of the data that is used to assess the system risk in this article.

162
The information required to compute the risk of the power grid are 1) the topology of the power grid along with its  For example, Feeder 1 is a grid with 293 nodes. Each of these nodes is assigned a random user by selecting a user from the 183 total 42 consumer database.

184
Historical temperature data: Using the geolocation and timestamp of the hourly load data points of each consumer from 185 the real-world consumer database, the historical hourly temperature data experienced by each consumer at the same time from 186 2015 to 2019 are collected from the nearest weather stations. This is collected using a weather API developed by the National 187 Oceanic and Atmospheric Administration (NOAA) [24]. 188 Using the geolocation and timestamp tagged load consumption and temperature data of different consumers from real-world 189 consumer databases and real-world weather stations of NOAA respectively, the conditional distributions of each consumer's 190 demand given the temperature (temperature-load response curves) are learned for different consumers. The learning process of 191 conditional distribution is further discussed in detail, later.

192
Climate model: Even though the historical temperature data from real-world weather stations of NOAA plays a crucial role 193 in data-driven analysis, it is not sufficient to forecast the vulnerabilities of the power grid due to the emission of greenhouse 194 gases i.e., global warming. To anticipate the vulnerability of the power grid accurately, a temperature forecasting model that 195 uses the information about global warming causing factors is necessary. To address this, two forecasting models are used, 196 one for the short-term and the other for the long-term temperature forecasting. In the short-term forecasting of temperature, 197 it is dominantly affected by the time of the day; hence NOAA weather API is used to obtain the short-term temperature 198 forecast data [24]. In the long-term forecasting of temperature, the temperature is significantly affected by global warming and 199 thereby the mean of the long-term data forecast changes. We use the representative concentration pathway (RCP) scenarios 200 that consider the impact of global warming to forecast the long-term temperature data as shown in Fig. 7 [23], [25]. Both of 201 the models obtained correspond to the geolocation of the real-world power grid locations in California. This spatiotemporal 202 analysis accurately models climate change both in short-term and long-term cases as a function of time and expected global 203 warming emissions respectively.

204
Electric vehicle charging data: The power grid is not only sensitive to the temperature dependent energy demand response 205 of the consumers, but it is also sensitive to the spatial effects. Here, we discuss the EV demand data that impacts the energy 206 demand of the power grid by considering factors such as the location of prominent communities with EVs, charging stations, 207   Here, we analyze the distribution of the EV charging profiles for various locations using zip codes. We also highlight the 213 spatial impact of EV loads.  is made possible by coalescing the power consumption data of unique users from real-world consumer databases, and the 230 temperature data from the weather API of NOAA using their hourly resolution timestamps.

231
For example, Fig. 9a is an error bar plot that shows a typical daily temperature versus load consumption response of a 232 consumer at a selected hour 2 (24-hour clock) in a year. We identify that most of the temperature response curves are parabolic 233 in nature. These temperature response functions are highly dictated by the usage of cooling and heating units as they are 234 major loads in households. We observe that these temperature response curves are "U" shaped and they are roughly symmetric 235 around an equilibrium temperature value as shown in Fig. 9a. From Fig. 9a, it can be observed that as the outside temperature   users in the available dataset use both cooling and heating units depending on whether they feel hot or cold respectively. Later 240 we use this knowledge of the temperature-load response function to predict the user's power consumption behavior based on 241 temperature.

242
It can be inferred by looking at the equilibrium temperature value of a consumer to know when the consumer stops using 243 the heating unit and starts using the cooling unit. For example, Fig. 9b presents the temperature response curve of two different 244 users, "A" and "B". It can be observed that the user "A" prefers to use the cooling unit than the heating unit when the 245 temperature is greater than 10 deg C. However, the user "B" prefers to use the cooling unit rather than the heating unit when 246 the temperature is greater than 15 deg C. Hence, we make sure to model the temperature-load response curve separately for 247 every user depending on the user's equilibrium temperature.

248
However, it is not sufficient to model the power demand only conditioned on the temperature since the power demand of a 249 consumer also depends on the time of the day. So, we identified to model the demand response functions of different consumers 250 considering both the temperature and the time of the day. For example, Fig. 10 shows the annual average power consumption 251 behavior of the user for the temperature and hour of the day in a year. It can be observed that at hours 12 to 15 (12pm to 252 3pm), there are no bar plots compared to other hours of the day. This is because the noontime is hot and the bar plots at hours 253 12 to 15 are only half of the "U" shape indicating the usage of only cooling units. Therefore, it is important to model the 254 demand response (P i ) curve for every unique user (say i) given both temperature (t) and hour of the day (h) i.e., P i |t, h.

255
The temperature values are divided into bins of width 3 deg C. The parameters of the distribution P |t, h = [(P 1 , P 2 , · · · , P n )|t, h]256 for n different consumers given a selected temperature bin value t and hour of the day h are learned by assuming that 257  P 1 , P 2 , · · · , P n are all conditionally independent given the temperature t and hour of the day h. Therefore, the parameters 258 of the conditional distributions for n different consumers are learned independently i.e., P |t, h = [(P 1 , P 2 , · · · , P n )|t, h] = 259 (P 1 |t, h) · (P 2 |t, h) · · · (P n |t, h). This conditional independence is assumed to handle the incomplete datasets corresponding to 260 the estimation of the joint distribution of total demand on power grid (P ) given a temperature and hour of the day. In other 261 words, the power consumption of two households can be estimated individually based on the outside temperatures they are 262 experiencing and the hour of the day. Finally, the parameters of each distribution of P i |t, h ∀i = 1, 2, · · · , n are estimated 263 by using a quadratic curve fitting (in a least square sense) since the shape of the temperature-load response function is "U" 264 shaped as shown in Fig. 9b.

265
Impact of real-world climate change induced demand growth on power grid vulnerability: Here, we evaluate the power 266 grid vulnerability using the learned conditional distribution (P |T, h) from the real-world consumers (Fig. 6), historical spa-267 tiotemporal temperature (Fig. 10, Fig. 9), and climate model databases (Fig. 7). For example, as shown in Fig. 11 to the climate change scenarios (from Fig. 7) for each node in the power grid (from Fig. 5). Using the IPCC temperature 283 predictions, the mean load consumption of the consumers at all the nodes in the grid are predicted. Finally using these predicted 284 load consumption values, the margin to collapse value is calculated for each bus using Φ from Methods. The color bar in 285 Fig. 14 indicates the margin to system collapse at each node in the real grid from Fig. 5. The size of each data point in the 286 figure is downscaled as the time (years) increases. It can be observed that the power grid's margin to collapse values decrease 287 as the years increase, due to the global warming-induced increase in energy demands. This deterioration is even worse in case 288 of RCP 8.5 climate change scenario from Fig. 14b when compared to the RCP 4.5 climate change scenario from Fig. 14a. 289 Therefore, it is crucial to identify the critical locations of the grid that makes it vulnerable, and the proposed framework easily 290 enables this identification.

291
The location of the intermittent energy resources and temperature-sensitive consumers can play an important role in the 292 power grid security. For example, Fig. 13 shows a layered abstraction of the power grid and how localization of loads impacts 293 the vulnerability of it. Specifically, Fig. 13a shows three layers 1) the satellite view (layer) of the distribution grid from Fig. 5 294 on the google map, 2) the schematic layer shows the one-line diagram of the actual topology of a selected Feeder 1, and 295 3) the graph layer shows the detailed nodes and edges of the selected Feeder 1. Each node in the graph layer is assigned a 296 consumer. Using each consumer's load and temperature profile data, the specific heat loss in KW/K (Kilo Watt per Kelvin) is 297 calculated. This KW/K metric helps to identify the most temperature-sensitive consumers using the historical temperature and 298 power consumption data. The specific heat loss corresponding to each consumer in the graph layers of Fig. 13a and Fig. 13b 299 are represented using a heat map. The difference between Fig. 13a and Fig. 13b is that the consumers are relocated at different 300 nodes (locations) in the grid. The margin to collapse values are presented for the nodes (nodes 55, 162, and 235) in Fig. 13a 301 and Fig. 13b using both climate change RCP 4.5 and RCP 8.5 scenarios (from Fig. 7). From Fig. 13a, we can observe that 302 the margin to collapse values for the nodes (node 162) that are closer to the substation have larger values when compared 303 to the nodes (nodes 55 and 235) that are farther from the substation. This shows that the vulnerability of power grids is not 304 only dependent on the amount of power consumption, but they are also dependent on the location of the loads/consumers. The 305 nodes 55, 162, 235 in Fig. 13b have lower heat-sensitive consumers when compared to that of the consumers in Fig. 13a. We 306 can also observe that the nodes with lower heat-sensitive consumers from Fig. 13b have higher margin to collapse values when 307 compared to that of the values in Fig. 13a. This indicates that the amount of power consumption also impacts the risk along 308 with the location of the loads/consumers. The risks being relatively high for the nodes that are farther from the substation can 309 also be observed here.

311
Here, a Monte Carlo simulation is proposed to assess the probability of a power grid's operating condition to violate either 312 the blackout boundary or a user-defined margin threshold to collapse due to global climate change. The global climate change 313 is modeled using the RCP scenarios [23], [25] which provide an estimate of change in temperature due to greenhouse gas 314 emissions.

315
The risk of the power grid operating in an unsafe/undesired region (α watts distance close to the system blackout event) at a given hour (h) of the day is given by p(T i |h) · p(Φ[(P 1 , P 2 , · · · , P n ) |T i , h] ≤ α), (1a) When α = 0, equation (1) represents the probability of power grid blackout.

316
The overall step-by-step procedure is as follows, first given that the temperature model (T ) due to climate change [23], 317 hour of the day (h) are known, the conditional distribution of the total demand on power grid (P |T, h) is learned. Where 318 P = [P 1 , P 2 , · · · , P n ], P i represents the power consumption of a consumer located at node i, and n represents the total number 319 of nodes on the power grid.

320
Once the distribution of (P |T, h) is learned, we can draw samples (N (P |T i , h)) of user's conditional power demand 321 distribution (P |T i , h) given a specific temperature (T i ) from a climate scenario of T = [T 1 , T 2 , · · · , T t ], and the hour of the 322 day (h) at which the risk has to be computed.

323
For every sample drawn from N (P |T i , h), Pareto front technique (Φ) from [19], [27] is used to estimate the distance 324 (margin to collapse) between the sample and the blackout boundary (for system-level margin to collapse values [27] is used 325 and for each bus margin to collapse values [19] is used). For example, when this margin to collapse Φ[P |T i , h] of a sample 326 is equal to zero then the power grid collapses i.e., Φ[P |T i , h] = 0.

327
The proposed framework enables the assessment of power grid security by giving the choice to the power system operators 328 to define a safe and unsafe operating region. For example, the power system operator may define the operating region of a 329 power grid to be safe when Φ[P |T, h] = m > α watts and unsafe when Φ[P |T, h] = m ≤ α watts where α is a threshold 330 in available watts that the grid can host without an outage event. Therefore, using Monte Carlo simulations, we can count the 331 number of samples drawn from the conditional demand distribution (P |T i , h) that drive the power grid to operate in an unsafe 332 region (N (Φ[P |T i , h] ≤ α)) and finally compute the risk as described in equation (1).

333
(a) Scenario when heat sensitive consumers are located at critical locations in the grid.
(b) Scenario when less heat sensitive consumers are located at critical locations in the grid.

334
All the authors contributed extensively to the work presented in this paper, including conception of the study, data identifi-335 cation and analysis, and drafting of the manuscript.

336
(a) Margin of each bus to system collapse corresponding to temperature projection (RCP 4.5 -blue line) from Fig. 7. (b) Margin of each bus to system collapse corresponding to temperature projection (RCP 8.5 -red line) from Fig. 7. Fig. 14: Margin to system failure for each node in the real distribution grid with respect to time. The climate change scenarios shown in Fig. 7 are used to estimate the temperature with respect to time. Using the estimated temperatures, the mean load consumption is predicted for all the consumers located in the modified real-world power grid. Finally using these predicted load consumption values, the margin to collapse scenario is calculated for each node using Φ from Methods. The color bar indicates the margin to system collapse at each node in the modified real-world power grid. As the (years) increases, the size of each data point in the figure is downscaled.