Model predictive control for on–off charging of electrical vehicles in smart grids

Massive integration of plug-in electric vehicles (PEVs) into a power grid potentially gives rise to its operating fluctuation, which is not easily controlled due to the unavailability of the PEVs’ power demand profiles prior to their random connections to the grid. The present paper considers the problem of joint coordination of PEV charging and grid power control to minimize both PEV charging cost and energy generation cost in meeting both residential and PEVs’ power demands and suppressing the potential impact of PEV integration. An on-off PEV charging strategy, under which individual PEVs either charge at a maximal power in on-charging-mode or do not charge at all in off-charging-mode at each of a sequence of time slots, is adopted to exploit its simple online implementation. Based on a recently developed model predictive control [1], which is free from assumptions about the probability distribution of PEV arrivals, knowledge of PEVs future demand, or the unlimited charging capacity of PEVs, a mixed integer nonlinear programming problem (MINP) in binary variables of the PEV charging strategy and continuous variables of the grid voltages is proposed at each time for implementing an online algorithm. Due to the large dimension of the continuous voltage variables and the binary on-off decision variables, this MINP is a large scale optimization problem, which is very computationally challenging. A major contribution of the paper is to develop a new solver for this MINP, which is practical for its online computation. Its capability for achieving an optimal solution is shown by numerical examples.


| INTRODUCTION
Significant breakthrough and innovation in battery and vehicle technology have driven an electric-vehicle (EV) boom over the past decade [1] to the extent that EVs are expected to account for 15-30% of new vehicles by 2030 [2]. There is no doubt that serving plug-in EVs (PEVs) is and will continue to be a critical function of tomorrow's smart grids to better leverage renewable energy, reduce power grid operation costs, and lower air pollution emissions [3,4]. However, the growing market penetration of PEVs presents a potential threat to existing power grid systems [5][6][7][8], in that unregulated PEV charging may result in higher peak loads and voltage violations. PEV charging coordination has already proven to be helpful in reducing the cost of power generation and shifting peak loads in the grid [9][10][11][12]. However, only a few of the existing solutions are robust enough to cater to a massive number of PEVs.
Moreover, the anticipated increases in PEV charging are enormous and must be balanced with the power grid's operational requirements.
Recently, PEV charging/discharging coordination has attracted more attention from researchers [13][14][15][16][17]. For example, the authors in Ref. [13] investigated a unidirectional vehicle-to-grid (V2G) network with the aim of maximising aggregator profit given different charging conditions. Reference [14] explored bidirectional power flows with a voltage control on distribution grids for coordinating PEV charging and discharging. Another approach was to flatten the total demand curve with bidirectional energy flow by a decentralised algorithm [15]. Reference [16] also proposed a decentralised charging algorithm. This approach is based on a mean-field game framework and is designed to suit a large number of PEVs, but it does not consider the operational requirements of the power grid, such as the power equation balance, voltage bounds, or line capacity. A new convex optimisation strategy considering battery voltage rise in PEV charging coordination was proposed in Ref. [17].
An on-off charging strategy for PEVs has renewed recent attention due to its simple control structure and efficient online implementation. Under this strategy, PEVs either charge at a fixed power in on-charging mode or do not charge at all in offcharging mode at each time slot [18][19][20]. Generally, the charging time slot varies from half an hour to one hour. Therefore, when PEVs are in the off-charging mode, they can be available for engaging other services. The authors in Ref. [18] studied the PEV on-off charging problem to minimise the overall costs of charging and power generation. By linearising the non-convex constraints of power flow equation, a mixed integer linear programing (MILP) was formulated for this problem. However, the compensation of modelling errors caused by the linearisation method was not analysed. Reference [19] developed a mixed integer non-linear programing (MINP) to address the PEVs' coordination by on-off charging strategy. The first-order Taylor expansion was used to linearise the MINP to a MILP. Such approximation may lead to infeasible results of the original MINP. A similar MILP considering V2G charging strategy was proposed in Ref. [13], where the negative impacts of PEVs integrated to power grid could be potentially suppressed by the bidirectional power source. All coordination methods in Refs. [13][14][15][18][19][20] are off-line mode, where the information of PEVs' arriving time, departing time and initial state of charging (SoC) should be known as a prior. It is not practical to use off-line charging algorithm in these applications. Actually, PEVs are randomly connected to the grid as such that information is hardly known beforehand.
Model predictive control (MPC) is regarded as an effective tool for PEV charging in smart grid. Reference [21] proposed an MPC-based PEV charging model, where the requirements of grid operation were ignored. In addition, PEVs were assumed to be fully charged within a single time slot. However, such an assumption is not possible due to the limit of current battery technology. A MILP for energy storage optimisation over a rolling time horizon was investigated in Ref. [22], which ignored the grid operation constraints as well.
Our previous work in Ref. [23] developed a novel MPC to handle the PEV coordination problem, where the overall costs of charging and power generation was minimised in meeting the requirements of PEV charging and power grid operation. Its distinct practicability is that there are no assumptions on PEVs' charging profiles including arrival/departure time, charging demand and battery SoC. The charging strategy in Ref. [23] is analogue in the sense that at each time slot, PEVs can be charged by any value of power within their battery capacity range. As such, it needs a mechanism to control this charge value, which is not always practical. The present paper adopts the aforementioned on-off charging for PEVs to exploit its easy and efficient online implementation, which also facilitates easy coordination and planning for other activities. However, in contrast to the analogue charging strategy in Ref. [23], which requires online computation of a large-scale non-convex problem on the continuous voltage and PEV charging variables, the on-off strategy requires online computation of a large-scale MINP on the continuous voltage variables and the binary PEV on-off charging decision variables. Thus, the bottleneck for implementing the onoff strategy is the online computation for the large-scale MINP, which is much more computationally difficult than the large-scale NP in Ref. [23]. To the authors' best knowledge, there are even no efficient off-line computations for largescale MINPs, which are the reason that all the previous works (see e.g. [18-20, 22, 24] and references therein) in different contexts must either linearise MINPs at the computation stage or utilise MILPs from the modelling stage to end up with MILPs with the aforementioned compromises. To accomplish the mission of PEV charging coordination with on-off charging strategy, novel techniques are developed to express computationally intractable binary constraints by computationally tractable continuous constraints and to measure the degree of satisfaction of the binary constraints. They are the principal ingredients in developing efficient algorithms for computational solution of this large-scale MINP. The main contributions of this paper are as follows: � A joint coordination solution to cater for massive PEV charging demand and comprehensive power grid control to stabilise fluctuating power demands that reduces total overall system costs. � An on-off strategy for PEV charging that is simple and convenient to implement online, which is crucial from a practical viewpoint. The strategy is based on a new MPC method and modelled with a novel large-scale MINP.
Notably, this method does not require any prior information about the PEVs' arrival and charging information. � A novel MINP solver to compute the joint coordination problem. The MINP model is large-scale, highly challenging computationally, and has no known solution. Hence, the solver is a significant advancement in computation. � Comprehensive simulations and analyses involving realistic charging scenarios with a range of PEV types that demonstrate the computational efficiency and practicality of the proposed methods.
Herein, the formulation of joint coordination of on-off PEV charging and power grid operation based on MINP-MPC is presented in Section 2, where its computational challenges are also analysed. The main technical contribution of the paper is in Section 3, where an efficient MINP solver is proposed. Section 4 considers the computation of a lower bound of this MINP. Numerical results and discussion are provided in Section 5, verifying the ability of the developed solver for seeking an optimal solution. The conclusion is presented in Section 6.
Notation. The imaginary unit is denoted by j, the Hermitian transpose of a vector/matrix A is denoted by A H . A⪰0 denote that A is positive semi-definite. The rank and trace of a matrix A is given by rank(A) and Tr(A), respectively. Rð⋅Þ and Ið⋅Þ, respectively represent the real and imaginary parts of a complex value. |N | denotes the cardinality of set N .

| JOINT COORDINATION OF ON-OFF PEV CHARGING AND POWER GRID OPERATION
A PEV coordination problem is studied, where the overall costs of charging and power generation was minimised in meeting the requirements of PEV charging and power grid operation. This section proposes an MPC for on-off PEV charging in a smart grid system and also analyses its principal computational challenges. This model is designed from the persecutive of the social organiser and distribution network operator. A general structure of smart grid is shown in Figure 1. Consider a power network with N nodes with N ≜ f1; 2; …; Ng denoted as the set of nodes, which are connected by flow lines L ⊆ N � N . Let N ðkÞ denote the set of incident nodes of node k. Suppose node k ∈ G as a generator node and node m ∉ G ⊆ N as a load node. If a node is to serve PEVs, then it is referred to a charging station (CS). The set of CSs is denoted by C.
A total of T time slots are assigned for the serving period that is T ≜ f1; …; T g. Each time slot is half an hour. Priceinelastic load varies according to the profile of residential power demand during the serving period.
In the power network, y km ∈ C is denoted as the admittance of line (k, m), V k (t 0 ) represents the voltage of node k at time slot t 0 while is a Hermitian symmetric matrix variable, whose entries are to replace the voltage operation V k ðt 0 ÞV � m ðt 0 Þ. Denote the active and reactive base demands as P l k ðt 0 Þ and Q l k ðt 0 Þ, respectively; and P g k ðt 0 Þ and Q g k ðt 0 Þ respectively as the real and reactive power generation by node k ∈ G. For t 0 ∈ T , the following constraints about power generation, voltage and phase balance, and line capacity are standard: IðW km ðt 0 ÞÞ ≤ RðW km ðt 0 ÞÞtanðθ max km Þ; ðk; mÞ ∈ L; ð4Þ where P g k , Q g k and � P g k , � Q g k are the limits of the real and reactive power generation, respectively; V k and � V k are the limits of the voltage magnitude, and θ max k;m is given for voltage phase balance, while S km is the upper bound of capacity for line (k, m).
Denote by H k the set of PEVs arriving at charging station k. Let t a;k n and t k n ;d be the arriving and departing time of a PEV k n , it should be fully charged before departing. Suppose C k n be the battery capacity and s 0 k n be the initial state of the battery. The charging power during each time slot of each battery is denoted by � P k n . Unlike Ref. [23], which allows each PEV to charge a power 0 ≤ P k n ðt 0 Þ ≤ � P k n ; the on-off charging strategy is proposed as follows. At each charging slot t 0 , PEV k n charges either with the fixed power (P k n ðt 0 Þ ¼ � P k n ) or zero power (P k n ðt 0 Þ ¼ 0). Figure 2 plots the on-off control for PEV charging in a serving period T ≜ f1; 2; …; 8g. Obviously, a PEV can be available for engaging other services when it is in off-chargingmode.
The following binary variables are introduced for the on-off charging mode. At each time slot t 0 , the charging power of PEV k n is P k n ðt 0 Þ ¼ τ k n ðt 0 Þ � P k n . To make PEV k n fully charged at its departure, the following constraint on binary variables τ k n ðt 0 Þ must be satisfied: -123 Here u h denotes the efficiency of the charging. Given (7) is equivalent to the linear equality constraint on the binary variables τ k n ðt 0 Þ: For any t 0 ∉ t k n ;a ; t k n ;b �, τ k n ðt 0 Þ ¼ 0, and for any k ∉ G, P g k ðt 0 Þ ≡ 0 and Q g k ðt 0 Þ ≡ 0. The supply energy at node k is defined as which is obviously a linear function of W k;m ðt 0 Þ¼ V k ðt 0 ÞV � m ðt 0 Þ though it is seen as a non-linear function of the complex voltage variables V k . For τ k ðt 0 Þ ¼ fτ k n ðt 0 Þg n∈H k , the demand energy at k ∈ C is defined as, leading to the following constraint which is a mixed integer non-linear constraint on the binary variables τ k n ðt 0 Þ and continuous complex voltage variables V k (t 0 ). On the other hand, the demand energy at k ∉ C is obviously defined as leading to the following constraint which is a non-linear constraint on the continuous complex voltage variables V k (t 0 ). By defining the continuous quantities W ≜ fW ðt 0 Þg; and the binary quantities the multi-objective function is defined as where f ðP g k ðt 0 ÞÞ represents power generation cost, β t 0 is the given PEVs' charging price, and � P avg is the averaged power demand over the serving period, which is estimated based on historical data. The first term of (11) is the summation of power generation cost and PEVs' charging cost, while the second term is used to flatten out the total power demand curve in stabilising the grid operations. γ > 0 is the weighting factor to trade-off the two conflicting objectives. The PEVs' charging price β t 0 used herein is the electricity price as PEVs are charged in residential areas. The PEVs' charging cost is customers' cost while the power generation cost is power stations' cost. The power generation cost f ðP g k ðt 0 ÞÞ is usually defined as a linear or quadratic function on the generated reactive power P g k ðt 0 Þ. The parameters of the linear or quadratic do not rely on the electricity price. Over the time horizon [1, T], the joint coordination of PEV charging problem to optimise the total costs of power generation and PEV charging while considering power demand stabilising can be formulated as follows: where the constraint (12c) is needed for legalising the nonlinear variable changes W km (t 0 ) = V k (t 0 )V m (t 0 ). Note that the conventional MPC [25,26] is not applicable as all equations in (12) are not given beforehand. Following Ref. [27], optimisation problem (12) is addressed by proposing the online predictive model at each time slot t 0 as follows.
Let C(t) denote the set of arrived PEVs, which are required to be charged by t. Suppose d k n ðtÞ is the remaining charging demand of PEV k n ∈ C(t) by its departure time t k n ;d . Hence, the following constraints must be imposed: where τ k n ðt 0 Þ ∈ f0; 1g; t 0 ∈ t; t k n ;d �; k n ∈ CðtÞ: ð14Þ Define ΨðtÞ ¼ max k n ∈CðtÞ t k n ;d , the prediction variable W P ðtÞ ≜ fW ðt 0 Þg t 0 ∈t;ΨðtÞ� ; R P ðtÞ ≜ fP g ðt 0 Þ; Q g ðt 0 Þg t 0 ∈t;ΨðtÞ� ; τ P ðtÞ ¼ fτ k n ðt 0 Þg k n ∈CðtÞ;t 0 ∈t;t kn;d � ; and the prediction objective The following MPC is solved over [t, Ψ(t)] at each time t. Only R(t), V(t) and τ(t) are employed to update the solution of (12): rankðW ðt 0 ÞÞ ¼ 1; for t 0 ∈ t; ΨðtÞ�: The difficulty of (16) is focused on the rank-one constraints (16d) and binary constraints (14). To cope with these difficulties, an efficient computational procedure exploiting only the solution of (16) at time slot t for updating online solution of (12) is proposed in the next section.

| TWO-STAGE OPTIMISATION-BASED SOLVER FOR MINP
During the computational procedure, only the solution of MINP (16) at time slot t is required for online updating. Thus, it is not necessary to handle the multiple non-convex matrix rank-one constraints (16d) for all t 0 ∈ [t, Ψ(t)]. In the following, we propose a two-stage optimisation scheme to tackle the computation of (16).
At the first stage, the rank-one constraints in (16d) is dropped and the optimisation (16) is relaxed to, min W P ðtÞ;R P ðtÞ;τ P ðtÞ F P ðR P ðtÞ; τ P ðtÞÞ s:t: ð16aÞ − ð16cÞ: ð17Þ Rðt 0 Þ and b τ k n ðt 0 Þ are the solution of optimisation (16). Otherwise, we go to the second stage and consider the following optimisation by substituting b τ k n ðtÞ into (16b): This non-linear optimisation problem, which involves only one rank-one constraint at t 0 = t in (18c), can be efficiently computed by our previously developed non-smooth optimisation algorithm [27][28][29]. To make the paper self-contained, this optimisation algorithm will be recalled in Section 3.2. The next subsection is devoted to the computation for MINP (17).

| Stage I: computational solution for MINP problem (17)
The key issue for solving MINP problem (17) is to cope with the binary constraint (14) in the optimisation problem (17). Our previous works [30][31][32] have shown that the exactly penalised optimisation, which simultaneously minimises the objective function and maximises the degree of satisfaction of the binary constraints, is appropriate for addressing the MINP SHI ET AL.
-125 (17). The computational efficiency and tractability of the exactly penalised optimisation are critically dependent on the function used to measure the degree of satisfaction of the binary constraints. Therefore, we develop a novel function to measure the degree of satisfaction of the binary constraints. Based on this function, a new path-following computational procedure which iteratively improves solution for the corresponding exactly penalised optimisation problem, is proposed as follows.
First, the equivalence of the binary constraint (14) is established with continuous constraints by the following lemma: The binary constraint (14) can be fulfiled with the following continuous constraints with linear constraint (15),  (20) given τ L k n ðt 0 Þ ¼ τ k n ðt 0 Þ, that is τ k n ðt 0 Þ ∈ f0; 1g, implying (14). Constraint (20) is called reverse convex as g(τ P (t)) is convex in τ P (t) [33]. With the decrease of L, g(τ P (t)) approaches the linear function  Nevertheless, choosing L close to 1 may not work as gðτ P ðtÞÞ − � τðtÞ will approach zero very fast with such setting.
To the authors' best knowledge, Lemma 1 is new. A particular result for L = 2 was obtained in our previous works [30,31].
Thanks to Lemma 1, L = 1.5 is set for Algorithm 1 to accelerate its convergence speed.
With the incorporation of the degree of satisfaction function g 1 into the objective of (17), the following penalised optimisation problem is obtained: min W P ðtÞ;R P ðtÞ;τ P ðtÞ ΦðR P ðtÞ; τ P ðtÞÞ ≜ F P ðR P ðtÞ; τ P ðtÞÞ þμg 1 ðτ P ðtÞÞ s:t: ð15Þ; ð16bÞ; ð16cÞ; ð19Þ; ð22Þ where μ > 0 is a penalty parameter. With a sufficiently large μ, this penalised optimisation problem is exact for (17). The solution of (22) is also an solution for (17) and thus satisfies g 1 (τ P (t)) = 0 [34]. To our best knowledge, using the function g 1 defined by (21) as a measurement to evaluate the satisfaction of bilinear constraint (14) instead of the conventional class � τðtÞ − gðτ P ðtÞÞ is quite new as well.
A path-following algorithm is developed to efficiently solve optimisation (22). First, a lower bound approximation for g (τ P (t)) is derived. Since g(τ P (t)) is a convex function, it is clear that at τ Hence, an approximation of the upper bounding for g 1 (τ P (t)) at the variable τ ðκÞ P ðtÞ can be easily obtained as This iterative procedure is summarised by a pseudo-code in Algorithm 1.

| Stage II: computational procedure for (18)
The main difficulty of optimisation problem (18) concentrates on the non-convex matrix rank-one constraint rank(W(t)) = 1, which can be handled by a non-smooth optimisation method proposed in our previous work [27][28][29]. To make the paper self-contained, this method is recalled here.
TrðW ðtÞÞ − λ max ðW ðtÞÞ ≥ 0; can be obtained directly, where λ max (W(t)) is the maximal eigenvalue of the matrix W(t). Thus, Tr(W(t)) − λ max (W(t)) can be used to measure the satisfaction of the matrix rank-one constraint rankðW ðtÞÞ ¼ 1: We incorporate this term into the objective function and obtain the following exactly penalised optimisation for ν > 0. where w ðκþ1Þ max ðtÞ is the normalized eigenvector of the largest eigenvalue λ max (W (κ+1) (t)) at time slot t. The optimisation (28) can be solved iteratively by the following convex problem: SHI ET AL.

| LOWER BOUNDING BY OFF-LINE COMPUTATION
To investigate the optimal performance of the MPC-based online computation in the previous section, we examine its off-line counterpart in this section, which requires the information of all PEVs including the arriving and departing time, initial SoC of the battery and future charging demand to be known beforehand. Of course, such off-line computation cannot be implemented in practice but it gives a lower bound for the practically implemented online computation.
Similarly, the off-line computation for (12) is of two following optimisation stages.
Stage I. The rank-one constraints in (12c) are dropped as thus (12) is relaxed to the following problem: and Stage II. Let c W , b R and b τ denote the obtained solution of MINP (31). If rankð c W ðtÞÞ ≡ 1 is fulfiled, t ∈ T then b R and b V is accepted as the solution of (12). Otherwise, for those t ∈ T with rankð c W ðtÞÞ > 1 one can substitute b τðtÞ to have (18) and use Algorithm 2 for its computation.

| Setup
Sedumi [35] solver in the framework of CVX [36] is applied to solve the convex optimisation problems (26), (30) and (32) with a Core i7-7600U CPU. The tested grid is a balanced network modified from IEEE 123 test feeder with the nominal voltage of 4.16 kV. Its details and data of system structure, physical limits and cost functions can be found in Ref. [37]. There are three distributed generators, which are respectively connected with node 16, 36, 56. Ten charging stations have been randomly placed on node 4, 7, 11, 17, 21, 28, 33, 40, 44 and 51. The charging period is set from 6 PM to 6 AM, divided into 24 30-min time slots. A truncated normal distribution (8, 1.5 2 ) is adopted to input the PEVs' arrival times independently. There are three types of PEVs: normal PEVs, which are required to be fully charged by 6 AM, median PEVs, which should be fully charged within six hours after their arrivals, and urgent PEVs, which must be charged immediately when connected to the grid. The number of normal PEVs, median PEVs and urgent PEVs for each CS are randomly generated using the following uniform distributions: U [10,15], U [3,6] and U [1,3], respectively. The randomly generated numbers of PEVs served at each CS are given in Table 1. The total number of three types of PEVs are 137, 39 and 23, respectively. According to the charging urgency, the energy price for normal PEVs, median PEVs and urgent PEVs are, respectively, defined as β t , 1.5β t and 2β t , where the energy price β t of on 17-20 May 2017 from 6 PM to 6 AM are given in Figure 3. The power capacity of the PEV battery is C k n ¼ 50 kWh. The initial SoC is set as 20%. The fixed charging power u h � P k n is set as 5 kWh to make that each PEV required eight time slots to be fully charged. Since there are more than 12!/8!4! = 990 feasible on-off charging selection for every PEV, the MINP (12) is not possible to be solved by any exhaustive search strategies. The price-inelastic demand P l k ðtÞ is defined according to Ref. [38] as P l k ðtÞ ¼ lðtÞ where � P l k T with the load demand � P l k indicates the total price-inelastic demand during the serving time period, while lðtÞ ∑ 24 t¼1 lðtÞ with the residential load demand l(t) at t indicates the proportion of price-inelastic demand at each time t. In our 128simulation, the data for � P l k is taken from [37], while the data for l(t) is taken from Ref. [39]. The residential load demand l (t) on 17-20 May 2017 from 6 PM to 6 AM are plotted in Figure 4. The weighting factor γ = 10 3 is set. ϵ = 10 −3 is used as the criteria to stop the proposed algorithm.

| Performance of the algorithms
The numerical results are summarised in Table 2. The effectiveness of online computation based on (18) is confirmed by observing that the cost of the power generation and PEV charging is almost the same to its counterpart computed by the off-line computation based on (12). The average CPU time of Algorithm 1 is provided in the sixth column of Table 2. Furthermore, Figure 5 indicates that Algorithm 1 for all profiles converge rapidly within several iterations. The average CPU time of Algorithm 2 for the four profiles is all within 1 min. It is worth noting that Algorithm 1 is developed to handle the challenging MINP problem (17) with large-scale binary charging variables and continuous voltage variables. To the best of our knowledge, no benchmarking methods can be used to efficiently handle this problem.
The total power demand, the real price-inelastic demand, and the charging demand under Profile one are plotted in Figure 6, while the three types of PEV charging demand under Profile 1 are shown in Figure 7. It can be seen the price-inelastic demand reaches a peak value at 6 PM but then decreases continuously till midnight and remains low values from 0 to 6 AM. The total power demand P tot (t) in real-time, which constitutes of the real price-inelastic demand and total charging demand in real time, is stable during the serving time period. The fluctuation rate of total power demand defined by maxf max t∈T P tot ðtÞ − � P avg � P avg ; � P avg − min t∈T P tot ðtÞ � P avg g is within 7%. One can see that, the magnitude of the total demand are all around 2100 kW, while the charging demand are within 0-700 kW. This means the PEV charging optimality has a substantial impact to the overall optimality. For other profiles, the results are similar. The voltage magnitude during the serving time period under Profile 1 is shown in Figure 8. The voltage magnitude starts to drop after 9 PM since most PEVs charge after that time but their values are always within the range of (1.165, 1.2] pu. Therefore, the negative impact of PEVs' integration to the grid has been successfully suppressed. Figure 9 plots the SoC of two normal PEVs and two median PEVs that arrives at different times randomly selected in Profile 1. The SoC of PEV keeps unchanged for several time slots as they don't charge at those time slots.
The computational performance of the off-line algorithm, which provides a lower bound for the online algorithm. The comparison of the solution between online charging and offline charging process are presented in the last three columns of Table 2. The gap between online charging and off-line charging cost is very light. Figure 10 plots SoC of a normal PEV and a median PEV randomly selected in Profile 1. The charging behaviour is obviously different between online and off-line charging even though their optimal costs are similar.
The performance of charging demand with various number of PEVs during the serving time period is also investigated.  CS #  4  7  11  17  21  28  33  40  44  51   Normal  14  12  14  14  15  12  13  15  13  15   Median  3  4  3  3  3  6  4  6  4  3   Urgent  3  2  2  3  3  2  3  1  1 Figure 11 shows the curves of charging demand corresponding to four different sets of PEVs under Profile 1. The number of PEVs connected to each CS is generated using the following uniform distributions: U [13,28], U [3,12] and U [2,6] for the three types of PEVs, respectively. The proposed algorithms still work well and converge rapidly. Notably, the charging demand climbs rapidly from 6 PM to 0 AM since PEVs are injected to the grid continuously during this period and then it maintains high values during 0 to 6 AM due to low residential demand and energy price during this period. It can also be seen that the total charging demand increases with the integration of more PEVs.

| CONCLUSIONS
The joint coordination solution for on-off PEV charging and controlling grid operations presented here is designed to meet the full charging demands of the massive influx of PEVs expected to enter the market between now and 2030 while maintaining grid operations within safety thresholds. This problem addressed is highly challenging because PEVs connect and disconnect from the grid at random and the binary on-off charging decisions are computationally intractable decision. The benefits of the proposed on-off PEV charging strategy include its efficient online implementation, the lack of a need for a mechanism to ensure analogue PEV charging values remain within battery capacity, and a solution that minimises the total overall costs of PEV charging and power generation.
This paper proposes a new MPC to address this task. At each time slot, the proposed MPC invokes computational solution of a large-scale MINP, which is solved by a novel and easily implemented online algorithm. Three types of PEVs: normal PEVs, median PEVs, and urgent PEVs have been included to stimulate the real charging scenarios. The efficiency and practicality of the system has been verified through comprehensive simulations, which show that the charging behaviour of PEVs and voltage behaviour of the power grid can be controlled very well even with a massive number of PEVs. Further, fluctuations in power demand are stabilised, and the overall cost of generating power and meeting charging demand is reduced. Although the proposed model is designed from the perspective of the social organiser and distribution network operator, it can also be utilised by an EV aggregator by modifying the objective function and keeping the constraints of the current model. Another extension of the proposed method to the problem of optimising demand response in smart grids is under study. -131 CONSTANTS � P avg Averaged power demand over the serving period. � P k n Fixed charging power of PEV n at charging station k ∈ C. β t 0 PEV charging price at time t 0 . C k n Battery capacity of PEV n at charging station k ∈ C. P l k ðt 0 Þ Active price-inelastic demands of node k at time slot t 0 . Q l k ðt 0 Þ Reactive price-inelastic demands of node k at time slot t 0 .

S km
Power capacity of line (k, m). u h Charging efficiency of PEV battery. y km Admittance of line (k, m). � P g k Upper bound of active power generated by node k ∈ G. � Q g k Upper bound of reactive power generated by node k ∈ G. � V k Upper bound of the amplitude of voltage injection at node k ∈ N . P g k Lower bound of active power generated by node k ∈ G. Q g k Lower bound of reactive power generated by node k ∈ G.

VARIABLES τ k n
Binary variable used to represent the on-off charging state of PEV n at charging station k ∈ C. P g k ðt 0 Þ Active power generation by generator node k ∈ G at time t 0 . Q g k ðt 0 Þ Reactive power generation by generator node k ∈ G at time t 0 . V k (t 0 ) Complex voltage variable of node k ∈ N at time t 0 . W(t 0 ) Semidefinite matrix variable introduced to transfer the nonconvex power flow constraints with voltage product.