Research on Flexibility Evaluation Method of Distribution System Based on Renewable Energy and Electric Vehicles

High penetration of renewable energy brings economic and environmental benefits to the grid, but also causes challenges to the distribution system flexibility due to its uncertainty. As a new variable load, electric vehicles (EVs) can increase system’s flexibility through interactions with the grid and promote the consumption of renewable energy. This paper studies the flexibility evaluation of distribution system considering the interaction between EVs and the grid. Firstly, the charging and discharging control strategies of EVs have been analyzed to reduce the impact of renewable energy output fluctuations on the distribution system flexibility, such as the G2V (Grid to vehicle) and V2G (Vehicle to grid) modes. Secondly, the battery capacity, driving plan, and transportation network are modeled to study the impact of the EVs connected to the grid on the distribution system flexibility. Thirdly, the flexibility evaluation method of distribution network is proposed based on the feasibility analysis of uncertain region. Finally, the analyses of the example show that the proposed method can realize the evaluation of the time series flexibility of the distribution system, which proves the effectiveness of the method proposed in this paper.


ST
time span, such as 24 a, b, c the operation cost coefficients of the thermal power unit P g,t the output of the thermal power unit at time t λ pv the cost coefficient of PV output curtailment P pvc,t the planned PV output curtailment x t , y t , z t the 0-1 variable representing the status, start and stop of the thermal power unit P g,min , P g,max the lower and upper limits of thermal power unit output R ramp the ramp rate of thermal power unit P pr pv,t the predicted output of PV P V ,t the charging power of the grid to the electric vehicle V at time t N ev the total number of EVs P Vmin,t , P Vmax,t the minimum and maximum charging power of electric vehicle V at time t The associate editor coordinating the review of this manuscript and approving it for publication was Ramazan Bayindir . P L,t the total load demand except for electric vehicles D N the set of distribution network nodes D A the set of power lines (j, n) ∈ L all lines with j as the starting node p ij the power flow of line ij p j the load at node j P pv,t the actual output of PV P up pv,t , P low pv,t the upper and lower limits of PV output µ pv,l , µ pv,u the upper and lower limits of PV output prediction error Cap r the battery capacity of EV D EV the actual driving distance of EV D AB the distance between two places in the gridding traffic network k the standard deviation parameter  nd the normalized deviation value x the decision variable in the optimization process P i , Q i injected active power and injected reactive power at node i respectively. U i , U j voltage amplitude of node i and j G ij , B ij conductance and susceptance of branch ij θ ij voltage phase angle difference between node i and j I l the branch l current amplitude I max l the maximum value of branch l current amplitude U i he voltage amplitude of node i U imax , U imin the upper and lower limits of node i voltage amplitude P REi,t , Q REi,t the active power and reactive power of renewable energy output at node i and time t S REi the capacity of renewable energy inverter at node i cosφ the lower limit of the power factor of PV output E SoC,t the residual capacity of energy storage at time t E SoC,min the minimum SoC E SoC,max the maximum SoC u c,t the charging flag at time t p c,t the actual charging power at time t p d,t the actual discharging power at time t p cmax the maximum charging power p dmax the maximum discharging power η c the charging efficiency of the energy storage battery η d the discharging efficiency t the time interval of charging and discharging X i,j (t), X i,j (t+1) the position update vectors of particle i in dimension j at the iteration t and t + 1 r 1 (t), r 2 (t), r 3 (t) the randomly generated coefficients X k,j (t), X i,j (t) the position vectors of particle i in dimension j at the iteration t X j (t) the average positions of all particles in dimension j at the iteration t T MCS the time span of simulation P IUF,t the value of insufficient upward flexibility at time t P IDF,t the value of insufficient downward flexibility at time t

I. INTRODUCTION A. MOTIVATION
Renewable energy has the advantages of environmental friendliness and economic convenience, such as wind power and solar energy. All countries in the world are vigorously promoting the construction of renewable energy power plants. The European Union, the United States and China have set targets for renewable energy generation to account for 100%, 80% and 60% of total power generation by 2050.
Although the advantages of renewable energy generation are significant, it also has the defects of fluctuation and uncertainty. Therefore, it is necessary to evaluate the flexibility of the distribution system in [1], and prepare measures for the supply and demand side safe operation of the distribution system within the uncertainty of renewable energy generation in advance. The flexibility of distribution system is mainly reflected in dispatching and operation, which is one of the important indexes to measure the operation performance of distribution system in [2]. Sufficient flexibility is helpful to alleviate or even eliminate the negative impact of renewable energy generation in distribution system, which is of great significance to the sustainable development of energy.

B. RELATED WORKS
The concept of distributon system flexibility was first formally proposed by the International Energy Agency (IEA) in [3] and the North American Electric Reliability Corporation (NERC) in [4]. At present, the references about flexibility evaluation are mostly from the perspective of distributon system. In [5], the ramp rate of thermal power unit is used to evaluate the flexibility sufficience of the distribution system generation side. Reference [6] proposes the definition of ''flexibility interval'' to evaluate the ability distribution system to cope with renewable energy fluctuations. In [7], robust optimization method is used to calculate the maximum fluctuation range of renewable energy that the distribution system can operate safely. But References [5]- [7] did not consider the impact of distribution network transmission constrains, which may lead to the evaluation of the flexibility is too optimistic. In [8], a distribution system flexibility evaluation method considering network constraints is proposed, but the network interconnection between multiple regions is not considered. To solve this problem, Reference [9] proposes a method to improve the flexibility and the consumption of renewable power output based on distribtuon network interconnection. In [10]- [12], the flexibility evaluation of distribution system with high proportion of renewable energy and potential problems are analyzed and discussed. Aiming at the problem of severe fluctuation of load in the distribution network with high proportion of grid-connected renewable energy, a micro network optimal scheduling model based on flexibility is proposed in [10]. An index system is proposed for the flexibility of distribution network in [11], [12], which comprehensively explains the indexes related to the flexibility. Electric vehicles (EVs) has attracted more and more attention due to its energy-saving and emission-reduction characteristics. At the same time, the grid-connected EVs has become a new research hotspot. The charge and discharge control strategy of EV is modeled in detail in [13]- [15], but the uncertainty of renewable energy generation is not taken into account. In order to improve the enthusiasm of EVs to participate in charge and discharge management, a series of studies on the formulation of charge and discharge price are carried out in [16]- [21]. As a variable load, EV is not only variable in time scale, but also uncertain in space scale due to the difference of driving plan. The traffic network is modeled in [22], [23], and the optimal driving path strategy of EVs is studied based on the coupling of traffic network and distribution network. Considering the constraints of EV and traffic network, the reliability of distribution system is studied in [24], [25], but the impact of high proportion of grid-connected renewable energy generation is not considered. At present, there are few studies on the evaluation of distribution system flexibility considering the interaction between EVs and distribution network in the world. Most of them have not involved the prediction error of uncertain variables such as renewable energy and EVs and the impact of the error on the time series operation or planning. In fact, a reasonable charge and discharge control strategy can not only reduce the uncertainty threat of EV as a variable load to the grid, but also improve the flexibility of the distribtuion system and increase the consumption of renewable energy output.
The prediction error of uncertain variables can be described by interval or region. Similar concepts have been used to describe the security work in the field of distribution network. A series of studies have been carried out on the existence proof, full quadrant characteristics and boundary observation of distribution network security region in [26]- [28]. The above studies ignore the constraints of node voltage in the distribution network, and do not consider the grid-connected EVs, high proportion of grid-connected renewable energy, and the prediction error of them.

C. CONTRIBUTION
In view of the above problems, flexibility evaluation method of distribution system is proposed in this paper based on renewable energy and EVs. The list of key contributions are as follows.
(1) The battery capacity, driving plan of EV and traffic network are modeled, the concept of distribution system flexibility considering EVs connected to the grid is defined, and the charging and discharging control strategy of EV in two modes of G2V and V2G is proposed, which can effectively reduce the impact of renewable energy output fluctuation on system flexibility.
(2) Based on the coverage of the feasible region to the uncertain region, the basic idea of flexibility evaluation is described. The concept of direction matrix is extended to the time series space, and the time series direction matrix is used to describe the specific direction of the uncertain region in the operation of distribution network, and analyze the periods of insufficient flexibility and uncertain variables that cause the insufficient flexibility.
(3) Based on the time series direction matrix, the uncertain variables are described. On this basis, the mathematical model of distribution network flexibility evaluation is established. The model can be used to quantify the dispatch capability of distribution network under the disturbance of multiple uncertain variables.
(4) The analyses of the example show that the proposed method can realize the evaluation of the time series flexibility of the distribution system, and can find the main period of the insufficient flexibility according to the value of the elements in the time series direction matrix. The resource allocation of the distribution system and its time series operation status will affect the flexibility index of the distribution network. At the same time, the interactions between distribution system and electric vehicles can improve distribution system flexibility greatly. Reasonable penetration of electric vehicles can maintain considerable flexibility of the system and improve consumption of renewable energy output.

D. ORGANIZATION
In this context, the definition of distribution system flexibility considering the renewable energy and EV is firstly proposed in Part II. Secondly, the charge and discharge control strategy of EV is modeled in detail, and the method to improve the flexibility of the system is studied in Part III. Thirdly, the traffic network, battery capacity and driving plan of EV are also modeled to obtain the information needed to implement the charge and discharge control strategy of EV in Part IV. Fourthly, for the scenario of high proportion of grid-connected renewable energy, considering the constraints such as node voltage in the time series operation of distribution network, a time series flexibility evaluation method of distribution network is proposed based on the adaptability of distribution network to uncertain variables in Part V. Finally, the effectiveness of the proposed method is verified by an example in Part VI. VOLUME 8, 2020

II. DEFINITION OF DISTRIBTUTION SYSTEM FLEXIBILITY
Based on the day ahead predictive value of PV output and the driving / charging plan reported by the EV, the day ahead unit commitment dispatching is carried out. The objective function is shown in Eq.1 and the constraints are shown in Eqs.2-8.
x t P g,min ≤ P g,t ≤ x t P g,max n|(j,n)∈L In formula, the time span ST is 24 hours. a, b and c are the operation cost coefficients of the thermal power unit. P g,t is the output of the thermal power unit at time t. λ pv is the cost coefficient of PV output curtailment. P pvc,t is the planned PV output curtailment. x t , y t , and z t are the 0-1 variable representing the status, start and stop of the thermal power unit. P g,min and P g,max represent the lower and upper limits of thermal power unit output. R ramp is the ramp rate of thermal power unit. P pr pv,t is the predicted output of PV. P V ,t is the charging power of the grid to the electric vehicle V at time t. N ev is the total number of EVs. P Vmin,t and P Vmax,t represent the minimum and maximum charging power of electric vehicle V at time t. P L,t is the total load demand except for electric vehicles. D N represents the set of distribution network nodes. D A is the set of power lines. (j, n) ∈ L represents all lines with j as the starting node. p ij is the power flow of line ij. p j is the load at node j.
Eqs.2-3 show the upper and lower limits of output and ramp rate constraints of thermal power unit. Eqs.4-5 show the start and stop state constraints of thermal power unit. Eq.6 is the power balance constraint of distribution system. Eq.7 limits the charging power of grid-connected EV. Eq.8 is the node power balance constraint.
By solving the unit commitment problem, the start and stop status, output of thermal power units and the charging plan of grid-connected EVs at all times can be obtained. Due to the randomness and fluctuation of PV output, there is a certain deviation between the actual output and the predicted value. In this paper, the uncertainty interval is used to characterize the accuracy of PV output prediction, as shown in Eq.9.
In formula, P pv,t is the actual output of PV, which can be obtained by Monte Carlo sampling method based on P pr pv,t and Eq.9. P up pv,t and P low pv,t are the upper and lower limits of PV output respectively. µ pv,l and µ pv,u are the upper and lower limits of PV output prediction error respectively.
In this paper, flexibility is defined as ''the ability of the system to use its own supply and demand side resource to cope with the uncertainty of renewable energy and EVs''. Flexibility can be divided into two categories: upward flexibility and downward flexibility. Upward flexibility refers to the power output that the system can increase when PV output suddenly decreases. Downward flexibility refers to the power output that the system can decrease when PV output suddenly increases. The demand for flexible resources of the system at time t is defined as P G,t , and the calculation method is shown as follows.
(1) When P pv,t > P pr pv,t , the system needs a downward flexibility resource P G,t = P pv,t − P pr pv,t to cope with the sudden increase of PV output.
(2) When 0 ≤ P pr pv,t − P pv,t ≤ P pvc,t , the reduction of PV output will not produce flexibility demand, and the actual PV output curtailment will change from the planned value P pvc,t to P pvc,t − P pr pv,t + P pv,t .
(3) When P pvc,t < P pr pv,t − P pv,t , an upward flexibility resource P G,t = P pr pv,t − P pvc,t − P pv,t is needed to cope with the sudden reduction of PV output.

III. CHARGE AND DISCHARGE CONTROL STRATEGY OF EV
As the link between the distribution system and the transportation system, EVs can provide flexible resources to the distribution network through reasonable charge and discharge control strategy, and improve the consumption of renewable energy output. The charging and discharging strategy of EV includes two modes, such as grid-to-vehicle (G2V) and vehicle-to-grid (V2G), which are shown as follows.

A. CONTROL STRATEGY OF G2V
G2V control mode refers to the distribtuion system changes the total load of demand side by adjusting the charging power of grid-connected EV at each time to balance the supply and demand. Define the state of charge (SoC) of EV i as SoC i shown in Eq.10. E i and E N ,i are the residual capacity and rated capacity of electric vehicles respectively.

1) DOWNWARD FLEXIBILITY PROVIDED BY G2V
When P pv,t > P pr pv,t , the thermal power unit can provide downward flexibility by reducing output. If P G,t ≤ min{R ramp , P g,t − P g,min }, then the ramp rate of thermal power units can meet the downward flexibility demand of the distribution system without the control strategy of G2V. If P G,t > min{R ramp , P g,t − P g,min }, then it is necessary to consume the excess PV output through the control strategy of G2V, and P G,t = P G,t − min{R ramp , P g,t − P g,min }. The control strategy of G2V is shown as follows.
(1) All the grid-connected EVs are arranged in ascending order according to SoC at time t, and then a set G2Vmax is formed, N G2Vmax is the number of elements in this set.
the optimization model is solved in Eq.11. After n G2V is solved, the first n G2V elements in set G2Vmax are taken to form a new set G2V , that is, the set of electric vehicles participating in the G2V control. The distribution system will supply power to these electric vehicles according to the maximum charging power P Vmax,t , which can fully meet the downward flexibility demand. ( then all electric vehicles in set G2Vmax will participate in the G2V control, and the distribution system will supply power to these electric vehicles according to the maximum charging power P Vmax,t . However, there will still be insufficient downward flexibility of the distribution system, and the lack of flexibility is When P pr pv,t − P pv,t > P pvc,t , the thermal power unit can provide upward flexibility by increasing output. If P G,t ≤ min{R ramp , P g,max − P g,t }, then the ramp rate of thermal power units can meet the upward flexibility demand of the distribution system without the control strategy of G2V. If P G,t > min{R ramp , P g,max − P g,t }, then it is necessary to meet the load demand through the control strategy of G2V, and P G,t = P G,t −min{R ramp , P g,max −P g,t }. The control strategy of G2V is shown as follows. (1) All the grid-connected EVs are arranged in descending order according to SoC at time t, and then a set ' G2Vmax is formed, N ' G2Vmax is the number of elements in this set. ( the optimization model is solved in Eq.12. After n' G2V is solved, the first n' G2V elements in set ' G2Vmax are taken to form a new set ' G2V , that is, the set of EVs participating in the G2V control. The distribution system will supply power to these EVs according to the minimum charging power P Vmin,t , which can fully meet the upward flexibility demand.
then all EVs in set ' G2Vmax will participate in the G2V control, and the distribution system will supply power to these EVs according to the maximum charging power P Vmin,t . However, there will still be insufficient upward flexibility of the distribution system, and the lack of flexibility is P G,t − ( V2G control means that the grid-connected EVs can transmit the power stored in its battery to the grid, together with the generator units to meet the power demand of the load side. Section A of Part III has introduced the G2V control strategy of EVs in detail. When all the grid-connected EVs participate in the G2V control strategy, and the upward flexibility of the distribution system is still insufficient, the V2G control strategy needs to be enabled, which is shown as follows. (1) All the grid-connected electric vehicles at time t are arranged in descending order according to SoC, and then a set V 2Gmax is formed, N V 2Gmax is the number of elements in this set.
(2) If The optimization model is solved in Eq.13 based on . P Vmaxf ,t is the maximum dis-charging power. After n V 2G is solved, the first n V 2G elements in set V 2Gmax are taken to form a new set V 2G , that is, the set of EVs participating in the V2G control. The EV will supply power to distribution system according to the maximum discharging power P Vmaxf ,t , which can fully meet the upward flexibility demand. (

then all EVs in set
V 2Gmax will participate in the V2G control, and the EV will supply power to distribution system according to the maximum discharging power P Vmaxf ,t . However, there will still be insufficient upward flexibility of the distribution system, and the lack of flexibility is P G,t −

IV. CHARACTERISTICS MODEL OF ELECTRIC VEHICLE
The interaction between EV and distribtuion network will improve the flexibility of distribution system, and the implementation of charge and discharge control strategy needs to obtain the information of the battery capacity characteristics of EVs. Therefore, this chapter models the battery capacity, driving plan of EVs and traffic network.

A. BATTERY CAPACITY MODEL
According to the research results on the European market of EV in [16], four types of EVs are summarized as follows. Type L7e is used for carrying goods, which maximum empty mass is 400kg or 550kg. Type M1 is used for carrying people, which is up to 9 passengers. Type N1 is used for carrying goods, which maximum empty mass is 3500kg. Type N2 is used for carrying goods, which maximum empty mass is between 3500kg and 12000kg. The probability distributions of battery capacity of these four types are shown in TABLE 1.

VOLUME 8, 2020
The probability expressions of battery capacity according to gamma distribution and normal distribution are shown in Eqs.14-15.
In formula, Cap r represents the battery capacity of EV. In order to simplify the calculation, this paper assumes that the single maximum driving distance of EV has a linear relationship with battery capacity.

B. DRIVING PLAN OF ELECTRIC VEHICLE
EVs are usually connected to the grid only at home, at work or at charging stations. In order to know how many times and how long the EV is connected to the distribution network every day, it is necessary to obtain its daily driving plan. Define the set of possible places that EVs may reach every day as dest = {home, working place, hospital, shopping mall, charging station, scenic spot}. Through investigation and statistics, the probability of EVs reaching these places can be obtained. The decision method of EV daily driving plan is shown as follows.
(1) After EV leaving home in the morning, it is assumed that the first destination of EV is the working place in this paper.
(2) It is assumed that the last destination of EV is home.
(3) According to the reaching probability of each palce, the daily driving plan of EV can be obtained by random sampling.

C. SIMULATION OF TRAFFIC NETWORK
In this paper, the method proposed in [24] is used to simulate the traffic network. The method is briefly described as follows.
(1) It is assumed that the traffic roads are parallel or vertical, as shown in Fig.1. The distance between any two points can be expressed by the total length of the lines connecting the two points.
(2) In fact, the traffic roads may not be parallel or vertical, and the EVs may not drive in the shortest distance. Therefore, the actual distance is calculated by normal distribution, as shown in Eq.16.
In formula, D EV is the actual driving distance of EV. D AB is the distance between two places in the gridding traffic network. k is the standard deviation parameter.
(3) When the residual battery capacity of EV is lower than the lower limit in TABLE 1, i.e. SoC< δ FC . δ FC is the ratio of lower limit to rate capacity of battery. EV will change the driving plan and drive to the nearest charging station.

V. FLEXIBILITY EVALUATION OF DISTRIBUTION NETWORK A. THEORETICAL ANALYSIS OF FLEXIBILITY EVALUATION
Because the adaptability of distribution network to uncertain variables has not been considered in the existing research of distribution network flexibility evaluation, it is necessary to propose the quantification indexes of flexibility evaluation. In order to describe any point in the uncertain region, an uncertain variable model based on time serires direction matrix is established.

1) DEFINITION AND QUANTIFICATION OF INDEXES
Considering the scenario of high proportion grid-connected renewable energy, the flexibility evaluation of distribution system is quantitative analysis of the feasible operation region of distribution network based on uncertain variables, such as renewable energy output forecast error, load forecast error, et al.
For the sake of simplicity, the two-dimensional uncertain variables are taken as an example to illustrate. The feasible region and the uncertain region can be expressed in a plane way. As shown in Fig.2, all points in the coordinate system correspond to a combination of uncertain variables, which is called the operation point. The red dotted curve is the boundary of the feasible region, and the internal area is the feasible region. The black rectangle of dash-dotted line represents the boundary of the uncertain region, and the internal area is called the uncertain region. For the combination of uncertain variables in the feasible region, the distribution system can dispatch the resources of supply and demand side to get the feasible operation scheme that does not violate the constraints. For the combination of uncertain variables outside the feasible region, the distribution system cannot get the feasible operation scheme that meets the constraints.
Point O in Fig.2 is the expected point, which corresponds to the expected value of uncertain variables combination. For any direction D starting from point O in the uncertain region, the feasible boundary point A and the uncertain boundary point B in this direction can be obtained. For any point C in this direction, the corresponding normalized deviation value is recorded as value nd , and the calculation method is shown in Eq.17.
In formula, || · || 2 is 2-norm. For this direction, the maximum feasible deviation is recorded as fd max , which is that the distance between the expected point and the boundary point of the feasible region is divided by the distance between the expected point and the boundary point of the uncertain region, i.e. the feasible maximum value nd based on the constraints. The calculation method is shown as follow.
For all directions in the uncertain region, the minimum δ D exists in only one direction, which is called the feasible maximum deviation as the flexibility index. The boundary point of the feasible region is the critical point, such as A c in Fig.2. The corresponding direction is the critical direction, and the boundary point of the uncertain region in the critical direction is the B c in Fig.2. When the uncertain variable is extended from two dimension to n dimension, the uncertain region is correspondingly extended to n dimension super rectangle in [29].

2) UNCERTAIN VARIABLE MODEL BASED ON TIME SERIES DIRECTION MATRIX
In fact, the operating point in the uncertain region can reflect the deviation degree of the uncertain variable from its expected value, which can be represented by two variables, such as the direction vector starting from the expected point and the normalized deviation value. The direction vector can be described by the direction matrix in [30]. The direction matrix is a diagonal matrix, each element on the diagonal represents an uncertain variable, and the value range of the element is [−1, 1]. However, it is possible for multiple direction matrices to represent the same direction, such as diag(0.5, 1) and diag(0.3, 0.6). To ensure that the direction matrix in the same direction is consistent, at least one element in the direction matrix is specified as 1 or -1. Specifically, if the direction matrix is divided by the largest absolute value of each element, then diag(0.3, 0.6) becomes diag(0.5, 1).
Because an uncertain variable in the direction matrix is only described by a single interval, it can only be used in scenarios where the fluctuation range of uncertain variable is small in different periods. However, for the time series operation of distribution system, the interval of uncertain variables in different periods is quite different. For example, because the output characteristics of PV are mainly affected by the light intensity, its output is zero at night every day, which is no uncertainty naturally. It is no need to consider the output of PV in the uncertain variables at night. However, the output of PV at noon generally deviates in a certain range of the predicted value, and It is necessary to consider the output of PV in the uncertain variables at noon. Therefore, the same type of uncertain variables also has a certain correlation in time series, such as photovoltaic output power. According to the above characteristics of uncertain variables, the concept of direction matrix is extended to time series direction matrix, such as in Eq.19.
In formula, T is the total number of periods. N is the total number of uncertain variables. DM is the time series direction matrix of T×N dimension, which represents a certain direction in the uncertain region. dm n t is the element value corresponding to the uncertain variable n at time t in the time series direction matrix, t=1, 2, . . . , T, n=1, 2, . . . , N. Since the boundary point of the uncertainty region in the critical direction is located at least one-dimensional boundary of the uncertainty region, such as point N b in Fig.2, one element in the time series direction matrix is 1 or −1 at least.
Any operating point in the uncertain region can be expressed as follow.
In formula, var is the uncertain variable. var exp is the expected value of the uncertain variable. var is the difference between the upper/lower limit of the uncertainty variable deviation and the expected value. ''•'' is the Hadamard product, which is the multiplication of the corresponding elements in the same position of matrixs.

B. FLEXIBILITY EVALUATION MODEL OF DISTRIBUTION NETWORK
Due to the flexibility evaluation needs to calculate the feasible maximum deviation and find the critical direction, the master problem and sub problem are set up in the research framework of the flexibility evaluation, and the mathematical model is established. VOLUME 8, 2020

1) RESEARCH FRAMEWORK OF FLEXIBILITY EVALUATION
According to the flexibility quantification indexe proposed in Section A of Part V, the flexibility index is essentially the feasible maximum deviation in the critical direction, so the flexibility evaluation needs to solve two problems as follows. One is to realize the calculation of the feasible maximum deviation in any direction in the uncertainty region, the other is to find the smallest feasible maximum deviation of all directions, so as to find the critical direction.
Based on the above analysis, this paper proposes a time series flexibility evaluation model including the master problem and the sub problem. The master problem seeks the critical direction from all directions in the uncertainty region, and the decision variable is the time series direction matrix. For each time series direction matrix of the master problem, the sub problem will be called to calculate the feasible maximum deviation of this direction. The objective of the sub problem is to solve the boundary of feasible region in a direction, i.e. the maximum deviation. The decision variables include the reactive power output of PV inverter and the charging and discharging power of energy storage. When the sub problem converges, the boundary of feasible region in this direction is returned to the master problem. Finally, the master problem can obtain the time series direction matrix of critical direction and the feasible maximum deviation of critical direction over many iterations, and the feasible maximum deviation of critical direction is the flexibility index.

2) MODEL OF MASTER PROBLEM a: OBJECTIVE FUNCTION
Eq.21 is the objective function of the master problem of flexibility evaluation. The critical direction is obtained in the uncertainty region, and the feasible maximum deviation of the critical direction is taken as the flexibility index.
In formula, F is the flexibility evaluation index. When F > 1, the distribution network has sufficient flexibility. On the basis of satisfying all the operation conditions in the uncertain region, the distribution system can also cope with a larger range of uncertain variables. When F = 1, the distribution network has precise flexibility, and the distribution system just satisfies all the operation conditions in the uncertain region. When F < 1, the distribution system has insufficient flexibility, and the node voltage or branch power will over-limit in some directions of the uncertain region.
In formula, dm n t is the element value corresponding to the uncertain variable n in the time series direction matrix DM at time t.

3) MODEL OF SUB PROBLEM a: OBJECTIVE FUNCTION
The objective function of the sub problem of flexibility evaluation is shown in Eq.23, which indicates that the feasible maximum deviation can be obtained in a direction of the time series direction matrix D of the uncertainty region.
In formula, x is the decision variable in the optimization process.

b: CONSTRAINTS OF DISTRIBUTION NETWORK
(1) power flow constraint In formula, P i and Q i are injected active power and injected reactive power at node i respectively. U i and U j are voltage amplitude of node i and j. G ij and B ij are conductance and susceptance of branch ij respectively. θ ij is voltage phase angle difference between node i and j.
(2) branch current constraint In formula, I l is the branch l current amplitude. I max l is the maximum value of branch l current amplitude.
(3) node voltage constraint In formula, U i is the voltage amplitude of node i. U imax and U imin are the upper and lower limits of node i voltage amplitude respectively.
(4) renewable energy generation constraints Reactive power is generated or absorbed by controlling the renewable energy inverter, which can effectively adjust the node voltage of distribution network. The constraints of PV include the capacity constraint of renewable energy inverter and the power factor constraint of renewable energy.
In formula, P REi,t and Q REi,t are the active power output and reactive power output of renewable energy at node i and time t. S REi is the capacity of renewable energy inverter at node i. cosφ is the lower limit of the power factor of renewable energy output. The reactive capacity of the renewable energy inverter is limited by the capacity and power factor of the inverter, and the renewable energy reactive capacity in the time series operation is determined by above two kinds of constraints.
(5) energy storage battery constraints The constraints of the energy storage battery include the upper and lower limit constraints of the state of charge (SoC), the power constraints of charging and discharging, and the recurrence relation constraints of SoC, as shown below.
Upper and lower limit constraints of SoC is shown as follow. In the process of charge and discharge of energy storage battery, overcharging or overdischarging will have a negative impact. Therefore, it is necessary to set the upper and lower limit of SoC of the energy storage battery.
In formula, E SoC,t represents the residual capacity of energy storage at time t. E SoC,min represents the minimum SoC. E SoC,max represents the maximum SoC.
Power constraints of charging and discharging is shown as follow.
In formula, u c,t is the charging flag at time t. When the energy storage battery is charging, u c,t = 1. When it is not charging, u c,t = 0. u d,t is the discharging flag at time t. When the energy storage battery is discharging, u d,t = 1. When it is not discharging, u d,t = 0.
In formula, p c,t is the actual charging power at time t. p d,t is the actual discharging power at time t. p cmax is the maximum charging power. p dmax is the maximum discharging power.
Recurrence relation constraints of SoC is shown as follow.
In formula, η c represents the charging efficiency of the energy storage battery. η d represents the discharging efficiency. t represents the time interval of charging and discharging.

C. SOLUTION ALGORITHM OF MODEL
According to the characteristics of the master problem and the sub problem of the time series flexibility evaluation of the distribution system, a solution algorithm of the dichotomy embedded social learning particle swarm optimization is proposed. The dichotomy is used to find the feasible maximum deviation of a direction in the uncertain region, and the social learning particle swarm optimization algorithm (SLPSO) is used to find the critical direction.

1) SOLUTION ALGORITHM OF SUB PROBLEM
The purpose of flexibility evaluation sub problem is to get the maximum feasible deviation in a direction, i.e. the boundary of feasible region in a direction. In this paper, the dichotomy is used to get the maximum feasible deviation in a direction, as shown below.
(1) In initialization, the expected point of the uncertainty region value nd,0 and search the upper bound of the deviation value nd,up,0 are set.
(2) Let t = 1, the two endpoints of the t-th iteration are value nd,t and value nd,up,t respectively, and the midpoint value nd,mid,t can be calculated at the same time.
(3) Judge the feasibility of midpoint value nd,mid,t can be judged as follows. If value nd,mid,t is feasible, the critical point in this direction is between value nd,mid,t and value nd,up,t , let value nd,t+1 = value nd,mid,t and value nd,up,t = value nd,up,t+1 . If value nd,mid,t is not feasible, the critical point in this direction is between value nd,t and value nd,mid,t , let value nd,t+1 = value nd,t and value nd,up,t+1 = value nd,mid,t .
In order to judge the feasibility, it is necessary to convert the normalized deviation value into the form of uncertain variables in Eq.20, and then the optimal power flow is calculated to determine whether there is a feasible solution.
(4) Iterations are completed until the convergence condition in Eq.32 is satisfied.
In formula, ε is the convergence factor to ensure the precise calculation of the boundary of the feasible region.

2) SOLUTION ALGORITHM OF MASTER PROBLEM
In this paper, SLPSO is used to solve the master problem. Traditional particle swarm optimization (PSO) only updates the particle position based on the historical information of population optimization and individual optimization, which leads to the loss of population diversity. PSO is easy to converge to local optimization, and is very sensitive to the initial value. The single dimension of each particle in SLPSO comes from different demonstration examples, which greatly enhances the population diversity and effectively avoids premature convergence of the algorithm. The solution flowchart of the master problem is detailed as follows.
(1) According to the search dimension, the population size pop and social influence factor sif are set. The direction matrix is expanded into one dimension, which is used as the position vector of particle swarm to initialize particle swarm.
(2) For all individuals in the population, the maximum feasible deviation is solved by dichotomy as the fitness value of all individuals, and all particles are ordered according to the fitness value.
(3) The positions of all particles are updated. The update vector of particle position consists of three parts. The first part is similar to the velocity component of PSO. The second part replaces the pbest of PSO, and the particle will learn from other particles whose fitness is better than itself. The third part replaces the gbest of PSO, and the particle will be affected by the average position of all particles in the population. The update expression of particle position is shown as follow.
In formula, X i,j (t) and X i,j (t+1) are the position update vectors of particle i in dimension j at the iteration t and t + 1 respectively. r 1 (t), r 2 (t) and r 3 (t) are the randomly generated coefficients respectively. X k,j (t) and X i,j (t) are the position vectors of particle i in dimension j at the iteration t.X j (t) is the average positions of all particles in dimension j at the iteration t.
In the solution process, the second part in Eq.33 can enhance the diversity of particles in the population and reduce the sensitivity of the algorithm to the initial population. The third part in Eq.33 can avoid the problem that When gbest is local optimum, the whole population will fall into local optimum.
According to the network structure, line and load parameters provided by the electric power company in practical application, the mathematical model of flexibility evaluation is established and solved by the proposed algorithm to obtain the the evaluation index of flexibility. Based on the evaluation index of flexibility, the adaptability degree of distribution network to uncertain variables can be evaluated. If the flexibility is insufficient, the constraints restricting the flexibility can be obtained according to the electric parameters at the critical point, such as node voltage. The flexibility can be improved by adjusting the operation strategy, adding energy storage and reactive power compensation devices, and other ways.

D. FLOWCHART OF FLEXIBILITY EVALUATION
The flowchart of flexibility evaluation is shown in Fig.4, as shown below.
(1) Parameters of distribution system and traffic system are initialized, such as PV predicted output, thermal power unit characteristics, total number of electric vehicles, battery capacity of electric vehicle, lower limit threshold δ FC of SoC, etc. The time interval is set to 1 hour. Set t=1, d=1, and the convergence threshold is δ CoV .
(2) Monte Carlo simulation method is used to make a one-year driving plan for each electric vehicle.
(3) The unit commitment models in Eqs.1-8 are solved, and the day ahead starting and stopping status and output plan of thermal power unit are obtained.
(4) Monte Carlo simulation method is used to extract the actual output of PV, and calculate the demand of distribution system for flexible resources P G,t is calculated.
(5) G2V and V2G control strategies is used to meet the flexible resources demand of distribution system.
(6) Each flexibility index is updated according to the method proposed in Sections A ∼C of Part IV.
(7) The SoC of the electric vehicle i is updated. (8) If the SoC of the electric vehicle i is lower than the predetermined threshold, the electric vehicle will stop the current driving plan and go to the nearest charging station for charging. The daily driving plan will be changed accordingly.
(9) If i < N ev , let i = i+1 and go back to (7), otherwise go to (10). (10) The end time of the day t day is set to 24. If t < t day , let t = t+1 and go back to step 4, otherwise go to step 11.
(11) The end time of the year d year is set to 365. If d < d year , let d = d+1, t=1 and go back to step 3, otherwise go to step 12.
(12) The flexibility evaluation index is calculated according to the method proposed in Sections A∼C of Part V, which is used to check whether the flexibility resources demand of distribution system can be met. If not, go back to (2), otherwise iterations are completed and results are output.

VI. EXAMPLE ANALYSIS A. DESCRIPTION OF EXAMPLE
IEEE 33 node radial distribution network is adopted in this example, relevant data of distribution network and traffic network are taken from [24]. The topological structure of distribution network is shown in Fig.5. The tie line is shown as the dotted line in Fig.5, which is normally open in this example. The reference voltage is 12.66kV. In this paper, the PV output power and load are set as uncertain variables. The prediction error of PV is 10% in [24], i.e. µ w,l =95% and µ w,u = 105%, as shown in Fig.6(a). The prediction error of load is 15% of the predicted value in [30], as shown in Fig.6(b). The PV output is much higher than load at noon. The distribution network has high proportion of grid-connected PV in this example. The nodes of PV connected to the grid are 6, 7, 13, 18, 28 and 33, each installation capacity is 1.2MW, and the minimum value of power factor is 0.95. Each PV rating is 1.2MW. There are 10000 vehicles in this system, and the penetration rate of EVs is 30%, i.e. 3000 EVs in this distirbution system. There are 750 EVs of L7e type, M1 type, N1 type and N2 type. The interval of SoC is [0.2, 0.9]. The battery capacity is obtained by random sampling according to TABLE 1, and the total EV rating is 110MWh.
Due to the strong correlation of renewable energy output in the distribution network in the same area in [31], for the sake of simplification, it is assumed that the correlation factors of all PVs are 1.0, and the correlation factors of all loads are 1.0. In this way, the uncertainty region and feasible region of a single time section can be expressed as a two-dimensional graph.
Contents discussed in the following sections are shown below. Firstly, through the flexibility evaluation of section, the influence of the distribution network time series operation state on the distribution network flexibility is illustrated. VOLUME 8, 2020 Secondly, through the result analysis of the time series flexibility evaluation, the function of the time series direction matrix and the feasibility of the proposed evaluation method are explained. Thirdly, through the calculation of the distribution network flexibility under different flexible resource allocation conditions, the effectiveness of the proposed evaluation model is shown. Fourthly, the comparison of different solution algorithms is made. Fifthly, the influence of gridconnected electric vehicles on system flexibility evaluation is detailed.

B. FLEXIBILITY EVALUATION OF SECTION
The reactive power of PV inverter is only determined by the active power output in this period. The capacity of storage energy batter is limited by the recurrence relation constraints of SoC, upper and lower limit constraints of SoC. The charge and discharge power of storage energy batter is affected by the charge and discharge power in the previous period. Therefore, the SoC of energy storage battery needs to be specified in the evaluation of section flexibility. 12:00 is taken as an example to analyze the influence of different SoC of energy storage battery on the flexibility evaluation results of section. Energy storage battery works in charging mode at noon, which is conductive to reduce the node voltage. When the SoC of energy storage battery is in the interval of [0.2, 0.6], the energy storage battery can be charged with the maximum power. Therefore, 0.6, 0.75 and 0.9 of SoC are considered for discussion, and the results are shown in Fig.7.  When SoC is 0.6 at 12:00, the energy storage battery in this period can be charged according to the maximum charging power, and the flexibility of the distribution system in this period can be evaluated. The direction matrix is diag(-0.2985, 1), which shows that PV has more influence on the flexibility evaluation index. The flexibility evaluation index is 0.8391, and distribution system has insufficient flexibility in this period. In the two-dimensional space, all directions in the two-dimensional uncertainty region are traversed based on the step length of 1 • , and the feasible region is shown in Fig.7(a). However, the feasible region still cannot completely cover the uncertainty region, and the flexibility evaluation index is less than 1, which leads to the insufficient flexibility of distribution system. When SoC is 0.75 at 12:00, the energy storage battery in this period is subject to the upper limit of SoC, and cannot be charged according to the maximum charging power. According to the flexibility evaluation of the distribution system in this period, the direction matrix is diag(-0.2985, 1), and the flexibility evaluation index is 0.6855. In this case, the feasible region and the uncertain region are shown in Fig.7(b). The feasible region cannot completely cover the uncertainty region, the flexibility evaluation index is less than 1, and this index is lower than that of SoC=0.6.

3) FLEXIBILITY EVALUATION RESULTS WHEN SOC IS 0.9 AT 12:00
When SoC is 0.9 at 12:00, the energy storage battery in this period is subject to the upper and lower limits of SoC, and cannot be charged. According to the flexibility evaluation of the distribution system in this period, the direction matrix is diag(-0.2985, 1), and the flexibility evaluation index is 0.5214. In this case, the feasible region and the uncertain region are shown in Fig.7(c). The feasible region is further reduced, the flexibility evaluation index is less than 1, and this index is lower than that of SoC = 0.6 and SoC = 0.75.
From the above analysis, two conclusions can be obtained as follows.
(1) Through the supply of flexibility resources, the value of flexibility evaluation index can be changed, but the direction matrix of time section cannot be changed. This is because only two uncertain variables in this period are involved in the flexibility evaluation of section. Increasing flexibility resources is equivalent to moving the feasible boundary towards the critical direction, and the direction matrix has not changed.
(2) When the flexibility evaluation of section is made, the current time series state of the flexible resources in distribution system will have a great impact on the flexibility evaluation index. It is necessary to formulate reasonable operation strategies to make the state of the flexible resources best, so as to realize the release of the resource flexibility potential in distribution system.

C. FLEXIBILITY EVALUATION OF TIME SERIES 1) RESULTS
The uncertain variables include load and PV power output, as shown in Fig.6. The simulation period T is set to 24, and the initial SoC of energy storage battery is set to 0.5. Through the simulation calculation, the flexibility evaluation index is 0.9322, which shows that the flexibility resources can not fully meet the flexibility demand of distribution system, i.e. distribution system has an insufficient flexibility, the feasible region and the uncertain region are shown in Fig.8. The time series direction matrix corresponding to the critical point is 24 × 2 dimension, as shown in TABLE 2. The time series direction matrix has been bold and reserved four decimal places.  Since PV has no power output in the period of 1:00-7:00 and 20:00-24:00, the element is 0 in the corresponding position of the direction matrix. At the same time, the elements in matrix corresponding to the load are all positive in above periods. Therefore, the distribution system is more likely to have node overload, node voltage and line transmission capacity over-limit in above periods.
As the same as the direction matrix in the flexibility evaluation result of single time section, the larger the element value in the direction matrix in the time series flexibility evaluation result, the greater the impact on the flexibility evaluation index. The element value of the direction matrix is large in the period of 11:00-15:00, especially the elements corresponding to PV, which shows that the elements of PV have a great impact on the results of time series flexibility evaluation. The main reason for this situation is that the higher VOLUME 8, 2020 power output of PV in above period, the more serious problem of node voltage over-limit will occur. The position of ''1'' element in the direction matrix corresponds to the PV at 13:00, and the problem of node voltage over-limit in this period is more serious than that of other periods.

2) CURVE OF NODE VOLTAGE AND VALUE OF DECISION VARIABLE AT CRITICAL POINT
Based on a set of uncertain variable values corresponding to the critical point, the distribution system is optimized with the maximum consumption of PV power output as the objective function. The node voltage of PVs is shown in Fig.9(a), reactive power output of distributed photovoltaic inverter is shown in Fig.9(b), and the SoC of energy storage battergy is shown in Fig.9(c).
As shown in Fig.9(a), the over-upper limit of node voltage occurred from 11:00 to 15:00, and the element value corresponding to PV is relatively large in above period. The over-lower limit of node voltage occurred at 8:00 and from 19:00 to 20:00, and the element value corresponding to load is relatively large in above period.
As shown in Fig.9(b), the reactive power output of the PV inverter is mainly concentrated from 11:00 to 15:00. The inverter absorbs the reactive power to reduce the node voltage and the degree of insufficient flexibility.
As shown in Fig.9(c), the energy storage battery is discharged at 6:00 to raise the node voltage of the distribution system. At the same time, the SoC is minimized to ensure that the energy storage battery can be charged to the upper limit capacity at noon to reduce the node voltage and improve the flexibility.

D. ANALYSIS OF FLEXIBLE RESOURCE CONFIGURATION
The impact of different flexible resource configuration on the flexibility evaluation index is analyzed to explain the effectiveness of the flexibility evaluation model. Four flexible resource configuration schemes are proposed to analyze the relationship between flexibility evaluation index and flexibility resource, as shown below. Scheme 1, no flexible resources in the distribution system. Scheme 2, PV is considered only. Scheme 3, both PV and energy storage battery are considered, the capacity of energy storage battery is 80kWh, the maximum charging power is 24kW, and the maximum discharging power is 38kW. Scheme 4, the distributed photovoltaic inverter and energy storage are considered, the capacity of energy storage is 140kWh, the maximum charging power is 42kW, and the maximum discharging power is 68kW.
The flexibility evaluation index of each scheme is shown in TABLE 3. The flexibility evaluation index is 0 in scheme 1. Without considering any flexible resources, the power flow calculation result has already violated the constraints based on the expected value of PV and load. The distribution system can not cope with more uncertain variables, so the result is 0. Considering the reactive power output of the PV inverter, the flexibility is greatly improved, and flexibility evaluation index is increased to 0.5207 in scheme 2. Considering the impact of energy storage battery charging and discharging on the flexibility, the flexibility evaluation index is increased to 0.9322 in scheme 3, but the distribution system is still of insufficient flexibility. Considering the increase of the energy storage battery capacity, the flexibility evaluation index is increase to 1.0392 in scheme 4. Therefore, the flexibility evaluation index may be greater than 1 in this research based on Eq. 21. Scheme 4 the flexibility evaluation index is greater than 1 (1.0392), which indicates the value of feasible boundary point A is larger than the value of uncertain boundary point B in each direction in Fig.2 That is, the uncertainty region is within the feasible region. The distribution network has sufficient flexibility. On the basis of satisfying all the operation conditions in the uncertain region, the distribution system can also cope with a larger range of uncertain variables.
Through the comparison of scheme 1, scheme 2 and scheme 3, it can be seen that the increase of flexible resource types can improve the flexibility. Through the comparison of scheme 3 and scheme 4, the increase of flexible resource capacity can also improve the flexibility.

E. COMPARISON OF DIFFERENT SOLUTION ALGORITHMS
The performance of SLPSO is detailed in this section. To ensure the fairness of algorithm comparison, the population size of SLPSO or PSO is set to 100. The learning factor of PSO is c 1 = c 2 = 2. PSO and SLPSO are respectively uesd to calculate the distribution system flexibility index of scheme 3, and the results are shown in TABLE 4. The calculation result of SLPSO algorithm is better than that of PSO. This is due to the single dimension of each particle in SLPSO comes from different demonstration examples, which greatly enhances the population diversity and effectively avoids premature convergence of the algorithm.
In formula, T MCS is the time span of simulation. P IUF,t is the value of insufficient upward flexibility at time t. P IDF,t indicates the value of insufficient downward flexibility at time t.
Calculation result of flexibility parameters is shown in TABLE 5, the IUFE of the distribution system is 0. When PV power output suddenly decreases, the system can guarantee the balance of supply and demand based on the dispatching of power plant ramp rate and electric vehicle charging and discharging control. Compared with the sufficient upward flexibility resources, the value of insufficient downward flexibility is 81.8257MWh. The total energy used by EVs for G2V and V2G control each year is 362.81MWh. The charge and discharge control strategy will enhance the flexibility of the distribution system. In addition, the energy of EVs participating in G2V control is greater than that of participating in V2G control. This is due to G2V mode can not only meet the demand for upward flexibility, but also satisfy the demand for downward flexibility. Moreover, G2V mode has a higher priority than V2G mode in response to the demand for upward flexibility of the system.
The increase of grid-connected EVs means that more EVs can participate in the control of G2V and V2G, so as to increase consumption of PV output and provide more flexible resources for the distribtuion system. The influence of penetration rate of EVs on flexibility of distribution system will be discussed as follow.
The curve of IUFE and IDFE with different penetration rates of EVs is shown in Fig.10(a). The comparation of the energy of EVs participating in the G2V and V2G mode is shown in Fig.10(b).
As shown in Fig.10(a), the value of IDFE is a decreasing function of EV penetration. When the penetration is 30-40%, the minimum value of IUFE is 0. When the penetration of EVs is 0, i.e. there is no EVs in the system, it is noted that the value of IUFE is less than that of IDFE. This is mainly because that the PV power output is the largest at noon. In order to provide more capacity for the consumption of PV output, the thermal power unit operates at the output close to the lower limit in the corresponding period, which results in that the downward ramp ability of the thermal power unit is weaker than the upward ramp ability, and the downward flexibility is more likely to be insufficient.
It can be seen from Fig.10(b) that with the increase of the penetration rate of EVs, the energy of EVs participating in the control of G2V and V2G will increase accordingly. However, the corresponding energy of EVs participating in G2V mode increases much faster than that of V2G mode. This is mainly because when the demand for upward flexible resources is needed to satisfied, the distribution system firstly selects G2V control strategy. Only when G2V control mode fails to fully balance supply and demand, V2G control mode will be used.

VII. CONCLUSIONS
This paper studies the evaluation of distribution system flexibility considering the renewable energy and EVs. The flexibility of distribution system is defined, the interaction between electric vehicle and distribution network is analyzed, and the control strategy of EV is described. Based on the time series direction matrix, the time series flexibility evaluation model of distribution system is established, and the solution algorithm of dichotomy embedded social learning particle swarm optimization is proposed. Through the example verification and analysis, the following conclusions are obtained.
(1) The flexibility evaluation model of distribution system is proposed in this paper, which can be used to evaluate the time series flexibility. The results show that the uncertain variables corresponding to the larger absolute value elements in the time series direction matrix limit the further increase of flexibility evaluation index, and the main period of insufficient flexibility can be judged accordingly.
(2) The flexible resource configuration of distribution system and its time series operation state will affect the flexibility evaluation index of distribution system. From the operation level, reasonable operation strategies can be formulated to release the potential of resource flexibility of distribution system. From the planning level, flexible resource configuration can be optimized to improve the flexibility of distribution system.
(3) Reasonable charging and discharging strategy of EV can greatly improve the flexibility of distribution system. Due to the priority of G2V control is higher than V2G, the energy of EVs participating in G2V mode is much higher than that of V2G mode. The expectation of insufficient downward flexibility of the system decreases with the increase of the penetration rate of EVs. When the penetration rate is 30-40%, the minimum value of IUFE is 0. When the penetration rate of EVs is 0, the upper and lower limit constrain of the generator unit output make that the upward flexibility of distribution system is better than the downward flexibility.
In the follow-up study, the coupling of uncertain variables can be considered to analyze its impact on the flexibility evaluation. The impact of flexibility demand on the optimization planning or resource configuration of distribution network can be analyzed. In addition, the setting of the optimal charge and discharge price to enhance the enthusiasm of the EV to participate in the flexible control will be discussed.