# Optimal Vehicle Charging in Bilevel Power-Traffic Networks via Charging
Demand Function

Yufan Zhang, Sujit Dey Yuanyuan Shi The authors are with the Department of
Electrical and Computer Engineering, University of California San Diego, San
Diego, California 92093, US.

###### Abstract

Electric vehicle (EV) charging couples the operation of power and traffic
networks. Specifically, the power network determines the charging price at
various locations, while EVs on the traffic network optimize the charging
power given the price, acting as price-takers. We model such decision-making
processes by a bilevel program, with the power network at the upper-level and
the traffic network at the lower-level. However, since the two networks are
managed by separate entities and the charging expense term, calculated as the
product of charging price and charging demand, is nonlinear. Solving the
bilevel program is nontrivial. To overcome these challenges, we derive the
charging demand function using multiparametric programming theory. This
function establishes a piecewise linear relationship between the charging
price and the optimal charging power, enabling the power network operator to
manage EV charging power independently while accounting for the coupling
between the two networks. With the derived function, we are also able to
replace the nonlinear charging expense term with a piecewise quadratic one,
thus guaranteeing solution optimality. Our numerical studies demonstrate that
different traffic demands can have an impact on charging patterns and the
power network can effectively incentivize charging at low-price nodes through
price setting.

Keywords: electric vehicles, bilevel program, coupled networks, charging
demand function

##  I Introduction

There is a growing trend towards electrifying a significant portion of the
transportation sector to decarbonize [1]. As of 2022, electric vehicle (EV)
sales accounted for 10% of all new car sales globally [2]. The increasing
number of EVs has raised concerns about their impact on the power network [3].
Without proper coordination, EV charging can compromise grid reliability. As
reported in [4], as little as 11% of heavy-duty vehicles charging
simultaneously in Texas can lead to significant voltage violations on the
transmission network. Therefore, careful management of EV charging is crucial.

Modeling a charging station as an aggregator is one of the approaches for
organizing charging smartly [5, 6]. This approach treats the charging demand
at a station as a flexible load that can be coordinated to meet the needs of
the grid. For example, it allows for EV charging to be delayed until there is
enough generation capacity on the grid, preventing high charging demand from
compromising grid reliability [7]. Previous studies have investigated the use
of aggregated EV charging for various purposes, such as frequency regulation
[8], ramping provision [9], and congestion relief [10], to name a few. To
achieve these goals, these studies treat EV driving behavior information such
as arrival time, charging duration, and energy demand as input parameters [6],
which can be obtained by forecasting approaches. From the perspective of
aggregators, the traveling plan of EVs and the charging management are viewed
as two separate decision-making processes. Although this assumption can
simplify the analysis, the coupling of those decision-makings is ignored,
which may lead to some unreliable decisions [11]. For example, if the EV
charging manager fails to consider the fact that EVs may travel to other
charging stations in response to electricity prices, it could result in a
charging demand hotspot where many EVs travel to a single station with low
electricity prices, overloading the charging station.

Linking vehicle travel plans with charging management has been advocated in
recent years [12]. For optimizing the collective driving behavior, an entity
represented by an independent traffic system operator (ITSO), is introduced
for settling the travel plans of vehicle drivers. Specifically, ITSO solves
the optimal traffic assignment problem, determining routes, charging
locations, and charging power for all vehicles in a way that minimizes total
traveling and charging costs. ITSO acts as a price-taker to the locational
marginal price (LMP) of electricity, which is determined by the optimal power
flow (OPF) problem solved by the independent distribution system operator
(IDSO). Research has investigated the joint optimization of ITSO’s traffic
assignment problem and IDSO’s OPF problem, demonstrating that the social
welfare of the coupled networks can be achieved under the LMP mechanism [11,
13]. Also, other studies have taken a game-theoretic perspective, with both
ITSO and IDSO seeking to optimize their own costs [14]. These studies have
shown that the power-traffic system game is a potential game, with the Nash
equilibrium coinciding with the joint optimization solution.

One challenge of coordinating the coupled power and traffic networks using the
joint optimization framework is that the two network operators belong to
separate entities, making it difficult to disclose the entire decision model
to each other. Recent works have resorted to distributed optimization for
solving the joint optimization problem, requiring multiple rounds of
information exchange between ITSO and IDSO [11, 13, 14]. However, such
iterative updates are not currently within the power network operation
structure, and iterative algorithms may not meet the time requirements of the
power network operation, as the lack of two-way real-time communication
between the systems makes multiple rounds of information exchange time-
consuming [15].

To maintain the existing operation framework of the power network, it is
desirable to solve the OPF problem independently and efficiently while
accounting for the coupling between power and traffic networks. One possible
approach is to model the interaction between the two networks as a bilevel
program [16], with the power network at the upper-level determining the LMP
and the traffic network at the lower-level determining the optimal charging
demand given the price. Interestingly, this bilevel formulation is equivalent
to the joint optimization and game-theoretic formulations proposed in previous
studies [11, 13, 14]. We will explore this equivalence further in Section
III-D.

However, solving the bilevel power-traffic optimization problem is not
trivial. Firstly, the charging expense that needs to be minimized is a product
of the charging price (an upper-level decision variable) and the charging
demand (a lower-level decision variable). The resulting bilinear term
introduces nonlinearity into the problem, making it difficult to solve with an
optimality guarantee. Secondly, solving the bilevel problem requires a
collection of the Karush-Kuhn-Tucker (KKT) conditions of the lower-level
problem, which in turn requires detailed knowledge of the ITSO decision model.

To address these challenges, we propose a novel approach for solving the
bilevel power-traffic optimization problem. Our approach is based on the
observation that the decision problems of the power and traffic networks are
coupled by the charging price and demand. If the relationship between the
charging price and the optimal charging demand can be obtained, the IDSO can
use such a function for solving the OPF problem independently, without the
need of knowing the detailed decision model of ITSO. Similar to the concept of
demand function [17] (which expresses the demand as a function of price), we
define such a relationship as the _charging demand function_.

As a critical step in our approach, we derive the charging demand function
using multiparametric programming theory [18, 19]. This function describes the
relationship between the charging price and the optimal charging power, which
is determined by the ITSO. By introducing the charging demand function, we are
able to “kill two birds with one stone”. On the one hand, the power network
can set the charging price and manage EV charging power by solving the OPF
problem independently, which aligns with the current operation framework of
the power network. Furthermore, the charging demand function is proven to be
piecewise linear, allowing us to transform the bilinear term in the charging
expense into a piecewise quadratic term. This transformation makes the bilevel
program for solving EV charging management with _optimality guarantees_.

Compared with existing studies, the main contributions of the paper are:

1) A theoretical derivation of the charging demand function, which describes
the relationship between the charging price and the optimal charging demand as
a piecewise linear function. This function allows for efficient vehicle travel
and charging management by the power network itself.

2) A solution strategy to the bilevel program that replaces the lower-level
ITSO problem with the charging demand function. This transformation converts
the bilinear term of the charging expense into a piecewise quadratic term,
enabling a computationally efficient solution of the bilevel power-traffic
program with optimality guarantees.

3) A complimentary understanding of different perspectives of modeling the
power and traffic networks interaction. We show our bilevel program is
equivalent to the models from the perspectives of joint optimization and game
theory, but using different solution techniques that do not require detailed
traffic models nor iterative information exchange.

The remaining parts of this paper are organized as follows. Section II
presents the bilevel setup and the decision-making problems of power and
traffic networks as preliminaries, whereas the details of charging demand
function derivation and solution strategy are given in Section III. Results
are discussed and evaluated in Section IV, followed by the conclusion and
future works.

Notation: The notation $\bm{x}_{\mathcal{I}}=[x_{i}]_{i\in\mathcal{I}}$
denotes a column vector comprised by elements in the set $\mathcal{I}$.
$|\cdot|$ denotes the cardinality of a set. $\bm{I}_{M}$ denotes an identity
matrix with the size of $\mathbb{R}^{M\times M}$. $\bm{1}_{M}$ is an all-one
column vector with the dimension of $\mathbb{R}^{M}$. $\bm{O}$ is an all-zero
matrix with a size to be defined by the use case. Given the column vectors
$\bm{x}_{1},\bm{x}_{2}$, we use the notation $[\bm{x}_{1};\bm{x}_{2}]$ to
denote the vertical stack of them.

##  II Preliminaries: Bilevel Setup and Networks Model

In this section, we first introduce the interaction between the power and
traffic networks as a bilevel framework in Section II-A. The decision problems
of traffic and distribution networks are described in Section II-B and Section
II-C, respectively.

###  II-A Bilevel Setup

Figure 1: Diagram of the bilevel setup. The problem of the ITSO is at the
lower-level, while the problem of IDSO is at the upper-level.

The setup of the problem is sketched in Fig. 1. We study the system-level
pricing strategy of non-profit IDSO, who manages the distribution network and
solves the OPF problem to minimize the generation cost for satisfying the
fixed load and EV charging demand. The resulting nodal electricity price,
namely LMP, affects the EVs’ routing and charging decisions. The traffic
network is managed by the non-profit ITSO, who optimally assigns the traffic
flow to minimize the travel cost brought by congestion and charging expenses.
The power and traffic networks are coupled by common infrastructure, namely
the EV charging stations. We follow the common assumption that the ITSO acts
as a price-taker in the decision problem, which optimizes EVs’ charging
behaviors given the charging price determined by IDSO. To model such a
relationship, a bilevel framework is considered here, where the decision
problem of IDSO solving the OPF is at the upper-level, and the optimal traffic
assignment problem of ITSO is at the lower-level. The two problems are coupled
by the decisions regarding the electricity price $\lambda_{i}$ and the
corresponding optimal charging demand $d_{i}^{*}$ at the node $i$ connected to
the charging station. The relationship between $\lambda_{i}$ and $d_{i}^{*}$
will be discussed in detail in Section III. For simplicity, we consider a
single-period problem but the formulation can be extended to the multi-period
case as future work.

###  II-B Traffic Network Model

The traffic network is modeled by a connected directional graph
$\mathcal{R}=(\mathcal{V},\mathcal{A})$, where $\mathcal{V}$ and $\mathcal{A}$
denote the sets of vertexes and arcs. We model the traveling time on arc $a$
as a linear function $\tau_{a}(\xi_{a})$ such that the traveling time
increases with the increase of the traffic flow $\xi_{a}$ on arc
$a\in\mathcal{A}$. ; see [11].

$\displaystyle\tau_{a}(\xi_{a})=\xi_{a}^{0}+\xi_{a}/R_{a},\forall
a\in\mathcal{A}$ (1)

where $\xi_{a}^{0},R_{a}$ are the parameters related to the shape of the
traveling time function.

Let $\mathcal{C}\in\mathcal{V}$ denote the set of nodes where the charging
stations are located in. We model the nodes of charging stations as “virtual
arcs” following [11] (see Fig. 2). At a node $i\in\mathcal{C}$, the EV driver
can choose the arc $l_{i}^{c}$ for charging or the arc $l_{i}^{n}$ for
skipping charging. The former decision incurs costs of time and electricity,
while the latter one incurs no cost. Concretely, the time costs for an EV on a
charging arc $l_{i}^{c}$ consist of two parts: waiting time and charging time.
The waiting time increases as the flow of EVs $\xi_{a},\forall a\in l_{i}^{c}$
increases and is modeled using the function in (1) [20] as well. The charging
time is associated with the average charging demand $e_{a}$, which is given by
$e_{a}/\rho_{a}$, where $\rho_{a}$ is the charging rate. The sum of waiting
and charging time constitutes the latency time spent in a charging station,
i.e.,

$\tau_{a}(\xi_{a})=e_{a}/\rho_{a}+\xi_{a}^{0}+\xi_{a}/R_{a},\forall
a\in\cup_{i\in\mathcal{C}}l_{i}^{c}$ (2)

For an EV on the virtual arc $l_{i}^{n}$, the time cost is zero, i.e.,

$\tau_{a}(\xi_{a})=0,\forall a\in\cup_{i\in\mathcal{C}}l_{i}^{n}$ (3)

With the notion of the virtual arcs, the extended traffic network is
introduced, with the set of arcs given by

$\mathcal{A}^{e}=\cup_{i\in\mathcal{C}}(l_{i}^{c}\cup
l_{i}^{n})\cup\mathcal{A}$ (4)

Figure 2: The illustration of virtual arcs to represent the charging behavior
in a charging station. An extra node (node #2) is introduced to form the
virtual arcs.

The ITSO optimizes the aggregated cost of total traffic system latency time
and charging expense which is given by

$\sum_{a\in\mathcal{A}^{e}}\gamma\cdot\xi_{a}\cdot\tau_{a}(\xi_{a})+\sum_{a\in\cup_{i\in\mathcal{C}}l_{i}^{c}}\lambda_{i}\cdot
e_{a}\cdot\xi_{a}\,,$ (5)

where $\gamma$ is the cost coefficient for one unit of travel time,
$\lambda_{i}$ is the electricity price at the charging station
$i\in\mathcal{C}$, which is obtained as the optimal dual solution of the OPF
problem solved by IDSO. From now on, we use the notation $d_{i}$ to denote the
charging demand in charging station $i$,

$d_{i}=e_{a}\cdot\xi_{a},\forall a\in\cup_{i\in\mathcal{C}}l_{i}^{c}$ (6)

Let the set $\mathcal{W}$ denote the set of finite origin–destination (O–D)
pairs. $\forall w\in\mathcal{W}$, there are a set of feasible routes
$\mathcal{R}_{w}$ that allows the driver to finish the trip on the extended
graph. Let $m_{w}$ be the travel demand of O-D pair $w$ and $f_{w}^{k}$ be the
flow on route $k\in\mathcal{R}_{w}$. We use
$\bm{f}_{w}\in\mathbb{R}^{|\mathcal{R}_{w}|}=[f_{w}^{k}]_{k\in\mathcal{R}_{w}}$
to denote the vector formed by the route flows of the O-D pair $w$, and define
$\bm{f}=[\bm{f}_{w}]_{w\in\mathcal{W}}$ as a collection of all the
$\bm{f}_{w}$. Let
$\bm{\xi}=[[\xi_{a}]_{a\in\mathcal{A}};[\xi_{a}]_{a\in\cup_{i\in\mathcal{C}}l_{i}^{n}};[\xi_{a}]_{a\in\cup_{i\in\mathcal{C}}l_{i}^{c}}]$
denote the arc flows on $\mathcal{A}^{e}$. The dimension of $\bm{\xi}$ is
$|\mathcal{A}^{e}|$, which is equal to $|\mathcal{A}|+2\cdot|\mathcal{C}|$. To
link the arc with the route, the link matrix
$\bm{A}_{w}\in\mathbb{R}^{|\mathcal{A}^{e}|\times|\mathcal{R}_{w}|}$ is
introduced for each O-D pair $w\in\mathcal{W}$, where the element in the
$a_{th}$ row is assigned with 1 if the arc $a$ is on the route belonging to
$\mathcal{R}_{w}$ and 0 otherwise. With these notations, the decision problem
of ITSO is formulated as

$\displaystyle\mathop{\min}_{\bm{\xi},\bm{f}}\sum_{a\in\mathcal{A}^{e}}\gamma\cdot\xi_{a}\cdot\tau_{a}(\xi_{a})+\sum_{i\in\mathcal{C}}\lambda_{i}\cdot
d_{i}$ (7a) $\displaystyle s.t.\ \ \eqref{6}$ (7b)
$\displaystyle\qquad\bm{1}_{|\mathcal{R}_{w}|}^{\top}\cdot\bm{f}_{w}=m_{w},\forall
w\in\mathcal{W}$ (7c)
$\displaystyle\qquad\bm{\xi}=\sum_{w\in\mathcal{W}}\bm{A}_{w}\cdot\bm{f}_{w}$
(7d) $\displaystyle\qquad 0\leq\bm{\xi}\leq\bm{\overline{\xi}}$ (7e)

where (7c) ensures that the sum of traffic flow on the routes of the O-D pair
$w$ equals the travel demand $m_{w}$. (7d) maps the flow on routes to the flow
on arcs. (7e) ensures the traffic flows on the arcs are non-negative and
within the limits. With the quadratic term introduced by the first term in
(7a), the decision problem of ITSO is a quadratic program (QP).

Remark 1: We note that the decision problem of ITSO in (7) can be extended to
the case considering EV discharging. For that, the battery degradation cost
$c(d_{i})$ will be considered in the objective, and modeled by a quadratic
function [21]. Also, the charging/discharging demand is modeled in a range,
i.e., $-e_{a}^{d}\cdot\xi_{a}\leq d_{i}\leq e_{a}^{c}\cdot\xi_{a},\forall
a\in\cup_{i\in\mathcal{C}}l_{i}^{c}$, where average discharging/charging rates
$e_{a}^{d},e_{a}^{c}$ are constants. When all the EVs in the traffic flow of
arc $a$ charge, the charging demand $d_{i}$ equals $e_{a}^{c}\cdot\xi_{a}$.
When all the EVs discharge, the charging demand $d_{i}$ equals
$-e_{a}^{d}\cdot\xi_{a}$. When part of the EVs charge and the remaining ones
discharge, the charging demand is in the middle of the range. To sum up,
considering discharging, the traffic network model is still a QP.

###  II-C Distribution Network Model

We use a graph $\mathcal{G}(\mathcal{N},\mathcal{L})$ to represent the power
distribution network, where $\mathcal{N}$, $\mathcal{L}$ are the sets of nodes
and lines, respectively. As charging stations are the common infrastructures
of distribution and traffic systems, we use $\mathcal{C}$ to denote the
charging stations. To satisfy the charging demand and fixed load, IDSO
optimally dispatches the generators under the network constraint and solves
the OPF problem. Let $\bm{g}=[g_{i}]_{i\in\mathcal{N}}$,
$\bm{v}=[v_{i}]_{i\in\mathcal{N}}$,
$\bm{\theta}=[\theta_{i}]_{i\in\mathcal{N}}$ and
$\bm{l}=[l_{i}]_{i\in\mathcal{N}}$ denote the vectors of generation, the
voltage amplitude and phase angle, and fixed load, respectively, and
$\bm{d}=[d_{i}]_{i\in\mathcal{N}}$ denote the vector of charging demand, where
$d_{i}=0,\forall i\in\mathcal{N}\backslash\mathcal{C}$. Given the generation
cost coefficient $\bm{c}=[c_{i}]_{i\in\mathcal{N}}$, the problem of IDSO
minimizes the total generation cost subject to the demand requirement and
power flow constraints,

$\displaystyle\mathop{\min}_{\bm{g},\bm{v},\bm{\theta}}\bm{c}^{\top}\bm{g}$
(8a) $\displaystyle s.t.\ \
0\leq\bm{g}\leq\overline{\bm{g}}:\underline{\bm{\tau}},\overline{\bm{\tau}}$
(8b) $\displaystyle\quad\
l_{i}+d_{i}+\sum_{j\in\Omega_{i}}K_{1ij}(v_{i}-v_{j})+$
$\displaystyle\quad\quad\quad\sum_{j\in\Omega_{i}}K_{2ij}(\theta_{i}-\theta_{j})=g_{i}:\lambda_{i},\forall
i\in\mathcal{N}$ (8c) $\displaystyle\quad\ -f_{ij}\leq
K_{1ij}(v_{i}-v_{j})+K_{2ij}(\theta_{i}-\theta_{j})\leq f_{ij}:$
$\displaystyle\qquad\qquad\qquad\qquad\quad\underline{\eta_{ij}},\overline{\eta_{ij}},\forall
i\in\mathcal{N},\forall j\in\Omega_{i}$ (8d) $\displaystyle\quad\
\underline{v}\leq
v_{i}\leq\overline{v}:\underline{\mu_{i}},\overline{\mu_{i}},\forall
i\in\mathcal{N},$ (8e)

where the dual variables are given after colons, and
$\underline{\bm{\tau}}=[\underline{\tau_{i}}]_{i\in\mathcal{N}},\overline{\bm{\tau}}=[\overline{\tau_{i}}]_{i\in\mathcal{N}}$.
Under the generation capacity $\overline{\bm{g}}$, (8b) limits the range of
generation. Let $\Omega_{i}$ denote the set of buses connected to node $i$.
(8c) is the nodal power balance constraint and linearized distribution power
flow model [22] is adopted, where
$K_{1ij}=\frac{x_{ij}r_{ij}}{r_{ij}^{2}+x_{ij}^{2}},K_{2ij}=\frac{x_{ij}^{2}}{r_{ij}^{2}+x_{ij}^{2}}$
are the parameters calculated by the line resistance $r_{ij}$ and line
reactance $x_{ij}$. (8d) and (8e) bound the range of power flow and voltage
amplitude, with $f_{ij}$ and $\underline{v},\overline{v}$ as the limits.

For the IDSO decision problem in (8), the charging demand $d_{i}$ serves as
the parameters and the electricity price $\lambda_{i}$ is the dual variable of
the constraint (8c). To clearly show the relationship between the problems of
IDSO and ITSO, the dual problem of (8) is derived, whose derivation is
achieved by restricting the derivative of Lagrange function of (8) to
$\theta_{i},v_{i},g_{i}$ to zero. We define the vectors
$\bm{\overline{\mu}}=[\overline{\mu_{i}}]_{i\in\mathcal{N}},\bm{\underline{\mu}}=[\underline{\mu_{i}}]_{i\in\mathcal{N}},\bm{\lambda}=[\lambda_{i}]_{i\in\mathcal{N}},\bm{\underline{\eta}}=[\underline{\eta_{ij}}]_{i\in\mathcal{N},j\in\Omega_{i}},\bm{\overline{\eta}}=[\overline{\eta_{ij}}]_{i\in\mathcal{N},j\in\Omega_{i}}$
as the collections of dual variables. The dual problem of (8) is given by,

$\displaystyle\mathop{\max}_{\bm{\overline{\tau}},\bm{\underline{\tau}},\bm{\overline{\mu}},\bm{\underline{\mu}},\bm{\lambda},\bm{\underline{\eta}},\bm{\overline{\eta}}}\sum_{i\in\mathcal{N}}-\overline{\mu_{i}}\cdot\overline{v_{i}}+\underline{\mu_{i}}\cdot\underline{v_{i}}-\overline{\tau_{i}}\cdot\overline{g_{i}}-$
$\displaystyle\sum_{i\in\mathcal{N}}\sum_{j\in\Omega_{i}}(f_{ij}\overline{\eta_{ij}}+f_{ij}\underline{\eta_{ij}})+\sum_{i\in\mathcal{N}}\lambda_{i}\cdot
l_{i}+\sum_{i\in\mathcal{C}}\lambda_{i}\cdot d_{i}$ (9a) $\displaystyle
s.t.\sum_{j\in\Omega_{i}}K_{2ij}(\lambda_{i}-\lambda_{j}+\overline{\eta_{ij}}-\overline{\eta_{ji}}-\underline{\eta_{ij}}+\underline{\eta_{ji}})=0,\forall
i\in\mathcal{N}$ (9b) $\displaystyle\quad\
\sum_{j\in\Omega_{i}}K_{1ij}(\lambda_{i}-\lambda_{j}+\overline{\eta_{ij}}-\overline{\eta_{ji}}-\underline{\eta_{ij}}+\underline{\eta_{ji}})+$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\overline{\mu_{i}}-\underline{\mu_{i}}=0,\forall
i\in\mathcal{N}$ (9c) $\displaystyle\qquad
c_{i}-\lambda_{i}-\underline{\tau_{i}}+\overline{\tau_{i}}=0,\forall
i\in\mathcal{N}$ (9d) $\displaystyle\qquad\overline{\mu_{i}}\geq
0,\underline{\mu_{i}}\geq 0,\overline{\tau_{i}}\geq 0,\underline{\tau_{i}}\geq
0,\overline{\eta_{ij}}\geq 0,\underline{\eta_{ij}}\geq 0.$ (9e)

Given $\underline{\tau_{i}}\geq 0$, (9d) can be further simplified as
$c_{i}-\lambda_{i}+\overline{\tau_{i}}\geq 0$. With the dual problem of IDSO
in (9), we model the interaction between IDSO and ITSO by a bilevel program in
the next Section.

##  III Optimal Electricity Pricing Strategy

In Section III-A, we formulate the interaction between IDSO and ITSO as a
bilevel program and point out the key difficulties of solving it. As one
critical step of the proposed solution strategy, we derive the charging demand
function in Section III-B. With it, we propose an effective solution strategy
for the bilevel program in Section III-C.

###  III-A Formulation of the Bilevel Power-Traffic Optimization

The studied interaction problem between distribution and traffic networks has
a bilevel structure, where the problem of ITSO in (7) is nested in the problem
of IDSO in (9). The bilevel program is formulated as follows

$\begin{split}&\mathop{\min}_{\bm{\overline{\tau}},\bm{\underline{\tau}},\bm{\overline{\mu}},\bm{\underline{\mu}},\bm{\lambda},\bm{\underline{\eta}},\bm{\overline{\eta}}}\sum_{i\in\mathcal{N}}\overline{\mu_{i}}\cdot\overline{v_{i}}-\underline{\mu_{i}}\cdot\underline{v_{i}}+\overline{\tau_{i}}\cdot\overline{g_{i}}+\\\
&\qquad\ \
\sum_{i\in\mathcal{N}}\sum_{j\in\Omega_{i}}(f_{ij}\overline{\eta_{ij}}+f_{ij}\underline{\eta_{ij}})-\sum_{i\in\mathcal{N}}\lambda_{i}\cdot
l_{i}-\sum_{i\in\mathcal{C}}\lambda_{i}\cdot d_{i}^{*}\\\
&s.t.\eqref{9(b)},\eqref{9(c)},\eqref{9(d)},\eqref{9(e)}\\\
&\quad\\{d_{i}^{*}\\}_{i\in\mathcal{C}}\in\mathop{\arg\min}_{\bm{\xi},\bm{f}}\sum_{a\in\mathcal{A}^{e}}\gamma\cdot\xi_{a}\cdot\tau_{a}(\xi_{a})+\sum_{i\in\mathcal{C}}\lambda_{i}\cdot
d_{i}\\\ &\qquad\qquad\quad\
s.t.\eqref{6},\eqref{7(c)},\eqref{7(d)},\eqref{7(e)}\end{split}$ (10)

where the maximization problem of IDSO in (9) is replaced with an equivalent
minimization problem by minimizing over the negative objective in (9a). Once
the bilevel program in (10) is solved, the electricity price for EV charging
is obtained from the optimal solution of the dual variable
$\lambda_{i}^{*},\forall i\in\mathcal{C}$.

Solving the bilevel optimization problem in (10) is typically done by
replacing the lower-level problem with its KKT conditions for a single-level
reduction [16]. However, this is not a trivial task. The presence of a
bilinear term $\sum_{i\in\mathcal{C}}\lambda_{i}\cdot d_{i}^{*}$ in the
objective function introduces nonlinearity, making the problem computationally
challenging to solve with the optimality guarantee. Additionally, the IDSO and
ITSO are separate entities, and the IDSO may not possess a detailed model of
the ITSO [11], which can make it difficult to derive the KKT conditions
necessary for solving (10) in practice.

Fortunately, the dual variable $\lambda_{i},\forall i\in\mathcal{C}$ in the
upper-level problem serves as the parameter in the lower-level problem of
ITSO, where the charging demand $d_{i},\forall i\in\mathcal{C}$ is the primal
variable. The relationship between the price $\lambda_{i},\forall
i\in\mathcal{C}$ and the optimal charging demand $d_{i}^{*},\forall
i\in\mathcal{C}$ is bijective and piecewise linear [23]. We will explicitly
derive this function in Section III-B. For now, let’s denote such bijection as
the charging demand function,

$\bm{d}^{*}_{c}=\pi(\bm{\lambda}_{c}),$ (11)

where
$\bm{d}^{*}_{c}=[d^{*}_{i}]_{i\in\mathcal{C}},\bm{\lambda}_{c}=[\lambda_{i}]_{i\in\mathcal{C}}$
are the vectors of charging demands and prices at the charging stations,
respectively. Once such a function is obtained, the bilinear term
$\sum_{i\in\mathcal{C}}\lambda_{i}\cdot d_{i}^{*}$ in the objective can be
transformed to a term regarding the charging prices $\bm{\lambda}_{c}$, i.e.,
$(\bm{\lambda}_{c})^{\top}\pi(\bm{\lambda}_{c})$. As such, IDSO can solve its
own problem in (9) to determine the electricity price, without the need of
knowing the detailed decision-making problem of ITSO. Therefore, deriving the
function of $\pi(\bm{\lambda}_{c})$ in (11) is the key step. The detailed
solution process will be discussed in the next subsection.

###  III-B Derivation of Charging Demand Function

The goal of this subsection is to derive the charging demand function between
the vector of electricity price $\bm{\lambda}_{c}$ and the optimal charging
demand $\bm{d}_{c}^{*}$ by multiparametric programming theory [23, 24]. Before
going into the details, we first rewrite the decision problem in (7) into a
compact form of QP. Let
$\bm{m}\in\mathbb{R}^{|\mathcal{W}|}=[m_{w}]_{m\in\mathcal{W}}$ be the vector
formed by the demands of O-D pairs, and
$\bm{E}\in\mathbb{R}^{|\mathcal{W}|\times\sum_{w\in\mathcal{W}}|\mathcal{R}_{w}|}$
be the coefficient matrix, we rewrite (7c) by $\bm{E}\bm{f}=\bm{m}$. Let
$\bm{A}\in\mathbb{R}^{|\bm{A}^{e}|\times\sum_{w\in\mathcal{W}}|\mathcal{R}_{w}|}=[\bm{A}_{1},...,\bm{A}_{|\mathcal{W}|}]$
denote the matrix formed by the horizontal stack of the link matrix of each
O-D pair in the set $\mathcal{W}$. Therefore, (7d) can be rewritten as
$\bm{A}\bm{f}=\bm{\xi}$. Denote
$\bm{G}\in\mathbb{R}^{2|\mathcal{A}^{e}|\times|\mathcal{A}^{e}|}$ as the
matrix formed by the vertical stack of the negative and positive values of the
identity matrix $\bm{I}_{|\mathcal{A}^{e}|}$.
$\bm{h}\in\mathbb{R}^{2|\mathcal{A}^{e}|}$ is the vertical stack of an all-
zero vector and $\overline{\bm{\xi}}$. The constraint (7e) can be rewritten as
$\bm{G}\bm{\xi}\leq\bm{h}$. To sum up, the compact QP form of (7) follows,

$\displaystyle\mathop{\min}_{\bm{\xi},\bm{f}}\frac{1}{2}\bm{\xi}^{\top}\bm{Q}\bm{\xi}+\bm{q}^{\top}\bm{\xi}$
(12a) $\displaystyle s.t.\ \bm{E}\bm{f}=\bm{m}:\bm{\psi}$ (12b)
$\displaystyle\quad\ \ \bm{G}\bm{\xi}\leq\bm{h}:\bm{\phi}$ (12c)
$\displaystyle\quad\ \ \bm{A}\bm{f}=\bm{\xi}:\bm{\delta},$ (12d)

where $\bm{q}$ is a column vector, whose first $|\mathcal{A}|$ elements equal
$\gamma\cdot\xi_{a}^{0},\forall a\in\mathcal{A}$, and the last $|\mathcal{C}|$
elements equal
$\gamma(\xi_{a}^{0}+e_{a}/\rho_{a})+e_{a}\cdot\lambda_{i},\forall
a\in\cup_{i\in\mathcal{C}}l_{i}^{c}$, while the remaining $|\mathcal{C}|$
elements equal zero. The detailed form of matrices
$\bm{Q},\bm{E},\bm{G},\bm{h}$ are given in Appendix A for completeness. The
charging prices $\bm{\lambda}_{c}$ are in the last $|\mathcal{C}|$ elements of
the coefficient $\bm{q}$ which lies in the linear part of the objective. The
charging demand is related to the flow $\xi_{a}$ on the virtual arc
$l_{i}^{c},\forall i\in\mathcal{C}$. We now proceed with the following
proposition to characterize the _local_ relationship between the price
$\bm{\lambda}_{c}$ and the optimal solutions.

######  Proposition 1.

Consider the convex QP problem (12). Let
$\bm{\xi}_{i},\bm{f}_{i},\bm{\psi}_{i},\bm{\phi}_{i},\bm{\delta}_{i}$ denote
the local mapping between the charging price $\bm{\lambda}_{c}$ and the
optimal solutions of primal and dual variables,
$\tilde{\bm{\xi}},\tilde{\bm{f}},\tilde{\bm{\psi}},\tilde{\bm{\phi}},\tilde{\bm{\delta}}$
denote the optimal primal and dual solutions under the given point
$\tilde{\bm{\lambda}}_{c}^{i}$. The calculation of
$\bm{\xi}_{i},\bm{f}_{i},\bm{\psi}_{i},\bm{\phi}_{i},\bm{\delta}_{i}$, in the
neighborhood of $\tilde{\bm{\lambda}}_{c}^{i}$, is given by,

$\begin{bmatrix}\bm{\xi}_{i}\\\ \bm{f}_{i}\\\ \bm{\psi}_{i}\\\
\bm{\phi}_{i}\\\
\bm{\delta}_{i}\end{bmatrix}=\begin{bmatrix}\tilde{\bm{\xi}}\\\
\tilde{\bm{f}}\\\ \tilde{\bm{\psi}}\\\ \tilde{\bm{\phi}}\\\
\tilde{\bm{\delta}}\end{bmatrix}+(\bm{M}_{0})^{-1}\bm{N}_{0}(\bm{\lambda}_{c}-\tilde{\bm{\lambda}}_{c}^{i}),$
(13)

where
$\bm{N}_{0}\in\mathbb{R}^{(4|\mathcal{A}^{e}|+\sum_{w\in\mathcal{W}}|\mathcal{R}_{w}|+|\mathcal{W}|)\times|\mathcal{C}|}$
and
$\bm{M}_{0}\in\mathbb{R}^{(4|\mathcal{A}^{e}|+\sum_{w\in\mathcal{W}}|\mathcal{R}_{w}|+|\mathcal{W}|)\times(4|\mathcal{A}^{e}|+\sum_{w\in\mathcal{W}}|\mathcal{R}_{w}|+|\mathcal{W}|)}$
that are defined as follows,

$\displaystyle\bm{N}_{0}=\begin{bmatrix}\bm{O}&-\bm{J}&\bm{O}&\bm{O}&\bm{O}&\bm{O}\end{bmatrix}^{\top};$
$\displaystyle\bm{M}_{0}=\begin{bmatrix}\bm{Q}&\bm{O}&\bm{O}&\bm{G}^{\top}&-\bm{I}_{|\mathcal{A}^{e}|}\\\
\bm{O}&\bm{O}&\bm{E}^{\top}&\bm{O}&\bm{A}^{\top}\\\
\bm{O}&\bm{E}&\bm{O}&\bm{O}&\bm{O}\\\
-\bm{I}_{|\mathcal{A}^{e}|}&\bm{A}&\bm{O}&\bm{O}&\bm{O}\\\
D(\tilde{\bm{\phi}})\bm{G}&\bm{O}&\bm{O}&D(\bm{G}\tilde{\bm{\xi}}-\bm{h})&\bm{O}\end{bmatrix}$

where $\bm{O}$ is an all-zero matrix, and
$\bm{J}\in\mathbb{R}^{|\mathcal{C}|\times|\mathcal{C}|}$ is a diagonal matrix,
whose diagonal element equals the average charging demand $e_{a}$, and the
operation $D(\cdot)$ creates a diagonal matrix from a vector.

The proof of Proposition 1 is given in Appendix B. Once (13) is obtained, we
can use the mapping of $\bm{\xi}_{i}$ to obtain $\pi_{i}(\bm{\lambda}_{c})$,
which is the local relationship between the optimal charging demand
$\bm{d}_{c}^{*}$ and price $\bm{\lambda}_{c}$, i.e.,

$\bm{d}_{c}^{*}=[\bm{O}\in\mathbb{R}^{|C|\times(|A|+|C|)},\bm{J}]\bm{\xi}_{i},$
(14)

where $\bm{O}$ is an all-zero matrix with the size of $|C|\times(|A|+|C|)$.

For a given price point $\tilde{\bm{\lambda}}_{c}^{i}$, in the neighborhood of
$\tilde{\bm{\lambda}}_{c}^{i}$, where (13) and (14) hold, a critical region
$\mathcal{B}_{i}(\tilde{\bm{\lambda}}_{c}^{i})$ is defined that is related to
a specific combination of the active constraints. Let
$\hat{\bm{G}},\hat{\bm{h}}$ denote the coefficients of inactive constraints in
(12c), and $\hat{\bm{\phi}}_{i}$ be part of the $\bm{\phi}_{i}$ corresponding
the active inequalities. By substituting $\bm{\xi}_{i}$ into the inactive
constraints, feasibility is ensured. The optimality is ensured by
$\hat{\bm{\phi}}_{i}\geq 0$. The critical region of
$\mathcal{B}_{i}(\tilde{\bm{\lambda}}_{c}^{i})$ is defined as

$\begin{split}&\tilde{\mathcal{B}}_{i}(\tilde{\bm{\lambda}}_{c}^{i})=\\{\bm{E}\bm{f}_{i}=\bm{m},\hat{\bm{G}}\bm{\xi}_{i}\leq\hat{\bm{h}},\bm{A}\bm{f}_{i}=\bm{\xi}_{i},\hat{\bm{\phi}}_{i}\geq\bm{0},\Lambda\\}\\\
&\mathcal{B}_{i}(\tilde{\bm{\lambda}}_{c}^{i})=\triangle\\{\tilde{\mathcal{B}}_{i}(\tilde{\bm{\lambda}}_{c}^{i})\\},\end{split}$
(15)

where $\Lambda$ denotes the given initial convex set of all feasible vectors
of charging prices, and $\triangle\\{\cdot\\}$ is defined as an operator which
removes redundant constraints.

Let
$\\{\mathcal{B}_{1}(\tilde{\bm{\lambda}}_{c}^{1}),...,\mathcal{B}_{N}(\tilde{\bm{\lambda}}_{c}^{N})\\}$
be the set of all critical regions, then

$\displaystyle\mathcal{B}_{i}(\tilde{\bm{\lambda}}_{c}^{i})\cap\mathcal{B}_{j}(\tilde{\bm{\lambda}}_{c}^{j})=\emptyset,\forall
i\neq j$ (16)
$\displaystyle\cup_{i=1}^{N}\mathcal{B}_{i}(\tilde{\bm{\lambda}}_{c}^{i})=\Lambda$

The goal is to compute all the critical regions via (15), and the
corresponding local policy by (14). The exploration begins from each facet of
the initial critical region, and all adjacent critical regions intersecting
the facet are enumerated. An illustration of such process is shown in Fig. 3.
In turn, each of these adjacent critical regions is considered and so on,
until there are no more regions to be discovered. The main steps are
summarized in Algorithm 1.

Figure 3:  The illustration of the enumeration of the critical regions, where
$\mathcal{B}_{1}(\tilde{\bm{\lambda}}_{c}^{1})$ is an initial critical region
whose one of the facets is $f$. The adjacent critical regions along this facet
form the set $\mathcal{A}$. Algorithm 1 The derivation of charging demand
function $\pi(\bm{\lambda}_{c})$ via multiparametric programming

1:The initial set of all feasible vectors of charging price $\Lambda$, the
initial given point $\tilde{\lambda}_{c}^{1}$, and the corresponding critical
region $\mathcal{B}_{1}(\tilde{\bm{\lambda}}_{c}^{1})$;

2:Calculate the local function by (13),(14) and initialize the set of
unexplored region $\mathcal{L}_{un}$ and the set of discovered region
$\mathcal{L}_{ds}$ by
$\mathcal{L}_{un}\longleftarrow\\{\mathcal{B}_{1}(\tilde{\bm{\lambda}}_{c}^{1})\\}$,
$\mathcal{L}_{ds}\longleftarrow\\{\mathcal{B}_{1}(\tilde{\bm{\lambda}}_{c}^{1})\\}$.

3:while $\mathcal{L}_{un}$ is not empty do.

4: Select and remove any region $\mathcal{R}$ in the set $\mathcal{L}_{un}$

5: for Each facet $f$ of $\mathcal{R}$ do

6: Compute set of critical regions $\mathcal{A}$ that along facet $f$,
$\mathcal{L}_{un}\longleftarrow\mathcal{L}_{un}\cup(\mathcal{A}\backslash\mathcal{L}_{ds})$,
$\mathcal{L}_{ds}\longleftarrow\mathcal{L}_{ds}\cup\mathcal{A}$

7: Calculate the local function by (13),(14) for each critical region in the
set $\mathcal{A}$

8: end for

9:end while

10:Collect all the solutions of local function and unify the critical regions
$\mathcal{L}_{ds}$

In a nutshell, the charging demand function $\pi(\bm{\lambda}_{c})$ on the
region of $\Lambda$ is a piecewise linear function, where the local function
$\pi_{i}(\bm{\lambda}_{c})$ is affine on the corresponding critical region
$\mathcal{B}_{i}(\tilde{\bm{\lambda}}_{c}^{i})$. With this representation of
$\bm{d}_{c}^{*}=\pi(\bm{\lambda}_{c})$, the price setting problem of IDSO can
be solved separately on each critical region with the local policy. The
solution process is detailed in the next subsection.

###  III-C Solution Strategy

With the charging demand function $\pi(\bm{\lambda}_{c})$, and the set of
critical regions, we use the half-space representation of a polytope and
rewrite each critical region $\mathcal{B}_{i}(\tilde{\bm{\lambda}}_{c}^{i})$
as $\bm{R}_{i}\cdot\bm{\lambda}_{c}\leq\bm{r}_{i}$, whose corresponding local
function is $\pi_{i}(\bm{\lambda}_{c})$. Thus, the problem of IDSO becomes a
QP problem in each local region, which is given by

$\displaystyle\mathop{\min}_{\bm{\overline{\tau}},\bm{\underline{\tau}},\bm{\overline{\mu}},\bm{\underline{\mu}},\bm{\lambda},\bm{\underline{\eta}},\bm{\overline{\eta}}}\sum_{i\in\mathcal{N}}\overline{\mu_{i}}\cdot\overline{v_{i}}-\underline{\mu_{i}}\cdot\underline{v_{i}}+\overline{\tau_{i}}\cdot\overline{g_{i}}+$
(17a)
$\displaystyle\qquad\sum_{i\in\mathcal{N}}\sum_{j\in\Omega_{i}}(f_{ij}\overline{\eta_{ij}}-f_{ij}\underline{\eta_{ij}})-\sum_{i\in\mathcal{N}}\lambda_{i}\cdot
l_{i}-(\bm{\lambda}_{c})^{\top}\cdot\pi_{i}(\bm{\lambda}_{c})$ $\displaystyle
s.t.\eqref{9(b)},\eqref{9(c)},\eqref{9(d)},\eqref{9(e)}$ $\displaystyle\quad\
\bm{R}_{i}\cdot\bm{\lambda}_{c}\leq\bm{r}_{i}$ (17b)

We solve the QP problem in (17) for every local critical region and return the
optimal solution of the problem, whose optimal objective is the smallest among
all the problems. The corresponding optimal solution of $\bm{\lambda}^{*}$
represents the optimal electricity price across all power network nodes, part
of which is the optimal charging price, i.e.,
$\bm{\lambda}_{c}^{*}=[\lambda_{i}^{*}]_{i\in\mathcal{C}}$.

Remark 2: Here, (17) is solved at single-period. For example, it can be solved
every 15 minutes, or one hour. We note that the proposed approach can be
extended to a multi-period problem, which minimizes the charging expenses,
i.e., the last term of 17(a), over time. For that, the charging demand at each
period is also replaced with the charging demand function at that period,
which is affine on its critical region. In this line, the total charging
expense over time is a quadratic term regarding the charging prices as well.
The resulting problem is still a QP and can be solved by many off-the-shelf
solvers.

Also, the OPF problem at the upper-level is replaced with its dual problem in
(17), which requires the OPF problem to be a convex one to make the strong
duality hold. In this line, more realistic distribution network models can be
considered, such as the one considering multiphase couplings and unbalanced
loads, once there is a convex relaxation for it.

###  III-D Relationship with Other Modeling Perspectives

Here, we discuss the relationship between this approach and others modeling
the interaction between IDSO and ITSO from collaborative [11, 13] and game-
theoretical [14] perspectives.

1) Collaborative Perspective: Ref. [11, 13] proposed to model the power and
traffic network interaction by an joint optimization approach, i.e.,

$\displaystyle\mathop{\min}\
\eqref{8(a)}+\eqref{7(a)}-\sum_{i\in\mathcal{C}}\lambda_{i}\cdot d_{i}$ (18)
$\displaystyle s.t.\text{IDSO constraints \eqref{8(b)}-\eqref{8(e)}}$ ITSO
constraints (7b)-(7e)

2) Game-theoretical Perspective: Assuming generators and EVs are price-takers
(as shown in Fig. 1), [14] proposed that the interaction between the power and
traffic networks can be modeled as a game with generators and EVs acting as
the players to maximize their own profits. The equilibrium in the coupled
power and traffic networks can be described by combining the KKT conditions of
the two problems, i.e.,

KKT conditions of IDSO problem (8) (19a) KKT conditions of ITSO problem (7)
(19b)

By proving the KKT conditions of (18) have the same form as that of (19), [14]
showed that the game-theoretical and collaborative perspectives are
equivalent. In this work, we will show that our bilevel program has the same
KKT conditions as that of (19), and therefore is equivalent to the model from
the game-theoretical perspective. The proof can be found in Appendix C.
Leveraging the relationship established in [14] that the game-theoretical and
collaborative perspectives are equivalent, our model is also equivalent to the
collaborative perspective. The relationship is summarized in Fig. 4.

Figure 4:  The relationship between bilevel, collaborative and game-
theoretical perspectives.

##  IV Case Study

This section presents the results of testing the effectiveness of the proposed
approach on a test system [25] consisting of the traffic network of Sioux
Falls city and the 33-bus distribution network [26], as shown in Fig. 5. The
coupled networks include four charging stations with specified flow capacity
$\overline{\xi_{a}}$ and charging demand capacity, as outlined in Table I. In
the traffic network, the O-D pairs of Nodes #1-#13 and Nodes #1-#20 are
considered for EVs and Nodes #2-#13 and Nodes #2-#20 are considered for
regular vehicles (RVs). Although our approach is applicable to the
consideration of more O-D pairs, we assume that there are only four pairs
having the traffic demand, while others don’t. The summary of routes in each
O-D pair can be found in [27]. The shape parameters in (5) and the cost of
unit time are given as $\xi_{a}^{0}=0,R_{a}=10^{4},\gamma=10^{3}$. The average
charging demand $e_{a}$ of an individual vehicle is set as 12 kWh, and the
charging rate $\rho_{a}$ is 200 kW. In the distribution network, the lower and
upper bounds of voltage are set as 0.9 p.u. and 1.1 p.u., respectively. There
are 5 distributed generators (G1 - G5) connected to the buses
#4,#13,#16,#19,#29, as shown in Fig 5. The capacity and the cost coefficient
of each generator are given in Table II. Also, the fixed loads of IEEE 33
network are scaled up to ten times.

Figure 5:  The Sioux Falls city traffic network and 33-bus distribution network coupled by charging stations (CSs). Table I:  The flow capacity and charging demand capacity for the four EV charging stations. |  Charging station 1 |  Charging station 2 |  Charging station 3 |  Charging station 4
---|---|---|---|---
Flow capacity (number of vehicles) | 100 | 250 | 250 | 100
Charging demand capacity (kWh) | 1200 | 3000 | 3000 | 1200
Table II: Generation capacity and cost for the five generators. |  G1 |  G2 |  G3 |  G4 |  G5
---|---|---|---|---|---
Generation capacity (kWh) | $10^{4}$ | $10^{4}$ | $3\times 10^{4}$ | $10^{4}$ | $10^{4}$
Cost coefficient ($/kWh) | 0.8 | 0.6 | 0.5 | 0.7 | 0.4

We conduct the experiments from three angles: (1) To investigate the impact of
O-D pair traffic demand on the charging demand; (2) To show the effectiveness
of the charging price on affecting the charging load patterns; (3) To show the
bilevel model is robust against the traffic demand forecast error to some
extent. All simulations are implemented on the laptop with Intel®CoreTM
i5-10210U 1.6 GHz CPU, and 8.00 GM RAM and based on the Multiparametric
Toolbox [28].

###  IV-A Impacts of Traffic Demands on Charging Demands and Prices

We first assign each O-D pair with equal traffic demand and test the proposed
method under different levels of it, namely 100, 200, and 300 vehicles. The
cost of IDSO, the number of critical regions, and the corresponding
computation time are shown in Table III. Unsurprisingly, the increase in
traffic demand results in an increase of charging demand, and therefore, the
operational cost of IDSO increases correspondingly. In addition, we observe
that as a parameter in the equality constraint (7c) of the ITSO decision
problem, the value of the traffic demand has impact on the traffic pattern and
results in different numbers of critical regions. Concretely, the number of
critical regions reduces as the traffic demand becomes larger. With the larger
demand, more inequality constraints regarding the arc flows are binding, which
results in a smaller combination of operation patterns, and therefore a
smaller number of critical regions. Since the proposed approach iterates over
all critical regions, the computation time increases when the number of
critical regions is larger. Furthermore, we also calculate the cost of IDSO by
solving the joint optimization in (18), whose results are listed in the last
row of Table III, and the cost of IDSO by solving (10) using single-level
reduction via KKT conditions, whose computation time is also reported.
Unsurprisingly, the costs of the IDSO by solving (18) and our bilevel program
in (10) are the same, as they are equivalent problems shown in Section III-D.
However, (18) requires the whole decision model of the ITSO, while our
approach doesn’t and still achieves optimality. For the comparison candidate
using single-level reduction via KKT conditions, it is transformed to a mixed
integer program. Therefore, its computation time is much longer than the
propose approach, which shows that the proposed approach is more
computationally efficient.

Table III:  Under different traffic demand levels, the cost of IDSO solved by the proposed approach for the bilevel program in (10), the number of critical regions, the computation time of the proposed approach, the computation time of solving (10) by single-level reduction via KKT conditions, the cost of IDSO of solving (10) by single-level reduction via KKT conditions, and the cost of IDSO by solving the joint optimization in (18). | $m_{w}=100$ | $m_{w}=200$ | $m_{w}=300$
---|---|---|---
The cost of IDSO ($) by the proposed approach | 19967 | 21478 | 23068
The number of the critical regions | 93 | 8 | 3
Computation time of the proposed approach (s) | 14 | 0.86 | 0.32
Computation time of solving (10) via KKT conditions (s) | 145 | 5640 | 44
The cost of IDSO ($) of solving (10) via KKT conditions | 19967 | 21478 | 23068
The cost of IDSO ($) by solving joint optimization [11, 13] | 19967 | 21478 | 23068

Figure 6: The charging demand at four charging stations under different
traffic demand levels.

We show the electricity price at each charging station and the charging demand
under different traffic demand levels in the first row of Table IV and Fig. 6,
respectively. Since the charging price under different traffic demand levels
remains the same, we show the results in a single line. As the line connecting
node #29 and node #25 is congested, the charging prices are different in
different charging stations. It is shown in Fig. 6 that the charging demands
in the traffic network can be affected by the charging price. With a lower
charging price, the charging demand in the charging station (CS) #2 is always
larger than the charging demand in CS #3, across different traffic demand
levels. Also, we observe that the charging price is not the only factor
affecting the charging demands. The traffic network structure, such as the
constituent routes of each O-D pair, has a significant impact as well. For
instance, among the 6 routes of the EV O-D pair Node #1-#13, 4 routes pass CS
#1. Therefore, it carries most of the charging demand of O-D pair Node #1-#13.
Hence, when the traffic demand of each O-D pair is 100 vehicles, the charging
demand in CS #1 is the largest. When the traffic demand of each O-D pair
increases to 200 and 300 vehicles, the demand in CS #1 remains operating at
the maximum capacity.

Table IV: The charging prices ($/kWh) at four charging stations after the decrease in generation cost ($/kWh) of G3. Generation cost of G3 |  Charging station 1 |  Charging station 2 |  Charging station 3 |  Charging station 4
---|---|---|---|---
0.5 | 0.7 | 0.62 | 0.71 | 0.5
0.3 | 0.8 | 0.6 | 0.82 | 0.3
0.1 | 0.74 | 0.48 | 0.75 | 0.1
0 | 0.78 | 0.47 | 0.8 | 0

In addition, we compare the proposed approach with a baseline charging
strategy in which all EVs choose to charge at the charging stations with the
lowest charging price, specifically at CS #2 and CS #4. We compare the costs
of IDSO, ITSO, and the operation cost of both networks defined in the
objective of (18), which is the negative value of the social welfare [11]. The
results of this comparison when the traffic demand is set to 100 vehicles are
presented in Table V.

Table V: The costs of IDSO, ITSO, and two networks’ operation cost under the proposed approach and the baseline approach |  IDSO cost ($) |  ITSO cost ($) |  Two networks’ operation cost ($)
---|---|---|---
The proposed approach | 19967 | 43258 | 63225
Baseline: charging EVs at the lowest LMP nodes | 19794 | 48294 | 68088

Although the baseline charging strategy results in the lower operating cost of
IDSO, the cost of ITSO is larger than the proposed approach, as the traffic
time on arcs increases when all EVs travel to the same stations with lower
prices for charging. The two networks’ operation cost of the baseline approach
increases by 7.6%, compared with the proposed bilevel optimization approach.

Moreover, to show the proposed approach is applicable to large scale
distribution network, we apply the proposed approach on the 85-bus [29] and
136-bus [30] distribution networks, as well. The generators and the charging
stations are connected to the nodes with the same indexes as the one shown in
Fig. 5. The traffic demand of 100 vehicles is considered. The charging demand
function is a piecewise linear function defined on 93 critical regions. We run
the program for ten times under 33-bus, 85-bus and 136-bus distribution
networks. The average computation time and its standard error are reported in
Table VI. The computation time increases as the scale of the network becomes
larger, as a larger number of variables and constraints are considered.
However, even for the case concerning the 136-bus network, the computation
time is acceptable for our optimization problem, which is often conducted on a
hourly basis.

Table VI:  The average and standard error of the computation time in ten runs under 33-bus, 85-bus and 136-bus distribution networks. | 33-bus network | 85-bus network | 136-bus network
---|---|---|---
Average | 14 s | 63 s | 160 s
Standard error | 2 s | 1 s | 6 s

###  IV-B Impacts of Electricity Prices on Charging Demands

To demonstrate the impact of charging prices on charging demand, we
investigate the effect of generation costs on charging demand patterns. We can
view the cost coefficient of the generator connected to a node as an
equivalent cost coefficient of a generation mix. Therefore, with the
increasing penetration of renewable sources, the cost of the generation mix
decreases. For this scenario, we test different values of the cost coefficient
of the generator G3 connected to node #29, where CS #4 is located, while
keeping the other coefficients the same. Specifically, we set the generation
cost of G3 as 0.5, 0.3, 0.1, and 0 $/kWh, respectively, and set the traffic
demand of each O-D pair to 100 vehicles. The charging prices are outlined in
Table IV. As the cost coefficient of G3 decreases, the price difference
between CS #1 and CS #4 becomes larger.

Fig. 7 illustrates the charging demands under different generation costs of
G3. The impact of charging prices on charging demand is significant. When the
cost coefficient of G3 is 0.5 $/kWh and 0.3 $/kWh, the charging demand at CS
#1 is higher than that at CS #4. However, as the cost coefficient of G3
decreases (and thus the price at CS #4 decreases), the charging demand at CS
#4 increases and eventually surpasses the charging demand at CS #1. This shows
that electricity prices can effectively alter the spatial distribution of
charging demand and incentivize EVs to charge at locations with lower
electricity prices, which may coincide with locations with plentiful renewable
sources and the marginal generation cost is low.

Figure 7: The charging demand at four charging stations under different
generation costs of G3.

###  IV-C Investigation on Inaccurate Traffic Demand Forecasts

It is widely acknowledged that obtaining high-quality estimates or forecasts
of traffic demand is challenging, and the actual traffic demand $m_{w},\forall
w\in\mathcal{W}$ may deviate from the forecast value. In this subsection, we
test the robustness of the bilevel model against the forecast error.
Specifically, we assume the accurate values of the traffic demand for each O-D
pair are all equal to 300 vehicles. We generate 30 samples of the possible
forecasts which are the integers ranging from 285 to 315 vehicles ($\pm 5\%$
deviations from the ground-truth traffic demand). The percentage of the
deviation between the IDSO’s costs under the forecast and the realization is
then calculated and shown in the histogram Fig. 8. The maximum cost deviation
percentage is less than 1.5%, which shows the IDSO’s costs under the forecast
and the realization are close. Since the deviation of the traffic demand is in
a relatively small range that doesn’t change the charging patterns, the
charging demand still reaches the maximum capacity of CS #1 and CS #4, and
only the demands in the CS #2 and CS #3 are different under different forecast
values. Also, as the charging price under different levels of charging demands
remains the same, the different charging demands in the two stations have a
small impact on the cost of IDSO. Therefore, the results show that the bilevel
model is robust against the forecast error of traffic demand to some extent.

Figure 8: The distribution of the deviation between the IDSO’s cost under the
forecast and the realization.

##  V Conclusion

In this work, we propose to manage EV charging via a bilevel optimization
approach, where the IDSO determines the charging price at the upper level, and
the ITSO decides the optimal charging power for the given price at the lower
level. The charging demand function depicting the theoretical relationship
between the charging price and the optimal charging demand is derived, which
is piecewise linear. With the charging demand function, the IDSO can solve the
optimal EV management problem without the need for a detailed model of ITSO,
while still considering the coupling of the power-traffic networks. In
addition, with the bilinear term replaced with a piecewise quadratic term, we
are able to solve the bilevel optimization efficiently with optimality
guarantees. Case studies reveal that the different levels of traffic demand
can affect the charging pattern in the charging stations. Also, charging
prices can effectively change the charging demand, where the charging stations
with lower price shoulder more charging demand than stations with higher
prices.

In the future, it would be interesting to extend our approach to multi-period
EV charging optimization. Besides, the single-phase distribution network model
is used here for simplicity. It is still required to explore more realistic
distribution network model, considering the multiphase couplings and
unbalanced loads. Additionally, our approach provides a novel solution to the
bilevel program by leveraging the relationship between the upper-level
variable and the optimal solution of the lower-level problem, allowing for the
upper-level problem to be solved independently with an optimality guarantee.
It is also interesting to apply the proposed approach to other tasks modeled
as a bilevel program.

## Acknowledgement

The authors would like to appreciate the Shanghai Jiao Tong University Grants.

## References

  * [1] L. Xie, C. Singh, S. K. Mitter, M. A. Dahleh, and S. S. Oren, “Toward carbon-neutral electricity and mobility: Is the grid infrastructure ready?” _Joule_ , vol. 5, no. 8, pp. 1908–1913, 2021.
  * [2] “Global ev sales reach 10 percent of new cars sold in 2022,” https://www.teslarati.com/electric-vehicle-sales-2022/.
  * [3] M. Muratori, “Impact of uncoordinated plug-in electric vehicle charging on residential power demand,” _Nature Energy_ , vol. 3, no. 3, pp. 193–201, 2018.
  * [4] R. El Helou, S. Sivaranjani, D. Kalathil, A. Schaper, and L. Xie, “The impact of heavy-duty vehicle electrification on large power grids: A synthetic texas case study,” _Advances in Applied Energy_ , vol. 6, p. 100093, 2022.
  * [5] N. Xu and C. Chung, “Challenges in future competition of electric vehicle charging management and solutions,” _IEEE Transactions on Smart Grid_ , vol. 6, no. 3, pp. 1323–1331, 2014.
  * [6] Z. Yi, Y. Xu, W. Gu, L. Yang, and H. Sun, “Aggregate operation model for numerous small-capacity distributed energy resources considering uncertainty,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 5, pp. 4208–4224, 2021.
  * [7] G. Lacey, G. Putrus, and E. Bentley, “Smart ev charging schedules: supporting the grid and protecting battery life,” _IET Electrical Systems in Transportation_ , vol. 7, no. 1, pp. 84–91, 2017.
  * [8] D. Lauinger, F. Vuille, and D. Kuhn, “Reliable frequency regulation through vehicle-to-grid: From eu legislation to robust optimization,” _arXiv preprint arXiv:2005.06042_ , 2020.
  * [9] S. I. Vagropoulos and A. G. Bakirtzis, “Optimal bidding strategy for electric vehicle aggregators in electricity markets,” _IEEE Transactions on Power Systems_ , vol. 28, no. 4, pp. 4031–4041, 2013.
  * [10] R. Li, Q. Wu, and S. S. Oren, “Distribution locational marginal pricing for optimal electric vehicle charging management,” _IEEE Transactions on Power Systems_ , vol. 29, no. 1, pp. 203–211, 2013.
  * [11] M. Alizadeh, H.-T. Wai, M. Chowdhury, A. Goldsmith, A. Scaglione, and T. Javidi, “Optimal pricing to manage electric vehicles in coupled power and transportation networks,” _IEEE Transactions on Control of Network Systems_ , vol. 4, no. 4, pp. 863–875, 2017.
  * [12] N. Chen, J. Ma, M. Wang, and X. Shen, “Two-tier energy compensation framework based on mobile vehicular electric storage,” _IEEE Transactions on Vehicular Technology_ , vol. 67, no. 12, pp. 11 719–11 732, 2018.
  * [13] F. Rossi, R. Iglesias, M. Alizadeh, and M. Pavone, “On the interaction between autonomous mobility-on-demand systems and the power network: Models and coordination algorithms,” _IEEE Transactions on Control of Network Systems_ , vol. 7, no. 1, pp. 384–397, 2020.
  * [14] Z. Zhou, S. J. Moura, H. Zhang, X. Zhang, Q. Guo, and H. Sun, “Power-traffic network equilibrium incorporating behavioral theory: A potential game perspective,” _Applied Energy_ , vol. 289, p. 116703, 2021.
  * [15] Y. Zheng and Qingshan, “A review of distributed optimization: Problems, models and algorithms,” _Neurocomputing_ , vol. 483, pp. 446–459, 2022.
  * [16] A. Sinha, P. Malo, and K. Deb, “A review on bilevel optimization: From classical to evolutionary approaches and applications,” _IEEE Transactions on Evolutionary Computation_ , vol. 22, no. 2, pp. 276–295, 2018\.
  * [17] D. S. Kirschen and G. Strbac, _Fundamentals of power system economics_. John Wiley & Sons, 2018.
  * [18] P. Tøndel, T. A. Johansen, and A. Bemporad, “An algorithm for multi-parametric quadratic programming and explicit mpc solutions,” _Automatica_ , vol. 39, no. 3, pp. 489–497, 2003.
  * [19] A. Grancharova and T. A. Johansen, _Multi-parametric Programming_. Berlin, Heidelberg: Springer Berlin Heidelberg, 2012, pp. 1–37.
  * [20] K. Davidson, “A flow travel time relationship for use in transportation planning,” in _Australian Road Research Board (ARRB) Conference, 3rd, 1966, Sydney_ , vol. 3, no. 1, 1966.
  * [21] Y. Bian, N. Zheng, Y. Zheng, B. Xu, and Y. Shi, “Predicting strategic energy storage behaviors,” _arXiv preprint arXiv:2306.11872_ , 2023.
  * [22] Z. Yi, Y. Xu, W. Gu, and W. Wu, “A multi-time-scale economic scheduling strategy for virtual power plant based on deferrable loads aggregation and disaggregation,” _IEEE Transactions on Sustainable Energy_ , vol. 11, no. 3, pp. 1332–1346, 2020.
  * [23] C. N. Jones and M. Morrari, “Multiparametric linear complementarity problems,” in _Proceedings of the 45th IEEE Conference on Decision and Control_ , 2006, pp. 5687–5692.
  * [24] T. Gal, _Postoptimal Analyses, Parametric Programming, and Related Topics_. Berlin, New York: De Gruyter, 2010\.
  * [25] X. Liu, C. B. Soh, T. Zhao, and P. Wang, “Stochastic scheduling of mobile energy storage in coupled distribution and transportation networks for conversion capacity enhancement,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 1, pp. 117–130, 2021.
  * [26] M. Baran and F. Wu, “Network reconfiguration in distribution systems for loss reduction and load balancing,” _IEEE Transactions on Power Delivery_ , vol. 4, no. 2, pp. 1401–1407, 1989.
  * [27] “The arc collection of routes in each o-d pair,” https://figshare.com/articles/figure/O-D_pair_data_xlsx/21762539.
  * [28] M. Herceg, M. Kvasnica, C. Jones, and M. Morari, “Multi-Parametric Toolbox 3.0,” in _Proc. of the European Control Conference_ , Zürich, Switzerland, July 17–19 2013, pp. 502–510, http://control.ee.ethz.ch/~mpt.
  * [29] D. Das, D. Kothari, and A. Kalam, “Simple and efficient method for load flow solution of radial distribution networks,” _International Journal of Electrical Power & Energy Systems_, vol. 17, no. 5, pp. 335–346, 1995.
  * [30] J. R. Mantovani, F. Casari, and R. A. Romero, “Reconfiguração de sistemas de distribuição radiais utilizando o critério de queda de tensão,” _Controle and Automacao_ , pp. 150–159, 2000.
  * [31] F. Alesiani, “Implicit bilevel optimization: Differentiating through bilevel optimization programming,” _arXiv preprint arXiv:2302.14473_ , 2023.

##  Appendix A The Coefficients of the Decision-making of ITSO in Compact QP
Form

$\bm{Q}$ is a diagonal matrix, whose the first $|\mathcal{A}|$ and the last
$|\mathcal{C}|$ principal diagonal elements equal $2\gamma/R_{a},\forall
a\in\mathcal{A}\cup(\cup_{i\in\mathcal{C}}l_{i}^{c})$, and the remaining ones
equal zero. That is, the diagonal elements are
$[\underbrace{\frac{2\gamma}{R_{1}},\ldots,\frac{2\gamma}{R_{|\mathcal{A}|}}}_{|\mathcal{A}|},\underbrace{0,\ldots,0}_{|\mathcal{C}|},\underbrace{\frac{2\gamma}{R_{l_{1}^{c}}},\ldots,\frac{2\gamma}{R_{l_{|\mathcal{C}|}^{c}}}}_{|\mathcal{C}|}]$.
Also,

$\displaystyle\bm{E}=\begin{bmatrix}\underbrace{1\ldots
1}_{|\mathcal{R}_{1}|}\underbrace{0\ldots 00\ldots
0}_{\sum_{i=2}^{|\mathcal{W}|}|\mathcal{R}_{i}|}\\\ \underbrace{0\ldots
0}_{|\mathcal{R}_{1}|}\underbrace{1\ldots
1}_{|\mathcal{R}_{2}|}\underbrace{0\ldots
0}_{\sum_{i=3}^{|\mathcal{W}|}|\mathcal{R}_{i}|}\\\ \ddots\\\ 0\ldots\ldots
0\underbrace{1\ldots
1}_{|\mathcal{R}_{|\mathcal{W}|}|}\end{bmatrix},\bm{G}=\begin{bmatrix}-\bm{I}_{|\mathcal{A}^{e}|}\\\
\bm{I}_{|\mathcal{A}^{e}|}\end{bmatrix},\bm{h}=\begin{bmatrix}\bm{0}_{|\mathcal{A}^{e}|}\\\
\overline{\bm{\xi}}\end{bmatrix}$

##  Appendix B Proof of Proposition 1

The Lagrangian of (12) is given by

$\begin{split}&L=\frac{1}{2}\bm{\xi}^{\top}\bm{Q}\bm{\xi}+\bm{q}^{\top}\bm{\xi}+\bm{\psi}^{\top}(\bm{E}\bm{f}-\bm{m})+\bm{\phi}^{\top}(\bm{G}\bm{\xi}-\bm{h})+\\\
&\qquad\qquad\quad\bm{\delta}^{\top}(\bm{A}\bm{f}-\bm{\xi})\end{split}$ (20)

When $\bm{\lambda}_{c}=\tilde{\bm{\lambda}}_{c}^{i}$, the optimal solution of
primal and dual variables becomes
$\tilde{\bm{\xi}},\tilde{\bm{f}},\tilde{\bm{\phi}},\tilde{\bm{\psi}},\tilde{\bm{\delta}}$.
The KKT conditions for stationarity, primal feasibility, and complementary
slackness are

$\displaystyle\bm{Q}\tilde{\bm{\xi}}+\bm{q}+\bm{G}^{\top}\tilde{\bm{\phi}}-\tilde{\bm{\delta}}=0$
(21a)
$\displaystyle\bm{E}^{\top}\tilde{\bm{\psi}}+\bm{A}^{\top}\tilde{\bm{\delta}}=0$
(21b) $\displaystyle\bm{E}\tilde{\bm{f}}-\bm{m}=0$ (21c)
$\displaystyle\bm{A}\tilde{\bm{f}}-\tilde{\bm{\xi}}=0$ (21d) $\displaystyle
D(\tilde{\bm{\phi}})(\bm{G}\tilde{\bm{\xi}}-\bm{h})=0,$ (21e)

where the operation $D(\cdot)$ creates a diagonal matrix from a vector. Taking
the differentials of these conditions gives the equations

$\displaystyle
d\bm{Q}\tilde{\bm{\xi}}+\bm{Q}d\bm{\xi}+d\bm{q}+d\bm{G}^{\top}\tilde{\bm{\phi}}+\bm{G}^{\top}d\bm{\phi}-\bm{I}_{|\mathcal{A}^{e}|}d\bm{\delta}=0$
(22a) $\displaystyle
d\bm{E}^{\top}\tilde{\bm{\psi}}+\bm{E}^{\top}d\bm{\psi}+d\bm{A}^{\top}\tilde{\bm{\delta}}+\bm{A}^{\top}d\bm{\delta}=0$
(22b) $\displaystyle d\bm{E}\tilde{\bm{f}}+\bm{E}d\bm{f}-d\bm{m}=0$ (22c)
$\displaystyle
d\bm{A}\tilde{\bm{f}}+\bm{A}d\bm{f}-\bm{I}_{|\mathcal{A}^{e}|}d\bm{\xi}=0$
(22d) $\displaystyle
D(\bm{G}\tilde{\bm{\xi}}-\bm{h})d\bm{\phi}+D(\tilde{\bm{\phi}})(d\bm{G}\tilde{\bm{\xi}}+\bm{G}d\bm{\xi}-d\bm{h})=0$
(22e)

We rewrite (22) into a compact matrix form

$\begin{split}&\begin{bmatrix}\bm{Q}&\bm{O}&\bm{O}&\bm{G}^{\top}&-\bm{I}_{|\mathcal{A}^{e}|}\\\
\bm{O}&\bm{O}&\bm{E}^{\top}&\bm{O}&\bm{A}^{\top}\\\
\bm{O}&\bm{E}&\bm{O}&\bm{O}&\bm{O}\\\
-\bm{I}_{|\mathcal{A}^{e}|}&\bm{A}&\bm{O}&\bm{O}&\bm{O}\\\
D(\tilde{\bm{\phi}})\bm{G}&\bm{O}&\bm{O}&D(\bm{G}\tilde{\bm{\xi}}-\bm{h})&\bm{O}\end{bmatrix}\\\
&\cdot\begin{bmatrix}d\bm{\xi}\\\ d\bm{f}\\\ d\bm{\psi}\\\ d\bm{\phi}\\\
d\bm{\delta}\end{bmatrix}=\begin{bmatrix}-d\bm{Q}\tilde{\bm{\xi}}-d\bm{q}-d\bm{G}^{\top}\tilde{\bm{\phi}}\\\
-d\bm{E}^{\top}\tilde{\bm{\psi}}-d\bm{A}^{\top}\tilde{\bm{\delta}}\\\
-d\bm{E}\tilde{\bm{f}}+d\bm{m}\\\ -d\bm{A}\tilde{\bm{f}}\\\
D(\tilde{\bm{\phi}})d\bm{h}-D(\tilde{\bm{\phi}})d\bm{G}\tilde{\bm{\xi}}\end{bmatrix}\end{split}$
(23)

The left-hand side coefficient matrix is the matrix $\bm{M}_{0}$. We wish to
compute the Jacobian
$\frac{\partial\bm{\xi}}{\partial\bm{\lambda}_{c}},\frac{\partial\bm{f}}{\partial\bm{\lambda}_{c}},\frac{\partial\bm{\psi}}{\partial\bm{\lambda}_{c}},\frac{\partial\bm{\phi}}{\partial\bm{\lambda}_{c}},\frac{\partial\bm{\delta}}{\partial\bm{\lambda}_{c}}$.
And take $\frac{\partial\bm{\xi}}{\partial\bm{\lambda}_{c}}$ for instance, it
can be calculated according to the chain rule,
i.e.,$\frac{\partial\bm{\xi}}{\partial\bm{\lambda}_{c}}=\frac{\partial\bm{\xi}}{\partial\bm{q}}\cdot\frac{\partial\bm{q}}{\partial\bm{\lambda}_{c}}$.
Therefore, we firstly calculate the Jacobian
$\frac{\partial\bm{\xi}}{\partial\bm{q}},\frac{\partial\bm{f}}{\partial\bm{q}},\frac{\partial\bm{\psi}}{\partial\bm{q}},\frac{\partial\bm{\phi}}{\partial\bm{q}},\frac{\partial\bm{\delta}}{\partial\bm{q}}$
by substituting $d\bm{q}=\bm{I}_{|\mathcal{A}^{e}|}$, and setting all other
differential terms in the right-hand side to zero. Therefore, the Jacobian is
given by

$\begin{split}&\begin{bmatrix}\frac{\partial\bm{\xi}}{\partial\bm{q}}&\frac{\partial\bm{f}}{\partial\bm{q}}&\frac{\partial\bm{\psi}}{\partial\bm{q}}&\frac{\partial\bm{\phi}}{\partial\bm{q}}&\frac{\partial\bm{\delta}}{\partial\bm{q}}\end{bmatrix}^{\top}\\\
&=\bm{M}_{0}^{-1}\cdot\begin{bmatrix}-\bm{I}_{|\mathcal{A}^{e}|}&\bm{O}&\bm{O}&\bm{O}&\bm{O}\end{bmatrix}^{\top}\end{split}$
(24)

The Jacobian of $\frac{\partial\bm{q}}{\partial\bm{\lambda}_{c}}$ equals

$\begin{bmatrix}\bm{O}&\bm{J}\end{bmatrix}^{\top}$ (25)

which is the stack of a zero matrix with the dimension of
$(|\mathcal{A}|+|\mathcal{C}|)\times|\mathcal{C}|$ and a diagonal matrix
$\bm{J}\in\mathbb{R}^{|\mathcal{C}|\times|\mathcal{C}|}$ whose principal
diagonal element equals the average charging demand $e_{a},\forall
a\in\cup_{i\in\mathcal{C}}l_{i}^{c}$. By applying the chain rule, the Jacobian
of
$\frac{\partial\bm{\xi}}{\partial\bm{\lambda}_{c}},\frac{\partial\bm{f}}{\partial\bm{\lambda}_{c}},\frac{\partial\bm{\psi}}{\partial\bm{\lambda}_{c}},\frac{\partial\bm{\phi}}{\partial\bm{\lambda}_{c}},\frac{\partial\bm{\delta}}{\partial\bm{\lambda}_{c}}$
is

$\bm{M}_{0}^{-1}\cdot\begin{bmatrix}-\bm{I}_{|\mathcal{A}^{e}|}&\bm{O}&\bm{O}&\bm{O}&\bm{O}\end{bmatrix}^{\top}\cdot\begin{bmatrix}\bm{O}&\bm{J}\end{bmatrix}^{\top}$
(26)

We use $\bm{N}_{0}$ to denote the multiplication of the last two matrices,
i.e.,

$\begin{bmatrix}\bm{O}&-\bm{J}&\bm{O}&\bm{O}&\bm{O}&\bm{O}\end{bmatrix}^{\top}$
(27)

Therefore, the Jacobian in (26) is replaced with $\bm{M}_{0}^{-1}\bm{N}_{0}$.

##  Appendix C The Derivation of the KKT Conditions of the Bilevel Program

Here, we derive the KKT conditions of the bilevel problem, where the IDSO
decision problem is at the upper-level, and the ITSO decision problem is at
the lower-level. We rewrite (10) by replacing the dual problem at the upper-
level with the primal problem in (8), i.e.,

$\displaystyle\mathop{\min}_{\bm{g},\bm{v},\bm{\theta}}\bm{c}^{\top}\bm{g}$
(28) $\displaystyle s.t.\eqref{8(b)},\eqref{8(c)},\eqref{8(d)},\eqref{8(e)}$
$\displaystyle\quad\
\mathop{\min}_{\bm{\xi},\bm{f}}\sum_{a\in\mathcal{A}^{e}}\gamma\cdot\xi_{a}\cdot\tau_{a}(\xi_{a})+\sum_{i\in\mathcal{C}}\lambda_{i}\cdot
d_{i}$ $\displaystyle\qquad\qquad\quad\
s.t.\eqref{6},\eqref{7(c)},\eqref{7(d)},\eqref{7(e)}$

We firstly derive the KKT conditions for the upper-level variables
$\bm{g},\bm{v},\bm{\theta}$ based on the Lagrangian function, i.e.,

$\displaystyle\sum_{j\in\Omega_{i}}K_{2ij}(\lambda_{i}-\lambda_{j}+\overline{\eta_{ij}}-\overline{\eta_{ji}}-\underline{\eta_{ij}}+\underline{\eta_{ji}})+$
$\displaystyle\qquad\lambda_{i}\frac{\partial
d_{i}}{\partial\theta_{i}}=0,\forall i\in\mathcal{N}$ (29a)
$\displaystyle\sum_{j\in\Omega_{i}}K_{1ij}(\lambda_{i}-\lambda_{j}+\overline{\eta_{ij}}-\overline{\eta_{ji}}-\underline{\eta_{ij}}+\underline{\eta_{ji}})+$
$\displaystyle\qquad\overline{\mu}_{i}-\underline{\mu}_{i}+\lambda_{i}\frac{\partial
d_{i}}{\partial v_{i}}=0,\forall i\in\mathcal{N}$ (29b) $\displaystyle
c_{i}-\lambda_{i}-\underline{\tau_{i}}+\overline{\tau_{i}}+\lambda_{i}\frac{\partial
d_{i}}{\partial g_{i}}=0,\forall i\in\mathcal{N}$ (29c) $\displaystyle
0\leq\bm{\underline{\tau}}\perp\bm{g}\geq 0$ (29d) $\displaystyle
0\leq\bm{\overline{\tau}}\perp(\bm{\overline{g}}-\bm{g})\geq 0$ (29e)
$\displaystyle
0\leq(K_{1ij}(v_{i}-v_{j})+K_{2ij}(\theta_{i}-\theta_{j})+f_{ij})\perp\underline{\eta_{ij}}\geq
0,$ $\displaystyle\qquad\forall i\in\mathcal{N},\forall j\in\Omega_{i}$ (29f)
$\displaystyle
0\leq(f_{ij}-K_{1ij}(v_{i}-v_{j})-K_{2ij}(\theta_{i}-\theta_{j}))\perp\overline{\eta_{ij}}\geq
0,$ $\displaystyle\qquad\forall i\in\mathcal{N},\forall j\in\Omega_{i}$ (29g)
$\displaystyle 0\leq(v_{i}-\underline{v})\perp\underline{\mu_{i}}\geq
0,\forall i\in\mathcal{N}$ (29h) $\displaystyle
0\leq(\overline{v}-v_{i})\perp\overline{\mu_{i}}\geq 0,\forall
i\in\mathcal{N}$ (29i) $\displaystyle\eqref{8(c)}$ (29j)

where (29a)-(29c) are the stationarity conditions (which include the gradient
of the upper-level variables $\bm{g},\bm{v},\bm{\theta}$ w.r.t. the lower-
level variables $\bm{d}$ [31]), and (29d)-(29i) are the complementarity
conditions, while (29j) ensures the primal feasibility.

As the lower-level problem is only associated with the dual variable
$\lambda_{i}$ and therefore $d_{i}$ is only a function of $\lambda_{i}$ (as we
show in (11)), the partial derivatives $\frac{\partial
d_{i}}{\partial\theta_{i}},\frac{\partial d_{i}}{\partial
v_{i}},\frac{\partial d_{i}}{\partial g_{i}}$ equal 0. Therefore, the
stationarity constraints in (29a),(29b),(29c) can be further simplified. With
the simplified stationarity constraints, the KKT conditions of the upper-level
problem are the same as the KKT conditions of IDSO problem in (8). Obviously,
the KKT conditions of the lower-level problem are the KKT conditions of ITSO
in (7). Therefore, the KKT conditions of the bilevel program (28) are in the
same form of (19), which ends the proof.