# Optimal location of EVs public charging stations based on a macroscopic urban electromobility model

Rémi Mourgues, Martin Rodriguez-Vega, Carlos Canudas de Wit

► To cite this version:

Rémi Mourgues, Martin Rodriguez-Vega, Carlos Canudas de Wit. Optimal location of EVs public charging stations based on a macroscopic urban electromobility model. CDC 2023 - 62nd IEEE Conference on Decision and Control, IEEE, Dec 2023, Singapore, Singapore. hal-04208196

HAL Id: hal-04208196

https://hal.science/hal-04208196v1

Submitted on 15 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Optimal location of EVs public charging stations based on a macroscopic urban electromobility model

Rémi Mourgues<fnref n="1" />, Martin Rodriguez-Vega<fnref n="1" />, Carlos Canudas-de-Wit<fnref n="1" />

**Abstract**— This paper introduces a graph-based dynamic model for electric vehicle (EV) mobility in urban areas. The model tracks EV state-of-charge (SoC) changes over time and space, along with power inputs from public charging stations (PCS). It considers driver behavior when deciding when and where to charge, accounting for factors like current SoC, distance to PCS, and charging cost. The model helps identify optimal PCS locations to enhance convenience for EV users and profitability for PCS owners. Additionally, an averaged version of the model is presented to reduce computational overhead while aiding in optimal PCS placement. Simulation results affirm the effectiveness of our model and optimization approach in identifying ideal charging station locations and enhancing EV charging infrastructure accessibility.

## I. INTRODUCTION

As the world transitions towards a decarbonized transportation sector, Electric Vehicles (EVs) have emerged as a promising alternative to conventional vehicles. To keep up with the predicted sharp increase of EVs, the expansion of charging infrastructure will need to be accelerated [1]. Although most of EV charging is currently done at private points at home or work [2], as the penetration rate increases, it can be expected that Public Charging Stations (PCS) will play a more important role: in some metropolitan areas around 50% to 60% of vehicles park at the roadside [3], and do not have ready access to private chargers.

The optimal placement of charging stations is crucial for maximizing the convenience of EV users and ensuring the efficient use of resources. Due to its importance, this problem has received growing attention in the literature. One of the most common approaches to solving the optimal PCS location problem is to maximize the area and number of vehicles covered by the PCS positions. In [4], the traffic network is modeled as a graph, where EV flows are defined on the links, and the objective is to maximize the sum of EV flows of links that have a PCS. [5] also uses a graph representation but optimizes the social cost instead, which takes into account recharging time, traffic congestion, and missed trips due to EVs' lack of energy. Other approaches, such as [6] and [7], consider instead the 2D area around each PCS, taking into account additional information and spatial statistics, such as population density, location of points of interest, and average occupancy of already existing PCS, to define the optimization criteria. More detailed reviews of the state of the art can be found in [8] from the point of view of

the transportation sector, and in [9] from the point of view of the power sector.

However, the existing solution approaches commonly neglect the time dimension, which is a critical factor in describing the formation of peaks of demand. This lack of adequate flow models, able to accurately forecast the spatiotemporal distribution of EVs charging demand, was identified in [8] as a key challenge.

This paper addresses the issue by presenting an enhanced Electromobility model, which delineates the spatiotemporal distribution of electric vehicles (EVs) and their average State of Charge (SoC) within an urban context. The primary contribution extends the model introduced in [10], establishing a connection between urban mobility patterns, EV flows, energy consumption, and power inputs from both private and public charging stations. The extension in this paper incorporates a driver behavior model in Section II, which calculates charging demand at each public charging station (PCS) based on factors like price, distance, EV SoC, and the attributes of competing PCS. Analytical properties of this model are examined, and a more manageable approximation is proposed in Section III.

The second contribution employs both models to formulate and address the optimal PCS location problem. Section IV's formulation considers EV recharging demand, benefits for PCS operators, and social equity in area coverage. The paper demonstrates the efficacy of this approach by solving optimization problems in a simulated scenario in Section V.

This paper offers tangible resources, rooted in authentic data, to empower policymakers, transportation planners, and industry stakeholders in their efforts to construct a more environmentally sustainable transportation network.

## II. ELECTROMOBILITY MODEL

In this section, we first introduce the original framework of the model, and then explain in more details the new contributions.

### A. EVs mobility and energy consumption model

We model the macroscopic movement of EVs in an urban metropolitan area by means of a graph-based framework, $G = \{N, L\}$ where the nodes $N$ correspond to the locations of interest where EVs move to/from, and the links $L \subset N \times N$ represent the real journeys between nodes. Following the framework introduced in [11], the nodes are divided in two types: Origins (e.g. municipalities of an urban area), and Destinations (e.g. workplaces, shopping areas, schools, leisure

<footnote n="1">¹ Univ. Grenoble Alpes, CNRS, INRIA, Grenoble INP, GIPSA-Lab, 38000 Grenoble, France {remi.mourgues@gipsa-lab.fr, martin.rodriguez-vega1@gipsa-lab.fr, carlos.canudas-de-wit@cnrs.fr}</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1: Illustration of the framework with a simple network. In orange the origin node, in purple the destination node, and in blue, the locations of two PCS. The arrows represent the links of the graph in each direction, with their corresponding flow and SoC losses. Only charging stations that are close enough to a link will provide power to it.</figcaption></figure>

<table><caption>TABLE I: Main notations of the mobility and energy consumption components of the Electromobility model</caption><thead><tr><th>Symbol</th><th>Description</th><th>Domain</th><th>Units</th></tr></thead><tbody><tr><td>ε<sub>i</sub>(t)</td><td>Avg. SoC of EVs in node i</td><td>[0, 1]</td><td>-</td></tr><tr><td>N<sub>i</sub>(t)</td><td>Number of EVs in node i</td><td>R<sub>+</sub></td><td>veh</td></tr><tr><td>φ<sub>i,j</sub>(t)</td><td>Flow of EVs from i to j</td><td>R<sub>+</sub></td><td>veh/h</td></tr><tr><td>Δε<sub>i,j</sub></td><td>SoC loss in the link (i, j)</td><td>[0, 1]</td><td>-</td></tr><tr><td>ℓ<sub>i,j</sub></td><td>Real path length of the link (i, j)</td><td>R<sub>+</sub></td><td>km</td></tr><tr><td>C</td><td>Average EVs battery capacity</td><td>R<sub>+</sub></td><td>kWh</td></tr><tr><td>ε<sup>k</sup><sub>i,j</sub></td><td>SoC at PCS k within the path (i, j)</td><td>[0, 1]</td><td>-</td></tr></tbody></table>

places). The main notations in this section are summarized in Table I.

The model states are: the number of EVs in each node $N_i(t)$, their average SoC $\varepsilon_i(t)$, and the EV flow between nodes $i \mapsto j$, $\phi_{i,j}(t)$. Each link $(i, j) \in \mathcal{L}$ connects Origin to Destination nodes (i.e. $\mathcal{G}$ is bipartite), and has road path length $\ell_{i,j}$. Each EV that traversing link $(i, j)$ losses $C\Delta\varepsilon_{i,j}$ energy, depending on the vehicle mass, average road grade, and aggregate resistive forces, as in [12]. Fig. 1 shows the main variables of interest and their interactions.

We consider the link EVs flows $\phi_{i,j}(t)$ as exogenous inputs given by the people mobility model from [11], depending on the number of people wanting to make each trip and a time profile function calibrated from real data (see [13]), coupled with a mode choice model from [10], calibrated with public survey data. The penetration rate $\eta$ of EVs is defined as the proportion of EVs in the total vehicle flow $\phi_V$.

The dynamics of $N_i$ is given by the conservation law,

$$ \dot{N}_i(t) = \sum_{j \in \mathcal{N}} (\phi_{j,i}(t) - \phi_{i,j}(t)). \quad (1) $$

The dynamics of the SoC are given by the energy balance from transported, lost, and PCS-provided energy, between nodes (see [10] for a more detailed derivation),

$$ \dot{\varepsilon}_i = \frac{1}{N_i} \sum_{j \in \mathcal{N}} \left( \phi_{j,i}(\varepsilon_j - \varepsilon_i - \Delta\varepsilon_{j,i}) + \frac{1}{C} \sum_{k=1}^{K} P_{j,i}^k \right), \quad (2) $$

where we assume there are a total of K PCS, and for each $k \in [1..K]$, $P_{j,i}^k$ is the power input from PCS k to EVs on the link $(j, i)$.

## B. Public charging stations with driver behavior

Consider that each PCS k has nominal total power $P_{tot}^k$, price per kilowatt-hour $\pi^k$, and spatial coordinates $(x^k, y^k) \in \mathcal{U}$, where $\mathcal{U}$ is the area of interest. The main notations for this section are summarized in Table II.

We assume that the EVs are primarily attracted by the PCS that are close enough to their paths. This is defined in terms of the “Attraction-law”, commonly used in human mobility applications [14],

$$ A_{i,j}^k = \frac{P_{tot}^k}{P_{ref}} \exp\left(-\left(d_{i,j}^k/\sigma\right)^2\right), \quad (3) $$

where $A_{i,j}^k$ is the "attraction-law" between the $(i, j)$-link and the k-PCS, $P_{ref}$ is a reference power parameter (e.g. the median power of charging points in many urban areas is 22 kW), $d_{i,j}^k$ is the shortest distance between the PCS and the path, and $\sigma$ is a parameter that specifies how far drivers are willing to go for a recharge. We say the link and the PCS are connected if $A_{i,j}^k$ is larger than a chosen threshold $A_{thresh}$.

We will now introduce the novel enhancements to the electromobility model, integrating driver behavior considerations. In this context, we presume that drivers assign a utility value to each charging station they can access along their route. They are then inclined to charge at stations with higher utility ratings. This addition is crucial because, without it, drivers would invariably opt to charge whenever their State of Charge (SoC) falls below 1 (as referenced in [10]). This discrete choice modeling approach has been extensively explored in the existing literature, see [15].

For that, we propose the following utility function $\mu: [0, 1] \times R_+ \times R_+ \rightarrow R$, which depends on current SoC of vehicles $\varepsilon$, the price $\pi$, and the distance to the charging station $d$,

$$ \mu(\varepsilon, \pi, d) = \gamma_{\varepsilon} \left( \frac{1}{\varepsilon} - \frac{1}{1-\varepsilon} \right) - \gamma_{\pi}\pi - \gamma_{d}e^{(d/\sigma)^2}, \quad (4) $$

where $\gamma_{\varepsilon}$, $\gamma_{\pi}$, and $\gamma_{d}$ are tunable parameters. The first term ensures that a high number of low-SoC EVs will want to charge, while only a few with high SoC will look for a PCS.

<table><caption>TABLE II: Notation summary of the charging at PCS components of the Electromobility model</caption><thead><tr><th>Symbol</th><th>Description</th><th>Domain</th><th>Units</th></tr></thead><tbody><tr><td>ℓ<sup>k</sup><sub>i,j</sub></td><td>Length along path (i, j) to PCS k</td><td>[0, ℓ<sub>i,j</sub>]</td><td>km</td></tr><tr><td>d<sup>k</sup><sub>i,j</sub></td><td>Shortest distance from k to path (i, j)</td><td>R<sub>+</sub></td><td>km</td></tr><tr><td>π<sup>k</sup></td><td>Cost of PCS k</td><td>R<sub>+</sub></td><td>€/kWh</td></tr><tr><td>P<sup>k</sup><sub>tot</sub></td><td>Total power of points in PCS k</td><td>R<sub>+</sub></td><td>kW</td></tr><tr><td>P<sub>ref</sub></td><td>Reference power for PCS attraction</td><td>R<sub>+</sub></td><td>kW</td></tr><tr><td>σ</td><td>Maximum deviation distance</td><td>R<sub>+</sub></td><td>km</td></tr><tr><td>μ<sup>k</sup><sub>i,j</sub></td><td>Utility of charging at k</td><td>R</td><td>-</td></tr><tr><td>μ&#773;<sup>k</sup><sub>i,j</sub></td><td>Utility of charging after k</td><td>R</td><td>-</td></tr><tr><td>γ<sub>ε</sub>, γ<sub>π</sub>, γ<sub>d</sub></td><td>Weight parameters for utility associated to SoC, price, and distance, respectively</td><td>R</td><td>-</td></tr><tr><td>β<sup>k</sup><sub>i,j</sub></td><td>Proportion of flow that wants to charge in PCS k</td><td>[0,1]</td><td>-</td></tr><tr><td>D<sup>k</sup><sub>i,j</sub></td><td>Charging demand at PCS k</td><td>R<sub>+</sub></td><td>kW</td></tr><tr><td>S<sup>k</sup><sub>i,j</sub></td><td>Charging supply of k into (i, j)</td><td>R<sub>+</sub></td><td>kW</td></tr><tr><td>P<sup>k</sup><sub>i,j</sub></td><td>Power provided by PCS k to link (i, j)</td><td>R<sub>+</sub></td><td>kW</td></tr></tbody></table>
The second and third terms specify the drivers' sensitivity to the price and distance to the PCS, respectively. The distance-related penalty is low for $d < \sigma$, but increases sharply afterwards. Other definitions for the utility could be proposed which consider other variables and interactions according to the application.

Denote by $\mu_{i,j}^k = \mu(\varepsilon_i, \pi^k, d_{i,j}^k)$ the utility of charging at PCS $k$ while traveling along the path $(i, j)$, and by $\bar{\mu}_{i,j}^k$ the utility of charging somewhere downstream, computed as

$$ \bar{\mu}_{i,j}^k = \max_{p|k<p} \{\mu_{i,j}^p\}, \quad (5) $$

where the notation $k \prec p$ means that $k$ precedes (is upstream of) $p$ along the same path, $\ell_{i,j}^k < \ell_{i,j}^p$. If there are no other PCS after $k$, we set $\bar{\mu}_{i,j}^k = 0$. The proportion of EVs that charge at $k$ is

$$ \beta_{i,j}^k = \frac{e^{\mu_{i,j}^k}}{e^{\mu_{i,j}^k} + e^{\bar{\mu}_{i,j}^k}}. \quad (6) $$

**Note:** $\beta_{i,j}^k \in [0, 1]$, $\lim_{\mu_{i,j}^k \to -\infty} \beta_{i,j}^k = 0$ and $\lim_{\mu_{i,j}^k \to \infty} \beta_{i,j}^k = 1$.

Thus the power provided by each station to each connected link is

$$ P_{i,j}^k = \min(D_{i,j}^k, S_{i,j}^k), \quad (7) $$

where the charging demand is

$$ D_{i,j}^k = C(1 - \varepsilon_i + \Delta\varepsilon_{i,j}^k)\phi_{i,j}\beta_{i,j}^k \prod_{p<k} (1 - \beta_{i,j}^p). \quad (8) $$

where $\Delta\varepsilon_{i,j}^k = \Delta\varepsilon_{i,j}^k \ell_{i,j}^k / \ell_{i,j}^k$, with $\ell_{i,j}^k$ the path length from $i$ to PCS $k$. The available PCS power supply is divided proportionally between all connected paths,

$$ S_{i,j}^k = \frac{D_{i,j}^k}{\sum_{(n,m)} D_{n,m}^k} P_{tot}^k. \quad (9) $$

**Example 1 (System with two-nodes and one CS).** Consider the example shown in Fig. 1 with only PCS 1 active. The system dynamics are:

$$ \begin{aligned} \dot{\varepsilon}_1 &= \frac{1}{N_1} \left( \phi_{2,1}(\varepsilon_2 - \varepsilon_1 - \Delta\varepsilon_{2,1}) + \frac{1}{C} P_{2,1}^1(\varepsilon_2) \right), & (10) \\ \dot{\varepsilon}_2 &= \frac{1}{N_2} \phi_{1,2}(\varepsilon_1 - \varepsilon_2 - \Delta\varepsilon_{1,2}). \end{aligned} $$

with $P_{2,1}^1(\varepsilon_2) = \min(C\beta_{2,1}^1\phi_{2,1}(1 + \frac{\ell_{2,1}^1}{\ell_{2,1}^1}\Delta\varepsilon_{2,1} - \varepsilon_2), P)$, and

$$ \begin{aligned} \dot{N}_1 &= \phi_{2,1}(t) - \phi_{1,2}(t), \\ \dot{N}_2 &= \phi_{1,2}(t) - \phi_{2,1}(t). \end{aligned} \quad (11) $$

### C. General model representation

Denote by $N(t) \in \mathbb{R}_+^{|N|}$, and $\varepsilon(t) \in [0, 1]^{|N|}$ the concatenated states of number of EVs $N_i(t)$ and SoC $\varepsilon_i(t)$, respectively, and by $\phi(t) \in \mathbb{R}_+^{|\mathcal{L}|}$ the concatenated vector of link flows $\phi_{i,j}(t)$.

Furthermore, denote by $u \in \mathbb{R}^{K \times 2}$ the location of K PCS, i.e.,

$$ u = [(x^1, y^1), (x^2, y^2), \dots (x^K, y^K)]^\top. $$

The model can be compactly written as

$$ \dot{N} = \mathcal{I}\phi(t), \quad (12) $$

$$ \dot{\varepsilon} = f(N(t), \phi(t), \varepsilon(t), u) = f(t, \varepsilon, u), \quad (13) $$

where $\mathcal{I} \in \{-1, 0, 1\}^{|N| \times |\mathcal{L}|}$ is the graph's incidence matrix, and $f$ is a vector function that gathers (2)–(7). As $\phi$ and $N$ do not depend on $\varepsilon$, $f$ can be considered as a function of time for ease of notation.

### III. MODEL PROPERTIES

In this section we study several model properties which are inherent to the physic of the system. In particular we analyze the average states and their stability, which are useful for simplifying the optimization problem and better understanding the aggregate system behavior.

#### A. Average system approximation

The proposed model was designed for the case of urban areas, which are characterized by periodic mobility patterns during a week. As a consequence, EV flows and numbers of EVs at each node resulting from solving (12) evolve periodically with a period $T$ (typically one week) i.e., $N(t) = N(t-T)$, $\phi(t) = \phi(t-T)$, Then, $\bar{N}$, and $\bar{\phi}$

$$ \bar{N} = \frac{1}{T} \int_{t-T}^{t} N(\tau) d\tau, \quad \bar{\phi} = \frac{1}{T} \int_{t-T}^{t} \phi(\tau) d\tau $$

describe the average vehicle number per node, and the average vehicle flow between nodes, respectively. Note that the periodicity of $N$ implies

$$ \bar{\phi}_{i,j} = \bar{\phi}_{j,i}, \quad \forall (i, j). \quad (14) $$

As a consequence, $f(t, \varepsilon, u)$ is periodic with respect to $t$, $f(t, \varepsilon, u) = f(t - T, \varepsilon, u)$. Hence, (13) can be used alone to study the model dynamic properties. For that, consider the following definition for the average SoC,

$$ \bar{\varepsilon}(t) = \frac{1}{T} \int_{t-T}^{t} \varepsilon(\tau) d\tau \quad (15) $$

which has dynamics

$$ \dot{\bar{\varepsilon}}(t) = \frac{1}{T} \int_{t-T}^{t} f(N, \phi, \varepsilon, u) d\tau. \quad (16) $$

**Lemma 1 (Approximated averaged equation).** The dynamic of $\bar{\varepsilon}(t)$ can be approximated using the Taylor expansion of $f$ around $(\bar{N}, \bar{\phi}, \bar{\varepsilon})$ by

$$ \dot{\bar{\varepsilon}}(t) \approx f(\bar{N}, \bar{\phi}, \bar{\varepsilon}(t), u). \quad (17) $$

**Proof.** Expanding $f(N, \phi, \varepsilon, u)$ around $(\bar{N}, \bar{\phi}, \bar{\varepsilon})$ and neglecting high-order terms, we get

$$ f(N, \phi, \varepsilon, u) \approx f(\bar{N}, \bar{\phi}, \bar{\varepsilon}, u) + \frac{\partial f}{\partial \varepsilon}(\varepsilon - \bar{\varepsilon}) + \frac{\partial f}{\partial \bar{N}}(N - \bar{N}) + \frac{\partial f}{\partial \bar{\phi}}(\phi - \bar{\phi}). $$

Noticing that, by construction, the integral of a period $T$ of the difference between the signals and their average do cancel, i.e.

$$ \int_{t-T}^{t} (\varepsilon - \bar{\varepsilon}) d\tau = \int_{t-T}^{t} (N - \bar{N}) d\tau = \int_{t-T}^{t} (\phi - \bar{\phi}) d\tau = 0, \\ \text{as a consequence } f(N, \phi, \varepsilon, u) \approx f(\bar{N}, \bar{\phi}, \bar{\varepsilon}, u). \quad \square $$

**Example 2 (Average approximation. Example 1-cont.).** With regard to the Example 1, and following Lemma 1, we obtain

$$ \begin{aligned} \dot{\varepsilon}_1 &= \frac{1}{\bar{N}_1} \left( \bar{\phi}(\bar{\varepsilon}_2 - \bar{\varepsilon}_1 - \Delta\bar{\varepsilon}_{2,1}) + \frac{1}{C} P_{2,1}(\bar{\varepsilon}_2) \right), & (18) \\ \dot{\varepsilon}_2 &= \frac{1}{\bar{N}_2} \bar{\phi}(\bar{\varepsilon}_1 - \bar{\varepsilon}_2 - \Delta\bar{\varepsilon}_{1,2}) \end{aligned} $$
<figure><img src="image_3.png" /><figcaption>Fig. 2: Example 1. Comparison between the full model (plain line) and the average model (dashed line).</figcaption></figure>

with $P_{2,1}$ given by:

$$P_{2,1}(\bar{\varepsilon}_2) = \min \left( C \bar{\phi}(1 - \bar{\varepsilon}_2)(1 - \bar{\varepsilon}_2 + \frac{1}{2} \Delta \varepsilon_{2,1}), P \right)$$

where $\beta_{2,1}$ in (6) is replaced by a linear approximation. Details are given in the Appendix I. The comparison between the original and average model is shown in Figure 2.

## B. Asymptotic-periodicity, Multiple equilibria and stability

In view of analyzing the dynamical properties of our model, let introduce the following definitions and properties.

**Definition 1 (Asymptotic-periodicity).** System (13) is said to be asymptotically periodic if its solutions satisfy

$$\lim_{t \to \infty} (\varepsilon(t) - \varepsilon(t-T)) = 0. \quad (19)$$

**Property 1 (Existence of Equilibria).** If System (13) is asymptotically periodic, then there exist one or more equilibrium states for its average counterpart $\bar{\varepsilon}(t)$,

$$\lim_{t \to \infty} \bar{\varepsilon}(t) = \bar{\varepsilon}^* \Leftrightarrow \lim_{t \to \infty} \dot{\bar{\varepsilon}}(t) = 0 \quad (20)$$

where $\bar{\varepsilon}^*$ is the average equilibrium state vector. Assuming that (17) holds exactly, then $\bar{\varepsilon}^*$ satisfies the property

$$f(\bar{N}, \bar{\phi}, \bar{\varepsilon}^*, u) = 0. \quad (21)$$

Note that $\bar{\varepsilon}^*$ will ultimately depend on the penetration rate $\eta$ so we may expect different positions and nature of equilibria depending on the total amount of EVs.

Besides, and because the structure of $P_{2,1}(\bar{\varepsilon}_2) = \min(D_{2,1}, S_{2,1})$, the system (18) (a similar behavior can be obtained for more complex systems) will exhibit different equilibrium and stability properties for the cases of: demand regime $D$, when $P_{2,1}(\bar{\varepsilon}_2) = D_{2,1}$, and the supply regime $S$, when $P_{2,1}(\bar{\varepsilon}_2) = S_{2,1}$, separated by the line $\bar{\varepsilon}_2 = \bar{\varepsilon}_{2c}$, where

$$\bar{\varepsilon}_{2c} := 1 + \frac{\Delta \varepsilon_{2,1}}{4} - \sqrt{\frac{\Delta \varepsilon_{2,1}^2}{16} + \frac{\eta_c}{\eta}(\Delta \varepsilon_{2,1} + \Delta \varepsilon_{1,2})}, \quad (22)$$

with the critical penetration rate $\eta_c$ given as

$$\eta_c := \frac{P}{\bar{\phi}_V C (\Delta \varepsilon_{1,2} + \Delta \varepsilon_{2,1})}. \quad (23)$$

it defines the transition between different systems regimes (as introduced by [10]). Its implications in the systems stability are given in Property 3.

**Definition 2 (Sustainable and Unsustainable demand regime).** We call sustainable (resp. unsustainable) demand regime (SDR, resp. UDR) the regime when the total power capacity of the PCS, is enough (resp. not enough) to satisfy

<figure><img src="image_4.png" /><figcaption>Fig. 3: Example 1. Evolution of the mean SoC of the system using the average model, for different $\eta$. The case $\eta > \eta_c$ reaches a plateau representing the fact that the cars are stuck when the SoC comes near the loss of the path.</figcaption></figure>

the total EVs energy demand. The SoC trajectories then reach a non-zero value (resp. converge towards zero).

We now study the stability and convergence of the system, starting with a general property.

**Property 2 (Stability).** If System (13) is asymptotically-periodic, then the average equilibrium state $\bar{\varepsilon}^*$ is stable, in the sense that the system converges to one of the values of $\bar{\varepsilon}^*$ satisfying $f(\bar{N}, \bar{\phi}, \bar{\varepsilon}^*, u) = 0$.

Asymptotic periodicity, and hence stability, can only be obtained for certain cases depending on the values of $\eta$. Indeed, as will be shown in Property 3, the system (18) exhibits a phase transition at the critical penetration rate $\eta_c$ between SDR and UDR.

**Property 3.** The system (18) has the following stability properties:

i) The following statements hold iff $\eta < \eta_c$:

a. There is an unique equilibrium point, $\bar{\varepsilon}_D^* \in D$ satisfying $f(\bar{N}, \bar{\phi}, \bar{\varepsilon}_D^*, u) = 0$,

b. $\bar{\varepsilon}_D^*$ is asymptotically stable: $\exists \delta$ s.t. $\forall \bar{\varepsilon}(0) \in B_\delta = \{\bar{\varepsilon} \in D : \|\bar{\varepsilon} - \bar{\varepsilon}_D^*\| \le \delta\}, \lim_{t \to \infty} \bar{\varepsilon}(t) = \bar{\varepsilon}_D^*$.

ii) Let $S^* := \{\bar{\varepsilon}^* \in S : f(\bar{N}, \bar{\phi}, \bar{\varepsilon}^*, u) = 0\}$, and $s(\bar{\varepsilon}) := \bar{\varepsilon}_1 - \bar{\varepsilon}_2 - \Delta\varepsilon_{1,2}$. Then:

a. $\forall \eta$, the trajectories $\bar{\varepsilon}(t)$ in $S$ converge exponentially towards $s(\bar{\varepsilon}) = 0$

b. $S^* \neq \emptyset$ iff $\eta = \eta_c$. Then, the set of equilibrium points, $\bar{\varepsilon}^* \in S^*$, belong to $s(\bar{\varepsilon}^*) = 0$.

iii) $\forall \eta \neq \eta_c, \forall \bar{\varepsilon}(0) \in S$, the trajectories $\bar{\varepsilon}(t)$ in $S$ reach in finite time:

a. Negative values by crossing the $\bar{\varepsilon}_2 = 0$ boundary if $\eta > \eta_c$;

b. $D$ by crossing the $S/D$ boundary if $\eta < \eta_c$.

**Proof.** See Appendix I.

**Remarks.** $\eta_c$ in (23) separates out power balance between the injected PCS power and power losses of the network. Case i) and ii) – b state that the SDR is reached for all $\eta \le \eta_c$. When $\eta = \eta_c$, an infinite number of equilibria exists characterized by $s(\bar{\varepsilon}^*) = 0$. ii) a indicates that the surface $s(\bar{\varepsilon}^*) = 0$ is always attractive. iii) – a corresponds to the UDR: the system loses more energy that the PCS can provide. In this case, the system trajectories reach in finite time the $\bar{\varepsilon}_2 = 0$ boundary. iii) – b corresponds to the
<figure><img src="image_5.png" /><figcaption>Fig. 4: Phase portraits of system (18) for different values of η. The circles represent the final state.</figcaption></figure>

SDR: the trajectories cross to the demand domain $\mathcal{D}$ and potentially reach an equilibrium $\bar{\varepsilon}_D^*$. Although those results do not include the study of potential oscillations between the two regimes, we can observe from the phase portraits in Fig. 4 that no cycles are observed at the vicinity of the $S/\mathcal{D}$ boundary when $\eta \le \eta_c$, and that the only equilibrium is indeed $\bar{\varepsilon}_D^*$ for $\eta < \eta_c$.

#### IV. PROBLEM FORMULATION OF THE OPTIMAL PLACEMENT OF CHARGING STATIONS

Consider the problem of finding the optimal locations, $u$, of $K$ PCS, with nominal powers $P_{\text{tot}}^k, k \in [1, K]$.

##### A. Objective function

Let $J(\varepsilon, u)$ be the cost function considering profits for the charge stations operators, user satisfaction, and social equity,

$$J(\varepsilon, u) = \underbrace{\lambda_1 J_1(\varepsilon, u)}_{\text{profit}} + \underbrace{\lambda_2 J_2(\varepsilon, u)}_{\text{satisfaction}} + \underbrace{\lambda_3 J_3(u)}_{\text{equity}} \quad (24)$$

where $\lambda_i > 0$ are weighting constants selected such that $0 \le J(\varepsilon, u) \le J_M$, with $J_i$ defined as follows:

* **Profit:** Highlights the priority of PCS operators in striving to maximize charging station profits through the minimization of the gap between total PCS capacity and the power actually supplied (sold)

$$J_1(\varepsilon, u) = \sum_k \left( P_{\text{tot}}^k - \sum_i \sum_j P_{i,j}^k(\varepsilon, u) \right). \quad (25)$$

* **Satisfaction:** describes the difference between the EV energy demand at a specific link and the expected power supply with PSC placement at location $u$. By minimizing this utility function, the satisfaction of EV drivers is maximized. $J_2$ is defined as

$$J_2(\varepsilon, u) = \sum_i \sum_j \left( D_{i,j}(\varepsilon) - \sum_k P_{i,j}^k(\varepsilon, u) \right), \quad (26)$$

where $D_{i,j}(\varepsilon)$ is the power demand on link $(i, j)$,

$$D_{i,j} = C \phi_{i,j} (1 + \Delta \varepsilon_{i,j} - \varepsilon_i) \quad (27)$$

* **Equity:** pertains to the geographical dispersion of PCS. Its purpose is to compel the optimizer to diversify PCS

locations, ensuring broader access for more users. This is done by minimizing

$$J_3(u) = \sum_{i \neq j} \sigma \left( \frac{P_{\text{tot}}^i P_{\text{tot}}^j}{\|u_i - u_j\|^2} \right), \quad (28)$$

where $\sigma: \mathbb{R}_+ \to [0, 1]$ is a sigmoid function.

##### B. Static optimization based on the average model

In this first formulation, we consider the averaged system $\bar{\varepsilon}(t)$, with the assumption that approximation (17) holds, and that the system is asymptotically periodic. Hence, from Property 1, system (17) will reach an average equilibrium. The optimization is then formulated as follows:

**Problem formulation 1 (Static optimization based on the average equilibrium).**

$$\begin{array}{ll} \min_{u \in \mathcal{U}} & J(\bar{\varepsilon}^*, u) \\ \text{s.t.} & f(\bar{N}, \bar{\phi}, \bar{\varepsilon}^*, u) = 0, \end{array} \quad (29)$$

where $\bar{N}, \bar{\phi}$ are pre-computed in advance.

##### C. Dynamic optimization based on the full model

Another alternative (more computationally involved), and in view to validate the previous simplified problem formulation, is to use the full model (12)-(13),

**Problem formulation 2 (Dynamic optimization based on the full model).**

$$\min_{u \in \mathcal{U}} \frac{1}{T} \int_{H-T}^{H} J(\varepsilon(t), u)dt \quad (30)$$

s.t. $\dot{\varepsilon}(t) = f(N(t), \phi(t), \varepsilon(t), u)$,

where $N(t)$ is pre-computed from (12) and stored during the whole time-optimization horizon $t = H$. The value of $H$ should be large enough such that $\varepsilon(t)$ is close to the asymptotic periodic behavior by the end.

#### V. OPTIMAL PLACEMENT STRATEGIES

Consider the network shown in Fig. 5-6, consisting of 4 origin and 5 destination nodes. Each origin has a population of 100 EVs. The EV flows of the network were set to emulate realistic daily commute cycles: they are shaped as
a Gaussian bell curve with a standard deviation of 1 hour, centered around 7h00 for origin-to-destination flows, and around 17h00 for destination-to-origin. The amplitudes are such that 95% of the total population moves during the day, and distributed randomly over all links. In these figures, the thickness of the graph links is proportional to the peak flows.

## A. Optimal placement strategies for a single PCS

1) Brute force approach and heatmap computation using the average model: To benchmark the solutions of the optimization Problems 1-2, we first compute the cost function $J$ for all possible locations of a single PCS. This allows to plot a heatmap of $J$ in the considered area, and also to identify the "true" optimal location. The network area is spatially discretized using a rectangular grid $\mathcal{U}$. Each rectangle is labeled by the index $k$, for all $k \in [1..K]$. Using the time-dependent flows as inputs we compute $\bar{\phi}_{i,j}$ for each $(i, j) \in \mathcal{L}$ and $\bar{N}_i$ for $i \in \mathcal{N}$. The objective function is computed for each of the candidate locations, and the optimal position $u_{opt}$ is set at the location with smallest $J$, as described by Alg. 1. Figure 5 shows the heatmap of the objective function for each position of the PCS, and the position of $u_{opt}$.

<caption>**Algorithm 1** Brute force approach and heatmap computation
- Average model</caption>

**Inputs:** $\phi(t)$
**Outputs:** $\{J(\bar{\varepsilon}_k^*, u_k)\}_{k=1...K}, J_{opt}, u_{opt}$
Initialize $J_{opt} \leftarrow J_M, u_{opt} \leftarrow (0, 0)$
From (12) compute $\bar{\phi}_{i,j}, \forall (i, j) \in \mathcal{L}$, and $\bar{N}_i, \forall i \in \mathcal{N}$
**for** $k = 1$ **to** $\mathcal{K}$ **do**
Let $u_k = (x_k, y_k)$
Find $\bar{\varepsilon}_k^*$ from $f(\bar{N}, \bar{\phi}, \bar{\varepsilon}_k^*, u_k) = 0$
**if** Solution exist **then**
Compute $J(\bar{\varepsilon}_k^*, u_k)$ using (24)
**else**
$J(\bar{\varepsilon}_k^*, u_k) \leftarrow J_M$
**end if**
**if** $J(\bar{\varepsilon}_k^*, u_k) < J_{opt}$ **then**
$J_{opt} \leftarrow J(\bar{\varepsilon}_k^*, u_k), u_{opt} \leftarrow u_k$
**end if**
**end for**

2) Optimization using the average model: We can solve the Problem 1 using the generic Matlab optimizer `fmincon`, which takes as input the objective function $J$, the constraints $f = 0$, and an initial point $u_0$ to compute the optimal location $u_{opt}$ using an interior-point method (only `fmicon` was used in this paper, other solvers will be studied in a future work.) Since Problem 1 is non-convex, it is possible to end up in a local minimum. To improve the optimization, we repeat the optimization $K$ times using random initial points and select the output that provides the best objective, as explained in Alg. 2. Fig. 5 shows the obtained optimal location of the PCS using $K = 5$, which in this case coincides with the "true" optimal location found with the brute force benchmark.

<figure><img src="image_6.png" /><figcaption>Fig. 5: Optimal location of 1 PCS for an arbitrary network. The background colors correspond to a heatmap of the objective function computed by locating the PCS at each coordinate in the domain. Whiter colors shown smaller values for $J$. The green square is the optimal location computed by an off-the-shelf solver.</figcaption></figure>

<caption>**Algorithm 2** Optimization with solver - Average model</caption>

**Inputs:** $\phi(t)$
**Outputs:** $J_{opt}, u_{opt}$
Initialize $J_{opt} \leftarrow J_M, u_{opt} \leftarrow (0, 0)$
From (12) compute $\bar{\phi}_{i,j}, \forall (i, j) \in \mathcal{L}$, and $\bar{N}_i, \forall i \in \mathcal{N}$
**for** $k = 1$ **to** $\mathcal{K}$ **do**
Choose $u_0 = (x_0, y_0) \in \mathcal{U}$ randomly
($u_k, J_k) \leftarrow \text{fmincon}(J, f(\bar{N}, \bar{\phi}, \cdot, \cdot), u_0)$
**if** $J_k < J_{opt}$ **then**
$J_{opt} \leftarrow J_k, u_{opt} \leftarrow u_k$
**end if**
**end for**

3) Brute force approach and heatmap computation using full model: Similarly to the average model analysis, we first study the solutions of the optimization Problem 2 using a brute force approach for the location of a single PCS, with the same spatial grid discretization $\mathcal{U}$. The network under study has a period $P = 1$ day, and we used a simulation horizon $H = 1$ week. In the full case, the input flows $\phi(t)$ are used to compute the number of EVs $N(t)$. To build the heatmap, for each candidate location of the PCS, the trajectories $\varepsilon(t)$ for $t \in [0, H]$ are computed, and the objective value is the average of the last period of the simulation. The optimal location $u_{opt}$ is the place with the smallest objective value $J$. For this, we used Alg. 3. The resulting heatmap is shown in Fig. 5.

4) Optimization using full model: To compare with the brute-force benchmark given previously, the `fmincon` optimizer is used as described in Alg. 4, which has a similar
**Algorithm 3** Brute force approach and heatmaps computation - full model

**Inputs:** $\phi(t), T, H$

**Outputs:** $\{J(\bar{\varepsilon}_k^*, u_k)\}_{k=1...K}, J_{opt}, u_{opt}$

Initialize $J_{opt} \leftarrow J_M, u_{opt} \leftarrow (0, 0)$

Compute $N(t)$ for $t \in [0, H]$ using (12)

**for** $k = 1$ to $K$ **do**

Let $u_k = (x_k, y_k)$

Compute $\varepsilon(t)$ for $t \in [0, H]$ from (13) using $u_k$

Compute $J_k = \frac{1}{T} \int_{H-T}^{H} J(\varepsilon(t), u_k)dt$ with (24)

**if** $J_k < J_{opt}$ **then**

$$J_{opt} \leftarrow J_k, u_{opt} \leftarrow u_k$$

**end if**

**end for**

structure to Alg. 2. Alg. 4 uses of the function $\circ obj$ to compute the value of the dynamic objective function for a given PCS position $u$. The resulting optimal PCS location $u_{opt}$ is shown in Fig.5, which coincides with the global optimal obtained with the brute force approach.

**Algorithm 4** Optimization with solver - full model

**Inputs:** $u, \phi(t), N(t), T, H$

**Outputs:** $J_{opt}, u_{opt}$

Initialize $J_{opt} \leftarrow J_M, u_{opt} \leftarrow (0, 0)$

Compute $N(t)$ for $t \in [0, H]$ using (12)

**for** $k = 1$ to $K$ **do**

Choose $u_0 = (x_0, y_0) \in \mathcal{U}$ randomly

$(u_k, J_k) \leftarrow \text{fmincon}(\circ obj(\cdot, \phi, N, T, H), u_0)$

**if** $J_k < J_{opt}$ **then**

$$J_{opt} \leftarrow J_k, u_{opt} \leftarrow u_k$$

**end if**

**end for**

**function** $\circ obj(u, \phi(t), N(t), T, H)$

Compute $\varepsilon(t)$ for $t \in [0, H]$ from (13)

Compute $J = \frac{1}{T} \int_{H-T}^{H} J(\varepsilon(t), u)dt$ with (24)

**return** $J$

**end function**

## B. Optimal location of multiple PCS

For multiple PCS, it is no longer feasible to use the brute force approaches in Alg. 1-3, as the number of computation complexity increases exponentially with the number of PCS. Therefore, we use Alg. 2 for the average model and Alg. 4 for the full model to locate 5 PCS. For the comparison to be fair, we use the same random sequence for the optimal search. The results are shown in Fig. 6, showing that both algorithms provide similar results. This was expected since the average formulation is designed to simplify the asymptotic-periodicity of the dynamic system, and contains the same spatial and vehicle demand properties of the full model. However, there could be situations where both models provide slightly different results, especially when the flows dynamics variance are large (e.g. the peak flows in different links occur during very different times of day). In terms of computation time, the averaged model optimization took 1.4 minutes to provide the results, whereas the full model

<figure><img src="image_7.png" /><figcaption>(a) Average model</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) full model</figcaption></figure>

Fig. 6: Optimal location of 5 PCS for an arbitrary network, using a solver. The same random initial conditions have been used. Note that the methods give similar results, but the first is 40 times faster.

required 53.5 minutes. This suggests that optimization based on the averaged model can be adopted when applying the methods to large-scale networks.

## VI. CONCLUSIONS

We studied the urban EV PCS location problem, by first refined the Electromobility model [10] adding a driver behavior model, enabling a more realistic assessment of driver demand. Then, we also analyzed the dynamic model properties and proposed a simplified, time-independent average model based on urban mobility patterns. We formulated two optimal placement problems using both models, considering EV user satisfaction, PCS operator profitability, and social equity in the cost function. Simulation results showed that both approaches yield similar satisfactory results, suggesting the averaged model improves computation efficiency.

# APPENDIX I
## PROOF OF PROPERTY 3

### A. Approximation of $\beta_{2,1}$

First, $\beta_{2,1}$ in (6) is approximated as

$$\beta_{2,1}(\bar{\varepsilon}_2) = \frac{1}{1 + \exp\left(-\frac{\gamma_{\varepsilon}}{\bar{\varepsilon}_2} + \frac{\gamma_{\varepsilon}}{1-\bar{\varepsilon}_2} + \gamma_{\pi}\pi\right)} \approx 1 - \bar{\varepsilon}_2 \quad (31)$$

<figure><img src="image_9.png" /><figcaption>Fig. 7: Shape of $\beta_{2,1}$ and its proposed linear approximation.</figcaption></figure>
where $\gamma_{\pi} = 0$ and $\gamma_{\varepsilon} = 1$ as the electricity price does not play any role in the single PCS case. The quality of this approximation is shown on the Figure 7. This approximation leads to the $P_{2,1}(\bar{\varepsilon}_2)$ given in Example 2.

**Proof of i): Stability in the demand regime.** It is straightforward to show from (18) that, in $\mathcal{D}$, the equilibrium is given by

$$ \begin{aligned} \bar{\varepsilon}_1^* &= \bar{\varepsilon}_2^* + \Delta\varepsilon_{1,2}, \\ \bar{\varepsilon}_2^* &= 1 + \frac{\Delta\varepsilon_{2,1}}{4} - \sqrt{\frac{\Delta\varepsilon_{2,1}^2}{16} + \Delta\varepsilon_{1,2} + \Delta\varepsilon_{2,1}}. \end{aligned} \quad (32) $$

In this case the solution is uniquely defined, and denoted $\bar{\varepsilon}_D^*$. It does not depend on the value of $\eta$, however in the case $\eta > \eta_c$, $\bar{\varepsilon}_D^* < \bar{\varepsilon}_{2c}$ so $\bar{\varepsilon}_D^* \in \mathcal{S}$.

Consider now the local approximation of system (18) around $\bar{\varepsilon}_D^* \in \mathcal{D}$, when $\delta\bar{\varepsilon}$ is a perturbation:

$$ \dot{\delta\bar{\varepsilon}} = \nabla f(\bar{\varepsilon}^*) \cdot \delta\bar{\varepsilon} \quad (33) $$

where $\nabla f$ is the Jacobian matrix of $f$:

$$ \begin{aligned} \nabla f(\bar{\varepsilon}^*) &= \eta \bar{\phi}_V \\ & \left( -\frac{1}{\bar{N}_1}, \frac{1}{\bar{N}_1} - \frac{2}{\bar{N}_1} \sqrt{\frac{\Delta\varepsilon_{2,1}^2}{16} + \Delta\varepsilon_{1,2} + \Delta\varepsilon_{2,1}} - \frac{1}{\bar{N}_2} \right). \end{aligned} \quad (34) $$

It can be shown from it that both eigenvalues are strictly negative, therefore the system is stable.

**Proof of ii): Stability in the supply regime.** From (18), for $\bar{\varepsilon} \in \mathcal{S}$, the equilibrium exists iff $\eta = \eta_c$, and it is characterized by

$$ s(\bar{\varepsilon}^*) = \bar{\varepsilon}_1^* - \bar{\varepsilon}_2^* - \Delta\varepsilon_{1,2} = 0. \quad (35) $$

To show that $s = s(\bar{\varepsilon}_1, \bar{\varepsilon}_2)$ is attractive in $\mathcal{S}$, consider the system solutions in $\mathcal{S}$:

$$ \begin{aligned} \bar{\varepsilon}_1(t) &= \alpha t - \frac{\rho}{a\bar{N}_1}(1 - e^{-at}) + \bar{\varepsilon}_1(0), \\ \bar{\varepsilon}_2(t) &= \alpha t + \frac{\rho}{a\bar{N}_2}(1 - e^{-at}) + \bar{\varepsilon}_2(0), \end{aligned} \quad (36) $$

with constants defined in the Table III. $s(t)$ becomes:

$$ s(t) = s(0)e^{-at} $$

and hence $s(t) \to 0$, regardless the value of $\eta$.

When $\eta = \eta_c$, $\alpha = 0$ and the system converges exponentially towards the solution:

$$ \begin{aligned} \bar{\varepsilon}_1^* &= \bar{\varepsilon}_1(0) - \rho/(a\bar{N}_1), \\ \bar{\varepsilon}_2^* &= \bar{\varepsilon}_2(0) + \rho/(a\bar{N}_2). \end{aligned} \quad (37) $$

Note that some trajectories can cross the $\mathcal{S}/\mathcal{D}$ boundary, however in that case $\bar{\varepsilon}_{D2}^* = \bar{\varepsilon}_{2c}$ so $\bar{\varepsilon}_D^* \in \mathcal{S}^*$. Moreover, as said in the main part, no cycle have been noticed and all trajectories converge towards their respective equilibrium.

When $\eta \neq \eta_c$, it is straightforward to see that no equilibrium exists.

**Proof of iii): Stability in the supply regime ($\eta \neq \eta_c$).**

In the case $\eta \neq \eta_c$, there is no equilibrium in $\mathcal{S}$, although the $s = 0$ hyperplane is still attractive. It can be shown that $\text{sign}(\alpha) = \text{sign}(\eta_c - \eta)$, from (36) the system trajectories will diverge in finite time towards:

* negative values for $\alpha < 0$, reaching the UDR;

* the $\mathcal{S}/\mathcal{D}$ boundary for $\alpha > 0$, then will follow the demand regime dynamic. As discussed previously, they then seem to converge towards $\bar{\varepsilon}_D^*$.

<table><caption>TABLE III: Constants of the dynamic evolution of the system for the supply regime.</caption><thead><tr><th>Constant</th><th>Value</th></tr></thead><tbody><tr><td>a</td><td>η̄φ<sub>V</sub>(1/̄N<sub>1</sub> + 1/̄N<sub>2</sub>)</td></tr><tr><td>b</td><td>Δε<sub>2,1</sub>̄φ<sub>V</sub>(η<sub>c</sub> - η) + Δε<sub>1,2</sub>̄φ<sub>V</sub>(η<sub>c</sub>/N<sub>1</sub> + η/N<sub>2</sub>)</td></tr><tr><td>α</td><td>η̄φ<sub>V</sub>/N<sub>2</sub>(b/a - Δε<sub>1,2</sub>)</td></tr><tr><td>ρ</td><td>η̄φ<sub>V</sub>(ε̄<sub>1</sub>(0) - ε̄<sub>2</sub>(0) - b/a)</td></tr></tbody></table>

REFERENCES

[1] European Federation for Transport and Environment, “Recharge EU: how many charge points will Europe and its Member States need in the 2020s,” Transport and Environment, p. 68, 2020.

[2] S. Wolff and R. Madlener, “Charged up? Preferences for Electric Vehicle Charging and Implications for Charging Infrastructure Planning,” SSRN Electronic Journal, no. 3, 2019.

[3] CEREMA and S. mixte des transports en commun de l'agglomération grenobloise, “Enquête Ménages Déplacements, Grenoble / Grande région grenobloise - 2010,” 2010.

[4] Y. W. Wang and C. C. Lin, “Locating multiple types of recharging stations for battery-powered electric vehicle transport,” Transportation Research Part E: Logistics and Transportation Review, vol. 58, pp. 76–87, 2013.

[5] F. He, Y. Yin, and J. Zhou, “Deploying public charging stations for electric vehicles on urban road networks,” Transportation Research Part C: Emerging Technologies, vol. 60, pp. 227–240, 2015.

[6] B. J. Mortimer, C. Hecht, R. Goldbeck, D. U. Sauer, and R. W. De Doncker, “Electric Vehicle Public Charging Infrastructure Planning Using Real-World Charging Data,” World Electric Vehicle Journal, vol. 13, no. 6, pp. 1–18, 2022.

[7] G. Dong, J. Ma, R. Wei, and J. Haycox, “Electric vehicle charging point placement optimisation by exploiting spatial statistics and maximal coverage location models,” Transportation Research Part D: Transport and Environment, vol. 67, no. November 2018, pp. 77–88, 2019.

[8] M. Kchaou-Boujelben, “Charging station location problem: A comprehensive review on models and solution approaches,” Transportation Research Part C: Emerging Technologies, vol. 132, no. November, p. 103376, 2021.

[9] F. Ahmad, A. Iqbal, I. Ashraf, M. Marzband, and I. Khan, “Optimal location of electric vehicle charging station and its impact on distribution network: A review,” Energy Reports, vol. 8, pp. 2314–2333, 2022.

[10] M. Rodriguez-vega, C. Canudas-de Wit, G. D. Nunzio, and B. Othman, “A Graph-Based Mobility Model for Electric Vehicles in Urban Traffic Networks : Application to the Grenoble Metropolitan Area,” in 21st European Control Conference (ECC), Bucharest, Romania, 2023.

[11] M. U. B. Niazi, C. Canudas-de-Wit, A. Kibangou, and P.-A. Bliman, “Optimal Control of Urban Human Mobility for Epidemic Mitigation,” in 2021 60th IEEE Conference on Decision and Control (CDC). IEEE, dec 2021, pp. 6958–6963.

[12] C. Fiori, K. Ahn, and H. A. Rakha, “Power-based electric vehicle energy consumption model: Model development and validation,” Applied Energy, vol. 168, pp. 257–268, 2016.

[13] U. Pratap, C. Canudas-de-Wit, and F. Garin, “Where, when and how people move in large-scale urban networks : the Grenoble saga,” Univ. Grenoble Alpes, CNRS, INRIA, Grenoble INP, Gipsa-lab, Tech. Rep., 2022.

[14] H. Barbosa, M. Barthelemy, G. Ghoshal, C. R. James, M. Lenormand, T. Louail, R. Menezes, J. J. Ramasco, F. Simini, and M. Tomasini, “Human mobility: Models and applications,” Physics Reports, vol. 734, pp. 1–74, 2018.

[15] H. Bouscasse, I. Joly, and J. Peyhardi, “A new family of qualitative choice models: An application of reference models to travel mode choice,” Transportation Research Part B: Methodological, vol. 121, pp. 74–91, 2019.