Robust Integration of Electric Vehicles Charging Load in Smart Grid’s Capacity Expansion Planning

Abstract

Battery charging of electric vehicles (EVs) needs to be properly coordinated by electricity producers to maintain the network reliability. In this paper, we propose a robust approach to model the interaction between a large fleet of EV users and utilities in a long-term generation expansion planning problem. In doing so, we employ a robust multi-period adjustable generation expansion planning problem, called R-ETEM, in which demand responses of EV users are uncertain. Then, we employ a linear quadratic game to simulate the average charging behavior of the EV users. The two models are coupled through a dynamic price signal broadcasted by the utility. Mean field game theory is used to solve the linear quadratic game model. Finally, we develop a new coupling algorithm between R-ETEM and the linear quadratic game with the purpose of adjusting in R-ETEM the uncertainty level of EV demand responses. The performance of our approach is evaluated on a realistic case study that represents the energy system of the Swiss “Arc Lémanique” region. Results show that a robust behaviorally-consistent generation expansion plan can potentially reduce the total actual cost of the system by 6.2% compared to a behaviorally inconsistent expansion plan.

Appendices

Appendices

Formulation of the ETEM Energy Model

In this appendix, we present the formulation of the deterministic ETEM energy model. We mostly use a similar notation as in Aliakbarisani et al. [1]. However, we slightly depart from Aliakbarisani et al. [1] where we use a regional demand response variable $\varvec{V}$. Let us start with the definition of the sets. Set ${\mathbb {T}}$ indicates time periods, and ${\mathbb {S}}$ is the set of time slices inside a period. Set ${\mathbb {C}}$ is for energy commodities, and ${\mathbb {P}}$ identifies technologies. Input-output energy flows are shown in set ${\mathbb {F}}$. Finally, set ${\mathbb {L}}$ identifies buses in different geographical zones. A full nomenclature is presented in Table 3.

Table 3 Nomenclature used in ETEM formulation

$$\begin{aligned}&\min \sum _{t,l} \rho ^t \cdot \left( \sum _{p} \alpha _{t,p}\varvec{C}_{t,l,p} + \pi _{t,p} ( \sum _{k =0}^{l_p - 1}\varvec{C}_{t-k,l,p} + \varOmega _{t,l,p}) +\sum _{s,c} \nu _{t,p}\varvec{P}_{t,s,l,p,c} \right. \nonumber \\&\qquad \left. +\sum _{s,c} \left( \lambda _{t,s,c} \varvec{I}_{t,s,l,c} - \lambda ^{\prime }_{t,s,c} \varvec{E}_{t,s,l,c}\right) + \sum _{s,l^{\prime },c} \lambda ^{''}_{t,s,l,l^{\prime },c} \varvec{T}_{t,s,l,l^{\prime },c} \right) , \end{aligned}$$

(32a)

$$\begin{aligned}&\text {s.t.}\nonumber \\&\left( \sum _{p \in {\mathbb {P}}^P_c} \varvec{P}_{t,s,l,p,c} + \varvec{I}_{t,s,l,c}\right) \eta _c + \sum _{l^{\prime } \ne l} \left( \eta _c \varvec{T}_{t,s,l^{\prime },l,c} - \varvec{T}_{t,s,l,l^{\prime },c}\right) \nonumber \\&\quad \ge \sum _{p \in {\mathbb {P}}^C_c} \varvec{P}_{t,s,l,p,c} + \varvec{E}_{t,s,l,c} \quad \forall \; t,s,l,c \in {\mathbb {C}} / {\mathbb {C}}^{{\mathcal {D}}} \end{aligned}$$

(32b)

$$\begin{aligned}&\Big (\sum _{l,p \in {\mathbb {P}}^{c}}\theta _{p}^c \cdot \beta _{t,s,p} \left( \sum _{k =0}^{l_p-1}\varvec{C}_{t-k,l,p} + \varOmega _{t,l,p}\right) + \sum _{l,p \in {\mathbb {P}}^P_c/ {\mathbb {P}}_c} \theta _{p}^c \cdot \varvec{P}_{t,s,l,p,c} + \varvec{I}_{t,s,l,c} \Big )\cdot \varrho _{t,s,c} \nonumber \\&\quad \ge \sum _{l,p \in {\mathbb {P}}^C_c} \varvec{P}_{t,s,l,p,c} + \varvec{E}_{t,s,l,c} \quad \forall t,s \in {\mathbb {S}}^{{\mathcal {G}}},c \in {\mathbb {C}}^{{\mathcal {G}}} \end{aligned}$$

(32c)

$$\begin{aligned}&\eta _{c} \sum _{p \in {\mathbb {P}}^{P_{c}}} \varvec{P}_{t,s,l,p,c} = \sum _{p \in {\mathbb {P}}^{C_{c}},s^{\prime } \in {\mathbb {S}}^s} \varvec{P}_{t,s^{\prime },l,p,c} \quad \forall t,s,l,c \in \mathbb {C}\mathbb {S} \end{aligned}$$

(32d)

$$\begin{aligned}&\sum _{p \in {\mathbb {P}}^P_c} \varvec{P}_{t,s,l,p,c} \ge \varTheta _{t,l,c} \varvec{V}_{t,s,l,c} \quad \forall t,s,l,c \in {\mathbb {C}}^{{\mathcal {D}}} \end{aligned}$$

(32e)

$$\begin{aligned}&\upsilon _{t, s,l, c} (1-\nu _{t, s,l, c}) \le \varvec{V}_{t, s,l, c} \le \upsilon _{t, s,l, c} (1+\nu _{t, s,l, c}) \quad \forall t,s,l,c \in {\mathbb {C}}^{{\mathcal {D}}} \end{aligned}$$

(32f)

$$\begin{aligned}&\sum _{s \in {\mathbb {S}}_j} \varvec{V}_{t, s,l, c} = \sum _{s \in {\mathbb {S}}_j} \upsilon _{t, s,l, c} \quad \forall t,l,c \in {\mathbb {C}}^{{\mathcal {D}}},j \in {\mathbb {J}} \end{aligned}$$

(32g)

$$\begin{aligned}&\sum _{c:p\in {\mathbb {P}}^P_c} \varvec{P}_{t,s,l,p,c} \le \beta _{t,s,p} \left( \sum _{k=0}^{l_p-1}\varvec{C}_{t-k,l,p} + \varOmega _{t,l,p}\right) \quad \forall t,s,l,p \notin {\mathbb {P}}^R, \end{aligned}$$

(32h)

$$\begin{aligned}&\sum _{c \in {\mathbb {C}}^p_m} \varvec{P}_{t,s,l,p,c} = \beta _{t,s,p} \left( \sum _{k=0}^{l_p-1}\varvec{C}_{t-k,l,p} + \varOmega _{t,l,p}\right) \quad \forall t,s,l,p \in {\mathbb {P}}^R \end{aligned}$$

(32i)

$$\begin{aligned}&\sum _{c \in {\mathbb {C}}_{f^{'}}} \varvec{P}_{t,s,l,p,c} = \eta _{f,f^{'}}^t \cdot \sum _{c \in {\mathbb {C}}_f } \varvec{P}_{t,s,l,p,c} \quad \forall f \in \mathbb {F}\mathbb {I}_p, f^{'}\in \mathbb {F}\mathbb {O}_p, t,s,l,p \end{aligned}$$

(32j)

$$\begin{aligned}&(\varvec{P},\varvec{I},\varvec{E},\varvec{T}) \in {\mathcal {Y}} \; , \; \varvec{C} \in {\mathcal {X}} \; , \; \varvec{P}\ge 0,\varvec{I}\ge 0,\varvec{E}\ge 0,\varvec{T}\ge 0,\varvec{C}\ge 0,\varvec{V}\ge 0\, . \end{aligned}$$

(32k)

Equation (32a) is the objective function which minimizes a discounted sum of all costs of the system over all regions $(l \in {\mathbb {L}})$ and time periods $(t \in {\mathbb {T}})$. The total cost includes investment, fixed and operational costs of technologies, import, export and transmission costs of energy commodities. Variables $\varvec{C},\varvec{P},\varvec{I},\varvec{E},\varvec{T}$ are capacity installation, energy production, import, export and regional transmission of energy, respectively. Parameters $\alpha _{t,p}$, $\nu _{t,p}$, $\lambda _{t,s,c}$, $\lambda ^{\prime }_{t,s,c}$, $\lambda ^{''}_{t,s,l,l',c}$, and $\pi _{t,p}$ are unit costs of the associated variables.

A commodity balance constraint is presented in (32b). It assures that during each period t and time-slice s, the regional procurement of energy commodity c is more or equal to the overall consumption. Total regional procurement, the left-hand side of the constraint, includes i) total production of commodity c in region l by all technologies producing it (${\mathbb {P}}^P_c$), ii) import of commodity c, and iii) net transmission of commodity c into region l. It is worth mentioning that while import and export refers to the transfer of energy from outside of the boundary of the energy system, energy transmission is the amount of energy which is produced in one region and transmitted to another region inside the energy system. Input energy c is multiplied by parameter $\eta _c$, which is the efficiency of the network with respect to commodity c, e.g., the efficiency of electricity transmission lines. On the right-hand side, the overall consumption of commodity c is equal to internal consumption by technologies ${\mathbb {P}}^C_c$, added to the amount of commodity that is exported.

In addition to the energy balance constraint (32b), for commodities in set ${\mathbb {C}}^{{\mathcal {G}}}$ constraint (32c) provides a safety margin, during peak time-slices $s \in {\mathbb {S}}^{{\mathcal {G}}}$, to protect against random events not explicitly represented in the model. Parameter $\varrho _{t,s,c} \in [0,\;1]$ represents the fraction of reserved capacity needed to ensure covering the peak load. The left-hand side of this constraint models the maximum amount of commodity $c \in {\mathbb {C}}^{{\mathcal {G}}}$ that can be procured in period t and time-slice s. This amount is equal to the sum of i) the maximum production capacity of commodity c by technologies that produce it as their main output (${\mathbb {P}}_c$), ii) the production of commodity c by technologies that produce c as their by-product of their main activity, and iii) the import of commodity c. The left-hand side is the total consumption similar to constraint (32b). Parameter $\theta _{p}^c$ is the proportion of technology production that can be used during the peak period. Constraint (32d) is a balance constraint for storage commodities. Namely, depending on the efficiency of the storage technology i.e., $\eta _c$, a fixed portion of the amount of storage at time-slice s can be restored at subsequent time-slice $s^{\prime }$. Constraints (32e) - (32g) together are modeling the notion of demand response. Specifically, constraint (32e) ensures that for each commodity c in the set of useful demand ${\mathbb {C}}^{{\mathcal {D}}}$, the total production accounted for by variable $\varvec{P}$ must satisfy the demand in each time slice and region. Parameter $\varTheta _{t,l,c}$ is the total demand and variable $\varvec{V}_{t,s,l,c}$ attempts to optimally distribute the total demand of period t into all time-slices $s \in {\mathbb {S}}$ inside period t and region l. Constraint (32f) limits the demand response to vary within a certain margin $\nu _{t,s,l,c}$ from the nominal value, i.e., $\upsilon _{t,s,l,c}$. In addition, since the shift of the demand is only possible between time slices in a day, constraint (32g) limits the sum of the demand response to be equal to the sum of the nominal values in each season.

Constraint (32h), known as capacity factor constraint, limits the maximum activity of each technology to the available capacity of that technology. Parameter $\beta _{t,s,p}$ is the capacity factor and simply represents the fraction of the total capacity which is available at each time-slice. In addition, constraint (32i) enforces that the production by renewable technologies take their maximum possible value considering the total available capacity of these technologies. In ETEM, the total installed capacity of a technology p in period t equals to the aggregated newly installed capacity in the last $l_p$ periods (the life duration of technology p) plus the already available capacity of that technology from before the beginning of the planning horizon. Mathematically speaking, the total available capacity of technology p is $\sum _{k=0}^{l_p-1}\varvec{C}_{t-k,l,p} + \varOmega _{t,l,p}$ (see the right-hand side of equations (32h), or (32i)), where $l_p$ is the life duration of technology p, $\varvec{C}_{t,l,p}$ is the newly installation of technology p in region l and period t, and $\varOmega _{t,l,p}$ is the already available capacity of technology p in period t from before the beginning of the planning horizon. In other words, the installed capacity of technology p in period t always remains in the system until $l_p$ years, and it retires (will not appear in the summation of newly installed technology, i.e., $\sum _{k=0}^{l_p-1}\varvec{C}_{t-k,l,p}$) after $l_p$ years. Constraint (32j) defines the efficiency of technology p. Parameter $\eta _{f,f{'}}^t$ is the efficiency of technology p with output flow $f^{'} \in \mathbb {F}\mathbb {O}_p$ and input flow $f \in \mathbb {F}\mathbb {I}_p$.

Finally, space ${\mathcal {X}}$ and ${\mathcal {Y}}$ represent operational, technical, and economical constraints that define the structure of the energy network and a desirable space for capacity and procurement decisions. For example, constraints on $\hbox {CO}_2$ emissions, technology market penetration, technology ramping, total energy import and export are among the constraints that form up the space ${\mathcal {X}}$ and ${\mathcal {Y}}$. Since these constraints do not affect our analysis in this paper, we omit to report them and refer the reader to Babonneau et al. [4] for a complete list of constraints.

Differences Between R-ETEM and Babonneau et al. [5]

In Babonneau et al. [5], the authors propose to robustify the deterministic version of ETEM presented in (1) by replacing constraint (3) with the robust constraint:

$$\begin{aligned} \varDelta _{l,c}^{{ \text{ min }}}\le \sum _{t\in {\mathbb {T}}}\sum _{s\in {\mathbb {S}}} \varTheta _{t,l,c} \varvec{V}_{t,s,l, c} \le \varDelta _{l,c}^{{ \text{ max }}}\,,\,\forall \,l\in {\mathbb {L}},c\in {\mathbb {C}}^{{\mathcal {D}}} \end{aligned}$$

(33)

with

$$\begin{aligned} \varDelta _{l,c}^{{ \text{ min }}}:=\sup _{\xi \in {\hat{\varXi }}} \sum _{(t,s)\in {\mathbb {T}}\times {\mathbb {S}}}\varDelta ^{l,c,-}_t + ({\bar{\varDelta }}^{l,c}_t-\varDelta ^{l,c,-}_t)\xi _t \end{aligned}$$

and

$$\begin{aligned} \varDelta _{l,c}^{{ \text{ max }}}:=\inf _{\xi \in {\hat{\varXi }}} \sum _{(t,s)\in {\mathbb {T}}\times {\mathbb {S}}}\varDelta ^{l,c,+}_t + ({\bar{\varDelta }}^{l,c}_t-\varDelta ^{l,c,+}_t)\xi _t \end{aligned}$$

where ${\hat{\varXi }}$ is a budgeted uncertainty set in the non-negative orthant, while $[\varDelta ^{l,c,-}_t,\,\varDelta ^{l,c,+}_t]$ is a confidence region for the demand response at time t, and ${\bar{\varDelta }}^{l,c}_t$ is the average value. We first note that (33) is a relaxation (therefore an optimistic approximation) of:

$$\begin{aligned}&\varDelta ^{l,c,-}_t + ({\bar{\varDelta }}^{l,c}_t-\varDelta ^{l,c,-}_t)\xi _t \le \varTheta _{t,l,c} \varvec{V}_{t,s,l, c}\,,\,\forall \,\xi \in {\hat{\varXi }},\,\forall \,(s,t,l,c)\in {\mathbb {S}}\times {\mathbb {T}}\times {\mathbb {L}}\times {\mathbb {C}}^{{\mathcal {D}}} \end{aligned}$$

(34a)

$$\begin{aligned}&\varTheta _{t,l,c} \varvec{V}_{t,s,l, c}\le \varDelta ^{l,c,+}_t + ({\bar{\varDelta }}^{l,c}_t-\varDelta ^{l,c,+}_t)\xi _t \,,\,\forall \,\xi \in {\hat{\varXi }},\,\forall \,(s,t,l,c)\in {\mathbb {S}}\times {\mathbb {T}}\times {\mathbb {L}}\times {\mathbb {C}}^{{\mathcal {D}}}\,, \end{aligned}$$

(34b)

which imposes a robust set of upper and lower bonds on the demand response plan in each time-slice. Thus, we can interpret any feasible solution of constraint (34) (and approximately (33)) as identifying a demand response that has a robust potential of being plausible with respect to the distribution that generated the statistics captured in $(\varDelta ^{l,c,-}_t,{\bar{\varDelta }}_t,\varDelta ^{l,c,+}_t)$. In particular, the constraint offers no protection against the demand response deviations captured by $(\varDelta ^{l,c,-}_t,{\bar{\varDelta }}_t,\varDelta ^{l,c,+}_t)$. This is especially noticeable when considering that $\varDelta _{l,c}^{{ \text{ min }}}$ and $\varDelta _{l,c}^{{ \text{ max }}}$ are both non-increasing and non-decreasing with respect to $\varDelta ^{l,c,-}_t$ and $\varDelta ^{l,c,+}_t$, respectively. This implies that if the observed demand response becomes more uncertain, constraint (33) actually becomes less restrictive for $\varTheta _{t,l,c}$.

Numerical Algorithm to Obtain the Average Optimal State of Charge and Demand

To compute the representative EV user’s average optimal state of charge and average demand, one needs to numerically solve FBODEs (14–15) to obtain $(\eta _\tau , w_\tau )$. Then, the solution, $(\varPsi _\tau , \varPi _\tau )_{\tau \in [0,T]}$ is obtained by solving the FBODEs (21–22). Finally, we the average optimal state of charge and average demand is calculated with the feedback formula (28) and (27). This procedure is done with a simple implementation of the Euler Scheme. Algorithm 2 presents the detail of this implementation.

Babonneau et al. [5] Returns the Solution of BI-DR Model (1)

In the following, we identify weak conditions under which the approach proposed in Babonneau et al. [5] converges in one iteration simply recommending the solution of the behaviorally inconsistent demand response model (1).

Assumption 1

The aggregate mean demand profile ${\mathbb {E}}[u_\tau ^*]$ observed at equilibrium in the mean field game integrates to the total discharge captured by the discharge profile, i.e.,

$$\begin{aligned} \int _0^T \varDelta _{\tau } \mathrm{d}\tau = \int _0^T v_\tau \mathrm{d}\tau . \end{aligned}$$

(35)

This assumption is for instance satisfied when the mean field game captures the stationary behavior of agents over a repeated period of length T. Indeed, in this context the optimal average state of charge at the beginning and end of the horizon must be equal, i.e., ${\mathbb {E}}[X_0^*]={\mathbb {E}}[X_T^*]$. Specifically,

$$\begin{aligned} 0&={\mathbb {E}}[X_T^*]-{\mathbb {E}}[X_0^*]={\mathbb {E}}\left[ x_0+\int _0^T(u_\tau ^*-v_\tau )\mathrm{d}\tau + \varepsilon W_T\right] -{\mathbb {E}}[X_0^*]\\&=\int _0^T{\mathbb {E}}[u_\tau ^*] \mathrm{d}\tau -\int _0^T v_\tau \mathrm{d}\tau \Rightarrow \int _0^T\varDelta _\tau \mathrm{d}\tau = \int _0^T v_\tau \mathrm{d}\tau \end{aligned}$$

Proposition 2

Let Assumption 1 be satisfied and the discharge profile of all linear quadratic game models be consistent with the total demand of the season that it describes in ETEM, namely:

$$\begin{aligned} \sum _{s\in {\mathcal {S}}}\varTheta _{t,l,c}\upsilon _{t,s,l,c}=\int _0^T v_\tau ^{l,c,t} \mathrm{d}\tau \,,\, \forall t \end{aligned}$$

(36)

Then, the solution ${\overline{V}}$ of the non-robust ETEM (1) satisfies constraint (33) formulated using the equilibrium of the mean field games.

Proof

We will first show that the lower bounding part of constraint (33) is satisfied. Namely,

$$\begin{aligned} \varDelta ^{{ \text{ min }}}_{l,c}&= \sup _{\xi \in {\hat{\varXi }}} \sum _{(t,s)\in {\mathbb {T}}\times {\mathbb {S}}}\varDelta ^{l,c,-}_{t,s} + ({\hat{\varDelta }}^{l,c}_{t,s}-\varDelta ^{l,c,-}_{t,s})\xi _{t,s} \end{aligned}$$

(37a)

$$\begin{aligned}&= \sup _{\xi \in {\hat{\varXi }}} \sum _{(t,s)\in {\mathbb {T}}\times {\mathbb {S}}} {\hat{\varDelta }}^{l,c}_{t,s} - ({\hat{\varDelta }}^{l,c}_{t,s}-\varDelta ^{l,c,-}_{t,s})(1-\xi _{t,s}) \end{aligned}$$

(37b)

$$\begin{aligned}&= \sum _{(t,s)\in {\mathbb {T}}\times {\mathbb {S}}} {\hat{\varDelta }}^{l,c}_{t,s} - \inf _{\xi \in {\hat{\varXi }}} \sum _{(t,s)\in {\mathbb {T}}\times {\mathbb {S}}} ({\hat{\varDelta }}^{l,c}_{t,s}-\varDelta ^{l,c,-}_{t,s})(1-\xi _{t,s}) \end{aligned}$$

(37c)

$$\begin{aligned}&\le \sum _{(t,s)\in {\mathbb {T}}\times {\mathbb {S}}} {\hat{\varDelta }}^{l,c}_{t,s} \end{aligned}$$

(37d)

$$\begin{aligned}&=\sum _{t\in {\mathbb {T}}}\left( \int _0^T\varDelta _\tau ^{l,c,t} \mathrm{d}\tau \right) \end{aligned}$$

(37e)

$$\begin{aligned}&= \sum _{t\in {\mathbb {T}}}\left( \int _0^T v_\tau ^{l,c,t} \mathrm{d}\tau \right) \end{aligned}$$

(37f)

$$\begin{aligned}&= \sum _{t\in {\mathbb {T}}}\left( \sum _{s\in {\mathbb {S}}}\varTheta _{t,l,c}\upsilon _{t,s,l,c}\right) \end{aligned}$$

(37g)

$$\begin{aligned}&=\sum _{t\in {\mathbb {T}}}\left( \sum _{s\in {\mathbb {S}}} \varTheta _{t,l,c}{\overline{V}}_{t,s,l,c}\right) \end{aligned}$$

(37h)

where ${\hat{\varXi }}:=\{\xi \in [0,\,1]^{{\mathbb {T}}\times {\mathbb {S}}}|\sum _{(t,s)\in {\mathbb {T}}\times {\mathbb {S}}}\xi _{t,s}\le \varGamma \}$, and which would confirm that ${\overline{V}}$ satisfies the lower bound. In detail, the inequality in (37d) follows from the infimum being greater or equal to zero due to ${\hat{\varXi }}\subseteq [0,\,1]^{{\mathbb {T}}\times {\mathbb {S}}}$ and ${\hat{\varDelta }}^{l,c}_{t,s}\ge \varDelta ^{l,c,-}_{t,s}$. The equality in (37d) follows from $\int _0^T\varDelta _\tau ^{l,c,t} \mathrm{d}\tau =\sum _{s\in {\mathbb {S}}} {\hat{\varDelta }}^{l,c}_{t,s}$ based on the definition of ${\hat{\varDelta }}^{l,c}_{t,s}$. Equalities (37f) and (37g) follow from Assumption 1 and (36), respectively. Finally, (37g) follows from the fact that ${\overline{V}}$ satisfies constraint (4).

A similar analysis can be done to prove that ${\overline{V}}$ also satisfies the upper bound in constraint (33). $\square$

The following corollary follows directly from Proposition 2 given the nature of the procedure described in Babonneau et al. [5], which consists in 1) solving model (1); 2) communicating the price information to the LQG model; 3) identifying the equilibrium using the limiting MFG; 4) adding (or updating) constraint (33) to model (1) based on MFG solution; 5) iterating until convergence.

Corollary 1

Given that Assumption 1 is satisfied and that the discharge profile of all linear quadratic game models be consistent with the total demand of the season that it describes in ETEM, then the procedure described in Babonneau et al. [5] will converge in one iteration prescribing the expansion plan that minimizes model (1).

We note that Appendix A.2 of Aliakbarisani et al. [1] also recently demonstrated that the same behavior occurs for Babonneau et al. [5]’s approach when the behavioral model is such that its expected demand response matches the demand response prescribed by model (1). Comparatively speaking, Assumption 1 is weaker as it only requires that this match occurs for the total demand response over all the time slices of each period t.