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

Abstract

Smart charging of electrical vehicles (EVs) is critical to provide the secure and cost-effective operation for distribution systems. Three model objective functions which are minimization of total supplied power, energy costs and maximization of profits are formulated. The conventional household load is modeled as a ZIP load that consists of constant power, constant current and constant impedance components. The imbalance of distribution system, constraints on nodal voltages and thermal loadings of lines and transformers are all taken into account. Utilizing the radial operation structure of distribution system, an extended iterative method is proposed to greatly reduce the dimensions of optimization variables and thus improve calculation speed. Impacts of the conventional household load model on the simulation results are also investigated. Case studies on three distribution systems with 2, 14, and 141 buses are performed and analyzed. It is found that the linear constrained convex quadratic programming model is applicable at each iteration, when the conventional household load is composed of constant power and constant impedance load. However, it is not applicable when the conventional household load consists of constant current load. The accuracy and computational efficiency of the proposed method are also validated.

1. 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.

Research on coordinated charging of EVs is generally focused on distributed and centralized coordinated charging technology. The distributed coordinated charging technology mainly uses fuzzy mathematics theory [4,5] and sensitivity analysis [6]. The centralized coordinated charging technology generally adopts sensitivity analysis [7,8], simulated evolutionary algorithms [9,10], or optimization techniques [11,12,13,14,15,16,17,18,19].

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.

2. Smart Charging Models of Electrical Vehicles (EVs)

2.1. 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.

In this figure, ${\dot{V}}_N^{abc}$ is the three-phase voltage; ${\dot{V}}_N^n$ is the neutral voltage of node N; and ${\dot{V}}_g$ is the ground voltage. $I_{N,N+1}^a$, $I_{N,N+1}^b$, $I_{N,N+1}^c$, and $I_{N,N+1}^n$ are the currents of phases $a$, $b$, $c$, and $n$ flowing from node $N$ to node $N+1$, respectively. $Z_N^{aa}$, $Z_N^{bb}$, $Z_N^{cc}$, and $Z_N^{nn}$ are self-impedances for phase $a$, $b$, $c$ and $n$ on the line between node $N+1$ and node $N$, respectively. $S_{LN}^{abc}$ 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:

$$A_{ij}=\{\begin{array}{l}1\\ -1\\ 0\end{array}\begin{array}{c}i=m\\ i=n\\ i\ne m,i\ne n\end{array}$$

(1)

Suppose the vector of currents injected into nodes in α (α = a, b, c) phase (with the power supply node excluded) is represented as $I_N^{\alpha }$, and the vector of currents for the branch of α phase is represented as $I_L^{\alpha }$. Then it has:

$$I_N^{\alpha }=A{I}_L^{\alpha }$$

(2)

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^{\alpha }$ can be calculated as:

$$I_L^{\alpha }=A^{-1}{I}_N^{\alpha }$$

(3)

(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.

2.2. 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)

$$J_1=\min {\displaystyle {\int }_{t_1}^{t_{\max }}{J}_0}\left(t\right)dt$$

(4)

$$J_0\left(t\right)={\displaystyle \sum _{n=1}^N{\displaystyle \sum _{\alpha =1}^3\left(P_{\mathrm{P},n,\alpha ,t}+P_{\mathrm{I},n,\alpha ,t}+P_{\mathrm{Z},n,\alpha ,t}\right)}}+{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{\alpha =1}^3{R}_{l,\alpha }\left(I_{l,\alpha ,t}^{\mathrm{r}2}+I_{l,\alpha ,t}^{\mathrm{i}2}\right)}}+{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{\alpha =1}^3{P}_{\mathrm{EV}k,\alpha ,t}+{\displaystyle \sum _{m=1}^M{P}_{\mathrm{EV}m,\beta ,t}}}}$$

(5)

where t₁ 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_{l,\alpha ,t}^r+\mathrm{j}{I}_{l,\alpha ,t}^i$ 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:

$$J_1=\min {\displaystyle \sum _{t=t_1}^{t_{\max }}{J}_0}\left(t\right)\Delta t$$

(6)

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):

$$J_2=\min {\displaystyle \sum _{t=t_1}^{t_{\max }}\rho \left(t\right)J_0\left(t\right)}\Delta t$$

(7)

where $\rho \left(t\right)$ 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:

$$J_3=\max {\displaystyle \sum _{t=t_1}^{t_{\max }}\varphi \left(t\right)\left({\displaystyle \sum _{k=1}^K{\displaystyle \sum _{\alpha =1}^3{P}_{\mathrm{EV}k,\alpha ,t}+{\displaystyle \sum _{m=1}^M{P}_{\mathrm{EV}m,\beta ,t}}}}\right)}\Delta t-{\displaystyle \sum _{t=t_1}^{t_{\max }}{J}_0\left(t\right)\Delta t-{\displaystyle \sum _{j=1}^J\beta E_j^{sh}}}$$

(8)

where $\varphi \left(t\right)$ is defined as $1.0\times 2^{6+24\mathrm{H}\left(t-18\right)-t}$. λ is the penalty price that the DSO pays to the EVs holder, because the EVs are not fully charged. Suppose $E_j^{t_{\mathrm{end}}}$ is the energy of jth EV when the optimization time is over and $E_j^{\mathrm{Cap}}$ is the capacity of the battery for the jth EV. $E_j^{sh}=E_j^{\mathrm{Cap}}-E_j^{t_{\mathrm{end}}}$ 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. $\varphi \left(t\right)$ 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:

$$0\le P_{\mathrm{EV}k,\alpha ,t}\le P_{\mathrm{EV}k,\max }$$

(9)

$$P_{\mathrm{EV}k,a,t}=P_{\mathrm{EV}k,b,t}=P_{\mathrm{EV}k,c,t}$$

(10)

Constraint on the charging power of each EV with single-phase charging mode is:

$$0\le P_{\mathrm{EV}m,\beta ,t}\le P_{\mathrm{EV}m,\max }$$

(11)

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:

$$\eta {\displaystyle \sum _{\alpha =1}^3{\displaystyle \sum _{t=t_{k\mathrm{s}}}^{t_{k\mathrm{e}}}{P}_{\mathrm{EV}k,\alpha ,t}\Delta t}}=E_k^{\mathrm{cap}}-E_k^{\mathrm{ini}}$$

(12)

$$\eta {\displaystyle \sum _{t=t_{m\mathrm{s}}}^{t_{m\mathrm{e}}}{P}_{\mathrm{EV}m,\alpha ,t}}\Delta t=E_m^{\mathrm{cap}}-E_m^{\mathrm{ini}}$$

(13)

where η is the charging efficiency. $E_k^{\mathrm{ini}}$, $E_k^{\mathrm{cap}}$, $E_m^{\mathrm{ini}}$, and $E_m^{\mathrm{cap}}$ 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:

$$\eta {\displaystyle \sum _{\alpha =1}^3{\displaystyle \sum _{t=t_{k\mathrm{s}}}^{t_{k\mathrm{e}}}{P}_{\mathrm{EV}k,\alpha ,t}}}\Delta t\le E_k^{\mathrm{cap}}-E_k^{\mathrm{ini}}$$

(14)

$$\eta {\displaystyle \sum _{t=t_{m\mathrm{s}}}^{t_{m\mathrm{e}}}{P}_{\mathrm{EV}m,\alpha ,t}}\Delta t\le E_m^{\mathrm{cap}}-E_m^{\mathrm{ini}}$$

(15)

Constraint on nodal voltage of the distribution system is:

$$U_{\min }\le U_{n,\alpha ,t}\le U_{\max }$$

(16)

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:

$$S_{l,\alpha ,t}\le S_{l,\max }$$

(17)

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:

$$S_{\mathrm{T}n,\alpha ,t}\le S_{\mathrm{T}\max }$$

(18)

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.

3. Extended Iterative Method

Nodal power balance equations are:

$${\dot{U}}_{n,\alpha ,t}{\dot{I}}_{n,\alpha ,t}^{\ast }={\dot{S}}_{n,\alpha ,t}$$

(19)

$${\dot{S}}_{n,\alpha ,t}=\left(P_{\mathrm{P},n,\alpha ,t}+P_{\mathrm{I},n,\alpha ,t}+P_{\mathrm{Z},n,\alpha ,t}+P_{\mathrm{EV}n,\alpha ,t}\right)+\mathrm{j}\left(Q_{\mathrm{P},n,\alpha ,t}+Q_{\mathrm{I},n,\alpha ,t}+Q_{\mathrm{Z},n,\alpha ,t}\right)$$

(20)

$$P_{\mathrm{I},n,\alpha ,t}=P_{\mathrm{I}0,n,\alpha ,t}\left(\frac{U_{n,\alpha ,t}}{U_{0,n,\alpha ,t}}\right)$$

(21)

$$P_{\mathrm{Z},n,\alpha ,t}=P_{\mathrm{Z}0,n,\alpha ,t}{\left(\frac{U_{n,\alpha ,t}}{U_{0,n,\alpha ,t}}\right)}^2$$

(22)

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:

$${\dot{I}}_{n,\alpha ,t}={\left({\dot{S}}_{n,\alpha ,t}/{\dot{U}}_{n,\alpha ,t}\right)}^{\ast }\begin{array}{c}\forall n,\end{array}\forall \alpha ,\forall t$$

(23)

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₀ 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:

$${\dot{U}}_{\mathrm{pre}(i)}^{\mathrm{abc}}={\dot{U}}_{\mathrm{pred}(\mathrm{pre}(i))}^{\mathrm{abc}}-Z_{\mathrm{pre}(i),\mathrm{pred}(\mathrm{pre}(i))}{\dot{I}}_{\mathrm{pre}(i),\mathrm{pred}(\mathrm{pre}(i))}^{\mathrm{abc}} \forall i,$$

(24)

where ${\dot{U}}_{\mathrm{pre}(i)}^{\mathrm{abc}}$ and ${\dot{U}}_{\mathrm{pred}(\mathrm{pre}(i))}^{\mathrm{abc}}$ are the three-phase predicted complex voltage vector of node Pre(i) and its parent node Pred(Pre(i)), respectively. $Z_{\mathrm{pre}(i),\mathrm{pred}(\mathrm{pre}(i))}$ and ${\dot{I}}_{\mathrm{pre}(i),\mathrm{pred}(\mathrm{pre}(i))}^{\mathrm{abc}}$ 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.

Suppose the charging power of each EV $P_{\mathrm{EV}j,\alpha ,t}^{\left(k\right)}$ is known, the voltage of each node per phase at time t, ${\dot{U}}_{n,\alpha ,t}^{\left(k+1\right)}$ can be calculated using power flow calculation. In this paper, the back/forward sweep method is used for power flow calculation. The node voltage ${\dot{U}}_{n,\alpha ,t}^{\left(k+1\right)}$ 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_{\mathrm{EV}j,\alpha ,t}^{\left(k+1\right)}$ is obtained.

The power flow is performed again with the charging power of each EV $P_{\mathrm{EV}j,\alpha ,t}^{\left(k+1\right)}$ used as the input information. The convergence criterion of the whole program is to maximize $\left|P_{\mathrm{EV}j,\alpha ,t}^{\left(k+1\right)}-P_{\mathrm{EV}j,\alpha ,t}^{\left(k\right)}\right|\le \epsilon$. 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.

4. Simulation Case Study

4.1. 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/\sqrt{3}$ kV, respectively. The voltage source is 1.05 p.u. The impedance matrix of the transmission lines is equal to:

$$Z_{\mathrm{Line}}=\left[\begin{array}{rrr}00276+0.0124\mathrm{i}& -0.0056+0.0060\mathrm{i}& -0.0056+0.0060\mathrm{i}\\ -0.0056+0.0060\mathrm{i}& 0.0276+0.0124\mathrm{i}& -0.0056+0.0060\mathrm{i}\\ -0.0056+0.0060\mathrm{i}& -0.0056+0.0060\mathrm{i}& 0.0276+0.0124\mathrm{i}\end{array}\right]$$

(25)

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₁ and different load models, the optimization results of the proposed iterative method denoted as PM.0 and PM.1 are shown in

Table 1. Optimization results for different load models.

Constant Current Load Model
Different Method	PM.0	PM.1	f.0	f.1
$U_{\mathrm{A}}$/p.u.	0.9687	0.9845	0.9696	0.9836
$U_{\mathrm{B}}$/p.u.	0.9544	0.9699	0.9552	0.9691
$U_{\mathrm{C}}$/p.u.	0.9550	0.9705	0.9558	0.9697
$P_{\mathrm{CHA}}$/p.u.	0.7734	1.1016	0.7475	1.1275
$P_{\mathrm{CHB}}$/p.u.	1.1442	1.4808	1.1203	1.5047
$P_{\mathrm{CHC}}$/p.u.	1.3281	1.6719	1.3068	1.6932
$J_1$/p.u.	15.2264		15.2263
Constant impedance load model
Different method	PM.0	PM.1	f.0	f.1
$U_{\mathrm{A}}$/p.u.	0.9634	0.9921	0.9652	0.9902
$U_{\mathrm{B}}$/p.u.	0.9500	0.9779	0.9516	0.9761
$U_{\mathrm{C}}$/p.u.	0.9502	0.9782	0.9519	0.9765
$P_{\mathrm{CHA}}$/p.u.	0.9827	0.8923	0.9277	0.9473
$P_{\mathrm{CHB}}$/p.u.	1.3460	1.2783	1.2954	1.3296
$P_{\mathrm{CHC}}$/p.u.	1.5262	1.4738	1.4793	1.5207
$J_1$/p.u.	14.9703		14.9696
Constant power load model
Different method	PM.0	PM.1	f.0	f.1
$U_{\mathrm{A}}$/p.u.	0.9747	0.9762	0.9748	0.9761
$U_{\mathrm{B}}$/p.u.	0.9595	0.9610	0.9596	0.9609
$U_{\mathrm{C}}$/p.u.	0.9604	0.9619	0.9605	0.9618
$P_{\mathrm{CHA}}$/p.u.	0.5593	1.3157	0.5569	1.3181
$P_{\mathrm{CHB}}$/p.u.	0.9337	1.6913	0.9312	1.6938
$P_{\mathrm{CHC}}$/p.u.	1.1210	1.8790	1.1184	1.8816
$J_1$/p.u.	15.5020		15.5020

, 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.

4.2. 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 parameters of distribution system.

Line	L(m)	${\mathit{R}}_1(\Omega )$	${\mathit{X}}_1(\Omega )$	${\mathit{R}}_0(\Omega )$	${\mathit{X}}_0(\Omega )$	${\mathit{I}}_{\mathbf{N}}\left(A\right)$
MV	1000	20.8	4	10	12	1000
1–2	190	0.0032	0.014	0.095	0.041	510
2–3	27.5	0.008	0.002	0.024	0.006	368
3–4	85	0.024	0.006	0.073	0.018	368
4–5	97.5	0.028	0.007	0.084	0.021	368
5–6	154	0.062	0.011	0.185	0.033	300
4–7	119	0.048	0.009	0.143	0.026	300
2–8	32.5	0.009	0.002	0.028	0.007	368
8–9	59	0.017	0.004	0.051	0.013	368
9–10	106	0.030	0.008	0.091	0.023	368
10–11	95	0.027	0.007	0.082	0.021	368
9–12	217	0.087	0.016	0.261	0.047	368

. 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/\sqrt{3}$ kV and $0.4/\sqrt{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.

As for the objective function J₁ and different load models, the optimization results are shown as

Table A1. Constant impedance load model.

Different Method		PM.0	PM.1	f.0	f.1
1	UA/p.u.	0.9952	1.0097	0.9943	1.0105
	UB/p.u.	0.9983	1.0129	0.9976	1.0136
	UC/p.u.	0.9974	1.0119	0.9966	1.0126
2	UA/p.u.	0.9545	0.9781	0.9530	0.9796
	UB/p.u.	0.9623	0.9860	0.9611	0.9873
	UC/p.u.	0.9607	0.9842	0.9594	0.9856
3	UA/p.u.	0.9498	0.9744	0.9482	0.9760
	UB/p.u.	0.9582	0.9829	0.9569	0.9842
	UC/p.u.	0.9565	0.9809	0.9551	0.9824
4	UA/p.u.	0.9407	0.9668	0.9386	0.9689
	UB/p.u.	0.9466	0.9735	0.9450	0.9750
	UC/p.u.	0.9451	0.9716	0.9434	0.9732
5	UA/p.u.	0.9359	0.9630	0.9336	0.9654
	UB/p.u.	0.9448	0.9735	0.9434	0.9749
	UC/p.u.	0.9327	0.9600	0.9307	0.9619
6	UA/p.u.	0.9364	0.9655	0.9342	0.9677
	UB/p.u.	0.9434	0.9740	0.9421	0.9754
	UC/p.u.	0.9083	0.9353	0.9055	0.9381
7	UA/p.u.	0.9395	0.9671	0.9375	0.9691
	UB/p.u.	0.9280	0.9547	0.9257	0.9569
	UC/p.u.	0.9453	0.9734	0.9438	0.9749
8	UA/p.u.	0.9492	0.9739	0.9476	0.9755
	UB/p.u.	0.9577	0.9825	0.9563	0.9837
	UC/p.u.	0.9561	0.9806	0.9546	0.9820
9	UA/p.u.	0.9428	0.9685	0.9408	0.9705
	UB/p.u.	0.9493	0.9758	0.9479	0.9771
	UC/p.u.	0.9481	0.9741	0.9464	0.9757
10	UA/p.u.	0.9377	0.9645	0.9354	0.9667
	UB/p.u.	0.9474	0.9758	0.9461	0.9771
	UC/p.u.	0.9349	0.9617	0.9329	0.9638
11	UA/p.u.	0.9380	0.9657	0.9359	0.9679
	UB/p.u.	0.9466	0.9758	0.9453	0.9770
	UC/p.u.	0.9244	0.9510	0.9219	0.9535
12	UA/p.u.	0.9408	0.9692	0.9389	0.9710
	UB/p.u.	0.9150	0.9415	0.9127	0.9438
	UC/p.u.	0.9486	0.9774	0.9472	0.9789
13	UA/p.u.	1.0069	1.0175	1.0062	1.0182
	UB/p.u.	1.0096	1.0201	1.0090	1.0207
	UC/p.u.	1.0087	1.0192	1.0080	1.0198
PCH3A/p.u.		0.1023	0.0852	0.0952	0.0923
PCH4A/p.u.		0.0830	0.1043	0.0915	0.0960
PCH5A/p.u.		0.0832	0.1043	0.0910	0.0965
PCH8A/p.u.		0.1000	0.0875	0.0945	0.0928
PCH9A/p.u.		0.0834	0.1041	0.0913	0.0962
PCH10A/p.u.		0.0831	0.1044	0.0913	0.0962
PCH7B/p.u.		0.2208	0.2480	0.2300	0.2387
PCH12B/p.u.		0.2223	0.2464	0.2292	0.2395
PCH6C/p.u.		0.2212	0.2475	0.2290	0.2398
PCH11C/p.u.		0.2208	0.2480	0.2302	0.2385
J1/p.u.		6.6278		6.6275
Calculation time/s		0.736		167

,

Table A2. Constant current load model.

Different Method		PM.0	PM.1	f.0	f.1
1	UA/p.u.	0.9963	1.0064	0.9957	1.0070
	UB/p.u.	0.9955	1.0105	0.9980	1.0110
	UC/p.u.	0.9978	1.0092	0.9972	1.0097
2	UA/p.u.	0.9570	0.9717	0.9559	0.9727
	UB/p.u.	0.9627	0.9819	0.9618	0.9829
	UC/p.u.	0.9613	0.9795	0.9603	0.9805
3	UA/p.u.	0.9525	0.9675	0.9514	0.9686
	UB/p.u.	0.9587	0.9785	0.9577	0.9795
	UC/p.u.	0.9571	0.9760	0.9560	0.9771
4	UA/p.u.	0.9434	0.9594	0.9420	0.9608
	UB/p.u.	0.9476	0.9682	0.9464	0.9694
	UC/p.u.	0.9457	0.9661	0.9444	0.9673
5	UA/p.u.	0.9387	0.9553	0.9372	0.9569
	UB/p.u.	0.9455	0.9684	0.9444	0.9695
	UC/p.u.	0.9338	0.9535	0.9323	0.9550
6	UA/p.u.	0.9386	0.9582	0.9372	0.9596
	UB/p.u.	0.9437	0.9691	0.9426	0.9701
	UC/p.u.	0.9108	0.9268	0.9088	0.9289
7	UA/p.u.	0.9419	0.9599	0.9406	0.9612
	UB/p.u.	0.9307	0.9473	0.9290	0.9490
	UC/p.u.	0.9454	0.9683	0.9443	0.9694
8	UA/p.u.	0.9519	0.9670	0.9508	0.9682
	UB/p.u.	0.9581	0.9782	0.9571	0.9792
	UC/p.u.	0.9568	0.9755	0.9557	0.9767
9	UA/p.u.	0.9455	0.9613	0.9442	0.9626
	UB/p.u.	0.9498	0.9710	0.9487	0.9721
	UC/p.u.	0.9491	0.9684	0.9478	0.9697
10	UA/p.u.	0.9404	0.9570	0.9389	0.9584
	UB/p.u.	0.9475	0.9712	0.9465	0.9723
	UC/p.u.	0.9368	0.9548	0.9353	0.9563
11	UA/p.u.	0.9404	0.9585	0.9390	0.9598
	UB/p.u.	0.9465	0.9713	0.9455	0.9723
	UC/p.u.	0.9272	0.9429	0.9254	0.9448
12	UA/p.u.	0.9430	0.9621	0.9418	0.9634
	UB/p.u.	0.9176	0.9339	0.9157	0.9357
	UC/p.u.	0.9489	0.9723	0.9478	0.9734
13	UA/p.u.	1.0080	1.0148	1.0076	1.0152
	UB/p.u.	1.0100	1.0182	1.0095	1.0187
	UC/p.u.	1.0092	1.0169	1.0088	1.0174
PCH3A/p.u.		0.0856	0.1019	0.0827	0.1048
PCH4A/p.u.		0.0742	0.1133	0.0781	0.1094
PCH5A/p.u.		0.0738	0.1137	0.07871	0.1088
PCH8A/p.u.		0.0851	0.1024	0.0833	0.1042
PCH9A/p.u.		0.0744	0.1131	0.0792	0.1083
PCH10A/p.u.		0.0736	0.11389	0.0778	0.1097
PCH7B/p.u.		0.1942	0.2745	0.2007	0.2681
PCH12B/p.u.		0.2039	0.2648	0.2091	0.2596
PCH6C/p.u.		0.2027	0.2660	0.2084	0.2604
PCH11C/p.u.		0.1951	0.2737	0.2016	0.2671
J1/p.u.		6.7695		6.7693
Calculation time/s		0.812		185

and

Table A3. Constant power load model.

Different Method		PM.0	PM.1	f.0	f.1
1	UA/p.u.	0.9976	1.0025	0.9975	1.0026
	UB/p.u.	0.9989	1.0077	0.9988	1.0078
	UC/p.u.	0.9983	1.0060	0.9982	1.0061
2	UA/p.u.	0.9599	0.9643	0.9599	0.9643
	UB/p.u.	0.9633	0.9771	0.9630	0.9774
	UC/p.u.	0.9621	0.9740	0.9619	0.9742
3	UA/p.u.	0.9556	0.9597	0.9556	0.9597
	UB/p.u.	0.9594	0.9734	0.9591	0.9737
	UC/p.u.	0.9579	0.9703	0.9576	0.9705
4	UA/p.u.	0.9467	0.9509	0.9467	0.9509
	UB/p.u.	0.9487	0.9621	0.9484	0.9624
	UC/p.u.	0.9465	0.9596	0.9462	0.9600
5	UA/p.u.	0.9420	0.9465	0.9420	0.9465
	UB/p.u.	0.9463	0.9624	0.9460	0.9627
	UC/p.u.	0.9352	0.9460	0.9348	0.9464
6	UA/p.u.	0.9413	0.9497	0.9414	0.9496
	UB/p.u.	0.9442	0.9632	0.9439	0.9635
	UC/p.u.	0.9137	0.9169	0.9132	0.9175
7	UA/p.u.	0.9448	0.9515	0.9449	0.9515
	UB/p.u.	0.9337	0.9389	0.9334	0.9393
	UC/p.u.	0.9457	0.9622	0.9454	0.9625
8	UA/p.u.	0.9551	0.9592	0.9550	0.9592
	UB/p.u.	0.9586	0.9731	0.9583	0.9734
	UC/p.u.	0.9577	0.9697	0.9575	0.9699
9	UA/p.u.	0.9487	0.9530	0.9487	0.9530
	UB/p.u.	0.9505	0.9655	0.9502	0.9658
	UC/p.u.	0.9503	0.9618	0.9500	0.9621
10	UA/p.u.	0.9436	0.9484	0.9436	0.9484
	UB/p.u.	0.9478	0.9659	0.9475	0.9662
	UC/p.u.	0.9390	0.9467	0.9387	0.9471
11	UA/p.u.	0.9432	0.9501	0.9433	0.9501
	UB/p.u.	0.9467	0.9660	0.9463	0.9664
	UC/p.u.	0.9304	0.9336	0.9300	0.9340
12	UA/p.u.	0.9457	0.9540	0.9457	0.9540
	UB/p.u.	0.9204	0.9252	0.9199	0.9258
	UC/p.u.	0.9494	0.9662	0.9492	0.9664
13	UA/p.u.	1.0094	1.0116	1.0093	1.0117
	UB/p.u.	1.0105	1.0160	1.0103	1.0161
	UC/p.u.	1.0099	1.0143	1.0098	1.0144
PCH3A/p.u		0.0671	0.1204	0.0682	0.1193
PCH4A/p.u		0.0650	0.1225	0.0632	0.1243
PCH5A/p.u		0.0642	0.1233	0.0651	0.1224
PCH8A/p.u		0.0686	0.1189	0.0701	0.1174
PCH9A/p.u		0.0653	0.1222	0.0658	0.1217
PCH10A/p.u		0.0639	0.1236	0.0628	0.1247
PCH7B/p.u		0.1658	0.3029	0.1667	0.3021
PCH12B/p.u		0.1841	0.2846	0.1858	0.2829
PCH6C/p.u		0.1829	0.2859	0.1845	0.2843
PCH11C/p.u		0.1676	0.3011	0.1687	0.3000
J1/p.u.		6.9322		6.9322
Calculation time/s		0.645		131

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 power between the proposed method and optimization functions is somewhat remarkable. The difference of the optimization results of the voltages, charging power, objective function between the proposed methods and the optimization function, fmincon, of MATLAB for constant impedance load is greater than that for the constant current load which is greater than that for the constant power load. This is because the charging power will cause the voltage drop, resulting in a decline in the power of conventional load. When the conventional household load is constant impedance load, the power of the conventional household load is proportional to the square of the voltage. The power decline is large. When the conventional household load is constant current load, the power of the conventional household load is proportional to the voltage. The power decline is small. When the conventional household load is constant power load, the power of the conventional household load is independent of voltage. The power declined is zero. Therefore, as for the proposed method, the sensitivity of the objective function to the charging power for constant impedance conventional household load is lower than that for constant current load which is lower than that for constant power load. As for constant impedance load, it can be seen that the difference caused by the sensitivity reduction is not much, as shown as Table A1, Table A2 and Table A3 in the Appendix A. The maximum difference of the charging power is less than 8%. The maximum difference of the voltage is less than 0.5%. The maximum difference of the objective function is less than 0.05%. However, the calculation time of the proposed method is much less than that of fmincon of MATLAB.

4.3. Simulation Case 3

4.3.1. Simulation Conditions

As shown in Figure 5, an actual distribution system with 141 nodes is introduced to test the 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:

(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 J₁, J₂, J₃, 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.

(7)

As for objective J₂ and J₃, during the optimization periods, the charging price is 0.77, 0.75, 0.70, 0.70, 0.66, 0.60, 0.48, 0.36, 0.34, 0.30, 0.28, 0.27, 0.27, 0.29 Y/kWh, and β is set to be 100.

(8)

The model of the conventional household load is assumed to be constant impedance load.

(9)

The single phase based power and voltage of the system is 160/3 kVA, $10/\sqrt{3}$ kV, and $0.4/\sqrt{3}$ kV, respectively.

4.3.2. 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 Figure 7, Figure 8 and Figure 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₁, 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 J₂, 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₃, 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₁, 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₂, 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₃, 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%.

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.

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.

4.3.3. Comparison with the Selected Method

The optimization results of the proposed method compared with the method in [15] are shown as

Table 3. Optimization results compared with [15].

Different Methods	Different Models	Objective Function	Minimal Voltage	Calculation Time
The proposed method	$J_1$/p.u.	49.9163	0.9002	19 s
	$J_2$/Yuan	1039.6	0.9000	21 s
	$J_3$/Yuan	31,331.7	0.9000	22 s
The method in [15]	$J_1$/p.u.	49.9266	0.9005	723 s
	$J_2$/Yuan	1040.9	0.9000	746 s
	$J_3$/Yuan	31,331.9	0.9000	757 s

. 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.

5. 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.