Electric Vehicle Charging and Discharging Coordination on Distribution Network Using Multi-Objective Particle Swarm Optimization and Fuzzy Decision Making

This paper proposed a optimal strategy for coordinated operation of electric vehicles (EVs) charging and discharging with wind-thermal system. By aggregating a large number of EVs, the huge total battery capacity is sufficient to stabilize the disturbance of the transmission grid. Hence, a dynamic environmental dispatch model which coordinates a cluster of charging and discharging controllable EV units with wind farms and thermal plants is proposed. A multi-objective particle swarm optimization (MOPSO) algorithm and a fuzzy decision maker are put forward for the simultaneous optimization of grid operating cost, CO2 emissions, wind curtailment, and EV users’ cost. Simulations are done in a 30 node system containing three traditional thermal plants, two carbon capture and storage (CCS) thermal plants, two wind farms, and six EV aggregations. Contrast of strategies under different EV charging/discharging price is also discussed. The results are presented to prove the effectiveness of the proposed strategy.


Introduction
Electric vehicles (EVs) will play a vital role in the future's sustainable transportation systems, since this technology is promising for environment, energy security, and improved fuel economy. Certain issues will need to be addressed in the event that the number of EVs on the road increases. One vital issue is the method by which these vehicles will be charged and if today's grid can sustain the increased demand due to more EVs. Although EVs' growing energy demand seems to be a heavy burden to the power gird, they can actually benefit the grid if we control the charging and discharging behavior of them properly. For example, nowadays, a growing quantity of renewable energy, such as wind and solar generation, is connected to the grid [1]. Due to the discontinuity of the renewable energy power generation, an energy storage system is needed to assist the grid to absorb the volatility of renewable energy power. Electric vehicles are considered to be the energy storage device with the most potential to absorb the volatility of renewable energy [2]. When a large number of EVs are aggregated, the huge total battery capacity is sufficient to stabilize the disturbance of the grid caused by the distributed power grid interconnection [3].
In earlier research, Amory Lovins first proposed the concept "V2G" in 1995, which is an abbreviation of "vehicle to grid", meaning EVs serve in discharge mode to support the grid. This

System Concept
In order to better schedule the charging and discharging behavior of EVs, since the capacity of a single EV's battery is hard to have a measurable influence on the transmission grid. The concept "EV aggregator" is present in this paper to represent a number of grid-connected controllable EVs. By aggregating EVs, a conceptual framework for a thermal-wind-EV system is built up as in Figure 1 shows: In this conceptual system, at the top layer, the supervisory control and data acquisition system (SCADA) system is used for a double purpose: one is to collect information, such as EV battery status, EV number, predicted wind power, predicted load, etc., which are feedback from EV aggregator, wind farms, and thermal plants; another is to make a schedule for the above units to generate/absorb power to/from the transmission grid.
At the middle layer, wind farms, thermal plants, and EV aggregator receive a working schedule transferred from the top layer and make a detailed work plan for each wind turbine, thermal power generating unit, and EV unit inside the system. The middle layer then send out power instructions to dispatch the input or output power of the aforementioned power units.
At the bottom layer, while a large number of EVs are connected to the transmission grid via power electronics devices, Each EV unit follows the power instruction from the middle layer and feeds its current status back constantly to the aggregator so that the system is able to adjust its power distribution strategy according to the real-time model.

Electric Vehicles Aggregator Model
The electric quantity equation of EV i is given as follows: where E i,k is the electric quantity of EV i at the time period k. d b is battery self-discharge rate, η C and η D are the battery charging and discharging efficiency factor, respectively, while P C i,k and P D i,k are the charging and discharging power of EV i at time period k, respectively. ∆t is the time step, given in hour.
The cost equation of the EV aggregator is as follows: where γ C i,k is the EV's charging price, while γ D i,k is the discharging cost for V2G services. Equations (1)-(3) are subject to constraints as follows: (1) Power constraints. P where P (C,min) i,k and P (C,max) i,k are the lower and upper limits of EV's charging power P C i,k respectively, P (D,min) i,k and P (D,max) i,k are the lower and upper limits of EV's discharging power P D i,k , respectively. (2) Electric quantity constraints.
E min i,k and E max i,k are the lower and upper limits of E i,k , which are decided by EV's battery capacity B and its depth of discharge (DoD): (3) EV will not charge and discharge at the same time period.
where SOC i,k is the battery state of charge of EV i , which is decided by the equation: SOC D i,d min is the lower bound of state of charge (SOC) when EV is discharging, while SOC C i,c max is the upper bound of SOC when EV is charging. These bounds are set by EV users.

Thermal Plant Model
In this thermal plant model, we consider two kinds of thermal generating plants, including traditional coal-fired (TC) power plants and coal-fired power plants with carbon capture and storage (CCS). The consumption function of traditional coal-fired power plants can be approximated by quadratic function: where P U i,k is the active power output of TC power plant i at time period k, while a i , b i , and c i are consumption characteristic factors of the function.
The emission equation of traditional coal-fired power plant is as follow: where α i and β i are emission characteristic factors of the function. CCS power plants need to consume a large amount of energy to capture CO 2 , which leads to an increase in generating cost. Assuming that the capture rate of CO 2 is ω, then the consumption function F C i and emission function C C i of CCS power plants are as follows: where P ccs i,k is the active power output of CCS power plant i at time period k. The amount of captured CO 2 is calculated as follows: Equations (12) - (16) are subject to constraints as follows: (1) Power constraints P where P For the safety of thermal plant operation, the power output of TC plants and CCS plants must be subject to ramping rate constraints as follows: where U G and D G are the highest increase and decrease value of the generating unit power output between adjacent time period.

Wind Farm Model
When wind power generation is allowed to be integrated into the grid, the power output of thermal plants decreased if the wind power generation increases. Assuming that the upper bound of wind generation is equal to the negative peak load regulation ability of thermal plants [25]: where P k min,i is the lower limit of thermal generating units, which includes TC and CCS units; P k L is the system load; P k loss is the system loss; P C i,k and P D i,k are EV i 's charging and discharging power, respectively. Consequently, the generating capacity of wind farm i at time period k can be calculated as follows: where P k wa,i is the available generating capacity, P k wc,i is the forecasted wind generation capacity, and P wr,i is the rated capacity that of Wind farm i.
The cost function of wind farms is as follows: where ε i is the operating cost factor of wind farms. Equations (20) - (22) are subject to constraint as follows:

Multi-Objective Optimization Model
According to the aforementioned EV aggregator model, thermal plant model, and wind farm model, the EV-wind-thermal coordination multi-objective optimization problem can be formulated as follows: where T is the scheduling time; N G is the number of TC units; N CC is the amount of CCS units; N W is the number of wind farms; N V is the total number of electric vehicles; π e i is the emission costs; π s i is the carbon capture costs. J2 = Minimize the global CO 2 emission J4 = Minimize the EV user s cost where ζ is the batteries degradation cost of EV users, which reflects the influence of deep charge and discharge on the EV's battery life. Equations (24)- (27) are subject to constraints as follows: (1) Power balance In Equation (28), power balance does not takes the power losses into account. The calculation of power losses is an open question for later research.
(2) Up-and down-spinning reserve constraints [26,27] where L% is the load reserve coefficient, which presents the uncertainty of load prediction, while w% confirm italics is the wind reserve coefficient, which presents the uncertainty of wind prediction.

Multi-Objective Particle Swarm Optimization Algorithm and Fuzzy Decision Maker
A multi-objective optimization problem (MOP) [18][19][20] can be described as follows: where f i (·) is objective function, k is the number of objective functions, and x = [x 1 ,x 2 , · · ·, x n ] T is the vector of independent variables. The aim of the multi-objective optimization is to determine a particle set of values x * = [x * 1 ,x * 2 , · · ·, x * n ] T from the feasible solution space which is closest to the ideal value of all the objective functions.
In single objective optimization problem (SOP), the global best G best is the solution of the optimal problem. However, in MOP, it is impossible to find a single G best that would be optimal for all objective functions simultaneously. Hence, a set of compromised solutions called "Pareto front", which consists of a set of solutions wherein no existing dominated solution, is utilised to solve the MOP.
Since the Pareto front usually contains a large amount of individual solutions, but a real practical problem only needs one solution, it is essential to choose a single best compromise solution out of the Pareto front by using decision maker (DM).
Based on the above analysis, the process of solving a MOP can be divided into following two main steps: (1) Find out the Pareto front of the MOP (2) Extract out the best compromise solution from the Pareto front

Multi-Objective Particle Swarm Optimization Algorithm
Particle swarm optimisation (PSO) [19] is an evolutionary soft computational optimisation technique developed by Kennedy and Eberhart. The system initially has a population of randomly-selected solutions. Each solution, which is called a particle, is given a random velocity and is flown through the problem space. The particles keep track of their best position and their corresponding fitness values. Among these best positions, the particle with the best fitness is the G best of the swarm.
In MOPSO [19][20][21], a set of particles are initialised in the D dimension decision space at random. For each particle i, a position x i = [x i1 , x i2 , · · ·, x id ] T in the decision space and a velocity v i = [v i1 , v i2 , · · ·, v id ] T are assigned. The particles change their positions and move towards the so-far best-found position p i = [p i1 , p i2 , · · ·, p id ] T . Besides, the so-far best-found position of the whole particles set p g = [p g1 , p g2 , · · ·, p gd ] T is kept. for every single particle, its evolution function in the D dimension decision space is as follows. The so-far best-found non-dominated solutions from the last generations are kept in an archive, which is an external population.
where c 1 and c 2 are the personal learning coefficient and global learning coefficient of the equation, respectively, by which the particle is able to make self summary and learn from outstanding individuals inside the whole particles set. ω is an inertia weight, which reflects the impact of a particle's formal speed to its current speed. In this paper, a self-adaptive inertia weight tuning method is adopted to enhance the searching ability and convergence speed of the optimal algorithm: where ω max is the upper limit of inertia weight, ω min is the lower limit of inertia weight, t max is the maximum iterations. The pseudo code of the proposed MOPSO algorithm is shown as follows: for i = 1 to N do 4: Individual best position p id update

Fuzzy Optimality Decision Making
In the real world, people prefer to describe the importance of an object and make decisions based on concepts such as "very important" or "less important" rather than using precise numbers. It is considered that, in the real world, even the most experienced expert can't precisely judge the importance of different objects by using certain numbers such as "0.7" or "0.5" [22]. It is considered that fuzzy sets or fuzzy numbers [22] are more appropriate than numbers as weight value used to describe the importance of different objects [25]. Hence, in this paper, a FODM approach is present to extract out the best compromised solution from the Pareto front, where the individual of the highest fuzzy optimality is the best compromised solution of the MOP.
In this paper, a triangular fuzzy numberÃ can be defined asÃ = (l, m, r), the membership function of which is defined as: where l and r are the lower and upper limits of the fuzzy numberÃ, and m is the modal value. The membership function of the triangular fuzzy number (Ã 1 ,Ã 2 , · · ·,Ã n ) can be illustrated as in Figure 2.   Where the red thick lines in Figure 3 represent the membership function of the maximum fuzzy subset max betweenÃ 1 andÃ 2 . DefineÃ max = max(Ã 1 ,Ã 2 , · · ·,Ã n ) as the maximum fuzzy subset among (Ã 1 ,Ã 2 , · · ·,Ã n ), the membership function ofÃ max is defined as: Hamming distance [24] d H (Ã i ,Ã j ) is used here to measure the distance between two fuzzy numbersÃ i andÃ j : The hamming distance d H (Ã 1 , max) and d H (Ã 2 , max) are illustrated in Figure 3, where s 1 ,s 2 ,s 3 , and s 4 are the geometric metric of hamming distance. d H (Ã 1 , max) = s 2 + s 3 , d H (Ã 2 , max) = s 1 + s 4 .
In this paper, the pseudo code of the proposed FODM algorithm is shown as follows: ,M k = max(r 1k ,r 2k , · · · ,r mk ) 6: Compute hamming distance: 7: Sort D i from small to large, the smallest is the best compromise solution.

Case Study
In this work, a 30 node grid system as shown in Figure 3 is adopted here to verify the effectiveness of the proposed optimal strategy, which contains two CCS thermal plants, three TC thermal plants, two wind farms, and six EV aggregators. Where the capacity of TC thermal plants is 150 MW, the capacity of CCS thermal plant is 150 MW, the capacity of a wind farm is 100 × 2 MW. There are 1000 EV in each EV aggregator, while the capacity of each EV's battery is 200 kWh. The upper limit of charging power is 40 kW, while the upper limit of discharging power is also 40 kW. Assuming the SOC of the 1000 EV is in accord with Gaussian distribution, set SOC D i,d min = SOC C i,c max = 50%. Assuming EV's charging price γ C i,k = 0.1 $/kWh, EV's discharging cost γ D i,k = 0.2 $/kWh, EV's battery degradation cost ζ = 0.13 $/kWh [2]. Simulation Parameters of the used thermal plants are listed as shown in Table 1    In this paper, the wind power model is predicted as Figure 4a shows, while the wind farm operating cost factor ε i is 50 $/MW, and the reserve coefficient of wind turbines, w%, is 20%.
In order to reduce the risk of incorporating uncertain wind forecasts into system scheduling, wind uncertainty should be taken into consideration before conducting further calculation. It is assumed that the wind forecast error is likely normally distributed [28]. The error of wind generation forecast is referred to asẽ. The forecasting error is estimated with a level of confidence α% [29], which means the probability of forecasting error being greater or equal toẽ is less than (100 − α)%. The wind generation capacity counted in the optimal schedule is calculated by: where P k wp,i is the predicted wind power generation of wind farm i.
Since we are concerned more with overestimation of the wind power (or that the power supply might not satisfy the demand), a one-side distribution curve is considered. By specifying the level of confidence, the value ofẽ can be estimated as follows: where µ e and σ e are the estimated mean value and standard deviation from sampling of the error of the historical forecast, respectively. The value of z α for 90%, 95%, and 99% confidence level can be chosen by Table 2 [29]. In the simulation setup of this paper, the confidence level is selected as 95%, the error mean µ e is 5 MW, and the standard deviation σ e is 3 MW.
Assuming that the system load model is as Figure 4b shows while there are no EVs in the EV aggregator. Assuming that the reserve coefficient of system load L% is 10%.  MOPSO parameters are initialized as follows: both the maximum number of iterations and population size are set to 200, while repository size is set to 50, inertia weight ω is 0.5, the personal learning coefficient c 1 is 1, while the global learning coefficient c 2 is 1.

Results and Analysis
According to the above curve and installment, the data of system operating status at 18:00, 1:00, and 4:00 are extracted out respectively to calculate the result of the simulation case. The cost of TC plants and CCS plants is calculated by using Equations (12) and (14), respectively, the emission of TC plants and CCS is calculated by using Equations (13) and (15). Equation (16) is used to calculate the amount of captured amount of CO 2 of CCS plants. Equation (22) is used to calculate the wind farms' cost, while EV's charging and discharging cost is calculated by using Equations (2) and (3). Equations (24)- (27) are used to calculate the global cost, CO 2 emission, wind curtailment, and EV users' cost, respectively. These equations are subject to the aforementioned constraints.
The Pareto front of the three scenes are shown as shown in Figure 5.  Setting up the fuzzy weight of the four objective functions by using triangle fuzzy members "Middle" (ω 1 = (0.3, 0.4, 0.5) ), "High" (ω 2 = (0.4, 0.5, 0.6)), "High" (ω 3 = (0.4, 0.5, 0.6)) and "Middle" (ω 3 = (0.3, 0.4, 0.5)) (the way to calculate it is shown as pseudo code in Section 4.2) respectively, the individual of the highest fuzzy optimality which is also the best compromised solution is selected from the Pareto front as Table 3 shows. By contrasting Figure 5 and Table 3, it can be seen that: From the above case it can be seen that when the system load is relatively high, for example at 18:00 (due to the fact that wind energy is relatively low by that time), the thermal plants need to go all-out to satisfy the power needs of the load in the system. However, when the EV aggregator is controlled to discharge its power to the grid, it is capable of sharing part of the thermal plants' work, hence the CO 2 emission can be relatively cut down. Also, EV users can receive economy benefit from sending their EVs to attend the V2G activity at the same time.
As a contrast, when the system load is relatively low, for example in the 1:00 scenario (due to the relatively high availability of wind energy at that time), while the thermal plants still need to maintain spinning, a large amount of wind energy will be abandoned for keeping the power balance of the grid. However, this power was saved because the EV aggregator is controlled to charge from the grid by the time. Hence, the wind curtailment of the system is reduced. Table 4 shows the result of a contrast simulation of two strategies under different discharging prices in the 18:00 scenario. One does not take the EV user's cost into consideration (strategy A), while another is the proposed strategy in this paper (strategy B). Table 5 shows the result of a contrast simulation of strategies A and the proposed strategy in this paper (strategy B) under different charging prices in the 4:00 scenario.  From Table 4 it can be seen that when the EV users' cost is not considered, EVs discharge more energy when the discharging price γ D i,k is lower, and discharge less energy when the discharging price is higher. This is because a higher discharging price means the grid company needs to pay more money to EV users, which will increase the grid operating cost. However, if the discharging price is too low, even below the charging price, EVs might discharge energy to charge itself while the battery degradation cost is not taken into consideration. This will definitely lead to an increase on EV user's cost, which is unacceptable for EV users since the they wish to achieve economy benefits from attending V2G activity. The advantage of strategy B (proposed in this paper) is, because the user's cost is taken into consideration while conducting system optimization, a compromised solution between system operating cost and EV users' cost will be chosen by using FODM. The system will not reduce the EVs' discharging power when the discharging price is relatively high (e.g., in the peak load period), and will limit EVs' discharging power when the discharging price is too low. As a result, CO 2 emission is reduced and the willing of EV users to attend the V2G activity is enhanced.
From Table 5 it can be seen that EVs' charging power is relatively high while EV users' cost is not taken into consideration. This is because a higher charging price means more income from EV users to the grid company, the system operating cost will also drop if the charging price is increasing. However, it is unacceptable for EV users if the charging price is too high. By using strategy B, system will drop part of the wind power and limit EV's charging power when the charging price is too high. Hence, the interest of EV users is protected. Otherwise, due to strategy B taking battery degradation cost into consideration inside its objective function, the charging power of strategy B will not be as high as that of strategy A.

Conclusions
This paper presents a dynamic environmental dispatch model which coordinates a cluster of charging and discharging controllable EV units with large-scale wind power farms and large thermal plants. A MOPSO algorithm for the optimization of system operating costs, CO 2 emissions, wind curtailment, and EV user's cost is presented to obtain the Pareto front of the optimization problem and a FODM method is proposed then to extract out the best compromise solution from the Pareto optimal set. The simulation results further demonstrate that the best compromised solution calculated by the proposed strategy is able to schedule charging and discharging of EV aggregator to balance supply and demand for active power, according to grid status. System operating costs, CO 2 emissions, and wind curtailment can be relatively reduced and the objective of EV users' cost is demonstrated to be able to ensure the interest of EV users attending V2G activity. The contribution of the proposed method can be described as follows: (1) establishing a mathematical model of wind farms, EV aggregators, and thermal plants while considering many operating constrains and security requirements; (2) using MOPSO as the solver to deal with the proposed MOP; (3) proposed a FODM-based algorithm to extract out the compromised best solution of the Pareto Front; (4) take EV users' cost into consideration in the optimization model, thus a best compromised solution between system operating cost, EV users' cost, CO 2 emission, and wind curtailment is obtained.