Hurricanes can cause devastating damage to overhead distribution lines leading to largepoweroutagesinelectricgrids. Poweroutagepredictionmodelscanhelputilities to plan for an expedited power recovery by identifying the extent of power disruptions before the arrival of a hurricane. These models often use multiple input parameters, including early warning forecasts of hurricane characteristics, and environmental and demographic information. We propose a quasi-binomial regression model to advance poweroutagemodelsandovercometheirexistinglimitations,suchasunboundedoutage predictions, limited extrapolation, and high uncertainties at low and high winds. This paper shows that the quasi-binomial model allows us to better capture the mechanics of power system failures due to hurricanes. We fitted our model to power outage data for four historical hurricanes, Harvey (2017), Michael (2018), Isaias (2020), and Ida (2021). WevalidatedourmodelfortheoutagesinFloridaduringHurricaneIan(2022). The quasi-binomial model outperformed existing random forest and negative binomial regression models with a 7%error versus 50%and76%, respectively. To demonstrate the quasi-binomial model’s good performance more comprehensively, we also tested a new beta regression model for outages. We show the quasi-binomial model had a smallercross-validationRootMeanSquaredErrorof0.23comparedto0.28forthebeta model. Finally, we show that our model also captures that grids with more redundant components can be more resilient to hurricane-caused outages. Thus, our proposed quasi-binomial model advances the state of the art for power outage predictions. INTRODUCTION Extreme weather events such as hurricanes can cause large-scale devastation to electric infrastructure (Smith 2020). Often, millions of customers lose power from hurricanes in the United States (US), e.g., (a) Hurricane Ida (2021) left 1.2 million 1 Arora, August 18, 2023 customers without power in Louisiana (AJOT 2021), (b) Hurricane Ian (2022) left more than 2.7 million customers without power in Florida (Florida Public Service Commission 2022; Giulia Carbonaro 2022; Cortes et al. 2022), (c) Hurricane Isaias (2020)causedmorethan3millionacrossfivestatesinUS(twitter.com2020;Aroraand Ceferino 2023a). These long-lasting power outages can cause permanent detrimental effectsforvulnerablecommunitiessincepowerisoftenrequiredtooperateothercritical infrastructures, e.g., hospitals, to respond to post-disaster emergency operations after disasters (Ceferino et al. 2020). The US Government has recognized the importance of security and resilience of critical infrastructure (Executive Office of the U.S. President 2022). For example, the US Senate passed the Grid Research and Security Act (2020) (Congress.gov 2020) to securethepowergrid. Additionally,theUSDepartmentofEnergy(DOE)hasalwaysput forwardthebenefitsofincreasingtheresilienceofpowergridinfrastructuretoextreme weather events (Office and August 2014). To achieve resilience, utilities must identify the vulnerabilities of the power systems to prepare the grid to withstand changing conditions, such as extreme weather events. However, they are often unprepared and struggletorespondtosuddenhurricane-inducedlargeblackouts,resultinginprolonged outageswhichhadhistoricallylastedonaverageforthreedays(Zimmermanetal.2017). Poweroutagemodelscanbecombinedwithearlywarninghurricaneforecasts(Can- gialosi et al. 2020) to provide utilities with outage and grid damage predictions three to seven days in advance. These predictions can be critical for resilience since they can inform the pre-impact planning and an early onset of coordination of repair crews fortherapidrecoveryfrompoweroutages(NationalAcademiesofSciences,Engineer- ing, and Medicine 2017). Thus, researchers have developed probabilistic and machine learning (ML) power outage prediction models to forecast outages before the arrival of a hurricane, e.g., from three hours to three days (Liu et al. 2005; Liu et al. 2007; Han et al. 2009a; Guikema et al. 2014; McRoberts et al. 2018). These models utilize informationonforecasthurricanewinds,environmentalconditions,demographics,and power system components to predict hurricane-induced outages. Previously, researchers have developed generalized linear models (GLM), general- izedadditivemodels(GAM),andrandomforest-basedpoweroutagemodels(Liuetal. 2005; Han et al. 2009a; Guikema et al. 2014). However, limited attention has been paid to understanding the compatibility of outages with Negative Binomial GAM and random forest models with the mechanics of power infrastructure failures (Arora and Ceferino 2023b). This paper presents a quasi-binomial GAM to predict power outages and overcome the above limitations. Unlike most previous models, the quasi-binomial distribution allows to model the fraction of customers without power in its appropriate range, from 0to100%. Additionally, the quasi-binomial model allows flexible variance modeling, i.e.,to represent high certainty in predictions for extremes at low and high winds. We also capture the effect of storm surges on power systems (Cruse and Kwasinski 2021) through a proxy for coastal flood exposure: distance of city to landfall (Pachev et al. 2023). We use data from Northeast, Southeast, and Southwest states on hurricane- induced power outages to make the power outage model generalizable across multiple cities at risk from hurricanes in the US. 2 Arora, August 18, 2023 LITERATURE REVIEW Liuetal.(2005)introducedanegativebinomialgeneralizedlinearmodel(GLM)to forecast storm-induced outages in North Carolina and South Carolina. They included four input parameters: maximum gust wind speeds, number of transformers, company indicator,andhurricaneindicator. Includingcompanyandhurricaneindicatorslimited the use of the forecast model to a particular region. Liu et al. (2008) also developed spatialgeneralizedlinearmixedmodels(GLMM)toaccountforanyspatialcorrelation in power outages but did not observe a significant improvement in the predictions. Later, Han et al. (2009b) worked on negative binomial GLM with outage data in the gulf coast and included extensive information on underground lines, overhead lines,switches,transformers,numberofcustomers,3-swindgusts,durationwindsover 20𝑚/𝑠,hurricaneindicators,pressure,timesincelasthurricanelandfall,precipitation, soil moisture, and land cover. Han et al. (2009a) also proposed a model without any storm or company indicator to make the model generalizable across the gulf coast region. The same group also developed negative binomial generalized additive models (GAMs) to forecast storm-induced power outages (Han et al. 2009a). GAMs showed improvement in outage predictions compared to GLMs, as GAMs can account for non-linearity between input variables and predicted outages. Guikema et al. (2010) and Quiring et al. (2011) used decision tree models with topologyandsoilparametersinformationtopredictoutages. Researchershaverecently usedothernon-parametricrandomforestmodelstohandlenon-linearityininputparam- eters(Guikemaetal.2014;McRobertsetal.2018;Shashaanietal.2018). Therandom forest model grows multiple parallel trees to reduce the variance and make predictions robust to outliers and noise (Breiman 2001). Guikema et al. (2014) used 3-s wind gusts, duration of winds over 20𝑚/𝑠, and population density to develop utility and hurricane-independentpredictionsofthefractionofcustomerswithoutpower. Nateghi et al. (2014) used the additional information on tree-trimming practices to improve the power outage predictions. However, the model had limited applicability since infor- mation on tree-trimming practices is generally not publicly available. Madrigano et al. (2015) included the information on tree species to further improve the power outage predictions with a random forest model. Tonn et al. (2016) predicted power outages at different confidence intervals using quantile random forests. Wanik et al. (2017) used lidar-derived tree data to develop another random forest power outage model. McRoberts et al. (2018) developed a two- stagemodeltopredictpoweroutages. Thefirststagewouldmakeabinarypredictionto identify outages or no outages in a region. The second stage would predict the fraction of customers without power for the regions with at least one outage. Shashaani et al. (2018)developedathree-stagepoweroutagemodelwiththefirststagetoidentifyoutage or no-outage; the second stage to classify outage regions into low, moderate, and high class; the third stage to predict the total number of outages. Researchers have also used artificial neural networks to predict outages (Xie et al. 2020). Haseltine and Eman (2017) used a neural network to predict the failure of indi- vidual power grid components before the arrival of a storm. However, high-resolution data on each grid component is generally not publicly available. Sun et al. (2016) used Twitter data to predict real-time outage, and Jaech et al. (2018) used repair logs 3 Arora, August 18, 2023 to forecast repair times. However, information from Twitter and repair logs are not availablebeforethearrivalofastormandcouldnotbeemployedforpre-stormplanning and coordination of resources. Arora and Ceferino (2023b) studied the limitations of the state-of-the-art power outagemodels. ThefirstlimitationwasthatnegativebinomialGAMscouldoverpredict power outages. For example, they showed negative binomial could predict as high as 16 times more outages than the actual number of customers. Second, random forest regressionhaslimitedcapabilitytoextrapolateforoutagesathighwindsduetosparsity of power failure data, i.e., predictions saturated at 70%of customers without power at wind speeds around 70𝑚/𝑠. Third, both the Negative binomial GAM and the randomforestdidnotcapturethemechanisticbehaviorofinfrastructurefailuresintheir uncertainty estimates. Negative binomial GAM showed high variance when predicted outages were close to 100% at high winds, and random forest showed high variances when predicted outages were close to 0% at low winds. However, we expect close to0%variance when outage predictions are close to 0% and 100%. For example, at very low winds, we expect little impacts to the power grid with high certainty, and at catastrophic wind levels, we expect large impacts also with high certainty. In addition, power system components near coastlines are heavily exposed to significant flooding, e.g.,due to storm surges (Cruse and Kwasinski 2021). Still, little attention has been given to incorporating this damage mechanism in power outage modeling. This paper presents the quasi-binomial outage model to overcome the aforemen- tioned limitations. We also included the proxy variable, the distance of a city to landfall, to model the damaging effects of storm surges on the power grid. DATA DESCRIPTION We included power outage data from five large historical hurricanes in developing the quasi-binomial regression power outage model. We acquired the power outage for Hurricane Harvey(2017) in Texas,Isaias (2020) inNew Jersey andNew York, andIda (2021) in Louisiana from PowerOuage, an independent organization that keeps track of power outages based on information from utilities in the US (PowerOutage 2022). In addition, we obtained the power outage data for Hurricane Michael (2018) and Ian (2022) from Florida Public Service Commission (2022). We used Hurricane Harvey, Michael, Isaias, and Ida’s data for model calibration and Ian’s data for validation. For training, we filtered outliers as PowerOutage’s reports are approximated in a few cases (PowerOutage2022). Weusedtheinterquartilerangecriteriaontheratioofhouseholds and the number of customers without power to remove 255 out of the 2577 total data points (cities) from Hurricane Harvey, Michael, Isaias, and Ida for calibration. Poweroutagepredictionsrelyoninputparameters,includinghurricane,environmen- tal, power system, and socio-demographic features. The input parameters are typically available at different scales. Thus, we adjusted the input data (e.g., aggregated or in- terpolated) to match the same spatial resolution of the outage data, as in our previous research (Arora and Ceferino 2023b). We developed the power outage model at the city-level spatial resolution consistent with the recent power outage models. Future research could develop an outage model with a finer resolution if input parameters and historical power outages at a finer spatial resolution are available. 4 Arora, August 18, 2023 Intotal,thisstudyuseddataon4.62millionoutagesin2322citiesacrossfourstorms covering 11 utilities in New York and 5 in New Jersey, five in Florida, 39 in Texas, and nine in Louisiana. Figure 1 shows the number of outages across New Jersey after Hur- ricaneIsaias(2020)todevelopthequasi-binomialmodel. Additionally,Supplementary Figures S1, S2, S3, and S4 show the power outages in Texas during Hurricane Harvey (2017), Florida during Hurricane Michael (2018), New York during Hurricane Isaias (2020), and Louisiana during Hurricane Ida (2021). We further describe the response and input variables, also summarized in Table 1. Response Variable We developed the power outage model to predict the fraction of customers without poweratthecitylevel. Wechosethefractionofcustomerswithoutpowerastheresponse variable to quantify the impact of power outages across cities having variable numbers of customers. Fig. 1.Distribution of power outages at quantile breakpoints in the aftermaths of Hurricane Isaias (2020) in New Jersey, as reported by PowerOutage Input Parameters Researchers have used different explanatory variables to predict hurricane-induced power outages (Liu et al. 2005; Han et al. 2009b; Han et al. 2009a; Guikema et al. 2010; Guikema et al. 2014; Shashaani et al. 2018). The explanatory variables include hurricane winds, whose intensity and extent can determine the scale of damage to power systems (Bjarnadottir et al. 2013). Other environmental conditions ( e.g.,land use,precipitation,soilmoisture)caninformaboutmorevulnerablesectionsofthepower 5 Arora, August 18, 2023 grid. For instance, power poles or trees in wet soil conditions are more susceptible to being torn out from soil under high wind conditions (McRoberts et al. 2018). Key power infrastructure data ( e.g.,the number of transformers) can also inform on the redundancy of power systems. In grids with redundant components, people are more likely to access electricity even if a component fails, e.g., from strong hurricane winds (Brown 2002). Finally, socio-demographic information ( e.g.,population density) can informonthecustomerdensityandpotentiallayoutofthepowergridinurbanandrural areas. Thus, we divided the input parameters into four categories: hurricane features, environmental features, power system information, and socio-demographic features, which are described next. Hurricane Features We included 3-s gust wind speed and duration of strong winds, i.e., duration of windsover 20𝑚/𝑠consistentwithpreviousresearch(Liuetal.2005;Hanetal.2009b; Hanetal.2009a;Guikemaetal.2010;Guikemaetal.2014;Shashaanietal.2018). We usedBestTrackdataforthehurricanepath(Knappetal.2010)andacompletetropical cyclonewindfieldmodelbyChavasetal.(2015)todeterminetheaxisymmetricwinds and background wind model byLin et al. (2012) to obtain the complete wind structure. Often, background winds are neglected, which can lead to underestimated maximum winds. For example, we calculated that the maximum 3-s gust wind increased from 37.5𝑚/𝑠to49.5𝑚/𝑠(32%more) in New Jersey after including background winds during Hurricane Isaias (2020). We consider the 3-s gust wind speed and duration of strong winds at the centroid of a city. We used the log of 3-s gust winds as an input parameter to better model power outages. We show the 3-s gust winds in New Jersey during Hurricane Isaias in Figure 2. Historically, apart from high hurricane winds, significant storm surges can also impact the power systems. Poudyal et al. (2022) showed that the compound effect of hurricane winds and storm surge could lead to more power loss compared to damages by hurricane winds alone in Texas. Not only a hurricane can cause a significant storm surge when it transitions from sea to land. A hurricane can also cause flooding when transitioningfromlandtoseaduetohighwaves(Buccietal.2022). Thus,weincluded the minimum distance between the city and the points where the hurricane crosses the coastal line as an input parameter to capture the compound risk to power systems from hurricane winds and storm surges. We normalized the distance with the radius of maximum winds to account for different hurricane sizes. Other Environmental Features We included various environmental features for outage modeling, e.g., land cover information,topography,soilconditions,andtreedensity. Theseenvironmentalfactors can help enhance outage predictions as described in the previous research (Liu et al. 2005; Han et al. 2009b; Han et al. 2009a; Guikema et al. 2010; Guikema et al. 2014; Shashaani et al. 2018; Arora and Ceferino 2023b). WeacquiredlandcoverdataavailableasNationalLandCoverData(NLCD)fromthe US Geological Survey (USGS) on Multi-Resolution Land Characteristics Consortium (Jin et al. 2021). The land cover classes can inform about distribution grid layouts, resulting in different outage patterns. For example, rural areas mostly have a radial 6 Arora, August 18, 2023 Fig. 2.Distribution of 3-s gust wind speeds during Hurricane Isaias (2021) in New Jersey. The wind gusts include both circulating and background components of the winds. layoutforthepowergridwhereeachcustomerisconnectedtoasinglefeeder. Thus,the failure of a single feeder in a radial grid leads to the shutting down of all downstream lines until the repair is completed for downed power lines, which can result in more outagepropagationduetofailedpowersystemcomponentsthaninurbanareas(Petersen 1982). USGSoriginallyclassifiedNLCDdatainto20classes,availableinrasterformat withcellsofsize 30𝑚×30𝑚. Wereclassifiedtheoriginal20NLCDclassesintonine major classes: developed, water, barren land, forest, scrub, grasslands, pasture land, cultivated cropland, and wetlands. We used following pythonpackages: rasterio,gdal, shapely,fiona,geopandas ,zonal_statistics to perform calculations on raster data. We calculated the percent of area occupied by each class of land cover in a city. Han et al. (2009b), Nateghi et al. (2014), McRoberts et al. (2018) have highlighted the importance of including precipitation and soil moisture conditions in power outage predictions. Many distribution lines are close to the trees, posing an imminent threat to failure of distribution lines from trees falling due to strong hurricane winds. The deviations in precipitations and soil moisture from standard conditions might increase the likelihood of trees uprooting from high winds resulting in more outages. Trees are morelikelytouprootduetostronghurricanewindsunderwetsoilconditions(Hanetal. 2009b; Nateghi et al. 2014). Additionally, persistent dry conditions can create gaps in soil,weakeningtreerootsandincreasingthechanceoftreeuprootingfromhighwinds. We acquired the information on soil moisture and precipitation at an hourly scale 7 Arora, August 18, 2023 for the past 42 years, starting in January of 1971 from Xia et al. (2012). Xia (2012) developed the Variable Infiltration Capacity (VIC) model available on National Land Data Assimilation System Phase 2 (NLDAS2). NLDAS2 data is available in a grid pattern with points spaced apart 0.125◦×0.125◦. We obtained the precipitation and soil moisture at the centroid of a city using nearest neighbor interpolation. We extracted the soil moisture available for three depths: 0−10𝑐𝑚,10−40𝑐𝑚, and40−100𝑐𝑚from NLDAS2 data. We converted the hourly soil moisture to daily soil moisture by taking an average of hourly soil moisture over a day. Soil moisture mayvarysignificantlyfromonelocationtoanotherbecauseofambientconditions. We normalized the soil moisture by computing the percentiles at each location after fitting the daily soil moisture data to Pearson Type III distribution using maximum likelihood estimate (MLE). Then, we evaluated the percentile of soil moisture a day before the storm would hit the city, denoted as CDF1: percentile of soil moisture for 0−10𝑐𝑚 depth, CDF2: percentile of soil moisture for 10−40𝑐𝑚depth, and CDF3: percentile of soil moisture for 40−100𝑐𝑚depth. Precipitation conditions are represented as standard precipitation index (SPI) (Wu et al. 2007; Guttman 1998; Casey 2016). SPI can represent the wet and dry conditions before a storm hits. We calculated the SPI for one month (SPI1), three months (SPI3), six months (SPI6), and 12 months (SPI12). SPI is obtained from precipitation time series data in three steps. First, we fitted the time series data to Pearson Type III distribution. Second, we calculated the percentile of precipitation for the month before the storm’s impact. In the final third step, we obtained SPI using the inverse of cumulative standard normal distribution on the percentiles computed in the previous step. We also considered the expected precipitation after the hurricane as an input parameter, as heavy precipitation can lead to weakened roots and more trees uprooting (McRoberts et al. 2018). We included the percentage of trees in a city to capture the falling tree hazard on electric lines and poles. We obtained the area covered by trees from the National Insect and Disaster Risk Maps (Krist Jr. et al. 2014) created by the US Department of Agriculture (USDA) available at a resolution of 240𝑚×240𝑚. We use zonal statics in ArcGIS to determine the total area covered by trees in a city. The depth of the soil from which plants and trees can effectively derive water and nutrients for their growth is defined as the root zone (RZ) depth. Trees with more RZ depth can withstand strong hurricane winds (McRoberts et al. 2018). We obtained the RZ data from USDA under Gridded Soil Survey Geographic (Soil Survey Staff 2021) as raster data with a resolution of 30𝑚×30𝑚. We used zonal statistics in ArcGIS to obtainthemeanRZdepthforacity. WeincludedtheRZdepthasaninputparameterto poweroutagepredictivemodelsasRZdepthcouldadditionallyindicatefallinghazards of trees on distribution lines and poles from strong hurricane winds. Chapman (2000), Miller et al. (2013), Guikema et al. (2010), Quiring et al. (2011), McRobertsetal.(2018)havefoundvariationinhurricanewindspeedswithvaryingsur- facetopography. Therefore,acity’smeanandmedianelevationaretypicaltopographic input parameters in outage predictions. We obtained mean and median elevation from a Digital Elevation Model (DEM) at a 30𝑚resolution in the Global Multi-Resolution Terrain Elevation Data (GMTED2010) by USGS (Danielson and Gesh 2011). 8 Arora, August 18, 2023 Feature Abbreviation Data Source Outages Outages (PowerOutage 2022) log of 3-s Gust Wind Speed logvmaxa,b(Chavas et al. 2015; Knapp et al. 2010) Duration of Strong Winds Duration Distance to Hurricane crossing coast Normalized with Radius to Maximum windsDistHurra,b Percent Developed Area Developed (Jin et al. 2021) Percent Water Area Water Percent Barren Area Barrena Percent Forest Area Forest Percent Scrub Area Scrub Percent Grassland Area Grasslanda Percent Pasture Area Pasture Percent Crops Cultivated Area Crops Percent Wetlands Area Wetlandsa,b Standard Precipitation Index 1 month SPI1 (Xia et al. 2012) Standard Precipitation Index 3 months SPI3 Standard Precipitation Index 6 months SPI6a,b Standard Precipitation Index 12 months SPI12 Soil Moisture 1𝑠𝑡Layer CDF1a,b(Xia et al. 2012) Soil Moisture 2𝑛𝑑Layer CDF2 Soil Moisture 3𝑟𝑑Layer CDF3a 7 day precipitation precipa,b(Xia et al. 2012) Root Zone Depth Rzonea(Soil Survey Staff 2021) Percent Treed Area Treesa,b(Krist Jr. et al. 2014) Mean Elevation Mean_Elea(Danielson and Gesh 2011) Median Elevation Median_Ele Substations Substations (data.gov 2022) Substation Normalized to Population Sub_Norma,b Road Length Normalized to Substations R_over_S Population Density Pop_Dena,b(American Community Survey 2019) Table 1. Parameters to build the power outage prediction models: all variables are rescaledatthecitylevel. Parametersaregroupedintocategoriesseparatedbyhorizontal lines. We selected one variable from each category from each group to minimize correlation across parameters. a-variables selected after filtering of correlated features b-finally selected variables after 5-fold cross-validation Power System Information Previously, Liu et al. (2005), Liu et al. (2007), Han et al. (2009a), Han et al. (2009b), Guikema et al. (2010), Quiring et al. (2011) have included information on power systems, such as density of transformers, type of protection device, number of 9 Arora, August 18, 2023 poles,numberoftransformers,lengthofoverheadlines,andlengthoftheunderground line, in outage predictions from hurricanes. The low number of protective devices and transformers, can indicate sparsity and low redundancy in power systems. The sparsity in power system protection devices, e.g., sectionalizing switches to separate the damaged part of power lines, profoundly affects power delivery to customers. The failureofoneprotectiondeviceinasparsepowernetworkcanleadtolargerandlonger power interruptions for many customers (Brown 2002). Also, National Academies of Sciences, Engineering, and Medicine (2017) reports that overhead distribution lines are the most vulnerable components of a power grid to strong hurricane winds. Power systems with many wood poles are particularly at high risk as these poles are generally designed to moderate wind speeds only up to 20𝑚/𝑠(IEEE 2007). Thus, including this information about power systems can help better capture the propagationofpoweroutages. However,duetosecurityconcerns,theinformationonthe powersystem’svariablesisnotgenerallypubliclyaccessible(StatesandAccountability 2020). Toovercomethisproblem,Birchfieldetal.(2017),Pahwaetal.(2014),Schultz etal.(2014),Schweitzeretal.(2010),Valenzuelaetal.(2019),Pisanoetal.(2019),Zhai etal.(2021)proposedgeneratingsyntheticpowergridsusingopenpublicgeoreferenced datae.g.,roads,buildings. Yet,formulatingandimplementingtheseprocedurescanbe a separate study. To keep the analysis simple and still include key power system information in our model,weadoptedafewsimplifiedtechniquesofsyntheticgridgeneration. Totrackthe numberofpolesexposedtohurricanewinds,weapproximatedthelengthofdistribution lines with the length of local roads in a city since most lines are along roads (Winkler etal.2010;Valenzuelaetal.2019). Inaddition,wegatheredinformationonthelocation and number of transmission substations from the Homeland Infrastructure Foundation Level Database (data.gov 2022) to use it as a proxy for the power systems’ sparsity or redundancy. To illustrate why a small number and density of substations can be an important indicator of lack of redundancy, Figure 3 shows the power flow in a radial power grid. The power flow in a grid starts at generation substations. Then transmission lines deliver power to transmission substations. Sub-transmission power lines deliver power to commercial customers or distribution substations. Primary distribution lines deliver powertodistributiontransformers,andthesecondarydistributionsystemdeliverspower to feeders. Finally, each customer is connected to a single feeder for power delivery. The failure of a feeder will result in a power outage for all downstream customers. According to the flow of power supply in a radial grid as depicted in Figure 3, more transmission substations would indicate more distribution transformers (Brown 2002). Thus,weincludedthenumberofsubstations,thenumberofsubstationsnormalized with population, and the number of transformers per unit road length (an indicator of power lines) in the power outage prediction model. Since some cities have no transmission substation, we calculated the variables representing power systems at a coarser resolution of the county level and assigned their values to all the cities within the county. 10 Arora, August 18, 2023 Fig. 3.Depiction of the flow of electricity in a radial power grid. First, power is deliveredfromagenerationsource( e.g.,coalplant)toatransmissionsubstation. Then, high-voltage transmission power lines transmit power to the distribution substation. Finally,powerisdeliveredtocustomersthroughdistributionpowerlines. Eachcustomer receivespowerfromasinglefeeder. Inasparseradialgridwithfewerfeeders,thefailure of a single feeder would result in outages for multiple downstream customers. Pictures for power plant, substation, and overhead powerlines are extracted from (USGS 2004), (Emergible Consulting Services, Accessed Online: 08/10/2023 ), and (JustEnergy, Accessed Online: 08/10/2023 ), respectively. Density Information Population density could also indicate the extent of lines and poles exposed to hur- ricane winds (Arora and Ceferino 2023b). Rural areas with a low population density generallyhavesparselylocatedpowergridcomponentswherethefailureofonecompo- nentcouldresultinapoweroutageformultiplecustomerscomparedtourbanareaswith redundant power components (Brown 2002). To capture this effect, population density has been used for power outage modeling (McRoberts et al. 2018). Thus, we also used populationdensityfrom5-yearestimatesinAmericanCommunitySurvey(2019)asan input parameter. POWER OUTAGE MODELING We develop the power outage model to predict the fraction of customers without power in a city. Beta, and binomial regressions (Olkin and Liu 2003) have been used 11 Arora, August 18, 2023 in many fields, ranging from financial to ecological studies (Yee 2012; Olkin and Liu 2003; Douma and Weedon 2019), when the response variable is a fraction distributed fromzerotoone,suchastheproportionofcustomerswithoutpowerintheaftermathof ahurricane. Unlikethelinearregressionmodel,thebeta,binomial,andquasi-binomial regression models do not assume homoscedasticity i.e.,constant variance for output variable with change in input parameters, and can model the variable dispersion in power outage predictions across different cities. Previously, researchers have applied methodsotherthanlinearregressiontomodelpoweroutages. Still,thosemodelsfailto capture the mechanics of power infrastructure failures as stated in Arora and Ceferino (2023b). Thus, we explore beta regression and quasi-binomial regressions as models to en- hance power outage predictions, better capturing the physics of power system failures. We present the quasi-binomial instead of the binomial regression because it better models large overdispersion, a feature extensively reported to be critical previously in outage prediction (Dunn and Smyth 2018; Han et al. 2009a). Also, notice that the quasi-binomial regression is a hierarchical model of the binomial regression; thus, it inherits its modeling strengths. Beta Regression Sincethefractionofcustomerswithoutpower 𝑦inacityisacontinuousvariablewith values ranging from zero to one, we can model 𝑦as beta distributed. The probability density function is given by 𝑓(𝑦|𝛼,𝛽)=𝑦𝛼−1(1−𝑦)𝛽−1 𝐵(𝛼,𝛽)(1) where 𝐵(𝛼,𝛽)=Γ(𝛼)Γ(𝛽) Γ(𝛼+𝛽)(2) and𝛼and𝛽are the parameters to be estimated, and Γ(.)is the gamma function. Correspondingly,themeanandvarianceforthefractionofcustomerswithoutpowerare 𝐸[𝑦]=𝛼 𝛼+𝛽(3) 𝑣𝑎𝑟[𝑦]=𝛼𝛽 (𝛼+𝛽)2(𝛼+𝛽+1)(4) Generally, the beta distribution is represented with an alternative form with param- eters𝜇and𝜙so that the mean 𝐸[𝑦]and variance 𝑣𝑎𝑟[𝑦]of beta distribution are given as: 𝐸[𝑦]=𝜇 (5) 𝑣𝑎𝑟[𝑦]=𝜇(1−𝜇) 1+𝜙(6) 12 Arora, August 18, 2023 This representation is known as mean-precision parameterization, where the previ- ously defined 𝜇is the mean, and 𝜙is the precision parameter (inverse of dispersion). Thus, we can rewrite the probability distribution in Eq. 1 with 𝜇and𝜙as 𝑓(𝑦;𝜇,𝜙)=Γ{𝜙} Γ{𝜇𝜙}Γ{(1−𝜇)𝜙}𝑦𝜇𝜙−1(1−𝑦)(1−𝜇)𝜙−1(7) where 0< 𝜇 < 1, and𝜙 > 0. The wide range of 𝜇and𝜙provides excel- lent flexibility to fit different probability distribution shapes, e.g.,skewed, symmetric, bell-shaped. Beta regression parameters are determined using Maximum Likelihood Estimates (MLE) similar to Generalized Linear Models (GLMs) (Dunn and Smyth 2018). The mean parameter 𝜇is related to the input parameters via a link function 𝑔(𝜇)=𝛽𝑋 (8) where𝛽is the vector of the unknown parameters determined using MLE, and 𝑋is thematrixwithinputparameters. Thelinkfunction 𝑔(𝜇)mustbemonotonicandtwice differentiable,and 𝜇isboundedbetween 0and1. (FerrariandCribari-Neto2004). The most common choice for the link function is a logit function. 𝑔(𝜇)=𝑙𝑜𝑔𝜇 1−𝜇 ;𝜇=𝑒𝛽𝑋 1+𝑒𝛽𝑋(9) The dispersion parameter 𝜙can also be linked to 𝑋similar to𝜇in Eq.8. However, for the scope of this paper, we keep 𝜙constant. Notice that under this assumption, the beta distribution still provides enough flexibility to fit different probability distribution shapes (e.g.,skewed, symmetric, bell-shaped). Finally, the log-likelihood function is maximizedtoestimateunknownparameters 𝛽and𝜙for𝑛observations. Wemaximized the log-likelihood for the outage data obtained for Hurricane Harvey (2017), Michael (2018), Isaias (2020), and Ida (2021). So far, the beta distribution proposed in Eq. 7 and Eq. 8 is defined on an open interval, i.e., 0< 𝑦 < 1. However, the fraction of customers without power in a city couldalsobeequaltoexactlyzeroorone. Fittingexactvaluesofzerosandonesleadsto poor estimates of 𝛽and𝜙. Previously researchers have introduced ways to model such boundary conditions for beta regressions (Smithson and Verkuilen 2006; Ospina and Ferrari 2012; Liu and Kong 2015; Swearingen et al. 2012). For example, Ospina and Ferrari (2012) nudged the output variables so that 𝑦would never reach values of zero and one. Furthermore, Ospina and Ferrari (2012), Liu and Kong (2015), Swearingen et al. (2012) developed zero, and one-inflated beta regressions to model the boundary conditions separately, but such models require additional parameters to be determined. Instead,weincludeaslackvariable 𝜖toboundthefractionofcustomerswithoutpower between open intervals (Smithson and Verkuilen 2006). ˜𝑦=  𝑦+𝜖 𝑖𝑓 𝑦=0 𝑦 𝑖𝑓 0< 𝑦 < 1 𝑦−𝜖 𝑖𝑓 𝑦=1(10) 13 Arora, August 18, 2023 We chose𝜖=10−8to ensure no significant changes in the accuracy of predicted outages. For example, New York City (NYC) is the most populous city in the US, with 7.8 million people. If NYC has no outages, adding a value of 10−8to the fraction of customerswithoutpowerwouldstillresultinlessthanoneoutage,whichisroundedoff to zero outages. The GLM-like formulation in Eq. 8 models a linear relationship between the logit of the mean of the fraction of customers without power and input parameters. However,thisrelationshipcanbenon-linear,asdemonstratedinpreviousresearch(Han et al. 2009a), which can be modeled with Generalized Additive Models (GAMs), an extension of GLMs. GAMs formulation using smoothing functions can capture this non-linear relationship. 𝑔(𝜇)=𝛽0+𝑝∑︁ 𝑗=𝑖𝛽𝑗𝑓𝑗(𝑥𝑗) (11) where𝑔(𝜇)is the link function in Eq. 9, 𝛽0is the intercept, and 𝑓𝑗are smoothing splines functions which are fitted to each parameter 𝑗, where𝑗∈[1,...,𝑝], and𝑝is thetotalnumberofinputparameters. Themostcommonexamplesofsmoothingsplines include P-splines, B-splines, and S-splines. Splines for each parameter could be fitted with polynomials of any degree. Thus, we modeled all parameters as polynomials of theS-spline,exceptwindgust(logvmax)anddistanceofthecitytolandfall(DistHurr) (Table 1). We modeled wind gusts with a polynomial of degree one because higher winds are expected to increase damage monotonically, e.g., zero damage at zero wind gusts and complete collapse of power systems with wind gusts from hurricanes of category 5(>90𝑚/𝑠). Also, cities closer to the point of the hurricane-crossing land tendtoexperiencehigherstormsurges. Hence,wealsomodeledtheparameterDistHurr as a degree one polynomial. We used mgcvlibrary in the R program to perform MLE estimates of unknown parameters (Wood 2017). Quasi-Binomial Regression The binomial regression is another popular class of regressions to model fractions arisingfromcountdata,suchasthefraction 𝑦ofoutagesinacitywith 𝑚totalcustomers. Quasi-binomialregressioninheritsthesebinomialregression’smodelingstrengthsand, inaddition,canhandleoverdispersionintheoutagepredictionsbetter(DunnandSmyth 2018; Consul 2010). Quasi-binomial distribution is a simplified form of beta-binomial distribution. The beta-binomial distribution assumes that the mean rate of failure ( 𝜇in Eq. 12)isbetadistributed(DunnandSmyth2018). Thus,theprobabilitymassfunction of beta-binomial distribution is given as 𝑃(𝑦;𝑚,𝛼,𝛽)=𝑚 𝑚𝑦𝑦𝛼−1(1−𝑦)𝛽−1 𝐵(𝛼,𝛽)(12) where𝐵()is the function from Eq. 2, and 𝑚is the number of customers in a city. The parameters 𝛼and𝛽are determined using MLE. Unlike beta regression, MLE for a quasi-binomial regression is weighted based on the number of customers in a city. Correspondingly,themeanandvarianceforthefractionofcustomerswithoutpowerare 14 Arora, August 18, 2023 𝐸[𝑦]=𝛼 𝛼+𝛽(13) 𝑣𝑎𝑟[𝑦]=𝑚𝛼𝛽(𝛼+𝛽+𝑚) (𝛼+𝛽)2(𝛼+𝛽+1)(14) Similartobetadistribution,wecanparametrizethebeta-binomialintermsof 𝜇and 𝜙, resulting in the typical form of the quasi-binomial distribution (Dunn and Smyth 2018). In that case, 𝐸[𝑦]=𝜇and𝑣𝑎𝑟[𝑦]=𝜙𝜇(1−𝜇). Here,𝜙is the dispersion parameter that allows handling varying overdispersion in the proportions of outages. Quasi-binomialregressioncanmodelthefractionsinaclosedinterval,i.e.,including zeros and ones. Thus, we do not need to apply the transformation in Eq. 10. Similar to the beta model, we used GAMs representation from Eq. 11 and logit link function from Eq. 9 for the quasi-binomial model. We use mgcvlibrary in the R program to fit the quasi-binomial model for the fraction of customers without power. FEATURE SELECTION Initially, we considered a large number of input parameters, which can challenge developing meaningful regression models (Kira and Rendell 1992). High-dimensional input data can result in degraded performance with overfitting, i.e.,small errors for training data but high error on test data. Thus, feature selection is crucial in model development,especiallyforsparsedatasetssuchasthecasefornaturalhazards(Kumar and Minz 2014). Thus, we selected the most appropriate and sufficient number of featurestopredictoutagesaccurately(Provost1999),asinpreviousresearchforoutage modeling (Han et al. 2009b). Feature selection also helps address multi-collinearity, i.e., when highly correlated input features introduce spurious causation because of similar importance coefficients in the regression. One method used in previous research to avoid multi-collinearity is Principal Component Analysis (Han et al. 2009b; Han et al. 2009a; Markhvida et al. 2018). However, we did not use it to make the model generalizable by keeping the featuresintheoriginalspace,whichalsoenabledustostudytheeffectsofkeyhurricane characteristics, such as winds, on power outages. Instead,wefilteredoutthehighly-correlatedinputfeatures. WeprovidethePearson correlation𝜌for all the variables in Supplementary Figure S5. First, we selected only one feature with highest 𝜌out of sets of highly-correlated input features with Pearson correlation coefficient ( 𝜌)>0.5, from Table 1. Then, we filtered out the input features with a correlation to the output (fraction of customers without power) below 0.1 to keep only the most important variables. For example, we filtered out the variable "Substations" as it had |𝜌|<0.1with fraction of outages, but we still kept Substations normalizedtopopulationssince |𝜌|=0.23. So,ourmodelwillstillincludetheeffectsof powersystemvariablesonoutages. Afterfilteringofvariables,weremoved14variables (out of 28) from further model development. The remaining variables for subsequent analysisare1)logvmax,2)DistHurr,3)Barren,4)Grassland,5)Wetlands,6)SPI6,7) CDF1,8) CDF3,9)precip,10) Rzone,11)Treed_area, 12)Mean_Ele,13)Sub_Norm, 14) Pop_Den. 15 Arora, August 18, 2023 MODEL SELECTION We first conducted a qualitative assessment of the performance of beta and quasi- Bionmial GAM models. We computed their standardized quantile-residuals to assess modelqualitybyevaluatinghowclosetheyaretoastandardnormaldistribution(Dunn andSmyth1996;Pereira2019;KlarandMeintanis2012). Figure4showsQ-Q(quantile- quantile)plotsforthequantileresidualsinbothmodels. Readerscanfindtheinformation on calculating quantile residuals in the supplementary material. We observe high deviationsforbetaregression’squantileresidualsatthetails,suggestingthatthismodel cannot account for extreme outages. We hypothesize these deviations are high because betaregressionsarenotsuitedtoincludezerosandones,asdiscussedearlier. Incontrast, quantile residuals of binomial regression show a closer match to the standard normal distribution. We hypothesize that quasi-binomial regressions perform better because they can handle predictions across cities of various sizes, including zeros and ones. Quasi-binomial regressions assign different weights according to the cities’ number of customers, unlike the beta regressions that assign the same weight for all (Douma and Weedon 2019). Fig. 4.Q-Q plots for quantile residuals, A: beta regression model, B: quasi-binomial regression model. Quantile residuals should be close to standard normal distribution for a good fit. However, plot A for beta regression shows high deviations for quantile residuals from standard normal distribution suggesting a poor fit for power outage modeling. QuantileresidualsinplotBforquasi-binomialregressionhaveadistribution closer to the standard normal distribution, suggesting a better fit for power outage modeling. Next, we conducted a quantitative comparison of the models’ performance through 5-fold cross-validation (Hastie et al. 2002), comparing the mean cross-validation root mean squared errors (RMSE) (Entekhabi et al. 2009). RMSE is given by 𝑅𝑀𝑆𝐸=√︂ 1 𝑛∑︁ (𝑦−ˆ𝑦)2 (15) 16 Arora, August 18, 2023 where𝑦is the observed fraction of customers without power, ˆ𝑦is the predicted fraction of customers without power, and 𝑛is the total number of observations. We obtained a 22%lower cross-validation RMSE for the quasi-binomial than for the beta regression (0.23 versus 0.28), further showing that the quasi-binomial regression can better model power outages. Hence, we use it for the rest of the paper. Finally, we filtered out additional variables to control for overfitting on the quasi- binomial regression. We fitted different regressions with the variables selected previ- ously and then dropped one variable at a time. We performed 5-fold cross-validation (Hastie et al. 2002) and compare the mean cross-validation RMSE for all the models (Figure 5). To obtain the feature importance, first, we subtracted the RMSE of each model and the RMSE of the base model with all variables from feature selection. Re- movingthemostimportantfeaturewillresultinthehighestincreaseinRMSE.Finally, we normalized the RMSE with the largest difference across all the models. Fig.5.Featureimportanceforallthevariablesselectedafterfilteringinfeatureselection. Distancetothelocationofahurricanecrossingthecoastlinenormalizedbytheradiusof maximumwindsandthelogarithmofmaximumwindsobservedarethemostimportant features to determine power outages. Features with red colored bars were filtered out as removal of these features resulted in reduced cross-validation RMSE. Features with blue colored bars were selected for further analysis. We show that distance to the hurricane crossing the coastline(DistHurr) and wind gusts(logvmax)arethemostessentialvariablesastheirexclusionresultedinanRMSE increase of 8.4% and 24.6 % (Supplementary Table S1). The increase in RMSE after dropping DistHurr and logvmax underlines the importance of these variables directly correlated to the significant damage that winds and storm surges can cause to power grids. Incontrast,excludinggrasslandarea,rootzone,meanelevation,andsoilmoisture 17 Arora, August 18, 2023 forthedepthof60-100cmresultedinlowerRMSEthanamodelwithallthevariables. Thus, we removed them. Previous studies have shown that the reliability of power systems can be lower in inland water areas (Ankit et al. 2021), as water can corrode the electric equipment, and the flow of debris coming from disposals during the hurricane can further damage electricalequipment. Weconfirmedthesearerelevantfeatures,asexcludinglandcover parameter wetlands leads to a higher estimate of RMSE. Also, population density and systeminformationvariablescanindicateurbanareaswithmorepowersystemresources (Brown 2002). We also observed that these are critical variables in our analysis, as removing those two parameters led to a higher RMSE. We observed an RMSE value of 0.235 for 5-fold cross-validation, filtering out grassland area, rootzone, mean elevation, and (Supplementary Table S1). We obtained 261.5 for 𝜙, which represents an overdispersion in outage data, a value that is consistent with previous research (Liu et al. 2005; Han et al. 2009a; Han et al. 2009b; Arora and Ceferino 2023b). The final selected variables are 1) logvmax, 2) DistHurr, 3) Barren, 4) Wetlands, 5) SPI6, 6) CDF1,7)precip,8)Treed_area,9)Sub_Norm,10)Pop_Den. Alltheselectedvariables are significant at a p-value of 0.01. VALIDATION OF PROPOSED MODEL We further validated the quasi-binomial model with outage data in Florida from Hurricane Ian (2022) (Florida Public Service Commission 2022). In addition, we alsocomparedthepredictionsofquasi-binomialregressiontotwootherstate-of-the-art modelsforoutagepredictions: negativebinomialregression(Liuetal.2005;Hanetal. 2009a; Han et al. 2009b) and random forest model (McRoberts et al. 2018; Shashaani et al. 2018). Figure 6, presents the predictions with quasi-binomial, negative binomial, random forest, and actual outages in Florida from Hurricane Ian (2022). Since the resolution of our developed model is at the city level, and the actual outages report is available at the county level, we aggregated the predicted outages at the county level to compare them with the actual outages. For quasi-binomial and random forest regression models, we obtained the mean fraction of customers without power by estimating the population-weighted mean of the percentage of outages at the city level. We obtained the predictions for negative binomial as the count of outages. Hence, for negative binomial predictions, we summed the outages at the city level to obtain the total outages at the county resolution and calculated percentage outages by taking the ratio with total customers in a county. We report the total outages predicted with quasi-binomial, random forest, and negative binomial regression in Table 2. We also reported the percent error on the total outage counts estimated as 𝐸𝑟𝑟𝑜𝑟 𝑜𝑛𝑡𝑜𝑡𝑎𝑙 𝑜𝑢𝑡𝑎𝑔𝑒𝑠 = Í𝑚𝑖ˆ𝑦𝑖−Í𝑚𝑖𝑦𝑖Í𝑚𝑖𝑦𝑖 ×100 (16) whereÍ𝑚𝑖ˆ𝑦𝑖are the total outages predicted by a regression model andÍ𝑚𝑖𝑦𝑖are the total outages in Florida from Hurricane Ian (2022). We also reported the RMSE on the fraction of customers without power in counties. The quasi-binomial model predicted the outages with the least error ( 7.04%error on total outages and RMSE of 0.13 on the fraction of outages) compared to Negative Binomial ( 61.2%error on total 18 Arora, August 18, 2023 Fig.6.Predictedoutages(a,b,andc),andactualoutages(d)inaftermathsofhurricane Ian (2022) in Florida. a.Quasi-binomial regression model predicted outages with an RMSE of 0.13 for percent outages in a county. b. Random forest regression model predicted outages with an RMSE of 0.15 c. The negative binomial regression model predicted outages with an RMSE of 0.28. outages and RMSE of 0.28) and Random Forest Regression models ( 49.7%error on total outages and RMSE of 0.15). Thus, the quasi-binomial regression showed the best model performance. Still, this 19 Arora, August 18, 2023 proposed model can be improved as we observed outage over-predictions in Osceola, BrevardCounty,andnearbyareas. AfterHurricaneIrma(2018),utilitieshavehardened the grid by undergrounding the vulnerable distribution lines (Pallone 2020). However, ourmodelhasmainlybeenfittedtopoweroutagesinregionswithoverheadlines. Inthe future,researcherscanenhancepoweroutagepredictionwithmorespecificinformation on the deployment of underground lines. Regression Model TotalOutages RMSE Error on Total Outages Quasi-Binomial 2,886,413 12.32% 7.04% Random Forest 4,649,134 14.85% 49.72% Negative Binomial 735,768 28.12% 76.30% Actual Outages 3,105,294 - - Table2.ReportedResultsfordifferentmodels. Quasi-binomialregressionmodelbetter predictedtheoutagesacrossFloridaduringHurricaneIan(2022)comparedtoRandom Forest and Negative Binomial regression models. Percent error on total outages and RMSEonpercentoutagesinacountyarecalculatedusingEqs. 16and15respectively. OVERCOMING LIMITATIONS OF EXISTING OUTAGE MODELS AroraandCeferino(2023b)discussedthelimitationsoftheexistingstate-of-the-art poweroutagemodels. Here,weshowthatthequasi-binomialmodelcanovercomethose limitations. Unbounded predictions from negative binomial GAMs AroraandCeferino(2023b)showedthatNegativeBinomialcanoverpredictoutages to values even larger than the number of customers since it had no upper bound for Hurricane Isaias (2020) in New Jersey. Similarly, negative binomial GAM predicted more outages than the number of customers in 7 cities in Florida during Hurricane Ian (2022),withanaverageoverpredictionratioof5.7. Forexample,thenegativebinomial regression model predicted 10,538 outages in Boca Grande, Florida, 7.6 times more than the number of customers, which is just 1,385. Conversely, the quasi-binomial model makes bounded predictions between 0 and 1 for the fraction of customers without power. Therefore, we did not observe any outage predictions larger than the number of customers. Lack of extrapolability with random forest We expect more damage to the power systems with increased wind speeds, i.e., a monotonic increase in outages with wind speeds. Arora and Ceferino (2023b) showed thatrandomforestmodelscanfailtocapturethesemonotonicincreases,leadingtopoor extrapolation of outage predictions for high and low winds when trained with limited outage data. To study this issue, we plotted how the random forest and quasi-binomial models capture the relationship between winds and outage predictions (Figure 7). We obtainedthepartialdependenceplotbyvaryingthewindspeedsfrom 0𝑚/𝑠to120𝑚/𝑠 whilekeepingallothervariablesconstantandtakingtheaverageoutagepredictionsover 20 Arora, August 18, 2023 all the observations in the training data. In this case, the wind data for training have a 5th and 95th percentile of 11.4𝑚/𝑠and62.5𝑚/𝑠. For the random forest model, we observe an almost monotonic increase in the fraction of customers without power within this range. However, we observe poor extrapolations outside. For example, we do not see a decrease in outages below winds of 11.4𝑚/𝑠, and predictions for wind speeds of 0𝑚/𝑠is above 0. Also, the maximum outages saturate at 0.85with the randomforestat 62.5𝑚/𝑠. ItisraretoobservehighwindeventsintheUS,aseventhe most wind-affected areas, such as Florida, have a return period of 100 years for winds of60𝑚/𝑠. The sparseness of the power outage data from hurricanes makes it difficult for the non-parametric random forest model to follow the simple mechanics of power infrastructure failures with higher winds. On the contrary, we observed a monotonic increase in risk for outages with wind speedsinFigure7aforthequasi-binomialmodel. Theoutagesconvergetozeroatwind speedsof 0𝑚/𝑠andcontinuetogrowwithincreasedwindspeedsforthequasi-binomial model. Fig. 7.a. Partial dependence of outages with winds for random forest and quasi- binomialmodels,b. Distributionofthewindspeedsinthetrainingdata. Randomforest fails to extrapolate for outages and winds outside the range of training data. Uncertainty Quantification in Outage Predictions AroraandCeferino(2023b)alsodiscussedtheinabilityofnegativebinomialGAMs and random forest models to predict uncertainties in power outage predictions that 21 Arora, August 18, 2023 follow the mechanics of infrastructure failures. Extreme events such as Hurricane Ian (2022),acategory5hurricanewithwindgustsover70m/s,arequitelikelytodevastate power systems. When Hurricane Ian (2022) made landfall, close to 100% customers in Charlotte County lost power. Another similar example is Hurricane Maria (2018) which destroyed power systems in Puerto Rico (Campbell 2018). On the contrary, we will likely observe no damage to power infrastructure at wind speeds closer to zero. The power system is structurally designed to resist at least wind speeds of 20𝑚/𝑠 (Bjarnadottir et al. 2013; IEEE 2007). Thus, we expect outage fractions of zero and close to one at low (close to zero) and high winds ( >70𝑚/𝑠), respectively, with very high certainty. We explore how uncertainty in outage predictions is captured by quasi-binomial, randomforest,andnegativebinomialmodelsbyanalyzingpredictedstandarddeviations for outages at different wind speeds in New Jersey (Figure 8). We vary winds keeping othervariablesconstantforacleancomparison,andshowthe 90%confidenceintervals on predictions on fractions of customers without power in Tuckerton Borough, New Jersey. Fig. 8.Uncertainty quantification in power outages for different models. Solid black represents the mean fraction of customers without power, and blue dotted lines repre- sent the 95% confidence intervals on outage predictions. a. Quasi-binomial regression method: shows lower uncertainty in outage predictions at low and high winds. b. Negative binomial regression method: shows lower uncertainty at low winds but large uncertaintiesathighwinds. c.Randomforestregressionmethod: showshighuncertain- ties at low and high winds. Werescaledthenegativebinomialpredictionstocomparethefractionofcustomers without power for all models. We truncated the predictions on the outage fraction for the negative binomial at 1 as there is no intrinsic upper bound in the model. For the negativebinomialfunction,varianceinoutagesgrowsasafunctionofthesquareofthe meanoutages(AroraandCeferino2023b)(Figure8b). Asdiscussedearlier,weexpect small uncertainties (e.g., small confidence interval ranges) for predictions of outage fractions close to 0 (e.g., at winds close to 0𝑚/𝑠) and 100% (e.g., at high winds of >70𝑚/𝑠). For negative binomial, the 90%confidence interval’s range is close to zero atlowwindspeeds,consistentwiththeexpectedlittletonodamagefrequentlyobserved in power systems at low wind speeds. However, we observed a very high variability 22 Arora, August 18, 2023 withincreasingwindspeeds,e.g.,the 90%confidenceintervalrangesfrom0.12to0.62 forwindspeedsof 80𝑚/𝑠. Thehighvariabilityathighwindsfornegativebinomialdoes not agree with the mechanics of infrastructure failures as complete damage is expected at high winds (>70m/s) with high certainty. Weusedquantilerandom forest(QRF)regressiontoobtainconfidenceintervals on predictions for a fraction of outages. QRF gives confidence intervals purely from the data because the random forest is non-parametric and does not assume any underlying probability distribution. Random forest predicts the outage with a 90%confidence intervalwithawiderange(almostfrom0to1)atwindsclosetozeroandarangeof0.52 athighwinds(from0.48to1). Thus,randomforestfailstoaccountforthephysicsofthe structuralbehaviorofpowersystemsatlowandhighwinds. Itisalsoimportanttonote from Figure 8c risk of outages does not grow monotonically. Outage fraction saturates at0.8fortherandomforestregression,showingthelimitedcapabilitytocapturetherisk correctly at high winds, consistent with observations in Arora and Ceferino (2023b). In contrast, quasi-binomial predictions show a small range in the 90%confidence intervalforoutagepredictionsatlowandhighwinds. Also,thepredictionsconvergeto1 withahighprobabilityatlargewinds. Wealsoobservedamonotonicincreaseintherisk of outages with wind speeds. Since quasi-binomial predictions are bounded between 0 and1,predictedoutagesneverexceed 100%ofcustomerswithoutpower. Therefore,the quasi-binomialmodelovercomesthe threelimitationsofstate-of-artoutagemodelsfor power outage predictions following closer the physics of infrastructure failures under hurricanes. RESILIENCE WITH REDUNDANT POWER INFRASTRUCTURE We investigated whether the quasi-binomial model captures how the power system configuration affects its vulnerability to hurricanes. To cleanly study the effect of redundancy, we explore two neighboring cities in New Jersey with similar population densities but different numbers of substations. We studied Upper Freehold Township, with 17 substations in the county, with 0.27 substations per 10,000 population and a population density of 0.06 /1000𝑚2, and Hamilton Township, with 25 substations in its county, with 0.67 substations per 10,000 population and a population density of 0.83/1000𝑚2. TheseareneighboringcitieswheretheDistHurrparameter(stormsurge effect parameter) will be almost similar, and thus, variation in winds will determine the variation in the extent of damage to power systems. We plotted the fragility for UpperFreeholdTownshipandHamiltonTownshipinFigure10. Weobservethatmore substations improve the resilience of power systems to hurricanes. For example, at a 100-year return wind gust of 32𝑚/𝑠for these cities, the fraction of outages is 0.27 (40.7% less) versus 0.51 for Hamilton Township than for Upper Freehold Township (25 substations versus 17). More substations generally indicate more transformers and feederssupplyingpowertocustomers,whichimpliesthattheHamiltonTownships’grid hasmoreredundancy. Weassumefourfeedersarebeingservedthroughonesubstation (Brown2002). Thus,aboutonefeederandthreefeedersper10,000customersinUpper Freehold Township (0.27 substations per 10,000 customers) and Hamilton Township (0.67 substations per 10,000 customers), respectively. To illustrate this point better, we draw the analogy to an example circuit shown in 23 Arora, August 18, 2023 Figure 9 with feeders connected to poles (yellow circles) to supply power to customers (red circles). In Figure 9a, all the customers receive power from one feeder. Thus, a 100%outage will occur when the feeder is inoperable. Similarly, 100%outage will occur for the circuit shown in Figure9b when all three feeders are inoperable. We obtainedtheprobabilityoffailureofonepolegivenwindspeedfromBjarnadottiretal. (2013). If𝑝is the probability of failure of one pole, and the failure is independent of otherpoles,thentheprobabilityoffailureof 𝑘polesis𝑝𝑘. Wepresentedtheprobability of outages for the example circuit in Figure 10b. We observe that redundancy in a circuit with more feeders reduces the circuit’s vulnerability to outages. For example, a community with one feeder has 24%probability of 100%outages at 32𝑚/𝑠compared to0%chanceswiththreefeeders. NotethatbothTownshipsinFigure10ahavenon-zero outagesat 20𝑚/𝑠winds,asduringahurricane,otherfactorssuchastreescanalsocause poweroutages. Thus,ourpresentedpoweroutagemodelcancaptureredundancyeffects even with limited information on the power systems, e.g., the number of substations. Fig. 9.Example Circuit to study the change in risk of hurricane-induced outages to communities by introducing redundancy in power systems a. power is supplied to customers through one feeder b. power is supplied to customers through three feeders CONCLUSION This paper presented a quasi-binomial regression-based hurricane-induced power outagepredictionmodeltoimprovepre-hurricaneresponseplanningtooutages. Wealso testedthebetaregressionmodelforthepoweroutagemodel,butthebetamodelhad 22% higher cross-validation RMSE than the quasi-binomial model. We sought to capture the underlying mechanics of power infrastructure failures with the presented quasi- binomialregressionmodel. Thepresentedoutagemodeladdressestheshortcomingsof existingpoweroutagemodels: (1)unboundedoutagepredictions,(2)poorextrapolation for higher winds, (3) variance estimates in outage predictions unrepresentative of the mechanics of infrastructure failures, and (4) limited ability to account for storm surge vulnerability to power systems. The presented outage model can predict the outage for largepowerinfrastructuresatthecitylevelwithalowerboundofzerotoanupperbound ofoneforthefractionofcustomerswithoutpower. Thequasi-binomialregressionmodel 24 Arora, August 18, 2023 Fig.10.Fragilitycomparisonwiththedifferentnumberofpowercomponents: a. Upper Freehold Township with 17 substations in its county, and Hamilton Township with 25 substationsinitscounty. b. fragilitycurveforoneandthreenumbersoffeedersforthe example circuits in Figure 9 alsoshowedamonotonicincreaseintheriskofoutageswithgustwindsspeeds. Finally, the quasi-binomial regression model shows small uncertainty for outage predictions at low and high winds. We utilize the historical outage datasets from the states of New Jersey, New York, Florida,Texas,andLouisianaintheUnitedStatestodevelopanoutagemodelapplicable to multiple cities prone to hurricane-induced outages. The outage modeling depends ondifferentinputparameters,includinghurricane,environmental,demographiccondi- tions, and power system information. Since many input parameters are only available after the hurricane, we studied outage predictions with limited parameters to enable outage predictions during hurricanes’ early warning, including only wind speed, the minimum distance of city to hurricane crossing land normalized with radius to maxi- mum winds, precipitation, barren area, wetlands area, soil moisture, standard precip- itation index, percent treed area, substation information, and population density. The mostimportantvariablesfortheoutagepredictionswerewindspeedsandtheminimum distance of city to hurricane crossing land normalized with radius to maximum winds, as these two factors are correlated to damage effects to power systems from hurricane winds and storm surges. We compared the performance of the proposed quasi-binomial regression outage modelwithstate-of-the-artoutagemodels,i.e.,thenegativebinomialandrandomforest models. We compared the models’ performance for Hurricane Ian (2022) in Florida. While the negative binomial model and random forest predicted outages with an error of76.30%and49.72%,respectively,theQuasi-Binomialregressionmodelhadanerror of7.04%,predicting2,886,413outagesforHurricaneIaninFlorida,closetotheactual 3,105,294 outages. Weillustratedhowourmodelcapturestheeffectoftheredundancyofpowersystems onhurricane-inducedoutages. Weinvestigatedthevariation ofriskofoutagesincities 25 Arora, August 18, 2023 withdifferentnumbersofsubstationsasaproxyforredundancy. Weobservedthatcities with more substations have fewer customers without electricity. For example, Mercer County has 25 substations, and Monmouth County has 17. At 32𝑚/𝑠of wind speed, Hamliton Township in Mercer County will have ∼41%lesser percent of customers withoutpowerthanUpperFreeholdTownshipinMonmouthCounty. Thus,powergrids withredundantcomponentssuchasparallelormeshedgrids(Brown2002)canincrease the resilience of power systems to extreme weather events. SUPPLEMENTARY INFORMATION Supplementary Information for the article is provided in a separate file. DATA AVAILABILITY All codes supporting this study’s findings are available from the corresponding author upon reasonable request. Power outage data is obtained from PowerOutage (PowerOutage2022). Thesourceforallotherdataisprovidedinthearticle,anddatais available from authors upon reasonable request. ACKNOWLEDGEMENTS The authors are thankful for the financial support provided by the NYU Tandon School of Engineering Fellowship and NYU CUSP Dissertation Fellowship. This research is also partially supported by Coalition for Disaster Resilient Infrastructure Fellowship (Grant No. 201924669).