# Multiple-scale analysis of the simplest large-delay differential equation

Gregory Kozyreff  Optique Nonlinéaire Théorique, Université libre de Bruxelles
(U.L.B.), CP 231, Belgium

###### Abstract

A delayed term in a differential equation reflects the fact that information
takes significant time to travel from one place to another within a process
being studied. Despite de apparent similarity with ordinary differential
equations, delay-differential equations (DDE) are known to be fundamentally
different and to require a dedicate mathematical apparatus for their analysis.
Indeed, when the delay is large, it was found that they can sometimes be
related to spatially extended dynamical systems. The purpose of this paper is
to explain this fact in the simplest possible DDE by way of a multiple-scale
analysis. We show the asymptotic correspondence of that linear DDE with the
diffusion equation. This partial differential equations arises from a
solvability condition that differs from the ones usually encountered in
textbooks on asymptotics: In the limit of large delays, the leading-order
problem is a map and secular divergence at subsequent orders stem from forcing
terms in that map.

##  I Introduction

In mathematical modelling, to be able to describe a physical, chemical, or an
automated process by a lumped-element model, i.e., by a finite set of time-
dependent variables, is a promising starting point for fruitful analysis.
Occasionally, the interaction between elements of the model takes place with a
time delay that cannot be neglected. In such a case, the appropriate
mathematical problem to be solved typically consists of one ore more Delay
Differential Equations (DDEs), also called functional differential equations,
as opposed to ordinary ones. The consequences can be dramatic: first order
DDEs can have periodic or even chaotic solutions [1, 2, 3, 4, 5, 6, 7]. Hence,
despite their similarity in writing, DDEs are fundamentally different from
ODEs and require specific consideration. It is easy to observe, for instance,
that linear, constant-coefficient DDEs generally admit an infinity of
independent exponential solutions. They may thus be regarded as ODEs of
infinite order [8]. In 1996, Giacommelli and Politi went one step further by
showing the asymptotic equivalence of some DDEs to partial differential
equations when the delay is large [9, 10]. This result followed earlier
numerical simulations [11] pointing to the relevance of a 2D representation of
the solutions of large-delay DDEs, see Figure 1(a). This “spatio-temporal
equivalence” has been confirmed in several subsequent works [12, 13, 14] and
the purpose of this paper is to discuss it in the simplest possible setting,
namely the equation

$\displaystyle y^{\prime}(t)+y(t)$ $\displaystyle=ry(t-T),$ $\displaystyle T$
$\displaystyle\gg 1,$ (1)

with initial data $y(t)=\psi(t)$ in the range $t\in[-T,0]$. The interpretation
of this equation is straightforward: it models a linear system whose internal
dynamics is governed by the left hand side, with a simple exponential decay,
and which is subjected to a delayed feedback with strength $r$. Here, $T$ is
the delay normalised by the eigenvalue of the isolated system. Equation 1
appears in recent developments as the deterministic part of a Langevin
equation

$y^{\prime}(t)+y(t)=ry(t-T)+\mu\xi(t),$ (2)

where $\xi(t)$ is a noisy forcing term. The latter can model Brownian motion
with a memory effect [15, 16] or delayed control [17]. More recently, it was
proposed to study chaotic diffusion mediated by a nonlinear DDE [18]. Equation
1 is also a special case of

$y^{\prime}(t)=ay(t)+by(t-T)$ (3)

which has been studied in detail as one of the simplest DDE [19, 20, 21].

Concretely, we will derive from eq. 1 the multiple-scale asymptotic
approximation

$y\sim r^{s}e^{-\ln(r)^{2}z/2}\times Y(s,z),$ (4)

with $s\sim t/(T+1)$, $z=t/T^{3}$ and

$\displaystyle\frac{\partial Y}{\partial z}$
$\displaystyle=\frac{1}{2}\frac{\partial^{2}Y}{\partial s^{2}},$
$\displaystyle Y(s,z)$ $\displaystyle=Y(s-1,z).$ (5)

This calculation serves three pedagogical purposes:

  * •

to highlight the general interest of investigating the large-$T$ limit in
DDEs,

  * •

to detail the inner workings of the method of multiple scales in that
particular framework,

  * •

to provide instructors with a new, easily workable application of that method,
beyond the usual perturbed harmonic oscillator, in the spirit of [22].

###  I.1 Application of the multiple-scale approximation

Before deriving eqs. 4 and 5, let us demonstrate their usefulness as a tool to
analyze eq. 1. Recall first that an initial condition $\delta(s+s_{0})$ at
$z=z_{0}$ evolves under the above diffusion equation as [23]

$Y(s,z)=\frac{e^{-(s+s_{0})^{2}/2(z-z_{0})}}{\sqrt{2\pi(z-z_{0})}}.$ (6)

From this, we may deduce, when $r=1$, that a peaked initial condition will
recur periodically with a flattening profile and a peak value that decreases
in time as $1/\sqrt{t+c}$, where $c$ is an appropriate constant. Indeed, a
gaussian initial condition

$e^{-k(t+t_{0})^{2}}$

in eq. 1 translates into $e^{-k(T+1)^{2}(s+s_{0})^{2}}$ in eq. 5. It may thus
be assimilated, up to an appropriate factor, to a delta function in the $s$
variable in the limit $T\to\infty$. By the same token, any similarly localized
initial condition can be asymptotically regarded as proportional to a delta
function on the $s$-scale. Therefore, eq. 6 applies. Furthermore, $t=-t_{0}$
corresponds to $z_{0}=-t_{0}/T^{3}$. Hence, evaluating eq. 6 at $s=-s_{0}$, we
deduce the following envelope for the maxima of $y(t)$:

$\frac{\text{const}}{\sqrt{z-z_{0}}}\propto\frac{1}{\sqrt{t+t_{0}}}.$ (7)

When $r\neq 1$ the exponential factor in eq. 4 is applied. One thus obtains
the envelope

$\text{const}\times\frac{r^{t\left(1/(T+1)+\ln(r)/T^{3}\right)}}{\sqrt{t+t_{0}}}.$
(8)

We compare the above formulas with numerical solutions of eq. 1 with $r=1$ and
$r=1.1$ respectively in figs. 1 and 2. fig. 1(a) confirms the relevance of the
“spatio-temporal” representation of the solution, in which $\mod(t,T+1)\propto
s$, $\lfloor x\rfloor$ denotes the integer part of $x$ and $\lfloor
t/(T+1)\rfloor$ is a discrete variable on which the evolution is so slow that
it is asymptotically equivalent to $z$. Next, fig. 1(b) and fig. 2 show the
quantitative agreement of eqs. 7 and 8 with the numerical simulations.
Especially striking is the non-monotonous behavior seen in the latter, which
conforms to the intuition brought by eq. 5. Indeed, diffusion promotes an
initial flattening and attenuation of the peaks, before the amplification
factor $r>1$ takes over.

Note that eq. 4 and eq. 6 actually make up an asymptotic approximation of the
Green function of eq. 1, in an alternative way to the exact solution given in
[15]. This Green function is used to study the noisy extension eq. 2 of eq. 1.

Figure 1: Numerical solution of eq. 1 with $r=1$, $T=30$, and initial
condition $y(t)=20e^{-(t+25)^{2}}$ in the range $-T<t<0$. (a) Pseudo spatio-
temporal plot, demonstrating the $T+1$ period of recurrence and the diffusive
behavior. (b) Usual representation as a function of $t$. Figure 2: Numerical
solution of eq. 1 with $r=1.1$, $T=30$, and initial condition
$y(t)=20e^{-(t+T/2)^{2}}$ in the range $-T<t<0$. The orange curve is the
envelope of the peaks predicted by eq. 8.

The rest of the paper is devoted to the derivation of eqs. 4 and 5. It starts
in section II with a standard linear analysis of eq. 1. The characteristic
equation is derived and then simplified in the large-$T$ limit. This gives us
insight on the relevant time scales for the multiple scale analysis, which we
develop in section III. We carry out the calculation twice, as the method can
be implemented in two slightly different but equally instructive ways. For the
sake of simplicity, we initially focus on $r=1$ and deal with general $r$ as a
later step. Finally, we present our conclusions in section IV.

##  II Linear spectrum

Let us look for a solution of (1) in the form of $y(t)=\exp\lambda t$.
Assuming $r=1$ for simplicity, this yields the characteristic equation

$\lambda+1=e^{-\lambda T}.$ (9)

Rearranging terms, we have

$\displaystyle\lambda$ $\displaystyle=-1+\frac{u}{T},$ $\displaystyle ue^{u}$
$\displaystyle=Te^{T}.$ (10)

The equation $ye^{y}=x$ possesses a discrete infinity of complex solutions,
denoted $W_{n}(x)$ and called “ProductLog[n,x]” in Mathematica. $W_{n}$ is the
$n^{\text{th}}$ branch of the complex Lambert function [24]. Hence, the
spectrum of eq. 1 is, exactly,

$\lambda_{n}=-1+\frac{W_{n}\left(Te^{T}\right)}{T}.$ (11)

Figure 3: The spectrum of eq. 1 for $T=100$ (blue) compared to the
approximation eq. 15 (orange circles). The separation between roots along the
imaginary axis asymptotes to $2\pi/(T+1)$.

Equation 11 is useful to draw the spectrum with a symbolic software but not
very enlightening to anyone who is not an expert of the Lambert function.
Fortunately, we can make significant progress in the limit $T\to\infty$.
Writing $\lambda T=\sigma+i\omega$, eq. 9 yields

$\displaystyle\frac{\sigma}{T}+1-e^{-\sigma}\cos\left(\omega\right)$
$\displaystyle=0,$
$\displaystyle\frac{\omega}{T}+e^{-\sigma}\sin\left(\omega\right)$
$\displaystyle=0.$ (12)

This suggests the expansions
$\sigma\sim\sigma_{0}+T^{-1}\sigma_{1}+T^{-2}\sigma_{2}+\dotsb$ and
$\omega\sim\omega_{0}+T^{-1}\omega_{1}+\dotsb$. At leading order, we get

$\displaystyle 1-e^{-\sigma_{0}}\cos\left(\omega_{0}\right)$
$\displaystyle=0,$ $\displaystyle e^{-\sigma_{0}}\sin\left(\omega_{0}\right)$
$\displaystyle=0,$ (13)

which implies $\sigma_{0}=0$ and $\omega_{0}=2n\pi$, where $n$ is an integer.
It is important to remark that $n$ must be assumed to be $O(1)$, so that
$\omega/T$ can be treated as $O(1/T)$ in eq. 12. To carry out the calculation
to higher orders presents no difficulty and can be proposed as an exercise.
One obtains

$\displaystyle\sigma_{1}$ $\displaystyle=0,$ $\displaystyle\omega_{1}$
$\displaystyle=-\omega_{0},$ $\displaystyle\sigma_{2}$
$\displaystyle=\frac{-\omega_{0}^{2}}{2},$ $\displaystyle\omega_{2}$
$\displaystyle=\omega_{0},$ $\displaystyle\dotsb$ (14)

Eventually,

$\lambda\sim\frac{1}{T}\left[-\frac{2n^{2}\pi^{2}}{T^{2}}+2in\pi\left(1-\frac{1}{T}+\frac{1}{T^{2}}\right)+\dotsb\right]\sim\frac{1}{T}\left(-\frac{2n^{2}\pi^{2}}{T^{2}}+\frac{2in\pi}{1+1/T}\right)=-\frac{2n^{2}\pi^{2}}{T^{3}}+\frac{2in\pi}{T+1}.$
(15)

The exact and approximate spectra are depicted in fig. 3 and are found to be
in good agreement as long as the imaginary part is small, i.e for mode numbers
satisfying $2\pi n\ll T$.

On the other hand, note that the diffusion equation eq. 5 has exponential
solutions $\exp(\kappa z+i\Lambda s)$ provided that

$\kappa=-\Lambda^{2}/2,$ (16)

while the periodic boundary condition in $s$ imposes $\Lambda=2n\pi$, with
integer $n$. Using the definitions of $s$ and $z$, one thus obtains
$\exp\\{[-2n^{2}\pi^{2}/T^{3}+2in\pi/(T+1)]t\\}$, which is consistent with eq.
15.

The imaginary part of $\lambda$ points to an oscillatory evolution with
approximate period $T+1$. Notice that time in eq. 1 is rescaled to the
intrinsic time scale of the isolated ($r=0$) dynamical system. With a more
general time unit, this intrinsic time scale would numerically be given by
$\bar{t}_{i}$, the delay by $\bar{T}$, and the period of damped oscillations
by $\bar{T}+\bar{t}_{i}$. In view of this, the periodicity is fixed by the
time required for information to be fed back into the dynamical system plus
the time required to internally process it. Such a conclusion would be
difficult to draw from the contemplation of eq. 11 alone.

The analysis for $r\neq 1$ follows the same pattern. We leave it as an
exercise to show that, in the general case,

$\lambda\sim\frac{1}{T}\left[-\frac{\ln
r}{2T^{2}}-\frac{2n^{2}\pi^{2}}{T^{2}}+\left(\ln
r+2in\pi\right)\left(1-\frac{1}{T}+\frac{1+\ln
r}{T^{2}}\right)+\dotsb\right].$ (17)

##  III Multiple-scale analysis

The linear stability analysis with $r=1$ reveals that the time scales for
oscillations and damping become infinitely separated as $T\to\infty$, which
suggests a multiple-scale analysis. In this section, we propose two slightly
distinct implementations of the method. The first one is a two-time expansion
of the solution,

$y(t)\sim Y(s,z),$ (18)

in which the fast time is a strained coordinate:

$s=\left(1+\frac{a_{1}}{T}+\frac{a_{2}}{T^{2}}+\dotsb\right)\frac{t}{T},$ (19)

and $z=t/T^{3}$. In that implementation, we construct a solution whose
periodicity is strictly $1$ in $s$. The constants $a_{1}$ and $a_{2}$ are
determined in the course of resolution, but they can in fact be guessed from
the results of the linear stability analysis.

The alternative is to pose a three-time ansatz

$y(t)\sim Y(\tau,\eta,z),$ (20)

where $\tau=1/T$ and $\eta=t/T^{2}$ with no preconception of how $f$ must
depend on each time scale. Both approaches have their advantages and
disadvantages. To introduce multiple time scales and convert a total
derivative into a partial differential operator can be a deterring prospect to
students who are exposed to the method for the first time. This speaks in
favor of minimizing the number of time scales and, hence, of choosing the
ansatz (18). On the other hand, the strained $s$ coordinate may appear pulled
out of a hat. This can make (20) more easily acceptable to those who prefer a
step-by-step approach. Eventually, given the lightness of the calculations, it
may be useful for the instructor to present the two variants of the
calculation and, in so doing, to demonstrate the robustness of the method. In
order to keep the calculation down to its essential details, we focus first on
the case $r=1$. We treat the case $r\neq 1$ as a later step.

###  III.1 Two-time calculation ($r=1$)

Assuming that $y$ is asymptotically given by the ansatz (18), differentiation
with respect to $t$ yields, by the chain rule,

$y^{\prime}(t)\sim\left(1+\frac{a_{1}}{T}+\frac{a_{2}}{T^{2}}+\dotsb\right)\frac{1}{T}\frac{\partial
Y}{\partial s}+\frac{1}{T^{3}}\frac{\partial Y}{\partial z}.$ (21)

On the other hand, the delayed term becomes

$y(t-T)\sim
Y\left(s-1-\frac{a_{1}}{T}-\frac{a_{2}}{T^{2}}-\dotsb,z-\frac{1}{T^{2}}\right)\\\
\sim Y\left(s-1,z\right)-\frac{a_{1}}{T}\frac{\partial}{\partial
s}Y\left(s-1,z\right)+\frac{1}{T^{2}}\left[\frac{a_{1}^{2}}{2}\frac{\partial^{2}}{\partial
s^{2}}-a_{2}\frac{\partial}{\partial s}-\frac{\partial}{\partial
z}\right]Y\left(s-1,z\right)+\dotsb$ (22)

Hence, eq. 1 is transformed into

$\left(\frac{1}{T}+\frac{a_{1}}{T^{2}}\right)\frac{\partial Y(s,z)}{\partial
s}+Y(s,z)\sim Y\left(s-1,z\right)-\frac{a_{1}}{T}\frac{\partial}{\partial
s}Y\left(s-1,z\right)\\\
+\frac{1}{T^{2}}\left[\frac{a_{1}^{2}}{2}\frac{\partial^{2}}{\partial
s^{2}}-a_{2}\frac{\partial}{\partial s}-\frac{\partial}{\partial
z}\right]Y\left(s-1,z\right)+O\left(\frac{1}{T^{3}}\right).$ (23)

Expanding $Y$ as $Y_{0}+T^{-1}Y_{1}+T^{-2}Y_{2}+\dotsb$, one obtains, at
$O(1)$,

$Y_{0}(s,z)=Y_{0}\left(s-1,z\right).$ (24)

so that $Y_{0}$ is periodic with period $1$ in $s$. Next, collecting $O(1/T)$
terms and using eq. 24, we find

$Y_{1}(s,z)=Y_{1}\left(s-1,z\right)-\left(1+a_{1}\right)\frac{\partial}{\partial
s}Y_{0}(s,z).$ (25)

The forcing term, being periodic, causes a secular divergence of $f_{1}$:

$Y_{1}(s+j,z)=Y_{1}\left(s,z\right)-j\left(1+a_{1}\right)\frac{\partial}{\partial
s}Y_{0}(s,z).$ (26)

Hence, the asymptotic ordering $T^{-1}Y_{1}\ll Y_{0}$ breaks down for
$j=O(T)$, unless we set

$a_{1}=-1.$ (27)

At $O(T^{-2})$, we then find

$Y_{2}(s,z)=Y_{2}\left(s-1,z\right)+\left(1-a_{2}\right)\frac{\partial
Y_{0}}{\partial s}+\frac{1}{2}\frac{\partial^{2}Y_{0}}{\partial
s^{2}}-\frac{\partial Y_{0}}{\partial z},$ (28)

where $Y_{0}$ is evaluated at $(s,z)$. Following the same reasoning as at the
previous order, the solvability condition is

$\frac{\partial Y_{0}}{\partial z}=\left(1-a_{2}\right)\frac{\partial
Y_{0}}{\partial s}+\frac{1}{2}\frac{\partial^{2}Y_{0}}{\partial s^{2}},$ (29)

At this stage, it may seem that the constant $a_{2}$ is still free. However,
if we insist that the periodic component of the solution is strictly $1$ in
$s$, then $Y_{0}$ is a combination of functions of the form $\exp(\kappa
z+2in\pi s)$, where, according to eq. 29,
$\kappa=2in\pi(1-a_{2})-2n^{2}\pi^{2}$. Now, the imaginary part of $\kappa$
perturbs the period of the solution and to avoid this, we make it vanish:

$a_{2}=1,$ (30)

which agrees with eq. 15. We thus obtained eq. 5.

###  III.2 Three-time calculation

If we assume the multi-time ansatz eq. 20, then

$y^{\prime}(t)\sim\frac{1}{T}\frac{\partial
Y}{\partial\tau}+\frac{1}{T^{2}}\frac{\partial
Y}{\partial\eta}+\frac{1}{T^{3}}\frac{\partial Y}{\partial z}$ (31)

and

$y(t-T)\sim Y\left(\tau-1,\eta-T^{-1},z-T^{-2}\right)\\\ \sim
Y\left(\tau-1,\eta,z\right)-\frac{1}{T}\frac{\partial}{\partial\eta}Y\left(\tau-1,\eta,z\right)-\frac{1}{T^{2}}\left(\frac{\partial}{\partial
z}-\frac{1}{2}\frac{\partial^{2}}{\partial\eta^{2}}\right)Y\left(\tau-1,\eta,z\right).$
(32)

We thus get to solve

$\left(1+\frac{1}{T}\frac{\partial}{\partial\tau}+\frac{1}{T^{2}}\frac{\partial}{\partial\eta}+\frac{1}{T^{3}}\frac{\partial}{\partial
z}\right)Y\left(\tau,\eta,z\right)=\left(1-\frac{1}{T}\frac{\partial}{\partial\eta}-\frac{1}{T^{2}}\frac{\partial}{\partial
z}+\frac{1}{2T^{2}}\frac{\partial}{\partial\eta^{2}}+\dotsb\right)Y\left(\tau-1,\eta,z\right).$
(33)

Expanding $Y$ again as $Y_{0}+T^{-1}Y_{1}+\dotsb$, we obtain

$Y_{0}(\tau,\eta,z)=Y_{0}(\tau-1,\eta,z),$ (34)

i.e., the same map as before, with $s$ replaced by $\tau$. At $O(T^{-1})$, we
have

$Y_{1}(\tau,\eta,z)-Y_{1}(\tau-1,\eta,z)=-\frac{\partial}{\partial\tau}Y_{0}(\tau,\eta,z)-\frac{\partial}{\partial\eta}Y_{0}(\tau-1,\eta,z)=-\left(\frac{\partial}{\partial\tau}+\frac{\partial}{\partial\eta}\right)Y_{0}(\tau,\eta,z).$
(35)

Since the right hand side is a period-1 function in $\tau$, solvability
requires that it vanishes:

$\frac{\partial Y_{0}}{\partial\tau}+\frac{\partial Y_{0}}{\partial\eta}=0.$
(36)

The general solution is, simply,

$Y_{0}(\tau,\eta,z)=Y_{0}(\tau-\eta,z),$ (37)

or, equivalenty,

$\displaystyle Y_{0}(\tau,\eta,z)$ $\displaystyle=Y_{0}(S,z),$ $\displaystyle
S=\tau-\eta=\frac{t}{T}\left(1-\frac{1}{T}\right).$ (38)

The equation for $Y_{1}$ then implies that it is periodic in $\tau$. Next, at
$O(T^{-2})$, and taking into account the periodicity of $Y_{0}$ and $Y_{1}$,
the problem for $Y_{2}$ is

$Y_{2}(\tau,\eta,z)-Y_{2}(\tau-1,\eta,z)=-\left(\frac{\partial}{\partial\tau}+\frac{\partial}{\partial\eta}\right)Y_{1}(\tau,\eta,z)+\left(\frac{\partial}{\partial
S}+\frac{1}{2}\frac{\partial}{\partial S^{2}}-\frac{\partial}{\partial
z}\right)Y_{0}(S,z).$ (39)

Secular divergence of $Y_{2}$ is avoided by making the right hand side vanish:

$\frac{\partial Y_{1}}{\partial\tau}+\frac{\partial
Y_{1}}{\partial\eta}=\frac{\partial Y_{0}}{\partial
S}+\frac{1}{2}\frac{\partial^{2}Y_{0}}{\partial S^{2}}-\frac{\partial
Y_{0}}{\partial z}.$ (40)

This is an equation for $Y_{1}$, with general solution

$Y_{1}(\tau,\eta,z)=Y_{1}(S,z)+\eta\left(\frac{\partial Y_{0}}{\partial
S}+\frac{1}{2}\frac{\partial^{2}Y_{0}}{\partial S^{2}}-\frac{\partial
Y_{0}}{\partial z}\right).$ (41)

There is therefore a secular divergence in $\eta$ of $Y_{1}$ and to avoid
this, we have a solvability condition on the solvability condition:

$\frac{\partial Y_{0}}{\partial z}-\frac{\partial Y_{0}}{\partial
S}=\frac{1}{2}\frac{\partial^{2}Y_{0}}{\partial S^{2}}.$ (42)

Finally, if we write

$\displaystyle Y_{0}(S,z)$ $\displaystyle=\phi(s,z),$ $\displaystyle
s=S+z=\tau-\eta+z=\frac{t}{T}\left(1-\frac{1}{T}+\frac{1}{T^{2}}\right),$ (43)

we obtain

$\frac{\partial\phi}{\partial z}=\frac{1}{2}\frac{\partial^{2}\phi}{\partial
s^{2}},$ (44)

i.e. the desired result. Note that the fact that $Y_{0}(\tau,\eta,z)$ is of
period $1$ in $\tau$ means that $\phi$ is of period $1$ in $s$, in full
consistency with what precedes.

Remark One may be tempted to purely and simply set $Y_{1}=0$, after noting (i)
that the initial problem is linear and (ii) that eq. 35 is homogenous after
applying the solvability condition. Indeed, this would expedite the derivation
of eq. 44. However, it would conceal the structure of nested solvability
conditions appearing at $O(T^{-2})$. While there is nothing wrong in seeking
the shortest route to the answer, especially in the eyes of an applied
mathematician, one should bear in mind that $Y_{1}$ should, in all generality,
be retained. This correction may be required to properly handle $O(T^{-1})$
terms in the initial condition or weak nonlinearities. Interestingly, while
there is only a single equation to solve at $O(T^{-2})$, more than one
solvability conditions can be extracted from it: the main one, eq. 40, and the
secondary one, eq. 42.

###  III.3 The case $r\neq 1$

We now revise the calculation in the more general case $r\neq 1$. Here again,
the two implementations yield the same result. We limit ourselves to the “two-
time” calculation and assume

$y\sim f(s,z),$ (45)

where $s$ is given by eq. 19. This time, eq. 1 is transformed into

$\left(\frac{1}{T}+\frac{a_{1}}{T^{2}}\right)\frac{\partial f(s,z)}{\partial
s}+f(s,z)\sim
r\bigg{\\{}f\left(s-1,z\right)-\frac{a_{1}}{T}\frac{\partial}{\partial
s}f\left(s-1,z\right)\\\
+\frac{1}{T^{2}}\left[\frac{a_{1}^{2}}{2}\frac{\partial^{2}}{\partial
s^{2}}-a_{2}\frac{\partial}{\partial s}-\frac{\partial}{\partial
z}\right]f\left(s-1,z\right)\bigg{\\}}+O\left(\frac{1}{T^{3}}\right).$ (46)

Expanding $f$ as $f_{0}+T^{-1}f_{1}+T^{-2}f_{2}+\dotsb$, the $O(1)$ problem is
now

$f_{0}(s,z)=rf_{0}\left(s-1,z\right).$ (47)

Posing $f_{0}(s,z)=r^{s}F_{0}(s,z)$, this yields

$F_{0}(s,z)=F_{0}\left(s-1,z\right),$ (48)

so that $F_{0}$ is periodic with period $1$ in $s$. Next, collecting $O(1/T)$
terms and using eqs. 47 and 48, we find

$f_{1}(s,z)=rf_{1}\left(s-1,z\right)-\left(1+a_{1}\right)r^{s}\left[\ln(r)F_{0}(s,z)+\frac{\partial}{\partial
s}F_{0}(s,z)\right].$ (49)

The terms between bracket, being periodic, cause a secular divergence of
$f_{1}$. Indeed, letting $f_{1}(s,z)=r^{s}F_{1}(s,z)$, one rapidly finds that

$F_{1}(s+j,z)=F_{1}\left(s,z\right)-\left(1+a_{1}\right)j\left[\ln(r)F_{0}(s,z)+\frac{\partial}{\partial
s}F_{0}(s,z)\right].$ (50)

Hence, irrespective of the factor $r^{s}$, $T^{-1}f_{1}$ ceases to be small
compared to $f_{0}$ when $j=O(T)$. To prevent this, we thus set

$a_{1}=-1.$ (51)

At $O(T^{-2})$, we have

$f_{2}(s,z)=rf_{2}\left(s-1,z\right)+r^{s}\Bigg{\\{}\left[\frac{\ln(r)^{2}}{2}+\ln(r)\left(1-a_{2}\right)\right]F_{0}(s,z)+\left(1+\ln(r)-a_{2}\right)\frac{\partial
F_{0}}{\partial s}+\frac{1}{2}\frac{\partial^{2}F_{0}}{\partial
s^{2}}-\frac{\partial F_{0}}{\partial z}\Bigg{\\}}.$ (52)

The solvability condition is now

$\frac{\partial F_{0}}{\partial z}=l_{0}F_{0}(s,z)+l_{1}\frac{\partial
F_{0}}{\partial s}+\frac{1}{2}\frac{\partial^{2}F_{0}}{\partial s^{2}},$ (53)

where $l_{0}=\frac{\ln(r)^{2}}{2}+\ln(r)\left(1-a_{2}\right)$ and
$l_{1}=1+\ln(r)-a_{2}$. The term multiplied by $l_{1}$ causes a change in the
periodicity of the solution and is therefore set to zero:

$a_{2}=1+\ln(r).$ (54)

This agrees with eq. 17. Hence, $l_{0}=\frac{-\ln(r)^{2}}{2}$. It only remains
to make a small change of variable, namely to write $F_{0}=e^{l_{0}z}Y(s,z)$
to obtain LABEL:{intro:diffusion}.

##  IV Discussion

The derivation presented in this paper hints to the great generality of the
diffusion equation eq. 5 as the linear backbone of large-delay differential
equations. Not all DDE develop the multiple-scale structure presented here
(see below) but when they do, diffusion is to be expected from the Taylor
expansion of the delayed term with the appropriate scales: see the $O(T^{-2})$
terms in the right hand side of eq. 32.

An interpretation of the delay in a DDE like eq. 1 is that the output of a
given system undergoes some kind of propagation before being fed back. The
form of the delayed term, $ry(t-T)$, is strongly suggestive of a hyperbolic
PDE as a mediator of this feedback, a point of view that was emphasized in
[25]. Why, then, should a parabolic PDE such as eq. 5 arise out of a
hyperbolic one? The answer to this question lies in the fact that the system
subjected to feedback is itself dissipative, being described in isolation by
the left hand side of eq. 1. Loosely speaking, a complete feedback cycle
includes a dispersion-less propagation of duration $T$ followed by an
attenuation of unit duration. In the spectral domain, the left hand side of
eq. 1 is a low-pass filter. This is where information is degraded, as happens
in diffusive processes. In this regard, one should bear in mind that,
numerically, $T$ is the ratio of the delay to the internal dissipation time of
the system. As we wrote earlier, in a general unit system, the delay is
numerically given by $\bar{T}$ and the internal dynamics is characterized by
$\bar{t}_{i}$, with $T=\bar{T}/\bar{t}_{i}$. It follows that the parameter $T$
does not appear as a result of the delay alone. It exists thanks to both the
delay _and_ the internal time scale. Even if the latter is short compared to
the former, the dissipation process occurring on the internal timescale is not
to be neglected. Another facet of this question is given by the picture of the
linear spectrum in fig. 3. Here, the plot of the real part of the eigenvalues
as a function of their imaginary part is analogous to dispersion relations
found in pattern formation, where the exponential growth rate of a
perturbation is plotted as a function of its wave number [26, 27, 28]. In that
context, the existence of a small band of wave numbers with near-zero, maximum
growth rate generically leads to diffusive terms in amplitude equations.
Presently, the spectrum is discrete but tends to a continuum as $T\to\infty$.
Hence, with a trained eye, diffusion can directly be anticipated at the sight
of fig. 3.

A general technical feature exemplified by the present calculation is that the
leading order problem is a map, see eqs. 24, 34 and 47. This is because the
main time scale is asymptotically set by the delay, and time derivatives of
functions evolving on this time scale are $O(T^{-1})$. In problems where the
period is not asymptotically given by the delay, a map is _not_ expected as
the leading order problem, even if the delay is large. Consider for example
Minorsky’s equation [29]

$\displaystyle y^{\prime\prime}(t)+\epsilon y^{\prime}(t)+\Omega^{2}y(t)$
$\displaystyle=-by^{\prime}(t-T)+\epsilon cy^{\prime}(t-T)^{3},$
$\displaystyle\epsilon$ $\displaystyle\ll 1,$ $\displaystyle
T=O\left(\epsilon^{-2}\right),$ (55)

which displays a Hopf bifurcations with an $O(1)$ frequency. Here $\Omega$
sets the oscillation period and the leading-order problem of the multiple-
scale analysis is the familiar harmonic oscillator. The linear spectrum is
again densified by the largeness of $T$ but it now displays a maximum with
near vanishing $\mbox{Re}(\lambda)$ in the vicinity of
$\mbox{Im}(\lambda)=\Omega$. This is akin to a Turing bifurcation in spatially
extended dynamical system. Following the preceding discussion, diffusion is
therefore expected in the amplitude equation. Taking nonlinear terms into
account one eventually obtains a Ginzburg-Landau equation (see [30] for a
detailed derivation.) In any case, whether the leading order problem is a map
or a harmonic oscillator, the general procedural idea of the multiple scale
analysis is the same: identify and kill secular divergences.

The present paper may give the reader the false impression that multiple-scale
analysis is the method of choice to study all DDE in the large-delay limit.
This is not the case and the following equation is a simple counter-example:

$-y^{\prime}(t)+y(t)=y(t-T)+y^{3}(t).$ (56)

Compared to eq. 1, we have simply changed the sign of the time derivate and
added a nonlinear term to avoid blow up. A numerical simulation is shown on
fig. 4. Independently of the initial condition, the solution of this equation
asymptotes to sustained square wave oscillations of period $2T$ that display
abrupt switching between approximately $\sqrt{2}$ and $-\sqrt{2}$. This
dynamical behavior is obviously not compatible with eq. 5. To study such a
limit cycle, the method of matched asymptotic expansions appears more
appropriate [31, 32].

Figure 4: Numerical solution of eq. 56 with $T=10$ and
$y(t)=0.1\sin(t)-0.02\cos(3t)$ for $t<0$ as an example of a non-diffusive
long-term behavior.

###### Acknowledgements.

G.K. is supported by the Fonds de la Recherche Scientifique - FNRS (Belgium.)

## References

  * Mackey and Glass [1977] M. C. Mackey and L. Glass, Oscillation and chaos in physiological control systems, Science 197, 287 (1977).
  * Ikeda [1979] K. Ikeda, Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system, Optics Communications 30, 257 (1979).
  * Farmer [1982] J. D. Farmer, Chaotic attractors of an infinite-dimensional dynamical system, Physica D: Nonlinear Phenomena 4, 366 (1982).
  * Erneux  _et al._ [2004] T. Erneux, L. Larger, M. W. Lee, and J.-P. Goedgebuer, Ikeda Hopf bifurcation revisited, Physica D: Nonlinear Phenomena 194, 49 (2004).
  * Sprott [2007] J. Sprott, A simple chaotic delay differential equation, Physics Letters A 366, 397 (2007).
  * Müller  _et al._ [2018] D. Müller, A. Otto, and G. Radons, Laminar chaos, Physical Review Letters 120, 10.1103/physrevlett.120.084102 (2018).
  * Banerjee and Biswas [2018] T. Banerjee and D. Biswas,  _Time-Delayed Chaotic Dynamical Systems_ (Springer International Publishing, 2018).
  * Bellman  _et al._ [1963] R. Bellman, K. L. Cooke, _et al._ ,  _Differential-difference equations_ (Academic press, 1963).
  * Giacomelli and Politi [1996] G. Giacomelli and A. Politi, Relationship between delayed and spatially extended dynamical systems, Physical Review Letters 76, 2686 (1996).
  * Giacomelli and Politi [1998] G. Giacomelli and A. Politi, Multiple scale analysis of delayed dynamical systems, Physica D: Nonlinear Phenomena 117, 26 (1998).
  * Arecchi  _et al._ [1992] F. T. Arecchi, G. Giacomelli, A. Lapucci, and R. Meucci, Two-dimensional representation of a delayed dynamical system, Physical Review A 45, R4225 (1992).
  * Wolfrum and Yanchuk [2006] M. Wolfrum and S. Yanchuk, Eckhaus instability in systems with large delay, Physical Review Letters 96, 10.1103/physrevlett.96.220201 (2006).
  * Wolfrum  _et al._ [2010] M. Wolfrum, S. Yanchuk, P. Hövel, and E. Schöll, Complex dynamics in delay-differential equations with large delay, The European Physical Journal Special Topics 191, 91 (2010).
  * Yanchuk and Giacomelli [2017] S. Yanchuk and G. Giacomelli, Spatio-temporal phenomena in complex systems with time delays, Journal of Physics A: Mathematical and Theoretical 50, 103001 (2017).
  * Budini and Cáceres [2004] A. A. Budini and M. O. Cáceres, Functional characterization of linear delay langevin equations, Physical Review E 70, 10.1103/physreve.70.046104 (2004).
  * Munakata  _et al._ [2009] T. Munakata, S. Iwama, and M. Kimizuka, Linear stochastic system with delay: Energy balance and entropy production, Physical Review E 79, 10.1103/physreve.79.031104 (2009).
  * Ando  _et al._ [2017] H. Ando, K. Takehara, and M. U. Kobayashi, Time-delayed feedback control of diffusion in random walkers, Physical review E 96, 012148 (2017).
  * Albers  _et al._ [2022] T. Albers, D. Müller-Bender, L. Hille, and G. Radons, Chaotic diffusion in delay systems: Giant enhancement by time lag modulation, Physical Review Letters 128, 10.1103/physrevlett.128.074101 (2022).
  * Hayes [1950] N. D. Hayes, Roots of the transcendental equation associated with a certain difference-differential equation, Journal of the London Mathematical Society s1-25, 226 (1950).
  * Stepan [1989] G. Stepan,  _Retarded Dynamical Systems: Stability and Characteristic Functions_ , Pitman research notes in mathematics series (Longman Scientific & Technical, 1989).
  * Michiels and Niculescu [2007] W. Michiels and S.-I. Niculescu,  _Stability and Stabilization of Time-Delay Systems_ (Society for Industrial and Applied Mathematics, 2007).
  * Edwards [2000] D. A. Edwards, An alternative example of the method of multiple scales, SIAM Review 42, 317 (2000).
  * Strauss [2007] W. A. Strauss,  _Partial differential equations: An introduction_ (John Wiley & Sons, 2007).
  * [24] Wikipedia, Lambert W function.
  * Bünner  _et al._ [1997] M. J. Bünner, T. Meyer, A. Kittel, and J. Par, On the correspondence of time-delay and spatially extended systems, Zeitschrift für Naturforschung A 52, 573 (1997).
  * Cross and Hohenberg [1993] M. C. Cross and P. C. Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. 65, 851 (1993).
  * Hoyle [2006] R. Hoyle,  _Pattern formation : an introduction to methods_ (Cambridge University, 2006).
  * Pismen [2006] L. Pismen,  _Patterns and interfaces in dissipative dynamics_ (Springer, 2006).
  * Pinney [1958] E. Pinney,  _Ordinary difference-differential equations_ (Univ of California Press, 1958).
  * Erneux [2009] T. Erneux,  _Applied delay differential equations_, Vol. 3 (Springer Science & Business Media, 2009).
  * Fowler [1982] A. C. Fowler, An asymptotic analysis of the delayed logistic equation when the delay is large, IMA Journal of Applied Mathematics 28, 41 (1982).
  * Nizette [2004] M. Nizette, Stability of square oscillations in a delayed-feedback system, Physical Review E 70, 10.1103/physreve.70.056204 (2004).


# Cross-Entropy-Based Approach to Multi-Objective Electric Vehicle Charging
Infrastructure Planning ††thanks: †These authors contributed equally.
††thanks: *Corresponding author: Hao Wang. ††thanks: This work was supported
in part by the Australian Research Council (ARC) Discovery Early Career
Researcher Award (DECRA) under Grant DE230100046.

Jinhao Li1†, Yu Hui Yuan1†, Qiushi Cui2, Hao Wang1,3∗ 1Department of Data
Science and AI, Faculty of IT, Monash University, Australia
2School of Electrical Engineering, Chongqing University, Chongqing, China
3Monash Energy Institute, Monash University, Australia
Emails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Pure electric vehicles (PEVs) are increasingly adopted to decarbonize the
transport sector and mitigate global warming. However, the inadequate PEV
charging infrastructure may hinder the further adoption of PEVs in the large-
scale traffic network, which calls for effective planning solutions for the
charging station (CS) placement. The deployment of charging infrastructure
inevitably increases the load on the associated power distribution network.
Therefore, we are motivated to develop a comprehensive multi-objective
framework for optimal CS placement in a traffic network overlaid by a
distribution network, considering multiple stakeholders’ interested factors,
such as traffic flow, PEV charging time cost, PEV travel distance, and the
reliability of the distribution network. We leverage a cross-entropy-based
method to solve the optimal CS placement and evaluate our method in a real-
world 183-node traffic network in Chengdu, China, overlaid by a 26-region
distribution network. It is demonstrated that our work provides various viable
planning options favoring different objectives for the stakeholders’ decision-
making in practice.

###### Index Terms:

Electric vehicle, charging infrastructure planning, multi-objective
optimization, cross-entropy method.

##  I Introduction

The substantial increase in greenhouse gas emissions in the past decades has
led to a steady rise in global temperature, creating worldwide awareness of
the adverse impacts of fossil fuel usage [1]. In particular, energy-related
carbon emissions generated by conventional combustion energy vehicles
contribute to approximately $15\%$ of global greenhouse gas emissions [2].
Mitigation efforts for global warming must be accelerated more urgently and
rapidly. Hence, promoting the adoption of pure electric vehicles (PEVs) has
been widely recognized as the most promising and effective solution, aligning
with the goal of net-zero transition [3].

Large-scale PEV deployment faces significant challenges due to the lack of
publicly available charging infrastructures, even though the government is
providing both financial and policy incentives. The placement problem of
charging stations (CS) in the traffic network is complex and influenced by
multi-variate factors, including the traffic throughput, the geospatial
location of candidate sites, and the accessibility of charging services for
PEV owners [4]. Moreover, the CS placement is closely related to the interests
of major stakeholders, e.g., PEV owners and the traffic network operator.
Specifically, PEV owners, being the largest group of stakeholders, prioritize
accessible charging services within the shortest time in the occurrence of the
charging demand as the ideal PEV charging infrastructure planning. The traffic
network operator prefers the CS placement at traffic nodes with higher traffic
throughput, such that the charging infrastructure is fully utilized. In
addition to the stakeholders’ interests, the impact of charging demand on the
power distribution network brought by associated with the CS placement must be
considered, as increasing charge demand can lead to significant voltage
deviations and threaten the reliability of the distribution network [5].

PEV charging infrastructure planning has drawn increasing attention in recent
studies, but is still in the early stage of development. Masoum et al. [6]
proposed a smart load management approach for coordinating PEV chargers in
distribution feeders, but neglected the impacts of the traffic network.
Hajimiragha et al. [7] employed a robust optimization method to primarily
address the positive environmental effects of PEV adoption in reducing carbon
emissions. Jia et al. in [8] facilitated CS placement by minimizing integrated
investment and operation costs of the charging infrastructure. Yao et al. [9]
developed a user equilibrium-based traffic assignment model, which maximizes
captured traffic flow under energy-related cost constraints. The developed
model requires specific information of PEV owners, but such information is
difficult to obtain in real practice. Wu et al. [10] proposed a stochastic
flow-capturing local model for optimal CS placement in the traffic network.

While the above studies tend to focus on either the traffic network or the
distribution network, recent works by [11, 12] proposed flow-based models
along with both network constraints to determine the optimal locations of
charging infrastructure. However, these approaches did not consider the
interests of PEV owners, as their primary focus was on maximizing the charging
support to the overall traffic flow in the traffic network. Improving the
accessibility of charging services for PEV owners is often overlooked but must
be considered to facilitate further adoption of PEVs for transportation
electrification.

As discussed above, existing studies lack a comprehensive framework covering
all important factors, e.g., traffic flow, the accessibility of charging
services, and the reliability of the distribution network, for the optimal CS
placement for charging infrastructure planning. We are motivated to bridge
such a research gap by formulating the PEV charging infrastructure planning as
a multi-objective optimization problem, including traffic-flow-oriented
charging support maximization, PEV charging time cost minimization, PEV travel
distance minimization, and distribution network reliability. We use the cross-
entropy method, known for its robustness and fast convergence [13], to derive
the optimal CS placement solution. Importantly, unlike previous studies that
usually provided a single optimal solution for charging infrastructure
planning, our work offers various planning options, as stakeholders often need
a set of solutions to tradeoff multiple objectives or factors in their
decision-making. The main contributions of our work are summarized as follows.

  * •

Multi-Objective Modeling in the Traffic Network Coupled with Distribution
Network: We model the PEV charging infrastructure planning in a traffic
network coupled with a distribution network as a multi-objective optimization
problem, considering multiple essential factors, such as traffic flow, PEV
charging time cost, PEV travel distance, and the reliability of the
distribution network.

  * •

Cross-Entropy-Based Solution Method and Real-World Simulations: We leverage
the cross-entropy method to solve the multi-objective CS placement problem,
which is validated in the real-world $183$-node traffic network located in
Chengdu, China, overlaid by a $26$-region distribution network.

  * •

Providing Multiple Viable Solutions in Practice: Our method offers multiple
viable CS placement options based on different preferences of the
stakeholders, allowing stakeholders to explicitly tradeoff multiple objectives
and facilitate better decision making in real practice.

The remainder of this paper is organized as follows. We formulate the multi-
objective PEV charging infrastructure planning model in Section II. Section
III introduces the cross-entropy solution method for solving the optimal CS
placement problem. A case study based on real-world data is presented in
Section IV. Section V concludes this paper.

##  II System Model

We show the system model for PEV charging infrastructure planning in Fig. 1.
We consider a traffic network with a set of traffic nodes denoted by
$\mathcal{N}$. The traffic network is overlaid by a multi-region distribution
network, and we denote the region set as $\mathcal{R}$. Our optimal PEV
charging infrastructure planning seeks to deploy a set of CSs, denoted by
$\mathcal{S}$, at traffic nodes, aiming to optimize multiple objectives, e.g.,
maximize charging support to traffic flow, reduce PEV charging time, minimize
PEV travel distance, and mitigate the effects of CS placement on the
reliability of the distribution network. We present the detailed models for
each objective from Section II-A to II-D, respectively.

###  II-A Maximizing Charging Support to Traffic Flow

One of the key objectives of the PEV charging infrastructure planning is to
fulfill PEVs’ charging demand, such that PEVs’ travel needs can be well
satisfied without worrying about the lack of charging infrastructure. Since
charging activities are more likely to occur at traffic nodes carrying
significant traffic flows [14], our infrastructure planning needs to, to the
utmost extent, consider the traffic flow and place CSs at the nodes with high
traffic. Specifically, we define the traffic flow at the $n$-th traffic node
as $f_{n}$, representing the total number of trips that both originate from
and end at the $n$-th node. The objective of the traffic-flow-oriented
charging support maximization problem can be formulated as

$J^{\text{flow}}=\sum_{n\in\mathcal{N}}\mathbb{I}_{n}f_{n},\\\ $ (1)

where $\mathbb{I}_{n}=\\{0,1\\}$ is the indicator variable to decide whether a
CS should be placed at $n$-th traffic node.

Figure 1: The system model of PEV charging infrastructure planning.

###  II-B PEV Charging Time Cost Minimization

In addition to providing maximum charging support to the overall traffic flow,
effective PEV charging infrastructure planning also needs to consider the
interests of PEV owners, address the known range anxiety barrier of charging
services, and ensure high accessibility of the charging sites [11]. As a
result, when deploying CSs, careful attention should be paid to the time
needed for charging services, including the PEV traveling time, queuing time,
and actual charging time.

####  II-B1 PEV Traveling Time

Given the limited battery capacity of the PEV, it is essential to place CSs at
the most appropriate nodes to ensure convenience for the PEV owners with the
shortest traveling time. The PEV traveling time for each traffic node can be
estimated using its average traveling time to all other nodes in the traffic
network [15]. We let $v$ denote the constant for average travel speed and
$d_{nm}$ denote the geospatial distance between the $n$-th and the $m$-th
nodes. To monitor the actual charging activity, traffic congestion is also
considered by incorporating the normalized traffic flow between the $m$-th and
the $n$-th nodes, i.e., $\hat{f}_{n}+\hat{f}_{m}$. The normalized traffic flow
is defined as $\hat{f}_{n}=f_{n}/F_{\text{max}}$, where $F_{\text{max}}$ is
the maximum traffic flow in the traffic network. The PEV traveling time for
the $n$-th traffic node can be defined as

$t^{\text{travel}}_{n}=\frac{1}{v\left(|\mathcal{N}|-1\right)}\sum_{m\in\mathcal{N}\backslash\\{n\\}}d_{mn}\left(\hat{f}_{n}+\hat{f}_{m}\right).$
(2)

####  II-B2 PEV Queuing Time

For a CS with multiple chargers, high charging demand at a CS-placed traffic
node may lead to queuing. To model this, we assume that the service system of
a CS follows a multiple-server queuing system with Poisson arrival and
exponential service times [16]. Specifically, the arrival of the PEVs is
described as a Poisson process and the interval time of arriving follows the
negative exponent distribution. The average arrival rate of PEVs at the $n$-th
traffic node can be estimated using historical data, which can be formulated
as

$\lambda_{n}=f_{n}\tau\omega\frac{\int_{t}^{t+\Delta t}f_{n}}{\Delta
t\sum_{m\in\mathcal{N}}\int_{t}^{t+\Delta t}f_{m}}\int_{t}^{t+\Delta
t}\hat{f}_{n},$ (3)

where $\tau$ is the charging frequency of a PEV, $\omega$ is the average ratio
of daily PEV charging [11], and $\Delta t$ is the historical time horizon. The
utilization of a CS, also known as the proportion of time that the CS is busy,
can be further calculated as

$\rho_{n}=\frac{\lambda_{n}}{C\mu},$ (4)

where $\mu$ is the average service rate of the charging device and $C$ is the
number of chargers equipped at a CS. The average number of PEVs in the queue
at the $n$-th node can be calculated using metrics defined in Eq. (3) and (4)
and is formulated as

$L^{\text{queue}}_{n}=\frac{\rho_{n}p_{n}^{\text{ava}}\left(\frac{\lambda_{n}}{\mu}\right)^{C}}{C!\left(1-\rho_{n}\right)^{2}},$
(5)

with the probability of the CS being available $p_{n}^{\text{ava}}$ formulated
as

$p_{n}^{\text{ava}}=1/\left[\sum_{c=0}^{C-1}\frac{\left(C\rho_{n}\right)^{c}}{c!}+\frac{\left(C\rho_{n}\right)^{C}}{C!\left(1-\rho_{n}\right)}\right].$
(6)

Combining the average arrival rate and the average waiting number of PEVs, we
can define the average queuing time at the $n$-th traffic node (if placed with
a CS) as

$t_{n}^{\text{queue}}=\frac{L_{n}^{\text{queue}}}{\lambda_{n}}.$ (7)

####  II-B3 PEV Charging Event Time

The charging event time here refers to the time required for a PEV to be fully
charged when connected to an available charger at a CS. We let
$SoC^{\text{ch}}$ denote the PEV’s state of charge (SoC) at the start of the
charging service, $E^{\text{PEV}}_{\text{max}}$ denote the storage capacity of
the PEV, $\eta^{\text{ch}}$ denote the charging efficiency of the PEV battery,
and $P_{\text{max}}^{\text{CS}}$ denote the rated power output of the utilized
charging device. The charging event time at a CS-placed traffic node can be
formulated as

$t_{n}^{\text{ch}}=\frac{\left(SoC_{\text{max}}-SoC^{\text{ch}}\right)E^{\text{PEV}}_{\text{max}}}{P^{\text{CS}}_{\text{max}}\eta^{\text{ch}}},$
(8)

where $SoC_{\text{max}}$ is the maximum allowed SoC of the PEV.

Integrating the PEV traveling time, queuing time, and charging event time, we
formulate the objective of the PEV charging time cost minimization as

$J^{\text{ch}}=\sum_{n\in\mathcal{N}}\mathbb{I}_{n}\left(t_{n}^{\text{travel}}+t_{n}^{\text{queue}}+t_{n}^{\text{ch}}\right).$
(9)

###  II-C PEV Travel Distance Minimization

Effective PEV charging infrastructure planning must consider the accessibility
of candidate sites for PEV owners. Specifically, the PEV travel distance to
CS-placed traffic nodes should be minimized to ensure high accessibility of
charging services. We let $\hat{f}_{mn}$ denote the normalized traffic flow
captured on the route from the $m$-th node to the $n$-th node. The objective
of PEV travel distance minimization can be formulated as

$J^{\text{dis}}=\sum_{m\in\mathcal{N}}\sum_{n\in\mathcal{N}\backslash\\{m\\}}\mathbb{I}_{n}d_{mn}\hat{f}_{mn}.$
(10)

###  II-D Optimizing Distribution Network Reliability

The CS deployment in the traffic network can increase the load in the coupled
distribution network, causing voltage deviations and threatening the
reliability of the distribution network [5]. In this work, we investigate how
to mitigate such negative effects of CS placement on a region-based
distribution network overlaid over the traffic network. We assume that CSs
within each region share the same substation for providing charging services.
We let $r$ denote one specific region of the distribution network and
$n^{\text{CS}}_{r}$ denote the number of traffic nodes selected as candidate
CS sites in the $r$-th region of the distribution network. To avoid excess CS
placement in one certain region of the distribution network, we introduce a
penalty coefficient $\pi$ with the aim of alleviating each region’s load
burden and ensuring the reliability of distribution networks. The objective of
minimizing the effect of CS placement on the distribution network is
formulated as

$J^{\text{DN}}=\pi\sum_{r\in\mathcal{R}}n_{r}^{\text{CS}}\sum_{n\in
r\cap\mathcal{N}}\frac{\left(SoC_{\text{max}}-SoC^{\text{ch}}\right)E_{\text{max}}^{\text{PEV}}}{\eta^{\text{ch}}}\hat{f}_{n}.$
(11)

The advantage of our model is that, instead of acquiring the topology and
parameters of the distribution network, we estimate the effect of CS placement
on each CS-placed traffic node using PEV charging demand and nodal-based
traffic flow as an effective approximation.

###  II-E Multi-Objective Optimization Formulation of PEV Charging
Infrastructure Planning

As presented from Section II-A to II-D, the PEV infrastructure planning
problem has four objectives: traffic-flow-oriented charging support
maximization, PEV charging time cost minimization, PEV travel distance
minimization, and optimizing distribution network reliability. To solve this
multi-objective optimization problem, we assign different weight coefficients
to each objective, denoted by $\alpha^{\text{flow}}$, $\alpha^{\text{ch}}$,
$\alpha^{\text{dis}}$, and $\alpha^{\text{DN}}$, respectively. This helps
transform the multi-objective optimization into a single-objective
optimization problem formulated as

$\displaystyle\min\hskip 2.5pt$
$\displaystyle\alpha^{\text{flow}}\bar{J}^{\text{flow}}+\alpha^{\text{ch}}\bar{J}^{\text{ch}}+\alpha^{\text{dis}}\bar{J}^{\text{dis}}+\alpha^{\text{DN}}\bar{J}^{\text{DN}}$
(12) s.t. $\displaystyle\sum_{n\in\mathcal{N}}\mathbb{I}_{n}\leq K,$ (13)
$\displaystyle\qquad\mathbb{I}_{n}=\\{0,1\\},\quad\forall n\in\mathcal{N},$
(14)

where objectives $J^{\text{flow}}$, $J^{\text{ch}}$, $J^{\text{dis}}$, and
$J^{\text{DN}}$ are normalized into range $[0,1]$ and then denoted by
$\bar{J}^{\text{flow}}$, $\bar{J}^{\text{ch}}$, $\bar{J}^{\text{dis}}$, and
$\bar{J}^{\text{DN}}$, respectively. In particular, since the objective
$J^{\text{flow}}$ aims to maximize charging support with respect to the
traffic flow, we normalize it unlike the other three minimized objectives,
which can be formulated as

$\bar{J}^{\text{flow}}=\frac{J^{\text{flow}}_{\text{max}}-J^{\text{flow}}}{J^{\text{flow}}_{\text{max}}-J^{\text{flow}}_{\text{min}}}.$
(15)

The constraint in Eq. (13) indicates that the number of CSs planned in the
traffic network cannot exceed the budget of the deployed CSs, which is denoted
by $K$. Eq. (14) describes the decision variable of the CS placement problem,
i.e., $\mathbb{I}_{n}$ to determine whether to place a CS at the $n$-th
traffic node.

##  III Cross-Entropy-Based Solution Method

To solve the PEV charging infrastructure planning problem formulated in
Section II-E, we propose a cross-entropy solution method that balances
different objectives to achieve optimal CS placement. Compared to classical
and evolutionary algorithms, the cross-entropy approach has several
advantages, including fast convergence, strong robustness, insensitivity to
initialized points, and most importantly, better interpretability in
optimization results [13].

We denote the objective value of the planning problem defined in Eq. (12) as
$\gamma=J(\mathcal{S})$. To minimize the objective $\gamma$, the cross-entropy
method first randomizes a family of probability distribution functions (PDFs)
for the traffic nodes, indicating the probability of being placed with a CS.
These PDFs are denoted as $f(\mathcal{S};\bm{p})$, where $\bm{p}$ is the
probability vector defined as $\bm{p}=[p_{1},p_{2},\cdots,p_{|\mathcal{N}|}]$.
The initialized PDFs are then updated in an adaptive manner to derive the
optimal probability for each traffic node to be placed with a CS.

Specifically, in the $t$-th iteration, we first set an adaptively updating
parameter denoted by $\gamma^{\text{ada}}_{t}$. The probability of a possible
planning scheme whose objective value is lower than such a threshold
$\gamma^{\text{ada}}_{t}$ can be formulated as

$l(\gamma)=P_{\bm{p}_{t}}\left\\{J\left(\mathcal{S}\right)\leq\gamma^{\text{ada}}_{t}\right\\}=\mathbb{E}_{\bm{p}_{t}}\mathbb{I}\left(J\left(\mathcal{S}\right)\leq\gamma^{\text{ada}}_{t}\right),$
(16)

where $\mathbb{E}_{\bm{p}_{t}}$ is the expectation of the PDFs in the $t$-th
iteration and $\mathbb{I}(J(\mathcal{S})\leq\gamma^{\text{ada}}_{t})$ is an
indicator variable that is equal to $1$ only if
$J(\mathcal{S})\leq\gamma^{\text{ada}}_{t}$. The primary objective of the
cross-entropy method is to minimize the cross entropy or variance of the
defined probability in Eq. (16). To achieve this, the cross-entropy method
updates the PDFs using the elite solutions, which are part of candidate
solutions with the lowest objective values, to guide the search towards the
global optimum in subsequent iterations [17]. Additionally, the objective
value of the worst solution among the elite solutions is set as the adaptively
updating parameter $\gamma_{t}^{\text{ada}}$ in each iteration.

We denote the candidate population size, i.e., the number of generated CS
placement plans, as $N^{\text{CE}}$. The proportion of elite solutions is
denoted by $\delta$, and the set of elite solutions is denoted by
$\bm{\mathcal{S}}^{\text{elite}}_{t}$. The updating process of the probability
parameter for the $n$-th traffic node in the $t$-th iteration can be
formulated as

$p_{t,n}=\frac{1}{\delta
N^{\text{CE}}}\sum_{\bm{\mathcal{S}}^{\text{elite}}_{t}}\mathbb{I}\left(J\left(\mathcal{S}\right)\leq\gamma_{t}^{\text{ada}}\right)\mathbb{I}_{n}.$
(17)

We set the maximum iteration number as $T^{\text{CE}}_{\text{max}}$ and assume
that the PDF in the cross-entropy method follows the Bernoulli distribution
$\mathcal{B}(\bm{p})$ with the success probability vector $\bm{p}$. The
detailed algorithmic procedure of our cross-entropy-based planning method is
presented in Algorithm. 1. By executing the cross-entropy method, we can
obtain the optimal solution for PEV infrastructure planning in the elite
solutions with the lowest objective value. More importantly, the iteration
process of the cross-entropy method clearly reveals how our algorithm
dynamically balances different objectives to achieve optimal planning,
providing insights in practice, which are analyzed in detail in Section IV.

Algorithm 1 The Cross-Entropy-Based Solution Method for PEV Charging
Infrastructure Planning

Initialize the probability vector $\bm{p}_{0}$.

for $t=1,\cdots,T^{\text{CE}}_{\text{max}}$ do

Generate $N^{\text{CE}}$ CS placement plans, where each node is selected based
on its corresponding Bernoulli distribution.

Calculate and rank the objective value $\gamma$ for all plans.

Update the adaptively updating parameter $\gamma_{t}^{\text{ada}}$.

Update the probability parameter for each node’s Bernoulli distribution
$\mathcal{B}(p_{t,n})$.

end for

##  IV Case Studies

###  IV-A Simulation Settings

We use realistic traffic flow data collected from a real-world $183$-node
traffic network located in Chengdu, China. The traffic network is coupled with
a simulated $26$-region distribution network. We illustrate traffic nodes
located in each region of the distribution network with a specific color, as
shown in Fig. 2. The time horizon of utilized traffic data is $7$ days. We
assume that the SoC of PEV when connected to a charger, denoted as
$SoC^{\text{ch}}$, follows a uniform distribution $\mathcal{U}(0.05,0.95)$. We
also set a stopping criterion for the cross-entropy method, specifying that
the iteration is terminated when the objective difference of the optimal
solutions in two consecutive iterations is smaller than $10^{-5}$. Our
algorithm is implemented using Python on an Intel i7-11800H 2.30 GHz/32 GB
laptop. The initialized parameters of the cross-entropy-based method are
presented in Table I.

Figure 2: The layout of the traffic network and distribution network in Chengdu. TABLE I: The initialized parameters. $|\mathcal{N}|$ | $128$ | $|\mathcal{R}|$ | $26$ | $v$ |  $30$ km/h
---|---|---|---|---|---
$\tau$ | $0.65$ | $\omega$ | $0.2$ | $C$ | $12$
$\mu$ | $40$ | $SoC_{\text{max}}$ | $0.95$ | $\eta^{\text{ch}}$ | $0.95$
$E_{\text{max}}^{\text{PEV}}$ |  $60$ KWh | $P_{\text{max}}^{\text{CS}}$ |  $60$ KW | $\pi$ | $0.5$
$K$ | $40$ | $\delta$ | $0.1$ | $N^{\text{CE}},T_{\text{max}}^{\text{CE}}$ | $100$

The weight coefficients of objectives, namely$\alpha^{\text{flow}}$,
$\alpha^{\text{ch}}$, $\alpha^{\text{dis}}$, and $\alpha^{\text{DN}}$, are all
set to a default value of $0.25$, indicating their equal importance. However,
in practice, the stakeholders, e.g., traffic network operator, distribution
network operator, and PEV owners, may have different preferences for specific
aspects of multi-objectives in the PEV charging infrastructure planning
decision-making. For example, different weights can be assigned to multi-
objectives, e.g., supporting more PEVs with charging demand, improving the
charging service with shorter charging time cost, or focusing on the
reliability of the distribution network. In our study, we analyze the optimal
CS placement based on stakeholders’ preferences, which are categorized into
six scenarios summarized in Table II. The default weight setting, i.e., $0.25$
for each objective, is the baseline Scenario for comparison. In contrast, we
assign larger weights, such as $0.5$ and $0.7$, to $\alpha^{\text{flow}}$,
$\alpha^{\text{ch}}$, and $\alpha^{\text{DN}}$, to prioritize the
corresponding objectives, respectively.

TABLE II: The objective weight settings in different simulation scenarios. Scenario Number |  Favored Objective | $\alpha^{\text{flow}}$ | $\alpha^{\text{ch}}$ | $\alpha^{\text{dis}}$ | $\alpha^{\text{DN}}$
---|---|---|---|---|---
$1$ | $J^{\text{flow}}$ | $\bm{0.5}$ | $0.5/3$ | $0.5/3$ | $0.5/3$
$2$ | $J^{\text{flow}}$ | $\bm{0.7}$ | $0.1$ | $0.1$ | $0.1$
$3$ | $J^{\text{ch}}$ | $0.5/3$ | $\bm{0.5}$ | $0.5/3$ | $0.5/3$
$4$ | $J^{\text{ch}}$ | $0.1$ | $\bm{0.7}$ | $0.1$ | $0.1$
$5$ | $J^{\text{DN}}$ | $0.5/3$ | $0.5/3$ | $0.5/3$ | $\bm{0.5}$
$6$ | $J^{\text{DN}}$ | $0.1$ | $0.1$ | $0.1$ | $\bm{0.7}$
Figure 3: The optimal CS placement in the baseline scenario. TABLE III: Simulation results derived from Scenario $1$ to $6$ and the baseline Scenario. Scenario | $\gamma_{\text{min}}$ | No. of CSs | Traffic Flow Supported | Average Charging Time |  Charging Demand in Distribution Network | Average Travel Distance
---|---|---|---|---|---|---
Baseline | $0.248$ | $21$ | $42.91\%$ |  $99$ mins |  $8,152$ MWh |  $38$ km
$1$ | $0.255$ | $37$ | $78.47\%$ |  $102$ mins |  $19,646$ MWh |  $30$ km
$2$ | $0.162$ | $39$ | $\bm{83.66\%}$ |  $127$ mins |  $26,312$ MWh |  $\bm{29}$ km
$3$ | $0.166$ | $13$ | $9.92\%$ |  $89$ mins |  $328$ MWh |  $46$ km
$4$ | $0.101$ | $20$ | $0.25\%$ |  $\bm{72}$ mins |  $824$ MWh |  $43$ km
$5$ | $0.167$ | $10$ | $0.03\%$ |  $129$ mins |  $\bm{1}$ MWh |  $47$ km
$6$ | $0.100$ | $10$ | $0.02\%$ |  $142$ mins |  $\bm{1}$ MWh |  $99$ km

(a) Scenario 1

(b) Scenario 2

(c) Scenario 3

(d) Scenario 4

(e) Scenario 5

(f) Scenario 6

Figure 4: The optimal CS placement solutions in Scenarios $1$ to $6$.

###  IV-B Charging Infrastructure Planning Solutions in Different Scenarios

The simulation results obtained in the aforementioned seven scenarios,
including the baseline scenario and additional six scenarios, are presented in
Table III. The optimal objective value $\gamma_{\text{min}}$, the number of
candidates for CS placement, the captured traffic flow, additional power
consumption for PEV charging services, and average PEV travel distance are
included in Table III. The corresponding optimal CS placement solutions in
seven scenarios are depicted in Fig. 3 (for the baseline scenario) and Fig. 4
(for Scenarios $1$ to $6$) for cross comparisons. The results in Table III
reveal that the optimal planning solutions from each scenario successfully
achieve their corresponding favored objectives, e.g., supporting more traffic
flow for Scenarios $1$ and $2$, while the baseline Scenario seems to balance
each objective for the optimal CS placement. The detailed analysis is
presented below.

####  IV-B1 Favoring Traffic-Flow-Oriented Charging Support Maximization

When favoring the traffic flow objective $J^{\text{flow}}$, our PEV
infrastructure planning tends to select more candidate sites ($37$ and $39$
CS-placed traffic nodes in Scenarios $1$ and $2$, respectively, as shown in
Table III) compared to other scenarios. The optimal planning in Scenarios $1$
and $2$ aims at supporting more PEVs with potential charging demand, i.e.,
$78.47\%$ and $83.66\%$ of traffic flow associated charging demand is well
supported in these two traffic-flow-favored scenarios, respectively.
Consequently, the selected candidate nodes tend to be densely clustered near
the city center of Chengdu (indicated as the black star in Fig. 4), as
depicted in Fig. 4a and 4b.

Moreover, by measuring the pair-wise geospatial distances between non-selected
traffic nodes and the candidate sites, we can observe that the traffic-flow-
favored planning solutions offer the minimum average distances among all
simulation scenarios, i.e., $30$ and $29$ kilometers for Scenarios $1$ and
$2$, respectively. However, such dense planning increases charging demand in
distribution network regions near the city center, with considerably
additional power consumption (more than $10,000$ MWh). Moreover, favoring
traffic flow associated charging demand support result in low accessibility
for potential PEV owners residing in outer city regions.

####  IV-B2 Favoring Charging Time Cost Minimization

Compared to other scenarios favoring $J^{\text{flow}}$ or $J^{\text{DN}}$, the
CS placement results in Scenarios $3$ and $4$ reveal that the candidate sites
tend to be evenly distributed in the traffic network, with relatively fewer
traffic nodes being selected, e.g., $20$ and $13$ candidates, respectively.
Such planning solutions, though capturing lower traffic flow, can
significantly reduce the PEV charging time, meeting the increasing charging
demand of PEV owners within a considerably shorter charging time cost ($89$
and $72$ minutes, respectively, as presented in Table III). As a result,
quality of service provided by the CS is improved.

####  IV-B3 Favoring Distribution Network Reliability

The optimal CS placement solutions in Scenarios $5$ and $6$ generate the
minimum additional power consumption among all simulated scenarios, with only
$10$ traffic nodes selected as candidate sites sparsely distributed outside
Chengdu City. However, such planning solutions present poor performance in
supporting traffic flow associated charging demand (e.g., $0.03\%$ and
$0.02\%$), leading to the average charging time of $129$ and $142$ minutes and
the average travel distance of $47$ and $99$ kilometers, in Scenarios $5$ and
$6$, respectively. These simulation results reveal that CS placement favoring
distribution network reliability is less beneficial for the majority of PEV
owners as well as the traffic network operator, indicating a clear interest
conflict among the stakeholders in PEV charging infrastructure planning.
Hence, achieving a trade-off between satisfying the increasing charging demand
and securing the distribution network is needed to balance the interests of
major stakeholders.

###  IV-C Trade-off among CS Planning Objectives

As discussed in Section IV-B, we see that compared to favoring one of the
objectives in the multi-objective CS placement problem, treating these
objectives with the same weight, i.e., the baseline scenario, appears to
strike an effective balance among the four objectives of the charging
infrastructure planning. As shown in Fig. 3, some of the selected candidates
are located around the city center while other CSs are placed relatively far
away from the city. Such a geospatial allocation pattern captures a fair
proportion of traffic flow on the one hand, while alleviating excessive load
stress on regions near the city center on the other hand.

It is worth noting that the main focus of our work is not to provide a single
optimal PEV charging infrastructure planning solution for a specific traffic
network coupled with a distribution network. Instead, our study aims to offer
a set of planning options considering multiple essential factors of the CS
placement and potentially conflicting interests in stakeholders’ decision-
making.

##  V Conclusion

In this paper, we proposed a cross-entropy-based method for optimal PEV
charging infrastructure planning. We formulated the CS placement problem in a
traffic network coupled with an electric distribution network as a multi-
objective optimization problem. We incorporated essential factors from the
perspectives of the traffic network operator, PEV owners, and the distribution
network operator, including traffic-flow-oriented charging support
maximization, PEV charging time cost minimization, PEV travel distance
minimization, and distribution network reliability optimization. The optimal
planning solutions were obtained using the cross-entropy algorithm. We
validated our method in a real-world $183$-node traffic network in Chengdu,
China, which is coupled with a $26$-region distribution network. Our method
offered various planning options for the stakeholders’ decision-making in
practice. The simulation results provided three noteworthy insights: 1)
overemphasizing the charging support for PEVs poses a challenge in managing
charging demand in regions with higher CS density for the distribution
network; 2) reducing the average charging time cost is crucial to enhance the
accessibility of charging infrastructures for PEV owners; 3) achieving a
trade-off should be reached to balance the interests of the distribution
network operator, traffic network operator, and PEV owners. Our work served a
comprehensive case study for solving the real-world PEV charging
infrastructure problem and highlighting the importance of effective CS
placement considering multiple stakeholders’ interests.

## References

  * [1] Z. Wu, M. Wang, J. Zheng, X. Sun, M. Zhao, and X. Wang, “Life cycle greenhouse gas emission reduction potential of battery electric vehicle,” _Journal of Cleaner Production_ , vol. 190, pp. 462–470, 2018.
  * [2] M. Metais, O. Jouini, Y. Perez, J. Berrada, and E. Suomalainen, “Too much or not enough? planning electric vehicle charging infrastructure: A review of modeling options,” _Renewable and Sustainable Energy Reviews_ , vol. 153, p. 111719, 2022.
  * [3] N. Melton, J. Axsen, and D. Sperling, “Moving beyond alternative fuel hype to decarbonize transportation,” _Nature Energy_ , vol. 1, no. 3, p. 16013, 2016\.
  * [4] M. Kavianipour, F. Fakhrmoosavi, H. Singh, M. Ghamami, A. Zockaie, Y. Ouyang, and R. Jackson, “Electric vehicle fast charging infrastructure planning in urban networks considering daily travel and charging behavior,” _Transportation Research Part D: Transport and Environment_ , vol. 93, p. 102769, 2021.
  * [5] S. Davidov and M. Pantoš, “Optimization model for charging infrastructure planning with electric power system reliability check,” _Energy_ , vol. 166, pp. 886–894, 2019.
  * [6] A. S. Masoum, S. Deilami, P. S. Moses, M. A. Masoum, and A. Abu-Siada, “Smart load management of plug-in electric vehicles in distribution and residential networks with charging stations for peak shaving and loss minimisation considering voltage regulation,” _IET generation, transmission & distribution_, vol. 5, no. 8, pp. 877–888, 2011.
  * [7] A. H. Hajimiragha, C. A. Canizares, M. W. Fowler, S. Moazeni, and A. Elkamel, “A robust optimization approach for planning the transition to plug-in hybrid electric vehicles,” _IEEE Transactions on Power Systems_ , vol. 26, no. 4, pp. 2264–2274, 2011.
  * [8] L. Jia, Z. Hu, Y. Song, and Z. Luo, “Optimal siting and sizing of electric vehicle charging stations,” in _2012 IEEE International Electric Vehicle Conference_ , 2012, pp. 1–6.
  * [9] W. Yao, J. Zhao, F. Wen, Z. Dong, Y. Xue, Y. Xu, and K. Meng, “A multi-objective collaborative planning strategy for integrated power distribution and electric vehicle charging systems,” _IEEE Transactions on Power Systems_ , vol. 29, no. 4, pp. 1811–1821, 2014.
  * [10] F. Wu and R. Sioshansi, “A stochastic flow-capturing model to optimize the location of fast-charging stations with uncertain electric vehicle flows,” _Transportation Research Part D: Transport and Environment_ , vol. 53, pp. 354–376, 2017.
  * [11] Y. Xiang, J. Liu, R. Li, F. Li, C. Gu, and S. Tang, “Economic planning of electric vehicle charging stations considering traffic constraints and load profile templates,” _Applied Energy_ , vol. 178, pp. 647–659, 2016.
  * [12] H. Zhang, S. J. Moura, Z. Hu, and Y. Song, “Pev fast-charging station siting and sizing on coupled transportation and power networks,” _IEEE Transactions on Smart Grid_ , vol. 9, no. 4, pp. 2595–2605, 2018.
  * [13] J. Bekker and C. Aldrich, “The cross-entropy method in multi-objective optimisation: An assessment,” _European Journal of Operational Research_ , vol. 211, no. 1, pp. 112–121, 2011.
  * [14] M. Daskin, “Network and discrete location: models, algorithms and applications,” _Journal of the Operational Research Society_ , vol. 48, no. 7, pp. 763–764, 1997.
  * [15] Q. Guo, S. Xin, H. Sun, Z. Li, and B. Zhang, “Rapid-charging navigation of electric vehicles based on real-time power systems and traffic data,” _IEEE Transactions on Smart Grid_ , vol. 5, no. 4, pp. 1969–1979, 2014. [Online]. Available: https://dx.doi.org/10.1109/tsg.2014.2309961
  * [16] Y. Zeng, L. Dong, and J. Ma, “Queuing modeling, analysis and simulation,” _Xi’an University of Electronic sciences and Technology Press_ , 2011.
  * [17] Z. I. Botev, D. P. Kroese, R. Y. Rubinstein, and P. L’Ecuyer, “Chapter 3 - the cross-entropy method for optimization,” in _Handbook of Statistics_ , ser. Handbook of Statistics, C. Rao and V. Govindaraju, Eds. Elsevier, 2013, vol. 31, pp. 35–59.