Learned Control Variates for Rare Grid-Outage Estimation at Electric Vehicle Charging Stations

# Learned Control Variates for Rare Grid-Outage Estimation at Electric Vehicle Charging Stations
Source: [https://arxiv.org/html/2610.27159](https://arxiv.org/html/2610.27159)
Sulyman IyandaAffiliation:Industrial and Systems Engineering Department King Fahd University of Petroleum and Minerals Dhahran, Saudi Arabia sulymaniyanda@gmail\.comMansur AriefAffiliation:Industrial and Systems Engineering Department King Fahd University of Petroleum and Minerals Dhahran, Saudi Arabia mansur\.arief@kfupm\.edu\.sa

###### Abstract

Public charging of electric vehicles \(EVs\) adds a large and synchronized load to the distribution feeders\. In most networks, it is rare for charging demand, weather, and local grid conditions to overload a feeder at the same time\. When it does happen, the consequences can be severe: an overloaded feeder stresses its transformers and can trigger cascading outages across a city\. Planners who design for high reliability therefore need to estimate very small outage probabilities\. To that end, plain Monte Carlo estimation with an expensive grid simulator can require a prohibitively large number of runs for this task\. To reduce the number of runs, we propose a control variate \(CV\) built from a learned classifier of the outage, trained on cheap observable features such as temperature and load forecasts\. We show that for any binary learner, the variance reduction is an explicit function of its recall and false\-positive rate\. Furthermore, we show that control variate based on a Gaussian proxy of the load loses its effect as outages become rarer, while a learner with non\-vanishing precision keeps removing a fixed fraction of the variance even in the extreme rare\-event regime\. We verify these results in numerical experiments using a feeder\-overload simulator across three orders of magnitude of outage probability and obtain that even weak learner with a single feature can remove a quarter of the variance, while a boosted\-tree learner is able to remove about 60%\. We then apply the estimator to the feeder network of Surabaya, Indonesia, in which the Gaussian load model leaves every feeder three to seven times the one\-day\-in\-ten\-years blackout target\. At the same simulation budget, planning the reserve with the learned CV and a 95% confidence limit reduces the city\-wide blackout days by a third compared with Monte Carlo point estimates\. These results show that learned control variates can help planners design distribution grids that stay reliable against rare outages\.

###### Index Terms:

electric vehicle charging, grid reliability, rare\-event simulation, control variates, variance reduction, machine learning

## IIntroduction

Electric vehicle \(EV\) charging is becoming one of the largest new loads on distribution networks\. Uncoordinated charging raises the peak demand and could push feeders and transformers beyond their ratings on the few days when charging coincides with a high ambient temperature, a demand forecast miss, or a storm\[[1](https://arxiv.org/html/2610.27159#bib.bib4),[2](https://arxiv.org/html/2610.27159#bib.bib5),[3](https://arxiv.org/html/2610.27159#bib.bib6),[4](https://arxiv.org/html/2610.27159#bib.bib31)\]\. Hosting\-capacity studies of distribution networks with EV charging stations show that the binding constraint is the loading on a small number of such extreme days\[[5](https://arxiv.org/html/2610.27159#bib.bib35)\]\. In cities where the grid already operates with little reserve, one such day could result in a local outage, in which public charging stations are among the first loads to be shed or to lose supply\[[6](https://arxiv.org/html/2610.27159#bib.bib30),[7](https://arxiv.org/html/2610.27159#bib.bib32)\]\.

Traditionally, the planning literature on charging stations treats an outage as a service failure, in which vehicles arrive, the station is out of service, and revenue is lost\[[8](https://arxiv.org/html/2610.27159#bib.bib1)\]\. The consequences of a feeder overload, however, go beyond the station\. An overloaded feeder that trips its protection moves its load onto the neighbouring feeders, while transformers that run hot age faster\. A sequence of such local events is the mechanism through which large blackouts begin\[[9](https://arxiv.org/html/2610.27159#bib.bib7),[10](https://arxiv.org/html/2610.27159#bib.bib8)\]\. For example, the August 2003 blackout in North America started from a small number of overloaded lines and cost several billion US dollars\[[11](https://arxiv.org/html/2610.27159#bib.bib3)\]\. Hence, for a candidate charging station, the grid operator needs to know the probability that its charging load, together with the other loads on the feeder, causes an overload on a given day\. Reliability planning works with such probabilities at levels of10−310^\{\-3\}to10−610^\{\-6\}per period\. The classical loss\-of\-load criterion of one day in ten years, for instance, corresponds to2\.7×10−42\.7\\times 10^\{\-4\}per day\[[12](https://arxiv.org/html/2610.27159#bib.bib9)\]\.

Estimating such a small probability is a rare\-event simulation problem\. A feeder model with thermal limits, protection logic and stochastic EV arrivals is expensive to run\. The Monte Carlo \(MC\) estimate usingnnsamples with a target probabilityθ\\thetahas relative standard error\(1−θ\)/\(n​θ\)\\sqrt\{\(1\-\\theta\)/\(n\\theta\)\}, which means a10%10\\%relative error at95%95\\%confidence requires about400/θ400/\\thetaruns, i\.e\., four million runs atθ=10−4\\theta=10^\{\-4\}\. Earlier work on robust charging\-station design\[[8](https://arxiv.org/html/2610.27159#bib.bib1)\]is able to reduce this cost with a control variate \(CV\) built from a Gaussian proxy of the electricity load\. This construction, however, assumes that the modeler already knows a Gaussian variable that is strongly correlated with the outage\. In practice however, it is often impractical to obtain such proxy\. Furthermore, as we show in this paper, such a proxy loses its effect as outages become rarer\.

To that end, we propose in this paper to replace the assumed proxy by a learned one\. Cheap features that are observable before the expensive simulation is run, such as the temperature anomaly, the day\-ahead load forecast and a storm flag, are often available\. A classifier trained on a modest number of simulator runs can then predict, though imperfectly, whether the simulator will report an overload\. We show that such a classifier only needs to be a weak learner\[[13](https://arxiv.org/html/2610.27159#bib.bib24)\], whose value as a control variate is determined by its recall and precision only\. We also show that, with a suitable training design, the learner keeps its effect at every rarity level\.

In sum, our main contributions are as follows\.First, we provide an exact variance formula for the CV estimator built from any binary learner, together with the rare\-event limit1−recall×precision1\-\\mathrm\{recall\}\\times\\mathrm\{precision\}of its variance ratio \(Theorem[1](https://arxiv.org/html/2610.27159#Thmtheorem1), Corollary[1](https://arxiv.org/html/2610.27159#Thmcorollary1)\)\.Second, we derive a ceiling on the variance reduction of any learner, which is set by the features, and a proof that a Gaussian proxy loses its effect as the outage becomes rarer \(Theorem[2](https://arxiv.org/html/2610.27159#Thmtheorem2), Proposition[1](https://arxiv.org/html/2610.27159#Thmproposition1)\)\.Third, we show the exact variance of the estimator with an estimated coefficient and an unlabelled pool, the cost\-optimal split between simulator runs and feature draws \(Theorem[3](https://arxiv.org/html/2610.27159#Thmtheorem3)\), and a training design that keeps the learner informative in the tail \(Proposition[2](https://arxiv.org/html/2610.27159#Thmproposition2)\)\.Finally, we demonstrate using a simulation study that compares plain MC, the Gaussian proxy, a weak learner, a strong learner and the Bayes\-posterior ceiling at outage probabilities from10−310^\{\-3\}to10−610^\{\-6\}based on a case study on the Surabaya network\[[8](https://arxiv.org/html/2610.27159#bib.bib1)\], which quantifies how the load model and the precision of the estimate change the reserve on each feeder and the blackout days the city would experience\.

## IIRelated Work

### II\-AEV Charging Load and Grid Outages

The effect of EV charging on distribution networks has been well documented\. Uncoordinated residential charging raises the evening peak and the loading of distribution transformers\[[1](https://arxiv.org/html/2610.27159#bib.bib4),[2](https://arxiv.org/html/2610.27159#bib.bib5)\]\. Furthermore, the choice of where and when vehicles charge could change the grid impact of a given fleet by a large factor\[[3](https://arxiv.org/html/2610.27159#bib.bib6)\]\. Reviews of the distribution\-level impacts list voltage drops, transformer ageing, losses and protection trips as the main consequences of EV charging, with demand management, storage and network reinforcement as the main remedies\[[4](https://arxiv.org/html/2610.27159#bib.bib31),[7](https://arxiv.org/html/2610.27159#bib.bib32)\]\. Hosting\-capacity studies, including recent work which shows that the capacity of a feeder for PV and charging stations is set by the loading on a small number of extreme days\[[5](https://arxiv.org/html/2610.27159#bib.bib35)\]\. Studies in African cities also show that a grid with little reserve shifts from occasional to frequent interruptions as charging grows\[[6](https://arxiv.org/html/2610.27159#bib.bib30)\]\. The consequences of a feeder overload, however, go beyond the feeder itself\. The analysis of cascading failures and of the record of large blackouts in North America shows that a sequence of local overloads is the mechanism through which large blackouts begin\[[9](https://arxiv.org/html/2610.27159#bib.bib7),[10](https://arxiv.org/html/2610.27159#bib.bib8)\], often at a large economic cost\[[11](https://arxiv.org/html/2610.27159#bib.bib3),[14](https://arxiv.org/html/2610.27159#bib.bib2)\]\. For this reason, reliability planning works with small outage probabilities per period\[[12](https://arxiv.org/html/2610.27159#bib.bib9)\]\. These studies describe the mechanism and the stakes of an overload\. However, they do not address how a utility could estimate a probability of order10−310^\{\-3\}to10−610^\{\-6\}for a specific feeder with a simulator that it can only run a limited number of times\.

### II\-BCharging\-Station Planning under Uncertainty

Charging\-infrastructure planning is usually formulated as a facility\-location problem\[[15](https://arxiv.org/html/2610.27159#bib.bib45)\]that trades investment and operating cost against service quality\. Early models were deterministic\[[16](https://arxiv.org/html/2610.27159#bib.bib41),[17](https://arxiv.org/html/2610.27159#bib.bib38)\], while later work added stochastic demand\[[18](https://arxiv.org/html/2610.27159#bib.bib42),[19](https://arxiv.org/html/2610.27159#bib.bib43)\]\. A smaller number of studies considered disruptions such as floods\[[20](https://arxiv.org/html/2610.27159#bib.bib39)\]or outages mitigated by on\-site storage\[[21](https://arxiv.org/html/2610.27159#bib.bib40)\]\. The review in\[[22](https://arxiv.org/html/2610.27159#bib.bib34)\]calls for models that quantify randomness explicitly\. In the facility\-location literature on unreliable facilities\[[23](https://arxiv.org/html/2610.27159#bib.bib44)\], the failure probabilities are treated as known inputs\. Earlier work on robust charging\-station design\[[8](https://arxiv.org/html/2610.27159#bib.bib1)\]is, to our knowledge, the first to estimate station reliabilities by simulation inside the planning loop, using a Gaussian\-proxy control variate to make this affordable\. However, such a proxy must be known in advance\. Furthermore, the Gaussian load model behind it is exactly the model that, as we show in this paper, misjudges the tail\. In this paper, we leave the location model aside and focus on the estimation step, which all of these models take as given\.

### II\-CRare\-Event Simulation and Variance Reduction

Estimating a small probability by plain MC is inefficient, since the relative error grows as1/n​θ1/\\sqrt\{n\\theta\}whenθ\\thetadecreases\[[24](https://arxiv.org/html/2610.27159#bib.bib14),[25](https://arxiv.org/html/2610.27159#bib.bib12)\]\. Two classical remedies are importance sampling and control variates\. Importance sampling changes the sampling distribution, whereas control variates keep it and subtract a correlated correction with zero mean\[[26](https://arxiv.org/html/2610.27159#bib.bib10),[27](https://arxiv.org/html/2610.27159#bib.bib16),[28](https://arxiv.org/html/2610.27159#bib.bib17),[29](https://arxiv.org/html/2610.27159#bib.bib15)\]\. Control variates for probability estimation in particular were studied in\[[30](https://arxiv.org/html/2610.27159#bib.bib18)\]\. Learned classifiers have also been used to build certifiable importance\-sampling proposals for black\-box rare events\[[31](https://arxiv.org/html/2610.27159#bib.bib26)\]\. With a learned proposal, however, the quality of the estimate depends on how well the proposal covers the rare\-event set, which has to be certified\. A learned control variate, in contrast, keeps the estimator unbiased for any classifier\. Several related ideas use a cheap model to correct an expensive one\. Multi\-fidelity MC uses a low\-fidelity model as a control for a high\-fidelity one and allocates the budget between the fidelities\[[32](https://arxiv.org/html/2610.27159#bib.bib20),[33](https://arxiv.org/html/2610.27159#bib.bib21)\]\. Prediction\-powered inference combines a machine\-learning predictor on a large unlabelled set with a small labelled set\[[34](https://arxiv.org/html/2610.27159#bib.bib22),[35](https://arxiv.org/html/2610.27159#bib.bib23)\]\. Control functionals build controls with a known mean from the score of the sampling density\[[36](https://arxiv.org/html/2610.27159#bib.bib19),[37](https://arxiv.org/html/2610.27159#bib.bib33)\]\. The estimator in Theorem[3](https://arxiv.org/html/2610.27159#Thmtheorem3)\(b\) is algebraically the same as the two\-fidelity estimator and the prediction\-powered estimator\. None of these works, however, considers the case in which the quantity of interest is a rare binary event\. In particular, they do not show how the variance reduction depends on the classifier, whether it survives as the event becomes rarer, or how the classifier should be trained when a nominal training set contains no event at all\.

### II\-DOutage Prediction and Weak Learners

Data\-driven outage prediction models use weather, load and infrastructure features to forecast where and how many outages a storm will cause\[[38](https://arxiv.org/html/2610.27159#bib.bib28),[39](https://arxiv.org/html/2610.27159#bib.bib29)\]\. Probabilistic forecasts of EV charging load and of ultra\-short\-term feeder load are also becoming common tools for operational reliability in smart grids\[[40](https://arxiv.org/html/2610.27159#bib.bib36),[41](https://arxiv.org/html/2610.27159#bib.bib37)\]\. These models are built and tuned for prediction, with the accuracy, the error or the area under the ROC curve as the criterion\. In machine learning, a weak learner is a classifier whose balanced error is only slightly below one half\. Boosting shows that such learners can be combined into a strong one\[[13](https://arxiv.org/html/2610.27159#bib.bib24),[42](https://arxiv.org/html/2610.27159#bib.bib25)\]\. In this paper, we use the notion for a different purpose\. We show that a weak learner is already a useful control variate, whose value is set by its recall and precision on the rare days \(Section[IV](https://arxiv.org/html/2610.27159#S4)\)\. These criteria differ from those used for prediction\.

## IIIProblem Setting and Estimators

### III\-AFeeders, Charging Stations and Outages

We consider a distribution feeder that supplies, among other loads, a public EV charging station\. The feeder has a rating, i\.e\., the thermal limit of its conductors or of the transformer behind it, and a protection scheme that disconnects it when the load exceeds the rating for long enough\. On an ordinary day, the peak loading stays well below the rating\. However, in a few days, several effects coincide: a high ambient temperature raises the air\-conditioning load and reduces the capacity of the cables, the day\-ahead forecast underestimates the demand, a storm shifts load onto the feeder, and many vehicles arrive at the station at once\. The loading then reaches the rating and trips\. We call such a day an*outage day*\.

An outage day has two kinds of cost\. For the charging operator, the station is out of service and revenue is lost\. For the grid operator, the tripped feeder moves its load onto its neighbors, whose transformers then run hot and age faster\. If the neighboring feeders are themselves close to their ratings, the disturbance could propagate as a cascading failure\[[9](https://arxiv.org/html/2610.27159#bib.bib7)\]\. This second cost is borne by every customer on the affected feeders and is much larger than the first\[[14](https://arxiv.org/html/2610.27159#bib.bib2)\]\. Hence, the approval of a charging station depends on the probability of an outage day on its feeder, which the grid operator requires to be small\. For example, the criterion of one loss\-of\-load day in ten years\[[12](https://arxiv.org/html/2610.27159#bib.bib9)\]corresponds to2\.7×10−42\.7\\times 10^\{\-4\}per day\. In general, the probabilities that a planner needs to certify may lie between10−310^\{\-3\}and10−610^\{\-6\}\.

The utility has three sources of information about a feeder\.First, it knows the current rating and has an outage record, from which the current outage probability is known with reasonable accuracy\. In the Surabaya data of\[[8](https://arxiv.org/html/2610.27159#bib.bib1)\], for instance, the reported reliabilities of the eleven feeders lie between0\.9740\.974and0\.9830\.983\.Second, it has cheap observable features for every past or simulated day, such as the temperature anomaly, the day\-ahead load\-forecast anomaly, a storm flag, the day of the week and the season\.Third, it has a simulator of the feeder that takes the features of a day, draws its own internal randomness \(vehicle arrivals, equipment state and protection behaviour\), computes the loading and reports whether the feeder trips\. One simulated day is expensive compared with one feature draw, since it involves a power\-flow or thermal model and a stochastic arrival process\. The planning question is which reserve, through a feeder upgrade or a cap on the charging load, brings the outage probability below the target\. This work focuses on the estimation part of this question: how to estimate the outage probability of a feeder at a given reserve with as few simulator runs as possible, together with a reliable and meaningful confidence interval for such a critical task\.

### III\-BNotation and Model

LetY=𝟏​\{outage day\}Y=\\mathbf\{1\}\\\{\\text\{outage day\}\\\}, and letθ=ℙ⁡\(Y=1\)\\theta=\\mathbb\{P\}\(Y=1\)be the outage probability of the feeder at the reserve under study\. The features of a day form a vectorX∈𝒳X\\in\\mathcal\{X\}with distributionFXF\_\{X\}, which can be drawn at costcXc\_\{X\}\. The simulator returnsYYgivenXXat costcY≫cXc\_\{Y\}\\gg c\_\{X\}\. We writeη⁡\(x\)=ℙ⁡\(Y=1∣X=x\)\\eta\(x\)=\\mathbb\{P\}\(Y=1\\mid X=x\)for the conditional outage probability given the features\. We assume a simulator, which can return a continuous severitySSand the standardized peak loading of the day, henceY=𝟏\{S≥κ\}Y=\\mathbf\{1\}\\\{S\\geq\\kappa\\\}for a ratingκ\\kappa\. A higher reserve thus corresponds to a higherκ\\kappa\. In the analysis that follows, we treats the simulator as a black box that simply returnsYY\.

Planners often associate the outage probabilityθ\\thetawith the following: the expected number of outage days per year \(i\.e\.,365​θ365\\,\\theta\), the probability of at least one outage day in a horizon ofTTdays \(i\.e\.,1−\(1−θ\)T≈T​θ1\-\(1\-\\theta\)^\{T\}\\approx T\\theta\), and and the station reliability \(i\.e\.,1−θ1\-\\theta\)\. All three are linear inθ\\thetato first order, hence a variance reduction obtained forθ\\thetaapply to them directly\.

### III\-CMonte Carlo and Control Variates

Withnnindependent simulated days,θ^MC=Y¯n\\hat\{\\theta\}^\{\\mathrm\{MC\}\}=\\bar\{Y\}\_\{n\}is unbiased withVar⁡\(θ^MC\)=θ⁡\(1−θ\)/n\\mathrm\{Var\}\(\\hat\{\\theta\}^\{\\mathrm\{MC\}\}\)=\\theta\(1\-\\theta\)/n\. The number of runs required for a relative half\-widthε\\varepsilonat confidence1−α1\-\\alphais

nMC​\(θ,ε,α\)=z1−α/22​\(1−θ\)/\(θ​ε2\),n\_\{\\mathrm\{MC\}\}\(\\theta,\\varepsilon,\\alpha\)=z\_\{1\-\\alpha/2\}^\{2\}\(1\-\\theta\)/\(\\theta\\,\\varepsilon^\{2\}\),\(1\)which grows as1/θ1/\\theta\. For a10%10\\%relative error at95%95\\%confidence, this is about400/θ400/\\thetaruns, i\.e\.,4×1054\\times 10^\{5\}runs atθ=10−3\\theta=10^\{\-3\}and4×1084\\times 10^\{8\}at10−610^\{\-6\}\. In other words, the planning has to observe about400400outage days in the Monte Carlo simulation to achieve this level of relative error\. To that end, a control variate is a random variableCCthat is simulated together withYY, is correlated withYY, and has a known meanγ=𝔼⁡\[C\]\\gamma=\\mathbb\{E\}\[C\]\. For any coefficientβ\\beta, the estimator

θ^βCV=Y¯n−β⁡\(C¯n−γ\)\\hat\{\\theta\}^\{\\mathrm\{CV\}\}\_\{\\beta\}=\\bar\{Y\}\_\{n\}\-\\beta\\,\(\\bar\{C\}\_\{n\}\-\\gamma\)\(2\)is unbiased\. Furthermore, withρ=Corr⁡\(Y,C\)\\rho=\\mathrm\{Corr\}\(Y,C\)andβ⋆=Cov⁡\(Y,C\)/Var⁡\(C\)\\beta^\{\\star\}=\\mathrm\{Cov\}\(Y,C\)/\\mathrm\{Var\}\(C\), its variance is\[[27](https://arxiv.org/html/2610.27159#bib.bib16),[29](https://arxiv.org/html/2610.27159#bib.bib15)\]

Var⁡\(θ^β⋆CV\)=\(1−ρ2\)​Var​\(Y\)/n\.\\mathrm\{Var\}\(\\hat\{\\theta\}^\{\\mathrm\{CV\}\}\_\{\\beta^\{\\star\}\}\)=\(1\-\\rho^\{2\}\)\\,\\mathrm\{Var\}\(Y\)/n\.\(3\)Thus,ρ2\\rho^\{2\}is the fraction of variance removed,1/\(1−ρ2\)1/\(1\-\\rho^\{2\}\)is the factor by which the number of simulator runs can be reduced, and\(1−ρ2\)​\(1−θ\)/\(n​θ\)\\sqrt\{\(1\-\\rho^\{2\}\)\(1\-\\theta\)/\(n\\theta\)\}is the relative standard error\.

In an earlier study\[[8](https://arxiv.org/html/2610.27159#bib.bib1)\], researcher use the control variateC=𝟏\{P≥t\}C=\\mathbf\{1\}\\\{P\\geq t\\\}for a Gaussian proxyPPof the load, withγ=Φ¯​\(t\)=1−Φ​\(t\)\\gamma=\\bar\{\\Phi\}\(t\)=1\-\\Phi\(t\)\. In this paper, we useC=g⁡\(X\)C=g\(X\)for a learnerggtrained on a separate set of simulator runs\. Its meanγ\\gammais either known, whenggthresholds a feature with a known distribution, or estimated from a pool of unlabelled feature draws, such as the features of every day in the utility’s weather archive\. Algorithm[1](https://arxiv.org/html/2610.27159#alg1)summarizes the procedure\.

Algorithm 1Learned control variate \(LCV\)Input:feature law

FXF\_\{X\}\(cost

cXc\_\{X\}per draw\); simulator

Y\|XY\\mid X\(cost

cYc\_\{Y\}\); training design

qq; sizes

ntrn\_\{\\mathrm\{tr\}\},

nn,

mm; learner class

𝒢\\mathcal\{G\}
1Draw

ntrn\_\{\\mathrm\{tr\}\}training features from

qq, run the simulator on them, and fit

g∈𝒢g\\in\\mathcal\{G\}with importance weights

FX/qF\_\{X\}/qwhere the fitting criterion requires them \(Section[IV\-D](https://arxiv.org/html/2610.27159#S4.SS4)\);

2Draw

nnfresh features from

FXF\_\{X\}, run the simulator, obtain

Y1,…,YnY\_\{1\},\\dots,Y\_\{n\}; set

Cl=g⁡\(Xl\)C\_\{l\}=g\(X\_\{l\}\);

3if*γ=𝔼⁡\[g⁡\(X\)\]\\gamma=\\mathbb\{E\}\[g\(X\)\]is known*then

4use

γ\\gamma;

5else

6draw

mmfurther features \(no simulator runs\) and set

γ^=1n\+m​\(∑l≤ng⁡\(Xl\)\+∑l≤mg⁡\(Xl′\)\)\\hat\{\\gamma\}=\\frac\{1\}\{n\+m\}\\big\(\\sum\_\{l\\leq n\}g\(X\_\{l\}\)\+\\sum\_\{l\\leq m\}g\(X^\{\\prime\}\_\{l\}\)\\big\);

7end if

8

β^←Cov^​\(Y,C\)/Var^​\(C\)\\hat\{\\beta\}\\leftarrow\\widehat\{\\mathrm\{Cov\}\}\(Y,C\)/\\widehat\{\\mathrm\{Var\}\}\(C\)from the

nnpairs;

θ^←Y¯n−β^​\(C¯n−γ\)\\hat\{\\theta\}\\leftarrow\\bar\{Y\}\_\{n\}\-\\hat\{\\beta\}\(\\bar\{C\}\_\{n\}\-\\gamma\)\(or

γ^\\hat\{\\gamma\}\);

Output:

θ^\\hat\{\\theta\}and its standard error from \([3](https://arxiv.org/html/2610.27159#S3.E3)\) or \([7](https://arxiv.org/html/2610.27159#S4.E7)\)\.

## IVAnalytical Results

Throughout this section,Y∈\{0,1\}Y\\in\\\{0,1\\\}withθ=ℙ⁡\(Y=1\)∈\(0,1\)\\theta=\\mathbb\{P\}\(Y=1\)\\in\(0,1\), andX∈𝒳X\\in\\mathcal\{X\}is the feature vector\.

### IV\-ABinary Learners

###### Definition 1\.

A binary learner is a measurableg:𝒳→\{0,1\}g:\\mathcal\{X\}\\to\\\{0,1\\\}with recallτ=ℙ⁡\(g⁡\(X\)=1∣Y=1\)\\tau=\\mathbb\{P\}\(g\(X\)=1\\mid Y=1\), false\-positive rateϕ=ℙ⁡\(g⁡\(X\)=1∣Y=0\)\\phi=\\mathbb\{P\}\(g\(X\)=1\\mid Y=0\), alarm rateγ=ℙ⁡\(g⁡\(X\)=1\)=θ​τ\+\(1−θ\)​ϕ\\gamma=\\mathbb\{P\}\(g\(X\)=1\)=\\theta\\tau\+\(1\-\\theta\)\\phiand precisionπ=ℙ⁡\(Y=1∣g⁡\(X\)=1\)=θ​τ/γ\\pi=\\mathbb\{P\}\(Y=1\\mid g\(X\)=1\)=\\theta\\tau/\\gamma\. Following\[[42](https://arxiv.org/html/2610.27159#bib.bib25),[13](https://arxiv.org/html/2610.27159#bib.bib24)\],ggis a*weak learner*with edgeε\>0\\varepsilon\>0if its balanced error rate12​\(\(1−τ\)\+ϕ\)\\tfrac\{1\}\{2\}\(\(1\-\\tau\)\+\\phi\)is at most12−ε\\tfrac\{1\}\{2\}\-\\varepsilon, i\.e\.,τ−ϕ≥2​ε\\tau\-\\phi\\geq 2\\varepsilon\. Note that no requirement is placed on accuracy, calibration or the area under the ROC curve\.

With this definition, we have the following:

###### Theorem 1\(Variance reduction from a binary learner\)\.

For a binary learner with0<γ<10<\\gamma<1andC=g⁡\(X\)C=g\(X\),

Cov⁡\(Y,C\)=θ⁡\(1−θ\)​\(τ−ϕ\),ρ2=\(τ−ϕ\)2​θ​\(1−θ\)γ⁡\(1−γ\)\.\\mathrm\{Cov\}\(Y,C\)=\\theta\(1\-\\theta\)\(\\tau\-\\phi\),\\quad\\rho^\{2\}=\\frac\{\(\\tau\-\\phi\)^\{2\}\\,\\theta\(1\-\\theta\)\}\{\\gamma\(1\-\\gamma\)\}\.\(4\)Hence, the CV estimator has a smaller variance than MC if and only ifτ≠ϕ\\tau\\neq\\phi\. A weak learner with edgeε\\varepsilonguaranteesρ2≥4​ε2​θ​\(1−θ\)/\(γ⁡\(1−γ\)\)\>0\\rho^\{2\}\\geq 4\\varepsilon^\{2\}\\theta\(1\-\\theta\)/\(\\gamma\(1\-\\gamma\)\)\>0\.

###### Proof\.

𝔼⁡\[Y​C\]=ℙ⁡\(Y=1,g⁡\(X\)=1\)=θ​τ\\mathbb\{E\}\[YC\]=\\mathbb\{P\}\(Y=1,g\(X\)=1\)=\\theta\\tauand𝔼⁡\[Y\]​𝔼​\[C\]=θ​γ\\mathbb\{E\}\[Y\]\\mathbb\{E\}\[C\]=\\theta\\gamma, soCov⁡\(Y,C\)=θ⁡\(τ−γ\)=θ⁡\(1−θ\)​\(τ−ϕ\)\\mathrm\{Cov\}\(Y,C\)=\\theta\(\\tau\-\\gamma\)=\\theta\(1\-\\theta\)\(\\tau\-\\phi\)\. Dividing its square byVar⁡\(Y\)​Var​\(C\)=θ⁡\(1−θ\)​γ​\(1−γ\)\\mathrm\{Var\}\(Y\)\\mathrm\{Var\}\(C\)=\\theta\(1\-\\theta\)\\gamma\(1\-\\gamma\)givesρ2\\rho^\{2\}\. ∎

Sinceρ2\\rho^\{2\}depends on\(τ−ϕ\)2\(\\tau\-\\phi\)^\{2\}, a learner that is worse than chance is as useful as one that is better by the same margin, and only the sign ofβ⋆\\beta^\{\\star\}changes\. This means, the only unusable learner is one whose alarms are independent of the outage\. More interestingly, we have the following limit:

###### Corollary 1\(Rare\-event limit\)\.

Forτ\>0\\tau\>0,ρ2=τ​π​\(1−ϕ/τ\)2​\(1−θ\)/\(1−γ\)\\rho^\{2\}=\\tau\\pi\\,\(1\-\\phi/\\tau\)^\{2\}\(1\-\\theta\)/\(1\-\\gamma\)\. Ifθ→0\\theta\\to 0while the precision stays bounded away from zero, thenγ,ϕ→0\\gamma,\\phi\\to 0and

Var⁡\(θ^β⋆CV\)Var⁡\(θ^MC\)→1−recall×precision\.\\frac\{\\mathrm\{Var\}\(\\hat\{\\theta\}^\{\\mathrm\{CV\}\}\_\{\\beta^\{\\star\}\}\)\}\{\\mathrm\{Var\}\(\\hat\{\\theta\}^\{\\mathrm\{MC\}\}\)\}\\;\\to\\;1\-\\mathrm\{recall\}\\times\\mathrm\{precision\}\.\(5\)If insteadϕ\\phiis held fixed asθ→0\\theta\\to 0, thenρ2=O⁡\(θ\)→0\\rho^\{2\}=O\(\\theta\)\\to 0\.

###### Proof\.

Substituteθ/γ=π/τ\\theta/\\gamma=\\pi/\\tauinto \([4](https://arxiv.org/html/2610.27159#S4.E4)\)\. Every factor other thanτ​π\\tau\\pitends to one whenγ,ϕ→0\\gamma,\\phi\\to 0\. The fixed\-ϕ\\phicase follows from \([4](https://arxiv.org/html/2610.27159#S4.E4)\) withγ→ϕ\\gamma\\to\\phi\. ∎

Corollary[1](https://arxiv.org/html/2610.27159#Thmcorollary1)gives the operational meaning of a weak learner for rare events\. Any weak learner reduces the variance\. However, the reduction stays bounded away from zero as the outage becomes rarer only if the false\-positive rate shrinks in proportion toθ\\theta, i\.e\., only if the learner has non\-vanishing precision\. High recall is easy to obtain by predicting an outage every day, whereas precision has to be learned\. Hence, a forecaster tuned to accuracy or ROC criteria is not tuned to the quantity that is meaningful for estimation\.

### IV\-BGeneral Learners and the Ceiling Set by the Features

With regard to more general learners and features, we have the following:

###### Theorem 2\(Factorisation and ceiling\)\.

Letη⁡\(x\)=ℙ⁡\(Y=1∣X=x\)\\eta\(x\)=\\mathbb\{P\}\(Y=1\\mid X=x\), and letg:𝒳→ℝg:\\mathcal\{X\}\\to\\mathbb\{R\}have0<Var⁡\(g⁡\(X\)\)<∞0<\\mathrm\{Var\}\(g\(X\)\)<\\infty\. Then

ρ2​\(Y,g⁡\(X\)\)=ρ2​\(η⁡\(X\),g⁡\(X\)\)⋅ρmax2\\rho^\{2\}\(Y,g\(X\)\)=\\rho^\{2\}\(\\eta\(X\),g\(X\)\)\\cdot\\rho^\{2\}\_\{\\max\}\(6\)whereρmax2=Var⁡\(η⁡\(X\)\)/θ⁡\(1−θ\)\\rho^\{2\}\_\{\\max\}=\\mathrm\{Var\}\(\\eta\(X\)\)/\\theta\(1\-\\theta\)\. Hence,ρ2​\(Y,g⁡\(X\)\)≤ρmax2\\rho^\{2\}\(Y,g\(X\)\)\\leq\\rho^\{2\}\_\{\\max\}for every learner, with equality if and only ifg⁡\(X\)g\(X\)is almost surely an affine function ofη⁡\(X\)\\eta\(X\)\.

###### Proof\.

By the tower property, we haveCov⁡\(Y,g⁡\(X\)\)=Cov⁡\(η⁡\(X\),g⁡\(X\)\)\\mathrm\{Cov\}\(Y,g\(X\)\)=\\mathrm\{Cov\}\(\\eta\(X\),g\(X\)\)\. Therefore,

ρ2​\(Y,g\)=\[Cov​\(η,g\)2/\(Var⁡\(η\)​Var​\(g\)\)\]⋅\[Var⁡\(η\)/Var⁡\(Y\)\],\\rho^\{2\}\(Y,g\)=\[\\mathrm\{Cov\}\(\\eta,g\)^\{2\}/\(\\mathrm\{Var\}\(\\eta\)\\mathrm\{Var\}\(g\)\)\]\\cdot\[\\mathrm\{Var\}\(\\eta\)/\\mathrm\{Var\}\(Y\)\],and Cauchy–Schwarz inequality bounds first factor by 1\. ∎

Three consequences follow from this theorem\. First, calibration does not matter, since an affine rescaling of a score leavesρ2\\rho^\{2\}unchanged\. Second, a thresholded score is a binary learner and is governed by Theorem[1](https://arxiv.org/html/2610.27159#Thmtheorem1)\. Third, the ceiling is a property of the features\. The information that a utility collects about each day therefore bounds what any learner, weak or strong, could give in terms of the estimation efficiency\.

For Gaussian proxy, we have the following result:

###### Proposition 1\(A Gaussian proxy loses its effect for rare events\)\.

Let\(S,P\)\(S,P\)be bivariate normal with standard margins and correlationr∈\(−1,1\)r\\in\(\-1,1\),Y=𝟏\{S≥κ\}Y=\\mathbf\{1\}\\\{S\\geq\\kappa\\\}andCt=𝟏\{P≥t\}C\_\{t\}=\\mathbf\{1\}\\\{P\\geq t\\\}\. Thensuptρ2​\(Y,Ct\)→0\\sup\_\{t\}\\rho^\{2\}\(Y,C\_\{t\}\)\\to 0asκ→∞\\kappa\\to\\infty\.

###### Proof sketch\.

For a binary control with nonnegative covariance, \([4](https://arxiv.org/html/2610.27159#S4.E4)\) givesρ2≤τ​π/\(\(1−θ\)​\(1−γ\)\)\\rho^\{2\}\\leq\\tau\\pi/\(\(1\-\\theta\)\(1\-\\gamma\)\)\. Taker≥0r\\geq 0, since the caser<0r<0follows by replacingPPwith−P\-P, and letc=\(1\+r\)/2<1c=\(1\+r\)/2<1\. Fort≥c​κt\\geq c\\kappa, the recallℙ⁡\(P≥t∣S≥κ\)\\mathbb\{P\}\(P\\geq t\\mid S\\geq\\kappa\)vanishes asκ→∞\\kappa\\to\\infty\. Indeed, conditionally onS=s≈κS=s\\approx\\kappa, the proxy isN⁡\(r​s,1−r2\)N\(rs,1\-r^\{2\}\), andt−r​κ≥κ⁡\(1−r\)/2t\-r\\kappa\\geq\\kappa\(1\-r\)/2\. For0≤t<c​κ0\\leq t<c\\kappa, the precision is at mostΦ¯​\(κ\)/Φ¯​\(c​κ\)\\bar\{\\Phi\}\(\\kappa\)/\\bar\{\\Phi\}\(c\\kappa\), which tends to zero by the Mills ratio\. Fort<0t<0, we use the complementary control1−Ct1\-C\_\{t\}, whose recall vanishes in the same way\. The three bounds are uniform intt\. ∎

The mechanism behind Proposition[1](https://arxiv.org/html/2610.27159#Thmproposition1)is the asymptotic tail independence of the bivariate normal distribution\[[43](https://arxiv.org/html/2610.27159#bib.bib27)\]\. The proxy and the load agree on moderate days but disagree on the extreme days that determine the outage\. By Corollary[1](https://arxiv.org/html/2610.27159#Thmcorollary1), a learner whose precision does not vanish keeps a constant fraction of the variance\. Whether such a learner exists depends on the mechanism of the outage\. If extreme loads are driven by observable extremes \(heat waves, forecast misses, storms\), the features stay informative in the tail\. If they are driven by unobserved noise, no estimator using those features is helpful\.

### IV\-CEstimated Coefficient, Unlabelled Pool, and Budget

To make our results practical, we extend our analysis to the estimated coefficient regime and report the following:

###### Theorem 3\(Estimated coefficient and unlabelled pool\)\.

Let thennlabelled pairs and themmunlabelled feature draws be independent, withggfixed,Var⁡\(C\)\>0\\mathrm\{Var\}\(C\)\>0andρ=Corr⁡\(Y,g⁡\(X\)\)\\rho=\\mathrm\{Corr\}\(Y,g\(X\)\)\.

1. \(a\)Withβ^\\hat\{\\beta\}estimated from the samennpairs andγ\\gammaknown,𝔼⁡\[θ^\]=θ\+O⁡\(1/n\)\\mathbb\{E\}\[\\hat\{\\theta\}\]=\\theta\+O\(1/n\)andn​\(θ^−θ\)→N⁡\(0,\(1−ρ2\)​Var​\(Y\)\)\\sqrt\{n\}\(\\hat\{\\theta\}\-\\theta\)\\to N\(0,\(1\-\\rho^\{2\}\)\\mathrm\{Var\}\(Y\)\)\.
2. \(b\)Withγ^\\hat\{\\gamma\}as in Algorithm[1](https://arxiv.org/html/2610.27159#alg1),θ^β\\hat\{\\theta\}\_\{\\beta\}is unbiased for every fixedβ\\beta, and Var⁡\(θ^β⋆\)=Var⁡\(Y\)n​\(1−ρ2​mn\+m\)\.\\mathrm\{Var\}\(\\hat\{\\theta\}\_\{\\beta^\{\\star\}\}\)=\\frac\{\\mathrm\{Var\}\(Y\)\}\{n\}\\Big\(1\-\\rho^\{2\}\\frac\{m\}\{n\+m\}\\Big\)\.\(7\)
3. \(c\)With budgetB=n⁡\(cY\+cX\)\+m​cXB=n\(c\_\{Y\}\+c\_\{X\}\)\+mc\_\{X\}, \([7](https://arxiv.org/html/2610.27159#S4.E7)\) is minimised at\(n\+m\)/n=cY​ρ2/\(cX​\(1−ρ2\)\)\(n\+m\)/n=\\sqrt\{c\_\{Y\}\\rho^\{2\}/\(c\_\{X\}\(1\-\\rho^\{2\}\)\)\}whenever this value exceeds 1\. The variance ratio to MC at the same budget isVar⋆/VarMC=\(\(1−ρ2\)​cY\+ρ2​cX\)2/\(cY\+cX\)\\mathrm\{Var\}^\{\\star\}/\\mathrm\{Var\}^\{\\mathrm\{MC\}\}=\\big\(\\sqrt\{\(1\-\\rho^\{2\}\)c\_\{Y\}\}\+\\sqrt\{\\rho^\{2\}c\_\{X\}\}\\big\)^\{2\}/\(c\_\{Y\}\+c\_\{X\}\), which tends to1−ρ21\-\\rho^\{2\}ascX/cY→0c\_\{X\}/c\_\{Y\}\\to 0\.

###### Proof sketch\.

\(a\) is the classical result for control variates with estimated coefficients\[[27](https://arxiv.org/html/2610.27159#bib.bib16),[28](https://arxiv.org/html/2610.27159#bib.bib17)\]\. \(b\) WriteC¯n−γ^=m/\(n\+m\)⋅\(C¯n−C¯m′\)\\bar\{C\}\_\{n\}\-\\hat\{\\gamma\}=m/\(n\+m\)\\cdot\(\\bar\{C\}\_\{n\}\-\\bar\{C\}^\{\\prime\}\_\{m\}\), where the pool meanC¯m′\\bar\{C\}^\{\\prime\}\_\{m\}is independent of the pairs\. Then expand the variance and substituteβ⋆\\beta^\{\\star\}\. \(c\) Withs=\(n\+m\)/ns=\(n\+m\)/n, the variance equalsVar⁡\(Y\)​f​\(s\)/B\\mathrm\{Var\}\(Y\)f\(s\)/Bwithf⁡\(s\)=\(cY\+s​cX\)​\(1−ρ2\+ρ2/s\)f\(s\)=\(c\_\{Y\}\+sc\_\{X\}\)\(1\-\\rho^\{2\}\+\\rho^\{2\}/s\)\. This function is convex and is minimised at the stated value\. ∎

Part \(b\) coincides with the two\-fidelity estimator\[[32](https://arxiv.org/html/2610.27159#bib.bib20)\]and with prediction\-powered inference with a tuned power parameter\[[35](https://arxiv.org/html/2610.27159#bib.bib23)\]\. Combined with Theorems[1](https://arxiv.org/html/2610.27159#Thmtheorem1)and[2](https://arxiv.org/html/2610.27159#Thmtheorem2), it gives the variance in closed form for a rare binary outcome\. In practice, the pool is the utility’s weather archive, on which the control is evaluated at no simulation cost\. Part \(c\) then gives the number of archive days that is worth one simulator run\.

### IV\-DTraining Design

Finally, we address the training nature of the learned CVs\. In the rare regime, a learner trained on features drawn from the nominal distributionFXF\_\{X\}sees almost no outage and almost no extreme day\. The following proposition shows that the training features may instead be drawn from any design\.

###### Proposition 2\(Validity under any training design\)\.

Letggbe any function of a training sample whose features are drawn from a distributionqqon𝒳\\mathcal\{X\}and labelled by the simulator, independently of thennestimation pairs drawn fromFXF\_\{X\}\. Then, conditionally ongg, Theorem[3](https://arxiv.org/html/2610.27159#Thmtheorem3)applies to the estimator of Algorithm[1](https://arxiv.org/html/2610.27159#alg1), withρ2\\rho^\{2\}computed underFXF\_\{X\}\. Furthermore, ifFX≪qF\_\{X\}\\ll qwith weightsw=d​FX/d​qw=dF\_\{X\}/dq, the weighted training\-sample ratiosτ^=∑lwl​Yl​g​\(Xl\)/∑lwl​Yl\\hat\{\\tau\}=\\sum\_\{l\}w\_\{l\}Y\_\{l\}g\(X\_\{l\}\)/\\sum\_\{l\}w\_\{l\}Y\_\{l\}andπ^=∑lwl​Yl​g​\(Xl\)/∑lwl​g​\(Xl\)\\hat\{\\pi\}=\\sum\_\{l\}w\_\{l\}Y\_\{l\}g\(X\_\{l\}\)/\\sum\_\{l\}w\_\{l\}g\(X\_\{l\}\)are consistent for the nominal recall and precision\. Hence,τ^​π^\\hat\{\\tau\}\\hat\{\\pi\}is a consistent estimate of the limit \([5](https://arxiv.org/html/2610.27159#S4.E5)\)\.

###### Proof\.

Conditionally ongg, the estimation pairs are i\.i\.d\. from the nominal distribution andggis fixed, so the classical results apply\. The weighted ratios are importance\-sampling estimates\[[26](https://arxiv.org/html/2610.27159#bib.bib10),[44](https://arxiv.org/html/2610.27159#bib.bib11)\]\. ∎

In practice, thentrn\_\{\\mathrm\{tr\}\}expensive training runs are placed where the outages are, by drawing the training features from a shifted design that over\-represents hot days and high forecasts\. Since the simulator returns the severitySS, the label𝟏\{S≥κ\}\\mathbf\{1\}\\\{S\\geq\\kappa\\\}can be formed afterwards for any ratingκ\\kappa\. One training set therefore serves every rarity level\. Score learners are fitted toSSby regression, while the threshold of a binary learner is chosen to maximiseτ^​π^\\hat\{\\tau\}\\hat\{\\pi\}\. For a utility, the simulation planning used to train the control is thus a one\-time investment, which can be reused when the ratings, the reliability targets or the stations under study change, which further emphasize the practical value of our proposed framework\.

## VNumerical Experiments

### V\-ASimulator and Experiment Setup

Each simulated day has a temperature anomalyx1∼N⁡\(0,1\)x\_\{1\}\\sim N\(0,1\), a day\-ahead load\-forecast anomalyx2∼N⁡\(0,1\)x\_\{2\}\\sim N\(0,1\)and a storm flagx3∼Bernoulli⁡\(0\.15\)x\_\{3\}\\sim\\mathrm\{Bernoulli\}\(0\.15\)\. The simulator returns the standardised peak feeder loading

S\\displaystyle S=a1​x1\+a2​max⁡\(x1,0\)2\+a3​x2\\displaystyle=a\_\{1\}x\_\{1\}\+a\_\{2\}\\max\(x\_\{1\},0\)^\{2\}\+a\_\{3\}x\_\{2\}\+a4​x3​max⁡\(ε2,0\)\+a0​ε1,\\displaystyle\+a\_\{4\}x\_\{3\}\\max\(\\varepsilon\_\{2\},0\)\+a\_\{0\}\\varepsilon\_\{1\},\(8\)with internal noiseε1,ε2∼N⁡\(0,1\)\\varepsilon\_\{1\},\\varepsilon\_\{2\}\\sim N\(0,1\), andY=𝟏\{S≥κ\}Y=\\mathbf\{1\}\\\{S\\geq\\kappa\\\}\. The coefficients\(a1,a2,a3,a4,a0\)=\(0\.45,0\.20,0\.55,0\.40,0\.30\)\(a\_\{1\},a\_\{2\},a\_\{3\},a\_\{4\},a\_\{0\}\)=\(0\.45,0\.20,0\.55,0\.40,0\.30\)represent a nonlinear heat effect, a linear forecast effect, a storm effect that only sometimes materialises, and unobserved noise\. The forecast anomaly thus acts as a Gaussian proxy of the load, with latent correlation0\.610\.61\. The ratingκ\\kappais calibrated by quadrature so thatθ∈\{10−3,10−4,10−5,10−6\}\\theta\\in\\\{10^\{\-3\},10^\{\-4\},10^\{\-5\},10^\{\-6\}\\\}exactly\. Since the noise structure is known,η⁡\(x\)\\eta\(x\)is also available by quadrature, and everyρ2\\rho^\{2\}is evaluated exactly by importance sampling over the features\. It is worth noting that our estimators never use this information\.

We compare plain MC with five controls: \(1\) the Gaussian proxy \(GP\)𝟏\{x2≥t\}\\mathbf\{1\}\\\{x\_\{2\}\\geq t\\\}, withttchosen to maximise its exactρ2\\rho^\{2\}; \(2\) a stump \(ST\)𝟏\{x1≥t\}\\mathbf\{1\}\\\{x\_\{1\}\\geq t\\\}, i\.e\., a one\-feature weak learner withttchosen by the weightedτ^​π^\\hat\{\\tau\}\\hat\{\\pi\}criterion; \(3\) a deliberately misspecified linear learner \(LL\), i\.e\., least squares ofSSonXX, with scoreΦ¯​\(\(κ−s^​\(x\)\)/σ^\)\\bar\{\\Phi\}\(\(\\kappa\-\\hat\{s\}\(x\)\)/\\hat\{\\sigma\}\); \(4\) gradient\-boosted trees \(BT\) regressingSSonXX\[[45](https://arxiv.org/html/2610.27159#bib.bib13)\], a standard strong learner for tabular data, with the same score; and \(5\) the Bayes posterior \(BP\)η⁡\(X\)\\eta\(X\), which attains the ceiling\. The learners are trained once onntr=104n\_\{\\mathrm\{tr\}\}=10^\{4\}runs from the tail\-focused designx1∼N⁡\(2,1\.52\)x\_\{1\}\\sim N\(2,1\.5^\{2\}\),x2∼N⁡\(1,1\.52\)x\_\{2\}\\sim N\(1,1\.5^\{2\}\), and are reused at everyθ\\theta\. For comparison, they are also trained on10410^\{4\}nominal runs\. The control means are computed once by importance sampling \(relative standard error below0\.3%0\.3\\%\)\.

### V\-BEstimation Results

Fig\. 1:\(a\) Variance reductionρ2\\rho^\{2\}of each control as the outage becomes rarer\. \(b\) Running estimate atθ=10−4\\theta=10^\{\-4\}as the sample grows: median \(line\) and 95% range \(shaded\) over 400 samples, for Monte Carlo \(black, dashed\) and the boosted\-tree control on the same days \(solid\); The true value is shown as the dotted line\. \(c\) Simulator runs needed to achieve a 10% relative error\.TABLE I:Replication study at about 20 outages per replicationFig\.[1](https://arxiv.org/html/2610.27159#S5.F1)\(a\) and Table[I](https://arxiv.org/html/2610.27159#S5.T1)present the main result\. The Gaussian proxy removes less than0\.2%0\.2\\%of the variance atθ=10−3\\theta=10^\{\-3\}and nothing beyond\. Even a hypothetical proxy with latent correlation0\.980\.98decays from0\.540\.54to0\.380\.38over the range \(grey curves\), as predicted by Proposition[1](https://arxiv.org/html/2610.27159#Thmproposition1)\. The learned controls, in contrast, keep their effect at every level\. The stump keepsρ2≈0\.25\\rho^\{2\}\\approx 0\.25\(1\.31\.3times fewer runs\), and the misspecified linear learner keeps about0\.400\.40\(1\.71\.7times\)\. The boosted trees keep0\.600\.60to0\.640\.64\(2\.52\.5to2\.82\.8times\), which is within5%5\\%of the ceiling of0\.630\.63to0\.650\.65set by the features\. Table[I](https://arxiv.org/html/2610.27159#S5.T1)repeats each estimator over 400, 300, 120 and 60 replications at the four rarity levels, withn=20/θn=20/\\thetasimulated days per replication, i\.e\., about 20 outages\. For each estimator, it reports the exactρ2\\rho^\{2\}, the equivalent\-runs multiplier1/\(1−ρ2\)1/\(1\-\\rho^\{2\}\), the empirical relative RMSE and its theoretical value\(1−ρ2\)​\(1−θ\)/\(n​θ\)\\sqrt\{\(1\-\\rho^\{2\}\)\(1\-\\theta\)/\(n\\theta\)\}from \([3](https://arxiv.org/html/2610.27159#S3.E3)\), the relative bias, and the coverage of the plug\-in95%95\\%interval\. In Table[I](https://arxiv.org/html/2610.27159#S5.T1), every empirical relative RMSE is within two standard errors of its theoretical value, e\.g\.,0\.210\.21to0\.240\.24for MC against0\.2240\.224\. The boosted\-tree control reduces it to0\.130\.13to0\.140\.14, while the95%95\\%intervals keep a coverage between0\.920\.92and0\.980\.98\. Atθ=10−6\\theta=10^\{\-6\}, a10%10\\%relative error requires3\.8×1083\.8\\times 10^\{8\}runs with MC and1\.5×1081\.5\\times 10^\{8\}with the boosted\-tree control, for a one\-time training cost of1\.2×1041\.2\\times 10^\{4\}runs\.

Fig\.[1](https://arxiv.org/html/2610.27159#S5.F1)\(b\) translates this result into more practical insights\. Atθ=10−4\\theta=10^\{\-4\}, i\.e\.,0\.0370\.037expected blackout days per year, the95%95\\%range of the MC estimate over 400 plannings is still±0\.021\\pm 0\.021days per year after10510^\{5\}simulated days and±0\.007\\pm 0\.007after10610^\{6\}\. The control variate reaches the same ranges after about5×1045\\times 10^\{4\}and4×1054\\times 10^\{5\}days, and neither estimator shows a bias\. Very short replications are the exception\. With only about two observed outages,β^\\hat\{\\beta\}is fitted to a handful of positive days, so the control may increase the variance and its plug\-in interval undercovers\. From about twenty outages onward, the coverage is close to nominal\.

Furthermore, it appears That the training design plays crucial role as much as the learner\. The same boosted learner trained on10410^\{4\}nominal runs hasρ2=0\.34\\rho^\{2\}=0\.34atθ=10−3\\theta=10^\{\-3\},0\.050\.05at10−410^\{\-4\}and zero below, since its training set contains no extreme day and tree ensembles extrapolate as a constant\. The tail\-focused design, in contrast, leaves hundreds of outage days in the training set even atθ=10−6\\theta=10^\{\-6\}\. Finally, with the control mean taken from a pool ofm=10​nm=10nunlabelled draws, the loss relative to a known mean is about a tenth ofρ2\\rho^\{2\}, as \([7](https://arxiv.org/html/2610.27159#S4.E7)\) predicts\.

### V\-CCase Study: Blackout Map of Surabaya

We now apply the estimator to a planning decision on the Surabaya network\[[8](https://arxiv.org/html/2610.27159#bib.bib1)\], which has 11 candidate feeders and 31 sub\-districts\. For each feederjj, the earlier work reports a masked daily peak load with meanμj\\mu\_\{j\}and standard deviationσj\\sigma\_\{j\}, a ratingzjz\_\{j\}, and a reliabilitypjp\_\{j\}between0\.9740\.974and0\.9830\.983at the current headroomh0​j=\(zj−μj\)/σj≈2h\_\{0j\}=\(z\_\{j\}\-\\mu\_\{j\}\)/\\sigma\_\{j\}\\approx 2\. We place a charging station on every feeder\. The planner then chooses the headroomhjh\_\{j\}, through a feeder upgrade or a cap on the charging load, so that the daily overload probabilityθj​\(hj\)\\theta\_\{j\}\(h\_\{j\}\)meets the one\-day\-in\-ten\-years targetθ⋆=2\.7×10−4\\theta^\{\\star\}=2\.7\\times 10^\{\-4\}\. The daily peak load is modelled asμj\+σj​S~j\\mu\_\{j\}\+\\sigma\_\{j\}\\tilde\{S\}\_\{j\}, whereS~j\\tilde\{S\}\_\{j\}is the severity \([8](https://arxiv.org/html/2610.27159#S5.E8)\) standardised so thatℙ⁡\(S~j≥h0​j\)=1−pj\\mathbb\{P\}\(\\tilde\{S\}\_\{j\}\\geq h\_\{0j\}\)=1\-p\_\{j\}\. The simulator therefore reproduces the reported reliability of every feeder and differs from the Gaussian model only in the tail\.

We compare three design procedures: \(a\) the Gaussian load model, which requiresh=Φ¯−1​\(θ⋆\)=3\.46h=\\bar\{\\Phi\}^\{\-1\}\(\\theta^\{\\star\}\)=3\.46for every feeder; \(b\) imprecise estimates, i\.e\., plain MC withn=104n=10^\{4\}simulated days per feeder \(the largest sample of the earlier work\) and the reserve sized to the point estimate; and \(c\) precise estimates, i\.e\., the boosted\-tree control variate on the same days, with the reserve sized to the smallest headroom whose95%95\\%upper confidence limit meets the target\. Since the simulator returns the severity, one replication gives the whole curveθj​\(h\)\\theta\_\{j\}\(h\)\.

![Refer to caption](https://arxiv.org/html/2610.27159v1/fig_map.png)Fig\. 2:Expected blackout days per year after a charging station is added and the reserve is sized for one day in ten years by \(a\) the Gaussian load model of the earlier work, \(b\) Monte Carlo with10410^\{4\}days and the point estimate, and \(c\) the learned control variate on the same days with a 95% upper limit\.Fig\.[2](https://arxiv.org/html/2610.27159#S5.F2)shows the true expected blackout days per year of each design\. Table[II](https://arxiv.org/html/2610.27159#S5.T2)summarises 200 replications per feeder by the city total of feeder\-blackout days per year, the share of feeders whose design ends above the target and above twice the target, the 90th percentile of the blackout days per year of a feeder \(averaged over the feeders\), and the share of the population served by a feeder above the target\. Under the Gaussian model, every feeder is under\-designed\. The required headroom is4\.054\.05to4\.464\.46, and the true blackout rate is0\.310\.31to0\.660\.66days per year, i\.e\., three to seven times the target, for a city total of5\.75\.7feeder\-blackout days per year\. This is the design consequence of Proposition[1](https://arxiv.org/html/2610.27159#Thmproposition1): the Gaussian model agrees with the simulator at the current operating point but underestimates the reserve needed in the tail\. With imprecise estimates \(panel \(b\)\),46%46\\%of the feeders end above the target and8%8\\%above twice the target\. Moreover, the planner cannot tell which feeders these are\. In the example shown, for instance, some districts \(Darmo Permai, Kenjeran and Rayon Ploso\) are at0\.140\.14to0\.150\.15days per year, while other district \(Embong Wungu\) is over\-built at0\.020\.02\. With precise estimates \(panel \(c\)\),15%15\\%of the feeders remain above the target\. The city total falls from1\.171\.17to0\.770\.77feeder\-blackout days per year, while the share of the population on a feeder above the target falls from46%46\\%to15%15\\%\. Table[II](https://arxiv.org/html/2610.27159#S5.T2)also breaks down the contribution the two components\. The confidence limit protects the design, since MC with the same limit reaches0\.780\.78city days\. The control variate then reduces the remaining failures, from24%24\\%to15%15\\%of the feeders above the target, and the 90th percentile from0\.140\.14to0\.110\.11days per year\.

TABLE II:City\-wide outcome over 200 replications of10,00010\{,\}000daysThe cost side of the planners is illustrated on the largest feeder \(Embong Wungu, withμ=28,979\\mu=28\{,\}979,σ=3,622\\sigma=3\{,\}622, rating36,22436\{,\}224\)\. The exact requirement is an upgrade of8,4158\{,\}415units \(29%29\\%of the base load\), whereas the Gaussian model asks for5,2745\{,\}274\(18%18\\%\)\. Under the confidence\-limit rule, MC over\-designs the feeder by1,4201\{,\}420,675675and330330units on average for samples of2×1042\\times 10^\{4\},10510^\{5\}and5×1055\\times 10^\{5\}days, against816816,459459and240240units for the control variate\. The control\-variate sample thus matches the decision quality of an MC sample about two to three times larger, close to the factor1/\(1−ρ2\)≈2\.71/\(1\-\\rho^\{2\}\)\\approx 2\.7from \([3](https://arxiv.org/html/2610.27159#S3.E3)\)\. Both estimators keep the true outage probability below the target in at least93%93\\%of the samples\.

## VIDiscussion

### VI\-AImplications for Grid Operators and City Planners

The results lead to a handful of implications for a distribution utility or a city that plans its EV infrastructure\.

First, reliability targets could be certified with far fewer simulations\. Since the planning quantities are linear inθ\\theta\(Section[III](https://arxiv.org/html/2610.27159#S3)\), the variance reduction transfers to them unchanged\. In the Surabaya case, the control variate reaches the decision quality of an MC sample two to three times larger\. Combined with an upper confidence limit, it also reduces the city\-wide blackout days by a third, and the population on an under\-designed feeder from46%46\\%to15%15\\%, compared with MC point estimates\. Furthermore, the load model plays important role as much as the estimator\. The Gaussian model of the earlier work reproduces the current outage rate of some districts \(including Embong Wungu\), yet it underestimates the required upgrade by about3,1003\{,\}100units\. The resulting risk is more than five times the accepted level\.

Next, investment in features should come before investment in models\. By Theorem[2](https://arxiv.org/html/2610.27159#Thmtheorem2), the information collected about each day bounds what any learner could achieve\. In our simulator, halving or doubling the unobserved noise moved the ceiling from0\.640\.64to0\.790\.79or0\.380\.38for every learner\. Hence, feeder loading histories, maintenance schedules, distributed\-generation output and the load forecasts that utilities already produce\[[40](https://arxiv.org/html/2610.27159#bib.bib36),[41](https://arxiv.org/html/2610.27159#bib.bib37)\]could raise the ceiling, while a more elaborate model on the same features could not\.

Furthermore, the training budget should be spent where the outages are\. A training sample placed on hot, high\-forecast days keeps its value at every reliability target and can be reused when the ratings or the targets change\. A sample placed on typical days, in contrast, is wasted in the rare regime\. Finally, the control depends only on the features of a day, so it can be evaluated on a decade of weather records at no simulation cost\. Theorem[3](https://arxiv.org/html/2610.27159#Thmtheorem3)\(c\) then gives the number of such days that is worth one simulator run\. Together, these points suggest a simple protocol for an operator: train once on a tail\-focused sample, evaluate the control on the archive, size the reserve with an upper confidence limit, and commit only after the sample contains at least about twenty outages\.

### VI\-BLimitations and Future Work

We also highlight a few limitations of this work and future work directions\. First, we are limited by the use of a synthetic simulator and synthetic features, which are chosen so that the exact posterior is available for evaluation\. In the Surabaya case, the simulator matches the reported reliability of every feeder, but the shape of its tail is assumed\. A feeder model driven by real weather and outage records is therefore the natural next step\. Furthermore, outages at different feeders are modelled independently, although a city\-wide storm correlates them, which we ignore in this study\. Finally, the tail\-focused design was also chosen by hand\. An adaptive design could alleviate this choice, which might allow us to combine our current learned control variate results with learned importance\-sampling proposal\[[31](https://arxiv.org/html/2610.27159#bib.bib26)\], which we plan to explore in our future work\.

## VIIConclusion

In this paper, we have shown that a classifier of grid outages trained on cheap observable features yields a valid control variate, which reduces the variance for every weak learner\. Its effect is governed by recall and precision through an exact formula and is bounded by a ceiling that the features alone determine\. Unlike a Gaussian proxy of the load, it does not vanish as the outage becomes rarer, provided that the learner is trained where the outages are\. On a feeder\-overload simulator, the learned controls removed between a quarter and nearly two thirds of the Monte Carlo variance at outage probabilities from10−310^\{\-3\}to10−610^\{\-6\}\. On the Surabaya network, the Gaussian load model yields every feeder three to seven times the one\-day\-in\-ten\-years target, while point\-estimate designs from plain Monte Carlo makes almost half of the feeders achieve above the target\. At the same simulation budget, the learned control variate with a confidence limit also reduces the city\-wide blackout days by a third\. Compared with Monte Carlo under the same rule, it also frees up the over\-design of the reserve by up to40%40\\%\. The proposed method thus makes the certification of rare overload probabilities affordable for a utility, at the cost of a one\-time training sample and a pool of historical feature data, which opens huge potential to design a more reliable grid against rare\-event outages\.

## References

- \[1\]\(2018\)Impact of uncoordinated plug\-in electric vehicle charging on residential power demand\.Nature Energy3\(3\),pp\. 193–201\.External Links:[Document](https://dx.doi.org/10.1038/s41560-017-0074-z)Cited by:[§I](https://arxiv.org/html/2610.27159#S1.p1.1),[§II\-A](https://arxiv.org/html/2610.27159#S2.SS1.p1.1)\.
- \[2\]K\. Clement\-Nyns, E\. Haesen, and J\. Driesen\(2010\)The impact of charging plug\-in hybrid electric vehicles on a residential distribution grid\.IEEE Transactions on Power Systems25\(1\),pp\. 371–380\.External Links:[Document](https://dx.doi.org/10.1109/TPWRS.2009.2036481)Cited by:[§I](https://arxiv.org/html/2610.27159#S1.p1.1),[§II\-A](https://arxiv.org/html/2610.27159#S2.SS1.p1.1)\.
- \[3\]S\. Powell, G\. V\. Cezar, L\. Min, I\. M\. L\. Azevedo, and R\. Rajagopal\(2022\)Charging infrastructure access and operation to reduce the grid impacts of deep electric vehicle adoption\.Nature Energy7\(10\),pp\. 932–945\.External Links:[Document](https://dx.doi.org/10.1038/s41560-022-01105-7)Cited by:[§I](https://arxiv.org/html/2610.27159#S1.p1.1),[§II\-A](https://arxiv.org/html/2610.27159#S2.SS1.p1.1)\.
- \[4\]I\. Nutkani, H\. Toole, N\. Fernando, and L\. P\. C\. Andrew\(2024\)Impact of ev charging on electrical distribution network and mitigating solutions–a review\.IET Smart Grid7\(5\),pp\. 485–502\.Cited by:[§I](https://arxiv.org/html/2610.27159#S1.p1.1),[§II\-A](https://arxiv.org/html/2610.27159#S2.SS1.p1.1)\.
- \[5\]J\. Sun, X\. Hu, H\. Liu, Z\. Zhang, H\. Peng, and Z\. Yang\(2025\)Hosting capacity assessment of distribution networks with PV and EV charging stations\.In2025 9th International Conference on Smart Grid and Smart Cities \(ICSGSC\),pp\. 155–159\.Cited by:[§I](https://arxiv.org/html/2610.27159#S1.p1.1),[§II\-A](https://arxiv.org/html/2610.27159#S2.SS1.p1.1)\.
- \[6\]J\. Lukuyu, R\. Shirley, and J\. Taneja\(2024\)Managing grid impacts from increased electric vehicle adoption in african cities\.Scientific Reports14\(1\),pp\. 24320\.Cited by:[§I](https://arxiv.org/html/2610.27159#S1.p1.1),[§II\-A](https://arxiv.org/html/2610.27159#S2.SS1.p1.1)\.
- \[7\]M\. O\. Khan, S\. Kirmani, and M\. Rihan\(2024\)Impact assessment of electric vehicle charging on distribution networks\.Renewable Energy Focus50,pp\. 100599\.Cited by:[§I](https://arxiv.org/html/2610.27159#S1.p1.1),[§II\-A](https://arxiv.org/html/2610.27159#S2.SS1.p1.1)\.
- \[8\]Y\. Akhra, I\. Vanany, and M\. Arief\(2026\)A robust electric vehicle charging station optimization under rare event disruptions\.Note:arXiv:2307\.05470vExternal Links:[Document](https://dx.doi.org/10.48550/arXiv.2307.05470)Cited by:[§I](https://arxiv.org/html/2610.27159#S1.p2.1),[§I](https://arxiv.org/html/2610.27159#S1.p3.1),[§I](https://arxiv.org/html/2610.27159#S1.p5.1),[§II\-B](https://arxiv.org/html/2610.27159#S2.SS2.p1.1),[§III\-A](https://arxiv.org/html/2610.27159#S3.SS1.p3.1),[§III\-C](https://arxiv.org/html/2610.27159#S3.SS3.p2.1),[§V\-C](https://arxiv.org/html/2610.27159#S5.SS3.p1.1)\.
- \[9\]I\. Dobson, B\. A\. Carreras, V\. E\. Lynch, and D\. E\. Newman\(2007\)Complex systems analysis of series of blackouts: cascading failure, critical points, and self\-organization\.Chaos17\(2\),pp\. 026103\.External Links:[Document](https://dx.doi.org/10.1063/1.2737822)Cited by:[§I](https://arxiv.org/html/2610.27159#S1.p2.1),[§II\-A](https://arxiv.org/html/2610.27159#S2.SS1.p1.1),[§III\-A](https://arxiv.org/html/2610.27159#S3.SS1.p2.1)\.
- \[10\]P\. Hines, J\. Apt, and S\. Talukdar\(2009\)Large blackouts in North America: historical trends and policy implications\.Energy Policy37\(12\),pp\. 5249–5259\.External Links:[Document](https://dx.doi.org/10.1016/j.enpol.2009.07.049)Cited by:[§I](https://arxiv.org/html/2610.27159#S1.p2.1),[§II\-A](https://arxiv.org/html/2610.27159#S2.SS1.p1.1)\.
- \[11\]U\.S\.–Canada Power System Outage Task Force\(2004\)Final report on the August 14, 2003 blackout in the United States and Canada: causes and recommendations\.Technical reportU\.S\. Department of Energy and Natural Resources Canada\.Cited by:[§I](https://arxiv.org/html/2610.27159#S1.p2.1),[§II\-A](https://arxiv.org/html/2610.27159#S2.SS1.p1.1)\.
- \[12\]R\. Billinton and R\. N\. Allan\(1996\)Reliability evaluation of power systems\.2nd edition,Plenum Press,New York\.Cited by:[§I](https://arxiv.org/html/2610.27159#S1.p2.1),[§II\-A](https://arxiv.org/html/2610.27159#S2.SS1.p1.1),[§III\-A](https://arxiv.org/html/2610.27159#S3.SS1.p2.1)\.
- \[13\]R\. E\. Schapire\(1990\)The strength of weak learnability\.Machine Learning5\(2\),pp\. 197–227\.Cited by:[§I](https://arxiv.org/html/2610.27159#S1.p4.1),[§II\-D](https://arxiv.org/html/2610.27159#S2.SS4.p1.1),[Definition 1](https://arxiv.org/html/2610.27159#Thmdefinition1.p1.1)\.
- \[14\]K\. H\. LaCommare and J\. H\. Eto\(2006\)Cost of power interruptions to electricity consumers in the United States \(US\)\.Energy31\(12\),pp\. 1845–1855\.External Links:[Document](https://dx.doi.org/10.1016/j.energy.2006.02.008)Cited by:[§II\-A](https://arxiv.org/html/2610.27159#S2.SS1.p1.1),[§III\-A](https://arxiv.org/html/2610.27159#S3.SS1.p2.1)\.
- \[15\]M\. S\. Daskin\(2013\)Network and discrete location: models, algorithms, and applications\.2nd edition,Wiley\.Cited by:[§II\-B](https://arxiv.org/html/2610.27159#S2.SS2.p1.1)\.
- \[16\]I\. Frade, A\. Ribeiro, G\. Gonçalves, and A\. P\. Antunes\(2011\)Optimal location of charging stations for electric vehicles in a neighborhood in lisbon, portugal\.Transportation Research Record2252\(1\),pp\. 91–98\.Cited by:[§II\-B](https://arxiv.org/html/2610.27159#S2.SS2.p1.1)\.
- \[17\]N\. Amilia, Z\. Palinrungi, I\. Vanany, and M\. Arief\(2022\)Designing an optimized electric vehicle charging station infrastructure for urban area: a case study from indonesia\.In2022 IEEE 25th International Conference on Intelligent Transportation Systems \(ITSC\),pp\. 2812–2817\.Cited by:[§II\-B](https://arxiv.org/html/2610.27159#S2.SS2.p1.1)\.
- \[18\]M\. Hosseini and S\. A\. MirHassani\(2015\)Refueling\-station location problem under uncertainty\.Transportation Research Part E: Logistics and Transportation Review84,pp\. 101–116\.Cited by:[§II\-B](https://arxiv.org/html/2610.27159#S2.SS2.p1.1)\.
- \[19\]B\. Yıldız, E\. Olcaytu, and A\. Şen\(2019\)The urban recharging infrastructure design problem with stochastic demands and capacitated charging stations\.Transportation Research Part B: Methodological119,pp\. 22–44\.Cited by:[§II\-B](https://arxiv.org/html/2610.27159#S2.SS2.p1.1)\.
- \[20\]Q\. Zhang, H\. Yu, G\. Zhang, and T\. Ma\(2023\)Optimal planning of flood\-resilient electric vehicle charging stations\.Computer\-Aided Civil and Infrastructure Engineering38\(4\),pp\. 489–507\.Cited by:[§II\-B](https://arxiv.org/html/2610.27159#S2.SS2.p1.1)\.
- \[21\]A\. Hussain, V\. Bui, and H\. Kim\(2020\)Optimal sizing of battery energy storage system in a fast ev charging station considering power outages\.IEEE Transactions on Transportation Electrification6\(2\),pp\. 453–463\.Cited by:[§II\-B](https://arxiv.org/html/2610.27159#S2.SS2.p1.1)\.
- \[22\]T\. Unterluggauer, J\. Rich, P\. B\. Andersen, and S\. Hashemi\(2022\)Electric vehicle charging infrastructure planning for integrated transportation and power distribution networks: a review\.ETransportation,pp\. 100163\.Cited by:[§II\-B](https://arxiv.org/html/2610.27159#S2.SS2.p1.1)\.
- \[23\]L\. V\. Snyder and M\. S\. Daskin\(2005\)Reliability models for facility location: the expected failure cost case\.Transportation Science39\(3\),pp\. 400–416\.Cited by:[§II\-B](https://arxiv.org/html/2610.27159#S2.SS2.p1.1)\.
- \[24\]S\. Asmussen and P\. W\. Glynn\(2007\)Stochastic simulation: algorithms and analysis\.Springer\.Cited by:[§II\-C](https://arxiv.org/html/2610.27159#S2.SS3.p1.1)\.
- \[25\]P\. Heidelberger\(1995\)Fast simulation of rare events in queueing and reliability models\.ACM Transactions on Modeling and Computer Simulation5\(1\),pp\. 43–85\.Cited by:[§II\-C](https://arxiv.org/html/2610.27159#S2.SS3.p1.1)\.
- \[26\]H\. Kahn and A\. W\. Marshall\(1953\)Methods of reducing sample size in Monte Carlo computations\.Journal of the Operations Research Society of America1\(5\),pp\. 263–278\.Cited by:[§II\-C](https://arxiv.org/html/2610.27159#S2.SS3.p1.1),[§IV\-D](https://arxiv.org/html/2610.27159#S4.SS4.p2.1.1)\.
- \[27\]S\. S\. Lavenberg and P\. D\. Welch\(1981\)A perspective on the use of control variables to increase the efficiency of Monte Carlo simulations\.Management Science27\(3\),pp\. 322–335\.Cited by:[§II\-C](https://arxiv.org/html/2610.27159#S2.SS3.p1.1),[§III\-C](https://arxiv.org/html/2610.27159#S3.SS3.p1.3),[§IV\-C](https://arxiv.org/html/2610.27159#S4.SS3.p2.1.1)\.
- \[28\]B\. L\. Nelson\(1990\)Control variate remedies\.Operations Research38\(6\),pp\. 974–992\.Cited by:[§II\-C](https://arxiv.org/html/2610.27159#S2.SS3.p1.1),[§IV\-C](https://arxiv.org/html/2610.27159#S4.SS3.p2.1.1)\.
- \[29\]P\. Glasserman\(2004\)Monte carlo methods in financial engineering\.Springer\.Cited by:[§II\-C](https://arxiv.org/html/2610.27159#S2.SS3.p1.1),[§III\-C](https://arxiv.org/html/2610.27159#S3.SS3.p1.3)\.
- \[30\]T\. C\. Hesterberg and B\. L\. Nelson\(1998\)Control variates for probability and quantile estimation\.Management Science44\(9\),pp\. 1295–1312\.Cited by:[§II\-C](https://arxiv.org/html/2610.27159#S2.SS3.p1.1)\.
- \[31\]M\. Arief, Y\. Bai, W\. Ding, S\. He, Z\. Huang, H\. Lam, and D\. Zhao\(2021\)Deep probabilistic accelerated evaluation: a robust certifiable rare\-event simulation methodology for black\-box safety\-critical systems\.InInternational Conference on Artificial Intelligence and Statistics \(AISTATS\),pp\. 595–603\.Cited by:[§II\-C](https://arxiv.org/html/2610.27159#S2.SS3.p1.1),[§VI\-B](https://arxiv.org/html/2610.27159#S6.SS2.p1.1)\.
- \[32\]B\. Peherstorfer, K\. Willcox, and M\. Gunzburger\(2016\)Optimal model management for multifidelity Monte Carlo estimation\.SIAM Journal on Scientific Computing38\(5\),pp\. A3163–A3194\.Cited by:[§II\-C](https://arxiv.org/html/2610.27159#S2.SS3.p1.1),[§IV\-C](https://arxiv.org/html/2610.27159#S4.SS3.p3.1)\.
- \[33\]B\. Peherstorfer, K\. Willcox, and M\. Gunzburger\(2018\)Survey of multifidelity methods in uncertainty propagation, inference, and optimization\.SIAM Review60\(3\),pp\. 550–591\.Cited by:[§II\-C](https://arxiv.org/html/2610.27159#S2.SS3.p1.1)\.
- \[34\]A\. N\. Angelopoulos, S\. Bates, C\. Fannjiang, M\. I\. Jordan, and T\. Zrnic\(2023\)Prediction\-powered inference\.Science382\(6671\),pp\. 669–674\.Cited by:[§II\-C](https://arxiv.org/html/2610.27159#S2.SS3.p1.1)\.
- \[35\]A\. N\. Angelopoulos, J\. C\. Duchi, and T\. Zrnic\(2023\)PPI\+\+: efficient prediction\-powered inference\.arXiv preprint arXiv:2311\.01453\.Cited by:[§II\-C](https://arxiv.org/html/2610.27159#S2.SS3.p1.1),[§IV\-C](https://arxiv.org/html/2610.27159#S4.SS3.p3.1)\.
- \[36\]C\. J\. Oates, M\. Girolami, and N\. Chopin\(2017\)Control functionals for Monte Carlo integration\.Journal of the Royal Statistical Society: Series B79\(3\),pp\. 695–718\.Cited by:[§II\-C](https://arxiv.org/html/2610.27159#S2.SS3.p1.1)\.
- \[37\]L\. F\. South, C\. J\. Oates, A\. Mira, and C\. Drovandi\(2023\)Regularized zero\-variance control variates\.Bayesian Analysis18\(3\),pp\. 865–888\.Cited by:[§II\-C](https://arxiv.org/html/2610.27159#S2.SS3.p1.1)\.
- \[38\]S\. D\. Guikema, R\. Nateghi, S\. M\. Quiring, A\. Staid, A\. C\. Reilly, and M\. Gao\(2014\)Predicting hurricane power outages to support storm response planning\.IEEE Access2,pp\. 1364–1373\.Cited by:[§II\-D](https://arxiv.org/html/2610.27159#S2.SS4.p1.1)\.
- \[39\]R\. Nateghi, S\. D\. Guikema, and S\. M\. Quiring\(2014\)Power outage estimation for tropical cyclones: improved accuracy with simpler models\.Risk Analysis34\(6\),pp\. 1069–1078\.Cited by:[§II\-D](https://arxiv.org/html/2610.27159#S2.SS4.p1.1)\.
- \[40\]X\. Hu, J\. Sun, H\. Liu, Z\. Zhang, F\. Qiu, and Z\. Yang\(2025\)Electric vehicle load forecasting based on probabilistic modeling and deep learning\.In2025 9th International Conference on Smart Grid and Smart Cities \(ICSGSC\),pp\. 456–461\.Cited by:[§II\-D](https://arxiv.org/html/2610.27159#S2.SS4.p1.1),[§VI\-A](https://arxiv.org/html/2610.27159#S6.SS1.p3.1)\.
- \[41\]K\. Chen, Y\. Zhang, Q\. Dai, and H\. Yu\(2025\)Improving the operational reliability of smart grids: an ultra\-short\-term load forecasting scheme based on anisotropic radial basis function kernel extreme learning machine\.In2025 9th International Conference on Smart Grid and Smart Cities \(ICSGSC\),pp\. 422–427\.Cited by:[§II\-D](https://arxiv.org/html/2610.27159#S2.SS4.p1.1),[§VI\-A](https://arxiv.org/html/2610.27159#S6.SS1.p3.1)\.
- \[42\]M\. Kearns and L\. Valiant\(1994\)Cryptographic limitations on learning Boolean formulae and finite automata\.Journal of the ACM41\(1\),pp\. 67–95\.Cited by:[§II\-D](https://arxiv.org/html/2610.27159#S2.SS4.p1.1),[Definition 1](https://arxiv.org/html/2610.27159#Thmdefinition1.p1.1)\.
- \[43\]M\. Sibuya\(1960\)Bivariate extreme statistics, I\.Annals of the Institute of Statistical Mathematics11\(2\),pp\. 195–210\.Cited by:[§IV\-B](https://arxiv.org/html/2610.27159#S4.SS2.p6.1)\.
- \[44\]H\. Shimodaira\(2000\)Improving predictive inference under covariate shift by weighting the log\-likelihood function\.Journal of Statistical Planning and Inference90\(2\),pp\. 227–244\.Cited by:[§IV\-D](https://arxiv.org/html/2610.27159#S4.SS4.p2.1.1)\.
- \[45\]F\. Pedregosa, G\. Varoquaux, A\. Gramfort, V\. Michel, B\. Thirion, O\. Grisel, M\. Blondel, P\. Prettenhofer, R\. Weiss, V\. Dubourg,et al\.\(2011\)Scikit\-learn: machine learning in Python\.Journal of Machine Learning Research12,pp\. 2825–2830\.Cited by:[§V\-A](https://arxiv.org/html/2610.27159#S5.SS1.p2.1)\.