Optimization of Energy Consumption Schedule ofResidential Loads and Electric VehiclesbyEnxin YaoM.S., Peking University, 2011B.S., Peking University, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Electrical and Computer Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)June 2016c© Enxin Yao, 2016AbstractIn the current electrical grid, utility companies have begun to use demand side management(DSM) programs and time-of-use (TOU) pricing schemes to shape the residential loadprofile. However, it is difficult for the residential users to respond to the pricing signal andmanually manage the operation of various household appliances. Hence, the autonomousenergy consumption scheduling of residential loads and electric vehicles (EVs) is necessaryfor the users to benefit from the DSM programs. In this thesis, we propose differentalgorithms to schedule the energy consumption of residential loads and EVs, and provideancillary services to the electrical grid. First, we study the DSM for areas with highphotovoltaic (PV) penetration. Since many rooftop PV units can be integrated in thedistribution network, the voltage rise issue occurs when the reverse power flow from thehouseholds to the substation is significant. We use stochastic programming to formulatean energy consumption scheduling problem, which takes into account the voltage rise issueand the uncertainty of the power generation from PV units. We propose an algorithm bysolving the formulated problem and jointly shave the peak load and reduce the reversepower flow. Subsequently, we study using the EVs to provide the frequency regulationservice. We formulate a problem to schedule the hourly regulation capacity of the EVsusing the probabilistic robust optimization framework. Our formulation takes into accountthe limited battery capacity of the EVs and the uncertainty of the automatic generationcontrol (AGC) signal. An efficient algorithm is proposed to solve the formulated problembased on duality. Last but not least, we study the market participation of an aggregatoriiAbstractwhich coordinates a fleet of EVs to provide frequency regulation service to an independentsystem operator (ISO). The two-settlement market system (i.e., the day-ahead market(DAM) and real-time market) is considered. We analyze two types of DAMs based on themarket rules of New York ISO and California ISO. We formulate a problem to determinethe bid for the aggregator in the DAM using stochastic program and conditional value atrisk. Efficient algorithms are proposed to tackle the formulated problem.iiiPrefaceChapters 2–4 encompass work which has been published or currently under review. Thecorresponding papers are under the supervision of Professor Vincent W.S. Wong. The pa-pers corresponding to Chapters 2–4 are also co-authored with Professor Robert Schober,who has provided constructive comments and suggestions on all the papers. Dr. PedramSamadi has co-authored the paper corresponding to Chapter 2 and provided helpful com-ments.For all chapters, I hereby declare that I am the first author of the corresponding papers.The following publications describe the work completed in this thesis.Journal Papers, Accepted• Enxin Yao, Pedram Samadi, Vincent W.S. Wong, and Robert Schober, “ResidentialDemand Side Management Under High Penetration of Rooftop Photovoltaic Units,”IEEE Trans. on Smart Grid, vol. 7, no. 3, pp 1597 - 1608, May 2016.• Enxin Yao, Vincent W.S. Wong, and Robert Schober, “Robust Frequency RegulationCapacity Scheduling Algorithm for Electric Vehicles,” accepted for publication inIEEE Trans. on Smart Grid, 2016.Journal Paper, Submitted• Enxin Yao, Vincent W.S. Wong, and Robert Schober, “Market Participation of Elec-tric Vehicles for Frequency Regulation Service Using an Aggregator,” submitted,2016.ivPrefaceConference Papers, Published• Enxin Yao, Pedram Samadi, Vincent W.S. Wong, and Robert Schober, “Tacklingthe Photovoltaic Integration Challenge in the Distribution Network with DeferrableLoad,” in Proc. of IEEE International Conference on Smart Grid Communications(SmartGridComm), Vancouver, Canada, Oct. 2013.• Enxin Yao, Vincent W.S. Wong, and Robert Schober, “A Robust Design of ElectricVehicle Frequency Regulation Service,” in Proc. of IEEE International Conferenceon Smart Grid Communications (SmartGridComm), Venice, Italy, Nov. 2014.• Enxin Yao, Vincent W.S. Wong, and Robert Schober, “Risk-Averse Forward Contractfor Electric Vehicle Frequency Regulation Service,” in Proc. of IEEE InternationalConference on Smart Grid Communications (SmartGridComm), Miami, FL, Nov.2015.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Residential Demand Side Management . . . . . . . . . . . . . . . . . . . . 21.1.1 DSM under Integration of Rooftop Photovoltaic Units . . . . . . . 31.2 Electric Vehicle Frequency Regulation Service . . . . . . . . . . . . . . . . 61.2.1 Frequency Regulation Service . . . . . . . . . . . . . . . . . . . . . 61.2.2 Characteristics of the EVs . . . . . . . . . . . . . . . . . . . . . . . 81.2.3 Performance-based Compensation Scheme . . . . . . . . . . . . . . 91.2.4 Two-settlement Market System . . . . . . . . . . . . . . . . . . . . 10viTable of Contents1.3 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.4 Summary of Results and Contributions . . . . . . . . . . . . . . . . . . . 131.5 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 DSM under High Penetration of Rooftop PV Units . . . . . . . . . . . 162.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.2 System Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.2.1 Voltage Rise Problem . . . . . . . . . . . . . . . . . . . . . . . . . 202.2.2 Cost Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.3 Energy Consumption Scheduling with Uncertain PV Power Generation . 262.3.1 Nested Stochastic Formulation . . . . . . . . . . . . . . . . . . . . 262.3.2 Sample Average Approximation . . . . . . . . . . . . . . . . . . . . 282.3.3 Convex Piecewise Linear Objective Function . . . . . . . . . . . . 302.3.4 Non-convex Piecewise Linear Objective Function . . . . . . . . . . 312.4 Performance Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352.4.1 Computational Complexity . . . . . . . . . . . . . . . . . . . . . . 382.4.2 Performance Gains for Voltage Rise Mitigation . . . . . . . . . . . 402.4.3 Benefits Regarding Electricity Bill and PAR . . . . . . . . . . . . . 402.4.4 Impact of Forecast Error . . . . . . . . . . . . . . . . . . . . . . . 432.4.5 Impact of PV Power Generation and Parameter η . . . . . . . . . . 462.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473 Frequency Regulation Capacity Scheduling for EVs . . . . . . . . . . . 503.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 503.2 System Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.2.1 Randomness of the AGC Signal . . . . . . . . . . . . . . . . . . . . 533.2.2 Performance-based Frequency Regulation Compensation . . . . . . 55viiTable of Contents3.3 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 583.3.1 Duality-based Problem Transformation . . . . . . . . . . . . . . . 653.3.2 Probability Bound and Selection of Parameter η . . . . . . . . . . 683.3.3 Trajectory of the AGC Signal in Hour h . . . . . . . . . . . . . . . 693.3.4 Stochastic Arrival and Departure of the EVs . . . . . . . . . . . . 723.3.5 Aggregate Capacity in the Day-ahead Market . . . . . . . . . . . 733.3.6 Algorithm and Implementation Issues . . . . . . . . . . . . . . . . 763.4 Performance Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 783.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 844 Market Participation of EVs for Frequency Regulation Service . . . . 874.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 874.2 System Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 894.2.1 Bid of the Aggregator . . . . . . . . . . . . . . . . . . . . . . . . . 914.2.2 DAM and RTM . . . . . . . . . . . . . . . . . . . . . . . . . . . . 924.2.3 Profit of the Aggregator . . . . . . . . . . . . . . . . . . . . . . . . 934.3 Problem Formulation and Algorithm . . . . . . . . . . . . . . . . . . . . 954.3.1 Risk-Averse Bid of the Aggregator . . . . . . . . . . . . . . . . . . 954.3.2 Aggregate Regulation Capacity of EVs . . . . . . . . . . . . . . . 994.4 Performance Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1034.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1085 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . 1095.1 Summary of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1095.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113viiiList of Tables2.1 List of deferrable loads . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36ixList of Figures1.1 The cumulative installed capacity of PV units in Canada. The data isobtained from [21]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.1 A schematic diagram of the proposed energy consumption scheduling scheme.The DNO performs the power flow analysis and provides the ECS with in-formation regarding the voltage rise (e.g., the sensitivity of the voltage withrespect to the household energy export) at the beginning of the day. . . . . 182.2 One line diagram of a 33 bus radial distribution test system. The householdsare connected to the buses. The voltage rise problem occurs when manyhouseholds have substantial positive energy export. The DSM programencourages households to shift loads to hours with high solar radiation tomitigate the voltage rise problem. . . . . . . . . . . . . . . . . . . . . . . . 212.3 Piecewise linear cost function for a household with PV unit. We considerthe case ps(t) ≤ pe(t) in this figure, i.e., the feed-in tariff does not exceedthe energy consumption price. On the other hand, Some utility companiespromote the rooftop PV installation aggressively and have ps(t) > pe(t). . . 252.4 Non-convex piecewise linear cost function. This case occurs when the feed-in tariff of the utility company is higher than the energy consumption price(i.e., ps(τ) > pe(τ)). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32xList of Figures2.5 Running time of Algorithm 2.1 with LP and MIP solvers under Price1 andPrice2, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 382.6 Comparison of the household voltage magnitude profile of the proposed al-gorithm and the HDA algorithm. The voltage magnitude is shown in pu,whereas the base value of pu is the standard voltage (e.g., 120 volt for house-holds). The voltage limit is 1.05 pu. In the proposed algorithm, the voltagerise is avoided by reducing the household energy export during hours withhigh solar power generation. . . . . . . . . . . . . . . . . . . . . . . . . . . 412.7 Comparison of the monthly electricity bill for different amount of dailyrooftop PV power generation. . . . . . . . . . . . . . . . . . . . . . . . . . 422.8 Comparison of the aggregate load PAR of the proposed algorithm and theHDA algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 432.9 The voltage rise versus the forecast error. The y-axis is the summationof the voltage rise over time slots and buses, which can be written as∑t∈T∑n∈B(Vn(t) − V¯ )+, where V¯ is the upper limit of the allowed volt-age variation and (·)+ = max{ · , 0}. . . . . . . . . . . . . . . . . . . . . . . 452.10 The saving on the monthly electricity bill versus the forecast error. They-axis is the percentage of the saving. The baseline of the percentage isthe saving when the energy generation and loads are perfectly known. Thesaving decreases as the forecast error increases. . . . . . . . . . . . . . . . . 462.11 The voltage rise versus the daily power generation of the rooftop PV unit.The y-axis is the summation of the voltage rise over time slots and buses,which can be written as∑t∈T∑n∈B(Vn(t)− V¯ )+. We choose η = 2× 10−4. 47xiList of Figures2.12 The impact of parameter η. The y-axis is the summation of the voltagerise∑t∈T∑n∈B(Vn(t) − V¯ )+. η is a parameter of the proposed algorithm.The performance of the cases when the ECS is not installed or the HDAalgorithm is used are not affected by η. . . . . . . . . . . . . . . . . . . . . 483.1 The EVs provide frequency regulation service by following an AGC signalissued by the ISO. The AGC signal is generated by the ISO in real-time andis random. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 533.2 Joint distribution of the hourly regulation up and regulation down com-ponents of an AGC signal. The distribution is obtained by analyzing theAGC signal data records from PJM for 2,208 hours. The figure reveals therandomness of the regulation up and regulation down components. . . . . . 553.3 Different trajectories of the AGC signal and the corresponding SOC of theEV battery. The initial SOC is set to be 0.5. . . . . . . . . . . . . . . . . . 713.4 A diagram of different components for the EV frequency regulation service. 773.5 A sample of the AGC signal obtained from PJM. . . . . . . . . . . . . . . 793.6 The average hourly prices of frequency regulation service from PJM in Dec.2014. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 803.7 Daily revenue versus the maximum hourly charged energy. . . . . . . . . . 833.8 (a) The revenue with respect to parameter η. (b) The regulation capacityand the performance score as a function of parameter η. . . . . . . . . . . 843.9 Probability Pi obtained via simulations and its lower bound obtained from(3.17). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 853.10 Revenue as a function of the price of the battery. Note that we assumed abidirectional charger in this figure. . . . . . . . . . . . . . . . . . . . . . . 864.1 Sample profiles of EVs’ connection with the power grid. . . . . . . . . . . . 88xiiList of Figures4.2 The aggregator coordinates the EVs to participate in the DAM and RTM,and provides frequency regulation service to the ISO. The aggregator needsto submit a bid to obtain a contract from the ISO in the DAM. . . . . . . 904.3 The prices in the DAM and RTM for the frequency regulation service inCAISO and NYISO. The y-axis is the sum of the prices of the regulation upand regulation down capacities. . . . . . . . . . . . . . . . . . . . . . . . . 1034.4 Expected profit versus maximum hourly charged energy. . . . . . . . . . . 1054.5 Expected profit of aggregator divided by number of EVs versus the numberof EVs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1064.6 CVaR of aggregator divided by number of EVs versus the number of EVs. 1064.7 Efficient frontier between the expected profit and the CVaR. This figureshows the tradeoff between the expected profit and the CVaR tuned byparameter β. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107xiiiList of AbbreviationsACE Area control errorAGC Automatic generation controlAMI Advanced metering infrastructureANSI American National Standards InstituteCAISO California ISOCLT Central limit theoremCPU Central processing unitCVaR Conditional value at riskDAM Day-ahead marketDER Distributed energy resourceDLC Direct load controlDNO Distribution network operatorDSM Demand side managementECS Energy consumption schedulerERCOT Electric Reliability Council of TexasESS Energy storage systemEV Electric vehicleFERC Federal Energy Regulatory CommissionGW GigawattHAN Home area networkxivList of AbbreviationsHDA HVAC load direct activationHVAC Heating, ventilation, and air conditioningIEEE Institute of Electrical and Electronics EngineersISO Independent system operatorLP Linear programmingMAPE Mean absolute percentage errorMIP Mixed integer programmingNYISO New York ISOPAR Peak-to-average ratioPJM Pennsylvania New Jersey Maryland InterconnectionPV Photovoltaicpu per unitRMCCP Regulation market capacity clearing priceRMPCP Regulation market performance clearing priceRTM Real-time marketSAA Sample average approximationSAE Society of Automotive EngineersSOC State of chargeTOU Time-of-useU.S. United StatesVaR Value at riskxvAcknowledgmentsFirst and foremost, I would like to express my gratitude to my supervisor, Professor VincentW.S. Wong, for his patience, encouragement, and invaluable advices during my Ph.D.study. I am lucky to have a supervisor who teaches me patiently and provides extensivesuggestions on my research. I greatly thank Professor Robert Schober for his constructiveadvice and comments through my research. I am grateful to Dr. Pedram Samadi who hashelped me at the start of my research.Also, I thank the members of my doctoral committee, Professor Lutz Lampe and Pro-fessor Juri Jatskevich, for the time and effort in evaluating my work and providing feedbackand suggestions.Many thanks go to all the colleagues and friends in room 4090 of Kaiser Building.I am always indebted to my wife and my parents for their love and support.This work was supported by the Natural Sciences and Engineering Research Council(NSERC) of Canada.xviDedicationTo my mother Suzhi KangxviiChapter 1IntroductionIn the current electrical grid, 20% of the generation capacity are used to meet peak load andare active in only 5% of the time. The demand side management (DSM) is a promisingmethod to enable more efficient utilization of the generation assets [1, 2]. Many util-ity companies have introduced DSM programs and time-of-use (TOU) pricing schemesto encourage users to shift their energy consumption from peak hours to off-peak hours[3, 4]. However, for the residential users, there are various appliances while the energyconsumption of each appliance is small. Hence, it is necessary to have autonomous energyconsumption scheduling algorithms to manage the operation of the household applianceswithout user intervention in order to allow the users to benefit from the DSM programs.The DSM programs can be used to encourage users to change load profile in order toachieve various objectives in the electrical grid. In the emerging smart grid, one of the keychallenges is how to support the rapid growth of the photovoltaic (PV) units and electricvehicles (EVs). The DSM programs can play an important role to help the integrationof PV units and EVs into the electrical grid. In areas which have high penetration ofrooftop PV units, the voltage rise problem has emerged as a major obstacle to integratethe PV units into the distribution network. In those areas, the DSM programs can beused to encourage the household users to shape load profile and mitigate the voltage riseproblem. On the other hand, the usage of the EVs increases rapidly and the EV becomesa promising resource to provide the frequency regulation service. The experiments in [5, 6]with Pennsylvania New Jersey Maryland Interconnection (PJM) demonstrate that EVs1Chapter 1. Introductionare capable to provide the frequency regulation service to an independent system operator(ISO) by changing EVs’ real-time charging or discharging power. Due to the limited batterycapacity of an EV, it is necessary to have aggregators to coordinate a large fleet of EVs inorder to meet the minimum capacity requirement to participate in the market of an ISO.The scheduling of the EVs’ frequency regulation capacity and the market participation ofthe aggregator have not been fully addressed.In this thesis, we propose different algorithms to schedule the energy consumption ofthe residential loads and EVs, and provide ancillary services to the electrical grid. Thestochastic optimization, robust optimization, and conditional value at risk (CVaR) areused in the design of the proposed algorithms. The rest of this chapter is organized asfollows. In Section 1.1, we provide an overview of the residential DSM and the voltage riseproblem for the integration of rooftop PV units. We introduce the EV frequency regulationservice in Section 1.2. We summarize the related works in Section 1.3. The contributionsmade in this thesis are summarized in Section 1.4, and the thesis organization is providedin Section 1.5.1.1 Residential Demand Side ManagementDSM can be used by the utility company to shape the load profile and reduce energyconsumption during peak hours. Peak shaving can reduce the energy generation cost andthe cost of plant investment [7]. There are generally two methods to shape the load, namelythe direct load control (DLC) and price-based DSM. In DLC, the utility company has acontract with the users and can remotely control the operation of users’ loads to reducethe consumption during a period of time [8]. On the other hand, the price-based DSMuses a dynamic pricing scheme for which the utility company sets a high price during peakconsumption hours to encourage the users to shift their load [9].2Chapter 1. IntroductionThe DSM programs have been widely used. The Hydro One cooperation in the provinceof Ontario in Canada has introduced TOU pricing schemes which are effective on weekdays,whereas the price of electricity varies from 17.5 cents to 8.3 cents for peak hours and off-peak hours, respectively [3]. The DSM programs are also used in the state of California[4] and the state of Illinois [10] in the United States (U.S.).Residential DSM has attracted significant attention from the research community [11–16]. An important challenge for residential DSM is that it is difficult for household usersto respond to the pricing signals [11]. With the development of the advanced meteringinfrastructure (AMI), smart meters can receive the dynamic prices (e.g., TOU, real-timeprices) from the utility company and have the potential to manage the energy consumptionof the household loads. In [11], an autonomous energy consumption scheduler (ECS) isproposed to help users make price-based control decisions. The ECS can be implementedwith the smart meter to receive price signals and control the loads. Moreover, the dynamicpricing scheme plays an important role in the DSM programs and the optimization of theresidential electricity retailing prices has attracted lots of attention [17–19]. An algorithmis proposed in [17] to determine the day-ahead TOU prices with two-stage stochastic opti-mization. The uncertain impact of the prices on the load profile is modeled with stochasticapproximation in [18]. The individual user’s responsiveness to the dynamic prices is takeninto account in [19] whereas online optimization is used to determine the real-time prices.1.1.1 DSM under Integration of Rooftop Photovoltaic UnitsIn recent years, the use of energy generation from PV units has proliferated. The cumula-tive installed capacity of PV units in the world has increased from 40 gigawatt (GW) to180 GW from year 2010 to 2014 [20]. Fig. 1.1 shows the rapid growth of installed capacityof PV units in Canada.3Chapter 1. Introduction0200400600800100012001400Installed Capacity (MW)Figure 1.1: The cumulative installed capacity of PV units in Canada. The data is obtainedfrom [21].The residential users are allowed to install PV units locally (i.e., the rooftop PV units).The rooftop PV units can reduce the electricity bills of households and may promote a lowcarbon footprint lifestyle. According to [22], around 784,000 homes and businesses in theU.S. have installed PV units by the middle of year 2015. In Canada, the installed capacityof grid-connected distributed PV units is around 30% of the total installed capacity of PVunits [21].The households with rooftop PV units may not only consume energy but also exportenergy to the electrical grid. Thereby, many utility companies have introduced net meteringprograms to determine the electricity bill for households with rooftop PV units [23, 24].Net metering programs typically have a feed-in tariff and allow households to sell theirextra energy to the electrical grid.The DSM under the net metering programs recently has attracted lots of attention[25–27]. The works in [25–27] focus on the case where a few households are equipped with4Chapter 1. Introductiondistributed energy resources (DER) and are encouraged to export their generation to thegrid. However, if a large number of users are equipped with PV units, the results in [25–27]may not be directly applicable because the reverse power flow in those areas may causethe voltage rise problem. Therefore, in areas with high penetration of PV units, new DSMprograms are necessary.The voltage rise problem is an important challenge for the integration of a large numberof rooftop PV units into the electrical grid [28, 29]. Traditional voltage control strategiesassume the unidirectional power flow from the substation to the households [30]. A substan-tial reverse power flow from households to the substation can cause the voltage magnitudeof some of the households to exceed the upper limit of the allowed voltage variation. Thisis referred to as the voltage rise problem. For example, in North America, the AmericanNational Standards Institute (ANSI) C84.1 standard [31] requires the voltage magnitudein the distribution network to be no more than 1.05 per unit (pu) of the normal state.In countries such as Germany which have high penetration of rooftop PV units, thevoltage rise problem has already emerged and different mechanisms have been proposed totackle the problem [32, 33]. The distribution network operators (DNOs) have upgraded thetransformers (e.g., with on-load tap changer transformers) and have enhanced the feederto host more PV units in some areas [32]. Moreover, the reverse power flow from PV unitscan be controlled by adjusting the active and reactive power of PV inverters [33]. The workin [34] proposes an optimal dispatch algorithm to set the active power and reactive powerof PV inverters. In particular, the German Renewable Energy Sources Act revision in 2012[35] enforces energy generation curtailment for rooftop PV units. The DNO measures thevoltage magnitude and sends control signals to the PV units for curtailment. The PVunits are required to either install remote controllers that curtail the generation based onthe control signal received from the DNO or permanently limit the active power feed-in to5Chapter 1. Introduction70% of the installed capacity [32].The voltage rise problem should be taken into account for the DSM programs in areaswith high penetration of PV units. The DSM programs can encourage users to shift load tohours with high solar radiation in order to reduce the reverse power flow from householdsto substations and mitigate the voltage rise issue [36, 37]. Compared to the generationcurtailment, the DSM program provides a promising alternative method to mitigate thevoltage rise issue without losing revenue of energy generation from PV units in the net-metering programs.1.2 Electric Vehicle Frequency Regulation ServiceIn recent years, there have been significant efforts to replace the combustion engine vehicleswith EVs. The cumulative sales of EVs in the U.S. have increased from 20,000 to 350,000from year 2011 to year 2015 [38]. The surge in EV sales can be partially attributed to therapidly falling cost of batteries. From year 2007 to year 2014, the cost of battery in EVs hasreduced from $1000 per kWh to $500 per kWh [39]. Most EVs (plug-in EVs and plug-inhybrid EVs) use the electricity drawn from the electrical grid to reduce the consumptionof fossil fuels. EVs can be regarded as distributed energy storage systems when they areconnected with the electrical grid. Thereby, EVs can be coordinated to provide frequencyregulation service to an ISO, such as the California ISO (CAISO).1.2.1 Frequency Regulation ServiceFrequency regulation service helps ISOs to keep the utility frequency around the nominalvalue (e.g., 50 Hertz or 60 Hertz) by compensating short term mismatch between generationand load. The real-time deviation of the utility frequency in the electrical grid dependson the power generation and the demand. That is, the utility frequency tends to increase6Chapter 1. Introductionwhen the generation is higher than the demand and tends to decrease vice versa. TheISOs monitor the real-time deviation of the utility frequency in their administrated areasand are responsible to maintain the utility frequency within a predefined range. Thereby,the ISO determines a real-time area control error (ACE) based on the frequency deviation,along with the power imports and exports with adjacent regions. Then, the ISO issues anautomatic generation control (AGC) signal based on the ACE every few seconds. The ISOsends the AGC signal to some power generators such that their real-time output power areadjusted according to the AGC signal.In the current electrical grid, the frequency regulation service is mostly provided bypower generators [40]. However, the ramping rate of the generators is typically slower thanthe rate required by the AGC signal, which results in high regulation capacity procurementsof the ISOs [41–43]. The pilot projects conducted by PJM [44] show that batteries andflexible loads can change their real-time energy consumption rapidly and follow the AGCsignal closely. The results from [45] show that resources with fast ramping rates suchas energy storage systems (ESS) and some of the loads can provide frequency regulationservice efficiently and suggest that CAISO can reduce its regulation capacity procurementby using fast ramping resources. Similar conclusion is obtained in the analysis of themarket participation of ESSs for the Electric Reliability Council of Texas (ERCOT) [42].The ESSs have begun to be used to provide the frequency regulation service by someISOs. Since their limited capacity, the ESSs are able to provide the frequency regulationservice only if their state of charge (SOC) are maintained within certain ranges. Theauthors in [46] propose to use filtered AGC signal to reduce the impact of the frequencyregulation service on the SOC of ESSs. Two contractual frameworks have been designedin [47] which constrain the ISOs to generate SOC-neutral and energy-neutral AGC signal,respectively. Recently, PJM has introduced a new filtered signal for ESSs, which is referred7Chapter 1. Introductionas RegD signal [48]. The hourly regulation up and regulation down components are almostidentical for the RegD signal. As a result, the RegD signal can reduce the impact of thefrequency regulation service on the SOC of ESSs. PJM allows the ESSs to follow the RegDsignal, whereas the power generators follow another signal (RegA) [49].Some of the flexible loads such as EVs can provide the frequency regulation service effi-ciently and the curtailment of the load is much faster than the ramping rate of the thermaland hydro power plants [50]. In year 2008, the Federal Energy Regulatory Commission(FERC) issued Order 719 [51], which required that the ISOs must accept bids from flexibleloads to provide frequency regulation service, besides other services, and must treat flexibleloads equally as other resources. The authors in [50] evaluate the potential of different loadsto provide the frequency regulation service. The heating, ventilation, and air conditioning(HVAC) loads in buildings are analyzed in [52], whereas the flexibility of the HVAC loadsis quantified and a model predictive control scheme is proposed to exploit the flexibilityto follow the supply. The analysis in [43] shows that the available frequency regulationcapacity of flexible loads can satisfy the current maximum regulation procurement (600MW) of CAISO. The EV is one of the flexible loads to provide the frequency regulationservice and has attracted lots of attention [5, 6, 53–55]. The experiments and analysis in[5, 6] show that EVs are able to follow the AGC signal by changing their real-time chargingand discharging rate quickly, and earn revenue for EV owners by providing the frequencyregulation service.1.2.2 Characteristics of the EVsEVs are intrinsically dispersed and each EV has limited regulation capacity. An aggregatoris typically used to coordinate a large fleet of EVs (e.g., thousands) in order to satisfy theminimum capacity (e.g., 0.1 MW) required to enter the wholesale market of an ISO. In8Chapter 1. Introduction[53], the role of the aggregator in market participation is discussed and several businessmodels are proposed. A virtual battery model is proposed in [54] for the aggregation of afleet of EVs and other residential loads. The authors in [55] evaluate the potential of EVto satisfy the frequency regulation demand in California.EVs can provide regulation service with both bidirectional chargers and unidirectionalchargers. The literature which considers the bidirectional chargers includes [5, 56–60]. Inthese papers, it is generally assumed that the EVs provide regulation up / down servicesby discharging / charging the EV batteries. On the other hand, a framework of unidirec-tional EV frequency regulation service is proposed in [61] to avoid the battery degradationcost occurred during discharging. The unidirectional EV frequency regulation service isconsidered to be the first step to use EVs to provide the frequency regulation service [61–64]. Thereby, the available regulation capacity of an EV depends on the allowed chargingrate of its charger. EV chargers have different maximum allowed charging rates and areclassified accordingly into level 1 (≤ 1.9 kW), level 2 (1.9 ∼ 19.2 kW), and level 3 (≥ 19.2kW) [65].1.2.3 Performance-based Compensation SchemeThe performance-based compensation scheme for the frequency regulation service is re-cently introduced by the ISOs in the U.S.. The FERC issued Order 755 [66] in Oct. 2011,which requires ISOs to introduce the performance-based compensation scheme in theirfrequency regulation market [67]. Under this scheme, the compensation for the frequencyregulation service depends on both the regulation capacity provided and the performanceof the service provider (e.g., EVs) in following the AGC signal, i.e., whether the EVstrack the AGC signal closely or not. The term performance-based compensation schemeis used to distinguish this new compensation scheme from the scheme previously used by9Chapter 1. IntroductionISOs where the revenue only depends on the capacity. The goal is to encourage the usageof fast ramping resources such as EVs and batteries in the frequency regulation market.Performance-based compensation schemes have been implemented by most ISOs in theU.S., e.g., PJM, CAISO, New York ISO (NYISO), and Midcontinent ISO. They provideeconomic incentives to encourage the resources to respond to the AGC signal quickly, ac-curately, and reliably. However, the significance of the performance-based compensationscheme on the frequency regulation service provided by the EVs and ESSs has not beeninvestigated.1.2.4 Two-settlement Market SystemThe day-ahead market (DAM) and real-time market (RTM) (i.e., the two-settlement sys-tem of the power markets) are widely used for the trade of frequency regulation servicein the U.S. [68–70]. On the day before the operation hours, the market participants (e.g.,power plants) submit their bid in the DAM to indicate their available capacity and askingprices. The ISO collects the bids and announces the market clearing prices of the regulationcapacity. A contract is awarded by the ISO if the aggregator is a winner in the bid.There are two types of DAMs with different market rules in the two-settlement mar-ket system. For the first type of DAM, the market participant needs to ensure that itsregulation capacity will meet the amount specified in the contract. This is the case inCAISO [68]. On the other hand, for the second type of DAM, if the capacity of the marketparticipant is insufficient compared to the amount specified in its contract, it can pay apenalty to the ISO to honor the contract. This is the case in NYISO [69].The aggregator can also participate in the RTM on the next day. If the aggregator hasextra capacity compared to the amount specified in its contract, the aggregator may sellextra capacity in the RTM. The market participation of aggregator in DAM and RTM is10Chapter 1. Introductiondiscussed in [62, 63, 71–73]. The participation of an aggregator in Iberian DAM is analyzedin [71]. In [72], an algorithm is proposed to leverage ESSs for EV charging under the marketparticipation. The bidding algorithm proposed in [63] takes into account the AGC signaland the market clearing prices. In [62], a multi-timescale control algorithm to determinethe charging rate of the EV under the market participation is proposed. However, the twotypes of DAMs have not been considered in the literature. We are interested in how doesthe aggregator determine its bid for the frequency regulation service in the DAM.1.3 Related WorkRecently, there has been a growing research interest on managing household appliancesand EVs to provide services to utility companies. In this chapter, we summarize some ofthe related work which generally fall into three categories.The first category contains the studies on using the load shifting and TOU prices toreduce the peak demand. The examples of the studies in this category are [7–9, 11–16].The authors in [12] categorize the residential loads into four types and formulate a utilitymaximization problem to schedule the operation of different types of loads. The householdappliances, such as washing machines and dryers are analyzed in detail in [13] and theresults show the potential of DSM programs on reducing the peak load. As estimatedin [13], the peak demand in the province of Ontario can be reduced by 125 MW if theoperation of some of household appliances can be delayed for a short period of time.The distributed wind power and the uncertainty of its power generation are consideredin the load scheduling algorithm in [14]. The work in [15] models the thermal dynamicsof HVAC loads and an algorithm is proposed to manage the operation of HVAC loadsunder the dynamic prices. The scheduling algorithm proposed in [16] takes into accountthe uncertainty of the load, the real-time prices, and inclining block rate tariffs. There are11Chapter 1. Introductionalso some other papers discussing how to design the optimal TOU prices or other dynamicpricing schemes, under the context of the load shifting, see [17–19]. This thesis is differentfrom the above works since peak shaving is jointly considered with the voltage rise issue.The second category of the related work is about the study of integration of PV unitsin the distribution network. This category contains [28, 29, 32–34, 36]. These referencesmainly focus on the challenges to integrate the rooftop PV units, especially the analysis ofthe voltage rise issue. Different approaches to handle the voltage rise issue are proposed in[28, 29, 32–34, 36]. Using the DSM to handle the voltage rise is considered in some of thereferences but has not been discussed in detail. On the other hand, there are some otherreferences [25–27] which study the DSM for households with DER such as PV units. In[25], a coordinated algorithm is proposed to minimize the users’ payments. The proposedalgorithm in [25] controls both the household loads and DERs. A game-theoretic approachis proposed in [26] which jointly controls the loads, DERs, and ESSs. In [27], an energyconsumption scheduling algorithm is proposed for a PV-based microgrid, where the hourlyprobabilities of different loads to be used are considered. However, the voltage rise issuehas been neglected in [25–27].The third category of the related work contains the studies of using EVs to providethe frequency regulation service. References [5, 6, 44, 46, 47, 50, 52–64, 71, 72] fall intothis category. There have been experiments on the EV frequency regulation service [5, 44]and also evaluation of the economical feasibility of the service [46, 47]. Different chargingcontrol strategies to provide the frequency regulation service are proposed in [58–62]. Thisthesis is different from [5, 6, 44, 46, 47, 50, 52–64, 71, 72] because we study the schedulingof regulation capacity under the performance-based compensation scheme. Moreover, themarket participation of the aggregator in the two-settlement market system is modeled inthis thesis.12Chapter 1. Introduction1.4 Summary of Results and ContributionsThis thesis proposes several algorithms to schedule the energy consumption of residentialloads and EVs, and provide ancillary services to the electrical grid. The results are dividedinto three chapters. The results and contributions in each chapter are as follows.1. In Chapter 2, we propose an algorithm to schedule the energy consumption of house-hold loads. The proposed algorithm can be used in the DSM programs in areas withhigh penetration of rooftop PV units and takes into account the voltage rise issuein those areas. We formulate a stochastic energy consumption scheduling problemto minimize an objective function which comprises both the electricity bill and anexternal cost for the voltage rise. We model the external cost using the sensitivity ofthe voltage magnitudes with respect to the household energy export and power flowanalysis. We solve the formulated problem using linear programming (LP) whenthe cost function is convex piecewise linear. Otherwise, when the cost function isnon-convex piecewise linear, we use mixed-integer programming (MIP) to solve theproblem. We evaluate the performance of the proposed algorithm using simulationswith a 33 bus radial distribution test system. The results confirm that the proposedalgorithm reduces the energy expenses of the users. The results also show that thevoltage rise problem can be mitigated with the proposed algorithm while the effect issubject to the forecast error, energy generation of PV units, and parameter selection.Part of the work of Chapter 2 has been presented in [37] whereas the full paper hasbeen accepted for publication [74].2. In Chapter 3, we study the EV frequency regulation service under the performance-based compensation scheme. This new compensation scheme has been implementedrecently by the ISOs in U.S. to encourage the usage of fast ramping resources such13Chapter 1. Introductionas batteries, flywheels, and EVs to participate in the market and provide the fre-quency regulation service. However, the effect of the performance-based compensa-tion scheme on the EV frequency regulation service has not yet been studied. Wemodel the revenue of EV frequency regulation service under this new compensationscheme and develop a problem formulation to maximize the revenue. Our formula-tion takes into account the limited battery capacity of the EVs and the uncertainty ofthe AGC signal. An efficient algorithm is developed by solving the formulated prob-lem. We evaluate the performance of the proposed algorithm by simulations basedon real AGC signal and prices from PJM. We compare with a benchmark algorithmwhere the uncertainty of the AGC signal is ignored. The results show that, underthe performance-based compensation scheme, the proposed algorithm can improvethe revenue of the EV frequency regulation service by around 10% compared to thebenchmark algorithm. We have presented part of the work of Chapter 3 in [59] andthe full paper has been accepted for publication [75].3. In Chapter 4, we study the market participation of an aggregator which coordinatesEVs to provide the frequency regulation service to an ISO. We develop a model forthe market participation of the aggregator by taking into account the two-settlementmarket system. We consider two types of DAMs based on the market rules of theCAISO and NYISO. Then, we formulate a problem to optimize the bid in the DAM.We use stochastic programing in the formulation to account for the uncertainty ofprices, along with random arrival and departure of EVs. Our formulation also in-corporates CVaR to model the financial risk caused by uncertainty of the regulationcapacity of the EVs. We develop efficient algorithms to solve the formulated prob-lem. We evaluate the performance of the proposed algorithms using simulations.We compare with an algorithm in the literature and use historical prices of frequency14Chapter 1. Introductionregulation service from CAISO and NYISO to perform the simulations. Results showthat the proposed algorithms improve the profit and reduce the financial risk. Partof the work of this chapter has been presented in [73] and the full paper is submittedfor publication in IEEE Trans. on Smart Grid.1.5 Thesis OrganizationThe rest of the thesis is organized as follows. In Chapter 2, we study the residential DSMfor areas with high penetration of rooftop PV units. We use an external cost to modelthe undesired voltage rise issue caused by the household energy export. We propose anenergy consumption scheduling algorithm for household appliances by jointly consideringTOU prices and the voltage rise issue. In Chapter 3, we study the coordination of the EVcharging to provide frequency regulation service to the ISOs. We develop a problem formu-lation to schedule the hourly regulation capacity of the EVs based on robust optimization.The performance-based compensation scheme is taken into account in the formulation. InChapter 4, we study an aggregator which coordinates EVs to participate in the DAM andRTM of an ISO. We formulate a problem based on stochastic programing to aggregate theuncertain regulation capacity of EVs and determine the bid in the DAM. Finally, this the-sis is concluded in Chapter 5 whereas some potential future works are provided. Chapters2-4 in this thesis are self-contained and included in separate journal or conference papers.The notations are defined separately for Chapters 2-4.15Chapter 2DSM under High Penetration ofRooftop PV Units2.1 IntroductionIn this chapter, we discuss the residential DSM in areas which have high penetration ofrooftop PV units. In those areas, the voltage rise issue is a major obstacle to integrate therooftop PV units in the distribution network and this issue has emerged in several countriessuch as Germany. The load control techniques and ESSs can be used to reduce the reversepower flow from households to substation during high solar radiation hours and mitigatethe voltage rise problem [36, 37, 76–78]. The work in [36] proposes a method to activatethe HVAC loads when the voltage rise is detected. The works in [76–78] focus on usingESSs to mitigate the voltage rise problem. ESSs can absorb the reverse power flow whenthe voltage rise is detected. However, the control of residential load and ESSs requires theapproval of household owners and the economic incentives are not fully addressed in theabove papers.The DSM under the high penetration of rooftop PV units needs to take into accountboth the economic incentives such as TOU pricing and the voltage rise problem. However,the challenge is how to jointly model the voltage rise problem and the energy consumptionscheduling problem. In this chapter, we propose to use the external cost and power flowanalysis to model the effect of the voltage rise issue. The contributions of this chapter are16Chapter 2. DSM under High Penetration of Rooftop PV Unitsas follows:• We propose a residential energy consumption scheduling algorithm for areas withhigh penetration of rooftop PV units. The proposed algorithm aims to reduce theenergy expenses of the users and mitigate the voltage rise problem. It shifts deferrableload from peak consumption hours to hours with high solar power generation.• We introduce a cost function to model the undesired voltage rise. The objectivefunction comprises a monetary cost for household energy consumption, the energygeneration revenue from PV units, and an external cost for the voltage rise. Theexternal cost is modeled with the sensitivity of the voltage magnitudes with respectto the household energy export, which can be obtained from power flow analysis.We formulate a stochastic energy consumption scheduling problem to minimize thecost function. We transform the problem into an LP problem when the feed-intariff does not exceed the energy consumption price and the cost function is convexpiecewise linear. Otherwise, when the cost function is non-convex, an MIP problemis formulated.• We compare our proposed algorithm with the algorithm in [36] on a 33 bus radialdistribution test system. Simulation results show that our proposed algorithm re-duces the energy expenses of the users. The results also show that our algorithm canmitigate the voltage rise problem for areas with high penetration of PV units andreduce the peak-to-average ratio (PAR) of the aggregate load.The rest of the chapter is organized as follows. The system model is introduced inSection 2.2. In Section 2.3, we present the problem formulation and the proposed algorithm.Numerical results are presented in Section 2.4, and a summary is given in Section 2.5.17Chapter 2. DSM under High Penetration of Rooftop PV UnitsSmart Meter(ECS Embedded)Information FlowDistribution Network TOU Pricing and Feed-in TariffUtility CompanyLoad+ -ESSPVPower Flow Analysis ResultsDistribution Network Operator (DNO)( )xbE tFigure 2.1: A schematic diagram of the proposed energy consumption scheduling scheme.The DNO performs the power flow analysis and provides the ECS with information re-garding the voltage rise (e.g., the sensitivity of the voltage with respect to the householdenergy export) at the beginning of the day.2.2 System ModelA block diagram of the proposed energy consumption scheduling scheme for a householdequipped with a PV unit and an ESS is shown in Fig. 2.1. The power flows from thedistribution network to the load when the power generated by the PV unit is insufficientto meet the demand. The direction of the power flow is reversed when the power generatedby the PV unit is higher than the load. Each household is equipped with an ECS [11], whichis embedded in the household smart meter. The ECS retrieves the pricing information fromthe utility company. A DNO performs the power flow analysis [30] for the local distributionnetwork and determines the hours when the households may have abnormally high voltagemagnitude if load shifting is not performed. We assume the DNO performs the power flow18Chapter 2. DSM under High Penetration of Rooftop PV Unitsanalysis based on historical records of demand and generation. The DNO sends the powerflow analysis results to the ECS at the beginning of each day. The ECS determines theoperational schedule of the deferrable load to minimize the electricity bill and reduce thereverse power flow for hours when voltage rise may happen. Examples of such deferrableload include clothes dryers and electric water heaters. Some other loads such as TV andcomputers are must-run load [16] because their energy consumption cannot be shifted. Wealso consider the ESSs in our model, e.g., battery systems owned by the household. Wedenote the sets of deferrable load, must-run load, and ESSs by Sd, Sm, and Ss, respectively.We note that the demand and PV energy generation are stochastic in practice and canbe different from the values obtained from the historical records. This may introduceinaccuracy when the DNO performs the power flow analysis. In Section 2.4, we analyzethe effect of forecast errors via simulations.We denote the set of operating time slots under consideration by T∆= {1, . . . , T}. Theselection of an appropriate length of the time slot depends on the formulated problem. Onthe one hand, a relatively short time slot is helpful to model the dynamics of the load andthe energy generation of PV units. Short time slots make the assumption that the activepower, reactive power, and phase on the buses are static within one time slot more realistic.On the other hand, longer time slots reduce the computational complexity of the algorithm.In this chapter, we choose the length of a time slot to be 15 minutes in order to achievea balance between the performance and the complexity of the algorithm. Furthermore,we use Edi (t), Emj (t), and Esl (t) to denote the energy consumption of the deferrable loadi ∈ Sd, the energy consumption of the must-run load j ∈ Sm, and the charged or dischargedenergy of the ESS l ∈ Ss at time slot t ∈ T , respectively. The energy generated from PVunits at time slot t is denoted by Eg(t). For a household connected to bus b, the exportedenergy at time slot t is denoted by Exb (t) and is given by19Chapter 2. DSM under High Penetration of Rooftop PV UnitsExb (t) = Eg(t)−∑i∈SdEdi (t)−∑j∈SmEmj (t)−∑l∈SsEsl (t). (2.1)The ECS adjusts Exb (t) by controlling the operation of deferrable loads and ESSs. Let etrepresent the vector of the energy consumption of the deferrable loads and the charged ordischarged energy of the ESSs at time slot t. et can be written aset =(Ed1 (t), . . . , Ed|Sd|(t), Es1(t), . . . , Es|Ss|(t)), (2.2)where |Sd| and |Ss| are the cardinalities of sets Sd and Ss, respectively. The rooftop PVunits are required to operate with unity power factor according to the Institute of Electricaland Electronics Engineers (IEEE) 1547 standard [79]. We restrict our model to considerthe active power generation (i.e., Eg(t)) of PV units.2.2.1 Voltage Rise ProblemThe voltage rise problem occurs when many households in the distribution network havelarge positive Exb (t). A 33 bus radial distribution network [80] is shown in Fig. 2.2.The households are connected to the buses. We denote the set of buses as B (e.g., B ={1, . . . , 33} in Fig. 2.2). We denote the voltage magnitude on bus n ∈ B as Vn(t) > 0.Note that Vn(t) is the root mean square of the voltage and t in Vn(t) denotes a time slot.We use b to specify the bus to which the considered household is connected, whereas busn is an arbitrary bus. The value of Vn(t) is determined using the power flow analysis. Weassume the active power, reactive power, voltage magnitude, and phase are static withinone time slot. We first consider the case when load shifting is not performed. The DNOexecutes the power flow analysis using historical records of demand and generation. Inthis case, Vn(t) may exceed an upper limit (e.g., 1.05 pu in North America) on some of the20Chapter 2. DSM under High Penetration of Rooftop PV Units12 3 87654 9 10 11 1615141312 1817212019 22323130 33282726 292423 25SubstationFigure 2.2: One line diagram of a 33 bus radial distribution test system. The householdsare connected to the buses. The voltage rise problem occurs when many households havesubstantial positive energy export. The DSM program encourages households to shift loadsto hours with high solar radiation to mitigate the voltage rise problem.buses. Hence, we obtain a set of buses with abnormally high voltage magnitude in timeslot t, which is denoted by Bv(t). We aim to decrease the voltage magnitude Vn(t) of busesn ∈ Bv(t) by decreasing Exb (t) via load shifting, when Bv(t) is not ∅.The sensitivity of the voltage magnitude Vn(t) of an arbitrary bus n with respect tohousehold energy export Exb (t) is characterized by the partial derivative∂Vn(t)∂Exb(t). We notethat Vn(t) depends on the power flow of all the buses. Moreover, a non-convex problem hasto be solved to obtain Vn(t). However, in order to reduce the computational complexity, alinear model of the external cost is adopted, and the partial derivatives ∂Vn(t)∂Exb(t)are used toanalyze the sensitivity of Vn(t) with respect to Exb (t).∂Vn(t)∂Exb(t)is obtained from power flowanalysis using the Newton-Raphson method [30]. That is, the Newton-Raphson methodis used to obtain the Jacobian matrix which contains the information about the partialderivatives of the active and reactive powers with respect to the voltage magnitude and21Chapter 2. DSM under High Penetration of Rooftop PV Unitsphase. Next, we calculate the inverse of the Jacobian matrix [81, 82], which containsthe partial derivatives of the voltage magnitude and phase with respect to the active andreactive powers. We use the partial derivative ∂Vn(t)∂Exb(t)obtained from the inverse of Jacobianmatrix to model the contribution of the energy import to the voltage rise. The exact valuesof ∂Vn(t)∂Exb(t)can be calculated using off-the-shelf software such as PowerWorld [83] or Matlab.We assume the DNO performs the power flow analysis based on the historical records ofdemand and generation. The DNO sends the results, i.e., ∂Vn(t)∂Exb(t)and Vn(t), to the ECS asthe input for the energy consumption scheduling algorithm.Note that in practice, it is possible that some utility companies have made investmentand upgraded the transformer with on load tap charger transformer (OLTC). In that case,the DSM program should be used along with the OLTC to mitigate the voltage rise. Werestrict our paper to focus on the DSM.2.2.2 Cost FunctionThe energy export of a household on an arbitrary bus b may contribute to the voltage riseproblem of other households in the distribution network (including the households on otherbuses). We model the contribution of the household energy export to the abnormally highvoltage magnitude of other households as an external cost. In economics, an externalityarises when a person engages in an activity that influences the well-being of a bystander[84]. We assume the household aims to minimize its own electricity bill and its externalcost.To tackle the voltage rise problem, the DNO can encourage households to reduce theenergy export during high solar radiation hours. We introduce a parameter h and assumean external cost occurs when the energy export is higher than h. The DNO determinesthe value of h by running the power flow analysis in an iterative manner. The power flow22Chapter 2. DSM under High Penetration of Rooftop PV Unitsanalysis is used to determine the bus voltage magnitude. In the power flow analysis, theDNO first sets the energy export from each household to be zero and runs the power flowanalysis. In this case, the voltage magnitude on each bus is lower than 1.05 pu. Then, theDNO gradually increases the energy export of all households simultaneously and repeatsthe power flow analysis to obtain the voltage magnitude for different energy exports. Asthe energy export increases, the voltage magnitude obtained in the power flow analysisincreases as well. When the voltage magnitude of any bus reaches 1.05 pu, the DNOstops increasing the energy export and the obtained value of the energy export can beused as the value of parameter h. The value of h depends on the voltage magnitude onthe secondary winding of the substation transformer, the line impedance, the line lengthbetween adjacent buses, and the number of households connected to each bus. Thesevalues are used as inputs in the power flow analysis to calculate the voltage magnitude.The value of h should be determined by the DNO for each distribution network accordingto its specific parameters.Note that there could be different methods to determine h while ensuring fairnessregarding the household energy export and voltage rise issue. For example, we may setdifferent h for households according to the size of their PV units. More importantly, thecontribution of the household energy export to the voltage rise is highly dependent on theline length between the substation and the household. Hence, we may design strategies toset h based on the line length. In this thesis, for simplicity, we assume the same parameterh to all households.Consider an arbitrary household connected to bus b ∈ B, we model the external cost ofthe household energy export ascx(et) = η∑n∈Bv(t)∂Vn(t)∂Exb (t)[Exb (t)− h]+, (2.3)23Chapter 2. DSM under High Penetration of Rooftop PV Unitswhere [x]+ = max{x, 0}.∑n∈Bv(t)∂Vn(t)∂Exb(t)is the contribution of one unit household energyexport to the abnormally high voltage magnitude of the buses n ∈ Bv(t). The voltagerise happens on some of the buses, i.e., the buses in set Bv(t). The parameter η is a non-negative coefficient. ∂Vn(t)∂Exb(t)is obtained using a sensitivity-based linear model of the voltagerise. The first order derivative ∂Vn(t)∂Exb(t)models the sensitivity of the voltage magnitude withrespect to the power flows. The sensitivity-based linear model has been widely used tostudy the voltage rise problem [81, 82, 85, 86]. The accuracy of the model may degradewhen power flows and voltage deviate from their expected values significantly.Assume the ECS aims to reduce the monetary cost (i.e., the electricity bill) and theexternal cost. The household purchases energy at price pe(t) when Exb (t) < 0 and sellsenergy with a feed-in tariff ps(t) when Exb (t) ≥ 0. The total cost at time slot t is given byc(et) =η∑n∈Bv(t)∂Vn(t)∂Exb (t)(Exb (t)− h)− ps(t)Exb (t), if Exb (t) ≥ h,−ps(t)Exb (t), if 0 ≤ Exb (t) < h,pe(t)|Exb (t)|, if Exb (t) < 0,(2.4)where |Exb (t)| represents the absolute value of Exb (t). Parameter η can be tuned to adjustthe weight between the external cost and the monetary cost. By increasing η, the ECStends to shift more load to the high solar radiation hours to reduce the external cost ofthe voltage rise problem. By decreasing η, the ECS puts more weight on minimizing theelectricity bill. When η = 0, the ECS only aims at minimizing the bill and neglects thepotential voltage rise problem. Fig. 2.3 shows that cost function c(et) is a convex piecewiselinear function [87] when ps(t) ≤ pe(t). However, it becomes a non-convex piecewise linearfunction when ps(t) > pe(t). Both cases will be considered in Section 2.3.The households schedule their load in a distributed manner. Distributed load control al-gorithms are preferred as they reduce the computational complexity and offer more privacy24Chapter 2. DSM under High Penetration of Rooftop PV UnitsCostenergy exportSlope isSlope isQuota hSlope is( ) ( )s ep t p t ( )sp t!( )ep t!( )tc e( )xbE t( )( )( )( )vsnxn t bV tp tE t"#$!$% Figure 2.3: Piecewise linear cost function for a household with PV unit. We considerthe case ps(t) ≤ pe(t) in this figure, i.e., the feed-in tariff does not exceed the energyconsumption price. On the other hand, Some utility companies promote the rooftop PVinstallation aggressively and have ps(t) > pe(t).compared to centralized algorithms. We investigate the energy consumption schedulingproblem for the ECS by an arbitrary household.We introduce the external cost in order to model the effect of the voltage rise issue.One promising method to enforce the households to consider the external cost is the mon-etization of the external cost [88], which has been used in several areas [89–92], such as thecarbon emission tax in the energy generation. This concept can be used for those areaswith high rooftop PV penetration. The external cost can be incurred to the user by meansof a carefully selected price function. For example, the utility company may collect anadditional charge from households with an energy export exceeding the threshold h. Someutility companies (e.g., Georgia Power) have already been collecting a monthly charge fromhouseholds which connect their rooftop PV units to the distribution network. With AMI,25Chapter 2. DSM under High Penetration of Rooftop PV Unitsthe utility company can monitor the household energy export in short time intervals anddetermine the external cost according to the energy export.2.3 Energy Consumption Scheduling with UncertainPV Power GenerationIn this section, we consider the energy consumption scheduling problem for minimizing theexpected payment and external cost of the user. We formulate the problem as a stochasticprogramming problem [93]. The uncertainty in the power generation from PV units istaken into account in our formulation. The PV power generation Eg(t) is intermittentand random in nature [94–96]. The ECS needs to determine et at the current time slot t,before the future PV power generation Eg(t + 1), . . . , Eg(T ) can be observed. In the nexttime slot, t+1, the power consumption vector et+1 will be updated as the new informationabout the PV power generation is received.For the rest of this chapter, we denote the current time slot by t, whereas τ representsan arbitrary time slot. We define Tt∆= {t, . . . , T} as the set of time slots from the currenttime slot t onwards.2.3.1 Nested Stochastic FormulationFor current time slot t, we define Eri (t) as the amount of energy required to finish theoperation of appliance i ∈ Sd. The energy requirement Eri (t) is updated asEri (t + 1) = max{Eri (t)− Edi (t), 0}, i ∈ Sd. (2.5)We denote the maximum energy consumption of the ith deferrable load within one time26Chapter 2. DSM under High Penetration of Rooftop PV Unitsslot by Emaxi,d . Moreover, Emaxl,s ≥ 0 represents the maximum energy that can be stored inESS l during one time slot. Eminl,s ≤ 0 is the maximum energy that can be derived fromESS l during one time slot. We denote the SOC of ESS l at time slot t by sl(t). Thebattery capacity of ESS l is denoted by bl. For ESS l, its SOC is updated assl(t) = sl(t− 1) +Esl (t)bl, l ∈ Ss. (2.6)At current time slot t, the ECS optimizes et for minimization of the cost from thecurrent time slot t onwards. The nested form of our stochastic problem isminimizeet∈χtc(et) + E[infet+1∈χt+1c(et+1) + · · ·+ E[infeT∈χTc(eT )]](2.7a)subject toT di∑τ=tEdi (τ) = Eri (t), i ∈ Sd, (2.7b)where T di is the deadline by which the operation of deferrable load i ∈ Sd has to be finished.E represents the expectation with respect to the uncertain PV power generation. The must-run load may also have uncertainty, depending on how accurately we can forecast themust-run load profile. We also assume that the household users determine the operationalconstraints of each appliance. The feasible set, χτ , τ ∈ Tt, in (2.7a) can be written asχτ = { eτ | 0 ≤ Edi (τ) ≤ Emaxi,d , ∀ i ∈ Sd,Eminl,s ≤ Esl (τ) ≤ Emaxl,s , ∀ l ∈ Ss,Smin ≤ sl(τ) ≤ 1, ∀ l ∈ Ss } . (2.8)The first two constraints in (2.8) ensure that the energy consumption of load i and chargedor discharged energy of ESS l at hour τ do not exceed their respective maximal values.27Chapter 2. DSM under High Penetration of Rooftop PV UnitsThe third constraint prevents the ESS from discharging to an overly low SOC, where Sminis the minimum SOC of the ESS, e.g., Smin = 0.3. Discharging an ESS to an overly lowSOC may degrade the life cycle of the ESS.Problem (2.7) is difficult to solve directly. Hence, we adopt the sample average approx-imation (SAA) technique [93] to solve problem (2.7).2.3.2 Sample Average ApproximationThe SAA technique generates scenarios for unknown future PV power generation andevaluates the objective function by averaging over different scenarios. Suppose we haveK scenarios {ω1, . . . ,ωK} based on the historical records of the PV power generation.Each scenario is a possible realization of (Eg(t + 1), . . . , Eg(T )). For the kth scenarioωk, the hourly PV power generation is denoted as (Eg(t + 1,ωk), . . . , Eg(T,ωk)), whereEg(τ,ωk), τ ∈ Tt+1, is the value of Eg(τ) under scenario ωk. We define K∆= {1, . . . , K} asthe set of scenarios. For an arbitrary scenario k ∈ K, let Edi (τ,ωk), Emj (τ,ωk), Esl (τ,ωk),Exb (τ,ωk), sl(τ,ωk), and e(τ,ωk) denote the energy consumption of the ith deferrable load,the energy consumption of the jth must-run load, the charged or discharged energy of thelth ESS, the household energy export, the SOC of the lth ESS, and the energy consumptionvector at time slot τ under scenario ωk, respectively.We estimate the expected cost by averaging the cost over different scenarios ωk, ∀k ∈ K.By using the SAA technique, problem (2.7) becomesminimizee(τ,ωk), τ ∈ Tt, k ∈ K∑k∈KP(ωk)∑τ∈Ttc(e(τ,ωk)) (2.9a)subject toT di∑τ=tEdi (τ,ωk) = Eri (t), i ∈ Sd, k ∈ K, (2.9b)0 ≤ Edi (τ,ωk) ≤ Emaxi,d , i ∈ Sd, τ ∈ Tt, k ∈ K, (2.9c)28Chapter 2. DSM under High Penetration of Rooftop PV UnitsEminl,s ≤ Esl (τ,ωk) ≤ Emaxl,s , l ∈ Ss, τ ∈ Tt, k ∈ K, (2.9d)Smin ≤ sl(τ,ωk) ≤ 1, l ∈ Ss, τ ∈ Tt, k ∈ K, (2.9e)e(t,ωk) = e(t,ωl), k, l ∈ K. (2.9f)The term e(τ,ωk) is the energy consumption vector of controllable loads and ESSs at timeslot τ under scenario ωk. c(e(τ,ωk)) denotes the cost function in time slot τ under scenarioωk.∑τ∈Ttc(e(τ,ωk)) in (2.9a) indicates the total cost under scenario ωk from time slot tto time slot T . The objective function (2.9a) is the cost from time slot t onwards, averagedacross the K scenarios. Constraints (2.9b)–(2.9e) are the extensions of the constraints in(2.7b) and (2.8) under scenario ωk. Constraint (2.9f) is the non-anticipativity constraintof stochastic programming [93]. In constraint (2.9f), t is the current time slot while τin problem (2.9) denotes an arbitrary time slot. Constraint (2.9f) reflects the fact thatthe ECS is required to make a deterministic decision in the current time slot before theunknown parameters are revealed. Constraint (2.9f) enforces the variables of the currenttime slot t to be equal under different scenarios, such that the obtained solution in thecurrent time slot t is deterministic, i.e., does not depend on a specific scenario. Problem(2.9) has an objective function which can be either convex or non-convex piecewise linear.For both the convex and non-convex cases, problem (2.9) will be studied in the followingsubsections.It should be noted that there are different approaches in the scenario-based analysis.In this chapter, we consider the case where we can minimize the expected cost usingSAA technique. On the other hand, in some cases, we are actually more interested in aparticular subset of the scenarios, e.g., the worst cases. Some other techniques such asrobust optimization or CVaR can be used for the analysis of the worst cases. We willexplore those techniques of the scenario-based analysis in the following chapters.29Chapter 2. DSM under High Penetration of Rooftop PV UnitsMoreover, we clarify that the proposed scheme actually has two steps. In the first step,the DNO calculates the sensitivity of the voltage rise with respect to the household energyexport in the power flow analysis and then sends the results to households. In the secondstep, the households schedule the energy consumption of their appliances based on thereceived information.2.3.3 Convex Piecewise Linear Objective FunctionFirst, we consider the case for which ps(τ) ≤ pe(τ), τ ∈ Tt. Fig. 2.3 illustrates the costfunction for this case. The cost function c(e(τ,ωk)) is a convex piecewise linear functionsince c(e(τ,ωk)) in (2.4) can be written asc(e(τ,ωk)) = max{− ps(τ)h +((η∑n∈Bv(τ)∂Vn(τ)∂Exb (τ)− ps(τ))(Exb (τ,ωk)− h)),− ps(τ)Exb (τ,ωk), pe(τ)(−Exb (τ,ωk))}. (2.10)By introducing auxiliary variables u(τ,ωk), τ ∈ Tt, k ∈ K, problem (2.9) can be rewrit-ten asminimizee(τ,ωk), u(τ,ωk),τ ∈ Tt, k ∈ K∑k∈KP(ωk)∑τ∈Ttu(τ,ωk) (2.11a)subject to u(τ,ωk) ≥ −ps(τ)h +((η∑n∈Bv(τ)∂Vn(τ)∂Exb (τ)− ps(τ))(Exb (τ,ωk)− h)), τ ∈ Tt, k ∈ K, (2.11b)u(τ,ωk) ≥ −ps(τ)Exb (τ,ωk), τ ∈ Tt, k ∈ K, (2.11c)u(τ,ωk) ≥ pe(τ)(−Exb (τ,ωk)), τ ∈ Tt, k ∈ K, (2.11d)constraints (2.9b)–(2.9e).30Chapter 2. DSM under High Penetration of Rooftop PV UnitsConstraints (2.11b)–(2.11d) are introduced to replace the piecewise linear objective func-tion c(e(τ,ωk)) in (2.10). Problem (2.11) can be solved efficiently using LP.2.3.4 Non-convex Piecewise Linear Objective FunctionFor the case in which ps(τ) > pe(τ), τ ∈ Tt, the cost function is a non-convex piecewiselinear function, cf. Fig. 2.4. The high feed-in tariff is typically used to encourage householdsto install rooftop PV units. For example, in the province of Ontario in Canada, the feed-intariff for rooftop PV units is 32.9 − 39.6 cents per kWh [24] while the residential energyconsumption price is 7.2− 12.9 cents per kWh. Moreover, in the state of Georgia in U.S.,the utility company has a residential energy retailing price of 4.8-9.6 cents per kWh whilepurchases energy from rooftop PV units at 17 cents per kWh [97]. The high feed-in tariffin the province of Ontario in Canada and Georgia state in U.S. are for the roof-top PVunits only.A non-convex piecewise linear objective function cannot be directly transformed intoan LP. To tackle problem (2.9) with a non-convex piecewise linear objective function, weintroduce new variables q0(τ,ωk), q1(τ,ωk), and q2(τ,ωk) to represent the components ofenergy export Exb (τ,ωk) in the three different pricing ranges, cf. Fig. 2.4. Then, we haveq0(τ,ωk) = min{ Exb (τ,ωk)−H, |H|}, (2.12)q1(τ,ωk) = min{max{Exb (τ,ωk), 0 }, h}, (2.13)q2(τ,ωk) = max{ Exb (τ,ωk)− h, 0 }, (2.14)Exb (τ,ωk) = H + q0(τ,ωk) + q1(τ,ωk) + q2(τ,ωk), (2.15)31Chapter 2. DSM under High Penetration of Rooftop PV UnitsCost Energy ExportQuotaConstanthHSlope is Slope isSlope is ! !, kc " #e !0 , kq " # !1 , kq " # !2 , kq " #( )sp "$( )ep "$( ) ( )s ep p" "%( , )xb kE " #( )( )( )( )vsnxbnVpE""& ""'($() Figure 2.4: Non-convex piecewise linear cost function. This case occurs when the feed-intariff of the utility company is higher than the energy consumption price (i.e., ps(τ) >pe(τ)).c(e(τ,ωk)) = pe(τ)(|H| − q0(τ,ωk))− ps(τ)q1(τ,ωk)−(ps(τ)− η∑n∈Bv(τ)∂Vn(τ)∂Exb (τ))q2(τ,ωk), (2.16)where H is a negative constant, cf. Fig. 2.4. The value of H is selected to be less than theminimum of Exb (τ,ωk). Now, problem (2.9) can be rewritten asminimizee(τ,ωk), q0(τ,ωk),q1(τ,ωk), q2(τ,ωk),τ ∈ Tt, k ∈ K∑k∈KP(ωk)∑τ∈Tt(pe(τ)(|H| − q0(τ,ωk))−ps(τ)q1(τ,ωk)−(ps(τ)− η∑n∈Bv(τ)∂Vn(τ)∂Exb (τ))q2(τ,ωk)) (2.17a)subject to 0 ≤ q0(τ,ωk) ≤ |H|, τ ∈ Tt, k ∈ K, (2.17b)0 ≤ q1(τ,ωk) ≤ h, τ ∈ Tt, k ∈ K, (2.17c)0 ≤ q2(τ,ωk), τ ∈ Tt, k ∈ K, (2.17d)q1(τ,ωk)(q0(τ,ωk)− |H|) = 0, τ ∈ Tt, k ∈ K, (2.17e)32Chapter 2. DSM under High Penetration of Rooftop PV Unitsq2(τ,ωk)(q1(τ,ωk)− h) = 0, τ ∈ Tt, k ∈ K, (2.17f)∑i∈SdEdi (τ,ωk) +∑j∈SmEmj (τ,ωk) +∑l∈SsEsl (τ,ωk) ≥ 0,τ ∈ Tt, k ∈ K, (2.17g)constraints (2.9b)–(2.9e), and (2.15).Constraints (2.17b)–(2.17d) ensure that q0(τ,ωk), q1(τ,ωk), and q2(τ,ωk) are selectedwithin their ranges. Because of constraint (2.17e), q1(τ,ωk) can have a non-zero value onlyif q0(τ,ωk) has reached its maximum. Similarly, because of constraint (2.17f), q2(τ,ωk) canbe non-zero only if q1(τ,ωk) has reached its maximum. Constraint (2.17g) ensures that,under high feed-in tariff, the households export energy from their PV units rather thantheir ESSs. Both Ontario Hydro and Georgia Power [24, 97] restrict that only the energyexport from certain types of renewable resources are eligible for high feed-in tariff. Other-wise, some households may misuse the high feed-in tariff by storing energy from a utilitycompany and sell the energy later. Under the current net metering programs of OntarioHydro and Georgia Power, the discharged energy from an ESS can be used to support theload but cannot be compensated if used to export energy. Some utility companies (e.g.,Georgia Power) require the households to install an additional meter for their PV units,such that the utility can determine whether the energy export comes from the PV unitsor not.The problem defined in (2.17) contains constraints (2.17e) and (2.17f) which are neitherlinear nor convex. To tackle this problem, we introduce new auxiliary binary variablesb1(τ,ωk) and b2(τ,ωk) to indicate whether q1(τ,ωk) and q2(τ,ωk) are strictly positive ornot, respectively. For example, we have b1(τ,ωk) = 1 if q1(τ,ωk) > 0, and b1(τ,ωk) = 0otherwise. Next, we transform problem (2.17) as33Chapter 2. DSM under High Penetration of Rooftop PV Unitsminimizee(τ,ωk), q0(τ,ωk),q1(τ,ωk), q2(τ,ωk),b1(τ,ωk), b2(τ,ωk),τ ∈ Tt, k ∈ K∑k∈KP(ωk)∑τ∈Tt(pe(τ)(|H| − q0(τ,ωk))−ps(τ)q1(τ,ωk)−(ps(τ)− η∑n∈Bv(τ)∂Vn(τ)∂Exb (τ))q2(τ,ωk)) (2.18a)subject to b1(τ,ωk) ≥q1(τ,ωk)h, τ ∈ Tt, k ∈ K, (2.18b)b2(τ,ωk) ≥q2(τ,ωk)G, τ ∈ Tt, k ∈ K, (2.18c)b1(τ,ωk) ≤q0(τ,ωk)|H|, τ ∈ Tt, k ∈ K, (2.18d)b2(τ,ωk) ≤q1(τ,ωk)h, τ ∈ Tt, k ∈ K, (2.18e)b1(τ,ωk), b2(τ,ωk) ∈ {0, 1}, τ ∈ Tt, k ∈ K, (2.18f)constraints (2.9b)–(2.9e), (2.15), (2.17b)–(2.17d), and (2.17g),where G is a positive constant and larger than the maximum of q2(τ,ωk). Constraints(2.18b), (2.18c), and (2.18f) yield b1(τ,ωk) = 1 and b2(τ,ωk) = 1 if q1(τ,ωk) and q2(τ,ωk)are strictly positive, respectively. Constraints (2.18b) and (2.18d) replace nonlinear con-straint (2.17e) in problem (2.17). They guarantee that q1(τ,ωk) is strictly positive onlyif q0(τ,ωk) is at its maximum. If q0(τ,ωk) < |H|, we have b1(τ,ωk) = 0 according toconstraints (2.18d) and (2.18f). Then, we can infer q1(τ,ωk) = 0 according to (2.18b).Similarly, constraints (2.18c) and (2.18e) replace nonlinear constraint (2.17f). They ensurethat q2(τ,ωk) is strictly positive only if q1(τ,ωk) is at its maximum. Problem (2.18) isan MIP problem and can be solved with an off-the-shelf optimization software such asMOSEK [98] or CPLEX [99]. We will discuss its computational complexity in Section 2.4.We present the proposed energy consumption scheduling algorithm in Algorithm 2.1.The algorithm is executed at the beginning of each hour t by the ECS. The ECS first34Chapter 2. DSM under High Penetration of Rooftop PV UnitsAlgorithm 2.1 Energy consumption scheduling algorithm executed at the beginning ofeach time slot t by the ECS1: Initialize Sd, Ss, Eri (t), Tdi , i ∈ Sd, Tt, and K.2: Retrieve pe(τ), ps(τ), τ ∈ Tt from the utility company through the communicationnetwork.3: Retrieve η, h, Vn(τ), and∑n∈Bv(τ)∂Vn(τ)∂Exb (τ)from the DNO.4: Construct scenarios ωk, k ∈ K, according to the historical records of power generationfrom the PV units.5: if ps(τ) ≤ pe(τ), τ ∈ Tt then6: Solve problem (2.11) to obtain the scheduling results.7: else8: Solve problem (2.18) to obtain the scheduling results.9: end if10: Control the operation of the deferrable load and ESSs according to the schedulingresults.11: Update Eri (t+ 1), sl(t + 1) according to (2.5) and (2.6), respectively.initializes the parameters and retrieves the pricing information and power flow analysisinformation from the utility company and the DNO, respectively (Lines 1-3). Then, theECS constructs scenarios ωk, k ∈ K, according to historical data records of PV powergeneration (Line 4). Subsequently, problem (2.11) is solved when the objective functionis a convex piecewise linear function, otherwise problem (2.18) is solved (Lines 5-9). Theresults are used to control the operation of the deferrable load and ESSs (Line 10). Finally,the ECS updates Eri (t + 1) and sl(t + 1) to prepare for the scheduling for the next hour(Line 11).2.4 Performance EvaluationIn this section, we present simulation results and assess the performance of the proposedenergy consumption scheduling algorithm. We compare our proposed algorithm with thealgorithm in [36], which we refer to as the HDA (HVAC load Direct Activation) algorithm.35Chapter 2. DSM under High Penetration of Rooftop PV UnitsTable 2.1: List of deferrable loadsLoad Energy Usage (kWh) Task DeadlineDish washer [0, 4] [6 pm, 8 pm]Water heater [0, 10] [7 pm, 12 am]Clothes dryer [0, 3] [10 pm, 12 am]While our proposed algorithm schedules household loads based on a cost function, theHDA algorithm is based on a real-time voltage signal. To be specific, the HDA algorithmassume that the the voltage magnitude is monitored and then households can activatethe idle HVAC load (e.g., water heater) when a voltage rise is detected. However, HDAdoes not take into account the TOU prices. We compare with the HDA algorithm sinceit considers how to use loads to handle the voltage rise problem, which shares the sameresearch interest of the work in this chapter. We also compare our proposed algorithmwith a benchmark algorithm, where the ECS is not deployed and the deferrable load is notcontrolled.We consider the 33 bus radial distribution test system [80]. We assume there are 10-40 households on each bus and each household has installed PV units and ESSs. Thehouseholds have deferrable loads as shown in Table 2.1. The daily consumption of themust-run load for each household is randomly selected from the interval [6, 18] kWh andthe hourly consumption profile is generated according to [100]. The installed capacity ofthe PV units for each household is randomly selected from the interval [1.2, 3.6] kW. Wedenote the percentage of deferrable load by β. We haveβ =∑τ∈T∑i∈Sd Edi (τ)∑τ∈T(∑i∈Sd Edi (τ) +∑j∈Sm Emj (τ)) × 100%.We consider the case where the household is equipped with a small-scale local battery36Chapter 2. DSM under High Penetration of Rooftop PV Unitssystem. The capacity of the battery is 4 kWh [78]. The maximum charged and dischargedenergy in one hour is 1 kWh and −1 kWh, respectively.We use two pricing schemes in order to analyze the two cases presented in Sections2.3.3 and 2.3.4. We denote the first pricing scheme by Price1, where we have ps(τ) = pe(τ).Price1 is similar to the TOU net metering program of Southern California Edison [101].The second pricing scheme is denoted by Price2, where the energy selling price is obtainedfrom [24]. Price2 represents the case when the utility company sets the feed-in tariff higherthan the price of energy consumption, as Ontario Hydro does.Our simulations cover a period of 30 days. The number of sunny, cloudy, and rainy daysis 23, 4, and 3, respectively. Different PV units may experience different weather condi-tions, especially when the distribution network covers a large geographical area. Moreover,because of the elevation angle of the sun, even on sunny days, some PV units may not beable to produce energy at their full capacity. Hence, for sunny days, we assume 85% of thePV units generate energy with their full capacity and the remaining 15% of the PV unitsproduce power at 60% of the capacity in each time slot. The hourly profile of PV powergeneration is obtained from [102].We first perform the power flow analysis with the Newton-Raphson method and obtainthe partial derivatives of the voltage magnitudes with respect to household energy export[30, 81, 82]. We assume η = 2 × 10−4 in (2.3) for simulation purpose. The value of η isselected such that the external cost is comparable to the monetary cost. We consider anarbitrary household on bus b = 18 in the 33 bus radial distribution test system (cf. Fig.2.2). The line impedances between buses are obtained from [80].37Chapter 2. DSM under High Penetration of Rooftop PV Units1000 2000 3000 4000 5000Number of variables01020304050Running time (seconds)LP, Price12000 3000 4000 5000 6000Number of variables01020304050Running time (seconds)MIP, Price2(a) Running time of LP solver (b) Running time of MIP solverFigure 2.5: Running time of Algorithm 2.1 with LP and MIP solvers under Price1 andPrice2, respectively.2.4.1 Computational ComplexityWe investigate the computational complexity of the proposed algorithm in Fig. 2.5. Weuse the optimization software MOSEK to solve LP problem (2.11) and MIP problem (2.18).The run time is evaluated with a desktop computer which has a quad-core processor and 16gigabyte memory. For an MIP problem, MOSEK uses the branch-and-cut method and isable to find an ǫ-optimal solution within a reasonable time. Here, ǫ denotes the optimalitygap (i.e., the gap between the obtained objective value and the optimal objective valuedoes not exceed a small value ǫ). The run time was measured for ǫ = 0.01.Fig. 2.5 shows the running time of Algorithm 2.1 under Price1 and Price2, respectively.LP problem (2.11) and MIP problem (2.18) are solved respectively under Price1 and Price238Chapter 2. DSM under High Penetration of Rooftop PV Unitsin Algorithm 2.1. We consider the problem formulation with different scales and simulatethe running time of the proposed algorithm versus the number of variables. The resultsare shown in Fig. 2.5. For the LP problem, the computational complexity increases inpolynomial time with respect to the number of variables. On the other hand, for the MIPproblem, the computational complexity increases exponentially with respect to the numberof variables. Hence, if the scale of MIP problem (2.18) further increases, the branch andcut method may not be applicable. In that case, problem (2.18) should be solved usingapproximation methods such as linear relaxation.We note that, in practice, a smart meter typically has much less computational re-sources compared to a desktop computer. For example, the Atmel smart meter [103] usesan ARM Cortex central processing unit (CPU) with a bus frequency of 120 MHz and thesize of the memory is 2 MB. Hence, for current smart meters, it is difficult to implementthe load scheduling algorithms which require to solve complicated optimization problems.On the other hand, we expect that future smart meters will have more computationalresources, considering that the prices of CPU and memory are constantly dropping. More-over, new computational infrastructure (e.g., cloud computing) can be leveraged to tacklethis problem. For an ECS, its task can be divided into three steps. First, the ECS needsto retrieve the values of required energy and task deadline from each appliance. Thesevalues are used as inputs of the optimization problem for load scheduling. Second, theoptimization problem is solved to determine the operating schedule of different appliances.Third, the ECS uses the results to control the appliances. A smart meter seems to be anideal candidate to perform the first and the third steps. However, it may not be necessaryfor the smart meter to perform the second step. Instead, the smart meter can offload thecomputational tasks to one or multiple servers in a cloud platform (e.g., Microsoft Azure[104]), which has abundant computational resources. After the servers have solved the op-39Chapter 2. DSM under High Penetration of Rooftop PV Unitstimization problem, the results can be sent to the ECS via a communication infrastructure.2.4.2 Performance Gains for Voltage Rise MitigationWe evaluate the performance of our proposed algorithm and the HDA algorithm [36] fortackling the voltage rise problem. We use the Newton-Raphson method to calculate thebus voltage magnitude for each time slot after load shifting. Fig. 2.6 shows the voltagemagnitude profile of bus 18 in Fig. 2.2. Bus 18 is selected as it has a longer distance fromthe substation compared to other buses and experiences a larger voltage variation (voltagedrop and voltage rise) when the power flow changes. As shown in Fig. 2.6, our proposedalgorithm avoids the voltage rise problem for households during the noon period (10:00 -14:00). The proposed algorithm outperforms the HDA algorithm for two reasons. First,we jointly considered the residential load and ESSs in our model. Second, for the HDAalgorithm, the tasks of some deferrable loads (e.g., dish washers, clothes dryers) may havealready been finished, when voltage rise happens. These loads may not be activated in thiscase. Note we consider Price2 in Fig. 2.6 but the proposed algorithm achieves a similarperformance in terms of voltage rise mitigation for Price1.2.4.3 Benefits Regarding Electricity Bill and PARFig. 2.7 shows the monthly electricity bill for different amount of daily rooftop PV powergeneration. When η = 2 × 10−4, the proposed algorithm reduces the user payment from$11.6 to $7.5 when the daily generation is 16 kWh under Price1. The user paymentdecreases from $1.3 to −$3.2 under Price2. For the proposed algorithm, we consider twocases, namely η = 2 × 10−4 and η = 0. The value of η can be tuned in the proposedalgorithm to adjust the weight on reducing the voltage rise magnitude and reducing theelectricity bill. When η = 0, the ECS neglects the potential voltage rise problem and only40Chapter 2. DSM under High Penetration of Rooftop PV Units2 4 6 8 10 12 14 16 18 20 22 24Time (hour)0.9811.021.041.061.081.1Voltage Magnitude (pu)Without ECS DeploymentHDA AlgorithmProposed AlgorithmVoltage Rise AvoidedVoltage LimitFigure 2.6: Comparison of the household voltage magnitude profile of the proposed algo-rithm and the HDA algorithm. The voltage magnitude is shown in pu, whereas the basevalue of pu is the standard voltage (e.g., 120 volt for households). The voltage limit is1.05 pu. In the proposed algorithm, the voltage rise is avoided by reducing the householdenergy export during hours with high solar power generation.aims to reduce the electricity bill. As shown in Fig. 2.7, the bill is reduced when η = 0,compared to the case when η = 2 × 10−4. In particular, the bill is decreased from $7.5to $6.1 when η is decreased from 2 × 10−4 to 0 under Price1. Under Price2, the bill isdecreased from −$3.2 to −$4.9. The saving of load shifting depends on the pricing schemeused by the utility company. For example, the Southern California Edison has recentlyintroduced a residential TOU pricing scheme where the price during peak hours is 4.18times higher than the price during the night hours [4]. The load shifting can achieve ahigher saving on electricity bill in this case. Moreover, if ECSs are widely installed, the41Chapter 2. DSM under High Penetration of Rooftop PV Units12 14 16 18 20Daily Rooftop PV Power Generation (kWh)02468101214161820UserPayment(DollarperMonth)Without ECS DeploymentHDA AlgorithmProposed Algorithm, η = 2× 10−4Proposed Algorithm, η = 012 14 16 18 20Daily Rooftop PV Power Generation (kWh)-15-10-5051015UserPayment(DollarperMonth)Without ECS DeploymentHDA AlgorithmProposed Algorithm, η = 2× 10−4Proposed Algorithm, η = 0(a) User Payment under Price1 (b) User Payment under Price2Figure 2.7: Comparison of the monthly electricity bill for different amount of daily rooftopPV power generation.total saving can be significant, considering the large numbers of households. In Fig. 2.7,the electricity bill is the highest when the ECS is not deployed because the loads are notshifted during the expensive peak hours. The HDA algorithm reduces the bill as someof the HVAC appliances finish their tasks during high solar radiation hours and consumeless during peak hours. Our proposed algorithm further reduces the bill as more loads areshifted from peak hours.Fig. 2.8 shows the PAR of the aggregate load versus the percentage of deferrable loadβ. The proposed algorithm reduces the PAR from 1.81 to 1.49 when β = 40%. As βincreases from 0% to 40%, the PAR decreases because more deferrable loads can be shiftedfrom peak hours to off-peak hours.42Chapter 2. DSM under High Penetration of Rooftop PV Units0 5 10 15 20 25 30 35 40Percentage of Deferrable Load β (%)1.41.451.51.551.61.651.71.751.81.851.9PARofAggregateLoadWithout ECS DeploymentHDA AlgorithmProposed AlgorithmFigure 2.8: Comparison of the aggregate load PAR of the proposed algorithm and theHDA algorithm.2.4.4 Impact of Forecast ErrorForecast error refers to the difference between the forecasted value and the actual value ofthe PV energy generation or the load. Forecast errors persist as both solar radiation andload are intrinsically stochastic. As forecasted values are used in the proposed algorithm,the forecast errors may degrade its performance. We use the mean absolute percentage43Chapter 2. DSM under High Penetration of Rooftop PV Unitserror (MAPE) to measure the forecast error, which is defined as12T∑t∈T(|Eˆg(t)− Eg(t)|Eg(t)+|∑j∈Sm Eˆmj (t)−∑j∈Sm Emj (t)|∑j∈Sm Emj (t)),where Eˆg(t) and∑j∈Sm Eˆmj (t) are the forecasted values of the energy generation of thePV unit and the energy consumption of must-run loads in time slot t ∈ T , respectively.We adopt different scenarios to model the uncertainty in energy generation. Therefore,the forecasted energy generation is calculated as Eˆg(t) =∑k∈K P(ωk)Eˆg(t,ωk), whereEˆg(t,ωk) is the forecasted energy generation under scenario ωk in time slot t. In Figs. 2.9and 2.10, we allocate the same MAPE for the energy generation and must-run load. This isbecause the gap between the energy generation and must-run load affects the load shifting.We are particularly interested in the performance of the proposed algorithm under differentforecast errors for this gap. Furthermore, we consider the summation of the voltage risemagnitude at different buses for different time slots, i.e.,∑t∈T∑n∈B(Vn(t) − V¯ )+, as ametric to evaluate the severity of the voltage rise problem. We have V¯ = 1.05 pu as weconsider the ANSI C84.1 standard [31] which is used in North America. Fig. 2.9 showsthe voltage rise versus the forecast error. As illustrated in Fig. 2.9, the performance ofthe proposed algorithm in mitigating the voltage rise problem degrades when the forecasterror is significant. However, our proposed algorithm is able to reduce the magnitude ofvoltage rise even in the presence of forecast error, since the ECS tends to shift the loadto high solar radiation hours. The voltage magnitude is affected by the load profile ofmultiple households while each household may overestimate or underestimate its energygeneration and consumption. The forecast errors of different households sometimes haveopposite impacts on the voltage magnitude.The forecast error may also degrade the ability of the proposed algorithm to reduce themonthly electricity bill of the users. We use the saving on the electricity bill as a metric of44Chapter 2. DSM under High Penetration of Rooftop PV Units0% 10% 20% 30% 40%Forecast Error (MAPE)00.511.522.5Voltage Rise (pu)Without ECS DeploymentHDA AlgorithmProposed AlgorithmFigure 2.9: The voltage rise versus the forecast error. The y-axis is the summation of thevoltage rise over time slots and buses, which can be written as∑t∈T∑n∈B(Vn(t) − V¯ )+,where V¯ is the upper limit of the allowed voltage variation and (·)+ = max{ · , 0}.performance. The saving is defined as the difference between the monthly bill of the userfor the proposed algorithm and the bill in the case when the ECS is not deployed. To havea baseline to compare with, we consider the case in which the energy generation and loadsare perfectly known. For the case with complete information of generation and demand, wecan achieve the maximum energy saving. In the following, we calculate the saving underdifferent forecast errors. The saving as a function of the forecast error is shown in Fig.2.10. Results show that the saving decreases as the forecast error increases. However, thesaving is still considerable in the presence of forecast error. This is because the ECS tendsto shift the load from peak hours when the energy consumption price is high to off-peakhours.45Chapter 2. DSM under High Penetration of Rooftop PV Units0% 10% 20% 30% 40%Forecast Error (MAPE)60%70%80%90%100%Saving on Electricity BillProposed Algorithm, Price1Proposed Algorithm, Price2Figure 2.10: The saving on the monthly electricity bill versus the forecast error. The y-axis is the percentage of the saving. The baseline of the percentage is the saving whenthe energy generation and loads are perfectly known. The saving decreases as the forecasterror increases.2.4.5 Impact of PV Power Generation and Parameter ηThe power generation of the PV units is important to the voltage rise problem. In thissection, we study the performance of the proposed algorithm for different daily powergenerations of the rooftop PV units and different values of parameter η. Fig. 2.11 depictsthe voltage rise magnitude versus the daily power generation of rooftop PV units underthe proposed algorithm. As shown in the figure, the voltage rise magnitude increases asthe power generation of rooftop PV units increases. Fig. 2.12 shows the magnitude of46Chapter 2. DSM under High Penetration of Rooftop PV Units0 4 8 12 16 20 24Daily rooftop PV power generation (kWh)00.511.522.533.54Voltage Rise (pu)Without ECS DeploymentHDA AlgorithmProposed AlgorithmFigure 2.11: The voltage rise versus the daily power generation of the rooftop PV unit.The y-axis is the summation of the voltage rise over time slots and buses, which can bewritten as∑t∈T∑n∈B(Vn(t)− V¯ )+. We choose η = 2× 10−4.the voltage rise versus parameter η. When parameter η is small, the ECS does not shiftload to mitigate the voltage rise. As parameter η increases, the voltage rise decreases sincethe ECS tends to put more weight on the voltage rise mitigation rather than reducing theelectricity bill. The installed capacity of the PV unit is assumed to be 3 kW in Fig. 2.12.2.5 SummaryIn this chapter, we proposed a residential energy consumption scheduling algorithm forareas with high penetration of rooftop PV units. The proposed algorithm aimed to jointlyshave the peak load and mitigate the voltage rise problem during high solar radiation47Chapter 2. DSM under High Penetration of Rooftop PV Units1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5x 10−400.511.52Parameter ηVoltage Rise (pu) Proposed AlgorithmFigure 2.12: The impact of parameter η. The y-axis is the summation of the voltage rise∑t∈T∑n∈B(Vn(t)− V¯ )+. η is a parameter of the proposed algorithm. The performance ofthe cases when the ECS is not installed or the HDA algorithm is used are not affected byη.hours. We modeled the adverse effect of excessive energy export on causing abnormallyhigh voltage magnitude as an external cost based on the power flow analysis. A stochasticprogram was formulated to minimize the household electricity bill and the external costunder the uncertainty in PV power generation. The objective function was convex piecewiselinear when the feed-in tariff did not exceed the energy consumption price. Otherwise, theobjective function was non-convex piecewise linear and we transformed the problem intoan MIP problem. Simulation results confirmed that the proposed algorithm reduced theelectricity bill and the PAR of the aggregate load. The proposed algorithm also mitigated48Chapter 2. DSM under High Penetration of Rooftop PV Unitsthe voltage rise problem caused by high penetration of rooftop PV units but the effect wassubject to the forecast error, generation from PV units, and parameter selection.49Chapter 3Frequency Regulation CapacityScheduling for EVs3.1 IntroductionIn this chapter, we study the EV frequency regulation service under the performance-basedcompensation scheme. As discussed in Chapter 1, the performance-based compensationscheme has been introduced by the ISOs in recent years but its significance on the EVfrequency regulation service has not been investigated in the literature [56, 60–64, 72,105–107]. As the hourly regulation capacity is purchased by the ISOs, we consider anhourly EV frequency regulation capacity scheduling problem under the performance-basedcompensation scheme.One method for improving the revenue of the EV frequency regulation service under theperformance-based compensation scheme is to ensure that the EVs follow the AGC signalreliably. However, an EV cannot follow the AGC signal to charge if it is fully charged andit cannot discharge if its battery is depleted of energy. On the other hand, the batteriesof the EVs may either get fully charged or depleted unexpectedly when the EVs followthe random AGC signal all the time. The authors in [46, 108] show that a filtered AGCsignal can improve the reliability of the EVs in following the AGC signal but can be usedonly by those ISOs (e.g., PJM) which require that the regulation up capacity be equalto the regulation down capacity. Since the performance-based compensation scheme has50Chapter 3. Frequency Regulation Capacity Scheduling for EVsbeen implemented by the ISOs, it is desirable to design new algorithms which improve therevenue of the EV frequency regulation service for this new compensation scheme. Themain contributions of this chapter can be summarized as follows:• We model the EV frequency regulation service under the performance-based com-pensation scheme. This new compensation scheme has been implemented recentlyby the ISOs in U.S. but its significance on the EV frequency regulation service hasnot yet been studied.• We develop a new problem formulation for scheduling the regulation capacity ofthe EVs. Thereby, the revenue under the performance-based compensation schemeis introduced in the objective function. Moreover, we use a robust optimizationframework in the formulation to encourage the EVs to follow the uncertain AGCsignal reliably most of the time. An efficient algorithm is developed to solve theformulated problem.• We evaluate the performance of the proposed algorithm and the EV frequency regu-lation service by simulations. To this end, real AGC signal and prices from PJM areused in our simulations. The results show that the proposed algorithm can improvethe revenue of the EV frequency regulation service by around 10% compared to abenchmark algorithm under the performance-based compensation scheme.We note that our work is different from [64] and [109]. First, we consider a performance-based compensation scheme whereas [64] and [109] focus on capacity-based compensation.Second, we use robust optimization to account for the uncertainty in the AGC signal. Theproposed robust optimization algorithm does not require an hourly dynamic update of thecapacity. On the other hand, the authors of [64] and [109] use Markov decision processand stochastic dynamic programming to model an aggregator which dynamically updates51Chapter 3. Frequency Regulation Capacity Scheduling for EVsits capacity at the beginning of each hour. The ISOs in U.S. have different market rulesand may or may not allow an hourly dynamic update of the capacity [49, 68].This chapter is organized as follows. The system model is introduced in Section 3.2.In Section 3.3, we formulate the frequency regulation capacity scheduling problem anddevelop an efficient algorithm. Numerical results are given in Section 3.4. A summary isgiven in Section 3.5.3.2 System ModelAn EV frequency regulation scheme is illustrated in Fig. 3.1. An aggregator coordinatesEVs to provide frequency regulation service to the ISO. First, the aggregator aggregates thehourly frequency regulation capacities of the EVs. The ISO purchases the capacities andthe aggregator enters into a contract with the ISO to provide frequency regulation service.Next, during the operation period, the aggregator retrieves the AGC signal issued by theISO every few seconds (e.g., every 2-6 seconds, depending on the ISO’s requirements)and broadcasts the AGC signal to the EVs. The EVs are obliged to provide frequencyregulation service by changing their real-time charging or discharging power based on theAGC signal. The information exchange between ISO, aggregator, and EVs is enabled bya two-way communication infrastructure.We denote the operation hours by H = {1, . . . , H} and the set of EVs by M ={1, . . . ,M}. In each hour h ∈ H, EV i ∈M has a baseline charging power xi(h), regulationup capacity vui (h), and regulation down capacity vdi (h). Our goal is to optimize the valuesof xi(h), vui (h), and vdi (h) to improve the frequency regulation revenue.52Chapter 3. Frequency Regulation Capacity Scheduling for EVsISOBatteryEVAGC signalTwo-way communication linkRegulation capacities1 2AggregatorNFigure 3.1: The EVs provide frequency regulation service by following an AGC signalissued by the ISO. The AGC signal is generated by the ISO in real-time and is random.3.2.1 Randomness of the AGC SignalThe AGC signal is generated by the ISO according to the real-time mismatch betweengeneration and load in the power grid. The AGC signal is updated in short intervals. Wedivide one hour into multiple time slots. Each time slot corresponds to the duration ofone interval of the AGC signal, i.e., one time slot lasts a few seconds. Let T = {1, . . . , T}denote the set of time slots in one hour. The AGC signal in time slot t ∈ T at hour h ∈ His denoted by q(h, t) ∈ [−1, 1].The AGC signal indicates the amount by which EV i ∈M should increase or decreaseits charging power, compared to the baseline charging power xi(h). A positive AGC signal(i.e., q(h, t) > 0) indicates that the power generation is lower than the load. In this case,EV i provides the regulation up service by multiplying the AGC signal with its regulationup capacity vui (h) and decreasing the charging power accordingly. A negative AGC signal53Chapter 3. Frequency Regulation Capacity Scheduling for EVs(i.e., q(h, t) < 0) indicates that the power generation is higher than the load. EV iprovides the regulation down service by multiplying the AGC signal with its regulationdown capacity vdi (h) and increasing the charging power accordingly. Note that for chargersthat comply with the Society of Automotive Engineers (SAE) J1772 standard [110], thecharging power can be changed by adjusting the duty cycle of the pulse width modulationin the charger’s pilot signal [105]. In this chapter, we assume that the charging power canbe continuously adjusted. On the other hand, in practice, some EVs may be equipped withsimple chargers which can only be turned on and off. These EVs are still able to providefrequency regulation service by turning on and off certain number of chargers [105].We analyze the AGC signal on an hourly basis as ISOs typically purchase hourly ca-pacities. We denote the regulation up component and the regulation down component ofthe AGC signal within hour h by fu(h) and f d(h), respectively. We havefu(h) =1T∑t∈T[q(h, t)]+ , h ∈ H, (3.1)f d(h) =1T∑t∈T[−q(h, t)]+ , h ∈ H, (3.2)where [x]+ = max{x, 0}. Let ei(h) denote the charged energy for EV i within hour h.Then, ei(h) can be expressed asei(h) =1T∑t∈T(xi(h)− vui (h) [q(h, t)]+ + vdi (h) [−q(h, t)]+ )= xi(h)− vui (h)fu(h) + vdi (h)fd(h). (3.3)The terms vui (h)fu(h) and vdi (h)fd(h) represent the discharged and charged energy dueto following the regulation up and regulation down components of the AGC signal duringhour h, respectively. Note that we assume EVs follow the AGC signal in every time slot54Chapter 3. Frequency Regulation Capacity Scheduling for EVs0 0.050.1 0.150.2 0.250.3 0.350.4 0.4500.050.10.150.20.250.30.350.40.4500.10.20.30.4Regulation downcomponent fd(h)Joint probabilitymass functionRegulation upcomponent fu(h)Figure 3.2: Joint distribution of the hourly regulation up and regulation down componentsof an AGC signal. The distribution is obtained by analyzing the AGC signal data recordsfrom PJM for 2,208 hours. The figure reveals the randomness of the regulation up andregulation down components.because EVs are fast ramping resources that have zero lost opportunity cost [46, 49].The random AGC signal leads to an uncertainty in the charged or discharged energyof the EVs, as shown in (3.3). Hence, we studied the statistical joint distribution of fu(h)and f d(h), i.e., P(fu(h), f d(h)), by analyzing the AGC signal data records of PJM, for theperiod from March 1, 2014 to May 31, 2014. The results are shown in Fig. 3.2 and revealthat fu(h) and f d(h) may deviate significantly from their expected values.3.2.2 Performance-based Frequency Regulation CompensationThe ISO makes two payments under the performance-based frequency regulation compen-sation scheme. First, the ISO purchases the regulation up and regulation down capacitiesat hour h at prices pu(h) and pd(h), respectively. Second, the ISO makes another pay-55Chapter 3. Frequency Regulation Capacity Scheduling for EVsment for EVs that follow the AGC signal (e.g., the regulation market performance clearingprice (RMPCP) in PJM). The corresponding performance price is denoted by pc(h). Letmu(h) and md(h) denote the summation of the absolute changes of the regulation up andregulation down elements of the AGC signal, respectively. We havemu(h) =∑t∈T∣∣∣ [q(h, t)]+ − [q(h, t− 1)]+ ∣∣∣, h ∈ H, (3.4)md(h) =∑t∈T∣∣∣ [−q(h, t)]+ − [−q(h, t− 1)]+ ∣∣∣, h ∈ H, (3.5)where |·| denotes the absolute value. We assume the communication links between the ISO,aggregator, and EVs are reliable. If EV i follows the AGC signal at hour h, it receivesa performance related payment of pc(h)(vui (h)mu(h) + vdi (h)md(h)). Let 1i,h denote anindicator function, which is equal to 1 when EV i follows the AGC signal in hour h, and isequal to 0 otherwise. We denote the price for purchasing energy at hour h by pe(h). Therevenue for EV i at hour h is denoted by ri(vui (h), vdi (h), xi(h)), and can be written asri(vui (h), vdi (h), xi(h))= −pe(h)(xi(h)− vui (h)fu(h) + vdi (h)fd(h))+ ξi( (pu(h)vui (h) + pd(h)vdi (h))+ 1i,hpc(h)(vui (h)mu(h) + vdi (h)md(h)) ), (3.6)where −pe(h)(xi(h)− vui (h)fu(h) + vdi (h)fd(h)), pu(h)vui (h) + pd(h)vdi (h), and pc(h)(vui (h)mu(h)+ vdi (h)md(h)) represent the cost for the charged energy, payment for the regulationcapacity, and the payment for following the AGC signal, respectively. The coefficientξi ∈ [0, 1] in (3.6) is referred to as the performance score in [49]. It is determined by theISO according to the performance of the EVs in following the AGC signal in past hours.If EV i fails to follow the AGC signal at hour h (i.e., 1i,h = 0), the ISO will degrade the56Chapter 3. Frequency Regulation Capacity Scheduling for EVsscore in future hours. Note that the ISO takes into account the performance over a longperiod (e.g., 100 hours in PJM) to calculate the performance score. Hence, it is difficult toexplicitly calculate the degradation of the revenue and performance score when EV i failsto follow the AGC signal at hour h. Instead, we aim to attain a performance score closeto 1 by ensuring that EVs follow the AGC signal most of the time.Equation (3.6) models the essence of the performance-based compensation scheme in[49]. In practice, for the calculation of the performance score, the ISO first comparesthe trajectory of the EVs’ charging power with the AGC signal and uses a correlationcoefficient to measure the accuracy of the EVs’ frequency regulation service for each hour[49]. The delay of the EVs in responding to the AGC signal is taken into account by theISO. The performance score is then calculated by the ISO according to the performanceof the EVs over the past hours. In order to have a tractable problem formulation, wemodel the revenue calculation method in [49] by equation (3.6), which reflects severalimportant characteristics of the performance-based compensation scheme. First, the ISOuses a performance score as a multiplier for the revenue. The EVs should track the AGCsignal quickly, accurately, and reliably to improve the performance score. Second, anadditional payment for following the AGC signal is introduced in the performance-basedcompensation scheme. Third, the payment for following the AGC signal is calculated basedon the absolute change of the AGC signal, which is referred to as the mileage of the AGCsignal in [49]. Although the ISOs have implemented different market rules under Order 755issued by the FERC, the above three characteristics are common to all performance-basedcompensation schemes.The EVs’ limited battery capacity makes it challenging to follow the AGC signal reli-ably. An EV cannot charge when its battery is fully charged, even if this is mandated bythe AGC signal. Similarly, it cannot discharge when its battery is depleted. In the next57Chapter 3. Frequency Regulation Capacity Scheduling for EVssection, we use a robust optimization framework to obtain the hourly regulation capacitiesand enable EVs to follow the AGC signal most of the time.3.3 Problem FormulationIn this section, we formulate an EV frequency regulation capacity scheduling problembased on the robust optimization framework in [111]. The formulated problem aims toselect vui (h), vdi (h), and xi(h) in order to maximize the revenue under the performance-based compensation scheme. The uncertain parameters of the AGC signal fu(h) andf d(h), h ∈ H are considered in the formulated problem to ensure that the EVs followthe AGC signal most of the time. We note that we do not aim to guarantee that EVsalways strictly follow the random AGC signal, because such a guarantee may lead to avery conservative solution (i.e., a small regulation capacity) and reduce the revenue of theEVs.We adopt the robust optimization framework in [111] to formulate our problem. We usethis framework because some of the constraints in our problem include multiple uncertainparameters and should be satisfied with a high probability. Note that the framework in[111] is different from most of the other robust optimization frameworks which considerhard constraints (e.g., [112]). A robust optimization framework with hard constraintsensures that the solution is always feasible when the uncertain parameters are within theiruncertainty sets. We use the framework in [111] instead of hard constraints to ensure theEVs follow the AGC signal in most of the time.The basic idea of the framework in [111] is that, although a single uncertain parametersometimes may take its worst case value, it rarely happens that all the parameters take theirworst case values simultaneously. Hence, the framework in [111] introduces a tunable designparameter to adjust the number of uncertain parameters which will take their worst case58Chapter 3. Frequency Regulation Capacity Scheduling for EVsvalues simultaneously. With an appropriate selection of the number of parameters assumingworst case values, the probability for the scenarios in practical systems to be worse than thescenarios considered in the formulation is small. In this case, the constraints which includeuncertain parameters are satisfied with a high probability. In this chapter, we borrow thebasic idea of framework in [111] and introduce an integer parameter η ∈ {0, 1, . . . , H}. Ourformulation ensures that the EVs follow an AGC signal where the unknown parametersfu(h) and f d(h) take worst case values in at most η hours and take their expected valuesin the remaining hours.Let τ represent an arbitrary hour in the set H(h) = {1, . . . , h}. We denote S ⊆ H(h)as the set of hours when fu(τ) and f d(τ) take their worst case values. The cardinality ofset S is denoted by |S|. The expected values of fu(τ) and f d(τ) are denoted by µu and µd,respectively. The uncertainty sets of fu(h) and f d(h) are denoted by fu(h) ∈ [0, fu,max],f d(h) ∈ [0, f d,max]. The constants fu,max and f d,max denote the maximum values of theregulation up and regulation down components of the AGC signal, respectively. Thevalues of fu,max and f d,max can be obtained by analyzing historical AGC signal data (seeSection 3.4). We consider two cases for fu(τ) and f d(τ). In the first case, fu(τ) andf d(τ) take worst case values in the set of hours S to increase the SOC of the EVs (i.e.,fu(τ) = 0, f d(τ) = f d,max, τ ∈ S). Let si(0) denote the initial SOC of EV i at the beginningof the operation hours. We denote the battery capacity of EV i as Bi. We assume thecharging efficiencies of both charger and battery are close to one. For EV i ∈M and hourh ∈ H, we have the following constraintsi(0) + max{S⊆H(h)∣∣ |S|≤η}(1Bi∑τ∈S(xi(τ)− vui (τ)0 + vdi (τ)fd,max)+1Bi∑τ∈H(h)\S(xi(τ)− vui (τ)µu + vdi (τ)µd) )≤ smax, (3.7)59Chapter 3. Frequency Regulation Capacity Scheduling for EVswhere smax is the maximum SOC of the EV battery (e.g., smax = 1). {S ⊆ H(h)∣∣ |S| ≤ η}are all subsets of H(h) where the number of elements is at most η. In the selected hoursτ ∈ S, the unknown parameters fu(τ) and f d(τ) take worst case values (i.e., fu(τ) =0, f d(τ) = f d,max) to increase the SOC of EV i, see the first term following the max in(3.7). In the remaining hours (i.e., τ ∈ H(h) \ S), fu(τ) and f d(τ) take the expectedvalues, see the last term on the left hand side of (3.7). Constraint (3.7) can be equivalentlyrewritten assi(0) + max{S⊆H(h)∣∣ |S|≤η}((1Bi∑τ∈S(xi(τ) + vdi (τ)fd,max)− 1Bi∑τ∈S(xi(τ)− vui (τ)µu + vdi (τ)µd))+(1Bi∑τ∈S(xi(τ)− vui (τ)µu + vdi (τ)µd)+ 1Bi∑τ∈H(h)\S(xi(τ)− vui (τ)µu + vdi (τ)µd) ))≤ smax.(3.8)Constraint (3.8) is equivalent to (3.7) as we remove and then add the same components inthe second and third lines. Constraint (3.8) can be simplified assi(0) + max{S⊆H(h)∣∣ |S|≤η}(1Bi∑τ∈Svui (τ)µu + vdi (τ)(fd,max − µd))+1Bi∑τ∈H(h)(xi(τ)− vui (τ)µu + vdi (τ)µd)≤ smax. (3.9)We replace the max operator in (3.9) using max g(x) = −min(−g(x)), where g(x) denotesan arbitrary function, and obtain60Chapter 3. Frequency Regulation Capacity Scheduling for EVssi(0)−1Bimin{S⊆H(h)∣∣ |S|≤η}(∑τ∈S−vui (τ)µu + vdi (τ)(µd − f d,max))+1Bi∑τ∈H(h)(xi(τ)− vui (τ)µu + vdi (τ)µd)≤ smax, i ∈M, h ∈ H. (3.10)Constraint (3.10) is equivalent to (3.7) as shown in (3.8) and (3.9).We now consider the second case when the unknown para-meters take worst case valuesto reduce the SOC of the EVs (i.e., fu(τ) = fu,max, f d(τ) = 0, τ ∈ S). The followingconstraint keeps the SOC above a minimumsi(0) +1Bimin{S⊆H(h)∣∣ |S|≤η}∑τ∈S(vui (τ)(µu − fu,max)− vdi (τ)µd)+1Bi∑τ∈H(h)(xi(τ)− vui (τ)µu + vdi (τ)µd)≥ smin, i ∈M, h ∈ H, (3.11)where smin is the minimum SOC of the EV battery (e.g., smin = 0). Constraints (3.10) and(3.11) confine the SOC of EV i to be within[smin, smax]at hour h. In this chapter, we useconstraints (3.10) and (3.11) to enable the EVs to follow the AGC signal most of the time.EVs can provide frequency regulation service in two modes, depending on whetherdischarging is allowed or not. If discharging is allowed, then there is an additional cost ofbattery degradation. In [113], a framework is proposed to estimate the battery degradationcost using a severity factor map. The severity factor map is a function which maps thetemperature and the SOC to a battery degradation factor. Let δ(θ, si(h)) denote thebattery degradation factor, where θ is the temperature and si(h) is the SOC of EV i athour h. In [114], the values of the battery degradation factor δ(θ, si(h)) are providedfor lithium-ion batteries, which are typically used in EVs. We denote the price of thebatteries in the EVs by pb. The unit of pb is $ per kW. The life cycle of the batteries isdenoted by l. Furthermore, let ai, di, and ψi denote the arrival time, departure time, and61Chapter 3. Frequency Regulation Capacity Scheduling for EVscharging demand requirement of EV i, respectively. We denote the expected values of thesummation of the absolute changes of the regulation up and regulation down componentsmu(h) and md(h) in (3.4) and (3.5) by λu and λd, respectively. The values of λu and λdcan be obtained based on (3.4) and (3.5), and the historical records of the AGC signal.We can formulate an EV frequency regulation capacity scheduling problem as follows:maximizevui (h), vdi (h), xi(h),i ∈ M, h ∈ H∑h∈H∑i∈M(E [pu(h)] vui (h) + E[pd(h)]vdi (h)+ E [pc(h)] (vui (h)λu + vdi (h)λd)− E [pe(h)] (xi(h)− vui (h)µu + vdi (h)µd)− pblδ(θ, si(h))([xi(h)]− + µu[xi(h)− vui (h)]−))(3.12a)subject to xi(h) + vdi (h) ≤ Emaxi , i ∈M, h ∈ H, (3.12b)xi(h)− vui (h) ≥ Emini , i ∈M, h ∈ H, (3.12c)vui (h), vdi (h) ≥ 0, i ∈M, h ∈ H, (3.12d)vui (h)− α∆ ≤ vdi (h) ≤ vui (h) + α∆, i ∈M, h ∈ H, (3.12e)−∆1ai≤h≤di ≤ xi(h) ≤ ∆1ai≤h≤di, i ∈M, h ∈ H, (3.12f)−∆1ai≤h≤di ≤ vui (h) ≤ ∆1ai≤h≤di, i ∈M, h ∈ H, (3.12g)−∆1ai≤h≤di ≤ vdi (h) ≤ ∆1ai≤h≤di , i ∈M, h ∈ H, (3.12h)di∑τ=ai(xi(τ)− vui (τ)µu + vdi (τ)µd) ≥ ψi, i ∈M, (3.12i)constraints (3.10) and (3.11),where [xi(h)]− denotes max{−xi(h), 0} and the termpblδ(θ, si(h))([xi(h)]− + µu[xi(h) −vui (h)]−)is the battery degradation cost. [xi(h)]− and µu[xi(h)−vui (h)]− are used to modelthe amount of the discharged energy due to a negative baseline charging power xi(h) and62Chapter 3. Frequency Regulation Capacity Scheduling for EVsthe regulation up service, respectively. Parameter ∆ is an arbitrarily large constant, e.g.,1010. α ∈ {0, 1} in (3.12e) specifies whether the ISO has separate regulation up andregulation down markets (α = 1) or requires the regulation up capacity to match theregulation down capacity (α = 0). On the other hand, E [pu(h)], E[pd(h)], E [pc(h)], andE [pe(h)] in (3.12a) denote the expected prices for regulation up capacity, regulation downcapacity, following the AGC signal, and the charged energy at hour h, respectively.Objective function (3.12a) represents an upper bound on the expected aggregate rev-enue of the EVs. An upper bound is considered because using (3.6) as the objectivefunction leads to an intractable problem. Note that performance score ξi in (3.6) can be anon-convex function with respect to the capacity. The gap between the upper bound andthe expected revenue is obtained by subtracting (3.6) from (3.12a) and can be rewritten as∑h∈H∑i∈M((1−ξi)(E [pu(h)] vui (h)+E[pd(h)]vdi (h))+(1−(ξiP(1i,h)))E [pc(h)] (vui (h)λu+vdi (h)λd)), where P(1i,h) represents the probability of 1i,h = 1. The gap approaches zeroas ξi → 1 and P(1i,h) → 1. In this chapter, we aim to ensure 1i,h = 1 most of the timeusing constraints (3.10) and (3.11).Constraints (3.12b) and (3.12c) guarantee that the real-time charging power in hour h iswithin the maximum Emaxi and minimum Emini hourly charged energy of EV i. Emini is 0 ifthe EV does not allow discharging. Otherwise, Emini can be negative. With current batterytechnology, EV owners may not allow discharging of their EVs because this may degradethe battery lifetime. However, in the future, as battery technology evolves and improves,we expect negative Emini to become a viable option. Constraint (3.12d) ensures that thefrequency regulation capacities have non-negative values. Constraint (3.12e) reflects thecharacteristics of two types of frequency regulation markets. Some ISOs (e.g., PJM, NewYork ISO) require the frequency regulation service provider to provide the same capacitiesfor regulation up and regulation down services. In this case, α = 0, i.e., constraint (3.12e)63Chapter 3. Frequency Regulation Capacity Scheduling for EVscan be rewritten as vui (h) = vdi (h). On the other hand, some other ISOs (e.g., CAISO,ERCOT) have separate regulation up and regulation down markets. In this case, α = 1,i.e., constraint (3.12e) is always satisfied as ∆ is very large. Constraints (3.12f) (3.12g),(3.12h), and (3.12i) model the different charging periods and demands of EVs. Constraints(3.12f), (3.12g), and (3.12h) ensure that the EVs provide the regulation capacity only whenthe EVs are connected with the power grid. Note that when h ∈ [ai, di], constraints (3.12f),(3.12g), and (3.12h) are always satisfied. Otherwise, when h 6∈ [ai, di], constraints (3.12f),(3.12g), and (3.12h) confine xi(h), vui (h), and vdi (h) to 0. On the other hand, constraint(3.12i) ensures that the EVs charge sufficient energy before their departure.Market prices pu(h), pd(h), pc(h), and pe(h) are unknown parameters. Although includ-ing the price uncertainty in the model can make the problem formulation more completeand may further improve the revenue, it may also make the formulation even more com-plicated. In this chapter, we focus on improving the ability of EVs to follow the AGCsignal under the uncertainty in the signal. The unknown prices also lead to uncertaintyin the revenue but may not prevent the EVs from following the AGC signal. Hence, ourmodel considers the expected values of the prices to reduce the computational complexity.Suboptimal solutions may result if the real values of the prices deviate significantly fromtheir expected values.Note that a linear model is used for the charging process of the battery system of theEVs in (3.10) and (3.11). We assume the maximum hourly charged energy is a constantand the charging efficiency is 1. This simplified linear model reduces the computationalcomplexity to solve the formulated problem. On the other hand, in practice, the maximumcharged energy is roughly a constant only when the SOC is within a certain range. Whenthe battery is nearly fully charged or depleted, the allowed charging rate of the EV can bedifferent.64Chapter 3. Frequency Regulation Capacity Scheduling for EVs3.3.1 Duality-based Problem TransformationProblem (3.12) is a non-convex problem because constraints (3.10) and (3.11) includecombinatorial optimization components. We first analyze constraint (3.10), which includesthe combinatorial component min{S⊆H(h)∣∣ |S|≤η}∑τ∈S(−vui (τ)µu + vdi (τ)(µd − f d,max)). Thiscomponent can be rewritten as follows [111]minimizew(τ), τ ∈ H(h)∑τ∈H(h)w(τ)(− vui (τ)µu + vdi (τ)(µd − f d,max) )(3.13a)subject to 0 ≤ w(τ) ≤ 1, τ ∈ H(h), (3.13b)∑τ∈H(h)w(τ) ≤ η, (3.13c)where w(τ), τ ∈ H(h), are the variables. Note that the optimal value of problem (3.13) isequal to the summation of the η smallest values of (−vui (τ)µu + vdi (τ)(µd − f d,max)) overhours τ ∈ H(h), which is equivalent to the component min{S⊆H(h)∣∣ |S|≤η}∑τ∈S(−vui (τ)µu +vdi (τ)(µd − f d,max)). Problem (3.13) is an LP problem as it optimizes variables w(τ).Furthermore, problem (3.13) is both feasible (e.g., w(τ) = 0, τ ∈ H(h)) and bounded(∑τ∈H(h)min{−vui (τ)µu+vdi (τ)(µd−f d,max), 0} is a lower bound of the objective function).From the strong duality theorem, the optimal values of problem (3.13) and its dual problemare the same. The dual problem of (3.13) can be written asmaximizey(τ), τ ∈ H(h), z−zη −∑τ∈H(h)y(τ) (3.14a)subject to y(τ) + z ≥ vui (τ)µu − vdi (τ)(µd − f d,max), τ ∈ H(h), (3.14b)y(τ), z ≥ 0, τ ∈ H(h), (3.14c)where y(τ) and z are the dual variables for constraints (3.13b) and (3.13c), respectively.65Chapter 3. Frequency Regulation Capacity Scheduling for EVsProblem (3.14) can be used to transform the combinatorial optimization componentsin constraints (3.10) and (3.11). For constraint (3.10), we replace the componentmin{S⊆H(h)∣∣ |S|≤η}∑τ∈S(−vui (τ)µu + vdi (τ)(µd − f d,max)) with the objective function in prob-lem (3.14) and add all constraints in (3.14) to problem (3.12). We use a similar ap-proach to convert constraint (3.11) in three steps. First, we convert the combinato-rial optimization component in constraint (3.11) into an LP problem by substituting−vui (τ)µu+ vdi (τ)(µd− f d,max) in problem (3.13) with vui (τ)(µu− fu,max)− vdi (τ)µd. Then,we convert the LP problem obtained in the first step into its dual problem. Finally, wesubstitute the combinatorial optimization component in constraint (3.11) with the objec-tive function of the dual problem obtained in the second step and add the constraints ofthe dual problem to problem (3.12). This leads to the following equivalent problemmaximizevui (h), vdi (h), xi(h),z1,i(h), z2,i(h),y1,i(h, τ), y2,i(h, τ)i ∈ M, τ ∈ H(h),h ∈ H∑h∈H∑i∈M(E [pu(h)] vui (h) + E[pd(h)]vdi (h)+ E [pc(h)] (vui (h)λu + vdi (h)λd)− E [pe(h)](xi(h)− vui (h)µu + vdi (h)µd)− pblδ(θ, si(h))([xi(h)]− + µu[xi(h)− vui (h)]−))(3.15a)subject to si(0) +1Bi∑τ∈H(h)(xi(τ)− vui (τ)µu + vdi (τ)µd)+1Bi(z1,i(h)η +∑τ∈H(h)y1,i(h, τ))≤ smax, i ∈M, h ∈ H, (3.15b)z1,i(h) + y1,i(h, τ) ≥ vui (τ)µu − vdi (τ)(µd − f d,max),i ∈M, τ ∈ H(h), h ∈ H, (3.15c)si(0) +1Bi∑τ∈H(h)(xi(τ)− vui (τ)µu + vdi (τ)µd)−1Bi(z2,i(h)η +∑τ∈H(h)y2,i(h, τ))≥ smin, i ∈M, h ∈ H, (3.15d)66Chapter 3. Frequency Regulation Capacity Scheduling for EVsz2,i(h) + y2,i(h, τ) ≥ vui (τ)(fu,max − µu) + vdi (τ)µd,i ∈M, τ ∈ H(h), h ∈ H, (3.15e)z1,i(h), z2,i(h), y1,i(h, τ), y2,i(h, τ) ≥ 0, i ∈ M, τ ∈ H(h), h ∈ H,(3.15f)constraints (3.12b)–(3.12i),where (z1,i(h), y1,i(h, τ)) and (z2,i(h), y2,i(h, τ)) are dual variables corresponding to the com-binatorial optimization components in constraints (3.10) and (3.11), respectively. Linearconstraints (3.15b) and (3.15d) replace constraints (3.10) and (3.11) by substituting theinvolved combinatorial optimization problems with the corresponding linear dual prob-lems. Constraints (3.15c), (3.15e), and (3.15f) are obtained from the constraints in thedual problems. However, the objective function (3.15a) contains a battery degradationcost −pblδ(θ, si(h))([xi(h)]− + µu[xi(h)− vui (h)]−), which makes problem (3.15) difficult tosolve. Note that the values of δ(θ, si(h)) can be measured experimentally and be recordedin a table [114]. To the best of our knowledge, there is no closed-form expression for thecalculation of δ(θ, si(h)) from θ and si(h). In order to make problem (3.15) tractable, weintroduce a parameter s¯i as the expected SOC during the charging period of EV i andauxiliary parameters ǫ1,i(h), ǫ2,i(h), i ∈ M, h ∈ H. Then, we rewrite problem (3.15) asfollows:maximizevui (h), vdi (h), xi(h),z1,i(h), z2,i(h),y1,i(h, τ), y2,i(h, τ)ǫ1,i(h), ǫ2,i(h)i ∈ M, τ ∈ H(h), h ∈ H∑h∈H∑i∈M(E [pu(h)] vui (h) + E[pd(h)]vdi (h)+ E [pc(h)] (vui (h)λu + vdi (h)λd)− E [pe(h)](xi(h)− vui (h)µu + vdi (h)µd)− pblδ(θ, s¯i)(ǫ1,i(h) + µuǫ2,i(h)))(3.16a)subject to ǫ1,i(h), ǫ2,i(h) ≥ 0, i ∈M, h ∈ H, (3.16b)67Chapter 3. Frequency Regulation Capacity Scheduling for EVsǫ1,i(h) ≥ −xi(h), i ∈M, h ∈ H, (3.16c)ǫ2,i(h) ≥ −xi(h) + vui (h), i ∈M, h ∈ H, (3.16d)constraints (3.12b)–(3.12i), (3.15b)–(3.15f).Constraints (3.16b), (3.16c), and (3.16d) ensure that ǫ1,i(h) ≥ [xi(h)]− and ǫ2,i(h) ≥[xi(h)− vui (h)]−. As the objective function in (3.16a) is decreasing with respect to ǫ1,i(h)and ǫ2,i(h), we have ǫ1,i(h) = [xi(h)]− and ǫ2,i(h) = [xi(h) − vui (h)]− when the objectivefunction is maximized.With the solution obtained in problem (3.16), the EVs are able to follow an AGC signalfor which parameters fu(h) and f d(h), h ∈ H, take their worst case values in at most ηhours and take their expected values in other hours. With an appropriate value of η,the probability that the AGC signal in a practical system is worse than the AGC signalconsidered in problem (3.12) is small (see Section 3.3.2). In this case, the constraints inproblem (3.12) are satisfied with a high probability with the AGC signal in a practicalsystem.3.3.2 Probability Bound and Selection of Parameter ηProblem (3.16) aims to enable the EVs to follow the AGC signal with a high probability.Let Pi denote the probability that EV i follows the AGC signal. Assume that the valuesof fu(h) and f d(h) in hour h are independent of their values in other hours. According to[111, Theorem 3], a lower bound on probability Pi is given byPi ≥1di − ai + 1di∑h=ai(1−12h−ai((1− µi)(h− ai⌊νi⌋)+h−ai∑l=⌊νi⌋+1(h− ail))), (3.17)68Chapter 3. Frequency Regulation Capacity Scheduling for EVswhere νi =η+h−ai2, µi = νi − ⌊νi⌋, and ⌊·⌋ denotes the floor function. Equation (3.17)provides a lower bound such that the probability that EV i follows the AGC signal is notless than the right hand side of (3.17).The lower bound on probability Pi in (3.17) is tuned by design parameter η. To bespecific, increasing η increases the lower bound and enables the EVs to follow the AGCsignal more reliably. On the other hand, decreasing η enables the EVs to provide morefrequency regulation capacity. The optimal choice of η depends on the characteristics ofthe AGC signal (especially whether the AGC signal will deviate in the same directionrepeatedly in multiple hours), and may vary from one ISO to another. As the statisticalcharacteristics of the AGC signal may not change frequently, a simulation study can beused to select an appropriate value of η, see Fig. 3.8 in Section 3.4.On the other hand, the lower bound on probability Pi in (3.17) is helpful for selectingan appropriate value of η. We first select a desired value for Pi, e.g., 95%. Then, as (3.17)can be calculated efficiently, we find the largest value of η (denoted by ηˆ) for which (3.17)is satisfied. The appropriate value of η should be selected from [0, ηˆ] because the righthand side of (3.17) is a lower bound on probability Pi. We use the above method in oursimulations in Section 3.4.3.3.3 Trajectory of the AGC Signal in Hour hIn this section, we verify that the uncertainty sets of fu(h) and f d(h) contain the worst-casetrajectories of the AGC signal in hour h. Note that two different AGC signal trajectoriesin hour h with the same values of fu(h) and f d(h) may have different effects on the SOCof the EV battery during hour h. On the other hand, it is difficult to directly model thetrajectory as there are a large number of possible trajectories. In this section, we studythe AGC signal trajectory on an hourly basis and have the following remark.69Chapter 3. Frequency Regulation Capacity Scheduling for EVsRemark 3.1 For an arbitrary hour h, if an EV is able to follow the AGC signal in thetwo extreme cases, namely the case when (fu(h) = fu,max, f d(h) = 0), and the case when(fu(h) = 0, f d(h) = f d,max), then the EV is able to follow the AGC signal under otherpossible trajectories within hour h, as long as fu(h) and f d(h) are in the ranges of fu(h) ∈[0, fu,max] and f d(h) ∈ [0, f d,max], respectively.We first consider two trajectories Ω1 and Ω2 of the AGC signal within hour h, for whichthe battery SOC of an EV has the largest deviation. For the first trajectory Ω1, the AGCsignal has all the regulation up component fu,max since the beginning of the hour to acertain time tˆ1 and regulation down component fd,max since tˆ1 to the end of the hour, c.f.Fig. 3.3(a). For the second trajectory Ω2, the AGC signal has all the regulation downcomponent f d,max since the beginning of the hour to a certain time tˆ2 and the regulationup component since tˆ2 to the end of the hour, c.f. Fig. 3.3(b). Ω1 and Ω2 are the worstcase trajectories of the AGC signal in hour h. The battery SOC of an EV has the largestdownward and upward deviations under Ω1 and Ω2, c.f. Figs. 3.3(c) and 3.3(d).For trajectory Ω1 in hour h, the SOC of an EV battery reaches a minimum at tˆ1, c.f.Fig. 3.3(c). In our problem formulation, we consider a case when fu(h) = fu,max andf d(h) = 0. A possible trajectory Ω3 of the AGC signal with fu(h) = fu,max and f d(h) = 0is illustrated in Fig. 3.3(e). The SOC of an EV battery at the end of hour h undertrajectory Ω3 is equivalent to the minimum SOC in hour h under trajectory Ω1, c.f. Figs.3.3(c) and 3.3(g). Note that an EV has the largest downward deviation of the battery SOCunder trajectory Ω1 in hour h. If the battery SOC of an EV remains above the minimumvalue under trajectory Ω3 at the end of hour h, the SOC will remain above the minimumvalue for other trajectories with fu(h) ∈ [0, fu,max] and f d(h) ∈ [0, f d,max]. Similarly, forthe second trajectory Ω2 in hour h, we show that the maximum SOC in hour h under Ω2is the same as the SOC at the end of hour h in a case when fu(h) = 0 and f d(h) = f d,max,70Chapter 3. Frequency Regulation Capacity Scheduling for EVs10 20 30 40 50 60−101Time (minutes)AGC signal Trajectory Ω110 20 30 40 50 60−101Time (minutes)AGC signal Trajectory Ω2(a) AGC signal trajectory Ω1 (b) AGC signal trajectory Ω25 10 15 20 25 30 35 40 45 50 55 6000.51Time (minutes)SOCSOC of EV battery under trajectory Ω15 10 15 20 25 30 35 40 45 50 55 6000.51Time (minutes)SOCSOC of EV battery under trajectory Ω2(c) SOC of EV battery under trajectory Ω1 (d) SOC of EV battery under trajectory Ω210 20 30 40 50 60−101Time (minutes)AGC signalTrajectory Ω310 20 30 40 50 60−101Time (minutes)AGC signal Trajectory Ω4(e) AGC signal trajectory Ω3 (f) AGC signal trajectory Ω45 10 15 20 25 30 35 40 45 50 55 6000.51Time (minutes)SOCSOC of EV battery under trajectory Ω35 10 15 20 25 30 35 40 45 50 55 6000.51Time (minutes)SOCSOC of EV battery under trajectory Ω4(g) SOC of EV battery under trajectory Ω3 (h) SOC of EV battery under trajectory Ω4Figure 3.3: Different trajectories of the AGC signal and the corresponding SOC of the EVbattery. The initial SOC is set to be 0.5.see Figs. 3.3(d) and 3.3(h).As explained in Remark 3.1 in the above, for a single hour h, the worst case trajectoryof the AGC signal in hour h is contained in the uncertainty sets of fu(h) and f d(h). Onthe other hand, note that the EVs need to follow the AGC signal in multiple hours. Ourproblem formulation considers the scenario where parameters fu(h) and f d(h), h ∈ H, ofthe AGC signal can take their worst case values in a limited number of hours (at most ηhours).71Chapter 3. Frequency Regulation Capacity Scheduling for EVs3.3.4 Stochastic Arrival and Departure of the EVsIn this section, we extend our problem formulation to a case where EVs have stochasticarrival and departure times. In Section 3.3.1, we have assumed that the arrival and de-parture times are known. In practice, the arrival and departure times of the EVs can beuncertain parameters.Two issues need to be considered when the EVs have uncertain arrival and departuretimes. The first issue concerns the entity to perform the scheduling of regulation capacity.When the arrival time is known, the aggregator schedules the regulation capacity in acentralized manner. However, if the arrival time is unknown, the regulation capacity needsto be determined by the EVs independently. Second, if the arrival time is unknown, theEV needs to consider different scenarios. A scenario is a possible realization of the futurearrival time, departure time, and charging demand requirement. The scenarios can begenerated based on the historical records of the EVs. Let ωk, k ∈ K denote a scenario whereK = {1, . . . , K} is the set of scenarios. The scenarios can be generated based on historicalrecords of EVs arrival and departure times. The probability for scenario ωk is denoted byP(ωk). We denote the arrival time, departure time, and charging demand requirement ofEV i for scenario ωk by ai(ωk), di(ωk), and ψi(ωk) respectively. Let vui (h, ωk), vdi (h, ωk),xi(h, ωk), z1,i(h, ωk), z2,i(h, ωk), y1,i(h, τ, ωk), y2,i(h, τ, ωk), ǫ1,i(h, ωk), and ǫ2,i(h, ωk) denotethe values of vui (h), vdi (h), xi(h), z1,i(h), z2,i(h), y1,i(h, τ), y2,i(h, τ), ǫ1,i(h), and ǫ2,i(h)under scenario ωk, respectively. For EV i ∈ M, the problem to schedule the regulationcapacity can be rewritten as72Chapter 3. Frequency Regulation Capacity Scheduling for EVsmaximizevui (h, ωk), vdi (h, ωk), xi(h, , ωk),z1,i(h, , ωk), z2,i(h, ωk),y1,i(h, τ, ωk), y2,i(h, τ, ωk),ǫ1,i(h, ωk), ǫ2,i(h, ωk),τ ∈ H(h), h ∈ H∑k∈KP(ωk)∑h∈H∑i∈M(E [pu(h)] vui (h, ωk) + E[pd(h)]vdi (h, ωk)+ E [pc(h)] (vui (h, ωk)λu + vdi (h, ωk)λd)− E [pe(h)](xi(h, ωk)− vui (h, ωk)µu + vdi (h, ωk)µd)− pblδ(θ, s¯i)(ǫ1,i(h, ωk) + µuǫ2,i(h, ωk))(3.18a)subject to −∆1ai(ωk)≤h≤di(ωk) ≤ xi(h, ωk) ≤ ∆1ai(ωk)≤h≤di(ωk),h ∈ H, k ∈ K, (3.18b)−∆1ai(ωk)≤h≤di(ωk) ≤ vui (h, ωk) ≤ ∆1ai(ωk)≤h≤di(ωk),h ∈ H, k ∈ K, (3.18c)−∆1ai(ωk)≤h≤di(ωk) ≤ vdi (h, ωk) ≤ ∆1ai(ωk)≤h≤di(ωk),h ∈ H, k ∈ K, (3.18d)di(ωk)∑τ=ai(ωk)(xi(τ, ωk)− vui (τ, ωk)µu + vdi (τ, ωk)µd) ≥ ψi(ωk),k ∈ K, (3.18e)constraints (3.12b)–(3.12e), (3.15b)–(3.15f), and (3.16b)–(3.16d).Constraints (3.18b)–(3.18e) extend constraints (3.12f)–(3.12i) under scenario ωk. Problem(3.18) can be solved to obtain the scheduling results under different scenarios.3.3.5 Aggregate Capacity in the Day-ahead MarketThe aggregator participates in the DAM of an ISO to sell the aggregate capacity of theEVs. In this section, we develop a method to determine the amount of regulation capacity73Chapter 3. Frequency Regulation Capacity Scheduling for EVssubmitted in the DAM. We aim to ensure that the EVs can provide the capacity on thenext day according to the amount submitted in the DAM. Let vu(h) and vd(h) denote theregulation up capacity and regulation down capacity submitted in the DAM, respectively.The aggregator needs to determine control variables vu(h) and vd(h) in order to participatein the DAM.The uncertain aggregate capacity of the EVs can be modeled as normal distributedbased on the Lyapunov central limit theorem (CLT) [115]. Note that the Lyapunov CLTdoes not require the involved random variables to have identical distributions, and hence,it is applicable to the regulation capacity of the EVs. According to the Lyapunov CLT,the aggregate capacities∑i∈M vui (h) and∑i∈M vdi (h) follow normal distributions when thenumber of EVs is large and both vui (h) and vdi (h) satisfy the following condition.Lyapunov condition: There exists a parameter ǫ > 0, such that the followingequalities hold [115, p. 362].limM→∞1(δv,u(h))2+ǫ ∑i∈ME[∣∣vui (h)− E [vui (h)] ∣∣2+ǫ] = 0, (3.19)limM→∞1(δv,d(h))2+ǫ ∑i∈ME[∣∣vdi (h)− E [vdi (h)] ∣∣2+ǫ] = 0, (3.20)where δv,u(h) =√∑i∈M(δv,ui (h))2 and δv,d(h) =√∑i∈M(δv,di (h))2. Here, δv,ui (h) andδv,di (h) are the standard deviation of vui (h) and vdi (h), respectively. The values of δv,ui (h)and δv,di (h) can be obtained based on the results of problem (3.18).Next, we show that both vui (h) and vdi (h) satisfy the Lyapunov condition. We introducea parameter F = maxi∈M(Emaxi −Emini ). F takes a finite value because the maximum hourlycharged energy Emaxi and minimum hourly charged energy Emini of an EV i ∈M are finite.74Chapter 3. Frequency Regulation Capacity Scheduling for EVsWe have 0 ≤ vui (h), vdi (h) ≤ F based on (3.12b) and (3.12c). Then, we have−F ≤ vui (h)− E [vui (h)] , vdi (h)− E[vdi (h)]≤ F. (3.21)According to [115, p. 362], since random variables vui (h) and vdi (h) are bounded as (3.21),the Lyapunov condition is satisfied for both vui (h) and vdi (h).Let N (0, 1) denote a standard normal distribution. The aggregate regulation up ca-pacity and regulation down capacity are modeled as normal distributions, which are givenby∑i∈Mvui (h) ∼ N(∑i∈ME[vui (h)],(δv,u(h))2), (3.22)∑i∈Mvdi (h) ∼ N(∑i∈ME[vdi (h)],(δv,d(h))2), (3.23)where E[vui (h)] =∑k∈KP(ωk)vui (h,ωk) and E[vdi (h)] =∑k∈K P(ωk)vdi (h,ωk). The valuesof vui (h,ωk) and vdi (h,ωk) are obtained by solving problem (3.18).The aggregator determines the capacities submitted in the DAM based on the distri-butions of∑i∈M vui (h) and∑i∈M vdi (h). In particular, the values of the regulation upcapacity vu(h) and regulation down capacity vd(h) submitted in the DAM are given byvu(h) =∑i∈ME[vui (h)]− Φ−1(β)δv,u(h), (3.24)vd(h) =∑i∈ME[vdi (h)]− Φ−1(β)δv,d(h), (3.25)where β is the desired confidence level (e.g., β = 0.99) and Φ−1 is the inverse of thecumulative distribution function of a standard normal distribution.75Chapter 3. Frequency Regulation Capacity Scheduling for EVsAlgorithm 3.1 Frequency regulation capacity scheduling algorithm executed by EV ibefore the operation hours1: Initialize α, η, M,H,K,E[pe(h)],E[pu(h)],E[pd(h)],E[pc(h)],h ∈ H, fu,max, fd,max, µu, µd, Emaxi , Emini , smax, smin, pb, and l2: Generate scenarios ωk, k ∈ K based on historical records of the arrival and departure timesof the EV3: Solve problem (3.18) to obtain vui (h, ωk), vdi (h, ωk), and xi(h, ωk), h ∈ H4: Send the values of vui (h, ωk), vdi (h, ωk), and xi(h, ωk), h ∈ H, k ∈ K to the aggregator3.3.6 Algorithm and Implementation IssuesAn algorithm for the EVs to schedule the frequency regulation capacity is presented inAlgorithm 3.1. In the algorithm, the EV first initializes the parameters (Line 1). Then,the EV generates scenarios based on its historical arrival and departure times (Line 2).Subsequently, problem (3.18) is solved and the results are sent to the aggregator (Lines 3– 4).The implementation of the proposed algorithm requires an autonomous scheduler, adata collecting and storage system, a charging rate controller, and user interface. A blockdiagram for these components is shown in Fig. 3.4. The arrival and departure times of theEVs should be collected and stored in a database as records. These records can be usedby an autonomous scheduler to generate scenarios for the proposed algorithm. On theday before the operation hours, the EVs schedule their capacity under different scenariosaccording to the proposed algorithm and send the results to an aggregator. The aggregatorcan generate a bid based on EVs regulation capacity using (3.24) and (3.25). The bid willbe submitted to the ISO in the market and the ISO will award a contract to allow the EVsto provide frequency regulation service. Then, during the operation hours on the next day,when an EV owner has parked the EV and plugged-in its charger, the user interface onthe charger will be shown and encourage the EV owner to specify the departure time andcharging demand requirement. If the EV owner skips this step, it means the departure76Chapter 3. Frequency Regulation Capacity Scheduling for EVsEVAutonomous SchedulerData StorageUser InterfaceCharging Rate ControllerAggregator ISOScheduling resultsScheduling resultsBidAGC signalHistorical arrival and departure timeThe actual arrival time, departure time, and charging demand specified by userContractFigure 3.4: A diagram of different components for the EV frequency regulation service.time is unknown and the EV should be charged as fast as possible. Otherwise, if the EVowner specifies the departure time and charging demand requirement, these values are sentto the data storage. The arrival time is sent to the data storage without user intervention.Then, the EV generates a new scenario based on the actual value of the arrival timeand the value of departure time specified by the user. Then, the EV runs the proposedalgorithm with the new scenario and sends the results to the charging rate controller.The regulation capacity of an EV is uncertain because of the EV’s random arrival anddeparture times. Equations (3.24) and (3.25) take into account the uncertainty of theregulation capacity and can be used to determine the values of the capacities submitted inthe DAM. Finally, the aggregator retrieves the AGC signal from the ISO and controls thereal-time charging rate of the EVs. When the AGC signal is positive at time slot t in hourh (i.e., q(h, t) > 0), the task of the aggregator is to decrease the charging rate of the EVs77Chapter 3. Frequency Regulation Capacity Scheduling for EVsby |q(h, t)|vu(h). Otherwise, when q(h, t) < 0, the task is to increase the charging rate ofthe EVs by |q(h, t)|vd(h). The aggregator assigns the charging rates to the EVs accordingto the regulation capacity of each EV. In particular, each EV i ∈ M should adjust itscharging rate to be xi(h)−vui (h)|q(h,t)|vu(h)∑i∈M vui (h)or xi(h) +vdi (h)|q(h,t)|vd(h)∑i∈M vdi (h), when q(h, t) is positiveor negative, respectively. In this manner, the aggregator provides exactly the amount ofregulation capacity submitted in the DAM.Next, we discuss the advantages and limitations of the proposed approach comparedto the works in the literature. The advantage of the proposed approach is that it takesinto account the performance-based compensation scheme and the uncertainty of the AGCsignal. The proposed approach improves the revenue of EV frequency regulation serviceunder the performance-based compensation scheme compared to the case when the un-certainty of the AGC signal is neglected, see Section 3.4. On the other hand, since theproposed approach leverages the probabilistic robust optimization framework [111], whichincludes a heuristic parameter η to control the level of conservativeness, the limitation ofthe proposed approach is that it cannot guarantee that the EVs always follow the AGCsignal. Instead, the proposed approach guarantees that the EVs follow the AGC signalwith a certain probability, c.f., Section 3.3.2.3.4 Performance EvaluationIn this section, we evaluate the performance of the proposed algorithm based on historicalrecords of the AGC signal from PJM. We conducted a statistical analysis of the AGC signaldata records [116] and obtained parameters fu,max = 0.249, f d,max = 0.302, µu = 0.134,µd = 0.145, λu = 9.336, and λd = 8.585. A sample of the AGC signal is shown in Fig. 3.5.The uncertainty of the AGC signal is taken into account in the proposed algorithm. Wecompare the proposed algorithm with a benchmark algorithm which has a deterministic78Chapter 3. Frequency Regulation Capacity Scheduling for EVs0 10 20 30 40 50 60Time (Minutes)-0.8-0.6-0.4-0.200.20.40.60.8The AGC SignalFigure 3.5: A sample of the AGC signal obtained from PJM.formulation. The benchmark algorithm is based on a similar problem formulation comparedto the formulation used in the proposed algorithm. The only difference in benchmarkalgorithm is that the values of fu(h) and f d(h) are set in a deterministic manner usingtheir expected values. In other words, the uncertainty of the AGC signal is ignored in thebenchmark algorithm. We compare with this benchmark algorithm to study the effect ofthe uncertain AGC signal on the EV frequency regulation service under the performance-based compensation scheme.The prices of the performance-based compensation scheme from PJM [116] are usedin the simulation. We averaged the regulation market capacity clearing price (RMCCP)and RMPCP of PJM from Dec. 1, 2014 to Dec. 31, 2014. The average prices are shownin Fig. 3.6. The RMCCP (i.e., pu(h) + pd(h)) is the price of providing the frequencyregulation capacity. The RMPCP (i.e., pc(h)) is introduced under the performance-basedcompensation scheme to reimburse the market participants (e.g., EVs) for following the79Chapter 3. Frequency Regulation Capacity Scheduling for EVs1 5 10 15 20 24Time (hours)0102030405060708090HourlyPrice($perMW)RMPCP×(λu(h) + λd(h))RMCCPFigure 3.6: The average hourly prices of frequency regulation service from PJM in Dec.2014.AGC signal. Note that the revenue for following the AGC signal is the product of RMPCP,the mileage of the AGC signal (i.e., λu + λd), and the capacity. On the other hand, therevenue of providing the capacity is the product of the RMCCP and the capacity. Hence,in Fig. 3.6, we present the price pc(h)(λu+λd) in order to compare with RMCCP fairly. Asshown in Fig. 3.6, the price of following the AGC signal is significant for the performance-based compensation scheme and needs to be considered in the problem formulation andthe simulations. Furthermore, we generate values of ξi using the following method tosimulate the performance score in the performance-based compensation scheme. First, wetest the algorithm outputs (i.e., xi(h, ωk), vui (h, ωk), vdi (h, ωk)) for the historical AGC signalfrom PJM. If EV i fails to follow the AGC signal in hour h under scenario ωk, we have80Chapter 3. Frequency Regulation Capacity Scheduling for EVs1i,h(ωk) = 0, where 1i,h(ωk) is the value of 1i,h in (3.6) under scenario ωk. We denote Di,kas the set of hours when 1i,h(ωk) = 0 and |Di,k| as its cardinality. As the performancescore is based on the average performance of following the AGC signal in historical hours[49], we simulate the performance score by generating values as followsξi = 1−∑k∈K |Di,k|HK.We consider a fleet of 1000 EVs. The maximum charged energy in one hour is 10 kWh [117].We first consider the case when the EVs have unidirectional chargers (charging only). Thebattery capacity of an EV is typically tens of kWh (e.g., 20 kWh for a Honda Fit and 85kWh for a Tesla Model S). We assume a battery capacity of 40 kWh for simulation purpose.We consider an overnight charging case where EVs charge at night and are used for drivingon the next day. The random arrival and departure times are generated according to [118]and there are 50 scenarios for each EV. The demand for charging energy is selected from[10, 30] kWh. We use β = 0.99 in (3.24) and (3.25) for simulation purpose.In Fig. 3.7, we compare the proposed algorithm with the benchmark algorithm. Asshown in Fig. 3.7, the proposed algorithm achieves a higher revenue than the benchmarkalgorithm under the performance-based compensation scheme. In particular, when themaximum hourly charged energy is 10 kWh, the average daily revenue increases from $133to $151. This is because the proposed algorithm yields a solution which enables EVsto follow the AGC signal most of the time. On the other hand, with the benchmarkalgorithm where the uncertainty of the AGC signal is ignored, the EVs sometimes stopfollowing the AGC signal when the EVs are fully charged. This will reduce the revenueunder the performance-based compensation scheme. Additionally, as can be observed fromFig. 3.7, for both algorithms, the revenue increases as the maximum hourly charged energyincreases. This is because the EVs can provide more regulation capacity when the chargers81Chapter 3. Frequency Regulation Capacity Scheduling for EVshave higher charging power.We study the effect of the tunable design parameter η on the proposed algorithm inFig. 3.8. We first use the probability bounds in Section 3.3.2 to determine that η ≤ 5 isneeded, given that the EVs need to follow the AGC signal with a probability of 95%. Therevenue as a function of parameter η is shown in Fig. 3.8(a). As can be observed fromFig. 3.8(a), when η increases, the revenue first increases when η is small and then decreaseswhen η is large. This is because of the tradeoff between the regulation capacity and theperformance score, which is shown in Fig. 3.8(b). As can be observed from Fig. 3.8(b), asη increases, the performance score increases while the regulation capacity decreases. Notethat under the performance-based compensation scheme, the revenue depends on both theregulation capacity and the performance score for following the AGC signal. Based on theresults in Fig. 3.8, parameter η is selected to be η = 1 in this chapter.In Fig. 3.9, we show the probability that an EV follows the AGC signal. Thereby,we compare Pi obtained via simulations with its lower bound obtained from (3.17). Foran arbitrary EV i, we assume parameters ai = 1 and di = 13. Fig. 3.9 reveals that theprobability that an EV follows the AGC signal is not less than the right hand side of (3.17).Next, we consider the situation when the EVs have bidirectional chargers which allowboth charging and discharging. There is a battery degradation cost for the discharging.We assume the battery lasts for 800 life cycles, which is inferred from a battery limitedwarranty [119]. The battery degradation severity factor is obtained from [114] where thetemperature is assumed to be 25 centigrade (◦C). According to [39], the prices of EVbatteries have dropped from $1000 per kWh to $500 per kWh from year 2007 to year 2014.For the prices of $500 per kWh, we found that the simulation results when the EVs havebidirectional chargers are actually similar to the results when the EVs use unidirectional82Chapter 3. Frequency Regulation Capacity Scheduling for EVs6 8 10 12 14 16Maximum hourly charged energy (kWh)8090100110120130140150160170180Revenue ($)Proposed algorithmBenchmark algorithm (deterministic)Figure 3.7: Daily revenue versus the maximum hourly charged energy.chargers. This is because the battery degradation cost is taken into account in the objectivefunction of our problem formulation. As the battery degradation cost is relatively highcompared to the gain obtained from discharging, the solution obtained from the proposedalgorithm will not instruct the EVs to discharge their batteries, if the price of the batteryis $500 per kWh. On the other hand, the prices of the battery are expected to drop furtherin the future. In Fig. 3.10, we present the expected daily revenue as a function of theprice of the battery. Fig. 3.10 reveals that, the revenue of the frequency regulation serviceincreases if the price of the battery can be decreased in the future. On the other hand,if the price of the battery is high, the EV batteries are not discharged for the frequency83Chapter 3. Frequency Regulation Capacity Scheduling for EVs0 1 2 3 4 5Parameter η2060100140180Revenue ($)(a) The revenue with respect to parameter η0 1 2 3 4 5Parameter η0246Capacity (MW)0.70.80.91Performancescoreξ iPerformance score ξiHourly average of regulation capacity(b) The tradeoff between the capacity and the performance score.Figure 3.8: (a) The revenue with respect to parameter η. (b) The regulation capacity andthe performance score as a function of parameter η.regulation service. As a result, the curve of the revenue becomes flat in Fig. 3.10 when theprice of the battery is high.3.5 SummaryIn this chapter, we studied the EV frequency regulation service under the performance-based compensation scheme. We first developed a model for the performance-based com-84Chapter 3. Frequency Regulation Capacity Scheduling for EVs0 2 4 6 8 10 12Parameter η0.20.30.40.50.60.70.80.91ProbabilityPiPi obtained via simulationsLower bound on Pi obtained from (3.17)Figure 3.9: Probability Pi obtained via simulations and its lower bound obtained from(3.17).pensation scheme by taking into account the market rules of the ISOs in U.S.. Then,we formulated a problem to schedule the regulation capacity of the EVs and maximizethe revenue under the performance-based compensation scheme. A robust optimizationframework was used in the formulation to enable the EVs to follow the uncertain AGCsignal most of the time. We performed numerical experiments using historical records ofthe AGC signal and prices from PJM. Simulation results showed that EVs can improvetheir revenue under the performance-based compensation scheme by taking into accountthe uncertainty of the AGC signal in their capacity scheduling.85Chapter 3. Frequency Regulation Capacity Scheduling for EVs100 200 300 400 500 600Price of battery ($ per kWh)100150200250300Expecteddailyrevenue($)Figure 3.10: Revenue as a function of the price of the battery. Note that we assumed abidirectional charger in this figure.86Chapter 4Market Participation of EVs forFrequency Regulation Service4.1 IntroductionIn this chapter, we study the market participation of the EVs to provide the frequencyregulation service to an ISO. The focus of this chapter is the two-settlement market system(i.e., DAM and RTM) of the ISOs. Moreover, there are two types of DAMs. For the firsttype of DAM, the amount of regulation capacity specified in the contract awarded in theDAM needs to be exactly fulfilled on the next day [68]. For the second type of DAM, amarket participant can settle a shortage of regulation capacity and honor the contract bypaying a penalty to the ISO, if necessary [69]. Due to EVs’ limited battery capacity, anaggregator can be used to coordinate a fleet of EVs to provide the frequency regulationservice. In this chapter, we study how does the aggregator determine its bid in the DAM.The available regulation capacity from a fleet of EVs is time-varying since the EVshave random plugged-in and unplugged times. In [57], a stochastic model for the EVs’regulation capacity is developed based on the hourly probability that an EV is connectedwith the power grid. An algorithm to assign regulation tasks to EVs with dynamic arrivaland departure times is proposed in [60]. In [120], an algorithm to estimate the regulationcapacity is proposed based on queuing theory. However, how does the aggregator determinethe bid in the DAM has not been fully addressed.87Chapter 4. Market Participation of EVs for Frequency Regulation Service1 2 3 4 5 6 7 8 9 10 11 12 13 14Time (days)unpluggedplugged-inProfile(a) Profile of EV 1.1 2 3 4 5 6 7 8 9 10 11 12 13 14Time (days)unpluggedplugged-inProfile(b) Profile of EV 2.1 2 3 4 5 6 7 8 9 10 11 12 13 14Time (days)unpluggedplugged-inProfile(c) Profile of EV 3.Figure 4.1: Sample profiles of EVs’ connection with the power grid.The regulation capacity of an aggregator is made up of the scattered, uncertain, andsmall-scale capacities of many EVs which makes the aggregator different from other mar-ket participants. Fig. 4.1 shows sample profiles of the stochastic connection periods ofEVs with the power grid. The data is collected in Vancouver, Canada [121]. As shown inFig. 4.1, the EV charging sessions are stochastic and different EVs tend to have differentcharging periods. Thus, it is necessary to account for the uncertainty of the arrival anddeparture times of the EVs in the bid. Moreover, financial risk may arise from the ran-domness of the available capacity of the EVs. The risk is the uncertainty that the profitis lower than the expected value or even becomes negative. Note that the profit can benegative when the aggregator needs to pay a penalty to honor the contract. Hence, as amarket participant, it is desirable to consider the financial risk and an effective measureof the financial risk is the CVaR [122, 123]. In summary, the aggregator is a new type ofmarket participant with unique characteristics and requires novel algorithms to determine88Chapter 4. Market Participation of EVs for Frequency Regulation Serviceits bid for the frequency regulation service. The main contributions of this chapter aresummarized as follows:• We model the market participation of an aggregator which coordinates the EVs toprovide frequency regulation service to an ISO. Thereby, we take into account twotypes of DAMs based on the market rules of CAISO and NYISO.• We formulate the problem of determining the bid of the aggregator in the DAM usingstochastic programming. Our formulation incorporates the CVaR to account for thefinancial risk caused by the uncertain capacity of the EVs. Efficient algorithms forsolving the formulated problem are provided.• We evaluate the performance of the proposed algorithm using simulations based onreal EV charging data, collected in Vancouver, Canada. The historical prices forthe frequency regulation service of CAISO and NYISO are used to study the twoconsidered types of DAMs. We compare the profit and financial risk of the aggregatorwhen it participates in the markets of CAISO and NYISO. Our results show thatthe impact of the uncertainty of the available capacity on the profit decreases as thenumber of EVs increases.The rest of the chapter is organized as follows. The system model is introduced inSection 4.2. In Section 4.3, we present the problem formulation and develop efficientalgorithms for tackling the problem. Simulation results are presented in Section 4.4, anda summary is given in Section 4.5.4.2 System ModelThe EVs can provide frequency regulation service to an ISO and participate in the DAMand RTM via an aggregator. Fig. 4.2 illustrates the operation of the aggregator which89Chapter 4. Market Participation of EVs for Frequency Regulation ServiceEVsAGC signalsignal to control the charging rateISO DAM Control CenterAggregatorcontractregulationcapacityRTM bidBidding Charging ControlFigure 4.2: The aggregator coordinates the EVs to participate in the DAM and RTM, andprovides frequency regulation service to the ISO. The aggregator needs to submit a bid toobtain a contract from the ISO in the DAM.aggregates the regulation capacity of EVs and sells the capacity to an ISO. The aggregatorfirst participates in the DAM and submits a bid to the ISO. A contract will be awardedby the ISO to the aggregator if the aggregator is one of the winners. According to thecontract, the aggregator will provide a certain amount of regulation capacity in each hourof the next day and the ISO will reimburse the aggregator for the capacity. As the EVs’regulation capacity is uncertain and can be different from the amount specified in thecontract, in the RTM on the next day, the aggregator may sell extra capacity beyond theamount submitted in the DAM. Finally, the ISO issues a real-time frequency regulationsignal to the aggregator. The aggregator sends a control signal to the EVs by dividingthe regulation signal proportionally according to the capacity of each EV. Then, the EVsprovide the frequency regulation service by changing their real-time charging rate accordingto the received signal.90Chapter 4. Market Participation of EVs for Frequency Regulation Service4.2.1 Bid of the AggregatorIn the DAM, the aggregator submits a bid to indicate its available capacity on the nextday and its asking price. In particular, the aggregator needs to specify its regulation upand regulation down capacities in the bid, i.e., the amount by which the EVs’ real-timecharging power can be decreased and increased, respectively. We use T = {1, . . . , T} todenote the set of operating hours on the next day. Let control variables vup(t) and vdn(t)denote the amount of regulation up capacity and regulation down capacity for hour t ∈ Tin the bid of the aggregator, respectively. We consider the case where the aggregator sets itsasking price honestly based on its costs. The calculation of costs can follow the guidelinesin [124]. The marginal cost to provide frequency regulation service mainly consists of thefuel cost and heat rate degradation cost [124]. The heat rate is the energy consumed by apower generator to produce 1 kWh of electricity. For resources such as EVs, these costs arezero [124, p. 54]. The authors of [63] also assume that the EVs’ marginal cost to providefrequency regulation service is zero.We model the aggregator as a price-taker, i.e., the aggregator accepts the prices an-nounced by the ISO and does not try to influence the prices. We assume that the aggregatorwins the bid and the ISO awards the aggregator a contract where the contract size is thesame as the capacity in the bid. The contract specifies that the aggregator is responsibleto provide vup(t) regulation up capacity and vdn(t) regulation down capacity at hour t ∈ Ton the next day. We denote the set of the EVs by N = {1, . . . , N}. The regulation upcapacity and regulation down capacity of EV i ∈ N at hour t ∈ T are denoted by vupi (t)and vdni (t), respectively. The values of vupi (t) and vdni (t) are uncertain when the aggregatorsubmits its bid in the DAM. In other words, the contract requires that the aggregatorprovides a certain amount of capacity on the next day while the available capacity of theEVs is uncertain.91Chapter 4. Market Participation of EVs for Frequency Regulation Service4.2.2 DAM and RTMWe study two types of DAMs. We use a binary parameter θ ∈ {0, 1} to distinguish betweenthese two types of DAMs. For the first type (θ=0), the aggregator is obliged to provide theamount of capacity specified in the contract with a probability close to one. We introducethe following chance constraints:P(vup(t) ≤∑i∈Nvupi (t))≥ γup(1− θ), t ∈ T , (4.1)P(vdn(t) ≤∑i∈Nvdni (t))≥ γdn(1− θ), t ∈ T , (4.2)where P(A) denotes the probability of event A. Parameters γup, γdn ∈ (0, 1) are the requiredconfidence levels with which the EVs can provide the capacity specified in the contract. γupand γdn take values close to one. Note that the aggregator may not be able to participate inthe market anymore, if the aggregator repeatedly fails to honor its contract [68, p. 87]. Asthe participants in electricity markets are usually assessed based on their reliability [125],we use chance constraints (4.1) and (4.2) to ensure that the EVs can provide the capacityspecified in the contract with a probability close to one when θ= 0. In contrast, for thesecond type of DAM (θ=1), a shortage of the available capacity (i.e., a positive value ofvup(t)−∑i∈N vupi (t) or vdn(t)−∑i∈N vdni (t)) is allowed and the aggregator can honor thecontract by paying a penalty to the ISO according to the amount of shortage. In this case,constraints (4.1) and (4.2) are always satisfied as θ=1. Note that the difference betweenthe two types of DAMs can be explained with economic theory [126]. The contract awardedby the ISO in the DAM can be regarded as a forward contract in economics because thecommodity (i.e., the regulation capacity) is delivered in the future (i.e., the next day).A forward contract can be settled in two ways, namely physical delivery and financial92Chapter 4. Market Participation of EVs for Frequency Regulation Servicesettlement [126]. The key difference between the two types of considered DAMs is that thecontract awarded in the first type of DAM can only be settled with physical delivery whilefinancial settlement is allowed in the second type of DAM.The aggregator may sell extra capacity (beyond the amount submitted in the DAM)in the RTM, regardless of the type of DAM of the ISO [68, 69]. Note that ISOs typicallypurchase capacity in the DAM according to their estimated demand on the next day, whichmay deviate from their real demand. Hence, ISOs purchase additional regulation capacityoccasionally in the RTM.On the other hand, for the second type of DAM, the aggregator may pay a penalty fora shortage of capacity. This penalty is used by the ISO to purchase regulation capacityfrom other market participants. In fact, this penalty is calculated based on the marketclearing prices in the RTM [69, p. 53], [127, p. 9]. Hence, the aggregator may have eitheran additional revenue or a penalty in the RTM, which needs to be considered to calculateits profit.4.2.3 Profit of the AggregatorThe profit of the aggregator has three components. The first component is the revenueearned in the DAM. Let pupdam(t) and pdndam(t) denote the prices of the regulation up capac-ity and regulation down capacity at hour t ∈ T in the DAM, respectively. The secondcomponent of the aggregator’s profit is the revenue or penalty in the RTM. We denote theaverage prices of the regulation up capacity and the regulation down capacity at hour t inthe RTM by puprtm(t) and pdnrtm(t), respectively. Note that a short time interval is typicallyused in the RTM, such as five minutes. Prices puprtm(t) and pdnrtm(t) can be obtained byaveraging the market clearing prices in the RTM over the time intervals within hour t. Weassume the aggregator is able to sell its capacity in the RTM when puprtm(t), pdnrtm(t) > 0.93Chapter 4. Market Participation of EVs for Frequency Regulation ServiceOn the other hand, if the ISO has no demand for regulation capacity in the RTM in hourt, then puprtm(t), pdnrtm(t) = 0 and the revenue in the RTM is zero. The third component ofthe aggregator’s profit are the payments by the aggregator to the EVs. We use pupev (t) andpdnev (t) to denote the prices that the aggregator pays for the regulation up capacity andregulation down capacity of an EV at hour t, respectively. The values of pupev (t) and pdnev (t)depend on the agreement between the aggregator and the fleet of the EVs. We assume thateach EV receives a payment according to its available capacity. Let r denote the profit ofthe aggregator, which is given byr =∑t∈T(pupdam(t)vup(t) + pdndam(t)vdn(t)+ (1− θ)(puprtm(t)(∑i∈Nvupi (t)− vup(t))+pdnrtm(t)(∑i∈Nvdni (t)− vdn(t))+)+ θ(puprtm(t)(∑i∈Nvupi (t)− vup(t))+ pdnrtm(t)(∑i∈Nvdni (t)− vdn(t)))−(pupev (t)∑i∈Nvupi (t) + pdnev (t)∑i∈Nvdni (t))), (4.3)where (x)+ = max{x, 0}. The first line in (4.3) is the revenue earned in the DAM. Thesecond and third lines in (4.3) denote the revenue in the RTM for the first type of DAM(θ=0). On the other hand, for the second type of DAM (θ=1), the revenue or penalty inthe RTM is given by the fourth and fifth lines in (4.3). The sixth line is the payment fromthe aggregator to the EVs. The capacities vup(t) and vdn(t) are control variables for theaggregator. The prices (pupdam(t), pdndam(t)) and (puprtm(t), pdnrtm(t)) are announced by the ISOin the DAM and the RTM, respectively. Hence, pupdam(t), pdndam(t), puprtm(t), and pdnrtm(t) areuncertain input parameters for the aggregator.94Chapter 4. Market Participation of EVs for Frequency Regulation Service4.3 Problem Formulation and AlgorithmIn this section, we formulate an optimization problem for the aggregator to determineits bid in the DAM based on an optimality criterion regarding its profit. Thereby, weconsider two metrics for the profit, namely the expected profit and the financial risk,where the risk is measured by the CVaR. In practice, the market participants not onlyconsider their expected profit but also tend to choose an option which leads to a lowerfinancial risk compared to other options. In economics, risk aversion is used to model theattitude of market participants who are reluctant to financial risk and prefer a steady payoff.Therefore, risk aversion is accounted for in our problem formulation. Our formulation alsotakes into account the random arrival and departure times of the EVs, and the uncertaintyof the prices.4.3.1 Risk-Averse Bid of the AggregatorWe consider an aggregator with risk aversion. Hence, both the expected profit and theCVaR are included in the objective function. Let α ∈ (0, 1) denote an arbitrary confidencelevel. The value at risk (VaRα) is the largest value of an auxiliary variable η for whichthe probability that the profit in (4.3) is less than η is less than or equal to 1−α, i.e.,VaRα = max{η | P(r<η)≤1− α}. The CVaRα is defined based on VaRα as the expectedprofit for the cases when the profit is not more than VaRα [122]. That is,CVaRα = Er≤VaRα(r), (4.4)where E denotes expectation with respect to random variables pupdam(t), pdndam(t), puprtm(t),pdnrtm(t),∑i∈N vupi (t), and∑i∈N vdni (t), t ∈ T . Note that both VaRα and CVaRα aremetrics for the financial risk and we use CVaRα in this chapter because it accounts for95Chapter 4. Market Participation of EVs for Frequency Regulation Servicethe profit beyond the confidence level α. In the DAM, the aggregator submits a bid whichcontains the capacity of the 24 hours of the next day. The problem for determining thebid is formulated as followsmaximizevup(t), vdn(t), t ∈ T(1− β) E(r) + β CVaRα (4.5a)subject to vup(t), vdn(t) ≥ 0, t ∈ T , (4.5b)constraints (4.1) and (4.2), (4.5c)where β ∈ [0, 1] is a tunable parameter which adjusts the weight of the expected profitE(r) and CVaRα. Increasing β puts more weight on CVaRα and improves the profit inthe worst cases. Decreasing β puts more weight on the expected profit and less weight onCVaRα. When β = 0, the aggregator aims to maximize the expected profit and neglects thefinancial risk. This is a special case which can be used to model a risk-neutral aggregator.Problem (4.5) has CVaRα in its objective function. The CVaRα is typically calculatedin a scenario-based manner [122, 123]. A scenario is a possible realization of the randominput parameters from the ISO (pupdam(t), pdndam(t), puprtm(t), pdnrtm(t)) and the EVs (∑i∈N vupi (t),∑i∈N vdni (t)). The historical prices from the past days from the ISO [128] can be used togenerate the values of pupdam(t), pdndam(t), puprtm(t), and pdnrtm(t) in different scenarios. On theother hand, the aggregator needs to estimate the available capacity of the EVs on the nextday. In particular, possible values of the regulation capacity of an EV can be obtained basedon its historical records of charging sessions. Then, the aggregate capacities∑i∈N vupi (t)and∑i∈N vdni (t) for different scenarios can be generated accordingly (c.f. Section 4.3.2).Let K denote the set of scenarios and |K| its cardinality. Let P(ωk) denote the probabilityof scenario ωk. We use ωk, k ∈ K, to denote the kth scenario. r(ωk) denotes the profitunder scenario ωk. We can rewrite problem (4.5) as follows [122]96Chapter 4. Market Participation of EVs for Frequency Regulation Servicemaximizeη, vup(t),vdn(t), t ∈ T((1− β) P(ωk)∑k∈Kr(ωk)+ β(η −P(ωk)1− α∑k∈K(η − r(ωk))+)) (4.6a)subject to constraints (4.1), (4.2), and (4.5b), (4.6b)where the expression η − P(ωk)1−α(∑k∈K (η − r(ωk))+)is used as a scenario-based estimateof CVaRα [122]. We introduce auxiliary variables φ(ωk), k ∈ K, to denote the value of(η − r(ωk))+. Problem (4.6) can be rewritten asmaximizeη, φ(ωk), k ∈ K,vup(t), vdn(t),t ∈ T((1− β) P(ωk)∑k∈Kr(ωk)+ β(η −P(ωk)1− α∑k∈Kφ(ωk))) (4.7a)subject to φ(ωk) ≥ η − r(ωk), k ∈ K, (4.7b)φ(ωk) ≥ 0, k ∈ K, (4.7c)constraints (4.1), (4.2), and (4.5b). (4.7d)Constraints (4.7b) and (4.7c) ensure that φ(ωk) ≥ (η− r(ωk))+. As the objective function(4.7a) is decreasing with respect to φ(ωk), the optimal solution of problem (4.7) is obtainedonly if either φ(ωk) = η − r(ωk) or φ(ωk) = 0. We study two cases for problem (4.7) andits chance constraints (4.1) and (4.2). First, when θ = 1, the right hand side of constraints(4.1) and (4.2) become 0 and the constraints are always satisfied. However, when θ = 0,i.e., for the first type of DAM, chance constraints (4.1) and (4.2) make problem (4.7)difficult to solve. The probabilities in constraints (4.1) and (4.2) are not amenable to anefficient solution and it is even difficult to validate whether a solution is feasible for chance97Chapter 4. Market Participation of EVs for Frequency Regulation Serviceconstraints. In the literature, chance constraints are usually tackled using approximationmethods [129–133]. The Bernstein approximation is used to tackle chance constraints in[129, 130]. It requires knowledge of the moment generating functions of the unknownparameters and these functions are hard to obtain for problem (4.7). In this chapter, weuse the framework in [131–133], which is a scenario-based method, to handle the chanceconstraints. Using the framework in [131–133], we obtain the following problem:maximizeη, φ(ωk), k ∈ K,vup(t), vdn(t),t ∈ T((1− β) P(ωk)∑k∈Kr(ωk)+ β(η −P(ωk)1− α∑k∈Kφ(ωk))) (4.8a)subject to vup(t) ≤∑i∈Nvupi,k(t), t ∈ T , k ∈ K, (4.8b)vdn(t) ≤∑i∈Nvdni,k(t), t ∈ T , k ∈ K, (4.8c)constraints (4.5b), (4.7b), and (4.7c), (4.8d)where vupi,k(t) and vdni,k(t) are the regulation up capacity and regulation down capacity ofEV i at hour t under scenario k. Problem (4.8) is a linear programming problem whichapproximates chance constraints (4.1) and (4.2) using scenario-based constraints (4.8b)and (4.8c), respectively. An interesting problem is how many scenarios are needed in orderto approximate the chance constraints. According to [131–133], if P(ωk) =1|K| , we canselect the number of scenarios as|K|≥⌈11− γ(B − 1 + ln1δ+√2(B−1)ln1δ+ln21δ)⌉, (4.9)where B is the number of variables in the chance constraint and γ = max{γup, γdn}.98Chapter 4. Market Participation of EVs for Frequency Regulation Serviceδ ∈ (0, 1] is a constant. If we select the number of scenarios according to (4.9), then thesolution obtained in problem (4.8) will satisfy chance constraints (4.1) and (4.2) with aprobability which is not less than 1−δ. We can set γup = γdn = 0.95, δ = 0.01 and generate185 scenarios such that chance constraints (4.1) and (4.2) are satisfied with a probabilitywhich is not less than 99%.4.3.2 Aggregate Regulation Capacity of EVsThe regulation capacity of an aggregator is the sum of the uncertain capacities of manyEVs. Possible values of the uncertain capacity of an individual EV can be derived basedon the historical records of its charging sessions. Let Li denote the set of recent recordsfor EV i ∈ N in the past J days. A record of charging session l ∈ Li contains the time tˆai,lwhen the EV charger is plugged-in, the time tˆdi,l when the charger is unplugged, the initialstate of charge (SOC) sai,l, and the SOC sdi,l when it departs. For simplicity of notation,we use term tˆ to denote the index of hour during J days whereas t is the index of hourin one day. That is, tˆ ∈ {1, . . . , 24 × J} and t ∈ {1, . . . , 24}. Let pe(tˆ) denote the pricesof charged energy at hour tˆ. We denote the regulation up component and the regulationdown component of the regulation signal at hour tˆ by fup(tˆ) = 1|Φ|∑τ∈Φmax{q(tˆ, τ), 0}and fdn(tˆ) = 1|Φ|∑τ∈Φmax{−q(tˆ, τ), 0}, respectively, where q(tˆ, τ) ∈ [−1, 1] denotes theregulation signal at time slot τ ∈ Φ at hour tˆ under record l. Here, Φ is the set of timeslots within one hour. On the other hand, the maximum hourly charged energy of EV i,the battery capacity of EV i, the SOC of EV i, and the price of the charged energy at hourtˆ are denoted by emaxi , bi, si(tˆ), and pe(tˆ), respectively. Let vˆupi,l (tˆ) and vˆdni,l (tˆ) denote theregulation up capacity and the regulation down capacity of EV i at hour tˆ given recordl, respectively. Constant µ > 0 denotes a coefficient of SOC at departure time, comparedto the monetary profit. The regulation capacity of EV i ∈ N for record l ∈ Ji can be99Chapter 4. Market Participation of EVs for Frequency Regulation Serviceobtained as(vˆupi,l (tˆ), vˆdni,l (tˆ), tˆ ∈ [tˆai,l, tˆdi,l])∗ =argmaxvˆupi,l(tˆ), vˆdni,l(tˆ),tˆ ∈ [tˆai,l, tˆdi,l]∑tˆ∈[tˆai,l,tˆdi,l](vˆupi,l (tˆ)pupev (tˆ) + vˆdni,l (tˆ)pdnev (tˆ)− pe(tˆ)(xi,l(tˆ)− fup(tˆ)vˆupi,l (tˆ) + fdn(tˆ)vˆdni,l (tˆ)))+ µ s(tˆdi,l) (4.10a)subject to xi,l(tˆ) + vˆdni,l (tˆ) ≤ emaxi , tˆ ∈ [tˆai,l, tˆdi,l], (4.10b)xi,l(tˆ)− vˆupi,l (tˆ) ≥ 0, tˆ ∈ [tˆai,l, tˆdi,l], (4.10c)vˆupi,l (tˆ), vˆdni,l (tˆ) ≥ 0, tˆ ∈ [tˆai,l, tˆdi,l], (4.10d)si(tˆai,l) = sai,l, (4.10e)si(tˆ+ 1) = si(tˆ) +1bi(xi,l(tˆ)− vˆupi,l (tˆ)fup(tˆ) + vˆdni,l (tˆ)fdn(tˆ)), tˆ ∈ [tˆai,l, tˆdi,l),(4.10f)s(tˆdi,l) ≥ sdi,l, (4.10g)where xi,l(tˆ) is the baseline charging rate of EV i at hour tˆ under record l. fup(tˆ)vˆupi,l (tˆ)and fdn(tˆ)vˆdni,l (tˆ) are the decrease and increase of the charged energy within hour tˆ due tothe regulation up service and the regulation down service, respectively. Hence, xi,l(tˆ)−fup(tˆ)vˆupi,l (tˆ)+fdn(tˆ)vˆdni,l (tˆ) denotes the charged energy in hour tˆ of EV i. The term µs(tˆdi,l)is introduced to model the effect of charging during the interval [tˆai,l, tˆdi,l] in the future.Constraints (4.10b) and (4.10c) ensure that the regulation capacity of the EVs are withintheir maximum and minimum charging rate. Constraint (4.10f) models the hourly updateof the SOC of EVs. Constraint (4.10g) ensures that the EV has charged sufficient energy atits departure. Note that problem (4.10) is formulated on an hourly basis because the hourlyregulation capacity is traded in the DAM. Problem (4.10) can be solved as a stochasticprogram [93] to obtain possible values of the regulation capacity of EV i for each scenario.100Chapter 4. Market Participation of EVs for Frequency Regulation ServiceAlgorithm 4.1 Bidding algorithm executed by the aggregator on the day before theoperation hours1: Initialize θ, γup, γdn, α, β, µ, T , N , and K2: Call cloud servers to execute Algorithm 4.2 for each EV and obtain matrices Gi, i ∈ N3: Construct scenarios ωk, k ∈ K, based on matrices Gi, i ∈ N , and the historical pricesin the past days from the ISO4: Solve problem (4.8) as a stochastic program to obtain vup(t) and vdn(t), t ∈ T5: Submit vup(t) and vdn(t), t ∈ T , to the ISOAn EV may have multiple or zero records in one day. Let J = {1, . . . , J} denote theset of J days during which the records were obtained. We introduce term j ∈ J to denotethe index of day. Let vupi,j (t) and vdni,j (t) denote the regulation up capacity and regulationdown capacity at day j ∈ {1, . . . , J}. Now we need to map the capacity for each record(vˆupi,l (tˆ), vˆdni,l (tˆ)) to the capacity for each day (vupi,j (t), vdni,j (t)):(vupi,j (t), vdni,j (t))=(vˆupi,l (t+24(j−1)), vˆdni,l (t+24(j−1))),if t+24(j−1) ∈ [tˆai,l, tˆdi,l],for j ∈ {1, . . . , J}, l ∈ Li,(0, 0), otherwise.(4.11)We use a J×48 matrix Gi to collect the values of the regulation capacity of EV i for eachhour in the past J days, i.e.,Gi =vupi,1(1) vdni,1(1) · · · vdni,1(24)............vupi,J(1) vdni,J(1) · · · vdni,J(24) . (4.12)When the aggregator has received matrices Gi, i ∈ N , it generates scenarios for problem(4.8) by randomly selecting values of vupi (t) and vdni (t) based on matrices Gi.The algorithms to determine the bid are presented in Algorithms 4.1 and 4.2. The101Chapter 4. Market Participation of EVs for Frequency Regulation ServiceAlgorithm 4.2 EV capacity scheduling algorithm executed in cloud servers for EV i ∈ Nwhen called by the aggregator1: Initialize emaxi , and bi2: Generate values of tˆai,l, tˆdi,l, sai,l, sdi,l, pe(tˆ), fup(tˆ), and fdn(tˆ), l ∈ Li, based on thehistorical records of the charging sessions of EV i3: for all l ∈ Li do4: Solve problem (4.10) as a stochastic program to obtain vupi,l (tˆ) and vdni,l (tˆ), tˆ ∈ [tˆai,l, tˆdi,l]5: end for6: Calculate Gi according to (4.12) and report Gi to the aggregatoraggregator executes Algorithm 4.1 and calls cloud servers to execute Algorithm 4.2 toobtain possible values of the regulation capacities of the EVs (i.e., the matrices Gi, i ∈ N ).The client-server model can be used for calculation of the possible EV regulation capacitiessuch that these values are available even when the EVs are not connected with the powergrid. The EVs report their historical records of charging sessions to the server in the cloudcomputing infrastructure (e.g., Amazon Web Services [134]) via communication systemsand these records are stored in the database of the servers. When called by the aggregator,the cloud servers execute Algorithm 4.2 for each EV to obtain matricesGi, i ∈ N , and sendthe results to the aggregator. Then, the aggregator continues the execution of Algorithm4.1, determines the values of vup(t) and vdn(t), and submits the values to the ISO. On theother hand, the aggregator may also participate in the RTM on the next day. In the RTM,the aggregator reports vup(t)−∑i∈N vupi (t) and vdn(t)−∑i∈N vdni (t) to the ISO. If thesevalues are positive, which means the aggregator has extra capacity, the ISO may purchasethis capacity and compensate the aggregator. If these values are negative and θ=1, theISO may charge a penalty to the aggregator, according to the market clearing prices in theRTM.102Chapter 4. Market Participation of EVs for Frequency Regulation Service1 24 48 72 96 120 144 168Time (hours)0204060Price($perMW) Prices in the DAMPrices in the RTM(a) The sample hourly price in a week (θ = 0, CAISO).1 5 10 15 20 24Time (hours)01020Price($perMW) Prices in the DAMPrices in the RTM(b) The average hourly price in a day (θ = 0, CAISO).1 24 48 72 96 120 144 168Time (hours)050100Price($perMW) Prices in the DAMPrices in the RTM(c) The sample hourly price in a week (θ = 1, NYISO).1 5 10 15 20 24Time (hours)0102030Price($perMW) Prices in the DAMPrices in the RTM(d) The average hourly price in a day (θ = 1, NYISO).Figure 4.3: The prices in the DAM and RTM for the frequency regulation service in CAISOand NYISO. The y-axis is the sum of the prices of the regulation up and regulation downcapacities.4.4 Performance EvaluationWe evaluate the performance of the proposed algorithms with real data collected in Van-couver, Canada [121]. The available data belongs to an “EV cloud” project in the province103Chapter 4. Market Participation of EVs for Frequency Regulation Serviceof British Columbia in Canada. Public EV charging stations are installed in the provincewhich are connected with a database via the cellular network. A small subset of detaileddata (which were collected in our campus) was available for our study. There were 2026records available for our study. Each record contains several tags including the time whenthe EV charger is plugged-in, the time when the EV charger is unplugged, and the amountof charged energy. These records are used in our simulation study to mimic the chargingsessions of the EVs. We consider a fleet of 1000 EVs. Each EV has a battery capacity of24 kWh [135]. The maximum hourly charged energy is assumed to be 6 kWh, which is atypical charging rate of a level 2 charger. The prices of the frequency regulation service inthe DAM and RTM in CAISO [136] and NYISO [128] are used in our simulations to studythe two considered types of DAMs. Fig. 4.3 shows samples of the prices in Jan. 2015from CAISO and NYISO. We used the prices from CAISO when θ=0 and the prices fromNYISO when θ=1, respectively, since CAISO uses the first type of DAM while the secondtype of DAM is used by NYISO. In our simulations, we assume pupev (t)=0.6E[pupdam(t)] andpdnev (t) = 0.6E[pdndam(t)] in (4.3). We set α=0.9, β=0.2, and µ=2 in our algorithm unlessstated otherwise.We simulated the expected profit of the aggregator when it participates in the marketsof CAISO and NYISO. The results are shown in Figs. 4.4. The aggregator is able to obtain ahigher profit in NYISO compared to the case when it participates in the market of CAISO.In particular, when the maximum hourly charged energy is 6 kWh, the daily profit is $45.2and $21.1 for NYISO and CAISO, respectively. The profit increases as the maximumhourly charged energy increases because a higher charging rate enables EVs to provide alarger frequency regulation capacity. The profit gradually saturates when the maximumhourly charged energy becomes very large because the frequency regulation capacity is alsolimited by the charging demand, when EVs are using unidirectional chargers.104Chapter 4. Market Participation of EVs for Frequency Regulation Service2 4 6 8 10Maximum Hourly Charged Energy (kWh)0102030405060ExpectedDailyProfit($)NYISOCAISOFigure 4.4: Expected profit versus maximum hourly charged energy.Next, we study the effect of the number of EVs on the market participation of theaggregator. Fig. 4.5 shows the expected profit of the aggregator divided by the numberof EVs as a function of the number of EVs. The CVaR of the aggregator divided bynumber of EVs is shown in Fig. 4.6. As shown in Figs. 4.5 and 4.6, as the number of EVsincreases, both the profit of the aggregator divided by number of EVs and the CVaR ofthe aggregator divided by number of EVs increases (i.e., the financial risk decreases). Theresults in Figs. 4.5 and 4.6 show that using an aggregator as an agent to represent a largenumber of EVs for market participation can improve the profit and reduce the financialrisk. This is because the uncertainty of the available regulation capacity is reduced as thenumber of EVs increases. We also observe that the effect of the uncertain capacity of theEVs on the profit and the CVaR is larger for an aggregator participating in the marketsof CAISO compared to the markets of NYISO. This is because CAISO requires that theaggregator must be able to fulfill the capacity submitted in DAM, which makes the bidding105Chapter 4. Market Participation of EVs for Frequency Regulation Service0 100 200 300 400 500 600 700Number of EVs00.010.020.030.040.05ProfitofAggregator/NumberofEVs($)NYISOCAISOFigure 4.5: Expected profit of aggregator divided by number of EVs versus the number ofEVs.0 100 200 300 400 500 600 700Number of EVs-0.01-0.00500.0050.010.015CVaRofAggregator/NumberofEVs($)CAISONYISOFigure 4.6: CVaR of aggregator divided by number of EVs versus the number of EVs.106Chapter 4. Market Participation of EVs for Frequency Regulation Service3.8 3.9 4 4.1 4.2 4.3 4.4 4.5Expected Profit ($)0.70.80.911.1CVaR($)β = 1 β = 0.8β = 0.6β = 0.4β = 0.2β = 0(a) Number of EVs is 100.11.5 12 12.5 13 13.5 14Expected Profit ($)22.533.5CVaR($)β = 1 β = 0.8 β = 0.6β = 0.4β = 0.2β = 0(b) Number of EVs is 300.38 39 40 41 42 43 44 45 46 47Expected Profit ($)789101112CVaR($)β = 1β = 0β = 0.8 β = 0.6β = 0.4β = 0.2(c) Number of EVs is 1000.Figure 4.7: Efficient frontier between the expected profit and the CVaR. This figure showsthe tradeoff between the expected profit and the CVaR tuned by parameter β.more sensitive to the uncertainty of the available capacity.We also analyze the tradeoff between the expected profit and the CVaR. The tradeoffcan be depicted with an efficient frontier [122], which comprises a set of efficient points. Anefficient point is a pair of expected profit and CVaR such that it is impossible to improveone of them without reducing the other one [122]. We are able to obtain the efficientfrontier when θ=1 (because the optimal solution of problem (4.5) can be obtained when107Chapter 4. Market Participation of EVs for Frequency Regulation Serviceθ=1). The efficient frontier is shown in Fig. 4.7. As can be observed from Fig. 4.7, thetradeoff between the expected profit and the CVaR can be tuned using parameter β. As βdecreases, the expected profit increases while the CVaR decreases. Moreover, the efficientfrontier depends on the number of the EVs. The simulation results of the efficient frontierwhen the number of EVs is 100, 300, and 1000, are shown in Figs. 4.7(a), 4.7(b), and4.7(c), respectively. We can select an appropriate value of β based on Figs. 4.7(a), 4.7(b),and 4.7(c) to achieve a desired tradeoff between expected profit and CVaR.4.5 SummaryIn this chapter, we studied an aggregator participating in the DAM and RTM to obtain acontract from an ISO for providing frequency regulation service. Two types of DAMs wereconsidered by taking into account the market rules of CAISO and NYISO. We formulateda stochastic optimization problem to determine the bid in the DAM. Risk managementwas taken into account in the formulation using the CVaR. An efficient algorithm wasproposed to solve the formulated problem. Numerical experiments were performed basedon real EV charging data collected in Vancouver, Canada. Simulation results showed thatan aggregator coordinating a large number of EVs (more than a few hundreds) was helpfulfor the market participation of EVs. As the number of EVs increased, the uncertainty ofthe available regulation capacity had less effect on the profit and financial risk.108Chapter 5Conclusions and Future WorkIn this final chapter, we summarize the results and highlight the contributions of this thesisin Section 5.1. We also propose ideas for future research topics in Section 5.2.5.1 Summary of ResultsThis thesis studies the energy consumption scheduling of the residential loads and EVs. Inthe following, we briefly review results and conclusions.In Chapter 2, we proposed an algorithm to schedule the energy consumption of resi-dential loads for areas with high penetration of rooftop PV units. The peak shaving andthe voltage rise issue were jointly considered in the proposed algorithm by shifting loadsfrom peak hours to hours with high solar radiation. We used the external cost and powerflow analysis to model the effect of the voltage rise. We formulated a stochastic optimiza-tion problem whereas the objective function was to minimize the electricity bill and theexternal cost. We evaluated the performance of the proposed algorithm with simulations.The results showed that the proposed algorithm can reduce the PAR of the load profileand also mitigate the voltage rise issue. Our conclusion was that the DSM scheme canbe a promising tool to handle the voltage rise issue but its effect was highly dependenton the electricity generation of PV units and other factors such as the forecast error andparameter selection.In Chapter 3, we developed a problem formulation to optimize the hourly frequency109Chapter 5. Conclusions and Future Workregulation capacity of the EVs. The performance-based compensation scheme was takeninto account. We developed a model for the performance-based compensation schemebased on the business practice manual of PJM. The uncertainty of the AGC signal was alsotaken into account in our formulation using a robust optimization framework. Numericalexperiments were performed based on data records of the AGC signal and prices from PJM.The results showed that the revenue of EV frequency regulation service can be improvedby taking into account the uncertainty of the AGC signal under the performance-basedcompensation scheme, compared to a benchmark case where the uncertainty of the AGCsignal is neglected. Our conclusion was that the performance-based compensation schemecan significantly affect the revenue of frequency regulation service and it should be takeninto account in the capacity scheduling and the economical evaluation. Moreover, underthe performance-based compensation scheme introduced by the ISOs, it became importantto take into account the uncertainty of the AGC signal when scheduling the capacity ofEVs in order to improve the revenue.In Chapter 4, we studied the market participation of an aggregator which coordinateda fleet of EVs to provide the frequency regulation service. We analyzed two types of DAMsbased on the market rules of CAISO and NYISO. The difference between these two typesof DAMs was whether the ISO allowed the contract awarded in the DAM to be settled bypaying a penalty or not. We formulated a stochastic optimization problem to determinethe bid in the DAM. Chance constraints and CVaR were used in the formulation to accountfor the uncertainty of EVs’ frequency regulation capacity and financial risks. We proposedefficient algorithms to solve the formulated problem. We evaluated the performance ofthe proposed algorithm via the simulation of the market participation of an aggregatorfor CAISO and NYISO. Our conclusion was that although the capacity of a single EV isuncertain, an aggregator coordinating a large number of EVs was able to participate in110Chapter 5. Conclusions and Future Workthe two-settlement market system and earn a profit from the market.5.2 Future WorkThe emerging smart grid needs to be efficient, reliable, and accommodate high penetrationof renewable resources. With the development of two-way communication systems in thesmart grid, the residential loads are expected to become controllable and the DSM will playan important role in the emerging smart grid. In the current research, most works havefocused on the exploration of the potential applications of the DSM to achieve differentgoals for the efficient and reliable operation of the electrical grid. In Chapters 2–4, we haveproposed algorithms to leverage residential loads and EVs to mitigate the voltage rise issueand to provide the frequency regulation service. However, the DSM is a vast research areaand many problems are still unsolved. In the following, we present some ideas for furtherresearch.1. Joint Operation of DSM and PV Inverter: In Chapter 2, we studied theDSM in areas with high penetration of PV units. We have shown that the voltagerise problem can be mitigated with DSM but the effect is subject to the forecasterror, the energy generation of rooftop PV units, and the parameter selection. Howto combine the DSM approach with the real-time operation of PV inverters, suchthat the voltage rise is guaranteed to be mitigated is an interesting problem worthexploration.2. EV Frequency Regulation Service with On-off Chargers: In Chapter 3, weproposed an algorithm to schedule the hourly frequency regulation capacity of theEVs. We have assumed that the EV’s charging rate can be continuously adjusted.This assumption is generally valid for the chargers which comply with the SAE J1772111Chapter 5. Conclusions and Future Workstandard [110]. An interesting extension of Chapter 3 is to consider the EV frequencyregulation service with on-off chargers. Note that in practice, some of the EVs areequipped with simple chargers which only allow on and off states. These EVs arestill able to provide the frequency regulation service. That is, when the AGC signalindicates regulation up, the aggregator should turn off some of the chargers. On theother hand, when AGC signal indicates regulation down, the aggregator should turnon additional EV chargers. However, the management of the on-off states of the EVchargers has not yet been fully addressed. It is an interesting problem to study howto manage the on-off states of the EV chargers to provide the frequency regulationservice while satisfying the charging demand.3. Impact of Distribution Network Constraints on EV Frequency Regula-tion Service: In Chapter 4, the market participation of an aggregator in thetwo-settlement market system is proposed. Since the aggregator needs to coordinatea fleet of EVs to adjust charging rate simultaneously, an interesting extension ofChapter 4 is to consider the constraints of the distribution network for the EV fre-quency regulation service. In particular, the frequency regulation down service canbe provided by increasing the charging power of the EVs simultaneously, which maycause significant voltage drop on the low-voltage distribution network. Thereby, thedistribution network allows a limited power flow to charge the EVs. An interestingproblem worth exploration is how to model the distribution network constraints forthe EV frequency regulation service.112Bibliography[1] A. Ipakchi and F. Albuyeh, “Grid of the future,” IEEE Power Energy Mag., vol. 7,no. 2, pp. 52–62, Mar. 2009.[2] Federal Energy Regulatory Commission, “Assessment of demand re-sponse and advanced metering,” Tech. Rep., Feb. 2011, Available at:http://www.ferc.gov/legal/staff-reports/2010-dr-report.pdf.[3] Ontario Hydro One Time-of-Use Pricing, Available at:http://www.hydroone.com/tou/Pages/Default.aspx, accessed: Sept. 2015.[4] Southern California Edison Residential Time-of-Use Rate Plan, Available at:https://www.sce.com/wps/portal/home/residential/rates/residential-plan/tou, ac-cessed: Sept. 2015.[5] W. Kempton, V. Udo, K. Huber, K. Komara, S. Letendre, S. Baker, D. Brun-ner, and N. Pearre, “A test of vehicle-to-grid (V2G) for energy storage and fre-quency regulation in the PJM system,” Tech. Rep., Nov. 2008, Available at:www.udel.edu/V2G/resources/test-v2g-in-pjm-jan09.pdf.[6] J. Tomic and W. Kempton, “Using fleets of electric-drive vehicles for grid support,”Journal of Power Sources, vol. 168, no. 2, pp. 459–468, Jun. 2007.[7] G. Strbac, “Demand side management: Benefits and challenges,” Energy Policy,vol. 36, no. 12, pp. 4419–4426, Dec. 2008.[8] N. Ruiz, I. Cobelo, and J. Oyarzaba, “A direct load control model for virtual powerplant management,” IEEE Trans. on Power Systems, vol. 24, no. 2, pp. 959–966,May 2009.[9] A. Conejo, J. Morales, and L. Baringo, “Real-time demand response model,” IEEETrans. on Smart Grid, vol. 1, no. 3, pp. 236–242, Dec. 2010.[10] Smart Pricing Program in the Illinois State of U.S., Available at:https://www.powersmartpricing.org/, accessed: May 2015.[11] A. H. Mohsenian-Rad, V.W.S. Wong, J. Jatskevich, R. Schober, and A. Leon-Garcia,“Autonomous demand side management based on game-theoretic energy consump-tion scheduling for the future smart grid,” IEEE Trans. on Smart Grid, vol. 1, no. 3,pp. 320–331, Dec. 2010.113Bibliography[12] N. Li, L. Chen, and S. H. Low, “Optimal demand response based on utility maxi-mization in power networks,” in Proc. of IEEE PES General Meeting, Detroit, MI,Jul. 2011.[13] P. Srikantha, C. Rosenberg, and S. Keshav, “An analysis of peak demand reductionsdue to elasticity of domestic appliances,” in Proc. of ACM e-Energy, Madrid, Spain,Jul. 2012.[14] C. Chen, K. G. Nagananda, G. Xiong, S. Kishore, and L. V. Snyder, “Acommunication-based appliance scheduling scheme for consumer-premise energymanagement systems,” IEEE Trans. on Smart Grid, vol. 4, no. 1, pp. 56–65, Mar.2013.[15] Z. Yu, L. McLaughlin, L. Jia, M. C. Murphy-Hoye, A. Pratt, and L. Tong, “Modelingand stochastic control for home energy management,” in Proc. of IEEE Power andEnergy Society (PES) General Meeting, San Diego, CA, Jul. 2012.[16] P. Samadi, H. Mohsenian-Rad, R. Schober, and V.W.S. Wong, “Tackling the load un-certainty challenges for energy consumption scheduling in smart grid,” IEEE Trans.on Smart Grid, vol. 4, no. 2, pp. 1007–1016, Jun. 2013.[17] L. Jia and L. Tong, “Optimal pricing for residential demand response: A stochasticoptimization approach,” in Proc. of Allerton Conference on Communication, Control,and Computing, Allerton, IL, Oct. 2012.[18] P. Samadi, H. Mohsenian-Rad, V.W.S. Wong, and R. Schober, “Real-time pricing fordemand response based on stochastic approximation,” IEEE Trans. on Smart Grid,vol. 5, no. 2, pp. 789–798, Mar. 2014.[19] S.-J. Kim and G. B. Giannakis, “Real-time electricity pricing for demand responseusing online convex optimization,” in Proc. of IEEE Innovative Smart Grid Tech-nologies Conference (ISGT), Washington, DC, Feb. 2014.[20] International Energy Agency, “Snapshot of global PV mar-kets,” Tech. Rep., Oct. 2014, Available at: http://www.iea-pvps.org/fileadmin/dam/public/report/technical/PVPS report -A Snapshot of Global PV - 1992-2014.pdf.[21] Y. Poissant and P. Luukkonen, “National survey report of PV power applica-tions in Canada,” International Energy Agency, Tech. Rep., Oct. 2013, Avail-able at: http://www.cansia.ca/sites/default/files/2014-14 rp-anu 411-pvuscl iea-pvps poissant luukkonen.pdf.[22] Solar Energy Industry Association, Available at: http://www.seia.org/research-resources/solar-industry-data, accessed: Sept. 2015.114Bibliography[23] Net Metering in the State of Illinois in U.S., Available at:www.illinoisattorneygeneral.gov/environment/netmetering.html, accessed: Feb.2016.[24] Feed-in tariff in Province of Ontario in Canada, Available at:http://fit.powerauthority.on.ca/what-feed-tariff-program, accessed: Jan. 2016.[25] A. Pedrasa, T. D. Spooner, and I. F. MacGill, “Coordinated scheduling of residentialdistributed energy resources to optimize smart home energy services,” IEEE Trans.on Smart Grid, vol. 1, no. 2, pp. 134–143, Nov. 2010.[26] I. Atzeni, L. G. Ordonez, G. Scutari, D. P. Palomar, and J. R. Fonollosa, “Demand-side management via distributed energy generation and storage optimization,” IEEETrans. on Smart Grid, vol. 4, no. 2, pp. 866–876, Jun. 2013.[27] C. O. Adika and L. Wang, “Autonomous appliance scheduling for household energymanagement,” IEEE Trans. on Smart Grid, vol. 5, no. 2, pp. 673–682, Mar. 2014.[28] R. A. Walling, R. Saint, R. C. Dugan, J. Burke, and L. A. Kojovic, “Summary ofdistributed resources impact on power delivery systems,” IEEE Trans. on PowerDelivery, vol. 23, no. 3, pp. 1636–1644, Jul. 2008.[29] R. Tonkoski, D. Turcotte, and T. El-Fouly, “Impact of high PV penetration onvoltage profiles in residential neighborhoods,” IEEE Trans. on Sustainable Energy,vol. 3, no. 3, pp. 518–527, May 2012.[30] J. D. Glover, M. S. Sarma, and T. Overbye, Power System Analysis and Design, 5thEdition. Cengage Learning, 2011.[31] American National Standard for Electric Power Systems and Equipment VoltageRatings, American National Standards Institute (ANSI) C84.1, 2006.[32] J. von Appen, M. Braun, T. Stetz, K. Diwold, and D. Geibel, “Time in the sun:The challenge of high PV penetration in the German electric grid,” IEEE Power andEnergy Magazine, vol. 11, no. 2, pp. 55–64, Mar. 2013.[33] T. Stetz, F. Marten, and M. Braun, “Improved low voltage grid-integration of pho-tovoltaic systems in Germany,” IEEE Trans. on Sustainable Energy, vol. 4, no. 2,pp. 534–542, Apr. 2013.[34] E. DallAnese, S. V. Dhople, and G. B. Giannakis, “Optimal dispatch of photovoltaicinverters in residential distribution systems,” IEEE Trans. on Sustainable Energy,vol. 5, no. 2, pp. 487–497, Apr. 2014.115Bibliography[35] “German Renewable Energy Sources Act (Erneuerbare-Energien-Gesetz),” Available at: http://www.erneuerbare-energien.de/EE/Navigation/DE/Gesetze/Das EEG/das eeg.html, accessed: May2015.[36] S. J. Steffel, P. R. Caroselli, A. M. Dinkel, J. Q. Liu, R. N. Sackey, and N. R.Vadhar, “Integrating solar generation on the electric distribution grid,” IEEE Trans.on Smart Grid, vol. 3, no. 2, pp. 878–886, Jun. 2012.[37] E. Yao, P. Samadi, V.W.S. Wong, and R. Schober, “Tackling the photovoltaic inte-gration challenge in the distribution network with deferrable load,” in Proc. of IEEESmartGridComm, Vancouver, Canada, Oct. 2013.[38] Hybrid Cars Market Dashboard, Available at: http://www.hybridcars.com/market-dashboard/, accessed: Sept. 2015.[39] B. Nykvist and M. Nilsson, “Rapid falling costs of battery packs for electric vehicles,”Nature Climate Change, vol. 5, no. 4, pp. 329–332, Apr. 2015.[40] B. Kirby, “Frequency regulation basics and trends,” Oakridge National Laboratory,Tech. Rep., Dec. 2004, Available at: http://www.osti.gov/bridge.[41] ——, “Ancillary services: Technical and commercial insights,” Tech. Rep.,Jul. 2007, Available at: http://www.consultkirby.com/files/Ancillary Services -Technical And Commercial Insights EXT .pdf.[42] Texas Energy Storage Alliance, “Storage participation in ERCOT,” Tech. Rep., Jan.2010.[43] H. Hao, B. M. Sanandaji, K. Poolla, and T. L. Vincent, “Frequency regulation fromflexible loads: Potential, economics, and implementation,” in Proc. of AmericanControl Conference (ACC), Portland, Oregon, Jun. 2014.[44] H. Ni, “PJM advanced technology pilots for system frequency control,” in Proc. ofIEEE Innovative Smart Grid Technologies (ISGT), Washington, DC, Jan. 2012.[45] Y. V. Makarov, J. Ma, S. Lu, and T. Nguyen, “Assessing the value ofregulation resources based on their time response characteristics,” Pa-cific Northwest National Laboratory, Tech. Rep., Jun. 2008, Availableat: http://www.pnl.gov/main/publications/external/technical reports/PNNL-17632.pdf.[46] C. Quinn, D. Zimmerle, and T. H. Bradley, “An evaluation of state-of-charge limita-tions and actuation signal energy content on plug-in hybrid electric vehicle, vehicle-to-grid reliability, and economics,” IEEE Trans. on Smart Grid, vol. 3, no. 1, pp.483–491, Mar. 2012.116Bibliography[47] D. Fooladivanda, C. Rosenberg, and S. Garg, “Analysis of energy- and SOC-neutralcontracts for frequency regulation with energy storage,” in Proc. of IEEE Smart-GridComm, Miami, FL, Nov. 2015.[48] PJM RegD Signal, http://www.pjm.com/committees-and-groups/task-forces/rpstf.aspx, accessed: Sept. 2015.[49] PJMManual 11: Energy and Ancillary Services Market Operations, Revision 72, Jan.2015, Available at: http://www.pjm.com/∼/media/documents/manuals/m11.ashx.[50] O. Ma, N. Alkadi, P. Cappers, P. Denholm, J. Dudley, S. Goli, M. Hummon, S. Kilic-cote, J. MacDonald, N. Matson, D. Olsen, C. Rose, M. D. Sohn, M. Starke, B. Kirby,and M. OMalley, “Demand response for ancillary services,” IEEE Trans. on SmartGrid, vol. 4, no. 4, pp. 1988–1995, Dec. 2013.[51] FERC Order 719: Wholesale Competition in Regions with Organized Elec-tric Markets, Oct. 2008, Available at: http://www.ferc.gov/whats-new/comm-meet/2008/101608/E-1.pdf.[52] M. Maasoumy, C. Rosenberg, A. Sangiovanni-Vincentelli, and D. S. Callaway, “Modelpredictive control approach to online computation of demand-side flexibility of com-mercial buildings HVAC systems for supply following,” in Proc. of American ControlConference (ACC), Portland, OR, Jun. 2014.[53] R. J. Bessa andM. A. Matos, “The role of an aggregator agent for EV in the electricitymarket,” in Proc. of 7th Mediterranean Conf. and Exhibition on Power Generation,Transmission, Distribution and Energy Conversion (MedPower), Agia Napa, Cyprus,Nov. 2010.[54] H. Hao, A. Somani, J. Lian, and T. E. Carroll, “Generalized aggregation and coordi-nation of residential loads in a smart community,” in Proc. of IEEE SmartGridComm,Miami, FL, Nov. 2015.[55] C. Goebel and D. S. Callaway, “Using ICT-controlled plug-in electric vehicles tosupply grid regulation in California at different renewable integration levels,” IEEETrans. on Smart Grid, vol. 4, no. 2, pp. 729–740, Jun. 2013.[56] S. Han, S. Han, and K. Sezaki, “Development of an optimal vehicle-to-grid aggregatorfor frequency regulation,” IEEE Trans. on Smart Grid, vol. 1, no. 1, pp. 65–72, Jun.2010.[57] ——, “Estimation of achievable power capacity from plug-in electric vehicles for V2Gfrequency regulation: Case studies for market participation,” IEEE Trans. on SmartGrid, vol. 2, no. 4, pp. 632–641, Dec. 2011.117Bibliography[58] W. Shi and V.W.S. Wong, “Real-time vehicle-to-grid control algorithm under priceuncertainty,” in Proc. of IEEE SmartGridComm, Brussels, Belgium, Oct. 2011.[59] E. Yao, V.W.S. Wong, and R. Schober, “A robust design of electric vehicle frequencyregulation service,” in Proc. of IEEE SmartGridComm, Venice, Italy, Nov. 2014.[60] S. Sun, M. Dong, and B. Liang, “Real-time welfare-maximizing regulation allocationin dynamic aggregator-EVs system,” IEEE Trans. on Smart Grid, vol. 5, no. 3, pp.1397–1409, May 2014.[61] E. Sortomme and M. El-Sharkawi, “Optimal charging strategies for unidirectionalvehicle-to-grid,” IEEE Trans. on Smart Grid, vol. 2, no. 1, pp. 131–138, Mar. 2011.[62] S. I. Vagropoulos, D. K. Kyriazidis, and A. G. Bakirtzis, “Real-time charging manage-ment framework for electric vehicle aggregators in a market environment,” acceptedfor publication in IEEE Trans. on Smart Grid, 2016.[63] S. I. Vagropoulos and A. G. Bakirtzis, “Optimal bidding strategy for electric vehicleaggregators in electricity markets,” IEEE Trans. on Power Systems, vol. 28, no. 4,pp. 4031–4041, Nov. 2013.[64] J. Donadee and M. D. Ilic, “Stochastic optimization of grid to vehicle frequencyregulation capacity bids,” IEEE Trans. on Smart Grid, vol. 5, no. 2, pp. 1061–1069,Mar. 2014.[65] SAE Charging Configurations and Ratings Terminology, Available at:http://www.sae.org/smartgrid/chargingspeeds.pdf, accessed: Sept. 2015.[66] FERC Order 755: Frequency Regulation Compensation in the Organized WholesalePower Markets, Oct. 2011, Available at: http://www.ferc.gov/whats-new/comm-meet/2011/102011/E-28.pdf.[67] A. D. Papalexopoulos and P. E. Andrianesis, “Performance-based pricing of frequencyregulation in electricity markets,” IEEE Trans. on Power Systems, vol. 29, no. 1, pp.441–449, Jan. 2014.[68] California ISO (CAISO) Business Practice Manual for Mar-ket Instruments, Version 38, Oct. 2015, Available at:http://bpmcm.caiso.com/Pages/BPMDetails.aspx?BPM=Market%20Instruments.[69] New York ISO (NYISO) Ancillary Services Manual, Version 4.2, Jun. 2015, Availableat: http://www.nyiso.com/public/webdocs/markets operations/documents/Manuals and Guides/Manuals/Operations/ancserv.pdf.[70] Electricity Reliability Council of Texas (ERCOT) Day-Ahead Mar-ket Operating Procedure Manual, Version 4.7, Apr. 2015, Available at:http://www.ercot.com/mktrules/bpm/.118Bibliography[71] R. J. Bessa, M. A. Matos, F. J. Soares, and J. A. P. Lopes, “Optimized bidding of aEV aggregation agent in the electricity market,” IEEE Trans. on Smart Grid, vol. 3,no. 1, pp. 443–452, Mar. 2012.[72] C. Jin, J. Tang, and P. Ghosh, “Optimizing electric vehicle charging with energystorage in the electricity market,” IEEE Trans. on Smart Grid, vol. 4, no. 1, pp.311–320, Mar. 2013.[73] E. Yao, V.W.S. Wong, and R. Schober, “Risk-averse forward contract for electricvehicle frequency regulation service,” in Proc. of IEEE SmartGridComm, Miami,FL, Nov. 2015.[74] E. Yao, P. Samadi, V.W.S. Wong, and R. Schober, “Residential demand side manage-ment under high penetration of rooftop photovoltaic units,” accepted for publicationin IEEE Trans. on Smart Grid, 2015.[75] E. Yao, V.W.S. Wong, and R. Schober, “Robust frequency regulation capacityscheduling algorithm for electric vehicles,” accepted for publication in IEEE Trans.on Smart Grid, 2016.[76] X. Liu, A. Aichhorn, L. Liu, and H. Li, “Coordinated control of distributed energystorage system with tap changer transformers for voltage rise mitigation under highphotovoltaic penetration,” IEEE Trans. on Smart Grid, vol. 3, no. 2, pp. 897–906,Jun. 2012.[77] M. J. E. Alam, K. M. Muttaqi, and D. Sutanto, “Distributed energy storage formitigation of voltage-rise impact caused by rooftop solar PV,” in Proc. of IEEE PESGeneral Meeting, San Diego, CA, Jul. 2012.[78] J. von Appen, T. Stetz, M. Braun, and A. Schmiegel, “Local voltage control strategiesfor PV storage systems in distribution grids,” IEEE Trans. on Smart Grid, vol. 5,no. 2, pp. 1002–1009, Mar. 2014.[79] IEEE 1547 Standard for Interconnecting Distributed Resources with ElectricPower Systems, http://grouper.ieee.org/groups/scc21/1547/1547 index.html, ac-cessed: Jan. 2016.[80] R. S. Rao, K. Ravindra, K. Satish, and S. V. L. Narasimham, “Power loss min-imization in distribution system using network reconfiguration in the presence ofdistributed generation,” IEEE Trans. on Power Systems, vol. 28, no. 1, pp. 317–325,Feb. 2013.[81] M. Brenna, E. D. Berardinis, F. Foiadelli, G. Sapienza, and D. Zaninelli, “Voltagecontrol in smart grids: An approach based on sensitivity theory,” Journal of Elec-tromagnetic Analysis and Applications, vol. 2, no. 1, pp. 669–678, Dec. 2006.119Bibliography[82] K. M. Rogers, R. Klump, H. Khurana, A. A. Aquino-Lugo, and T. J. Overbye, “Anauthenticated control framework for distributed voltage support on the smart grid,”IEEE Trans. on Smart Grid, vol. 1, no. 1, pp. 40–47, Jun. 2010.[83] Power system simulator PowerWorld, Available at: www.powerworld.com, accessed:Oct. 2015.[84] N. G. Mankiw, Principle of Economics, 6th edition. Cengage Learning, 2011.[85] L. Yu, D. Czarkowski, and F. de Leo´n, “Optimal distributed voltage regulation forsecondary networks with DGs,” IEEE Trans. on Smart Grid, vol. 3, no. 2, pp. 959–967, June 2012.[86] S. Bolognani and S. Zampieri, “A distributed control strategy for reactive powercompensation in smart microgrids,” IEEE Trans. on Automatic Control, vol. 58,no. 11, pp. 2818–2833, Nov. 2013.[87] S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge University Press,2004.[88] R. H. Coase, “The problem of social cost,” Journal of Law and Economics, vol. 3,no. 1, pp. 1–44, Oct. 1960.[89] A. D. Owen, “Renewable energy: Externality costs as market barriers,” Energy Pol-icy, vol. 34, no. 5, pp. 632–642, Mar. 2006.[90] P. Rafaj and S. Kypreos, “Internalisation of external cost in the power generationsector: Analysis with global multi-regional MARKAL model,” Energy Policy, vol. 35,no. 2, pp. 828–843, Feb. 2007.[91] X. Li and C. W. Yu, “Impacts of emission trading on carbon, electricity and renewablemarkets: A review,” in Proc. of IEEE PES General Meeting, Minneapolis, MN, Jul.2010.[92] Q. Chen, C. Kang, Q. Xia, and D. Kirschen, “Optimal flexible operation of a CO2capture power plant in a combined energy and carbon emission market,” IEEE Trans.on Power Systems, vol. 27, no. 3, pp. 1602–1609, Aug. 2012.[93] A. Shapiro, D. Dentcheva, and A. Ruszczynski, Lectures on Stochastic Programming:Modeling and Theory. SIAM, 2009.[94] A. Yona, T. Senjyu, A. Y. Saber, T. Funabashi, H. Sekine, and C. H. Kim, “Ap-plication of neural network to 24-hour-ahead generating power forecasting for PVsystem,” in Proc. of IEEE PES General Meeting, Pittsburgh, PA, Jul. 2008.120Bibliography[95] E. Lorenz, D. H. J. Hurka, and H. G. Beyer, “Irradiance forecasting for the powerprediction of grid-connected photovoltaic systems,” IEEE Journal of Selected Topicsin Applied Earth Observations and Remote Sensing, vol. 2, no. 1, pp. 2–10, Mar.2009.[96] P. Bachera, H. Madsena, and H. A. Nielsen, “Online short-term solar power forecast-ing,” Solar Energy, vol. 83, no. 10, pp. 1772–1783, Oct. 2009.[97] Solar purchase program in Georgia state in U.S., Available at:http://www.georgiapower.com/pricing/residential/residential-tariffs.cshtml, ac-cessed: Mar. 2015.[98] Optimization software MOSEK, Available at: http://www.mosek.com/, accessed:June 2013.[99] Optimization software CPLEX, Available at: http://www-01.ibm.com/software/commerce/optimization/cplex-optimizer/, accessed: Jan.2016.[100] J. A. Jardini, C. M. V. Tahan, M. R. Gouvea, S. U. Ahn, and F. M. Figueiredo, “Dailyload profiles for residential, commercial and industrial low voltage consumers,” IEEETrans. on Power Delivery, vol. 15, no. 1, pp. 375–380, Jan. 2000.[101] Southern California Edison Net Metering Program, Available at:https://www.sce.com/NR/sc3/tm2/pdf/ce158-12.pdf, accessed: Oct. 2015.[102] Solar Data from National Renewable Energy Laboratory (NREL), Available at:http://www.nrel.gov/rredc/solar data.html, accessed: Aug. 2013.[103] Atmel Smart Metering Solutions, Available at:http://www.atmel.com/applications/Smart Energy/default.aspx, accessed: Nov.2014.[104] Microsoft Azure: Cloud Computing Platform, Available at:http://azure.microsoft.com, accessed: Nov. 2014.[105] E. Sortomme and K. W. Cheung, “Intelligent dispatch of electric vehicles performingvehicle-to-grid regulation,” in Proc. of IEEE Int’l Electric Vehicle Conf. (IEVC),Greenville, SC, Mar. 2012.[106] J. Lin, K.-C. Leung, and V. O. K. Li, “Online scheduling for vehicle-to-grid regulationservice,” in Proc. of IEEE SmartGridComm, Vancouver, Canada, Oct. 2013.[107] C. Wu, H. Mohsenian-Rad, and J. Huang, “Vehicle-to-aggregator interaction game,”IEEE Trans. on Smart Grid, vol. 3, no. 1, pp. 434–442, Mar. 2012.121Bibliography[108] C. Quinn, “Evaluation of distributed energy storage for ancillary service provision,”Master’s thesis, Colorado State University, 2011.[109] J. Donadee and M. D. Ilic, “Stochastic co-optimization of charging and frequencyregulation by electric vehicles,” in Proc. of IEEE North American Power Symposium(NAPS), Champaign, IL, Jan. 2012.[110] SAE J1772 Standard on Electric Vehicle and Plug-in Hybrid Electric Vehicle Con-ductive Charge Coupler, Available at: http://standards.sae.org/j1772 201001/, ac-cessed: Jan. 2016.[111] D. Bertsimas and M. Sim, “The price of robustness,” Operations Research, vol. 52,no. 1, pp. 35–53, Feb. 2004.[112] A. Ben-Tal and A. Nemirovski, “Robust optimization methodology and applica-tions,” Mathematical Programming, vol. 92, no. 3, pp. 453–480, May 2002.[113] J. Donadee and M. D. Ilic, “Estimating the rate of battery degradation under a sta-tionary Markov operating policy,” in Proc. of IEEE PES General Meeting, NationalHarbor, MD, Jul. 2014.[114] S. Onori, P. Spagnol, V. Marano, Y. Guezennec, and G. Rizzoni, “A new life estima-tion method for Lithium-ion batteries in plug-in hybrid electric vehicles applications,”Intl. J. of Power Electronics, vol. 4, no. 3, pp. 302–319, June 2012.[115] R. B. Ash and C. A. Doleans-Dade, Probability and Measure Theory, Second Ed.Elsevier Science, 1999.[116] PJM Interconnection Market Based Frequency Regulation,http://www.pjm.com/markets-and-operations/ancillary-services/mkt-based-regulation.aspx, accessed: Sept. 2015.[117] Teslamotors outlet charging adapter guide, Available at:http://www.teslamotors.com/en CA/charging/outlet, accessed: Sept. 2015.[118] A. Santos, N. McGuckin, H. Nakamoto, D. Gray, and S. Liss, “Summary oftravel trends: 2009 national household travel survey,” U.S. Department of Trans-portation, Federal Highway Administration, Tech. Rep., Jun. 2011, Available at:http://nhts.ornl.gov/2009/pub/stt.pdf.[119] Tesla Motors model S vehicle warranty, Available at:https://www.teslamotors.com/sites/default/files/blog attachments/ms vehiclewarranty.pdf, accessed: Jan. 2016.[120] A. Y. S. Lam, K.-C. Leung, and V. O. K. Li, “Capacity estimation for vehicle-to-grid frequency regulation services with smart charging mechanism,” IEEE Trans. onSmart Grid, vol. 7, no. 1, pp. 156–166, Jan. 2016.122Bibliography[121] FleetCarma EV Cloud: Collecting electric vehicle charging data in the provinceof British Columbia, Canada, Available at: https://www.fleetcarma.com/evcloud,accessed: Sept. 2015.[122] A. J. Conejo, M. Carrion, and J. M. Morales, Decision Making Under Uncertaintyin Electricity Markets. Springer, 2010.[123] R. T. Rockafellar and S. Uryasev, “Optimization of conditional value-at-risk,” Jour-nal of Risk, vol. 2, no. 3, pp. 21–41, Apr. 2000.[124] PJM Manual 15: Cost Development Guidelines, Revision 26, Nov. 2014, Availableat: http://www.pjm.com/∼/media/documents/manuals/m15.ashx.[125] IEEE Standard Definitions for Use in Reporting Electric Gener-ating Unit Reliability, Availability, and Productivity, Available at:http://www.nerc.com/docs/pc/gadstf/ieee762tf/762-2006.pdf, accessed: Jan.2016.[126] S. J. Deng and S. S. Oren, “Electricity derivatives and risk management,” The In-ternational Journal of Energy, vol. 31, no. 6, pp. 940–953, May 2006.[127] New York ISO Rate Schedule 3, Nov. 2013, Available at:http://www.nyiso.com/public/webdocs/markets operations/committees/mc/meeting materials/2015-12-17/Agenda%2006 BTMNG%20MST%20Sec%2015 3%20Revisions.pdf.[128] NYISO Market Operation Data Records, Available at:http://www.nyiso.com/public/markets operations/market data/pricing data/index.jsp, accessed: Jan. 2016.[129] A. Nemirovski and A. Shapiro, “Convex approximation of chance constrained pro-grams,” SIAM Journal of Optimization, vol. 17, no. 4, pp. 969–996, Nov. 2006.[130] W. W. L. Li, A. Y. J. Zhang, A. M. C. So, and M. Z. Win, “Slow adaptive OFDMAsystems through chance constrained programming,” IEEE Trans. on Signal Process-ing, vol. 58, no. 7, pp. 3858–3869, Jul. 2010.[131] A. M. C. So and A. Y. J. Zhang, “Distributionally robust slow adaptive OFDMAwith soft QoS via linear programming,” IEEE Journal on Selected Areas in Commu-nications, vol. 31, no. 5, pp. 947–958, May 2013.[132] A. Nemirovski and A. Shapiro, “Scenario approximations of chance constraints,” inProbabilistic and Randomized Methods for Design under Uncertainty, G. Calafioreand F. Dabbene, Eds. London: Springer, 2005, ch. 1, pp. 3–47.123Bibliography[133] M. C. Campi and S. Garatti, “The exact feasibility of randomized solutions of uncer-tain convex programs,” SIAM Journal of Optimization, vol. 19, no. 3, pp. 1211–1230,Dec. 2008.[134] Amazon Web Services, Available at: https://aws.amazon.com/, accessed: Sept. 2015.[135] Nissan Leaf EV (Model S) Specifications, Available at:http://www.nissan.ca/en/electric-cars/leaf/versions-specs/version.s.html, accessed:Sept. 2015.[136] CAISO Open Access Same-time Information System (OASIS),http://oasis.caiso.com/, accessed: Sept. 2015.124