Accelerating System Adequacy Assessment using the Multilevel Monte Carlo Approach
aa r X i v : . [ s t a t . C O ] A p r Accelerating System Adequacy Assessmentusing the Multilevel Monte Carlo Approach
Simon Tindemans
Department of Electrical Sustainable EnergyDelft University of TechnologyDelft, The [email protected]
Goran Strbac
Department of Electrical and Electronic EngineeringImperial College LondonLondon, [email protected]
Abstract —Accurately and efficiently estimating system per-formance under uncertainty is paramount in power systemplanning and operation. Monte Carlo simulation is often usedfor this purpose, but convergence may be slow, especially whendetailed models are used. Previously published methods tospeed up computations may severely constrain model complex-ity, limiting their real-world effectiveness. This paper uses therecently proposed Multilevel Monte Carlo (MLMC) framework,which combines outputs from a hierarchy of simulators toboost computational efficiency without sacrificing accuracy. Itexplains which requirements the MLMC framework imposes onthe model hierarchy, and how these naturally occur in powersystem adequacy assessment problems. Two adequacy assessmentexamples are studied in detail: a composite system and a systemwith heterogeneous storage units. An intuitive speed metric isintroduced for easy comparison of simulation setups. Dependingon the problem and metric of interest, large speedups can beobtained.
Index Terms —adequacy assessment, computational efficiency,Monte Carlo methods, storage dispatch, time-sequential simula-tion
I. I
NTRODUCTION
Operational and planning problems in the power systemdomain often involve the assessment of (sub-)system perfor-mance across a range of probabilistically modelled scenarios.For all but the simplest power system models, this cannotbe done analytically, and Monte Carlo (MC) simulations areused instead. MC simulations are a powerful general purposecomputation method with a long tradition in power systemapplications [1], but convergence to the correct answer may beslow. A number of different variance reduction methods existto speed up convergence of Monte Carlo estimates, e.g. [1],[2]. One of these, importance sampling, has recently grown inpopularity for power system applications, especially in combi-nation with automatic tuning of model bias parameters usingthe cross-entropy approach [3], [4]. However, implementingimportance sampling typically requires deep insight into themodel, and limits the design freedom, e.g. for simulationsinvolving complex decision making or sequential actions.
This research was supported by the SMART-SAFE project, funded throughthe HubNet (Extension) programme (EPSRC EP/N030028/1).
The Multilevel Monte Carlo (MLMC) method was intro-duced in the context of computational finance to speed up av-eraging over sample paths, without compromising model detailor accuracy [5]. Initial applications involved the combinationof multi-resolution models (geometric sequences), but otherapplications have subsequently evolved. A good overview ofthe method and its applications is given in [5]. The MLMCapproach has recently been used in a reliability context tospeed up the estimation of the average mission time of largesystems in [6]. In [7], electrical distribution system risk metricswere estimated using MLMC, using a multi-scale approach tosimulate component failures and repairs.This paper considers how the MLMC framework [5] can beused to accelerate risk calculations, in particular in applicationsrelating to system adequacy assessment of complex systems.The contributions of this work are as follows.1) A concise overview of the MLMC approach to theestimation of risks is given. It is shown how the struc-ture required for MLMC simulation naturally occursin adequacy assessment problems, and can often beimplemented with minimal changes to the constituentmodels. Two examples of common model patterns aregiven.2) An intuitive speed metric is introduced that allowsfor fair comparison between Monte Carlo simulationapproaches, and across risk measures.3) Two case studies are presented, each representing oneof the common model patterns. The MLMC approachresults in large speedups, in one case speeding upsimulations by a factor 2000 compared to conventionalMonte Carlo sampling. The sensitivity of computationalspeed to the model stack is investigated.II. M
ETHODOLOGY
A. Mathematical problem statement
Power system performance indicators often take the formof risk measures q that are expressed as the expectation of a The framework of estimating expectation values is less limiting than itmay seem. For example, if one is interested in estimating the distribution of X in addition to the expectation E [ X ] , one can define a series of quantities X ( v ) := I X ≤ v , so that E (cid:2) X ( v ) (cid:3) = F X ( v ) . erformance indicator X (a random variable), i.e. q = E [ X ] .Formally, the random variable X may be seen as a function X : Ω → R that associates a numerical outcome with everysystem state ω ∈ Ω in a sample space Ω . The probabilisticbehaviour of the system, and therefore of X , is defined byassociating probabilities (sets of) states.In the context of system adequacy assessment, the proba-bilistic behaviour of a power system is typically specified usinga bottom up model that defines demand levels, componentstatus, generator output levels, etc. This model generates boththe sample space Ω (the set of all possible combinationsof component states) and the associated probabilities. Thefunction X deterministically evaluates any specific state ω ∈ Ω and computes a numerical performance measure for that state.The risk measure q = E [ X ] is then the (probability weighted)average of the function X over all states.For even moderately complex systems, it is not possibleto compute the quantity of interest q = E [ X ] analytically,nor can it be computed by enumeration of all states in Ω . In such cases, it is common to resort to Monte Carlosimulation, in which power system states ω ( i ) , i = 1 , , . . . are generated using the probabilistic bottom-up model andanalysed to provide relevant outcomes X ( ω ( i ) ) . It should benoted that at any time, multiple outcomes X ( a ) , X ( b ) , . . . (e.g.number of outages, energy not supplied) can be measuredsimultaneously, at little to no extra cost. In the mathematicalanalysis that follows, only a single risk measure q = E [ X ] isdiscussed, but the methods can trivially be applied in parallel. B. Conventional Monte Carlo
A brief summary of conventional Monte Carlo simulationis given in this section, as a point of reference for followingsections. In conventional Monte Carlo simulation, the quantity q = E [ X ] is approximated by the Monte Carlo estimator ˆ Q MC ≡ n n X i =1 X ( i ) , (1)where { X (1) , . . . , X ( n ) } represents a random sample from X , with each X ( i ) independent and identically distributed to X . Note that we distinguish the random variable X ( i ) thatrepresents the i -th random draw from X , and its realisation x ( i ) in a particular experiment or simulation run. The MCestimate for a simulation run is thus given by q ≈ ˆ q = 1 n n X i =1 x ( i ) . (2)We proceed to use the generic expression (1) to reason aboutthe convergence of the result. The error ∆ Q MC obtained inthis approximation is ∆ Q MC = ˆ Q MC − q. (3) We use the statistics convention that a sample is a set of sampled values,rather than the computational science convention where each x ( i ) is a sample. The MC estimator ˆ Q MC is unbiased, and, as a result of thecentral limit theorem, for a sufficiently large sample size n , ∆ Q MC is normally distributed, so that ∆ Q MC ∼ N (cid:16) , σ Q MC (cid:17) . (4)The variance of ˆ Q MC follows from the MC estimator (1): σ Q MC = σ X n . (5)As a result, the standard error σ ˆ Q MC equals σ X / √ n , indicat-ing the typical O ( n − / ) convergence of MC simulations.To quantify the computational efficiency of an MC simula-tion, we denote by τ the average time required to generate asingle realisation x ( i ) . The time spent to generate a sample ofsize n is then t MC = nτ. (6)Using this relation, the variance (5) can be expressed as σ Q MC ( t MC ) = σ X τt MC . (7) C. Multilevel Monte Carlo
For multilevel Monte Carlo (MLMC), we assume to haveat our disposal a hierarchy of models M , . . . , M L that gen-erate random outputs X , . . . , X L , the expectations of whichapproximate E [ X ] with increasing accuracy. Specifically, weconsider the case where the top level model M L is themodel of interest, i.e. q = E [ X L ] . The lower level models M , . . . , M L − are approximations of the top level modelthat are faster to evaluate but have a bias, i.e. E [ X l By comparing the expressions for the variance of the con-ventional and multilevel MC approaches, we can investigatethe potential speedup resulting from the MLMC approach.Let us consider the times ˜ t MC and ˜ t ML required to convergeto a given variance ˜ v = σ Q MC (˜ t MC ) = σ ∗ Q ML (˜ t ML ) . Then,combining (5) and (14) results in the expressionspeedup = ˜ t MC ˜ t ML = σ X √ τ P Ll =0 σ Y l √ τ l ! . (15) In practice, the variance of the lowest level is similar to that ofthe direct MC simulator, σ Y ≈ σ X , and the cost of evaluatingthe highest level pair is at least that of a direct evaluation of thehighest level, i.e. τ L ≥ τ . Considerable speedups are possibleif σ Y l √ τ l ≪ σ X √ τ for all l . Intuitively, this occurs when eachsimplified model M l − is much faster than the next level M l ,but returns very similar results for the majority of samples.Examples where this occurs naturally in the context of powersystem adequacy assessment will be discussed in Sections IVand V.In order to compare the compuational efficiency of var-ious implementations, we require an operational definitionof ‘computational speed’. Monte Carlo simulations are oftenrun with the goal to estimate the quantity q with a certainrelative accuracy, expressed using the coefficient of variation c q = σ Q /q . We note that both (7) and (14) can be broughtinto the form c q |{z} computational‘distance’ = z q |{z} speed × t |{z} time . (16)This implicitly defines the computation speed z q as z q := q tσ Q ( t ) . (17)This definition may be compared with the ‘figure of merit’used in [8]. The inclusion of the quantity q in (17) hasa number of advantages, provided that q = 0 . First, thespeed has dimensions / time , independent of the measure q .Second, speeds corresponding to different metrics are directlycomparable. For example, when z LOLE < z EENS , this indicatesthat the LOLE estimator is the limiting factor in achievingconvergence to a given coefficient of variation. And finally, thespeed metric and the implied computational distance are easilyinterpretable in terms of simulation outcomes. For example,in order to achieve a coefficient of variation of 1% (i.e. a‘distance’ 10,000) using a speed of s − , a simulation runof s is required.In the course of a simulation run, (17) can be used toestimate the computational speed, replacing q and σ Q bytheir empirical estimates. The speed z q for MC and MLMCestimation follow from (5) and (10) as z q,MC = ˆ q MC t MC ˆ σ X /n , (18) ˆ σ X = P ni =1 ( x ( i ) − n P nj =1 x ( j ) ) n − , (19)and z q,ML = ˆ q ML t ML P l ˆ σ Y l /n l , (20) ˆ σ Y l = P n l i =1 ( y ( i ) l − n l P n l j =1 y ( j ) l ) n l − . (21)II. C ONSIDERATIONS FOR IMPLEMENTATION A. Joint sample spaces The core of the MLMC algorithm is the joint generation ofsample pairs ( X ( l,i ) l , X ( l,i ) l − ) , used in (9b), in such a way thatthey are maximally correlated. The random variables X l and X l − have sample spaces Ω l and Ω l − , respectively, whichmust be combined into a joint sample space Ω ′ l . We highlighttwo common model patterns that naturally achieve this. 1) Pattern 1: component subsets: One common occurrencein system adequacy studies is that the lower level model M l − omits components that are present in the higher level model M l . As a result, the sample space Ω l can be written as aCartesian product Ω l = Ω l − × A l , (22)where A l is the sample space of components present in M l but not in M l − . We may then identify Ω ′ l and Ω l . In practicalterms this means that samples can be generated at the higherlevel l and unused elements are discarded for the simplermodels M l − . An example of this design pattern is exploredin Section IV. 2) Pattern 2: identical randomness: It is also easy toconceive of scenarios where M l and M l − have identicalsample spaces, so that Ω ′ l = Ω l = Ω l − . (23)This occurs when both models are driven by the same setof random inputs, but the higher level model performs morecomplex processing. An example of this model pattern is givenin Section V. B. Direct evaluation of expectations Occasionally, the base model M is sufficiently simple topermit direct computation of r = E [ X ] , either analyticallyor using a numerical approximation procedure. In those cases,the long run efficiency is enhanced by evaluating r directlyinstead of using its MC estimate. The standard deviation σ Y isthen equal to 0, or a value commensurate with the accuracy ofthe numerical approximation of r . Although direct evaluationof the lowest level is nearly always preferred, there maybe cases where the evaluation of E [ X ] is a comparativelytime-consuming operation and the optimal trade-off is morecomplex. In the examples that follow in Sections IV andV, direct evaluation is always possible, and results in fasterconvergence of the overall MLMC estimator.The use of an analytical result at the lowest level alsohighlights a connection between the MLMC method and thecontrol variate approach [5]. The control variate similarlymakes use of a simplified model for which an explicit solutioncan be calculated. It can therefore be considered as a specialcase of a bilevel MLMC procedure where the value E [ X ] isknown and the output X is scaled for optimal convergence.The control variate approach was used in [2] to speed upcomposite system adequacy assessment - a problem that isalso addressed in Section IV. C. Implementation Simulations were implemented in Python 3.7 and were runon an Intel i5-7360U CPU under macOS 10.14.6. A genericmultilevel sampler was developed with specialisations forparticular simulation studies. No effort was made to optimisethe execution speed of individual models, because the aim ofthis paper is not to maximise execution speed per se, but toinvestigate the relative speed between sampling strategies. Thecode used to generate the results in this paper is available [9].All MLMC simulations started with an exploratory run inwhich a sample with fixed size n (0) is taken at each level set Y l , in order to determine initial estimates of the evaluation cost ˆ τ l and variance ˆ σ Y l . This initial run is followed by a sequenceof follow-up runs, each parameterised by a target run time t ∗ .Given t ∗ , optimal sample sizes at each level were determinedusing (13) and the most up to date estimates of evaluationtimes ˆ τ l and variances ˆ σ Y l . For all results in this paper, 10runs with an estimated run time of 60 seconds (each) wereused, for a total run time of approximately 600 seconds.One practical concern with determining optimal samplesizes using (13) is that the values of σ Y l are estimated usingrelatively small data sets. In power system risk assessment,the simulation outputs X l often involve measurements of rareevents, so that there is a high probability that Y ( i ) l = 0 , andtherefore ˆ σ Y l ≪ σ Y l (or even ˆ σ Y l = 0 ). If the estimated valueis used naively in (13), this leads to undersampling of Y l ,thereby exacerbating the problem because fewer samples aregenerated that can correct the estimate of σ Y l . To mitigate thisrisk, the variance estimators were adjusted as follows. First, aconservative estimate for the variance of X was obtained as ˜ σ X = max l (ˆ σ X l ) . (24)Next, we assumed for the lowest level estimator that ˜ σ Y ≈ ˜ σ X , and that the ratio Y l +1 /Y l of variances of subsequentlevel contributions is lower-bounded by a factor α . Therefore,updated variance estimates are computed as ˜ σ Y l = max(ˆ σ Y l , α l ˜ σ X ) , (25)for those pairs l where E [ Y l ] is estimated by sampling. Forthe simulations, the value of α was heuristically set to 0.1.Finally, in simulations, multiple risk measures q ( a ) , q ( b ) , . . . were estimated in parallel. In determining optimal samplesizes, one of these was selected as the ‘target measure’ tooptimise for, so that its mean and variance estimates wereinserted in (13) to determine the optimal allocation of samplecounts n l .IV. C OMPOSITE SYSTEM ADEQUACY ASSESSMENT The first case study is a system adequacy assessment of thesingle area IEEE Reliability Test System (RTS) [10]. A two-level MLMC approach is used, where the upper level, i.e. thestudy of interest, is a hierarchical level 2 (HL2) study [1]: acomposite system adequacy assessment that takes into accounttransmission line outages and constraints. The lower levelHL1 is a single node assessment that omits the transmissionystem. This is in accordance with the subset model patternin Section III-A1. A. Models1) Model M - Composite system adequacy assessment(HL2): The RTS model defines outage probabilities of gen-erators and transmission lines, which were modelled as inde-pendent two state Markov models. Maintenance and transientoutages were not considered. Load levels were sampled byuniformly selecting an hour from the annual demand trace andassigning loads to each node in proportion to the maximumnodal demands.Therefore, at the upper level ( l = 1 ), a sampled systemstate ω ( i )1 consists of: (i) the nodal demand d ( i ) n for n ∈ N ,the set of nodes; (ii) the generator status γ ( i ) j ∈ { , } for j ∈ G , with G = S n ∈N G n , where G n is the set of generatorsin node n ; (iii) the line status λ ( i ) k ∈ { , } for k ∈ L , theset of transmission lines. Let generator and line flow limits begiven by g max j and f max k . Then, the amount of curtailment C is computed by the linear program C ( ω ( i )1 ) = min c |N| ,g |G| X n ∈N c n , (26)subject to ≤ c n ≤ d ( i ) n , ∀ n ∈ N ≤ g j ≤ γ ( i ) j g max j , ∀ j ∈ G− f max k ≤ X n ∈N M ( i ) kn [ X j ∈G n g j + c n − d ( i ) n ] ≤ f max k , ∀ k ∈ L X n ∈N [ X j ∈G n g j + c n − d ( i ) n ] , where the matrix M ( i ) = D ( i ) A ( A T D ( i ) A +1 / |N | ) − relatesbus injections and line flows. The directed line-node incidencematrix A has elements +1 for outgoing lines and − forincoming lines; the diagonal matrix D ( i ) has elements D ( i ) kk = λ ( i ) k /x k , where x k is the reactance of line k . The element-wiseconstant / |N | ensures invertibility, eliminating the need fora designated slack bus. In cases where line outages resulted inmultiple islands, problem (26) was formulated and solved foreach island independently and the curtailments were summedto obtain the total system curtailment. Linear optimisationwas performed using scipy.optimize.linprog , withthe revised simplex method. 2) Model M - Generation adequacy assessment (HL1): For HL1 assessment, a single-node generation adequacy anal-ysis is performed, without transmission line constraints andoutages. The lower level system state ω ( i )0 can thus be obtainedfrom ω ( i )1 by omitting the line status variables. For this HL1study, the curtailment is calculated as C ( ω ( i )0 ) = max , X n ∈N d ( i ) n − X j ∈G n γ ( i ) j g max j . (27) TABLE IC OMPOSITE SYSTEM ADEQUACY ASSESSMENT - AVAILABLE MODELS model description z PLC [1/s] z EPNS [1/s] direct evaluation M HL2 . 31 0 . no M HL1 . ∗ . ∗ optional ∗ : inherent estimation bias 3) Risk measures: Two common risk measures were com-puted: the probability of load curtailment (PLC) and expectedpower not supplied (EPNS). The related performance measures X q,l are defined in terms of the load curtailment (27) and (26)as X PLC ,l ( ω ) = C l ( ω ) > , (28) X EPNS ,l ( ω ) = max(0 , C l ( ω )) . (29) B. Results Throughout, Monte Carlo estimates of risk measures aregiven with the relevant number of significant digits, fol-lowed by the estimated standard error in parentheses. Thus, . × − stands for an estimate of . with astandard error of . . For all MLMC runs, an initialexploratory run with n (0) = 100 was used, followed by 10runs of approximately 60 seconds. The target risk measure forsample size optimisation was EPNS. Unless stated otherwise,thermal line ratings were scaled to 80% of the nominal values,to tighten network constraints.Table I compares the speed obtained with the individualmodels for the estimation of both PLC and EPNS risk mea-sures, and whether direct evaluation of the expectation ispossible (for an effective ‘sampling speed’ of z = ∞ ). Thesenumbers were estimated at the end of 10-minute conventionalMC runs. The HL1 model is over one hundred times fasterthan the HL2 model for both measures, but of course thiscomes at the cost of an estimation bias.Table II compares the results of three different estimators.The top row is the conventional Monte Carlo estimator thatdirectly performs the HL2 study. The middle row representsa two-level MLMC approach where HL2 sampling is com-bined with HL1 sampling, immediate leading to significantspeedups of 2.5 (for PLC) and 10 (for EPNS). In the thirdconfiguration (bottom row), further speedups are obtained byeliminating sampling of the lower level model, and com-puting the lower level estimates r PLC , = E [ X PLC , ] and r EPNS , = E [ X EPNS , ] directly by convolution using 1 MWdiscretisation steps.An interesting observation is that, for the regular MC sam-pler, the speed of PLC estimation (0.31 s − ) is larger than thatof EPNS estimation (0.17 s − ). However, the MLMC samplersees much more substantial speedups for EPNS estimationthan for PLC estimation. This is only partially caused bythe EPNS-focused sample size optimisation. The other factoris that the discontinuous loss-of-load indicator (28) is lessamenable to successive approximation [5]. ABLE IIC OMPOSITE SYSTEM ADEQUACY ASSESSMENT - COMPARISON OF APPROACHES PLC estimation EPNS estimationestimator models used run time [s] PLC z PLC [1/s] speedup EPNS [MW] z EPNS [1/s] speedupMC M . × − M , M . × − . M , M . × − . OMPOSITE SYSTEM ADEQUACY ASSESSMENT - MULTILEVELCONTRIBUTIONS term PLC EPNS [MW] τ l [ms] n l ˆ r . × − . . 93 158 ˆ r . × − . . . × − . Table III gives insight into the multilevel structure ofthe regular MLMC estimate. For both PLC and EPNS, therefinement term ˆ r is substantially smaller than the crudeestimate ˆ r . More importantly, sampling from the HL1 modelis substantially faster (0.023 ms per evaluation) than the HL2-HL1 difference term (5.4 ms per evaluation), due to the linearprogram (26) involved in the latter. The MLMC algorithmadapts to this cost difference by invoking the HL1 modelnearly 50 times as often.Finally, Table IV shows the impact on convergence speedof varying the thermal line ratings between 80% and 100% ofthe nominal values. Higher line ratings cause fewer constraints,which results in a slight reduction in speed for the regular MCsampler. On the other hand, the MLMC sampler experiencesvery large speedups as the difference between the results fromthe HL1 and HL2 models becomes smaller, so that fewer(expensive) HL2 evaluations are required. Once again, thegains in EPNS estimation speed exceed the gains in PLCestimation speed. TABLE IVC OMPOSITE SYSTEM ADEQUACY ASSESSMENT - THERMAL RATINGS relative PLC estimation EPNS estimationline rating z MC z ML speedup z MC z ML speedup0.8 0.31 1.04 3.3 0.17 2.54 150.9 0.26 1.38 5.3 0.14 4.69 341.0 0.25 2.11 8.6 0.12 16.7 143 V. D ISPATCH OF STORAGE The second example concerns the assessment of systemadequacy in the presence of energy-constrained storage units(e.g. batteries). The energy constraints couple decisions insubsequent time slots, thus necessitating the use of time-sequential Monte Carlo simulations. Convergence for time-sequential simulations tends to be much slower than for snap-shot problems, due to significant correlations in visited systemstates. An additional complication is deciding an appropriate dispatch strategy for energy storage units. A greedy EENS-minimising discharging strategy was recently proposed in [11],as a reasonable default dispatch strategy for adequacy studies. A. Models The Great Britain (GB) adequacy study from [11] is re-produced here, with an eye on speeding up estimation ofloss of load expectation (LOLE) and expected energy notsupplied (EENS) risks using the MLMC approach. Individualsimulations are run for a sequence of 8760 hours (1 year). Thesystem performance in a simulated year is driven entirely bythe net generation margin trace M t ( ω ( i ) ) = g ( i ) t + w ( i ) t − d ( i ) t , t ∈ { , . . . , } , (30)where the sampled state ω ( i ) consists of the demand trace d ( i ) t ,wind power trace w ( i ) t and conventional generation trace g ( i ) t .Annual demand traces are chosen randomly from historicalGB demand measurements for 2006-2015 (net demand, [12]).Annual wind traces are similarly sampled from a syntheticdata set for hypothetical GB wind power output for the period1985-2014, derived from MERRA reanalysis data and anassumed constant distribution of wind generation sites withan installed capacity of 10 GW [13]. Conventional generationtraces are generated using an assumed diverse portfolio ofthermal units. The portfolio of 27 storage units was based onstorage units contracted in the GB 2018 T-4 capacity auction.The reader is referred to [11] for further model details.We consider four different storage dispatch models. Theresulting storage dispatch (with sign convention that con-sumption is positive) is denoted by S t,l ( ω ) , and is entirelydetermined by the net generation margin M t ( ω ( i ) ) . All fourmodels are defined on the same sample space Ω , providingan example of the model pattern described in Section III-A2.However, the models differ tremendously in computationalcomplexity, as is clear from the descriptions below. TABLE VT IME - SEQUENTIAL SIMULATION WITH STORAGE - AVAILABLE MODELS directmodel description z LOLE [1/s] z EENS [1/s] evaluation M EENS-optimal . 105 0 . no M sequential . ∗ . ∗ no M average . ∗ . ∗ optional M ′ no storage . ∗ . ∗ optional ∗ : inherent estimation biasABLE VIT IME - SEQUENTIAL SIMULATION WITH STORAGE - MODEL COMPARISON run LOLE estimation EENS estimationestimator models used time [s] LOLE [h/y] z LOLE [1/s] speedup EENS [MWh/y] z EENS [1/s] speedupMC M . M , M , M ′ . M , M . M , M , M . 112 2 113 1) Model M - EENS-optimal dispatch: The storage dis-patch S t, ( ω ) is computed using the EENS-minimising algo-rithm given in [11]. It is sequential and requires complex logicfor each step. 2) Model M - Sequential greedy dispatch: The storagedispatch S t, ( ω ) is computed using a heuristic approximationof the EENS-minimising policy. Storage units s ∈ S are sortedby decreasing time to go (from full) e s /p s , where e s and p s are energy and discharge power ratings, respectively. Then,a sequential greedy dispatch is performed, charging whenpossible, and discharging only when required to avoid loadcurtailment. Evaluating this model requires one sequential passper storage unit, but the simulation steps are trivial. 3) Model M - Constant peak-shaving dispatch: The stor-age fleet is optimistically approximated by a single storageunit with e = P s e s and p = P s p s . A mean dailydemand profile ˜ d is computed by averaging demand overall historical days. This profile is used to compute a single daily dispatch pattern ˜ s that solves to following quadraticoptimisation problem to flatten the average total demandprofile ( ˜ d h + ˜ s h ) h =1:24 : ˜ s = arg min s ,e X h =1 ( ˜ d h + s h ) , (31)subject to − p ≤ s h ≤ p, h = 1 , . . . , ≤ e h ≤ e, h = 1 , . . . , e h +1 = e h + s h × , h = 1 , . . . , e = e + s × . This problem was solved using the Python quadprog pack-age. The resulting annual storage dispatch is obtained byrepeating the 24-hour dispatch pattern: S t, ( ω ) = ˜ s ( t mod 24) . (32)Because S t, is a deterministic load offset, risk measures forthis model can be computed by convolution. 4) Model M ′ - No storage: This alternative lowest levelmodel does not use storage at all, so that S t, = 0 . 5) Risk measures: The net generation margin M t ( ω ) andstorage dispatch S t,l ( ω ) result in a curtailment trace as follows C t,l ( ω ) = max[0 , − M t ( ω ) + S t,l ( ω )] , ∀ t. (33) The LOLE and EENS risk measures can be computed usingthe performance measures X LOLE ,l ( ω ) = X t =1 C t,l ( ω ) > , (34) X EENS ,l ( ω ) = X t =1 C t,l ( ω ) × h. (35) B. Results Table V compares the speed obtained with the individualmodels for the estimation of both LOLE and EENS riskmeasures, and whether direct evaluation of the expectation ispossible with each model. All numbers were estimated at theend of 10-minute conventional MC runs. Very large differencesin model speed are visible, with the detailed model M beingover 500 times slower than the crude model M ′ .For MLMC simulations, an exploratory run with n (0) = 20 was used, followed by 10 runs of 60 seconds, where samplesizes were optimised for the EENS risk measure. In all cases,the crude estimate r = E [ Y ] was evaluated using a convo-lution approach. Results are shown in Table VI, comparingthe performance of three MLMC architectures with directMC simulation. A three-layer architecture using model M ′ (without storage) as a bottom layer achieved speedups of 10(LOLE) and 30 (EENS), but much better results were obtainedwhen the daily average dispatch model M was used - evenwhen a two-layer MLMC stack was created by omitting theintermediate sequential greedy dispatch model.The results show that the MLMC performance is very sen-sitive to the choice of levels, but robust speedups are availableeven for sub-optimal model choices. The best performingarchitecture is further analysed in Table VII. It can be seenthat the contribution from the final refinement ˆ r is minimal,i.e. the heuristic model is very accurate, which is key to theobserved speedup of . The MLMC algorithm dynamically TABLE VIIT IME - SEQUENTIAL SIMULATION WITH STORAGE - MULTILEVELCONTRIBUTIONS term LOLE [h/y] EENS [MWh/y] τ l [ms] n l ˆ r − . ˆ r − . − ˆ r n/a n/asum . djusted sample sizes to generate more samples evaluating Y = X − X than on the costly evaluation of Y = X − X .Moreover, no samples are spent on the contribution ˆ r , whichcan be computed directly by convolution. As a result, the speed z EENS is able to exceed even that of the fastest model inTable V (for regular MC estimation).VI. C ONCLUSIONS AND FUTURE WORK This paper has set out how the MLMC approach canbe applied to power system risk analysis, and specificallyto system adequacy assessment problems. Common modelpatterns were identified that are particularly amenable toMLMC implementation, and a computational speed measure(17) was introduced to quantify simulation speed in a waythat is easily comparable across tools, Monte Carlo methodsand risk measures. Two case studies illustrate the potential forspeeding up estimation of risk measures, and the ability toapply the method to complex simulations.In future work, we will consider automatic selection ofoptimal model stacks, and explore the scope for the applicationof multi-index Monte Carlo schemes [14] .A CKNOWLEDGMENTS We thank Kate Ward and Iain Staffell for the provision ofGB demand and wind power data, and Michael Evans forsharing Python code for the EENS-minimising dispatch. Weare also grateful for insightful discussions with Chris Dent,Matthias Troffaes, Amy Wilson and Stan Zachary.R EFERENCES[1] R. Billinton and W. Li, Reliability Assessment of Electric Power SystemsUsing Monte Carlo Methods . Boston, MA: Springer US, 1994.[2] R. Billinton and A. Jonnavithula, “Composite system adequacy assess-ment using sequential Monte Carlo simulation with variance reductiontechniques,” IEE Proceedings - Generation, Transmission and Distribu-tion , vol. 144, no. 1, p. 1, 1997.[3] F. Belmudes, D. Ernst, and L. Wehenkel, “Cross-Entropy Based Rare-Event Simulation for the Identification of Dangerous Events in PowerSystems,” Probabilistic Methods Applied to Power Systems, 2008.PMAPS ’08. Proceedings of the 10th International Conference on , pp.1–7, 2008.[4] A. Leite da Silva, R. Fernandez, and C. Singh, “Generating CapacityReliability Evaluation Based on Monte Carlo Simulation and Cross-Entropy Methods,” IEEE Transactions on Power Systems , vol. 25, no. 1,pp. 129–137, feb 2010.[5] M. B. Giles, “Multilevel Monte Carlo methods,” Acta Numerica , vol. 24,pp. 259–328, may 2015.[6] L. J. Aslett, T. Nagapetyan, and S. J. Vollmer, “Multilevel Monte Carlofor Reliability Theory,” Reliability Engineering & System Safety , vol.165, pp. 188–196, sep 2017.[7] A. Huda and R. ˇZivanovi´c, “Improving distribution system reliabilitycalculation efficiency using multilevel Monte Carlo method,” Interna-tional Transactions on Electrical Energy Systems , vol. 27, no. 7, p.e2333, jul 2017.[8] P. Henneaux and P.-E. Labeau, “Improving Monte Carlo simulationefficiency of level-I blackout probabilistic risk assessment,” in , Durham, 2014.[9] S. Tindemans, “Code release for ‘Accelerating System Adequacy As-sessment using the Multilevel Monte Carlo Approach’,” Apr. 2020. DOI:10.5281/zenodo.3757332 [10] C. Grigg, P. Wong, P. Albrecht, R. Allan, M. Bhavaraju, R. Billinton,Q. Chen, C. Fong, S. Haddad, S. Kuruganty, W. Li, R. Mukerji,D. Patton, N. Rau, D. Reppen, A. Schneider, M. Shahidehpour, andC. Singh, “The IEEE Reliability Test System-1996. A report prepared bythe Reliability Test System Task Force of the Application of ProbabilityMethods Subcommittee,” IEEE Transactions on Power Systems , vol. 14,no. 3, pp. 1010–1020, 1999.[11] M. P. Evans, S. H. Tindemans, and D. Angeli, “Minimizing UnservedEnergy Using Heterogeneous Storage Units,” IEEE Transactions onPower Systems , vol. 34, no. 5, pp. 3647–3656, sep 2019.[12] I. Staffell and S. Pfenninger, “The increasing impact of weather onelectricity supply and demand,” Energy , vol. 145, pp. 65–78, feb 2018.[13] ——, “Using bias-corrected reanalysis to simulate current and futurewind power output,” Energy , vol. 114, pp. 1224–1239, nov 2016.[14] A.-L. Haji-Ali, F. Nobile, and R. Tempone, “Multi-index Monte Carlo:when sparsity meets sampling,”