Enhancing Distribution Resilience with Mobile Energy Storage: A Progressive Hedging Approach

Enhancing Distrib ution Resilience with Mobile Ener gy Storage: A Progressi ve Hedging Approach Jip Kim, Student Member , IEEE, and Y ury Dvorkin, Member , IEEE, Email: { jipkim, dvorkin } @nyu.ed u Abstract —Electrochemical energy storage (ES) u nits (e.g. bat- teries) hav e been field-validated as an effi cient back-up resource that enhance resilience of the distribution system in case of natural disasters. Howe ver , u sing these units fo r resilience is not sufficient to economically justify their installation and, theref ore, these units a re often installed in locations where they incur the greate st economic v alue during normal opera tions. Motivated by the recent progress in transportable ES technologies, i.e. ES units can be mo ved u sing public transportation ro utes, this paper p roposes to use this spatial flexibility to brid ge the gap between th e economically optimal locations during n ormal operations and disaster -specific locations where extra back-up capacity is necessary . W e propose a two-stage optimization model that optimizes inv estments in mobile ES un its in the first stage and can re-route the in st all ed mobile ES un its in th e second stage to a void the expected l oad sheddi ng caused by disaster fo recasts. Since the proposed model cann ot be solved efficiently with off-the-sh elf solvers, even for relativ ely small instances, we apply a pr ogre ssiv e hedging a lgorithm. The pro posed m odel and progr essive hedging algorithm ar e tested thro ugh two illustrativ e examples on a 15-bus radial distributi on test system. Index T erms —Energy storage, resilience, progr essive hedging I . I N T R O D U C T I O N The US p ower grid is vuln erable to n atural disasters (e.g. flooding , extreme winds, earthqu akes) as they increase the likelihood of critical equipmen t failures [1]. Distribution sys- tems are p articularly affected b y natural disasters due to the compou nding effect of line o utages, radial topolo gy , and lim- ited back -up resources. Furth ermore, large- scale d isasters can affect the power grid in mu ltiple locations cau sing th e domino effect, which may in turn spread power ou tages acro ss large geogra p hical are a s, e ven if some areas would no t otherwise be a ffected by disasters [ 2]. Th erefore , it is impo rtant to contain th e failures within distribution systems a n d p revent their fu rther pro p agation in the transmission system. Provided a time ly d isaster forecast, power gr id operators can pr ev ent these failure s, or at le a st mitigate their imp acts, by strategically pla c ing flexible back -up resou r ces in th e distribution system to strengthe n its weak n odes (buses) an d weak lin ks (lines) [3]. Since the disaster f o recasts ar e no rmally av ailable on a short notice (fo r instanc e , NO AA forec asts are av ailable 4 8-16 8 hou rs ahead [4]), o nly a f ew tech nologies are ph ysically suitable to be de p loyed, relocated or ac tivated within this time fr ame: distributed en ergy reso urces, includ ing portable diesel gen e r ators [5] a nd parked electric vehicles [3], adaptive microg rids [6], or preventi ve load she dding [7]. The o bjective of this p aper is to illustra te that mobile ES The authors are with the Department of Electrical and Computer Engineer- ing, T andon School of Engineering, New Y ork Uni versi ty . This work was supported in part by the US NSF Grant No. ECCS-1760540. units, i.e. those that can be m oved among d ifferent distribution system locations u sing regular transp ortation ro utes, are also econom ic a lly and physically suitable for this application . A mob ile ES u n it, of te n referred to as ‘storag e-on–wh eels’, is an emerging technolo gy that has been re cently developed in the fo rm of a trailer-mounted electroch emical b attery . Consol- idated E d ison of New Y ork is cur rently c onsidering installing such mobile E S units to redu ce th e impac t o f PV g eneration on th eir distribution system in New Y o rk City a n d d efer costly distribution u pgrade s [8]. As r eported in [8], each mo bile ES unit will featu r e a Lithium -ion battery that can sto r e up to 800 kWh of energy with the maximu m charging/d isch arging duration of 1 0 hou rs. Relativ e to o ther resilienc e reso u rces, the mobile ES un its have multiple advantages. Fir st, th ey are more en vironmen tally friend ly than p ortable diesel gen erators and can be used without un necessary noise and air pollution . Furthermo re, strict pollutio n standa rds in certain jurisdic tio ns, especially metr opolitan areas, pr o hibit the use of portab le diesel g enerators d uring no rmal op erations, which redu ce their value to the distribution system. Second , u nlike electric vehicles or adaptive m icrogrid s, the mobile E S units can be d irectly op erated by the power gr id oper a to r and do not require advanced commu nication in frastructu r e or engage m ent with electr icity co nsumers. Finally , rec ent ad vances in mater ia l sciences ar e expected to enable mobile ES units with just-in- time re-c o nfigura b le power and energy r atings (so - called flow batteries) tha t can f u rther enh ance the back -up flexibility th a t these units can p r ovide in case of natural disasters. Since mobile ES un its ca n b e used during no rmal op e r ations and f or d ealing with the effects of natural disasters, there is a need to integrate their mo bility in plan ning too ls for power g rid o perator s. A typ ical state-of- the-art plannin g tool for stationa ry ES un its, [9], [10], is rou tinely formu lated as a two-stage stochastic mixed-integer linear program (MILP) and co nsiders ES units a s stationar y resou rces. In th ese to ols, the first stage optimizes the ES locatio n s an d sizes, wh ile the second stage fixes the first-stage decisions and co-o ptimizes the operatio n o f existing resource s an d newly in stalled ES units. T o immu nize the first-stage inv estment d ecisions ag ainst a variety o f oper ating cond itions, the second stage n ormally considers multiple op e rating scen arios. If the ES units were mobile, the c o mplexity of the plann in g tool would incr ease. First, the E S mobility imp lies that the ES location is not fixed in the second stage and needs to be o p timized (so- called recourse decision s). Sec o nd, the recou rse d ecision on each mo bile ES u nit is binary , wher e it attains the value of 1, if that unit n e eds to be moved, or 0 o therwise. The two-stage stochastic MILPs with binary recou rse decision s are particularly comp utationally deman d ing, and o ften exist- ing solu tion strategies, e.g. Bende r s’ deco mposition, perfor m poorly when applied to such problems [11]. This compu ta - tional com plexity can b e overcome by using the pro gressive hedging (PH) algorithm [12]. Recently , the PH algorithm ha s gained atten tio n in the co ntext of a two-stage stochastic un it commitmen t proble m [13]. Unlike Benders’ decom position, the PH alg orithm do es n ot exploit the two-stage stru c tu re of the u nderlyin g optimization and d oes no t separate th e first- and second-stage decisions. Instead , the PH algo rithm partition s the or ig inal prob lem in a scenario-b ased fashion an d the first- and secon d-stage decisions are op timized for eac h scenario indepen d ently . In each scen ario-ba sed pro blem, the relax ation of the first-stage decision is p enalized with an exogen ous penalty coefficient. The algorith m iterates u ntil th e first-stag e decisions across all scen arios conver ge with a given tolerance . The scenario -based decom position u sed in the PH algor ithm is shown to be an effecti ve solu tio n strategy for two-stage MI LPs with b inary recourse d ecisions and, therefore, it is applica b le for p lanning too ls with mobile ES units. The pr imary focu s of th is pro blem is two-f old: 1) It ta kes the perspective of the distribution system ope r- ator and fo rmulates the optim ization prob lem to decide on the inv estments in mobile ES units. The pro posed optimization is a two-stage mixed-integer progr am with binary r e c ourse d ecisions, wh ich acco unt for th e reloca- tion of mobile ES units under each specific scenario. This o ptimization ach iev es the trade-o ff between the econom ic value of mobile ES units d uring normal opera- tions and their ab ility to en hance the d istribution system resilience in case of natur al disasters. 2) The proposed optimization is solved using the PH al- gorithm. The n umerical results demo nstrate that this method o u tperfor ms off-the-shelf solvers in terms of the co mputation al perform a nce a n d solution accuracy . The nume r ical experiments also su g gest that the PH perfor mance can be improved by tun in g extern alities (e.g. p enalty coefficients). The rem ainder of th is pape r is o rganized as follows. Sec- tion II and Sectio n III describ e th e propo sed op timization and the PH imp lem entation. Section IV d escribes the case stud y and dem onstrates the u sef ulness of mo bile ES u nits. Section V conclud e s the paper . I I . P L A N N I N G W I T H M O B I L E E S U N I T S In the following, we co n sider a distribution system la y out typical to the US power sector . The system has a radial topolog y , where sets B an d L , indexed b y b and l , rep r esent the distribution buses ( nodes) and lines (edges). The single distribution system o p erator is responsib le for all planning an d operating decisio n s, as well as incu rred costs. Th e d istribution power flows a r e m odeled using a second-o rder conic ( SOC) relaxation of ac power flows that is comp utationally tr a c table, [14]. In this case the objective of th e system operato r is to allocate th e m obile ES units with fixed power and en ergy ratings (the set of units is denoted by K , indexed by k ) in such a way th at these u nits are op e r ated as stationar y resour ces during th e nor mal operatio ns and can be moved to a different location in case of natur al disasters. Set S , ind exed by s , contains scena r ios of the normal op e r ations and emergency condition s. Each scen ario is a ssign ed prob a bility ω s and h as the num ber of time intervals denoted by set T , indexed by t . The o bjective fu nction of th is o ptimization is then g iv en as: min  γ · I C ( x k ) + X s ∈S ω s · O C s ( u kbts | x k )  (1) I C ( x k ) = X k ∈K  C P · P k + C E · E k  · x k (2) OC s ( u kbts | x k ) = X t ∈T ,b ∈B C g b · p g bts + X k ∈K ,b ∈B ,t ∈T C V oLL b · p ls bts + X k ∈K ,b ∈B ,t ∈T      h 100     C P  p ch kbts + p dis kbts   , ∀ s ∈ S (3) Eq. (1) minimizes the sum o f th e in vestment cost I C ( · ) , where γ is a daily capital recovery factor from [9] that prorates the in vestmen t cost on a d a ily b asis, and the expected expected daily o perating co st over all scenarios OC s ( · ) . Eq. (2) c o m- putes the inv estment cost o f in stalling mo b ile ES unit k with fixed power an d energy r atings P k and E k , eac h priced with parameters C P and C E as in [ 9], [10]. The installation d ecision is modelled by bina r y variable x k ∈  0 , 1  . If x k = 1 , mob ile ES u nit k is in stalled, othe r wise x k = 0 . Note that P k and E k are par ameters, but can be extended to decision variables. Eq. (3 ) c omputes the op erating cost for each scen ario, wh ich includes: operatin g cost of gen eration resour ces with the incr e- mental co st C g b produ cing p g bts , lo ad shed ding co st based o n p ls bts and value of lo st load C V oLL b , an d the ES degradatio n co st. The degradation cost is computed as explained in [15] and depend s on c harging ( p ch kbts ) and discharging ( p dis kbts ) decision variables and p arameter h , the degradation slop e taken from [15]. As descr ibed belo w , OC s ( u kbts | x k ) depends on the binary d ecision u kbts ∈  0 , 1  , w h ich is a re c o urse decision on the placeme n t of mobile ES unit k at bus b at time in te r val t under scenario s a nd that in turn dep ends on the scenario- and time-indep endent decision x k . The b in ary recou rse is deno te d in (1) as u kbts | x k and is im plicitly inter n alized in the rig ht- hand side of (3) via variables p ch kbts and p dis kbts . The mo bile ES units are oper a te d as ( ∀ k ∈ K , t ∈ T , s ∈ S ): e kts = e k,t − 1 ,s + X b ∈B  p ch kbts · ℵ ch − p dis kbts / ℵ dis  , (4) 0 ≤ e kts ≤ E k , (5) 0 ≤ p ch kbts · ℵ ch ≤ P k · u kbts , ∀ b ∈ B (6) 0 ≤ p dis kbts / ℵ dis ≤ P k · u kbts , ∀ b ∈ B (7) − K · p ch kbts ≤ q ch kbts ≤ K · p ch kbts , ∀ b ∈ B (8) − K · p dis kbts ≤ q dis kbts ≤ K · p dis kbts , ∀ b ∈ B (9) Eq. (4) relates the en ergy state-of-charge e kts and ch a rg- ing/discharging decisions with imperfect efficiency ℵ ch = ℵ dis < 1 . Eq . (5) limits the en ergy stored to E k . E q. (6)- (7) limit the ma x imum ch arging and discharging power to P k and binary variable u kbts ∈ { 0 , 1 } in d icates if mob ile ES u nit k is located a t bus b in time period t and scenario s . If u kbts = 0 , the energy state-of-ch arge remains unchan ged. Eq. (8) and (9) relate the reactive power injec tion of m obile ES units to th eir charging a nd discharging power via para m eter K gi ven by a desirable p ower factor . No te that (6)-(7) co ntain prod ucts of continuo us and b inary variables, which can be linearize d u sin g the big-M method ( see [1 0] for details). The distribution system is operated as ( ∀ t ∈ T , s ∈ S ): ( f p lts ) 2 + ( f q lts ) 2 ≤ ((1 − α lts ) S l ) 2 , ∀ l ∈ L (10) ( f p lts − a lts · R l ) 2 + ( f q lts − a lts · X l ) 2 ≤ (( 1 − α lts ) S l ) 2 , ∀ l ∈ L (11) v s ( l ) ,t,s − 2( R l · f p lts + X l · f q lts ) + a lts  R 2 l + X 2 l  = v r ( l ) ,t,s , ∀ l ∈ L , (12) ( f p lts ) 2 + ( f q lts ) 2 a lts ≤ v s ( l ) ,t,s , ∀ l ∈ L (13) − X l | r ( l )=0 ( f p lts − a lts · R l ) − p g 0 ,t,s + G 0 · v 0 ,t,s = 0 , (14) − X l | r ( l )=0 ( f q lts − a lts · X l ) − q g 0 ,t,s − B 0 · v 0 ,t,s = 0 , (15) f p bts − X l | r ( l )= b ( f p lts − a lts · R l ) − p g bts + P d bts − p l s bts + G b · v bts − X k ∈K p dis kbts + X k ∈K p ch kbts = 0 , ∀ b ∈ B , (16) f q bts − X l | r ( l )= b ( f q lts − a lts · X l ) − q g bts + Q d bts − q l s bts − B b · v bts − X k ∈K q dis kbts + X k ∈K q ch kbts = 0 , ∀ b ∈ B , (17) P g b ≤ p g bts ≤ P g b , ∀ b ∈ B G , (18) Q g b ≤ q g bts ≤ Q g b , ∀ b ∈ B G , (19) 0 ≤ p ls bts ≤ P d bts , ∀ b ∈ B , (20) 0 ≤ q ls bts ≤ Q d bts , ∀ b ∈ B , (21) V b ≤ v bts ≤ V b , ∀ b ∈ B , (22) Eq. (10)-( 13) are the SOC-based p ower flow mo del as in [14], where active and reactive power flows are f p lts and f q lts and squared magnitu de o f nod al voltages at sen ding and receiving buses of ea c h line l are v s ( l ) and v r ( l ) . The apparen t flow limit, r esistance an d reactance of each line are given by S l , R l and X l . Scenario -depen dent parameter 0 ≤ α lts ≤ 1 is used to emulate the imp act o f natural disasters on the distribution system lin es. If α lts = 1 , line l is tripped . On the other han d, if 0 < α lts < 1 , line l is operated with a reduced apparent flow limit. Finally , if α lts = 0 , line l is operated with the nor m al appar ent flow limit. Eq. (14)-(15) enforce the activ e and re activ e nodal power balance for the root bus (de n oted with index 0 ) of the distribution system. Th e nodal power balan ce for o ther distribution buses is enforc e d in (1 6) an d (17), where P d bts and Q d bts are the active and reactive power d emand. Active and r eactiv e power injection s of conventional g enerators are co nstrained in (18) an d (19) using their minimum a n d maximu m limits ( P g b , Q g b , P g b , Q g b ). The nodal lo ad she d ding is limited to the nod a l dem and in (20) and (2 1). The nodal voltages ar e to be kept with the u pper ( V b ) and lower ( V b ) limits enfor ced in (22). The mo bility of ES u n its is modeled as: X b ∈B u kbts ≤ x k , ∀ k ∈ K , t ∈ T , s ∈ S (23) X k ∈K u kbts ≤ N ES b , ∀ b ∈ B , t ∈ T , s ∈ S (24) u kbt,s 1 = u kb,t 1 ,s 1 , ∀ s ∈ S , k ∈ K , b ∈ B , t ∈ T \  t 1  (25) u kb,t 0 ,s = u kb,t 0 ,s 1 , ∀ k ∈ K , b ∈ B , s ∈ S (26) u kb 1 ts − u kb 1 ,t +1 ,s ≤ 1 − u kb 2 ,t + τ ,s , ∀ s ∈ S \  s 1  , k ∈ K , b 1 6 = b 2 ∈ B , t ∈ T , τ ∈  1 , · · · , min( T d b 1 ,b 2 ,t , N t − t )  (27) Eq. ( 23) relates u kbts to the inv e stme nt decision x k accounted for in (1). If x k = 0 , it follows fr o m ( 23) u kbts = 0 . If x k = 1 , u kbts is optimized for each scen ario. Eq. (24) limits the numb er of m o bile ES units ( N ES b ) that can be co n nected to bus b . Eq. (25) ensures that mobile ES u nits are operated as station a ry reso urces du ring the no rmal o perations (d enoted as scena r io s = 1 ). Eq. (2 7) m odels the transition delay on m oving mob ile ES unit fro m bus b 1 to bus b 2 , where T d b 1 ,b 2 ,t is a gi ven transition time between buses b 1 and b 2 . In practice, the value o f T d b 1 ,b 2 ,t can be d etermined b ased on the av ailability an d le n gth of transpo r tation routes. Note that (27) is somewhat equivalent of minimum up and down time constrain ts on the on /off status of conventional generators modeled in UC problems, see [13]. I I I . S O L U T I O N T E C H N I Q U E The optimization prob le m in (1)-(2 7) is a two-stage stoch as- tic mixed-integer progr a m with binary recourse decisions. In general, such prob le m s can be solved usin g o ff-the-shelf solvers, but their perfor mance is limited , especially fo r large networks. T o s olve this problem efficiently , we apply th e PH alg o rithm. Following [ 12], we decom pose th e o riginal problem in (1)-(2 7) in N s = car d ( S ) subpr oblems, wh ere ea c h subprob lem can be solved indepe n dently in a pa rallel fashion. Due to the pagination limit, we cond e nse the trad itional 9 steps of the PH algorithm, [12], in three steps: • Step 1 : The PH algo r ithm is initialized by setting the iteration coun ter i = 0 an d PH multiplier m ( i =0) s = 0 . • Step 2: Each of N s subprob lems is solved in parallel to obtain binary decisions u ( i =0) kbts . When all sub problem s are solved, we com p ute u as the weighted average of all subprob lem solutio ns. • Step 3: T he iteration co unter is upd ated. For each sub- problem we up date the value of PH mu ltiplier m ( i ) s using the optimal solutio n of the previous iteration u ( i − 1) kbts , the value of u and exo genous penalty c o efficient ρ . Then we resolve each subproblem , where the deviation of u ( i ) kbts from u is p e nalized. After the sub problem s are solved, the value o f u is r ecompu ted using u ( i ) kbts . T h e iterative process con tinues until the mismatch g ( i ) is b elow a given tolerance ε . T o avoid issues with convergence, we f ollow the r e c ommend a- tion o f [12] and set the value of ρ based on the cost co e fficients of respe ctiv e hedged variables. Algorithm 1: PHA for Installing M obile ES units Step 1 . i := 0 , m ( i =0) s := 0 Step 2 . for s ← 1 to N s do u ( i =0) kbts ← a rg min u γ · I C + OC s end u = P s ∈S ω s · u ( i =0) k,b,t =1 ,s Step 3 . do i ← i + 1 for s ← 1 to N s do m ( i ) s ← m ( i − 1) s + ρ ·  u ( i − 1) k,b,t =1 ,s − u  u ( i ) kbts ← arg min u γ · I C + OC s + m ( i ) s · u k,b,t =1 ,s + ρ 2 k u k,b,t =1 ,s − u k 2 end u = P s ∈S ω s · u ( i ) k,b,t =1 ,s g ( i ) ← P s ∈S ω s ·    u ( i ) k,b,t =1 ,s − u    while convergence : g ( i ) < ǫ return u ( i ) kbts I V . C A S E S T U DY The case study is performe d based on the 15-bus radial distribution test system described in [16]. Th e distribution system o perator aims to install one mobile ES units with E k =1MWh, P k =0.15MW and ℵ ch = ℵ dis =0.9. Based on [9], the capital co sts ar e C P =$1,00 0/kW and C E =$50/kWh and the expec ted lifetime is 1 0 years. The disaster f orecasts are represented by th e scenarios d escribed in T ab le I. Case 1 co n- siders the d isaster repr esented by one scenario ( s 2 ) and Case 2 considers the disaster rep resented by 5 scenar ios ( s 2 , · · · , s 6 ). In both c ases, the no rmal op erations are repr esented by one scenario ( s 1 ). Each disaster scenario is assumed to start at t = 6 hour and affect th e ca p acity of the d istribution lines, given in T able I, v ia parameter α lts during the rest of the time horizo n. The transition time is defined as T d b 1 ,b 2 ,t = min( | b 1 − b 2 | , d b 1 ,b 2 ) where d b 1 ,b 2 is the number of lines between buses b 1 and b 2 . T he value of lost load is C V oLL b = $5,0 00/MWh. All simulations have been carried out using Gu robi solver v7.5 .1 on Julia 0.6 [17] with an In tel Core i7 2.6 GHz pr ocessor with 16GB of memo ry . T ABLE I. D I S A S T E R S C E NA R I O S F O R C A S E 1 A N D C A S E 2 Case 1 Case 2 Scenari o s 1 2 1 2 3 4 5 6 Probabil ity ω s 0.9 0.1 0.9 0. 02 0.02 0.02 0.02 0.02 Line ∗ l - 4 - 1 12 4 8 9 ∗ Line numbers are assigned for the edge connecti ng bus b with its parent bus as illustrate d in Fig. 1, i.e. l = b . T ABLE II. O B J E C T I V E F U N C T I O N V A L U E S A N D C P U T I M E S F O R T H E B F A N D P H I M P L E M E N TA T I O N Objecti ve function, $ CPU time, s BF PH BF PH Case 1 1032.06 1032.07 35.18 7.99 Case 2 1787.18 1788.41 21,519 2,594 1) Computa tional performa nce: T ab le II compare s the propo sed PH implementa tio n and the bru te- force (BF) im- plementation (i.e. s olving the p ropo sed optimiza tion using Gurobi witho ut any algorith mic enhanc ement) in terms of their comp utational perf o rmance and op tim ality . In both cases, the PH and BF imple m entations attain near ly iden tical values of the o bjective functions and respective minimizers. On the o ther hand, the PH implemen tation solves th e proposed optimization 4x times faster in Case 1, where o n ly 2 scenarios are conside r ed. As the numb er o f scenarios increases to size as in Case 2, th e comp utational gains increase to 8x times. Th is compariso n reveals the comp utational ef ficiency of the PH implementatio n for solving instances with a larger numb er of scenarios. Giv en th e sup erior perfo r mance of the PH algo rithm described in T able I I, it is used to obtain all nume rical results in the rest o f the paper . 2) Effect of the ES mobility: The b enefits of mobile ES units can be revealed by comp a r ing the ir o perations with the case with stationar y ES u nits and withou t ES units. This T ABLE III. T O TA L L O S T L OA D A N D O B J E C T I V E F U N C T I O N F O R I N S TA N C E S W I T H O U T E S , W I T H S TA T I O NA RY E S A N D W I T H M O B I L E E S w/o ES Station ary ES Mobile ES Case T otal lost load, s 1 0 0 0 1 MWh s 2 0.44 0 0 Objecti ve functio n, $ 1185.62 1032.51 1032.07 s 1 0 0 0 T otal s 2 3.47 2.89 2.18 Case lost load, s 3 6.82 6.82 5.42 2 MWh s 4 0.44 0 0 s 5 0 0 0 s 6 0 0 0 Objecti ve functio n, $ 2036.20 1997.95 1788.41 6 7 9 10 8 5 4 3 2 0 1 12 1 1 T&D Interface 13 14 S2 Stationary ES Forecast disaster (t=[6:24]) Mobile ES t=[3:5] S2 Fig. 1: The 15-bus test system from [16 ] used in the case study with the storage placement decisions in t he disaster scenario of Case 1. T ABLE IV. B U S L O C AT I O N S A N D S TA T E - O F - C H A R G E ( S O C ) O F M O B I L E E S U N I T S I N C A S E 2 s Time interva l # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 T 2 2 2 T 1 1 1 1 1 T 0 0 0 T 1 1 1 1 1 T 0 0 Bus # 3 1 T 0 0 T 12 12 12 12 12 12 12 T 0 0 0 0 0 0 T 12 12 12 12 4 1 T 2 T T 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 1 1 T 2 T T 8 8 8 8 8 8 T T T 1 1 1 1 1 1 1 – 2 6 1 T T T T 9 9 9 9 9 T T T T 1 1 1 1 1 1 1 1 1 1 1 0.5 0.5 0.5 0.5 0.5 0.5 0 .49 0.49 0.47 0. 45 0. 42 0 .4 0.38 0.38 0.38 0.35 0.33 0. 33 0.3 0.3 0.3 0.3 0.3 0.3 2 0.54 T 0.69 0.84 1 T 0.9 0.75 0.6 0.45 0.3 T 0.45 0.6 0.75 T 0.6 0.45 0.3 0. 15 0 T 0.15 0.3 SoC, MW h 3 0.65 T 0.8 0.95 T 0.82 0.68 0.55 0.41 0.27 0.14 0 T 0.15 0.3 0.45 0.6 0.75 0.9 T 0.75 0.6 0.45 0.3 4 0.64 T 0.79 T T 0.78 0.76 0.74 0. 71 0.68 0.65 0.62 0.59 0.56 0.53 0.5 0.47 0.44 0.41 0.38 0.35 0. 33 0.31 0.3 5 0.35 0.2 T 0.14 T T 0.27 0.41 0.56 0.7 0.85 1 T T T 0.91 0.79 0.67 0.56 0.45 0.37 0.32 T 0.3 6 0.35 T T T T 0.49 0.62 0.75 0.87 1 T T T T 0.9 0.8 0.8 0.67 0.54 0.43 0.34 0.3 0.3 0.3 Labels ‘T ’ defines that the mobile ES unit is in transit. compariso n is given in T ab le III for Case 1 and Case 2. The effect of the mo bile ES units is rather negligible on the objective fu nction in Case 1 due to the small scenario set. Howe ver , th e use of either station ary or mobile ES u nits reduces the total lost load ( compu te d as P b ∈B P t ∈T p ls bts ) relativ e to the case without ES units. On the other h and, the effect of mobile ES u n its is evident in th e numer ical results for Case 2. First, u sing mobile ES units r e duces the total lost load across all scen arios conside r ed r elativ e to the simulations without ES units and with stationary ES units. Second, the mobile ES units ensure the least-cost op e ration, achieved in part due to the r educed total load shedd ing. If the ES unit were station a ry , as in [9] and [10], it would be installed at bus 6 . this decision chang es for mo bile ES u nits. Fig. 1 describes the ES placement in Case 1. In this case, the m obile ES is placed at bus 1 d uring normal operatio ns and moves to bus 4 in anticipation of the disaster, which is to occur at t = 6 hour and to affect line 4. Follo wing the disaster, th e m obile ES unit b egins to supp ly power to downstream buses 5 and 6. The stor ed energy is eno ugh to av oid load shedd ing u ntil period t = 24 hour s. In Case 2, whic h considers 5 disaster scenarios, the mob ile ES is r o uted an d dispa tch ed as d escribed in T able IV. Sim ilar ly to Case 1, the mo bile ES unit is placed at bus 1 d uring nor mal oper ations ( s 1 ) an d then has fiv e un iq ue routing trajecto ries fo r each disaster scenario ( s 2 · · · s 6 ). These trajectories vary in te r ms of the optimal locations an d d ispatch decisions on m obile ES units, wh ich unde r pins th e importan ce of care f ully calibrating d isaster scen arios (see [4]). V . C O N C L U S I O N This paper describes an app roach to op timize investments of the d istribution system operato r in mo bile ES units. The ability of mob ile ES un its to move b etween different loc ations is used to trade-off the least-c ost o peration s durin g norm a l oper ations and the n eed to en hance power g rid resilience in case of natural disasters. The prop osed optim ization is a two-stage mixed-integer pro gram with b inary reco urse decisio n s, wh ich account fo r the reloc a tio n of mo bile ES u nits und er specific disaster scen arios. Th e pro p osed optimiza tio n is solved using the PH algorithm. The n umerical experiments reveal that the mobile ES reduc e th e o perating c o sts and th e to tal a m ount of load sheddin g cau sed by n atural disasters re lati ve to the cases without E S units o r with stationary ES units. R E F E R E N C E S [1] J.-P . W atson and et al, “Conceptua l frame work for de vel oping resilienc e metrics for the elect ricity , oil, and gas sectors in the united states, ” Sandia Nat. Lab . , T ech. Rep. SAND2014-18019 , 2014. [2] S. Neumayer and E. Modiano, “ Assessing the effe ct of geograph icall y correla ted fa ilures on interc onnecte d power -communication network s, ” in 2013 IEEE SmartGridComm Conf. IEE E , 2013, pp. 366–371. [3] Y . W ang, C. Chen, J. W ang, and R. Baldick , “Researc h on resilience of po wer systems under natural disastersa re view , ” IE EE Tr ans. on P ower Systems , vol. 31, no. 2, pp. 1604–1613, 2016. [4] E. Rappapo rt and et al, “ Adva nces at the national hurricane center , ” W eather and F oreca sting , vol. 24, no. 2, pp. 395–419, 2009. [5] S. Lei, J. W ang, C. Chen, and Y . Hou, “Mobile emer genc y generator pre-posit ioning and real-ti m e allocatio n for resilient response to natural disasters, ” IE EE T rans. on Smart Grid , 2016. [6] Z. W ang and J. W ang, “Self-healing resilient distributio n systems based on sectionali zatio n into microgrids, ” IE EE T rans. on P ower Systems , vol. 30, no. 6, pp. 3139–3149, 2015. [7] N. Smith and R. McCann, “ Analysis of distrib uted generation sources and load shedding schemes on isolate d grids case study: The bahamas, ” in 2014 Int. Conf. R en. E n. Res. & A ppl. IEEE, 2014, pp. 301–306. [8] Con edison de veloping energy s torage system on wheels. [Online]. A vaila ble: https:/ /goo.gl/hMr fiA [9] H. Pa nd ˇ zi ´ c, Y . W ang, Y . Dvo rkin, and D. Kirschen, “Near -optimal method for siting and sizing of distribut ed storage in a transmission netw ork, ” IEEE T rans. Pwr . Syst. , vol. 30, no. 5, pp. 2288–2300, 2015. [10] R. Fernandez-Bla nco and et al, “Optimal ener gy storage s iting and sizing: A wecc case study , ” IEEE T rans. Sust. En. , vol. 8, no. 2, pp. 733–743, April 2017. [11] H. D. Sherali and X. Zhu, T wo-Stag e Stochastic Mixed-Inte ger Pr o- grams: Algorithms and Insights . Boston, MA: Springer US, 2009. [12] J.-P . W atson and D. L. W oodruf f, “Progressiv e hedging innov ations for a cl ass of stoch astic mixed-i nte ger resource al locat ion problems, ” Comp. Manag . Sc. , vol. 8, no. 4, pp. 355–370, Jul. 2010. [13] S. Ryan and et al, “T o ward scalable, parall el progressi ve hedging for stochasti c unit commitment, ” in P ower and Energy Society General Meeti ng (PES), 2013 IEEE . IEEE, 2013, pp. 1–5. [14] M. Fari var and S. Low , “Branch fl ow model: Relaxati ons and con vexi- fication , ” IEEE T rans. P wr . Syst. , vol. 28, no. 3, pp. 2554–64, 2013. [15] M. A. Ortega -V azquez, “Optimal schedulin g of elect ric vehicle char g- ing and vehi cle-t o-grid s ervices at household lev el including battery degra dation and price uncertaint y , ” IET Generation, T ransmission & Distrib ution , vol. 8, no. 6, pp. 1007–1016, 2014. [16] A. Pap av asiliou, “ Analysis of distribu tion locational marginal prices, ” IEEE T rans. on Smart Grid , 2017. [17] M. Lubin and I. Dunning, “Computing in operations research using julia , ” INFORMS J . on Comp. , vol. 27, no. 2, pp. 238–248, 2015.