Smart Charging of EVs in Residential Distribution Systems Using the Extended Iterative Method

Jian Zhang 1, Mingjian Cui 2,*, Hualiang Fang 3 and Yigang He 1,* 1 School of Electrical Engineering and Automation, Hefei University of Technology, Hefei 230009, China; 17775357967@163.com 2 Department of Mechanical Engineering, University of Texas at Dallas, Richardson, TX 75080, USA 3 School of Electrical Engineering, Wuhan University, Wuhan 430072, China; hlfang@whu.edu.cn * Correspondence: mingjian.cui@utdallas.edu (M.C.); 18655136887@163.com (Y.H.); Tel.: +1-469-805-1156 (M.C.); +86-186-551-36887 (Y.H.)


Introduction
The energy sector faces critical challenges worldwide with regard to the security of power supply, environmental impacts, and energy costs [1].Energy investments are trending towards innovations to improve the energy efficiency and the environmental friendliness.EVs present significant benefits over traditional vehicles with regards to their non-reliance on oil, reduced harmful gas emissions, and lowering fluctuations of renewable energy sources.Recently, a number of plans were put forward in China, e.g., the "Ten Cities, Thousand Vehicles" program, to promote the development of EVs, which have made the markets for EVs develop rapidly [2,3].However, uncoordinated charging of massive numbers of EVs could significantly increase the network losses, reduce the energy efficiency, lower voltages and overload distribution transformers or lines, while smart charging of EVs can greatly improve the economic benefits and ensure the safe operation of the distribution system.
In [7,8], a real-time smart load management strategy is proposed for the coordinated charging of EVs based on sensitivity analysis techniques.However, the control variables are the charging locations rather than the charging power of EVs.In [9,10], a genetic algorithm (GA) is used to optimize the coordinated charging of EVs.Due to the inherent characteristics of the GA it takes too much computational time for the coordinated charging of a large population of EVs.
Since the coordinated charging of EVs is a large scale optimization problem, many techniques are proposed to improve the calculation speed.In [11,12], a linear constrained convex quadratic programming method is used to correct the nodal voltage iteratively.However, there are some disadvantages in the aforementioned two methods.The objective function is the minimization of power losses.The conventional household load model is assumed to be constant power load.However, if the conventional household load is not constant power load and the objective function is associated with nodal voltages, the method utilized in [11,12] would not be applicable.Moreover, the imbalance of distribution systems, constraints on nodal voltages, and thermal loadings of lines and transformers are not taken into account.
In [13,14], a linear programming of the coordinated charging of EVs is proposed to maximize the total charging energy of EVs.Inequality constraints of nodal voltage and thermal loadings of transformers and lines are all considered and linearized in the model.However, this model cannot be applicable to the nonlinear objective function which is not linearly related to the charging power of EVs, such as the minimization of total power supply.In [15], a mixed integer linear programming of coordinated charging of EVs is proposed to maximize the revenue of power corporations.Since the charging power of EVs is not optimized, the results are not optimal.
In [16], a quadratic programming approach is proposed to optimize the charging and discharging power of EVs with the time-of-use power price and battery degradation costs considered.In [17], a coordination strategy for optimal charging of EVs considering the congestion of distribution system is proposed.In [18], quadratic programming is used to minimize the power losses of the distribution system.Three-phase photovoltaic inverters and EV chargers are adopted to transfer power from the highly loaded phase to the less loaded phase.In [19], load factor, load variance and network losses are proven to be equivalent under certain conditions.Minimizing network losses can be transformed into minimizing load factor or load variance so as to reduce calculation complexity.However, the constraints on nodal voltages or thermal loadings of transformers and lines are not taken into account in the aforementioned four models.When massive numbers of EVs are connected to the distribution network, the constraints on nodal voltages and/or thermal loadings of transformers and lines can really be a factor that limits the charging power of EVs.Neglecting the constraints on nodal voltages and/or thermal loadings of transformers and lines can greatly improve the calculation speed, but may result in the charging power of EVs being unfeasible.
In this paper, we extend the method in [11,12] to a more general situation to improve the computational efficiency.The objective function is not confined to the minimization of power losses and the conventional household load model is no longer confined to the constant power load.The imbalance of distribution systems, constraints on nodal voltages and thermal loading of transformers and lines are all considered.The accuracy and computational efficiency are compared with selected sophisticated methods.
The organization of this paper is as follows: three EV smart charging models are presented in Section 2. The method that extends the iterative technique to a more general situation is introduced in Section 3. The accuracy and computational efficiency of the proposed method are presented in Section 4. Section 5 concludes the paper.

Basic Assumptions
A schematic diagram of a simple distribution system is illustrated in Figure 1.This system can be reduced to the three-phase three-wire system using the Kcron's reduction.

Basic Assumptions
A schematic diagram of a simple distribution system is illustrated in Figure 1.This system can be reduced to the three-phase three-wire system using the Kcron's reduction.

S
is the three-phase apparent power of the load at node N .In order to reasonably simplify the degree of complexity, some basic calculation conditions are set as follows: (1) Under normal circumstances, the distribution system uses the radial operation structure.Thus, it is suitable to represent it with node branch incident matrix.Evidently, the number of nodes in the radial network is one more than that of the branches (the grounded branches are not considered).Suppose the element Aij in the node branch incident matrix A corresponds to the ith node, jth branch, which is directed from node m to n.Then Aij can be represented as: Suppose the vector of currents injected into nodes in α (α = a, b, c) phase (with the power supply node excluded) is represented as    , and the vector of currents for the branch of α phase is represented as    .Then it has: where N is the total number of nodes of the distribution network excluding the power supply nodes.L is the total number of branches of the distribution network excluding the grounded branches.Thus,    can be calculated as: (2) Suppose the charging location of each EV is fixed.Only the optimization of the charging power is considered while the optimization of the charging location is not considered.(3) The fluctuation of the external power grid is not taken into account.The capacity of external power grid is assumed to be large enough and the substation bus is taken as the slack node in every optimization period.V g is the ground voltage.I a N,N+1 , I b N,N+1 , I c N,N+1 , and I n N,N+1 are the currents of phases a, b, c, and n flowing from node N to node N + 1, respectively.Z aa N , Z bb N , Z cc N , and Z nn N are self-impedances for phase a, b, c and n on the line between node N + 1 and node N, respectively.S abc LN is the three-phase apparent power of the load at node N. In order to reasonably simplify the degree of complexity, some basic calculation conditions are set as follows: (1) Under normal circumstances, the distribution system uses the radial operation structure.Thus, it is suitable to represent it with node branch incident matrix.Evidently, the number of nodes in the radial network is one more than that of the branches (the grounded branches are not considered).Suppose the element A ij in the node branch incident matrix A corresponds to the ith node, jth branch, which is directed from node m to n.Then A ij can be represented as: Suppose the vector of currents injected into nodes in α (α = a, b, c) phase (with the power supply node excluded) is represented as I α N , and the vector of currents for the branch of α phase is represented as Then it has: where N is the total number of nodes of the distribution network excluding the power supply nodes.L is the total number of branches of the distribution network excluding the grounded branches.Thus, I α L can be calculated as: (2) Suppose the charging location of each EV is fixed.Only the optimization of the charging power is considered while the optimization of the charging location is not considered.
(3) The fluctuation of the external power grid is not taken into account.The capacity of external power grid is assumed to be large enough and the substation bus is taken as the slack node in every optimization period.(4) It is assumed that the output power of the charger can be adjusted continuously.

Objective Functions and Constraints
The objective function of the first model is the power supply minimization of the distribution system shown as Equations ( 4) and ( 5) where t 1 and t max are the optimization start and end time, respectively.P P,n,α,t , P I,n,α,t , and P Z,n,α,t are the active power of the conventional household constant impedance, constant current and constant power load of node n phase α at time t, respectively.R l,α and I r l,α,t + jI i l,α,t are the resistance and current of line l phase α at time t, respectively.P EVk,α,t and P EVm,β,t are the charging power of kth, mth EV with three-phase and single phase charging mode, located at phase α, β at time t, respectively.K and M are the total number of EVs with three-phase and single phase charging mode, respectively.
Due to the complexity of integral computation, in current research and practical control, the discretization technique is adopted.The total optimization time is divided into several periods.In each period, it is assumed that the conventional household load does not change and the charging power of each EV keeps constant.Thus, the objective function can be transformed as: where ∆t is the interval of sampling and control.The objective function of the second model is the minimization of energy costs that distribution system operator (DSO) purchases from the external power grid shown as (7): where ρ(t) is the power price that the operator (DSO) purchases from outside at time t.
The objective function of the third model is the maximization of profit for DSO shown as: where φ(t) is defined as 1.0 × 2 6+24H(t−18)−t .λ is the penalty price that the DSO pays to the EVs holder, because the EVs are not fully charged.Suppose E t end j is the energy of jth EV when the optimization time is over and is the capacity of the battery for the jth EV.
is the energy cut down for the jth EV.The first term is the profits that the EVs are charged at early evening which is encouraged by the DSO.H(t) is the Heaviside step function.φ(t) is a very large value at early evening.The more energy EVs are charged at early evening, the more profits the DSO can make.
Constraint on the charging power of each EV with three-phase charging mode is: P EVk,a,t = P EVk,b,t = P EVk,c,t (10) Constraint on the charging power of each EV with single-phase charging mode is: where P EVk,max and P EVm,max are the maximal charging power of the kth and mth EV with three-phase and single phase charging mode, respectively.Constraints on the power demand of each EV with three-phase and single phase charging mode for the first and second models are: where η is the charging efficiency.
, and E cap m are the initial energy and the battery capacity of the kth, mth EV with three-phase and single-phase charging mode, respectively.t ks , t ke , t ms , and t me are the charging start and end time for the kth, mth EV with three-phase and single-phase charging mode, respectively.
As for the third model constraints on power demand are: Constraint on nodal voltage of the distribution system is: where U n,α,t , U min , and U max are the voltage of node n, phase α at time t and its lower and upper limits, respectively.Constraint on thermal loadings of each line is: where S l,α,t and S l,max are the power of the line l phase α at time t and its maximum, respectively.Constraint on thermal loadings of each transformer is: where S Tn,α,t and S Tmax are the apparent power of the distribution transformer for phase α at time t and its maximum, respectively.As for the first and second models, the constraints are Equations ( 9)- (18).As for the third model, the constraints are Equations ( 9)- (18).
In this paper, the technique of sensitivity analysis presented in [13] is introduced to linearize the constraints on nodal voltages and thermal loadings of transformers and lines.

Extended Iterative Method
Nodal power balance equations are: . .
S n,α,t = (P P,n,α,t + P I,n,α,t + P Z,n,α,t + P EVn,α,t ) + j(Q P,n,α,t where Q P,n,α,t , Q I,n,α,t , and Q Z,n,α,t are the reactive power of ZIP load on node n phase α at time t, respectively.P I0,n,α,t and P Z0,n,α,t are the active power of constant current and power conventional household load when the voltage on node n phase α at time t is U 0,n,α,t , respectively.The current can be obtained from Equation (19), the formula for which is shown as: .
It can be seen from Equations ( 3), ( 20)-( 23) that the current of each line and its conjugate are the linear function of the charging power of EVs, assuming that the voltages are fixed to some known values.As the branch current magnitude square is equal to the current multiplied by its conjugate and power loss is a linear function of the current magnitude square, the power loss is the quadratic function of the charging power.However, when the conventional household load contains constant impedance and/or constant current load, the first term of J 0 in Equation ( 5) is associated with nodal voltage.The conventional iterative methods in [11,12] cannot be applicable any more.To solve this problem, an extended iterative method is formulated as follows: (P1) Using the slack node as the root node, the depth first search (DFS) program is used to calculate the forward node sequence Pre and the parent node sequence Pred.(P2) At each iteration, using the complex voltage vector of root node, the branch impedance matrix and the branch current complex vector, node sequence Pre, Pred, the predicted voltage of each node per phase at time t related to charging power can successively be calculated as: .
U abc pre(i) and .
U abc pred(pre(i)) are the three-phase predicted complex voltage vector of node Pre(i) and its parent node Pred(Pre(i)), respectively.Z pre(i),pred(pre(i)) and .I abc pre(i),pred(pre(i)) are the impedance matrix and three-phase complex branch current vector between node Pre(i) and its parent node Pred(Pre(i)), respectively.(P3) The square of the predicted voltage magnitude of each node per phase at each time t is calculated with the complex voltage multiplied by its conjugate.Voltage magnitude can be obtained by root calculation.(P4) The active and reactive power of the conventional household load for each node per phase at each time t can be calculated by substituting the predicted voltage magnitudes into Equations ( 21) and (22).
By the above procedures P1~P4, the iterative method can be applicable once more.Clearly, the predicted voltage complex vector and its conjugate are linear functions of the charging power, such as Equation (24).Therefore, the square of the predicted voltage magnitude is the quadratic function of the charging power.If the conventional household load is only composed of constant power and/or constant impedance load, the active power of conventional household load is the quadratic function of the charging power.While if the conventional household load contains constant current load, the active power of the conventional household load is no longer the quadratic function of the charging power.But it is still a convex function of the charging power.
The flow chart of the proposed method is shown in Figure 2. The basic input information includes the topological structure of distribution system, parameters of distribution transformer and lines, load curves of each node per phase, charging locations, charging power limit, battery capacity, initial state of charge (SOC), charging start and end time of each EV, etc.  .When the difference of the infinity norm of a two adjacent optimal charging power vector is less than a pre-given small positive number, the procedure is terminated and the results are output.

Simulation Case 1
A simple distribution system with two nodes is used to validate the precision of the proposed n,α,t can be calculated using power flow calculation.In this paper, the back/forward sweep method is used for power flow calculation.The node voltage .

U (k+1)
n,α,t is then used as the known quantity, and a linear constrained convex quadratic programming or convex programming model is formulated at each iteration, and P (k+1) EVj,α,t is obtained.The power flow is performed again with the charging power of each EV P (k+1) EVj,α,t used as the input information.The convergence criterion of the whole program is to maximize When the difference of the infinity norm of a two adjacent optimal charging power vector is less than a pre-given small positive number, the procedure is terminated and the results are output.

Simulation Case 1
A simple distribution system with two nodes is used to validate the precision of the proposed method.The single line diagram is shown in Figure 3.The base power and base voltage of each phase are 160/3 kVA and 10/ √ 3 kV, respectively.The voltage source is 1.05 p.u.The impedance matrix of the transmission lines is equal to: Energies 2016, 9, 985 8 of 20 Line 0.0276 + 0.0124i -0.0056 + 0.0060i -0.0056 + 0.0060i -0.0056 + 0.0060i 0.0276 + 0.0124i -0.0056 + 0.0060i -0.0056 + 0.0060i -0.0056 + 0.0060i 0.0276 + 0.0124i The optimization period is 12:00~14:00.During the optimization period, there is a 40 kW and 80 kW conventional household load connected to bus 1 during the periods of 12:00~13:00 and 13:00~14:00, respectively.The power factor is 0.95.There are 10, 14 and 16 EVs connected to phase A, B, and C of bus 1, respectively.The power demand of each EV is 10 kWh.The upper limit of charging power of each EV is 10 kW.The upper and lower limit of the voltage are 1.10 p.u. and 0.90 p.u., respectively.
As for the objective function J1 and different load models, the optimization results of the proposed iterative method denoted as PM.0 and PM.1 are shown in Table 1, where f.0 and f.1 represent the optimization functions of MATLAB, fmincon, for charging periods of 12:00~13:00 and 13:00~14:00, respectively.Evidently, the optimization results of voltages and objective function are very close between PM.0 and f.0, as well as between PM.1 and f.1.However, the difference of optimal charging power between the proposed method and optimization functions is somewhat remarkable for the constant impedance load model.The optimization period is 12:00~14:00.During the optimization period, there is a 40 kW and 80 kW conventional household load connected to bus 1 during the periods of 12:00~13:00 and 13:00~14:00, respectively.The power factor is 0.95.There are 10, 14 and 16 EVs connected to phase A, B, and C of bus 1, respectively.The power demand of each EV is 10 kWh.The upper limit of charging power of each EV is 10 kW.The upper and lower limit of the voltage are 1.10 p.u. and 0.90 p.u., respectively.
As for the objective function J 1 and different load models, the optimization results of the proposed iterative method denoted as PM.0 and PM.1 are shown in Table 1, where f.0 and f.1 represent the optimization functions of MATLAB, fmincon, for charging periods of 12:00~13:00 and 13:00~14:00, respectively.Evidently, the optimization results of voltages and objective function are very close between PM.0 and f.0, as well as between PM.1 and f.1.However, the difference of optimal charging power between the proposed method and optimization functions is somewhat remarkable for the constant impedance load model.

Simulation Case 2
A distribution system with 14 nodes is used to validate the precision and the calculation speed of the proposed method.The single line diagram is shown in Figure 4.The parameters of the lines are shown in Table 2.The capacity of the transformer is 400 kVA.Both the positive and zero sequence impedance of the transformer are 0.06 + i0.0125 p.u.The base power and voltage of each phase are 160/3 kVA, 10/ √ 3 kV and 0.4/ √ 3 kV, respectively.Node 14 is the slack node and its voltage is 1.05 p.u.The optimization period is 12:00~14:00.During the periods of 12:00~13:00 and 13:00~14:00, the load connected to each node per phase is 3 kW and 1.5 kW, respectively.The power factor is 0.95.There is one EV connected to node 3~5, 8~10 phase A, node 7, 12 phase B, node 6, 11 phase C, respectively.The power demand of each EV connected to phase A, B, C are 10 kWh, 25 kWh, and 25 kWh, respectively.The upper limits of charging power of each EV connected to phase A, B, and C are 10 kW, 25 kW, and 25 kW, respectively.The upper and lower limits of the voltage are 1.10 p.u. and 0.90 p.u., respectively.

Simulation Case 2
A distribution system with 14 nodes is used to validate the precision and the calculation speed of the proposed method.The single line diagram is shown in Figure 4.The parameters of the lines are shown in Table 2.The capacity of the transformer is 400 kVA.Both the positive and zero sequence impedance of the transformer are 0.06 + i0.0125 p.u.The base power and voltage of each phase are 160/3 kVA, 10 / 3 kV and 0.4 / 3 kV, respectively.Node 14 is the slack node and its voltage is 1.05 p.u.The optimization period is 12:00~14:00.During the periods of 12:00~13:00 and 13:00~14:00, the load connected to each node per phase is 3 kW and 1.5 kW, respectively.The power factor is 0.95.There is one EV connected to node 3~5, 8~10 phase A, node 7, 12 phase B, node 6, 11 phase C, respectively.The power demand of each EV connected to phase A, B, C are 10 kWh, 25 kWh, and 25 kWh, respectively.The upper limits of charging power of each EV connected to phase A, B, and C are 10 kW, 25 kW, and 25 kW, respectively.The upper and lower limits of the voltage are 1.10 p.u. and 0.90 p.u., respectively.Table 2.The parameters of distribution system.

Line L(m)
  As for the objective function J1 and different load models, the optimization results are shown as Tables A1-A3 in the Appendix A. The calculation results of voltages and objective functions are very close for PM.0 and f.0, as well as for PM.1 and f.1.However, the difference between optimal charging   capability of the proposed method.Length, impedance and rated current of lines are shown in Table 2.The transformer's rated capacity is 400 kVA with both positive and zero sequence impedance are 0.06 + i0.0125 p.u.As shown by the arrows in Figure 5, there are total 165 households in the distribution system.The lumped load is connected to phase A. It indicates 31 single phase household loads of another area.While the rest 134 household loads indicate three-phase symmetrical ones.As shown by the small black circles in Figure 5, there are 67 EVs distributed in the distribution system.Taking one EV per household for example, the permeability of EV is 50%.Charging locations of each EV is the same as the household load.As the load curve of each household is not available, typical summer load curves between 18:00~7:00 as shown in Figure 6 is assigned to each household.Different load curves may result in different optimal charging power curves for each EV.However, the capability of the algorithm is independent of load curves.Other simulation conditions are set as follows:

Line
(1) All EV owners are willing to participate in coordinated charging and charging power of each EV is fully controllable.Optimization period of time is between 18:00~7:00.(2) Maximal charging power of each EV is 4 kW.The battery capacity of each EV is 20 kWh.
Charging efficiency is 0.98.As for J1, J2, J3, initial SOC of each EV is set to be 0.25, 0.5, and 0.5, respectively.(3) Each EV adopts single phase charging mode.EVs located at area 1, 2 and 3 are connected to phase A, B, and C, respectively.(4) The optimization time interval is 1 h.( 5) The upper and lower voltage limit are 1.1 p.u. and 0.9 p.u., respectively.( 6) Node 141 is taken as the slack node, and its voltage is kept constant as 1.05 p.u.The rest of the nodes are taken as PQ nodes.

Simulation Results
Typically, the program converges within five iterations.Charging power of some EVs located at typical nodes for different objective functions are shown in Figures 7-9.It can be seen that the charging power of EVs located at different nodes vary from each other greatly.As for the objective function J 1 , in the evening peak time, the charging power of EVs is large to reduce the conventional household load power by lowering the voltage.As for EV located at node 15, the charging power is maximal in the evening peak time while it is minimal at late night trough time.This is because node 15 is near to the start of the distribution system, the network losses of additional current caused by this charging load is low.Charging with a higher power at the peak time can lower the voltage.Thus, the power of the conventional household load can be reduced.As for EVs located at node 91, 140, 30, 108, 126, the charging power at peak time is lower than that at trough time.This is because these nodes are far from the start of the distribution system.Increasing the charging power at trough time can effectively reduce the network losses of the additional current caused by the charging power.Charging power is also closely associated with the power price.As for the objective function J2, when the power price is low, the charging power is high.Conversely, when the power price is high, 18:00 20:00 22:00 24:00 2:00 4:00 6:00 0.7 0.8 0.9    Charging power is also closely associated with the power price.As for the objective function J2, 18:00 20:00 22:00 24:00 2:00 4:00 6:00 0.7 0.8 0.9 Charging power is also closely associated with the power price.As for the objective function J 2 , when the power price is low, the charging power is high.Conversely, when the power price is high, the charging power is low (even zero).As for objective function J 3 , the charging power of the EVs with priority such as those connected to node 15, 140, and 126 is very high during peak time in order to get fully charged as soon as possible for use in the early evening.As for the EVs without priority such as those connected to node 91, 30, 108, the charging power is high when the power price is low.Conversely, when the power price is high, the charging power of those EVs is low (even zero).
Total charging power of EVs, total conventional household load of household, total power losses and total power supply of the distribution system for different objective functions are shown in Figure 10.As for the objective function J 1 , the total charging power is somewhat uniform.At the peak time, it is slightly higher than that at the trough time.The total power supply is similar to that of the total conventional load.As for objective J 2 , the total charging load is concentrated on the time when the power price is low.In the peak power price time, the total charging power is very small, even zero.As for objective J 3 , the total charging load is concentrated in the early evening and the low power price time.Charging power is also closely associated with the power price.As for the objective function J2, when the power price is low, the charging power is high.Conversely, when the power price is high, the charging power is low (even zero).As for objective function J3, the charging power of the EVs with priority such as those connected to node 15, 140, and 126 is very high during peak time in order 18:00 20:00 22:00 24:00 2:00 4:00 6:00 0 0.5 to get fully charged as soon as possible for use in the early evening.As for the EVs without priority such as those connected to node 91, 30, 108, the charging power is high when the power price is low.Conversely, when the power price is high, the charging power of those EVs is low (even zero).
Total charging power of EVs, total conventional household load of household, total power losses and total power supply of the distribution system for different objective functions are shown in Figure 10.As for the objective function J1, the total charging power is somewhat uniform.At the peak time, it is slightly higher than that at the trough time.The total power supply is similar to that of the total conventional load.As for objective J2, the total charging load is concentrated on the time when the power price is low.In the peak power price time, the total charging power is very small, even zero.As for objective J3, the total charging load is concentrated in the early evening and the low power price time.Thermal loadings of transformer for different objective functions are shown in Figure 11.Clearly, it increases significantly during the periods when the charging power is high.But all of them are less than 80%.
18:00 20:00 22:00 24:00 2:00 4:00 6:00 0 0.5 Thermal loadings of transformer for different objective functions are shown in Figure 11.Clearly, it increases significantly during the periods when the charging power is high.But all of them are less than 80%.Thermal loadings of the main cable for different objective functions are shown in Figure 12.Clearly, it increases significantly during the periods when the charging power is high, but all of them are less than 85%.It is evident that neither the currents of transformer nor those of main cable are the binding constraints on this network.Clearly, the network equipment is more than well suited to accommodate the additional load required by the high penetration of EVs, assuming the proposed coordinated charging scheme is introduced.Thermal loadings of the main cable for different objective functions are shown in Figure 12.Clearly, it increases significantly during the periods when the charging power is high, but all of them are less than 85%.It is evident that neither the currents of transformer nor those of main cable are the binding constraints on this network.Clearly, the network equipment is more than well suited to accommodate the additional load required by the high penetration of EVs, assuming the proposed coordinated charging scheme is introduced.
Clearly, it increases significantly during the periods when the charging power is high, but all of them are less than 85%.It is evident that neither the currents of transformer nor those of main cable are the binding constraints on this network.Clearly, the network equipment is more than well suited to accommodate the additional load required by the high penetration of EVs, assuming the proposed coordinated charging scheme is introduced.The three-phase voltage profiles of node 6, which is located at the terminal of the distribution system are shown as Figure 13.As the lumped load is connected to phase A, the voltage of phase A for node 6 is slightly lower than that of phase B and C when charging power is low.Clearly, the high charging load results in a significant drop in voltage.If the voltage constraint is not applied, the voltage of phase C is much lower than that of the lower bound.However, if the voltage constraint is applied, the voltage of all the three phase is above that of the lower bound.Thus, the voltage is really a binding constraint when there are a large amount of EVs connected to the distribution network.
(a) (b)  (b) Thermal loadings of main cable for J 2 ; (c) Thermal loadings of main cable for J 3 .
The three-phase voltage profiles of node 6, which is located at the terminal of the distribution system are shown as Figure 13.As the lumped load is connected to phase A, the voltage of phase A for node 6 is slightly lower than that of phase B and C when charging power is low.Clearly, the high charging load results in a significant drop in voltage.If the voltage constraint is not applied, the voltage of phase C is much lower than that of the lower bound.However, if the voltage constraint is applied, the voltage of all the three phase is above that of the lower bound.Thus, the voltage is really a binding constraint when there are a large amount of EVs connected to the distribution network.
system are shown as Figure 13.As the lumped load is connected to phase A, the voltage of phase A for node 6 is slightly lower than that of phase B and C when charging power is low.Clearly, the high charging load results in a significant drop in voltage.If the voltage constraint is not applied, the voltage of phase C is much lower than that of the lower bound.However, if the voltage constraint is applied, the voltage of all the three phase is above that of the lower bound.Thus, the voltage is really a binding constraint when there are a large amount of EVs connected to the distribution network.

Comparison with the Selected Method
The optimization results of the proposed method compared with the method in [15] are shown as Table 3.Clearly, the calculation efficiency of the proposed method is much higher than the method described in [15].This is because at each iteration, optimization variables only consist of the charging power of EVs, which greatly reduce the number of optimization variables.Whereas, in [15] nodal voltages are included in the optimization variables.Although linearization techniques are applied and mixed integer linear programming is formulated, because the number of optimization variables is large, the calculation time increases greatly.From the simulation results, it can be seen that the proposed method has a rapid convergence rate in the optimization of large-scale coordinated charging of EVs, and can greatly improve the economic benefits and guarantee the safe operation of distribution system.

Conclusions
In this paper, three models for smart charging of EVs were formulated to minimize the total power supply, power costs and maximize the operation profits.The iterative method was extended to solve the developed models.The extended algorithm was accommodated to a mix of ZIP load of the conventional household with the three-phase imbalance of distribution network, constraints on nodal voltages and thermal loadings of transformers and lines taken into account.When the conventional household load was composed of constant impedance and/or constant power loads, a linearly constrained convex quadratic programming model could be formulated at each iteration.When the conventional household load consists of the constant current load, a linear constrained convex programming model could be formulated instead.
The simulation results indicate that: (1) Compared with the optimization results in [15], the developed method had much higher calculation efficiency.(2) This developed method could avoid the risk of the node voltage exceeding the lower limit, and provided reliable and economic operation of the distribution system when massive numbers of EVs penetrate into the grid.

Figure 1 .I
Figure 1.Schematic diagram of a simple distribution system.

Z
are self-impedances for phase a , b , c and n on the line between node 1 N  and node N , respectively.abcLN

Figure 1 .N
Figure 1.Schematic diagram of a simple distribution system.

Figure 2 .
Figure 2. Flow chart of the proposed algorithm for smart charging of EVs.

Figure 3 .
Figure 3. Single line diagram of a 2-node distribution system.

Figure 4 .
Figure 4. Single line diagram of a 14-node distribution system.

Figure 4 .
Figure 4. Single line diagram of a 14-node distribution system.

Figure 5 .
Figure 5.The single line diagram of an actual distribution system with 141 nodes.Figure 5.The single line diagram of an actual distribution system with 141 nodes.

Figure 5 . 20 Figure 6 .
Figure 5.The single line diagram of an actual distribution system with 141 nodes.Figure 5.The single line diagram of an actual distribution system with 141 nodes.Energies 2016, 9, 985 11 of 20

Figure 7 .
Figure 7. Optimal charging power of EVs at typical nodes for the objective function J1.

Figure 8 .
Figure 8. Optimal charging power of EVs at typical nodes for J2.

Figure 9 .
Figure 9. Optimal charging power of EVs at typical nodes for J3.

Figure 7 .Figure 7 .
Figure 7. Optimal charging power of EVs at typical nodes for the objective function J 1 .

Figure 8 .
Figure 8. Optimal charging power of EVs at typical nodes for J2.

Figure 9 .
Figure 9. Optimal charging power of EVs at typical nodes for J3.

Figure 8 .
Figure 8. Optimal charging power of EVs at typical nodes for J2.

Figure 9 .
Figure 9. Optimal charging power of EVs at typical nodes for J3.

Figure 10 .
Figure 10.Relationship between charging load and conventional household load for J1, J2, and J3.(a) Relationship between charging load and conventional household load for J1; (b) Relationship between charging load and conventional household load for J2; (c) Relationship between charging load and conventional household load for J3.

Figure 10 .
Figure10.Relationship between charging load and conventional household load for J 1 , J 2 , and J 3 .(a) Relationship between charging load and conventional household load for J 1 ; (b) Relationship between charging load and conventional household load for J 2 ; (c) Relationship between charging load and conventional household load for J 3 .

Figure 11 .
Figure 11.Thermal loadings of transformer for J 1 , J 2 , and J 3 .(a) Thermal loadings of transformer for J 1 ; (b) Thermal loadings of transformer for J 2 ; (c) Thermal loadings of transformer for J 3 .

Figure 12 .
Figure 12.Thermal loadings of main cable for J1, J2, and J3.(a) Thermal loadings of main cable for J1; (b) Thermal loadings of main cable for J2; (c) Thermal loadings of main cable for J3.

Figure 12 .
Figure12.Thermal loadings of main cable for J 1 , J 2 , and J 3 .(a) Thermal loadings of main cable for J 1 ; (b) Thermal loadings of main cable for J 2 ; (c) Thermal loadings of main cable for J 3 .

Figure 13 .
Figure 13.Voltage of node 6 under different scenarios for J 1 , J 2 , and J 3 .(a) Voltage of node 6 under different scenarios for J 1 ; (b) Voltage of node 6 under different scenarios for J 2 ; (c) Voltage of node 6 under different scenarios for J 3 .
Flow chart of the proposed algorithm for smart charging of EVs.
P  used as the input information.The convergence criterion of the whole program is to maximize

Table 1 .
Optimization results for different load models.Single line diagram of a 2-node distribution system.

Table 1 .
Optimization results for different load models.

Table 2 .
The parameters of distribution system.
(9) model of the conventional household load is assumed to be constant impedance load.(9)Thesingle phase based power and voltage of the system is 160/3 kVA, 10/ Optimal charging power of EVs at typical nodes for J 2 .
Optimal charging power of EVs at typical nodes for J 3 .