A Model Predictive Control Approach for
Low-Complexity Electric Vehicle Charging
Scheduling: Optimality and Scalability
Wanrong Tang Student Member, IEEE and Ying Jun (Angela) Zhang, Senior Member, IEEE

arXiv:1502.01456v3 [math.OC] 1 Apr 2016

Department of Information Engineering, The Chinese University of Hong Kong
Shatin, New Territories, Hong Kong
Abstract—With the increasing adoption of plug-in electric vehicles (PEVs), it is critical to develop efficient charging coordination
mechanisms that minimize the cost and impact of PEV integration to the power grid. In this paper, we consider the optimal PEV
charging scheduling, where the non-causal information about
future PEV arrivals is not known in advance, but its statistical
information can be estimated. This leads to an “online” charging
scheduling problem that is naturally formulated as a finitehorizon dynamic programming with continuous state space and
action space. To avoid the prohibitively high complexity of solving
such a dynamic programming problem, we provide a Model
Predictive Control (MPC) based algorithm with computational
complexity O(T 3 ), where T is the total number of time stages.
We rigorously analyze the performance gap between the nearoptimal solution of the MPC-based approach and the optimal
solution for any distributions of exogenous random variables.
Furthermore, our rigorous analysis shows that when the random
process describing the arrival of charging demands is first-order
periodic, the complexity of proposed algorithm can be reduced
to O(1), which is independent of T . Extensive simulations show
that the proposed online algorithm performs very closely to the
optimal online algorithm. The performance gap is smaller than
0.4% in most cases.

I. I NTRODUCTION
A. Background and Contributions
The massive deployment of PEVs imposes great challenges
to smart power grid, such as voltage deviation, increased power
losses, and higher peak load demands. It is critical to design
PEV charging mechanisms that minimize the cost and impact
of PEV integration. Previously, PEV charging coordination has
been extensively studied to minimize power loss, minimize
load variance, or minimize charging cost, etc [3]–[7]. Ideally,
the load demand can be flattened as much as possible if
the information about future charging demand is known noncausally when calculating the charging schedule. However, in
practice, a PEV charging station only knows the load demand
This work was presented in part as [1]. This work was supported in part
by the National Basic Research Program (973 program Program number
2013CB336701), and three grants from the Research Grants Council of
Hong Kong under General Research Funding (Project number 2150828 and
2150876) and Theme-Based Research Scheme (Project number T23-407/13N).
W. Tang is with the Department of Information Engineering, The Chinese
University of Hong Kong, Shatin, New Territories, Hong Kong (Email:
twr011@ie.cuhk.edu.hk).
Y. J. Zhang is with the Department of Information Engineering, The Chinese
University of Hong Kong. She is also with Shenzhen Research Institute, The
Chinese University of Hong Kong (Email: yjzhang@ie.cuhk.edu.hk).

of the PEVs that have arrived, but not that of the PEVs
coming in the future. Fortunately, the statistical information
of the future charging demands can often be acquired through
historic data, which benefits the control of the PEV charging
scheduling in practical scenarios.
In this paper, we consider the optimal PEV charging
scheduling, assuming that the future charging demand is not
known a priori, but its statistical information can be estimated.
In particular, we define the cost of PEV charging as a general
strictly convex increasing function of the instantaneous load
demand. Minimizing such a cost leads to a flattened load demand, which is highly desirable for many reasons [3]–[7]. The
online PEV charging scheduling problem is formulated as a
finite-horizon dynamic programming problem with continuous
state space and action space. To avoid the prohibitively high
complexity of solving such a dynamic programming problem,
we provide a Model Predictive Control (MPC) approach to
obtain a near-optimal solution. Instead of adopting the generic
convex optimization algorithms to solve the problem, we
propose an algorithm with computational complexity O(T 3 )
by exploring the load flattening feature of the solution, where
T is the total number of time stages. We rigorously analyze
the performance gap between the near-optimal solution of the
MPC-based approach and the optimal solution, and the result
applies to any distributions of exogenous random variables.
Specially, the performance gap is evaluated by the Value
of the Stochastic Solution (VSS), which represents the gap
between the solution of the approximate approach and that
of dynamic programming problem [8]–[10]. Furthermore, our
analysis shows that when the random process describing
the arrival of charging demands is first-order periodic, the
complexity of proposed algorithm can be reduced to O(1),
which is independent of T . Extensive simulations show that
the proposed algorithm performs very closely to the optimal
solution. The performance gap is smaller than 0.4% in most
cases. As such, the proposed algorithm is very appealing for
practical implementation due to its scalable computational
complexity and close to optimal performance.
The rest of the paper is organized as follows. A review of the
related work on the PEV charging scheduling with uncertain
load demand is presented in Section I-B. We introduce the
problem formulations of both offline and online PEV charging
problem in Section II. In Section III, we propose a MPC based
online algorithm and analyze its performance gap. The O(1)-

complexity algorithm is given when the arrival process is firstorder periodic in Section IV. Simulation results are presented
in Section V. Finally, the paper is concluded in Section VI.
B. Related Work
The works on the PEV charging scheduling with uncertain
PEV load demand include both simulation-based evaluations
[11], [12] and theoretical performance guarantees [13]–[17].
Meanwhile, MPC is one of most commonly approaches for
which has been widely adopted in recent studies [12]–[15].
[12] leverages the MPC based method to design a dynamic
charging and driving cost control scheme. Both [13] and [14]
apply MPC algorithms to minimize the load variation. [15]
proposes a plug and play MPC approach to minimize the
voltage fluctuations by assuming that the load demand is timeperiodic. Compared to [12]–[15], in this paper we analyze
the performance gap between the solution of MPC approach
and the optimal solution regardless of the distribution of the
load demand. Besides, we provide a more scalable algorithm
with O(1)− complexity as well as the optimality analysis
for the case when the load demand is first-order periodic.
Additionally, the objective functions in [12]–[15] are quadratic
forms of load demand. Whereas in this paper, the objective
function is a general strictly convex increasing function which
reflects both the charging cost and the load variance.
As to the amount of information needed, the EV charging
scheduling algorithms in [16] and [17] require the probability
distribution of the random PEV arrival process. In contrast,
the proposed algorithm in this paper only requires the firstorder moment, i.e., the expected values of the random demand
patterns. In practice, it is a lot easier to obtain the expected
values of random process than to obtain the probability distribution of a random process. Convex cost functions are also
considered in [18] and [19]. Both of them devise the online
algorithms for battery charging control problems, where there
is no charging deadline for the battery. The PEV charging
scheduling problem in this paper differs from stationary battery
charging in that each PEV has a demand to be satisfied before
a certain deadline.
II. P ROBLEM F ORMULATION
We consider the PEV charging scheduling problem, where
PEVs arrive at the charging station at random instants with
random charging demands that must be fulfilled before a
random departure time.
A. Optimal Offline PEV Charging Problem
For the ease of understanding, we first introduce an ideal
case, where all the non-causal information of based load and
the PEVs, including the arrival times, departure times and
charging demands are known to the charging station before
the system time. The entire system time is divided into T
equal-length time slots. Let N denote the set of PEVs that
arrive during the system time. Notice that for a given time
slot number T , N is itself a random set due to the random
arrival of PEVs. We denote by I(t) the set of PEVs that

(s)

are in the charging station during slot t. Denote by ti and
(e)
ti the arrival time slot and the departure time slot of PEV
i, respectively. di denotes the charging demand that PEV i
requires. The charging station needs to decide, the charging
rate xit , ∀i ∈ I(t). To satisfy the demand di , xit must satisfy
Pt(e)
i
(s) xit = di . Let st be the total charging rate of time slot
t=ti
t, i.e.,
X
xit , ∀t = 1, 2, · · · , T,

st =

(1)

i∈I(t)

which is also called charging load at time t. The total load
consists of both the charging load and the inelastic base load in
the same location. The base load, denoted by lt , represents the
load of other electricity consumptions at time t except
P for PEV
charging. Then, the total load at time t is given by i∈It xit +
lt . Suppose that the charging cost at time t is a strictly convex
increasing function of the total load, denoted by f (st + lt ).
The convexity and increasing property of f (st + lt ) reflects
the fact that each unit of additional power demand becomes
more expensive to obtain and make available to the consumer.
For example, in the wholesale market, the instantaneous cost
can be modeled as an increasing quadratic function of the
instant load [4]–[6]. On the other hand, the convexity of f (st +
lt ) also captures the intent of reducing the load fluctuation
over time [7]. Then the total cost over time T is computed
PT
(s) (e)
as t=1 f (st + lt ). In the ideal case, assume that lt , ti , ti ,
and di for all t = 1, · · · , T, i ∈ N are known non-causally at
the beginning of the system time. Then, the charging station
can solve (2) and obtain the optimal charging rate, denoted by
x∗it for all time t and the optimal total cost, denoted by Ψ1 .
Such a solution is referred to as an “optimal offline solution”.
T
X

Ψ1 = min
xit


f


X

xit + lt 

(2a)

xit = di , ∀i ∈ N ,

(2b)

t=1

i∈I(t)

(e)

ti
X

s.t.

(s)
t=ti

(e)

(s)

xit ≥ 0, ∀t = ti , · · · , ti , ∀i ∈ N .

(2c)

In particular, the optimal
rate, denoted by
P total charging
∗
s∗t , is defined as s∗t =
x
.
Note
that there are in
i∈I(t) it
total O(T |I(t)|) variables in (2), where |I(t)| denotes the
cardinality of the set I(t). This number can be quite large
when the number of cars present at each time slot, |I(t)|, is
large. Next, we propose an equivalent transformation of (2)
that drastically reduces the number of variables. In particular,
the following Theorem 1 shows that as long as we find the
optimal s∗t ∀t, the optimal x∗it ∀i, t can be obtained by earliest
deadline first (EDF) scheduling.
Theorem 1: If a set of st ’s satisfy the following inequality
for all n = 1, · · · , T
n
X
t=1

X
(e)

i∈{i|ti

di ≤
=t}

n
X
t=1

st ≤

n
X
t=1

X
(s)

i∈{i|ti

di ,

(3)

=t}

then there exists at least a set of xit ’s that is feasible to (2).
One such set of xit ’s can be obtained by EDF scheduling,
which charges the PEV i ∈ I(t) with the earliest deadline

at a rate st at each time t. Moreover, when st = s∗t , the set
of xit ’s obtained by EDF scheduling are the optimal solution,
x∗it , to (2).
Proof: Please see the detailed proof in Appendix A. To
see Theorem 1, note that (3) implies that the total energy
charged by any time slot n is no less than the total charging
demand that must be satisfied by time n and no more then
the total charging demand of PEVs which have arrived up
to time n. On the other hand, by EDF scheduling, PEVs
with earlier deadlines must be fully charged before those
with later deadlines can be charged. Thus, (3) guarantees
the fulfillment of the charging demands of each individual
PEV. With Theorem 1, we can transform (2) to the following
equivalent problem with T variables.
Ψ1 = min
st

s.t.

T
X
t=1
n
X

f (st + lt )
st ≥

t=1
n
X

n
X
j=1

st ≤

t=1

n
X
j=1

statistical information of the unknown future demands. The
charging rate sk , once determined, cannot be changed in the
future. Specifically, the remaining charging
Pk−1 demand of PEV
i at time k is given by dˆki = di − t=t(s) xit . Note that,
i
dˆki = di for all PEVs that have not yet arrived by time k − 1.
A close look at (4) suggests that the charging schedule st
only depends on the total charging demand that needs to be
finished before a certain time, but not the demand due to
individual
P PEVs. Thus, for notational simplicity, we define
d˜kt = i∈{i|t(e) =t} dˆki , ∀t = k, · · · , T, as the total unfinished
i
charging demand at time k that must be completed by time t.
With this, we define the state of system at time t as
Dt = [lt , d˜tt , d˜tt+1 , · · · , d˜tT ],

(4a)

(5)

d˜tt0

X
(e)

i∈{i|ti

di , ∀n = 1, · · · , T. (4b)
=j}

X

where lt is the base load at time t,
is the total unfinished
charging demand at time t that must be completed by time
t0 . Let ξ t represent the random arrival events at time t. ξ t is
defined as
t
ξt = [ιt , ηtt , ηt+1
, · · · , ηet t ],

di , ∀n = 1, · · · , T. (4c)

(s)
i∈{i|ti =j}

s∗t

to (4) has an interesting feature: it
The optimal solution
does not change with the cost function f (st + lt ), as long as
f is strictly convex. Moreover, s∗t also minimizes the variance
of total load subjecting to (4b) and (4c), where
the variance
PT
PT
t=1 st +lt 2
of total load is defined as t=1 (st + lt −
) [13],
T
[14]. This is proved in Theorem 2.
Theorem 2: The optimal solution s∗t to (4) does not change
with the cost function f (.), as long as f (.) is strictly convex.
Moreover, s∗t is a load flattening solution that minimizes the
variance of total load.
Proof: Please see the detailed proof in Appendix B.
Remark 1: In practice, a capacity constraint on st + lt is
present for each t due to the hardware limitations and security
concerns. The constraint is omitted in our formulation for
the following reason. Theorem 2 indicates that the optimal
solution s∗t to (4) minimizes the variance of total load. That is,
any other scheduling solution would have a higher peak load,
and therefore is more likely to violate the capacity constraint.
In this sense, the optimal solution s∗t to (4) is “capacity
optimal” in the sense that if the optimal solution to Problem
4 (or equivalently Problem 2) violates the capacity constraint,
then there does not exist any other scheduling solutions that
satisfy the capacity constraint.
B. Online PEV Charging problem
For the online PEV charging problem, the charging schedule only depends on the statistic information of future load
demand, the current based load and the remaining charging
demands and deadlines of the PEVs that have arrived so far.
In contrast to the ideal case in the last subsection, in practice
the charging station only knows the remaining charging demands and departure deadlines of the PEVs that have already
arrived, as well as the current base load. The online charging
scheduling algorithm computes the charging rate sk at each
time slot k based on the known causal information and the

(6)

where ιt is the base load at time t, ηtt0 is the total charging
demand that arrive at time t and must be fulfilled by time t0 ,
et is the latest deadline among the PEVs that arrive at time t.
Then, the state transition, defined as
Dt+1 := g(st , Dt , ξt+1 ),

(7)

is calculated as follows:
lt+1 = ιt+1

(8)

and




= d˜tt0 − st −
d˜t+1
t0

0
tX
−1

+ +
, ∀t0 = t+1, · · · , T. (9)
d˜tj   +ηtt+1
0

j=t

+

Here, [x] = max{x, 0}. With the above definitions of system
state and state transition, we are now ready to rewrite (4) into
the following finite-horizon dynamic programming problem.
Qk (Dk ) = min
sk

f (sk + lk ) + Eξk+1 [Qk+1 (g(sk , Dk , ξk+1 ))]
(10a)

s. t.

d˜kk ≤ sk ≤

T
X

d˜kt ,

(10b)

t=k

where Qk (Dk ) is the optimal value of the dynamic programming at time k. The left side of (10b) ensures all charging
demands to be satisfied before their deadlines. The right side
of (10b) implies that the total charging power up to a certain
time cannot exceed the total demands that have arrived up to
that time. By slight abuse of notation, in the rest of the paper
we denote the optimal solutions to both the online and offline
problems as s∗k , when no confusion arises. The actual meaning
of s∗k will be clear from the context. Suppose that s∗k is the
optimal solution to (10) at stage k. Then, the total cost at the
end of system time, denoted by Ψ2 , is provided by
Ψ2 =

T
X
k=1

f (s∗k + lk ).

(11)

Note that (10a) comprises nested expectations with respect
to the random PEV arrivals in the future time stages. Except
for few special cases, it is hard to provide the closed-form
of the optimal solution to (10). On the other hand, (10)
can be solved by the standard numerical methods, such as
backward reduction and the sample average approximation
(SAA) based on Monte Carlo sampling techniques [8]–[10],
[22]. These algorithms typically incur a computational complexity that grows exponentially with both the time span T
and the dimensions of state and decision spaces. Note that (10)
involves continuous state and decision spaces. Discretization
of these spaces leads to a curse of dimensionality, rendering
the computational complexity prohibitively high.
III. MPC- BASED O NLINE C HARGING A LGORITHM
In view of the extremely high complexity of standard
numerical methods, we are motivated to obtain a near-optimal
solution by solving a much simpler problem, which replace
all exogenous random variables by their expected values. This
is referred to as the expected value problem [8]–[10] or the
MPC approach [12]–[15] in the literature. Notice that the firstorder moment, i.e., expectation, of a random process is much
easier to estimate than the other statistics, e.g., variance or the
probability distribution. Thus, the assumption of knowing the
expected values is weaker than the assumptions in other EVcharging algorithms [20], which assume that the probability
distributions of the random processes are known.
Instead of solving the problem using generic convex optimization approaches, we propose a low-complexity online
Expected Load Flattening (ELF) algorithm by exploring the
load flattening feature of the optimal solution to the expected
value problem, as shown in Section III-A. Section III-B provides the theoretical analysis of the performance gap between
the optimal solution to the expected value problem and the
optimal solution to (10).

end of system time, denoted by Ψ3 , is defined as
Ψ3 =

min
sk

s. t.

f (sk + lk ) +

T
X

f (st + νt )

(12a)

t=k+1
j
X

st ≥

j
X

t=k

t=k

j

T
X

X
t=k

st ≤

t=k

d˜kt +

j
X

j
X

µm
n , ∀j

j

d˜kt +

X

em
X

0
d¯tt00 =

µm
n , ∀j = k, · · · , T. (12c)

In each time k, we solve problem (12) and obtain the optimal
charging solution s∗k . Then, problem (12) is resolved with the
updated d˜kt according to the realization of the PEVs arrived in
next time. So on and so forth, we obtain the optimal charging
solution s∗k for time stage k = 2, · · · , T . The total cost at the

(

0
d˜tt00 , for t00 = k, · · · , T, t0 = k,
0
µtt00 , for t00 = t0 , · · · , T, t0 = k + 1, · · · , T .

(14)

The key idea of online ELF algorithm is to balance the
charging load among all time slots k, · · · , T . Specifically, step
3 - 5 is to search the time interval [i∗ , j ∗ ] that has the maximum
load density among current time slots and record the maximum
load density. The time slots with maximum load density are
then deleted, and the process is repeated until the current time
slot k belongs to the maximum-density interval, i.e., i∗ = k.
Algorithm 1: Online ELF Algorithm

1
2
3

4

input : Dk , µt , t = k + 1, · · · , T
output: sk
initialization i = 0, j = 0;
repeat
For all time slot i = k, · · · , T, j = i, · · · , T , compute
Pj
P
0
( jt00 =t0 d¯tt00 + νt0 )
0
i∗ , j ∗ = arg max { t =i
}. (15)
k≤i≤j≤T
j−i+1
Set
y∗ =

5

m=k+1 n=m

(13)

where s∗k is the optimal solution to (12) at time stage k. The
solution to (12) is always feasible to (10) in the sense that
it always guarantees fulfilling the charging demand of the
current parking PEVs before their departures. This is because
the constraints of sk in (10) are included in (12).
Due to the convexity of f (·), the optimal solution is the one
that flattens the total load as much as possible. By exploiting
the load flattening feature of the solution, we present in
Algorithm 1 the online ELF algorithm that solves (12) with
complexity O(T 3 ). The online ELF algorithm have a lower
computational complexity than generic convex optimization
algorithms, such as the interior point method, which has a
complexity O(T 3.5 ) [21]. Notice that similar algorithms have
been proposed in the literature of speed scaling problems [23],
[24] and PEV charging problems [5]. The optimality and the
complexity of the algorithm have been proved therein, and
hence omitted here. The algorithm presented here, however,
paves the way for further complexity reduction to O(1) in
Section IV. For notation brevity, we denote in the online ELF
algorithm

= k, · · · , T, (12b)

m=k+1 n=m

f (s∗k + lk ),

k=1

A. Algorithm Description
Denote the expectation of ξ t as µt = [νt , µtt , · · · , µtT ],
where νt = E[ιt ], µtt0 = E[ηtt0 ], ∀t0 = t, · · · , T. Replacing ξ t
in (10) with µt , we obtain the following deterministic problem:

T
X

6
7

Pj ∗

Pj ∗

¯t0
t00 =t0 dt00
j ∗ − i∗ + 1

t0 =i∗ (

+ νt0 )

.

(16)

Delete time slot i∗ , · · · , j ∗ and relabel the existing time
slot t > j ∗ as t − j ∗ + i∗ − 1.
until i∗ = k;
Set sk = y ∗ − lk .

B. Optimality Analysis
In this subsection, we analyze the optimality of the solution to (12). Notice that MPC approximates the non-causal

random variables by their expected values regardless of their
distribution functions. As a result, such approximation may
lead to unacceptably large performance loss, depending on the
distribution of the random variables. Therefore, the MPC approximation is not always justifiable. A well-accepted metric,
Value of the Stochastic Solution (VSS) is adopted to evaluate
optimality gap between the optimal online solution and the
solution to the expected value problem [8]–[10]. Previous
work, e.g., [9], [10], mainly evaluates VSS using numerical
simulations. Whereas in our analysis, we show that VSS is
always bounded regardless of the distribution of the future
EV charging demands. This provides a strong theoretical
justification of adopting MPC approximation to solve our
problem.
Let Ξ denote a scenario, which is defined as a possible
realization of the sequence of random load demand [22],
Ξ = [ξ2 , ξ3 , · · · , ξT ].

(17)

Here, we treat ξ 1 as deterministic information since the
demand of PEVs arrived at the first stage is known by the
scheduler. Let Φ1 , Φ2 and Φ3 be the expectation of the optimal
value of the offline problem (4), the online problem (10)
and the expected value problem (12), respectively, where the
expectation is taken over the random scenarios. That is,
Φ1 = EΞ [Ψ1 (Ξ)] , Φ2 = EΞ [Ψ2 (Ξ)] , Φ3 = EΞ [Ψ3 (Ξ)] .

(18)

It has been proved previously [8], [9] that
Φ1 ≤ Φ2 ≤ Φ3 .

(19)

To assess the benefit of knowing and using the distributions
of the future outcomes, the VSS is defined as
VSS = Φ3 − Φ2 .

(20)

To show that the online ELF algorithm yields a bounded VSS,
we need to bound Φ3 and Φ2 . Generally, it is hard to calculate
Φ2 or analyze the lower bound of Φ2 directly [9], [10]. Thus,
we choose to analyze the lower bound of Φ1 instead, since
(19) shows that the lower bound of Φ1 is also the bound of
Φ2 . In what follows, we will show the lower bound of Φ1 in
Proposition 1 and the upper bound of Φ3 in Proposition 2.
Proposition 1:
Φ1 ≥ T f

!
Pe1 ˜1 PT Pet j PT
t=1 dt +
t=2
j=t µt +
t=1 νt
T

.

(21)

Proof: Please see the detailed proof in Appendix C.
Let O(t) be the set that O(t) = {(m, n)|em ≥ t, m =
1, · · · , t, n = t, · · · , em }. Then, we show that Φ3 is bounded
by Proposition 2.
Proposition 2: For any distribution of ξ t , t = 1, · · · , T ,
there is

Φ3 ≤ E 

T
X
t=1


f


X

ηnm + ιt  .

(22)

(m,n)∈O(t)

Proof: Please see the detailed proof in Appendix D.
Now, we are ready to present Theorem 3, which states that

the VSS is bounded for any distribution of random variables.
Theorem 3: For any distribution of random vector ξ t , t =
1, · · · , T, n = t, · · · , T, there is

VSS ≤ E 

T
X
t=1


f


X

ηnm + ιt  − T f

(m,n)∈O(t)



Γ
T


,

(23)

Pe1 ˜1 PT Pet j PT
dt + t=2 j=t µt + t=1 νt .
where Γ = t=1
Theorem 3 can be easily derived by Proposition 1 and
Proposition 2. In practice, the performance gap between the
online ELF algorithm and the optimal online algorithm is often
much smaller than the bound of VSS. This will be elaborated
in the numerical results in Section V.
IV. O NLINE ELF A LGORITHM UNDER F IRST- ORDER
P ERIODIC R ANDOM P ROCESSES
Notice that the complexity of O(T 3 ) of online ELF algorithm mainly comes from step 3, which exhaustively searches
the maximum-density period [i∗ , j ∗ ] over all subintervals
within [k, T ]. When the random arrival process is first-order
periodic stochastically 1 , we argue that the searching scope
can be reduced to one period from the whole system time
T . Thus, the complexity of step 3 is limited by the length
of a period instead of T . As a result, the complexity of
the algorithm reduces from O(T 3 ) to O(1), implying that
it does not increase with the system time T , and thus the
algorithm is perfectly scalable. In practice, the arrival process
of the charging demands are usually periodic stochastically.
For example, the arrival of charging demands at a particular location is statistically identical at the same time every
day during weekdays (or weekends). the National Household
Travel Survey (NHTS) 2009 gathers information about daily
travel patterns of different types of households in 2009, and
shows that the daily travel statistics (e.g., Average Vehicle
Trip Length, Average Time Spent Driving, Person Trips,
Person Miles of Travel) are very similar for each weekday or
weekend, but different between weekday and weekend [28]. In
Section IV-A, we investigate the case when the random load
demand process is first-order periodic. In Section IV-B, we
provide a closed-form solution to (12) for a special case when
the load demand process is the first-order stationary.
A. First-Order Periodic Process
In this subsection, we consider the case when the arrival
process is first-order periodic. Specifically, the first-order periodic process means that the first-order moment (i.e., mean)
of the random process is periodic. That is, at current time stage
k, for all t = k + 1, · · · , T, we have
µt = E[ξ t ] = µt+p ,

(24)

where ξ t is the random arrival events at time t, µt is the
expectation of ξ t , and p is the length of period. Then, instead
1 The first-order periodic stochastic process is defined as a stochastic process
whose first-order moment, i.e., mean is periodic. That is, the mean of the
random arrival events, i.e., µt is periodic. However, the actual realizations of
the arrival events ξt are uncertain and not periodic.

of considering µt for t = k + 1, · · · , T , we only need to
consider µt for one period, i.e., for t = k + 1, k + p:
k+1
k+1
µk+1 = [νk+1 , µk+1
k+1 , µk+2 , · · · , µk+e1 , 0, · · · , 0],
..
.

(25)

k+p
k+p
µk+p = [νk+p , µk+p
k+p , µk+p+1 , · · · , µk+ep , 0, · · · , 0].

Here, en ≤ T, n = 1, · · · , p is the maximum parking time
for PEVs arriving at time k + n. Specially, we define ê as
ê = max{ek+1 , ek+2 , · · · , ek+p }. We decompose the search
region {i, j|i = k, · · · , T, j = i, · · · , T } into three subregions, defined as Π1 = {i, j|i = k, j = k, · · · , k + ê}, Π2 =
{i, j|i = k, j = k + ê + 1, · · · , T } and Π3 = {i, j|i = k +
1, · · · , T, j = i, · · · , T }, respectively. We denote by X̂, Ŷ , Ẑ
the maximum densities of region Π1 , Π2 , Π3 , respectively.
Indeed, the largest of X̂, Ŷ , Ẑ is the maximum density of the
interval [i∗ , j ∗ ] ⊆ [k, T ] over all possible pairs i, j ∈ {i =
k, · · · , T, j = i, · · · , T }. Let [î1 , ĵ1 ], [î2 , ĵ2 ], [î3 , ĵ3 ] be the
intervals with the maximum density over region Π1 , Π2 , Π3 ,
respectively. By definition, î1 = î2 = k. Similar to the
stationary case, X̂ can be calculated by searching ĵ1 over
{k, · · · , k + ê}. That is,
˜k
n=k (dn

Pt
X̂ =

max

k≤t≤k+ê

P
P
+ νn ) + tn=k tm=n µn
m
.
n−k+1

(26)

Moreover, Lemma 1 shows that Ŷ and Ẑ can be calculated
once the maximum density of interval [k + 1, k + ê] has been
obtained. First, we introduce some definitions which help to
show Lemma 1. Let Π4 be a subset of Π3 , where Π4 is defines
as Π4 = {i, j|i = k + 1, · · · , k + ê, j = i, · · · , k + ê}, and
[ī, j̄] be the interval with
maximum
density of region Π4 , i.e.,
Pk+en n
Pj
m=n µm +νn )
n=i (
ī, j̄ = arg max
.
j−i+1
k+1≤i≤j≤k+ê

Lemma 1: The maximum densities of Π2 and Π3 are calculated by
Pĵ2
Ŷ =

n=k (

Pk+en
m=n

µn
m + νn )

ĵ2 − k + 1

Pĵ3
, Ẑ =

n=î3

(

Pk+en
m=n

µn
m + νn )

ĵ3 − î3 + 1

,

(27)

respectively, where î3 = ī,
ĵ2 =

and

(
max{j̄, k + ê + 1}, if j̄ < ī + p,
j̄ + (r − 1)p,
otherwise.

(
j̄,
if j̄ < ī + p,
ĵ3 =
j̄ + (r − 1)p, otherwise.

(28)

(29)

Proof: Please see the detailed proof in Appendix E.
Based on Lemma 1, we can modified the searching region
of step 3 of online ELF algorithm as follows:
• if j̄ < ī + p, the interval with the maximum density
during time stages [k + 1, T ] is [ī, j̄]. Then, in step 3
of the online ELF algorithm, the search region of i, j
is reduced from {i, j|i = k, · · · , T, j = i, · · · , T } to
{i, j|i = k, · · · , ī, j = i, · · · , ī, j̄}.
• If j̄ ≥ ī+p, the interval with the maximum density during
time stages [k + 1, T ] is [ī, j̄ + (r − 1)p]. Then, in step

3 of the online ELF algorithm, the search region of i, j
can be reduced from {i, j|i = k, · · · , T, j = i, · · · , T }
to {i, j|i = k, · · · , ī, j = i, · · · , ī, j̄ + (r − 1)p}.
As a result, the searching region of the online ELF algorithm
is only related to [k + 1, k + ê] instead of T . Thus, the
computational complexity of the online ELF algorithm is O(1)
instead of O(T 3 ) under first-order periodic process.

B. First-order Stationary Process
In this subsection, we show that the optimal solution to
(12) can be calculated in closed form if the arrival process is
first-order stationary. Due to the page limit, we only provide
the main results here. By first-order stationary, we mean that
the statistical mean of ξ t , i.e., νt and µtt0 , t0 = t, · · · , T only
depends on the relative time difference τ = t0 − t, but not the
absolute value of t. We can then replace νt by ν and replace
µtt0 by µτ , where τ = t0 − t. Then, µt is no longer a function
of t, and can be represented as
µ = [ν, µ1 , µ2 , · · · , µē , 0, · · · , 0],

(30)

where ē is the maximum parking time of a PEV. We denote
by X, Y, Z the maximum densities of region Π1 , Π2 , Π3 ,
respectively. Then, X is calculated by
X=
( Pn

t=k

max

P
d˜kt + n
j=1 (n − k − j + 1)µj + lk − ν
n−k+1

k≤n≤k+ē

)
+ν

,

(31)

and Y, Z are provided in Lemma 2.
Lemma 2: The maximum densities of Π2 and Π3 are
achieved by setting i2 = k, j2 = T, i3 = k + 1, j3 = T ,
and calculated by
Y =

Pk+ē ˜k Pk+ē
t=k dt +
j=1 (T − k − j + 1)µj + lk − ν
T −k+1
(T
−
k
−
j
+
1)µj
j=1
+ ν.
T −k

+ ν, (32a)

Pk+ē
Z=

(32b)

Proof: Please see the detailed proof in Appendix F. The
largest of X, Y, and Z is the maximum density of the
interval [i∗ , j ∗ ] ⊆ [k, T ] over all possible pairs i, j ∈ {i =
k, · · · , T, j = i, · · · , T }. Specially, if X or Y is the largest
one, then k is already contained in the maximum-density
interval, and thus X or Y is the optimal charging rate at time
k. On the other hand, if Z is the largest, then the maximumdensity interval, i.e., [k + 1, T ], does not include k. Following
Algorithm 1, we will delete the maximum-density interval and
repeat the process. Now, time slot k is the only remaining time
slot after deletion. This implies that all charging demands that
have arrived by time slot k should be fulfilled during time slot
k. These arguments are summarized in Proposition 3, which
provides the closed form solution to (12).
Proposition 3: When the random load demand process is
first-order stationary, the optimal charging schedule to (12) is

TABLE I
PARAMETER SETTINGS OF THE PEV TRAFFIC PATTERNS

(33a)
(33b)
(33c)

t=k

V. S IMULATIONS
In this section, we investigate the performance of the
proposed online ELF algorithm through numerical simulations.
All the computations are solved in MATLAB on a computer
with an Intel Core i3-2120 3.30GHz CPU and 8 GB of
memory. For comparison, we also plot the optimal solution to
(10), which is obtained by SAA method [22], and a heuristic
solution by the online AVG algorithm [6], which is obtained by
charging each PEV at a fixed rate, i.e., its charging demand
divided by its parking time. Let the expected cost of AVG
algorithm is denoted by Φ4 . Define the relative performance
loss of ELF and AVG compared with the optimal online
2
2
and Φ4Φ−Φ
, respectively. Similar to [5]
solution as Φ3Φ−Φ
2
2
[6], we adopt an increasing quadratic cost function in the
simulations, i.e., f (st + lt ) = (st + lt )2 . Note that the cost
function is increasing and strictly convex, since the load st +lt
is always non-negative.
A. Average Performance Evaluation
In this subsection, we evaluate the average performance of
the online ELF algorithm under three different traffic patterns,
i.e., light, moderate, and heavy traffics. In particular, the
system time is set to be 24 hours, and each time slot lasts
10 minutes. The PEV arrivals follow a Poisson distribution
and the parking time of each PEV follows an exponential
distribution [25]–[27]. The mean arrival and parking durations
of the three traffic patterns are listed in Table I. The main
difference lies in the arrival rates at the two peak hours, i.e.
12 : 00 to 14 : 00 and 18 : 00 to 20 : 00. The settings of the
peak hour match with the realistic vehicle trips in National
Household Travel Survey (NHTS) 2009 [28]. Specially, the
average number of total PEVs simulated in scenario 1, 2
and 3 are 104, 204 and 304, respectively. We choose the
base load profile of one day in the service area of South
California Edison from [7]. Each PEV’s charging demand is
uniformly chosen from [25, 35]kW h. Each point in Fig. 1,
Fig. 2 and Table II is an average of 105 independent instances
of the scenarios listed in Table I. In Fig. 2, the total loads
s∗k are plotted over time. We notice that the total load of
the online ELF algorithm follows closely to that of optimal
online algorithm, whereas that of the AVG algorithm has a
larger gap from the optimal online algorithm. The average
costs normalized by that of the optimal online algorithm
are plotted in Fig. 1. Moreover, the VSS and the relative
performance loss are listed in Table II. Both the figure and
the table show that ELF performs very close to the optimal
online algorithm. The VSS and the relative performance loss
are no more than 0.1536 and 0.38% respectively. In contrast,
the relative performance loss of AVG algorithm is up to 5.82%,

Arrival Rate (PEVs/hour)
S. 1 S. 2
S. 3

Time of Day
08 : 00 − 10 : 00
10 : 00 − 12 : 00
12 : 00 − 14 : 00
14 : 00 − 18 : 00
18 : 00 − 20 : 00
20 : 00 − 24 : 00
24 : 00 − 08 : 00

7
5
10
5
10
5
0

7
5
35
5
35
5
0

Mean Parking
Time (hour)

7
5
60
5
60
5
0

10
1/2
2
1/2
2
10
0

TABLE II
AVERAGE PERFORMANCE COMPARISON UNDER THREE TRAFFIC PATTERNS

Scenario

VSS

Φ3 −Φ2
Φ2

Φ4 −Φ2
Φ2

1
2
3

0.1178
0.1319
0.1536

0.19%
0.28%
0.38%

3.50%
4.46%
5.82%

which is more than 15 times of that of ELF algorithm. The
relative performance loss of the approximate online algorithm
reflects the percentage of the extra cost compared with the
optimal online algorithm. Obviously, the performance loss is
always the smaller the better. For example, from the report
of Rocky Mountain Institute, the average electricity rate is
11.2cents/kW h, and the average load for a charging station is
100kW [29]. Then, the expected electricity cost of a charging
station for one years is $967680. 6% relative performance loss
means AVG algorithm leads to $58060 extra cost, while the
proposed online algorithm with 0.38% relative performance
loss leads to $3677 extra cost, or an annual saving of $54383
compared with the AVG algorithm.
B. Complexity of The Online Algorithm ELF
In this subsection, we verify the computational complexity
of the online ELF algorithm and also compare the complexity
of online ELF algorithm with that of the optimal online
algorithm and online AVG algorithm, respectively. We still
adopt the SAA method as the optimal online algorithm. Since
the complexity of SAA is too high, truncation is often adopted
in the SAA to reduce the complexity at a cost of performance
loss [8]. As such, we also simulate the complexity of truncated
SAA with a truncation period of 3 hours. Each PEV’s charging
demand is uniformly chosen from [25, 35]kW h, the arrivals
follow a Poisson distribution and the parking time of each
1.07

Optimal online

Online ELF

Online AVG
1.0582

1.06
1.05

1.0446

1.04

Normalized cost

given by the following close-form:

X − lk ,
if X = max{X, Y, Z},




 Y − lk ,
if Y = max{X, Y, Z},
s∗k =
k+ē
X



d˜kt ,
otherwise.



1.035

1.03

1.02
1.01

1 1.0019

1

1.0028

1

1.0038

1
0.99
0.98
0.97

Scenario 1

Fig. 1.

Scenario 2

Scenario 3

Normalized costs of three algorithms in three scenarios.

120
Total load of optimal online
Total load of online ELF
Total load of online AVG

65
90

Total load of optimal online
Total load of online ELF
Total load of online AVG

110
100

60
80
Load (kW)

Load (kW)

Load (kW)

90
55

70

50

80
70

60
60
45

Total load of optimal online
Total load of online ELF
Total load of online AVG

50

50

40

40
40

8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 1 2 3 4 5 6 7 8
Time (hour)

8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 1 2 3 4 5 6 7 8
Time (hour)

(a) Scenario 1: light traffic

(b) Scenario 2: moderate traffic
Fig. 2.

(c) Scenario 3: heavy traffic

Base load and total load of three algorithms in three scenarios.

4

10

65
Optimal online without truncation
Optimal online with truncation
Online ELF
Online AVG

3

10

60

Load (kW)

Computional time (s)

8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 1 2 3 4 5 6 7 8
Time (hour)

2

10

55

50
1

10

45

0

10

Total load of optimal online
Total load of online ELF
Total load of online ORCHARD

40
1

Fig. 3.

2

3

4

5
6
Time (hour)

7

8

9

10

CPU computational time over the system time.

PEV follows an exponential distribution, where the arrival rate
and mean parking durations are the same as that in the peak
hours of scenario 2 in Section V.A. We simulate 10 cases in
total, where the system time are set to be 1, 2, · · · , 10 hours.
For each algorithm, we record the CPU computational time at
each time stage and calculate the average CPU computational
time as the sum of CPU computational times of all time
stages divided by the number of time stages. Each point in
Fig. 3 is an average of 100 independent instances. Fig. 3
shows that the CPU computational time of both the online
algorithm ELF and AVG almost grow linearly with the system
time. In contrast, for the SAA based online algorithm with
or without truncation, the average CPU computational times
grow very quickly as system time increases. We notice that
when system time increases to 4 hours, the optimal online
algorithm without truncation consumes more than 2 hours and
the optimal online algorithm with truncation consumes more
than 30 minutes. Meanwhile the proposed online algorithm
ELF only takes several seconds. It is foreseeable that the
computational complexity of optimal online algorithm will
become extremely expensive as we further increase the system
time.
C. Performance Comparison with Online Algorithm ORCHARD [5]
In this section, we compare the proposed online algorithm
ELF with the online algorithm ORCHARD proposed in refer-

8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 1 2 3 4 5 6 7 8
Time (hour)

Fig. 4.

Load comparison of ORCHARD and ELF

ence [5] on the properties of the optimality and computational
complexity. First, we evaluate the average performance of the
proposed online algorithm and online algorithm ORCHARD.
To facilitate the comparison, we also adopt the average performance of optimal online algorithm as a bench mark. For
the system parameters, we use the default settings of scenario
1 in Section V.A. We simulate 105 cases and plot the total
load (the sum of the base load and the PEV charging load)
over time for the three online algorithms in Fig. 4. In addition,
the average performance ratios normalized against the optimal
online solution are shown in Fig. 5 respectively. Fig. 4 shows
that compared with online algorithm ELF, the online algorithm
ORCHARD always produces a larger gap from the optimal
online solution. From Fig. 5, we can see that the proposed
online algorithm ELF achieves a much lower expected cost
than the online algorithm ORCHARD, which indicates that
online algorithm ELF owns a better average performance than
online algorithm ORCHARD.
To compare the computation complexity of online algorithm
ORCHARD and ELF, we adopt the similar case study in
Section V.D of reference [5]. Specifically, we simulate the
CPU computational time of online algorithms by varying the
arrival rates of the PEVs during one day. For the system
parameter, we use the same settings as scenario 1 in Section
V.A except the arrival rates, which are assumed to be the same
during 8 : 00 − 18 : 00 and 0 after 18 : 00. We vary the arrival

Optimal online

2

Online ORCHARD

Online ELF
1.73

1.8
1.6

Normalized cost

1.4
1.2

1.0019

1

1
0.8
0.6
0.4
0.2

problem. Instead of adopting the standard numerical methods with high complexity, we provide a MPC-based online
algorithm with O(T 3 )-complexity. We rigorously analyze the
performance gap between the solution of the MPC-based
approach and the optimal solution for any distribution of
exogenous random variables. Moreover, we show that the
proposed algorithm can be made scalable under the first-order
periodic process of load demand. Besides, our analyses are
validated through extensive simulations.

0

A PPENDIX
Fig. 5.
ELF

Average performance ratios of online algorithm ORCHARD and

A. Proof of Theorem 1:
5

Computional time (s)

10

Optimal online without truncation
Optimal online with truncation
Online ELF
Online ORCHARD

4

10

(s)

(e)

i∈{i|ti =1}

3

10

2

10
100

Fig. 6.

We use the inductive method to show that through EDF
scheduling, all the PEVs can be fulfilled charging before
deadlines. For n = 1, (3) becomes
X
X
di ≥ s1 ≥
di .
(34)

150

200

250
300
350
Number of PEVs

400

450

500

Thus, by EDF scheduling, we can first satisfy the demand of
PEVs whose deadline at time stage 1. That is, for any PEV
(e)
i ∈ {i|ti = 1}, we set
xi1 = di .

CPU computational time over the number of PEVs.

rate in 8 : 00−18 : 00 from 10 to 50 (PEVs/hour) that leads to
the average number of total PEVs during one day varies from
100 to 500. For each specified average number of PEVs, we
simulate the average performance of 108 independent instances
for the online ELF algorithm, the optimal online algorithm and
the online ORCHARD algorithm, respectively, and record the
average CPU computational times of each cases for the three
algorithms, respectively. The results are plotted in Fig. 6. Fig. 6
shows that the average CPU computational time of the optimal
online algorithm grows quickly as the number of total PEVs
increases, while the average CPU computational time of the
online ELF algorithm and the online ORCHARD algorithm
grows slowly as the number of total PEVs increases. When
the number of PEVs is 200, the average CPU computational
time of the optimal online algorithm without truncation is
more than 24 hours. Even for the the optimal online algorithm
with truncation, the average CPU computational time is about
100 minutes. Whereas, the proposed online algorithm ELF
only takes about 4 minutes. Fig. 6 also indicates that the
computational time of online ELF and online ORCHARD as
the number of the PEVs increases.
As a conclusion, the case study shows that the online
algorithm ELF has a better average performance than online
algorithm ORCHARD. It also indicates that the CPU computational time of online ELF and online ORCHARD are similar.
VI. C ONCLUSIONS
In this paper, we formulate the optimal PEV charging
scheduling problem as a finite-horizon dynamic programming

i∈{i|ti =1}

(35)

Assuming that for all time stage m, EDF scheduling can fulfill
charge all the PEVs which depart at or before time stage m,
i.e., there exists at least a set of xit ’s that satisfy
(e)

ti
X

(e)

xit = di , ∀i ∈ {i|ti

≤ m},

(36a)

(s)
t=ti

(s)

(e)

(e)

xit ≥ 0, ∀t = ti , · · · , ti , ∀i ∈ {i|ti
Since

m
X

st ≥

t=1

m
X
t=1

X

≤ m}.

di ,

(36b)

(37)

(e)
i∈{i|ti =t}

Pm
Pm P
then, t=1 st − t=1 i∈{i|t(e) =t} di represents the amount
i
of power which is outputted from the charging station during
time stage 1, · · · , m and charged to the PEVs with deadline
after time stage m. By EDF scheduling, once the PEVs which
depart at time m have been fulfilled charging, we will first
charge the PEVs which depart at time stage m + 1. Thus, if
m
X

st −

t=1

m
X

X

di ≥

t=1 i∈{i|t(e) =t}
i

X

di ,

(38)

(e)

i∈{i|ti =m+1}

we finish charging of PEVs with deadline m + 1, and then go
to charge the PEVs with deadline m + 2. If
m
X
t=1

st −

m
X

X

di <

t=1 i∈{i|t(e) =t}
i

X

di ,

(39)

(e)
i∈{i|ti =m+1}

then the
PEVs with
m + 1 have been charged as
Pm
Pmdeadline
P
power t=1 st − t=1 i∈{i|t(e) =t} di . At time stage m + 1.
i

Since

m+1
X

st ≥

t=1

m+1
X
t=1

X

di ,

(e)

i∈{i|ti

(40)

=t}

1) If x∗it1 = 0 for a particular PEV i at a time slot t1 ∈
(s)
(e)
{ti , · · · , ti }, then, by complementary slackness, we
have ωit1 > 0. From (47a),
f 0 (st1 + lt1 ) = λi + ωit1 .

then,

m
m
X
X
di − 
st −

X

sm+1 ≥

(e)

i∈{i|ti

t=1

=m+1}

t=1


X
(e)

i∈{i|ti


di  , (41)
=t}

which means all the PEVs with deadline m+1 can be fulfilled
charging. This is because we will charge the PEVs with
deadline m + 1 first by the EDF scheduling. Thus, there exists
at least a set of xit ’s that satisfy

(48)

2) If x∗it2 > 0 for PEV i during a time slot t2 ∈
(s)
(e)
{ti , · · · , ti }, we can infer from (47c) that ωit2 = 0.
Then,
f 0 (st2 + lt2 ) = λi .
(49)
On the other hand, since f (st +lt ) is a strictly convex function
of st + lt , then f 0 (st + lt ) is an increasing function. From the
above discussions, we get the following two conclusions:

(e)

ti
X

(e)

xit = di , ∀i ∈ {i|ti

= m + 1},

(42a)

1) If x∗it1 > 0, x∗it2 > 0, then by (49),

(s)

f 0 (st1 + lt1 ) = f 0 (st2 + lt2 ) = λi .

t=ti

xi,m+1 ≥ 0, ∀i ∈

(e)
{i|ti

= m + 1}.

(42b)

Combining (36) and (42), we get that all the PEVs whose
deadline at or before stage m + 1 can be fulfill charging, i.e.,
there exist at least a set of xit ’s that satisfy

Due to the monotonicity of f 0 (st ), we have s∗t1 + lt1 =
s∗t2 + lt2 .
2) If x∗it1 = 0, x∗it2 > 0, then by (48) and (49), there is
f 0 (st1 + lt1 ) = λi + ωit1 > f 0 (st2 + lt2 ) = λi . (51)

(e)
ti

X

(e)

xit = di , ∀i ∈ {i|ti

≤ m + 1},

(50)

Since f 0 (st ) is a increasing function, we have s∗t1 +lt1 ≥
s∗t2 + lt2 .

(43a)

(s)

t=ti

(s)

(e)

(e)

xit ≥ 0, ∀t = ti , · · · , ti , ∀i ∈ {i|ti

≤ m + 1}.

(43b)

Therefore, we can conclude that by EDF scheduling, there
always exists at least a set of xit ’s that is feasible to (2). This
completes the proof.


Consider two function fˆ(st + lt ) and f¯(st + lt ). Let x̂∗it and
denote the optimal solutions to (2) with f (st +lt ) replaced
by fˆ(st + lt ) and f¯(st + lt ), respectively. Define ŝ∗t , s̄∗t as
X
X
ŝ∗t =
x̂∗it , s̄∗t =
x̄∗it , t = 1, · · · , T,
(52)

x̄∗it

i∈I(t)

B. Proof of Theorem 2:
First, we show that if there exists a PEV parking in the
station at both time t1 and t2 , i.e.,
(s)

(e)

t1 , t2 ∈ {ti , · · · , ti },

(44)

respectively. Suppose that there exists a time slot t1 such that
ŝ∗t1 < s̄∗t1 .
Since

T
X

and
x∗it1 ≥ 0, x∗it2 > 0,
s∗t1 + lt1 ≥ s∗t2 + lt2 .

f 0(

(s)

(e)

xit + lt ) − λi − ωit = 0, i ∈ N , t = ti , · · · , ti ,

s̄∗t =

t=1

X

di ,

(54)

i∈N

ŝ∗t2 > s̄∗t2

(55)

ŝ∗t1 + ŝ∗t2 = s̄∗t1 + s̄∗t2

(56)

Thus, we can find a PEV i ∈ N such that

(47a)

x̂∗it1 < x̄∗t1 , x̂∗it2 > x̄∗it2 .

(57)

x̂∗it2 > 0

(58)

As a result,

(e)

λi (di −

T
X

and

i∈I(t)

ti
X

ŝ∗t =

(53)

there must exist another time slot t2 such that

(46)

The Karush-Kuhn-Tucker (KKT) conditions to the convex
problem (2) are
X

t=1

(45)

then the optimal total loads at time t1 and t2 must satisfy that

i∈I(t)

xit ) = 0, i ∈ N ,

(47b)

(s)
t=ti

ωit xit = 0, i ∈ N , t =

since
(s)
ti , · · ·

(e)
, ti ,

(47c)

where λ, ω are the non-negative optimal Lagrangian multipliers corresponding to (2b) and (2c), respectively. We separate
our analysis into the following two cases:

x̄∗it2

≥ 0. Based on (46), there is
ŝ∗t2 + lt2 ≤ ŝ∗t1 + lt1 .

(59)

Combining (53)(56)(59), we get
s̄∗t2 + lt2 < ŝ∗t2 + lt2 ≤ ŝ∗t1 + lt1 < s̄∗t1 + lt1 .

(60)

Since f¯(st + lt ) is a strictly convex function of st + lt , then,
based on (56) and (60), we have
f¯(s̄∗t1 + lt1 ) + f¯(s̄∗t2 + lt2 ) > f¯(ŝ∗t1 + lt1 ) + f¯(ŝ∗t2 + lt2 ). (61)

PT

t=2

Pet

j=t

µjt . Then, by Jensen’s inequality [30],

Ψ1 (E[Ξ]) ≥

T
X

f

PT

s +l

which indicates that P t=1T t t is a constant. Then, we see
T
PT
s +l
that t=1 (st + lt − t=1T t t )2 is a strictly convex function
of st + lt . This completes the proof.

C. Proof of Proposition 1:
First, we show that Ψ1 (Ξ) is a convex function of Ξ. For
any Ξ, define s∗t (Ξ) as the optimal solution that minimizes
Ψ1 (Ξ) subject to (4b) - (4c). For an arbitrary pair of Ξ0 and
Ξ00 , let Ξ000 = λΞ0 + (1 − λ)Ξ00 ∀λ ∈ [0, 1]. Then, the charging
schedule st (Ξ000 ) such that st (Ξ000 ) = λs∗t (Ξ0 )+(1−λ)s∗t (Ξ00 )
still satisfies (4b) - (4c) with Ξ = Ξ000 due to the linearity of
the constraints. Based on the convexity of f (st + lt ), we have
T
X
t=1

≤λ

T
X

f (s∗t (Ξ0 )

T
X

+ lt ) + (1 − λ)

t=1

(63)
f (s∗t (Ξ00 )

+ lt )

T
X

f (s∗t (Ξ000 ) + lt ) ≤

t=1

T
X

f (st (Ξ000 ) + lt ).

(64)

t=1

Combining (63) and (64), we have

T
X

T
X
n=t

d˜tn ≤

T
X

t=1

n=t−1

T
X

ηnt =

em
X X

ηnm ,

(68)

m∈M n=t

where the inequality holds since st−1 ≥ d˜t−1
t−1 and M =
{m|e
≥
t,
m
=
1,
·
·
·
,
t}
Thus,
we
have st ≤
m
P
m
,
where
O(t)
is
a
bounded
set
for t =
η
(m,n)∈O(t) n
hP
i
P
T
m
1, · · · , T . Therefore, E
f
(
+
ι
)
is an
η
t
t=1
(m,n)∈O(t) n
upper bound of Φ3 . This completes the proof.

E. Proof of Lemma 1:
We provide the proof by discussing the following two cases:
1)If j̄ ≥ ī + p, which means that [ī, j̄] and [ī + p, j̄ + p]
overlaps with each other, then the density of interval [ī, j̄ + p]
is higher than that of [ī, j̄], and the density of [ī, j̄ + 2p] is
higher than that of [ī, j̄ + p]. So on and so forth. Finally, we
see that the interval [ī, j̄ + (r − 1)p] has the maximum density
over region Π3 . Thus, we have î3 = ī, ĵ3 = j̄ + (r − 1)p, and
Pj̄+(r−1)p Pk+en n
( m=n µm + νn )
n=ī
.
j̄ + (r − 1)p − ī + 1

(69)

Likewise, for the region {i, j|i = k, j = k + ê + 1, · · · , T },
we have ĵ2 = j̄ + (r − 1)p, and
Ŷ =

Pj̄+(r−1)p Pk+en n
( m=n µm + νn )
n=k
.
j̄ + (r − 1)p − k + 1

(70)

2)If j̄ < ī + p, then, the density of interval [ī, j̄] is higher
than that of [ī, j̄ + p], and the density of [ī, j̄ + p] is higher
than that of [ī, j̄ + 2p]. So on and so forth. Finally, we see
that the interval [ī, j̄] has the maximum density over region
Π3 . Thus, we have î3 = ī, ĵ3 = j̄, and the corresponding
maximum density
n=ī

P n n
( k+e
m=n µm + νn )
.
j̄ − ī + 1

(71)

For the region Π2 , if j̄ ≤ k + ê + 1, then

f (s∗t (Ξ00 ) + lt ).

Thus, we have established the convexity of Ψ1 (Ξ) over the
set of Ξ. Therefore, we have
E [Ψ1 (Ξ)] ≥ Ψ1 (E[Ξ]).

et
X
n=t

(65)

t=1

On the other hand, based on
1, · · · , e1 , j = t, · · · , et , we

d˜t−1
− d˜t−1
n
t−1 +

Ẑ =

f (s∗t (Ξ0 ) + lt ) + (1 − λ)



D. Proof of Proposition 2:
PT
st ≤ n=t d˜tn holds for all stage t. On the other hand, we
have

Pj̄

f (s∗t (Ξ000 ) + lt )

t=1

≤λ

(67a)
!
Pe1 ˜1 PT Pet j PT
t=1 νt
j=t µt +
t=2
t=1 dt +
.
T
(67b)

This completes the proof.

t=1

for all λ ∈ [0, 1]. On the P
other hand, let s∗t (Ξ000 ) be the optimal
T
solution that minimizes t=1 f (st + lt ) subject to (4b) - (4c)
000
with Ξ = Ξ . Then,

T
X

= Tf

Ẑ =

f (st (Ξ000 ) + lt )

T

t=1

s̄∗t

This contradicts with the fact that the
is the optimal
total charging rate for objective function f¯(st + lt ). Therefore, the optimal charging solution s∗t is the same for any
strictly convex function f (st + lt ). Next, we show that optimal solution s∗t isPa load flattening solution that minimizes
T
PT
t=1 st +lt 2
) subjecting to (4b) and (4c).
t=1 (st + lt −
T
Based on the argument that s∗t is the same for any strictly
convex function fP(st + lt ), then it is equivalent to show that
T
PT
t=1 st +lt 2
) is a strictly convex function of
t=1 (st + lt −
T
st + lt . Since
P
PT
PT
i∈N di +
t=1 lt
t=1 st + lt
=
,
(62)
T
T

!
Pe1 ˜1 PT Pet j PT
j=t µt +
t=1 νt
t=2
t=1 dt +

Pk+ê+1 Pk+en n
n=k (
m=n µm + νn )
Ŷ =
.
ê + 2

If j̄ > k + ê + 1, then

(66)

the definition
of st ,P
d˜it , µjt , i
PT
e1 ˜1
have
t=1 st =
t=1 dt

Pj̄
Ŷ =

=
+

(72)

n=k (

Pk+en

This completes the proof.

µn
m + νn )
.
j̄ − k + 1
m=n

(73)



we have

F. Proof of Lemma 2:

Pj−i+2
µt
(j − i + 2 − t)µt + t=1
+ν
j−i+2
Pj−i+1
(82)
(j − i + 2 − t)µt
t=1
≥
+ν
j−i+1
= ρ(i, j),
Pj−i+1

ρ(i, j + 1) =

Let ρ(i, j) denote the maximum density of [i, j]. For any
i = k, j = k + ē + 1, · · · , T , the density of interval [i, j] is
given by
Pk+ē
Pj−k Pj−2k+1
µt + t=k d˜kt + lk + (j − k)ν
ρ(i, j) = k=1 t=1
j−k+1
Pk+ē
Pj−k
(j − k + 1 − t)µt + t=k d˜kt + lk − ν
+ ν.
= t=1
j−k+1
(74)
To prove that the maximum density is achieved by setting
j = T , we only need to show ρ(i, j) is a non-decreasing
function of j for each given i, i.e.,
ρ(i, j) ≤ ρ(i, j + 1), ∀k + ē + 1 ≤ j ≤ T − 1.
Since
Pj−k

t=1 (j

− k + 1 − t)µt +
j−k

Pk+ē ˜k
t=k dt

≤

j−k+1
X

µt ,

t=1

(75)

(76)

we have
ρ(i, j + 1)
Pj−k
P
˜k Pj−k+1 µt + lk − ν
(j − k + 1 − t)µt + k+ē
t=k dt +
t=1
+ν
= t=1
j−k+1
Pj−k
Pk+ē ˜k
(j − k + 1 − t)µt + t=k dt + lk − ν
+ν
≥ t=1
j−k
=ρ(i, j),
(77)

which implies (75). Hence, Y is the maximum density of
[k, j], j = k + ē + 1, · · · , T . Next, we show that Z is the
maximum density of [k + 1, T ]. For any k + 1 ≤ i ≤ j ≤ T ,
the density of interval [i, j] is given by
Pj−i+1 Pj−i+2−k
µt
t=1
ρ(i, j) = k=1
+ν
j−i+1
(78)
Pj−i+1
t=1 (j − i + 2 − t)µt
=
+ ν.
j−i+1
To prove that the maximum density is achieved by setting
i = k + 1, j = T , we only need to show ρ(i, j) is a nondecreasing function of j for each given i, i.e.,
(79)

and a non-increasing function of i for each given j, i.e.,
ρ(i, j) ≥ ρ(i + 1, j), ∀k + 1 ≤ i + 1 ≤ j ≤ T.

which implies (79). On the other hand, as
Pj−i
j−i+1
X
t=1 (j − i + 1 − t)µt
≤
µt , ∀k + 1 ≤ i ≤ j ≤ T,
j−i
t=1
(83)
then
Pj−i
(j − i + 1 − t)µt
ρ(i + 1, j) = t=1
+ν
j−i
Pj−i
Pj−i+1
µt
t=1 (j − i + 1 − t)µt +
t=1
+ν
≤
j−i+1
= ρ(i, j),
(84)
which implies (80). This completes the proof.

R EFERENCES

k + ē + 1 ≤ j ≤ T − 1,

ρ(i, j) ≤ ρ(i, j + 1), ∀k + 1 ≤ i ≤ j ≤ T − 1,

t=1

(80)

On one hand, since
Pj−i+1
j−i+2
X
t=1 (j − i + 2 − t)µt
≤
µt , ∀k + 1 ≤ i ≤ j ≤ T,
j−i+1
t=1
(81)

[1] W. Tang and Y. J. Zhang, “Online electric vehicle charging control
with multistage stochastic programming,” in IEEE 48th Annu. Conf.
Information Sciences and Systems (CISS), pp. 1-6, Mar. 2014.
[2] J. A. P. Lopes, F. J. Soares, and P. M. R. Almeida, “Integration of electric
vehicles in the electric power system,” Proc. of the IEEE, vol. 99, no. 1,
pp. 168-183, 2011.
[3] E. Sortomme, M. M. Hindi, S. D. J. MacPherson, and S. S. Venkata,
“Coordinated charging of plug-in hybrid electric vehicles to minimize
distribution system losses,” IEEE Trans. on Smart Grid, vol.2, no.1, pp.
198-205, 2011.
[4] Z. Ma, D. Callaway, and I. Hiskens, “Decentralized charging control of
large populations of plug-in electric vehicles,” IEEE Trans. on Control
Systems Technology, vol. 21, no. 1, pp. 67-78, 2013.
[5] W. Tang, S. Bi, and Y. J. Zhang, “Online coordinated charging decision
algorithm for electric vehicles without future information,” IEEE Trans.
on Smart Grid, vol. 5, no. 6, pp. 2810 - 2824, 2014.
[6] Y. He, B. Venkatesh, and L. Guan, “Optimal scheduling for charging and
discharging of electric vehicles,” IEEE Trans. on Smart Grid, vol. 3, no.
3, pp. 1095-1105, 2012.
[7] L. Gan, U. Topcu, and S. H. Low, “Optimal decentralized protocol for
electric vehicle charging,” IEEE Trans. on Power System, vol.28, iss. 2,
pp. 940-951, 2012.
[8] J. R. Birge and F. Louveaux, Introduction to Stochastic Programming,
New York: Springer, 1997.
[9] B. Defourny, D. Ernst, and L. Wehenkel, “Multistage stochastic programming: a scenario tree based approach to planning under uncertainty,” Decision Theory Models for Applications in Artificial Intelligence: Concepts
and Solutions, 2011.
[10] F. Maggioni and S. Wallace, “Analyzing the quality of the expected value
solution in stochastic programming,” Annals of Operations Research, pp.
37-54, 2012.
[11] R. Leou, C. Su, and C. Lu, “Stochastic analyses of electric vehicle
charging impacts on distribution network,” IEEE Trans. on Power System,
vol. 29, no. 3, pp. 1055-1063, May 2014.
[12] L. Rao and J. Yao, “SmartCar: smart charging and driving control
for electric vehicles in the smart grid,” IEEE Global Communications
Conference (GLOBECOM), pp. 2709-2714, Dec. 2014.
[13] N. Chen, L. Gan, S. H. Low, and A. Wierman, “Distributional analysis
for model predictive deferrable load control,” 53rd IEEE Conference on
Decision and Control (CDC), pp. 6433-6438, Dec. 2014.
[14] L. Gan, A. Wierman, U. Topcu, N. Chen, and S. H. Low, “Realtime deferrable load control: handling the uncertainties of renewable generation,”
in Proceedings of the fourth international conference on Future energy
systems (ACM e-Energy) , pp. 113-124, May 2013.

[15] S. Bansal, M. N. Zeilinger, and C. J. Tomlin, “Plug-and-play model
predictive control for electric vehicle charging and voltage control in
smart grids,” 53rd IEEE Conference on Decision and Control (CDC), pp.
5894-5900, 2014.
[16] G. Li and X. Zhang, “Modeling of Plug-in Hybrid Electric Vehicle
Charging Demand in Probabilistic Power Flow Calculations,” IEEE Trans.
on Vehicular Technology, vol. 63, no. 6, pp. 2600-2612, 2014.
[17] T. Zhang, W. Chen, Z. Han, and Z. Cao, “Charging scheduling of electric
vehicles with local renewable energy under uncertain electric vehicle
arrival and grid power price,” IEEE Trans. on Smart Grid, vol. 3, no.
1, pp. 492-499, 2012.
[18] I. Koutsopoulos, V. Hatzi, and L. Tassiulas, “Optimal energy storage
control policies for the smart power grid,” in Proc. of IEEE International
Conference on Smart Grid Communications (SmartGridComm), pp. 475480, 2011.
[19] L.
Huang, J.
Walrand, and K.
Ramchandran, “Optimal demand response with energy storage management,”
http://arxiv.org/pdf/1205.4297.pdf, 2012.
[20] W. Feller, An introduction to probability theory and its applications,
John Wiley & Sons, 2008.
[21] Y. Ye, Interior Point Algorithms: Theory and Analysis, WileyInterscience Press, 1997.
[22] A. Shapiro, D. Dentcheva, and A. Ruszczynski, Lectures on Stochastic
Programming: Modeling and Theory, MPS-SIAM, Philadelphia, 2009.
[23] F. Yao, A. Demers, and S. Shenker, “A scheduling model for reduced
cpu energy,” in Proc. IEEE Symp. Foundations of Computer Science, pp.
374-382, 1995.
[24] N. Bansal, T. Kimbrel, and K. Pruhs, “Speed scaling to manage energy
and temperature,” Journal of the ACM (JACM), vol. 54, no. 1, pp. 1-39,
2007.
[25] M. Alizadeh, A. Scaglione, J. Davies, and K. S. Kurani, “A scalable
stochastic model for the electricity demand of electric and plug-in hybrid
vehicles,” IEEE Trans. on Smart Grid, vol. 5, no. 2, pp. 848-860, 2014.
[26] X. Zhang and S. Grijalva, “An Advanced Data Driven Model for
Residential Electric Vehicle Charging Demand,” technique report, Georgia
Institute of Technology, 2015.
[27] S. Chen and L. Tong, “iEMS for large scale charging of electric vehicles
architecture and optimal online scheduling,” in Proc. IEEE Int. Conf.
Smart Grid Commun. (SmartGridComm), pp. 629-634, Nov. 2012.
[28] A. Santos, A. N. McGuckin, H. Y. Nakamoto, D. Gray, and S. Lis,
Summary of travel trends: 2009 national household travel survey, Federal
Highway Administration, Washington, DC, 2011.
[29] James W. May and Matt Mattila, “Plugging In: A Stakeholder Investment Guide for Public Electric-Vehicle Charging Infrastructure,” Rocky
Mountain Institute, 2009.
[30] S. Boyd and L. Vandenberghe, Convex optimization, Cambridge university press, 2004.

