Aggregative Game for Distributed Charging Strategy of PEVs in a Smart Charging Station

Abstract

This paper proposes a charging strategy for plug-in electric vehicles (PEVs) in a smart charging station (SCS) that considers load constraints and time anxieties. Due to the rapidly growing load demand of PEVs and the load capacity investments in infrastructure, PEV charging needs to be subject to overload limits, beyond which failures can occur. The time anxiety is presented to address some of the uncertainties that may arise while charging PEVs. Under an aggregative game framework, this paper constructs a price-driven charging model to minimize costs by choosing the optimal charging strategy. Meanwhile, since the driver information is an aggregated item in the PEV cost function, the drivers’ privacy can be protected. Then, a distributed reflected forward–backward (RFB) splitting method is developed to search for the generalized Nash equilibria (GNE) of the game. The convergence of the proposed algorithm and the effectiveness of the charging strategy are verified by the detailed simulation and results.

1. Introduction

Governments worldwide are promoting plug-in electric vehicles (PEVs) as a clean alternative to conventional gasoline vehicles due to the depletion of fossil fuel resources and environmental pollution, both of which are serious issues. Compared with the conventional gasoline vehicle, the PEV reduces carbon dioxide emissions and the overall operating cost [1]. However, with the rapid growth of PEVs, the total charging load of PEVs entering the electric grid has increased [2], which leads to an overload of charging stations. According to [3], uncoordinated PEV charging lowers the electrical grid’s power quality. However, the demand presented in [4] might actually flatten as a result of the coordinated charging of PEVs. To increase the electric grid’s operational effectiveness and security, we must formulate effective charging strategies to coordinate and control the charging behavior of PEV drivers.

Since PEVs are capable of storing electrical energy, some PEVs can be used as energy suppliers and transmit power in both directions, to and from a smart charging station (SCS), which can mitigate the effects of overcharging at SCS at peak times. Meanwhile, in order to improve the accuracy of the charging strategy, the driver’s time anxiety needs to be taken into account, in addition to the plug-in electric vehicle factor. The authors of Vatanparvar et al. [5] introduced the concept of demand–price elasticity to address the demand–response model based on the drivers’ behavior. A robust-index method was developed in [6] to resolve driver behavior uncertainties and minimize violations to comfort in household load scheduling. In [7], a distributed method was proposed and a means of reducing the impact of some uncertain events under a non-cooperative game was proposed.However, some of these driver behavior models were homogeneous and lacked theoretical justification. In order to improve the validity of the model, four different PEV driver behaviors were proposed to study the effect of time anxiety on plug-in electric vehicle charging. There are two main classes of control architecture in this charging strategy, namely the centralized and distributed approaches.

The authors of Yan et al. [8] proposed a four-stage optimal control method for electric vehicle charging stations (EVCSs) to reduce operating costs and maintain the supply–demand power balance. In [9], a centralized control approach based on multistage droop control was used to operate an island microgrid in the presence of high PEV penetration. Furthermore, the model in [10] utilized a centralized real-time charging scheme based on a convex relaxation method to coordinate PEV charging. However, in the centralized strategy, each driver, as well as the control center, needs to know the complete information, including cost function, local feasible decision sets and affine sharing constraints, which can easily compromise the drivers’ privacy. Meanwhile, in many cases, no central node can communicate with all drivers in both directions, and it is not possible to realize large-scale driver communication in a real-time manner. Therefore, a distributed strategy is used to compute the local decisions corresponding to PEVs and generalized Nash equilibria and communicate with the adjacent local decisions using the drivers’ local data, which reduces the computational and interaction burden and maintains the privacy of driver information. Based on a decentralized protocol, Li et al. [11] presented a charge control technique for large-scale PEVs under the Newton-type algorithm. The authors of Gan et al. [12] studied a decentralized algorithm to optimize PEV charging during off-peak hours. However, the distributed algorithm described above rarely considers the strategic interaction between multiple PEVs in the charging station.

Therefore, to study communication between multiple PEVs, we then focused on game theory, as this is a powerful tool for analyzing the interactions between multiple decision makers and can improve the model’s performance. Since the famous Cournot model was proposed, aggregative games have become an important type of game theory. In the description for the aggregative game in [13], each driver was not subject to a one-to-one interaction, but was subject to a number of aggregations across the charging strategy. Recent studies in [14,15,16,17] considered the linear aggregation functions and quadratic cost functions in such games, as well as their interaction with plug-in electric vehicle charging. Based on the literature, a distributed reflected forward–backward splitting method was proposed in [18] to find the generalized Nash equilibria of the aggregative game. Unlike the distributed algorithm presented above, this algorithm exploits the aggregated coupling structure in the cost function, meaning that each agent only needs to exchange and maintain an aggregated estimate, not including the estimate of the even multiplier, reducing the computational and communication burden of the model. In addition, since each PEV only shares its estimate of the total, it does not need to share the complete information, further protecting the privacy of the drivers. The main contributions of this paper are summarized as follows:

(1)

A new price-driven charging model combining time anxiety and load constraints is constructed to minimize the cost of an individual PEV driver within the framework of an aggregative game. In particular, as everyone only knows the final summation result, not the specific information, the aggregation game can better protect the privacy of drivers.

(2)

Load constraints are proposed to protect the safety of SCS. Then, four PEV driver behaviors are proposed based on different time anxiety states and load constraints. Meanwhile, the effects of time anxiety under four different driver behaviors are compared, and the effects of uncertain occurrence events are reduced by the charging strategy.

(3)

A distributed reflected forward–backward algorithm is designed to seek the generalized Nash equilibria of the model. The proposed algorithm seeks its optimal response charging strategy regarding the current load and time anxiety in the electric grid, thus preventing overload in the smart charging station and mitigating the impact of uncertain events that may occur at the PEV charging time. The algorithm obtains an improvement in significant convergence compared to the numerical values of the FB algorithm [19].

2. System Model and Problem Formulation

As this paper focuses on charging PEVs, the studied system is named SCS, which is considered to be part of the electric grid. In [20,21,22], SCS is composed of PEVs, which are powered by the electric grid, and the PEV can transmit power in both directions between the electric grid. Meanwhile, each PEV i ($i\in \mathcal{I}=\left\{1,2,\cdots ,N\right\}$) is the electric grid user in Figure 1. We define $d_t$ as the basic demand load of the system at time t. Let vector $x_i=col(x_{i,1},\dots ,x_{i,T})$ denote the charging curve of PEV i over all the time slots and $x_{i,t}$ define the charging power of PEV i at time t. In the SCS, the driver can query the electrical grid constraint described in Figure 1 to prevent the required cost from exceeding the load. Moreover, as different drivers have different needs, they may have a preference for a particular charging time, as they encounter various forms of uncertainty, i.e., events that may occur at a later charging time; therefore, they may feel an urge to finish charging the PEV earlier. Therefore, each PEV in SCS can be charged flexibly based on the available information, which can meet the demand and reduce the charging cost.

2.1. Feasible Charging Coordination Constraint Profiles

2.1.1. Battery Capacity Constraint for PEV i

We describe the battery dynamics of PEV i using a linear model (1), in which ${\pi }_i^t$ represents the state of charge of PEV i at time t. Meanwhile, ${\pi }_i^{min}$ denotes the lower limit of battery capacity and ${\pi }_i^{max}$ denotes the upper limit of battery capacity for PEV i:

$${\pi }_i^{t+1}={\pi }_i^t+x_{i,t},{\pi }_i^{min}\le {\pi }_i^t\le {\pi }_i^{max}$$

(1)

2.1.2. Charging Constraint for PEVs

For each PEV i, the total electrical energy that it obtains needs to meet its charging requirement needs. Based on this constraint, let $R_i$ be the required electrical energy. Therefore, the following equation holds:

$$\sum _{t=1}^T{x}_{i,t}=R_i$$

(2)

Moreover, each PEV i can choose to charge and discharge according to its situation, which depends on the PEV battery capacity. Note here that PEVs cannot be selected for charging or discharging at the same time; only one of them can be selected. $\underline{x_i}$ is defined as the minimum discharging power of a PEV i, and $\overline{x_i}$ is defined as the maximum charging power of a PEV i:

$$\underline{x_i}\le x_{i,t}\le \overline{x_i}$$

(3)

2.1.3. Overload Constraint for Charging PEVs

To avoid overloading the SCS, the charging demand of all the PEVs cannot exceed the maximum electrical energy supplied by the SCS. Let $C_{max}$ denote the total PEV load, and the basic demand load at moment t is set to $d_t$, which denotes the basic electrical load transmitted through the SCS at time t. As the constraint spatially couples the charging demand of all PEVs using the SCS at time t, constraint (4) is referred to as a joint constraint:

$$\sum _{i=1}^N{x}_{i,t}+d_t\le C_{max}$$

(4)

2.1.4. Feasible Charging Profiles

We assume that the set ${\mathcal{X}}_i$ is the feasible charging profile and ${\mathcal{X}}_i$ is a nonempty, compact convex set. Let $\mathcal{X}=\left\{\left(x_1,x_2,\dots ,x_N\right)\mid x_i\in {\mathcal{X}}_i,\forall i\in \mathcal{I}\right\}$ denote the set of charging profiles of all PEVs. Then, $\mathcal{X}$ is also nonempty, compact and convex. Therefore, set ${\mathcal{X}}_i$ can be written as follows:

$${\mathcal{X}}_i=\left\{x_i\mid \left(1\right)\left(2\right)\left(3\right)\left(4\right)\right\}$$

(5)

2.2. Cost Function of PEVs

We assume that there are N drivers, in which r ($r\le N$) start at the charging time and enter a state of anxiety after a period of time. For PEV i, $t_i^a$ is denoted as the arrival time to the SCS and $t_i^d$ as the departure time at which SCS will be left; the discharge is no longer considered after anxiety begins. Considering the aggregative game framework, we defined $\sigma \left(x\right)=\frac{1}{N}{\sum }_{i=1}^N{x}_i$, and the minimization cost function of PEV i under the corresponding energy consumption constraint is defined as follows:

$$\underset{x_i}{min}{f}_i\left(x_i,\sigma \left(x\right)\right)=\sum _{t=1}^T{S}_{i,t}{\int }_0^{x_{i,t}}{\rho }_t\left(d_t+{\mu }_{i,t}+\sum _{j\ne i}{x}_{j,t}\right)d{\mu }_{i,t}$$

(6)

in which $S_{i,t}$ denotes the driver’s anxiety influence at time t. Next, the pricing function ${\rho }_t\left(L_t\right)$ is defined

$$\begin{array}{cc}\hfill & {\rho }_t\left(L_t\right)={\alpha }_t{L}_t+{\beta }_t\hfill \\ \hfill & \underset{x_i}{min}{f}_i\left(x_i,\sigma \left(x\right)\right)=\sum _{t=1}^T{S}_{i,t}{\int }_0^{x^{i,t}}\left({\alpha }_t\left(d_t+{\mu }_{i,t}+\sum _{j\ne i}{x}_{j,t}\right)+{\beta }_t\right)d{\mu }_{i,t}\hfill \end{array}$$

(7)

where $L_t=d_t+{\mu }_{i,t}+{\sum }_{j\ne i}{x}_{j,t}$ represents the total load at time t. The pricing function is for a nonlinear form that varies with the increasing charging load of PEV i, while the models in [23,24] use ${\rho }_t$ as a linear function. As a consequence, the linearly decreasing marginal benefits [12,23,24] correspond to the PEV quadratic pricing function. In the meantime, the PEV i increases a load of ${\mu }_{i,t}$ to the electricity grid at time t; the driver of PEV pays an amount of ${\rho }_t\left(d_t+{\mu }_{i,t}+{\sum }_{j\ne i}{x}_{j,t}\right)\mathsf{\Delta }{\mu }_{i,t}$ for charging at the time t. When $\mathsf{\Delta }{\mu }_{i,t}$ tends toward zero and ${\mu }_{i,t}$ changes from zero to $x_{i,t}$, it then evolves into Formula (8):

$$\begin{array}{c}\hfill \underset{x_i}{min}{f}_i\left(x_i,\sigma \left(x\right)\right)={\left(d+N\sigma \left(x\right)-\frac{1}{2}{x}_i\right)}^{\mathrm{T}}{S}_i\alpha x_i+1_T^{\mathrm{T}}{S}_i\beta x_i\end{array}$$

(8)

2.3. Time Anxiety for Drivers

The design of the cost function (8) shows that a higher price leads to lower charging power. Furthermore, if the driver has time anxiety about charging and discharging the PEV, the charging power needs to be higher. Based on this conclusion, $S_{i,t}$ must be small to have a high charging power. Following that, at the time of charging, the value of $S_{i,t}$ gradually rises. Meanwhile, the local objective function was defined according to the problem of anxious and non-anxious drivers, so we obtained two different settings.

The first setting is non-anxious drivers, and we defined S as a constant. This means that the driver has a constant anxiety factor at the charging time, which does not change over time. In this model, we denoted $AT\left(i\right)$ as the arrival time of PEV i to SCS and denoted $DT\left(i\right)$ as the departure time of PEV i from SCS:

$$S_{i,t}=\left\{\begin{array}{cc}S,\hfill & \mathrm{A}T\left(i\right)\le t<DT\left(i\right)\\ 0,\hfill & \mathrm{otherwise}\end{array}\right.$$

(9)

and the next setting is anxious drivers. We defined $V_{i,t}$ as the ratio of the duration of the anxiety state to the overall anxiety time. Thus, the following equation holds:

$$V_{i,t}=\frac{t-t_i^a}{t_i^d-t_i^a},t_i^a\le t<t_i^d$$

(10)

With the time anxiety studied in this paper, at the time from $AT\left(i\right)$ to $t_i^{a,}$ the time anxiety is at its lowest value $S_{i,min}$, and from $t_i^d$ to $DT\left(i\right)$, it is at its highest value $S_{i,max}$. This paper will explore the behavior of PEV drivers based on [7] and propose four different behaviors:

(1)

Non-time-anxious driver (NTAD): This type of driver reaches an anxious time directly after entering a state of peak anxiety (Figure 2).

(2)

Less time-anxious driver (LTAD): This type of driver has anxiety values that rise quickly and then slowly after entering the anxious time. The rise is faster and then slower (Figure 3).

(3)

Mid-time-anxious driver (MTAD): This type of driver enters anxious time with anxiety values increasing at a uniform rate (Figure 4).

(4)

High-time-anxious driver (HTAD): This type of driver has anxiety values that rise slowly and then quickly after entering anxious time. The rise is slow and then fast (Figure 5).

Thus, we can obtain the following equation:

$$S_{i,t}=\left\{\begin{array}{cc}{S}_{i,max},\hfill & \mathrm{for}\mathrm{NTAD}\hfill \\ ln\left[V_{i,t}(e-1)+1\right]\times \left(S_{i,max}-S_{i,min}\right)+S_{i,min},\hfill & \mathrm{for}\mathrm{LTAD}\hfill \\ V_{i,t}\left(S_{i,max}-S_{i,min}\right)+S_{i,min},\hfill & \mathrm{for}\mathrm{MTAD}\hfill \\ \frac{e^{V_{i,t}}-1}{e-1}\left(S_{i,max}-S_{i,min}\right)+S_{i,min},\hfill & \mathrm{for}\mathrm{HTAD}\hfill \end{array}\right.$$

(11)

In Figure 2, Figure 3, Figure 4 and Figure 5, the value of $S_{i,max}-S_{i,min}$ indicates the PEV driver’s time anxiety, i.e., the depth of time anxiety. This indicates that the smaller the $S_{i,min}$, the greater the time anxiety of the PEV driver. Furthermore, a greater time anxiety will increase the willingness of the PEV driver to meet the charging demand earlier and closer to the departure time, thus reducing the amount of time anxiety.

According to the preceding discussions, two factors will influence PEV drivers’ time anxiety. The intensity of time anxiety is the first, and the curve’s shape is the second. To represent the effect of curve shape, the curve shapes of the four different behavior types regarding time anxiety are shown in Figure 2, Figure 3, Figure 4 and Figure 5. We chose NATD as the reference, i.e., we chose $S_{i,max}$ as the reference value. The difference between this and the impact of time anxiety caused by the suggested behavior is referred to as the impact difference. For example, $\mathsf{\Delta }{S}_{i,t}^{\mathrm{LTAD}}=S_{i,t}^{NATD}-S_{i,t}^{\mathrm{LTAD}}$ is defined as the impact difference of the LTAD behavior. Therefore, the following inequality (12) holds in Figure 2, Figure 3, Figure 4 and Figure 5:

$$\mathsf{\Delta }{S}_{i,t}^{\mathrm{NTAD}}\le \mathsf{\Delta }{S}_{i,t}^{\mathrm{LTAD}}\le \mathsf{\Delta }{S}_{i,t}^{\mathrm{MTAD}}\le \mathsf{\Delta }{S}_{i,t}^{\mathrm{HTAD}}$$

(12)

3. Distributed Charging Strategy

In this section, we propose a distributed reflected forward–backward splitting method to find a GNE of the function (8). Then, for the optimal response of each PEV driver, a distributed charging strategy is proposed. We assume that all PEVs in the SCS are selfish and each PEV is only allowed information about local problem data. Centralized control methods are typically unavailable in this situation.

3.1. Game Model

As the PEV charging and discharging problem is a generalized Nash equilibrium [24], we consider a group of agents $\mathcal{I}=\{1,\dots ,N\}$ that seek a GNE of the aggregative game with globally shared affine constraints and the gradient condition of KKT necessary optimality conditions can be then given by (13):

$$\begin{array}{cc}\hfill {\nabla }_{x_i}{f}_i\left(x_i,\sigma \left(x\right)\right)=& \left(d+N\sigma \left(x\right)+x_i\right)S_i\alpha +1_T{S}_i\beta \hfill \end{array}$$

(13)

Assumption 1. 

For each $i\in \mathcal{I}$, the function in (12) is differentiable and convex, and $\mathsf{\Omega }\subset R^n$ is a closed convex set.

We define $R_m$ as m-dimensional Euclidean space and each agent i chooses its local decision $x_i\in {\mathsf{\Omega }}_i\subset R^n$. We call $\mathrm{x}=col\left(x_1,\cdots ,x_N\right)=col{\left(x_i\right)}_{i\in \mathcal{I}}\in \mathsf{\Omega }\subset R^{Nn}$ the decision profile, i.e., the stacked vector ${\prod }_{i=1}^N{\mathsf{\Omega }}_i$$=\mathsf{\Omega }$. The aim of each agent i is to optimize its objective function, $f_i\left(x_i,\sigma \left(x\right)\right):\mathsf{\Omega }\to R$, within its feasible decision set. Note that $f_i\left(x_i,\sigma \left(x\right)\right)$ is nonlinearly coupled to the decisions of the other agents, but may not be explicitly coupled to the decisions of all other agents. We denote

$$\mathcal{X}:=\left\{x_i\in {\mathsf{\Omega }}_i\mid \sum _{i=1}^N{A}_i{x}_i\le \sum _{i=1}^N{b}_i\right\}$$

(14)

where $A_i\in R^{m\times n}$, $b_i\in R^m$ are local data from agent i. Thus, we can obtain the following formula:

$$\underset{x_i}{min}{f}_i\left(x_i,\sigma \left(x\right)\right)s.t.\sum _{i=1}^N{A}_i{x}_i\le \sum _{i=1}^N{b}_i$$

(15)

Furthermore, by considering games with affine sharing constraints $Ax\le b$, and supposing $x^{\ast }$ as a GNE of game (12), the optimal solution to the following convex optimization problem is defined as

$$\underset{x_i}{min}{f}_i\left(x_i^{\ast },\sigma \left(x^{\ast }\right)\right)s.t.x_i\in {\mathsf{\Omega }}_i,A_i{x}_i\le b-\sum _{j\ne i,j\in \mathcal{I}}{A}_j{x}_j^{\ast }$$

(16)

where $b={\sum }_{i=1}^N{b}_i\in R^m$. The set ${\mathsf{\Omega }}_i$ denotes the local decision set of agent i and the matrix $A_i$ defines how agent i is involved in the coupling constraint. However, the constraints of Equation (5) are clearly different from those of Equation (16). Therefore, the matrix $A_i$ and vector $b_i$ are then divided into two submatrices, $H_i$ and $W_i$, and subvectors, $P_i$ and $Q_i$. We define these as follows:

$$W_i=\left[\begin{array}{c}{0}_{(i-1)T\times T}\\ -E\\ 0_{(N-i)T\times T}\\ 0_{(i-1)T\times T}\\ E\\ 0_{(N-i)T\times T}\\ I\end{array}\right],P_i=\left[\begin{array}{c}{0}_{(i-1)}\\ R_i\\ 0_{(N-i)}\end{array}\right]$$

(17)

$$Q_i=\left[\begin{array}{c}{0}_{(i-1)T}\\ \left({\pi }_i^1-{\pi }_i^{min}\right)1_T\\ 0_{(i-1)T}\\ 0_{(i-1)T}\\ \left({\pi }_i^{max}-{\pi }_i^1\right)1_T\\ 0_{(i-1)T}\\ \frac{C_{max}}{N}{1}_T\end{array}\right],H_i=\left[\begin{array}{c}{0}_{(i-1)\times T}\\ 1_{1\times T}\\ 0_{(N-i)\times T}\end{array}\right]$$

(18)

where $0_T$ represents the zero vector of T dimension, I represents the unit matrix and $E=\left[\begin{array}{ccccc}1& & & & \\ 1& 1& & & & \\ ·& ·& ·& & & \\ ·& ·& & ·& & \\ ·& ·& & & ·& \\ 1& 1& ·& ·& ·& 1\end{array}\right]\in R^{T\times T}$. By definition, the constraints ${\sum }_{i=1}^N{W}_i{x}_i\le {\sum }_{i=1}^N{Q}_i,{\sum }_{i=1}^N{H}_i{x}_i={\sum }_{i=1}^N{P}_i$ are satisfied. Then, the constraints in Equation (5) are satisfied. Note that ${\sum }_{i=1}^N{H}_i{x}_i={\sum }_{i=1}^N{P}_i$ is an equation constraint, while Equation (16) contains only the inequality constraint. Therefore, we represent this Equation constraint using two inequality constraints, satisfying both ${\sum }_{i=1}^N{H}_i{x}_i\le {\sum }_{i=1}^N{P}_i$ and $-{\sum }_{i=1}^N{H}_i{x}_i\le -{\sum }_{i=1}^N{P}_i$. If using this approach, we need to rewrite $A_i=\left[\begin{array}{c}{W}_i\\ H_i\\ -H_i\end{array}\right]$ and $b_i=\left[\begin{array}{c}{Q}_i\\ P_i\\ -P_i\end{array}\right]$, and the constraints in Equation (5) are converted into the constraints in Equation (16).

Assumption. 2 

For each $i\in \mathcal{I}$, and for each $\xi \in \mathrm{E}$, the function $f_i(x_i,\sigma \left(x\right),\xi )$ is Lipschitz continuous, convex, and continuously differentiable. For $\sigma \left(x\right)$, the Lipschitz constant $\ell \left(\sigma \right(x),\xi )$ is integrable in ξ.

Among all possible generalized Nash equilibria, we are concerned with those solution sets that correspond to the set of solutions to an appropriate variational inequality. For this purpose, let us define the (pseudo) gradient mapping as

$$\mathrm{F}\left(x\right)=col\left({\left(\mathrm{E}\left[{\nabla }_{x_i}{f}_i\left(x_i,\sigma \left(x\right),{\xi }_i\right)\right]\right)}_{i\in \mathcal{I}}\right)$$

(19)

We define a local Lagrangian function for agent i as $\mathcal{L}(x,z,\lambda )=〈\mathrm{F}\left(x^{\ast },J{x}^{\ast }\right),x〉+{\ell }_{\mathsf{\Omega }}\left(x\right)+{\lambda }^{\mathrm{T}}(Ax-b)+{\ell }_{R_{+}^{Nm}}\left(\lambda \right)+z^{\mathrm{T}}{L}_{\lambda }\lambda$, where ${\lambda }_i\in R_{+}^m$ is a dual variable associated with the coupling constraint. When $x^{\ast }$ is an optimal solution to (16), the following Karush–Kuhn–Tucker (KKT) conditions are satisfied:

$$\forall i\in \mathcal{I}:\left\{\begin{array}{c}0\in \mathrm{F}\left(x^{\ast },J{x}^{\ast }\right)+{\mathrm{N}}_{\mathsf{\Omega }}\left(x^{\ast }\right)+A^{\mathrm{T}}{\lambda }^{\ast }\hfill \\ 0=L_{\lambda }{\lambda }^{\ast }\hfill \\ 0\in A{x}^{\ast }-b-{\mathrm{N}}_{R_{+}^{Nm}}\left({\lambda }^{\ast }\right)+L_{\lambda }{z}^{\ast }\hfill \end{array}\right.$$

(20)

In order to ensure that all the preceding signs are + and facilitate the operation, the third formula becomes $0\in -A{x}^{\ast }+b+L_{\lambda }{\lambda }^{\ast }+{\mathrm{N}}_{R_{+}^{Nm}}\left({\lambda }^{\ast }\right)-L_{\lambda }{z}^{\ast }$. Since $\widehat{u}=Mu$ is orthogonal to $Ju$, there is no consistent vector in the space of $\widehat{u}$ to make $L_u\widehat{u}=0$ when, and only when, $\widehat{u}=0$. Meanwhile, we introduce $L_u$ to implement a distributed estimation. If we use I instead of $L_u$, we need the central node to pass the average information. Based on this result, the extended KKT condition is as follows:

$$\left\{\begin{array}{c}0\in \mathrm{F}\left(x^{\ast },J^{\ast }+{\widehat{u}}^{\ast }\right)+{\mathrm{N}}_{\mathsf{\Omega }}\left(x^{\ast }\right)+A^{\mathrm{T}}{\lambda }^{\ast }\hfill \\ 0=c{L}_u{\widehat{u}}^{\ast }\hfill \\ 0=L_{\lambda }{\lambda }^{\ast }\hfill \\ 0\in -A{x}^{\ast }+b+L_{\lambda }{\lambda }^{\ast }+{\mathrm{N}}_{R_{+}^{Nm}}\left({\lambda }^{\ast }\right)-L_{\lambda }{z}^{\ast }\hfill \end{array}\right.$$

(21)

3.2. Distributed Algorithm

We assume that each driver only knows their local data, i.e., $f_i(x_i,\sigma \left(x\right))$, ${\mathsf{\Omega }}_i$, $A_i$ and $b_i$, which contains their own private information. Meanwhile, the shared affine coupling constraints are decomposed such that each driver knows only one local block of the constraint matrix. Note that $A_i$ describes how agent i participates in the coupling constraints (shared global resources), which is also assumed to be privately known by driver i. The globally shared constraint $Ax\le b$ then couples the set of feasible decisions of the agents, but is not known to any agent. Next, we describe the preprocessing process leading to the distributed iteration proposed in the algorithm. In this, $x_{i,k}$, $u_{i,k}$, $z_{i,k}$ and ${\lambda }_{i,k}$ are the state variables of agent i at iteration k and ${\tau }_i$, $v_i$, ${\alpha }_i$ are fixed constant step-sizes for driver i. We define the weighted adjacency matrix $W={\left[w_{i,j}\right]}_{i,j}\in R^{N\times N}$. The set of neighbors of PEV i is ${\mathcal{N}}_i^{\lambda }=\left\{j\mid w_{i,j}>0\right\}$. We define the operators as follows:

$$\begin{array}{c}\overline{\mathrm{A}}:\left(\begin{array}{c}x\hfill \\ u^{\perp }\hfill \\ z\hfill \\ \lambda \hfill \end{array}\right)\mapsto \left(\begin{array}{c}\mathrm{F}\left(x,Jx+u^{\perp }\right)\hfill \\ c{L}_u{u}^{\perp }\hfill \\ 0\hfill \\ b\hfill \end{array}\right)+\left(\begin{array}{c}0\\ 0\\ 0\\ L_{\lambda }\lambda \end{array}\right)\end{array}$$

(22)

$$\begin{array}{c}\overline{\mathrm{B}}:\left(\begin{array}{c}x\hfill \\ u^{\perp }\hfill \\ z\hfill \\ \lambda \hfill \end{array}\right)\mapsto \left(\begin{array}{c}{\mathrm{N}}_{\mathsf{\Omega }}\left(x\right)\hfill \\ 0\hfill \\ 0\hfill \\ N_{R_{+}^{\mathrm{Nm}}}\left(\lambda \right)\hfill \end{array}\right)+\left[\begin{array}{cccc}0& 0& 0& A^{\mathrm{T}}\\ 0& 0& 0& 0\\ 0& 0& 0& L_{\lambda }\\ -A& 0& -L_{\lambda }& 0\end{array}\right]\left(\begin{array}{c}x\hfill \\ u^{\perp }\hfill \\ z\hfill \\ \lambda \hfill \end{array}\right)\end{array}$$

(23)

in which $c\in R_{+}$. Meanwhile, the metric matrix $\mathsf{\Phi }$ is defined as follows. Note that the matrix $\mathsf{\Phi }$ is symmetric and positive definite, and the first term of $\mathsf{\Phi }$ is the antisymmetric matrix in the operator $\overline{\mathrm{B}}$:

$$\begin{array}{cc}\hfill \mathsf{\Phi }& =\left[\begin{array}{cccc}0& 0& 0& -A^{\mathrm{T}}\\ 0& 0& 0& 0\\ 0& 0& 0& -L_{\lambda }\\ -A& 0& -L_{\lambda }& 0\end{array}\right]+\left[\begin{array}{cccc}{\tau }^{-1}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& v^{-1}& 0\\ 0& 0& 0& {\alpha }^{-1}\end{array}\right]+\left[\begin{array}{cccc}-S_C& -S& 0& 0\\ -S& {\kappa }^{-1}I& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{cccc}{\tau }^{-1}-S_C& -S& 0& -A^{\mathrm{T}}\\ -S& {\kappa }^{-1}I& 0& 0\\ 0& 0& v^{-1}& -L_{\lambda }\\ -A& 0& -L_{\lambda }& {\alpha }^{-1}\end{array}\right]\hfill \end{array}$$

(24)

Through the RFB algorithm, we can obtain the following formula $\mathsf{\Phi }{\overline{v}}_{k+1}+\overline{\mathrm{B}}{\overline{v}}_{k+1}=\mathsf{\Phi }{\overline{v}}_k-\overline{\mathrm{A}}{v}_k$ and $v_{k+1}=2{\overline{v}}_{k+1}-{\overline{v}}_k$, which is in the form of a distributed reflected forward–backward splitting method to find zeros of (${\mathsf{\Phi }}^{-1}\mathrm{A}+{\mathsf{\Phi }}^{-1}\mathrm{B}$). We substitute the operators into $\mathsf{\Phi }{\overline{v}}_{k+1}+\overline{\mathrm{B}}{\overline{v}}_{k+1}=\mathsf{\Phi }{\overline{v}}_k-\overline{\mathrm{A}}{v}_k$, and then, through calculation, we can obtain the following equation:

$$\left\{\begin{array}{c}{\mathrm{N}}_{\mathsf{\Omega }}\left({\overline{x}}_{k+1}\right)+{\overline{x}}_{k+1}={\overline{x}}_k+\tau \left(-A^{\mathrm{T}}{\overline{\lambda }}_k-\mathrm{F}\left(x_k,J{x}_k+u_k^{\perp }\right)-c{L}_k{u}_k\right)\hfill \\ {\overline{u}}_{k+1}^{\perp }={\overline{u}}_k^{\perp }-\kappa c{L}_k{u}_k^{\perp }+M\left({\overline{x}}_{k+1}-{\overline{x}}_k\right)\hfill \\ {\overline{z}}_{k+1}={\overline{z}}_k-v{L}_{\lambda }{\overline{\lambda }}_k\hfill \\ {\mathrm{N}}_{R_{+}^{Nm}}\left({\overline{\lambda }}_{k+1}\right)+{\overline{\lambda }}_{k+1}={\overline{\lambda }}_k+\alpha \left(A\left(2{\overline{x}}_{k+1}-{\overline{x}}_k\right)+L_{\lambda }\left(2{\overline{z}}_{k+1}-{\overline{z}}_k\right)-b-L_{\lambda }{\lambda }_k\right)\hfill \\ x_{k+1}=2{\overline{x}}_{k+1}-{\overline{x}}_k\hfill \\ u_{k+1}^{\perp }=2{\overline{u}}_{k+1}^{\perp }-{\overline{u}}_k^{\perp }\hfill \\ z_{k+1}=2{\overline{z}}_{k+1}-{\overline{z}}_k\hfill \\ {\lambda }_{k+1}=2{\overline{\lambda }}_{k+1}-{\overline{\lambda }}_k\hfill \end{array}\right.$$

(25)

its initial condition is ${\overline{u}}_0={\overline{x}}_0$ and $u_0=x_0$.

Based on the above conclusions, the algorithm can find the GNE of the game, i.e., the strategy that finds the cost minimization for PEV i. Its convergence was proved in [18]. Meanwhile, by writing the above algorithm in distributed form, we can obtain the following Algorithm 1:

Algorithm 1: Distributed charging strategy with reflected forward-backward algorithm.

1: Initialization:$S_{i,min}=S_{i,max}=S$, $x_i\in {\mathsf{\Omega }}_i$, $u_i\in R^n$, ${\lambda }_i\in R^m$ and $z_i\in R^m$ 2: Task: slove (16) 3: For$k=1:k_{max}$ 4:         For $i=1:N$ 5: (1)Receives$x_{i,k}$ for $j\in {\mathcal{N}}_i^f$, $u_{j,k}$, $z_{j,k}$ and ${\lambda }_{j,k}$ for $j\in {\mathcal{N}}_i^{\lambda }$ then updates 6: ${\overline{x}}_{i,k+1}={proj}_{{\mathsf{\Omega }}_i}\left({\overline{x}}_{i,k}+{\tau }_i\left(-A_i^{\mathrm{T}}{\overline{\lambda }}_{i,k}-{\nabla }_{x_i}{f}_i\left(x_{i,k},u_{i,k}\right)-c{\sum }_{j=1}^N{\omega }_{ij}\left(u_{i,k}-u_{j,k}\right)\right)\right)$ 7: ${\overline{u}}_{i,k+1}={\overline{u}}_{i,k}-c\kappa {\sum }_{j=1}^N{\omega }_{ij}\left(u_{i,k}-u_{j,k}\right)+{\overline{x}}_{i,k+1}-{\overline{x}}_{i,k}$ 8: ${\overline{z}}_{i,k+1}={\overline{z}}_{i,k}-v_i{\sum }_{j=1}^N{\omega }_{ij}\left({\overline{\lambda }}_{i,k}-{\overline{\lambda }}_{j,k}\right)$ 9: ${\overline{\lambda }}_{i,k+1}={proj}_{R_{+}^m}({\overline{\lambda }}_{i,k}+{\overline{\alpha }}_i{A}_i(2{\overline{x}}_{i,k+1}-{\overline{x}}_{i,k}-{\overline{\alpha }}_i\left(b_i-{\sum }_{j=1}^N{\omega }_{ij}({\lambda }_{i,k}-{\lambda }_{j,k})\right)$ 10: $+{\overline{\alpha }}_i{\sum }_{j=1}^N{\omega }_{ij}\left(2\left({\overline{z}}_{i,k+1}-{\overline{z}}_{j,k+1}\right)-\left({\overline{z}}_{i,k}-{\overline{z}}_{j,k}\right)\right))$ 11: (2)Receives${\overline{x}}_{i,k}$, ${\overline{u}}_{i,k}$, ${\overline{z}}_{i,k}$, ${\overline{\lambda }}_{i,k}$then updates 12: $x_{i,k+1}=2{\overline{x}}_{i,k+1}-{\overline{x}}_{i,k}$ 13: $u_{i,k+1}=2{\overline{u}}_{i,k+1}-{\overline{u}}_{i,k}$ 14: $z_{i,k+1}=2{\overline{z}}_{i,k+1}-{\overline{z}}_{i,k}$ 15: ${\lambda }_{i,k+1}=2{\overline{\lambda }}_{i,k+1}-{\overline{\lambda }}_{i,k}$ 16:         end 17: end 18: if${\sum }_{t=t_i^a}^{t_i^d}{x}^{i,t}>p_i^c$then 19:        $S_{i,min}=S_{i,min}-\epsilon$ 20:        Go back to step 3 21: else 22:        break 23: end if

Step 3 can be explained as follows. We define a step-size $\epsilon$ and set a threshold value $p_i^c$ for time-anxious drivers. Since the initialization $S_{i,min}=S_{i,max}=S$, i.e., all drivers are not time-anxious, the algorithm is executed once, and a Nash equilibrium is solved. Then, we determine whether the sum of charging capacity ${\sum }_{t=t_i^a}^{t_i^d}{x}_{i,t}$ for PEV i at the time $t\in [t_i^a,t_i^d]$ is greater than the threshold value $p_i^c$. If ${\sum }_{t=t_i^a}^{t_i^d}{\overline{x}}_{i,t}>p_i^c$, let $S_{i,min}=S_{i,min}-\epsilon$. Then, the new $S_{i,min}$ can be substituted to solve the Nash equilibrium once more. Otherwise, it is straightforward to derive the Nash equilibrium solution.

4. Simulation and Numerical Results

In this section, the performance of the proposed algorithm with load constraints and time anxiety is evaluated by minimizing the charging cost for PEV i in the SCS. For further illustration, we consider the SCS with 10 PEVs in the residential area. We also investigated the charging power distribution under non-time-anxious and time-anxious conditions. The simulation configuration was set up as follows.

4.1. Overload Control for 10 PEVs

In this scenario, depending on the owner’s preferences and needs, 10 PEVs arrive at the SCS. The charging needs of each PEV were chosen between 50 and 60 KW, and the charging power was chosen between 10 and 15 KW/h.

Table 1. Constraint-related parameters.

PEV i	${\mathit{\pi }}_{\mathit{i}}^1$	${\mathit{R}}_{\mathit{i}}$	${\mathit{\pi }}_{\mathit{i}}^{min}$	${\mathit{\pi }}_{\mathit{i}}^{max}$	$\overline{{\mathit{x}}_{\mathit{i}}}$	$\underline{{\mathit{x}}_{\mathit{i}}}$	$\mathit{AT}$	$\mathit{DT}$
1	$7.5$	53	5	75	15	$-15$	12	17
2	6	56	5	80	10	$-10$	22	14
3	8	$52.5$	5	75	10	$-10$	3	12
4	$5.6$	$51.4$	5	70	15	$-15$	2	10
5	$6.7$	51	5	65	10	$-10$	15	21
6	$6.5$	56	5	75	15	$-15$	2	8
7	$6.1$	$52.4$	5	65	10	$-10$	4	20
8	9	51	5	80	10	$-10$	13	24
9	$9.2$	$50.8$	5	70	15	$-15$	10	22
10	$7.5$	53	5	75	15	$-15$	18	24

lists the PEV parameters, arrival time (AT), and departure time (DT). Usually, the daily peak charging demand occurs from 12:00 to 17:00, when people go to the SCS to charge their PEVs. Therefore, we define ${\alpha }_t=0.3$$/kWh during peak hours (i.e., from 18:00 to 6:00), ${\alpha }_t=0.2$$/kWh during off-peak hours (i.e., other times), and the initial electricity price is ${\beta }_t=0.3$$/kWh. We assume that the maximum power supplied by the SCS at time t is 130 KW, which is obtained from the value of [25].

As shown in Figure 6, the red dotted line shows the maximum capacity to support PEV charging at time t, i.e., $L_{max}=$ 130 KW. The green and purple solid lines denote the base load and the charging requirement load, respectively, for the 24 h. As shown in Figure 6, the SCS is severely overloaded at peak hours due to uncoordinated PEV charging, which may damage the SCS. Therefore, we present the overload control constraint in the framework of the game (16). After the overload control constraint, the load profile is shown in Figure 7. Compared with Figure 6, as the SCS in Figure 7 holds the charging capacity fixed at peak hours, the charging strategy shifts the excess charging capacity to free time, i.e., off-peak hours. The results show that the charging requirement load is always below the value of $L_{max}$. Therefore, the strategy ensures the safety of the SCS.

4.2. Time Anxiety for PEVs

The simulation provides a further charging strategy considering load constraints. According to the settings in Table 1, the charging and discharging powers of each PEV are shown in Figure 8.

In these four pictures, the blue line indicates a charging strategy that does not consider time anxiety, while the red-brown line indicates a charging strategy that does consider time anxiety (note that the two are only different charging strategies; the amount that the PEV is charged does not change but is simply shifted). Figure 9 shows the NATD driver’s behavior in our simulation. From the discussion in Section 2.3, it follows that the driver has no time anxiety; therefore, the blue and red-brown lines overlap, and we use one of the blue lines to represent the PEV.

Figure 10 shows the LATD driver’s behavior in our simulation. From the previous discussion in Section 2.3, it follows that this driver will have some time anxiety, i.e., there may be a delay in the charging of the electric vehicle due to something that occurs during the time anxiety; therefore, a charging strategy that considers time anxiety will move some of the charging within the time interval $[t_i^a,t_i^d]$ of time anxiety in the non-time anxious time interval, effectively avoiding the situation of missed EV charging due to unpredictable circumstances. This is represented in the graph by the transfer of the charge from anxious energy to shifted energy, which can be seen in Figure 10 as a reduction in the charge in the time anxious interval. The PEV driver is satisfied with the current charging method, as it is considered robust.

We simulated the behavior of the MATD driver, and, as can be seen in Figure 11, it shifted more of its charge to the rest of the time interval during the anxiety time interval than in Figure 10. The HATD driver, as can be seen in Figure 12, shifted more of its charge to the rest of the time interval during the anxious period than in Figure 11. This is because PEV drivers of various anxiety levels have a predefined threshold, which is lower if the drivers want to be more robust in response to uncertain events (e.g., HATD has the strongest anxiety). If the current anxiety energy exceeds the threshold, the PEV driver is dissatisfied with the current charge level, which would be insufficient to meet his or her charging needs in an uncertain event. Our algorithm is implemented in an iterative manner until the PEV driver succeeds in bringing his/her anxiety energy below its threshold, as indicated by the red-brown line. Our algorithm is iteratively implemented until the anxiety energy of the electric vehicle driver falls below a particular threshold.

The simulation of these four different driver behaviors leads to the same conclusion as discussed in Section 2.3: the greater the driver’s time anxiety, the greater that driver’s willingness to meet the charging demand earlier, i.e., more charging is transferred within the interval $[t_i^a,t_i^d]$ of time anxiety.

4.3. Convergence Analysis

In this simulation, we provide an iterative process performed by all the PEVs at 18:00. PEVs do not need a central node to be able to bidirectionally communicate with all other PEVs. To preserve privacy, each PEV computes the corresponding decision in a distributed manner using its cost function, feasible set and coupling constraints. At each iteration k, the PEVs update their decisions and their estimates. The iterative process $x_k$ denotes the decision variable with PEV i; $u_k$ denotes the aggregated estimate, which includes the parameters that affect the electricity price; ${\lambda }_k$ is used to ensure that the constraints hold; $z_k$ is used as an auxiliary variable to ensure that ${\lambda }_k$ is consistent. As shown in Figure 13, the PEVs all converge to their optimal charging strategies during the iterative process. As a result, the proposed distributed RFB algorithm under the aggregative game can solve the charging problem in SCS proficiently.

In addition, we simulated the convergence accuracy of the RFB algorithm and the FB algorithm. Compared to the FB algorithm, the RFB algorithm converged to ${10}^{-7.1}$ and the FB algorithm converged to ${10}^{-3.6}$ after 9000 iterations under the same conditions shown in Figure 14. The effectiveness of the algorithm is demonstrated by the fact that it only takes 3000 iterations to converge to ${10}^{-3.6}$ using the RFB algorithm. The faster convergence indicates that all PEVs are solved to arrive at the optimal strategy faster.

5. Conclusions

In this paper, a new distributed charging and discharging strategy for PEVs, based on time anxiety and load constraints, was proposed in the framework of the aggregative game. Time anxiety was proposed to mitigate the effects of some uncertain events that may occur during charging. The load constraint was proposed to make the PEV charging more coordinated and protect the safety of the SCS. Detailed case studies were presented, showing that the charging strategy with time anxiety and load constraint considered is more reasonable and reduces the total cost. The distributed reflected forward–backward algorithm was designed to seek the generalized Nash equilibria of the game model. The proposed algorithm achieved optimal driver response in a theoretically fast, distributed manner and protected the driver’s privacy. The effectiveness of the proposed algorithm was verified by example simulations. In the future, we will focus on the aggregative game for charging stations.