Development of robust participation strategies for an aggregator of electric vehicles in multiple types of electricity markets

The ﬂexible charging and discharging behaviours of electric vehicles (EVs) can provide a promising proﬁt opportunity for EV owners or aggregators in deregulated market. How-ever, a major barrier to overcome is the uncertain market prices. In this context, this study aims to establish a novel and robust model for EV aggregators to optimally participate in multiple types of electricity markets while considering the uncertainty of the prices in various markets. First, the electricity market framework and the obligation of an EV aggregator are introduced, and each EVaggregator is entitled to participate in day-ahead, real-time, and regulation service markets. To cope with forecasting errors of market prices and other uncertainties, a robust optimisation model is developed to schedule the charging/discharging of EVs in different markets. A stepwise bidding strategy is then proposed to formulate robust market bids, offering more ﬂexibility in handling volatile market prices in actual electricity markets. After the clearing of multiple types of electricity markets, a real-time dispatch model for EV aggregators is proposed for allocating the settlements and regulation orders to each EV. Finally, scenario-based simulation is employed to demonstrate the effectiveness of the proposed bidding strategy using actual market data and regulation signals.


INTRODUCTION
In the past decade, battery technology has developed significantly, and the number of electric vehicles (EVs) has increased rapidly around the globe, making modern power systems transform to sustainable and smart energy systems [1,2]. EVs have the potential of contributing to future power systems largely owing to their flexibilities and fast responses to system operation requirements, which enhance system operation performance in both security and economics. Meanwhile, the arbitrary and undirected charging behaviours of significant numbers of EVs could result in negative impacts on the systems [3][4][5]. Thus, proper methods must be adopted to control/direct the EV charging behaviours, such as through direct EV charging control [6] or market-based control [7]. Thus far, a popular approach for handling the scattered and uncertain EV charging demand is to introduce an aggregation agent, normally denoted as an EV aggregator (EVA) [8,9]. An EVA controls the charging schedules of individual EVs and interacts with system operators through control directives and market signals. Generally, the market-based mechanism is more effective as it provides economic incentives and protects the confidential information for EVs. The state-of-the-art study on multimarket bidding models can be found in [10], whereas bidding strategies for conventional power producers are reported in [11][12][13][14]. The bidding strategies for EVs and EVAs adhere to similar principles as those of conventional market participants. Existing studies have explored the participation and optimal bidding strategies of EVAs in electricity markets. For example, the bidding strategy in the day-ahead (DA) market, considering the market price uncertainties, is formulated as a bilevel optimisation in [15]; the bidding strategy in DA and secondary reserve market is optimised in [16,17]; the energy prices risks are assessed through stochastic optimisation in [18]; the impact of the uncertain reserve market price on the EVA bidding strategy is modelled using a scenario-based stochastic optimisation method in [19]; a stochastic bidding strategy of EV charging stations in a frequency containment reserve market is presented in [20]; the combined bidding in multiple ancillary services employing vehicle-to-grid (V2G) are proposed in [21][22][23][24]; regulation capacity bids are solved through stochastic optimisation in [25]; coordinated bidding and management of EVAs in combined energy and ancillary markets are developed in [26,27]; the uncertainties of energy and ancillary service market prices are modelled by the well-established conditional value-at-risk (CVaR) approach in [28][29][30]; the fuzzy optimisation technique is employed in [31] to manage uncertainty and coordinate strategic bidding problems for EVs in various ancillary service markets; the robust optimisation method is employed in [32] to develop the optimal bidding strategy with energy market price uncertainty considered; EVs are combined with the energy hub concept in [33] and the power, heating, and cooling systems are integrated into the developed strategic bidding model.
As validated by these studies, the most advantageous way of utilising an EVA is to optimise its charging schedules in energy markets based on market price predictions while reserving certain capacities for the ancillary markets. Thus, a reasonable coordinated bidding strategy in multimarket seems to be the desired solution to the market participation of the EVAs. However, the following concerns are yet to be addressed: • Although various approaches such as stochastic and robust optimisation methods have been applied in handling uncertainty factors associated with the strategic bidding problems, focus has been laid on the determination of the participated quantity in the electricity market while ignoring the bidding price or addressing it in simplified ways [34,35]. Therefore, existing bidding strategies for EVAs cannot appropriately accommodate stochastic market price fluctuations. • Aggregation of EVs is extensively employed in reducing the complexity of the bidding model of EVs to enable their participation in an electricity market. Disaggregation of the market clearing results (i.e. energy schedules and ancillary service capacities) and their subsequent allocation to individual EVs remain unsolved. In addition, dispatching the charging/discharging power of EVs by EVAs to provide the required services is also a problem to be addressed.
In this study, we aim to develop a complete solution for EVAs to participate in electricity markets as follows: (i) the interaction mechanisms between the EVAs and EV owners; (ii) the optimisation model for energy scheduling of EVAs in multimarkets; (iii) the robust approach for deciding the bid prices while considering uncertainties; (iv) the real-time dispatch model for EVAs for allocating cleared market results and ancillary service requirements.
The remainder of this paper is organised as follows: In Section 2, we discuss the electricity market frames, particularly the DA, real-time (RT) energy market, and regulation service (RS) market. The relationships between EVA and EV owners are presented in Section 3, along with the analysis of EVA costs and revenues. Robust optimisation is employed in Section 4 to solve the proposed multimarket energy scheduling model, and stepwise bidding functions are introduced to formulate the bid prices that are adaptive to market price fluctuations. In Section 5, we develop the dispatch model to allocate energy and RS schedules to individual EVs during the RT operations. Scenariobased simulation is adopted to validate the robustness of the proposed bidding strategy, and the design of test scenarios is demonstrated in Section 6. The case studies are conducted in Section 7, where the effectiveness of the proposed strategy is verified through test scenarios based on real-world data, and the results of EVAs with different risk preferences are compared. Conclusion and future work are discussed in Section 8.

Energy market
Generally, EVAs will purchase energy for EV charging through market transactions in a deregulated market environment. The transactions could be long-term energy contracts between market participants, such as bilateral contracts with steady prices and lower financial risks, or short-term trades in wholesale energy markets, such as DA and RT market with volatile prices and higher financial risks. In this study, we emphasise the interactions between the EVAs and the independent system operator (ISO) in short-term DA and RT markets, whereas the long-term transactions are not considered. According to existing market designs [10,11,36], the proposed market frame is discussed in the subsequent subsections. The flowchart of the market operation is provided in Figure 1. The proposed EVA bidding models and dispatch strategies are also listed in Figure 1 to clearly illustrate the connections among different sections/models. As shown in Figure 1, three models, namely, Model-D, Model-R, and Model-H are introduced to solve EVA optimisation problems at different market stages. Detailed descriptions and formulations of these models will be described in Sections 4 and 5.

DA market
The EVAs could simultaneously submit both the demand bids (battery charging) and the supply offers (battery discharging) to the ISO for energy market clearing. The bids/offers for the hourly DA energy market of operating day D should be submitted to the ISO before 12:00 on day D-1, implying that the EVAs should forecast their energy demand and the DA market clearing price (MCP). The bids/offers submitted to the ISO should contain both the power quantities and bid prices for all the operating hours of day D. The EVAs can submit a single power-price block (pair) bid/offer, or multiple power-price blocks (pairs) for each operating hour. In this study, it is assumed that the EVAs will submit their bids/offers in the form of stepwise functions that contain a series of power-price pairs for each trading hour. From 12:00 to 16:00 on day D-1, the ISO will assess the received DA market bids/offers for day D, and clear the DA market for maximum social welfare while considering security constraints. The DA MCP and successful bids/offers will be

RT market
The RT market for day D will open for bids/offers at 16:00-18:00 on day D-1, immediately after the clearance of the DA market. Market participants who are willing to revise their cleared DA schedules will rebid in the RT market. The RT market clears after 18:00, and the cleared hourly RT MCP and successful bids/offers will be published by the ISO as well. From 0:00 to 24:00 on the operating day D, the ISO may also organise hour-ahead balancing markets for market participants to revise their schedules with more accurate forecasts, which are not considered in this study. All the market participants are supposed to follow their successful bids/offers in the energy markets throughout the operating day D, as significant deviations from the cleared power quantities may be subjected to financial penalties by the ISO.

Ancillary service market
Among different ancillary service market frames around the globe, ancillary services generally comprise regulation, spinning/non-spinning reserve, operating reserve, and other services such as reactive power support and black-start capability [13,37,38]. It is assumed that the EVAs will bid for the RS only in the hourly ancillary services market because the RS normally enjoys the highest price, thus EVAs are likely to compensate for their costs for providing the service. The RS contains regulation up and down services, for which bidding can be done separately in some markets. In this study, an EVA will bid a unique price for a specific regulation capacity, which applies to both the regulation up and down services. In addition to the regulation capacity price, the EVA also gets paid for the actual dispatched regulation up/down energy based on the regulation up/down energy prices.
The RS market for day D opens at 16:00 on day D-1, after the clearance of DA energy market. The EVAs may bid for regulation capacities if they can increase/decrease their charging demands compared to the successful DA bids. Notably, the EVAs have to meet certain requirements before they are eligible to bid in the RS market. For example, the least regulation capacity in Pennsylvania-New Jersey-Maryland (PJM) ancillary market is 0.1 MW [36], implying that the EVAs with available capacities less than 0.1 MW are not eligible to bid in the RS market for that particular period. The RS market will be cleared after the market closes at 18:00 on day D-1. From 0:00 to 24:00 on the operating day D, the market participants with successful RS bids are supposed to increase/decrease their generation/consumptionbased on regulation signals from ISO and the volume of the accepted regulation capacity.
The difference between the accepted regulation capacity and the provided regulation power will result in financial penalties to the service providers as well. The EVAs can revise their regulation bids at least 1 h ahead of the operating hour, to handle unexpected circumstances or forecasted errors. The modification of accepted regulation capacity at hour H+1 on day D can be revised no later than the beginning of hour H on day D. However, only the modification that decreases the regulation capacity is eligible.

ROLE OF ELECTRIC VEHICLE AGGREGATORS
Thus far, several approaches are available for the EV owners to recharge their batteries when required, namely, distributed home charging, centralised charging, and battery swapping [39]. The charging requirements of EVs are generally affected by the battery characteristics, daily driving mileage statistics, and driving habits. Although the charging/discharging of an individual EV could be arbitrary and stochastic, the behaviour of an EVA normally shows statistical regularity because a relatively steady daily power consumption can be expected. Through smart metering and communication device, an EVA is capable of monitoring and controlling the charging/discharging of individual EVs when they are parked at home or work.

EVA responsibilities
The major function of an EVA is to fulfill the EV charging requirements of its clients. However, the actual EV charging demand of any single EV could vary daily, making it difficult for the EVA to determine whether the customers' charging need has been met or not. According to the vehicle travel statistics [40], a substantially large proportion of EVs are parked throughout the day, and the most frequent vehicle usages occur during the morning hours. Thus, the following criterion will be adopted to evaluate the charging services of an EVA: Criterion: The EVA is supposed to charge each EV battery to a predefined level before a predefined time-point of each day and financial penalties may apply if the EVA fails to fulfill the task. Meanwhile, if the EV owner's driving behaviour makes it impossible for the EVA to charge the battery to the agreed level at the specific time, then the EVA is expected to charge the battery as fully as it can to avoid penalties.
Generally, the amount of energy stored in batteries is indicated by the state-of-charge (SoC). Let Ω E and Ω T denote the set of EVs and time slots; the proposed criterion for a single EV (∀e ∈ Ω E ) during the considered time horizon (∀t ∈ Ω T ) can be described as follows: where S C e,t , p ch e,t , and p dis e,t denote the SoC level, charging power, and discharging power of e-th EV at time t, respectively; B E e , S C,min e , S C,max e , p ch,max e , and p dis,max e denote the battery capacity, minimum/maximum SoC levels, and the maximum charging/discharging power of the e-th EV, respectively; t 0 and T con e represent the initial time and the predefined EV charging check time-point, whereas S Con e denotes the predefined SoC level of the e-th EV; ch and dis denote the charging/discharging efficiency of the battery. u ch e,t is introduced because the battery cannot be charging and discharging simultaneously.
Based on the proposed criterion, the responsibility of the EVA can be quantified while mitigating the interferences introduced by the stochastic driving patterns of individual EVs. The check time-point of EVA charging services can normally be set during morning hours when most EV customers expect to start their daily usage. The EVA will charge the idle EV batteries at its convenience for the remainder of the day.
Apart from the EVAs that provide home charging services, EV charging stations and parking slots with EV charging facilities can also provide EV charging services, which can be analysed similarly. Thus, this study only emphasizes on the home charging EVAs.

Cost and profit of an EVA
As economical entities, the EVAs are expected to evaluate their costs and maximise their profits while performing their EV charging responsibilities. The profit of an EVA can be expressed as follows: where R profit a , R income a , and C cost a denote the net profit, total income, and total cost of the a-th EVA, respectively. Ω A represents the set of EVAs.
In this study, the EV customers are assumed to pay for the charging services only, whereas the EVA needs to pay for all the mandatory fees including utility tariffs. Thus, the costs of the a-th EVA is normally caused by: (i) depreciation of metering device and charging facilities; (ii) grid tariffs and network loss fees, denoted as C GT a ; (iii) battery degeneration cost in the case of V2G employments, denoted as C BD a ; (iv) cost of purchasing electricity in both the DA and RT energy markets, denoted as C PE a ; (v) all categories of penalties, such as failure to fulfill the EV charging requirement and failure to comply with successful bids in the electricity markets, denoted as C Penal a ; and (iv) Other miscellaneous expenses.
Ignoring the charging facility depreciation and miscellaneous expenses owing to difficulties in quantifying them, the daily cost of the a-th EVA (∀a ∈ Ω A ) can be separately calculated as follows: where loss and GT denote the expense ratios of network losses and grid tariffs, respectively; c dis denotes the battery degeneration cost of unit discharging energy and dis,com denotes the compensation rate such that the EV customers will gain profits for discharging their batteries; E net a and E dis a denote the net charged and discharged energy, respectively; DA t and RT t along with P DA,CH a,t and P RT,CH a,t denote the MCP and charging power of DA and RT market at time t, respectively; uc , dev , ru , and rd along with E uc a , E dev a , E ru,m a , and E rd,m a denote the penalty rate and energy mismatch for undercharging EV batteries, the deviations from the cleared bids, and the mismatches of regulation up/down services, respectively.
The income of the a-th EVA is mostly generated from charging services, which is denoted as R CH a . Moreover, the EVA also gets paid for its V2G operation and provision of RS, which are denoted as R DIS a and R REG a , respectively. The total income of the a-th EVA (∀a ∈ Ω A ) can be calculated as follows: The EVA may optimise its profit by maximising the objective function (7) while accommodating the constraints in (1)- (6), for all the EVs within its control. Statistically, the desired daily charging energy is relatively steady, thus, leaving the profit of an EVA mostly related to its success in the electricity markets.

EVA BIDDING STRATEGY IN MULTIMARKETS
According to the proposed market framework, the EVAs will have to bid in the DA, RT, and RS markets on day D-1 with forecasted market prices and uncertain EV driving patterns of day D. The EVAs are subject to significant risks if their bidding strategy will not handle these uncertainties appropriately. In this study, robust optimisation is employed to handle the uncertainties introduced by the fluctuating market prices.

Robust bidding model in DA energy market
Although the bidding in RT and RS markets starts after the clearance of the DA energy market, the EVAs still have to consider the price fluctuations in those markets while optimising their DA market bids because the clear results from the DA market lay the foundations for the bidding in RT and RS markets.
Theoretically, the charging constraints of every EV should be considered in optimising the EVA bidding strategies, which will increase the model complexity, parameter randomness, and computational burden. In this study, the EVAs ignore the individual EV charging/discharging characteristics and only consider the holistic charging requirements.
Considering the price uncertainties in DA, RT, RS markets as well as the regulation energy price uncertainties for regulation up/down services, a robust bidding model for the a-th EVA (∀a ∈ Ω A ) in the DA market is developed based on the existing studies in [41][42][43] and formulated as follows: P act,min Here, the objective (17) maximizes the profit of the a-th EVA on the operating day D with the selected robust coefficient denote the extra discharged and charged energies for providing regulation up/down services at time t, respectively, the values of which are calculated from (47)-(49). Equation (33) calculates the total stored energy in EV batteries at time t, denoted as P net a,t . Term ΔP net,REG a,t denotes the net energy changes for providing RS, which is calculated from (37). Equation (34) defines the daily charging requirements of all the EV batteries, andP total a denotes the daily energy demand based on historical data of the a-th EVA. Equation (35) constrains the total charged energy considering the EV customer usage requirement demonstrated in (1). Term P cr a,t c denotes the total cumulative charging demand based on historical data, whereas t c denotes the check time-points. Equation (36) calculates the discharged/charged energy for regulation up/down services, and RU t and RD t denote the estimated probabilities that the regulation up/down will be deployed at time t, respectively.
Constraints (38)-(39) denote the total charging and discharging power limits, where P CH,max a,t and P DIS,max a,t denote the maximum total charging/discharging power of the a-th EVA at time t, respectively. Equations (40)-(43) constrain the bid charging/discharging power in DA and RT markets, and the integer variable U CH a,t is introduced based on (4)- (6). Equation (44) constrains the quantities of P Rcap a,t , and integer variable U REG a,t is introduced to indicate the participation of EVA a in the RS market at time t. Term P Rcap,min denotes the minimum regulation capacity that the ISO accepts, whereas REG is a coefficient that represents the proportion of P CH,max a,t acceptable as the regulation capacity of the a-th EVA at time t.
The net consumption of the EVA may be negative (discharging power to the grid) if the accepted regulation capacity is larger than the accepted charging power in energy markets. In this case, the EVA has to discharge its available batteries, and extra battery life losses and costs will be introduced, which are calculated from (45)-(51). Term P act,min a,t denotes the net discharging power, and the integer variable U EXT a,t represents whether the EVA has to discharge its batteries during regulation up service at time t, respectively. Terms and are set to very low and very high positive values, respectively. The robust bidding model in the DA market is denoted as Model-D. Notably, Model-D employs the equivalent linear formulation from [44]. In summary, Model-D can be described as follows.
Objective (17) The optimised results of Model-D are identified with superscript * , for example, the optimised P DA,CH a,t is denoted as P * DA,CH a,t . Different bidding results will be obtained by modifying the value of Γ DA , representing the optimal bidding strategies of EVAs with various risk tolerance levels. Notably, Model-D only solves the quantities of EVA charging/discharging power throughout the operating day D, whereas the bidding prices remains unsettled.

Robust bidding model in RT and RS markets
After the ISO clears the DA market, the EVAs will receive both the DA MCP and their accepted bids, based on which they can alter their energy schedules by bidding in the RT market and provide RS by bidding in the RS market.

Bidding prices and stepwise bidding strategy
The EVAs will follow the optimised charging and discharging power quantities when bidding in the markets. Meanwhile, the bidding prices that are required for the success of market bidding also deserve careful formulation. The EVAs can bid the market cap prices for demand bids, or zero prices for supply offers and RS offers [27].
The single block bidding strategy is simple and easy to implement; however, it cannot adjust to the possible fluctuations of the market prices, and hence might not be economically efficient. The simplified bid price may create the following problems when facing fluctuating market prices (e.g. caused by renewable energy generation and forecasting errors): • If the market clearing price is significantly lower than the forecasted price of the aggregator and the energy originally scheduled at the same period is low to avoid the estimated high charging cost, only the conservative energy schedule attained by Model-D and Model-R will be cleared by the market, resulting in higher charging costs for the aggregator and its EVs. • If the market clearing price is significantly higher than the forecasted price of the aggregator and the energy originally scheduled at the same period is high to utilise the estimated low charging cost, the energy schedule attained by Model-D and Model-R will be cleared by the market, resulting in higher charging costs for the aggregator and its EVs.
In this study, a stepwise bidding strategy is adopted for the formulation of the bidding prices to enable these prices to accommodate unexpected market price fluctuations.

DA market bid prices
Model-D optimises the power quantities that should be accepted in the DA energy market. Nonetheless, the composition of both the DA bids and the RT bids represents the final bid quantities in the energy markets. Moreover, it could be more efficient if the bid quantities become even larger/smaller in the case of extremely low/high DA market prices. Assuming that the power quantities with positive values represent purchasing power from the market, the stepwise DA bidding function for each hour may comprise the following power blocks (PBs): , that is, the optimal DA market bids; , that is, the optimal bids of combined DA and RT market; , that is, the maximum charging power without reducing the optimal regulation capacity; , that is, the maximum charging power; , that is, the maximum discharging power without reducing the optimal regulation capacity; , that is, the maximum discharging power; PB7 0, that is, no charging/discharging bids.
To formulate an eligible stepwise bidding function, the following principle must be adhered to: Principle: The bid price of the x-th PB, denoted as PBx , must be higher than that of a PB with smaller power quantity and lower than that of a PB with larger power quantity. The PBs with the same quantities must share the same bid prices.
First, the price of PB1 should be settled. Here, two cases are developed as follows: Case I: PB1 ≥ 0, implying the optimised EVA schedule is to purchase power from the DA market. The bid price for PB1 should be higher than the highest possible DA MCP for acceptance. Thus, the bid price for PB1 can be formulated as follows: Case II: PB1 < 0, implying the optimised EVA schedule is to sell power to the DA market. The bid price for PB1 should be lower than the lowest possible DA MCP, as described by the following: The case when PB1 ≥ 0 is used as an example in formulating the prices of the remainder of PBs, and the case when PB1 < 0 can be handled analogically. PB2 t can be described as follows: , which is believed to not exceed the forecast minimum DA pricẽD A,min t based on (67), makes the accepted power quantity in the DA market PB2 on condition that the DA MCP is not higher than the PB2 t . Thus, the accepted power quantity in the DA market will automatically become larger to accommodate the low DA MCP, and the profit of the EVA is likely to increase in this manner. The case when PB2 < PB1 can be analysed similarly.
PB3 t will make sure that the EVA will not increase its DA market consumption to PB3 unless the DA market price becomes extremely low; thus, more profit for providing EV charging services can be guaranteed. Considering the daily EV charging requirement is relatively steady, the demand increase at a particular time slot will result in demand decline at other periods. Hence, the EVA has to evaluate the possible charging profit of the entire trading day before determining PB3 t as follows: where Δ inc denotes a predefined minimum bid price increment, whereas Λ() denotes the function that returns the average value of a series. PB4 denotes the maximum charging power of the EVA, where the optimised regulation capacity will be dismissed. PB4 t can be calculated as follows: wherẽR cap,max t ,̃R U,E,max t , and̃R D,E,max t denote the estimated maximum price for regulation capacity, regulation up, and down services, respectively. Term Δ 4 t denotes the revenue loss for reducing the regulation capacity from P * Rcap a,t to 0, which is calculated based on extra V2G cost, net energy consumption, and estimated maximum income for providing RS.
PB5 and PB6 both represent the discharging capability of the EVA. If the EVA discharges its batteries, it will cause the EVA extra cost to charge the batteries again. The extra cost for providing V2G can be calculated as follows: In this study, PB5 t , PB6 t and PB7 t are calculated as follows: t , and PB7 t can still be calculated with (67), (73)-(74), while the others can be calculated as follows: Moreover, the calculated bidding prices must be modified to accommodate the proposed principle; thus, the following constraints must be respected, ∀i, j ∈ {1,2,…,7}: The demonstration of the proposed bidding strategy when PB1 ≥ 0 is shown in Figure 2.

RT market bid prices
When the DA market clears at 16:00 on day D-1, the EVAs have another chance to modify their energy consumption schedules for day D in the RT market. In most cases, the EVAs are willing to strictly follow their optimal RT bidding results to prevent possible violations of EV charging requirements and constraints. Thus, the EVAs are assumed to submit single-block energy bids for RT market clearing, as they generally prefer not to change the bided quantities. The quantity of PB for RT market bidding, which is denoted as PBr, is equal to the optimised RT market schedule shown in PBr PBr

Charging power
Discharging power

Bid price
Bid price when PBr>0 Bid price when PBr<0

FIGURE 3
Bidding strategy for EVA in the RT market (81). The success of single block bids in the RT market will be guaranteed by the bid price as formulated in (82) as follows.
If PBr ≥ 0, a very high bid price, as shown in (82), makes sure that the bid block will be accepted in most of the time. If PBr < 0, a low bid price is adopted considering the cleared DA MCP, to enable the supply offer to be accepted in the RT market. From the proposed RT market bidding strategy, the EVAs are likely to win their desired bid while simultaneously avoiding losses caused by extreme RT MCP cases. The RT market bid curve is demonstrated in Figure 3.

RS market bid prices
Although the RT market is not cleared when deciding the RS market bids, the proposed RT market bidding strategy normally guarantees the acceptances of RT market bids. Therefore, the accepted DA and RT energy schedules will be considered fixed when bidding in the RS market. If the RS MCP becomes significantly higher or lower than those forecasted, the EVAs may get extra revenue by increasing or decreasing the optimised regulation capacity. Thus, a stepwise RS bidding function is also employed, which contains the following power quantities: , that is, the optimised regulation capacity; PBe2P Rcap,min , that is, the minimum eligible regulation capacity; PBe3 0, that is, EVA does not provide regulation service; PBe4P ⊕Rcap,max a,t , that is, the maximum regulation capacity that the EVA is able to provide.
The bid price for PBe1 should be such as to ensure that the optimised RS capacity is likely to be accepted, as shown in (83) as follows: where Δ e2 t denotes the revenue loss for reducing the P  The bid prices in the RS market must respect the principle as well, and the bid curve is shown in Figure 4.

EVA REAL-TIME DISPATCH MODEL
The EVAs submit their demand bids/supply offers to the ISO based on the proposed bidding strategy and the market timeline. The accepted bids/offers will be the guidelines for EVAs to dispatch the EV charging/discharging on operating day D. Ideally, the EVAs will simply follow their successful bids/offers when dispatching their EVs. However, the following issues may make the dispatch of EVs more complicated: i. The accepted bids/offers could not satisfy the EV charging requirements due to reasons such as less successful market bidding results, or a significant increase in EV charging requirement. ii. The available EV charging/discharging power could not meet the accepted energy quantities owing to unexpected EV driving patterns; for example, the number of available EVs is significantly lower than expected, thus making it impossible to consume the cleared power/provide the accepted regulation capacity. iii. The market clear results only define the total consumption/output of the EVAs, whereas the detailed dispatch schedules for individual EVs remain unsettled.
Thus, the EVAs need to allocate the accepted bids/offers to individual EVs according to their actual driving behaviours and make necessary decisions in case the accepted bids/offers or the EV charging requirements cannot be fulfilled.

Hourly EV dispatch model
The EV parking probability is employed in optimising the bidding model, but it does not contain information concerning the individual EV parking statuses, that is, which EV is parked or not at a specific time. Thus, the EVAs are not able to dispatch the charging/discharging schedule of each EV without actual driving data, which can only be obtained in real-time on day D.
For simplicity, the EV parking/driving statuses of all the EVs are assumed to be constant within each operating hour. For example, an EV that is parked at the beginning of hour H will remain parked until the end of H, and it may not drive away from home/parking slot before the beginning of H+1. The accuracy of EV real-time dispatch can be enhanced by similarly dividing each hour into multiple time slots. Without loss of generality, the detailed EV charging/discharging arrangements will be dispatched hourly in this study, and dispatch models with smaller time slots can be derived correspondingly.
The dispatch of individual EVs at hour H is done at the beginning of H, and the objective function of the a-th EVA with a set of EVs will aim at following the cleared bids/offers with its best effort, while voiding unnecessary penalties. The hourly EV dispatch model for hour H is shown as follows.
The objective function (92) minimises the possible penalties that an EVA may face throughout the operating day D. In (92), a weighted sum of EV charging power is introduced and the weights of EV charging powers will help the EVA to decide EV charging control schemes. Terms K a,e,k and P K a,e,k denote the weight and power quantity of power block k of the eth EV, respectively. Terms Ω E(a) and Ω K stand for the set of EVs within the control of the a-th aggregator and the weights, respectively. Constraint (93) denotes the charging/discharging power limits of a single EV, whereas p ev a,e,H and pk a,e,H denote the energy schedule and parking status of the e-th EV, respectively. Constraints (94)-(95) denote the SoC limits for the parked EVs. The total EV charging/discharging power of the a-th aggregator, denoted as P ev,act a,H , is calculated from (96). The EVA will attempt to keep the total power to the accepted DA and RT energy bids, as demonstrated in (97)  is positive, the EVA will have to pay for the energy mismatch penalty as discussed in Equation (12) in Section 3.2. Furthermore, the EVA will also pay/get paid for purchasing/selling power based on the RT MCP.
Equations (99)-(107) calculate the SoC levels of all the EVs at hour H concerning the EV charging model described in (1)-(6).
In (99), ΔS C+ a,e,H and ΔS C− a,e,H respectively denote the SoC increment/decrement if the EV is parked, whereas ΔS C,trip a,e,H denotes the SoC decline when the EV is not parked. If an EV is parked and available to be dispatched by the EVA, its SoC profile will be managed by the EVA through charging/discharging controls. Meanwhile, if an EV is not parked, its SoC will drop along with its mileage, and the entire process may take some time until the end of the trip. In this study, the EV SoC level is assumed to be uniformly decreasing during the non-parking period, whereas the total SoC drop is derived from the EV driving pattern.
Equations (108) Δp ev a,e,H,1 represents the power quantity that must be charged at this hour; otherwise, the EV will be undercharged at the agreed time even if it is later charged at its full charging power. Consequently, the EVA is subjected to the penalty for undercharging its EV if the charging power of the e-th EV is lower than Δp ev a,e,H,1 , as demonstrated in (111). Δp ev a,e,H,2 represents the average required charging power to meet the EV charging agreement on condition that the EV remains parked until the agreed time. Δp ev a,e,H,3 is the remainder of its maximum charging power. Based on these three divided parts, the EVAs will be able to appropriately allocate their successful bids to individual EVs. The first part must be satisfied to avoid EV undercharge penalties. The second part is the desired amount of EV charging power such that the EVA will not have to fully charge that battery and make it available for regulation services during the subsequent hours. The EVA will decide to charge/discharge the third part of EV charging power at its convenience. In case the accepted energy consumption is not sufficient to accommodate all the power parts with the same priority, for example, the sum of the first part charging power of all the idle EVs is greater than the successful energy bids in the energy markets, the EVA will allocate its purchased energy considering the weights of each charging power block, which can be calculated as follows: Generally, the EVs with lower SoC levels enjoy higher charging priorities. During the early morning, the weights of EV charging power blocks are calculated based on their SoC levels and the remaining charging time (assuming they will stay parked) before their predefined check time. Afterwards, both the first and the second parts are equal to zero, and the weights will depend on their SoC levels only. In summary, the EVA's hourly dispatch model, denoted as Model-H, can be described as follows:

Revision of regulation capacities
The proposed market frame enables the EVAs to reduce their regulation capacities at least an hour before the operating time. Regulation capacity reduction for H+1 can be calculated with the forecast EV driving patterns and the real-time EVA dispatch strategy at H. The following four steps will be employed: Step 1: Evaluate the initial EV SoC level at the end of period H, which is obtained from the EV real-time dispatch model (Model-H), denoted as S ⊕C a,e,H . S ⊕C a,e,H , will be used as the initial EV SoC to estimate the regulation capacity revision of the e-th EV for the period H+1.
Step 2. Estimate the EV parking statuses at H+1. The number of parked EVs at H+1 was originally estimated with the EV parking probability. Here, a modification is adopted considering the forecast error at H as follows: pk H +1 =̃p k,ori wherep ev,up a,e,H +1 andp ev,dw a,e,H +1 estimate the maximum power output/input of the e-th parked EV, respectively. Term̃p k a,e,H +1 denotes the estimated parking status of the e-th EV at H+1 based on step 2. If the number of parked EVs is believed to increase from H to H+1, the SoC profiles of the newly arrived EVs are unknown to the EVA. This can be handled conservatively by ignoring their capacities in providing regulation up service, whereas their full charging power is eligible for the regulation down service.
Step 4: The quantity of regulation capacity that should be modified, which is denoted asP or take risks by submitting modification bids with less or no decreased regulation capacities. In this study, the EVAs are assumed to be risk-averse entities, and they will submit the revision bids based on the calculatedP Rcap,mod a,H +1 .

5.3
Control strategy for regulation services During the RT system operation, the accepted RS suppliers will be frequently called to provide regulation up/down services by the ISO within every hour. Once an EVA receives a regulation signal, it is supposed to immediately reduce/increase the total charging power by the accepted RS capacity to fulfill the regulation up/down services. The developed weights of the EV charging power blocks also help the EVA to decide the EV to be controlled to provide the RS. The following procedures are employed in this study: i. The RS signal from the ISO is received, and the parked EVs screened. If the regulation down service is desired, then proceed to step ii, otherwise jump to step iii. ii. The EVA will increase its total charging power for the regulation down service. Thus, Δp ev a,e,H,1 has the highest priorities, whereas Δp ev a,e,H,3 enjoys the lowest. The PBs with lower priorities will not be considered unless the full charge of all the PBs with higher priorities fails to fulfill the accepted regulation capacity. If the full charge of all the PBs exceeds the accepted regulation capacity, those PBs with the lowest priorities and lowest charging weights will be dismissed for the RS, until the exact regulation capacity is achieved. iii. The EVA will decrease its total charging power, or even increase its total discharging power, for the regulation up service. The decision-making procedure is similar to that in step ii, but the EVA prefers to discharge the PBs with lower priorities and lower weights. If the discharge of all the charging PBs is not able to fulfill the accepted regulation capacity, the EVA will discharge some EV batteries. The batteries with the highest SoC levels will be adopted for discharging until the exact regulation capacity is achieved.
Through the proposed procedure, the EVA will quickly determine the detailed schemes for the RS. The flowchart of the proposed procedure is presented in Figure 5.

Regulation signal data set
Set size: Data element: 1. Regulation up data 2. Regulation down data

FIGURE 6
Simulation scenario data sets

CASE STUDY SETUP
The proposed bidding strategy is expected to formulate risktolerant solutions for the EVAs without perfect forecasts of electricity market prices and EV driving behaviours. Once a bidding strategy is decided, different market prices and EV driving patterns will lead to various market clearing and EV dispatch results. Therefore, a scenario-based simulation is adopted to verify the effectiveness of the proposed strategy.

Scenario design
The effectiveness of the proposed bidding strategy can be validated by changing the electricity market prices, regulation signals, and individual EV driving information, as demonstrated in Figure 6. The three data sets will be adopted to create several scenarios for verifying the proposed strategy. Electricity market prices dataset: Each data element contains three price series that cover a unique operating day, namely, the DA, RT, and RS market price data. Theoretically, a data ele-ment with price data from different operating days can be used to generate test scenarios, but the correlations among the three markets are not guaranteed.
Thus, to create a data set with N mp elements, market prices of N mp days will be used. In this study, the DA, RT, and RS market price data were obtained from the real PJM energy market data of April 2015, where N mp = 30 [36]. The RS market data represents the regulation capacity price, whereas the energy prices for regulation up and down services are assumed to follow two normal distribution N($90/MW, ($10/MW) 2 ) and N($15/MW, ($6/MW) 2 ), respectively.
Regulation signal data set: The real regulation signals in the PJM regulation market from 4 May to 10 May 2014 are adopted as the test regulation signal data in this study [36]; hence, N rs = 7. The time resolution of the regulation signal is 4 s; thus, there are 21,600 regulation signal data in each day.
EV driving data set: The initial SoC levels of all the EVs are assumed to follow the normal distribution of N (0.45, 0.05 2 ), and 0.8 is assumed to be their contracted SoC level with their aggregators. In terms of the contract delivery time (T con e ), half of all the considered EVs expect that charging finishes at 6:00 each day. The others are divided evenly into two groups that expect the charging to finish at 7:00 and 8:00 each day. The daily driving mileage of EVs is assumed to follow the normal distribution of N (60 mi, (5 mi) 2 ).
The EV daily driving/parking status data can be generated based on the data in [40] and [45]. Detailed EV driving/parking statuses are attained through the following steps: i. The EV parking statuses from 0:00 to 1:00 will be generated initially. A randomly generated number ranging from 0 to 1 is allocated to each EV. If the number is lower than the home parking probability from 0:00 to 1:00, the EV is assumed to be parked at home during that period; otherwise, the EV is not parked at home and, thus, is unavailable to the aggregators. ii. For the remainder of the operating hours of the day, N ev random numbers ranged from 0 to 1 will be gener- In this study, it is specified that N ev = 200 and N dr = 50, indicating that 50 sets of EV driving data will be generated based on the proposed steps and each set of data contains detailed driving information of 200 EVs.
For each optimised EV bidding strategy, a total of N mp × N rs × N dr scenarios, that is, 30 × 7 × 50 scenarios are employed in the scenario-based simulation. The size of the test scenarios is suitable for demonstrating the effectiveness of the proposed strategy.

Estimated data for robust optimisation
The EV aggregators have no access to the real market prices, regulation signals, and EV driving data generated for scenariobased simulation. Their bidding strategy will be built on the estimation of these parameters.

Electricity market prices
The average DA, RT, and RS market prices of the previous market price data set is assumed to be the aggregators' price prediction. Moreover, an estimation interval is added to the forecasted price at each hour to accommodate the robust bidding   Figure 7. The estimation accuracy is summarised in Table 1, where 'underestimate' and 'overestimate' represent the cases where the real price is higher or lower than the forecasted interval, respectively. Generally, the level of confidence is approximately 90%. The estimated regulation up/down energy prices are assumed to be fixed values of $100/MWh and $15/MWh throughout the day, and the estimated prices intervals are ($75/MWh, $125/MWh) and ($5/MWh, $25/MWh), respectively.

Regulation signals
The regulation signals determine the aggregators' actual regulation schedules, but it is impossible to estimate the precise regulation signal series after every 4 s. Instead, the called probabilities RU t and RD t are set to 25% and 20%, respectively, in this study as the aggregators' estimation. Notably, the regulation down service is called more frequently than the regulation up service based on the aforementioned regulation signal data.

EV driving data
Considering the EVA responsibility discussed in Section 3 and historical EV charging profiles, certain checkpoints are necessary to allocate the daily purchasing schedule. If the total daily energy consumption is 1, and R ev,cr a,t denotes the cumulative electricity consumption of aggregator a by the end of t. The adopted checkpoints are summarised in Table 2.  The checkpoints 13:00 and 19:00 are obtained from historical EV data to accommodate possible charging requirements during the day. These checkpoints can be omitted if required. In this study, the energy consumption rate of EV batteries is assumed to be 0.2 kWh/mi.

Parameters
Some parameters are employed in optimising the EV bidding strategy and the scenario-based simulation. GT + loss is set to $20/MWh, and the values of other parameters are summarised in Table 3. The presented models (i.e. Model-D, Model-R, and Model-H) are all solved by AMPL/CPLEX [46,47].

Market bidding
The EVA evaluates the forecasts and historical data, formulates the optimal energy schedule, and bids in the DA market. Table 4 demonstrates the impact of robust parameter Γ DA on the EVA's energy schedule and expected revenue, whereas denotes the constraints violation probability. As the RT market price is more likely to fluctuate severely compared to the DA market price, the EVA will reduce its involvement in the RT market to avoid risks as Γ DA increases. The expected revenue also decreases if the EVA requires lower . In this case study, five EVAs that represent different risk preferences are selected to validate the feasibility of the proposed strategy. The preferred of EVA1, EVA2, EVA3, EVA4, and EVA5 are assumed to be 20%, 10%, 5%, 1%, and 0.01%, respectively. The optimised DA market schedule and the final  Figure 8. Here, all the EVAs will not discharge electricity to the power system in the combined DA and RT energy market. However, the EVAs will prefer to discharge their batteries at certain periods in the DA market while purchasing more energy in the RT market instead. Thus, they may gain extra profits owing to the price differences between the DA and RT markets. Moreover, the EVAs may also schedule a 0.1 MW charging power in the energy market in the afternoon to accommodate the minimum RS capacity requirement of the ISO. Otherwise, the EVAs are not eligible to bid in the RS market during those periods, which is believed to reduce their revenue from the RS market. EVAs formulate their bids in the DA market based on the proposed stepwise bidding strategy. EVA4, for example, bids in the DA market of the next operating day as shown in Figure 9. The bid blocks of PB1-PB7 are depicted in red, green, yellow, blue, cyan, pink, and grey, respectively. Taking the market bidding curve at the first time slot as an example, the bid prices for non-positive charging power blocks PB5-PB7 (i.e. cyan, pink, and grey blocks) are very high (more than $200/MWh), indicating that this EVA will not act as a power producer unless the market clearing price is higher than the cost associated with battery discharging. Meanwhile, the EV charging power blocks PB3 and PB4 are low or even negative bid prices, indicating that the EVA does not prefer to increase power procurement unless the market price can compensate for the opportunity cost and the reduction in ancillary service market rewards. Bidding blocks PB1 and PB2 represent actual energy demands of the EVA, and thus their bid prices are in the normal range of energy market price time series. Hence, the proposed stepwise bidding strategy enables the EVAs to automatically adjust their EV charging schedules with varying market clearing prices.
The RT and RS market bidding strategies are based on the cleared results of the DA market, which can be obtained through scenario-based simulation in the latter sections. As a demonstration of the proposed bidding strategy, the PJM DA market price on 1 April 2015 is adopted, and the RT and RS market bids will be formulated based on the EVA's cleared DA market results. Figure 10 shows the bids of EVA4 in the RT and RS markets based on its cleared DA market results, and the acceptable in Model-R is also set to 1%. The black bars in Figure 10 represent the RT market bids, whereas the grey bars represent the RS capacity market bids.

Scenario-based simulation
In the scenario-based simulation, the cleared EVA energy market and RS market schedules are obtained by combining the  EVA bids and real electricity market price data. Thus, the bids of each EVA will have N mp cleared schedules. Afterwards, each cleared EVA schedule will be dispatched based on regulation signal and EV driving data, where N rs × N dr detailed EV charging/discharging schedules are formulated. Table 5 compares the market cleared schedules and the EVAs' optimal schedules. Case I indicates that the market cleared schedule is the same as the optimised schedule of the EVA, whereas cases II and III indicate that the cleared schedule is smaller or larger than the optimised schedule, respectively.  Considering the market price forecast accuracy, the proposed multi-market bidding strategy is effective. The RT market bids generally suffer from the lowest bidding accuracy owing to the highly volatile RT market price test data. Moreover, the average errors between optimised schedules and cleared results in the RT market are significantly higher than those of DA and RS market, and the error drops gradually with the decrease of . On the operating day D, the EVAs will formulate their hourly energy schedules based on the cleared market results and the proposed real-time dispatch model (Model-H), as they have access to the real-time EV driving data. Moreover, the EVAs also modify the original schedules and dispatch their available EVs based on the regulation signals from the ISO. The originally scheduled and dispatched hourly energy profiles of EVA4 are shown in Figure 11. The dispatched profile has a time resolution of 4 s, and the hourly data is the average consumption at that particular hour. As can be observed from Figure 11, the actual average energy consumptions are generally slightly higher than the original schedules, as the regulation down service is more likely to be called based on the adopted RS signals.
A detailed 4-s dispatch result of EVA4 with market price data and the RS signal is shown in Figure 12, including the average EV SoC profiles.
The cleared RS capacities may not be fulfilled at specific periods owing to the unpredictable EV driving behaviours. Table 6 FIGURE 12 Detailed EVA charging/discharging behaviour and average SoC profile based on the selected scenario lists some results about RS capacity modification. In Table 6, P mis and P mis denote the average and maximum mismatched RS capacity, whereas R mis and R mis denote the average and maximum mismatched RS capacity while considering the RS capacity modification, respectively. Revising the RS capacity before the operating hour will reduce P mis and possibly reduce the penalties for RS capacity mismatches. However, as the EVA cannot perfectly predict the driving patterns of all the EVs, RS capacity mismatches may still exist.

Profit and risk analysis
From the perspective of an EVA, the profit and risk for participating in the electricity markets are the primary concerns. The economic statistics of the tested EVAs are listed in Table 7.
Introduced as indexes to evaluate risks [48], max and mean denote the maximum deviation and standard deviation of EVA's R profit a obtained from test scenarios, respectively. Generally, the average profit of an EVA decreases slightly by reducing while the EVA enjoys a less fluctuating profit and alleviates risks. The penalties cannot be completely avoided, as errors always exist in market price and EV driving behaviour forecasts. However, the results indicate that the EVAs will alleviate some of these penalties by employing a more riskconservative bidding strategy. The profits in Table 7 are generally higher than their corresponding values in Table 4. This can be explained by (i) the conservativeness of robust optimisation, and (ii) the market price forecast error.
According to Table 7, it is difficult to conclude whether a universal optimal risk preference for the EVAs to bid in the multimarket exists. Based on the proposed test scenarios, the riskseeking strategy seems to outperform the risk aversion strategies, as the maximum and minimum profits of EVA1 are both higher than that of the other EVAs. However, the profits of EVA1 also suffer from the most severe fluctuations, especially when the market prices become volatile. Furthermore, the performances of EVA1 may not always be the best if other market prices and EV driving data are utilised. Thus, it is proposed that the EVAs adjust their risk settings to accommodate the volatile market price profiles. Additionally, improving the market price forecast accuracy will certainly help the EVAs when making market bidding decisions.

CONCLUSIONS
The participation of EVAs in energy and ancillary service markets is addressed in four phases. The first phase discusses the electricity market framework and the interactions between the EVA and EV customers, which are the foundations for developing bidding strategies. The second phase proposes a robust optimisation model for EVAs to schedule their consumption in multiple markets. The third phase proposes the stepwise bidding strategies based on the optimised schedules to incorporate with possible market price fluctuations. The last phase develops a real-time energy allocation strategy for EVAs to follow their successful market bids and regulation signals. All aspects pertaining to EVA participation in multiple kinds of markets are covered by these models. The effectiveness of the proposed strategies is demonstrated through scenario-based simulations. Some test scenarios are based on actual electricity market data, and simulation results indicate that the proposed bidding strategy is feasible for EVAs with different risk preferences for optimising their bid decisions in multiple electricity markets. The influence of unforeseeable EV driving behaviours and the regulation orders can be handled efficiently by the EVA dispatch model (Model-H). As is expected, case studies demonstrate that participation in the RS market could generate considerable and reliable revenues for the EVAs. Thus, coordinated optimisation of participation strategies for each EVA in the energy markets and RS market represents a promising method of increasing revenues.
The developed robust bidding strategy can accommodate risk preference of various EVAs; however, the forecast errors of market prices and stochastic EV driving behaviours still make the profit of each EVA suffer uncertainty to some extent. Future studies will focus on the strategies of stabilising the revenue of each EVA in participating in electricity market operation as follows: • considering the financial market and long-term contract for managing risks of volatile prices in energy markets and ancillary service markets; • developing a bidirectional communication interface and a control structure for EVAs and EV owners, thereby enhancing the dispatching accuracy of EVAs and the satisfaction degrees of EV owners; and • combining EVAs with other energy sources, such as wind farms, to form a virtual power plant or a market alliance and utilise the characteristics of each participant to pursue a more advantageous outcome in multiple types of electricity markets.

ACKNOWLEDGEMENT
This work is supported by National Natural Science Foundation of China (U1910216).