**4. Establish a Double-Layer Model**

In this section, a bilevel optimization model for the location of charging facilities is proposed. The upper level model determines the location of the charging facilities by selecting the top *p* links. The lower level model calculates users' generalized travel cost, and randomly allocates the flow demand between the OD pairs to the filtered paths through known charging facilities location.

## *4.1. Upper Level Problem*

The upper model is designed to maximize the total covered EV link flow by deploying a given number of charging facilities. The EV flow is covered when the charging facility is on the link. That is

$$\text{Max}\sum\_{a} v\_{ac} x\_a \tag{1}$$

$$\text{Subject to } \sum\_{a \in A} x\_a = p \tag{2}$$

Equation (1) is the objective function of the upper model, and Equation (2) is the budget constraint indicating the number of charging stations *p* in a given network.

### *4.2. Lower Level Problem*

The Bureau of Public Road (BPR) function in the lower model is employed [42].

$$t\_a(v\_{ag\_{"\prime}}, v\_{ac}) = t\_a^0 \left\{ 1 + 0.15 \times \left( \frac{v\_{ag} + v\_{ac}}{H\_a} \right)^4 \right\} \tag{3}$$

Since *U*, *t w ck*, *k* are parameters that are independent of the equilibrium flow of the link, the Jacobian matrix [43] in the lower layer problem is as follows:

$$\frac{\partial \overline{\epsilon}\_{k\emptyset}^{w}}{\partial v\_{ac}} = \frac{\partial \overline{\epsilon}\_{k\varepsilon}^{w}}{\partial v\_{a\emptyset}} = 0.6 \sum t\_a^0 \delta\_{ak}^w \frac{\left(v\_{a\emptyset} + v\_{a\emptyset}\right)^3}{H\_a 4} \tag{4}$$

It proves that the Jacobian matrix is symmetrical, and the lower layer model can be established as a convex function problem. It is assumed that all types of travelers make path selection in a random manner. According to the random utility theory, the probability that type-*i* traveler chooses the path *k* between OD pair *w* is:

$$p\_{ki}^{\
u} = \frac{\exp(-\theta\_i \overline{c}\_{ki}^{\overline{\mathbf{c}}^{\overline{\mathbf{c}}}})}{\sum\_{k \in \mathbb{R}\_w} \exp(-\theta\_i \overline{c}\_{ki}^{\overline{\mathbf{c}}^{\overline{\mathbf{c}}}})} \,\forall k, w, i \tag{5}$$

Assume that elastic travel demand *q <sup>w</sup> <sup>i</sup>* of type-*i* traveler is a strictly monotonically decreasing function of the expected minimum travel time between OD pair *w*, and with an upper bound.

$$q\_{\,\,i}^{\,w} = D\_{\,\,uni}(\,\mathbb{C}\_{\,\,uni}) \,\, \leq \overline{q}\_{i}^{w} \quad \forall w, i \tag{6}$$

According to the discrete selection theory [44,45], define *C wi* as:

$$\mathbb{C}\_{\text{wi}}(\boldsymbol{c}^{\text{w}}(\boldsymbol{\hat{x}}) = \mathbb{E}\left[\min\_{\mathbf{k} \in \mathbb{R}\_{\text{w}}} \left\{ \widetilde{\boldsymbol{c}}\_{\text{k}i}^{\text{w}} \right\} \boldsymbol{c}^{\text{w}}(\boldsymbol{\hat{x}}) \right] = -\frac{1}{\theta\_{i}} \ln \sum\_{\substack{\mathbf{k} \in \mathbb{R}\_{\text{w}} \\ \mathbf{k} \in \mathbb{R}\_{\text{w}}}} \exp(-\theta\_{i} \widetilde{\boldsymbol{c}}\_{\text{k}i}^{\text{w}}) \quad \forall \boldsymbol{w}, \boldsymbol{i} \tag{7}$$

For the type-*i* traveler, its SUE condition can be expressed as:

$$f\stackrel{w}{\;}\_{ki} = q\,\,^w\_i p\,^w\_{ki}\,\forall k, w, i \tag{8}$$

In the mixed network, SUE-ED problem can be described by the equivalent mathematical programming model:

$$\begin{split} \text{Min } Z(\mathbf{x}, f, q) &= \sum\_{\mathbf{z} \in \mathcal{A}} \int\_{0}^{v\_{\mathbf{z}}} t\_{\mathbf{z}}(\mathbf{w}) \text{d}w + \sum\_{i \in I} \frac{1}{\mathcal{O}\_{i}} \sum\_{w \in \mathcal{W}} \sum\_{k \in \mathbb{R}\_{w}} f^{w}\_{k}(\text{ln}\, f^{w}\_{k\mathbf{i}} - 1) \\ &- \sum\_{i \in I} \frac{1}{\mathcal{O}\_{i}} \sum\_{w \in \mathcal{W}} \sum\_{k \in \mathbb{R}\_{w}} q^{w}\_{i}(\text{ln}\, q^{w}\_{i} - 1) - \sum\_{i \in I} \sum\_{w \in \mathcal{W}} \int\_{0}^{q^{\mathbf{z}}} D^{-1}\_{\mathbf{w}i}(\mathbf{w}) \text{d}w \end{split} \tag{9}$$

$$\text{Subject to: } q\,\,^w\_{\,\,l} = \sum\_{k} f\,\,^w\_{\,\,ki'}\,\,\forall w\,\,i\tag{10}$$

$$\left|f\right|\_{\text{ki}}^w \ge 0, \ \forall k, w, i \tag{11}$$

$$\{q\:\:\:q\:\:^w\_i\geq 0,\;\:\forall w\,, i\tag{12}$$

$$\upsilon\_{\rm ai} = \sum\_{w \in \mathcal{W}} \sum\_{k \in R\_w} f\_{ki}^{\rm w} \delta\_{ak'}^{w} \,\,\forall a \,\, i \tag{13}$$

$$R\_{\mathbf{c}} - l\_k^{i\bar{j}} \ge 0, \ \forall k \tag{14}$$

Equation (10) is the flow conservation constraint. Equation (11) is the non-negative constraint of path flow for type-*i* travelers. Equation (12) is the non-negative constraint of type-*i* traveler's OD flow demand. Equation (13) is the correlation between link flow and path flow. The novelty of this problem is that the introduction of sub-paths in the Equation (14) can generate a feasible set of paths in advance from the finite paths between each OD pair. The superscripts *i* and *j* include the origin *r* and the destination *s* of all OD pairs in Equation (14).

It is necessary to prove the equivalence between the solution of the proposed program Equation (9) and the solution of the SUE-ED model. The generalized Lagrangian function [46] of the mathematical model is constructed as follows:

$$\mathcal{L} = \mathcal{Z} + \sum\_{i} \sum\_{w} \lambda\_{\bar{w}i} (q\_i^w - \sum\_{k \in \mathcal{R}\_w} f\_{\bar{k}i}^w) - \sum\_{i} \sum\_{w} \sum\_{k} \mu\_{\bar{k}i}^w f\_{\bar{k}i}^w - \sum\_{i \in \mathcal{I}} \sum\_{w \in \mathcal{W}} \nu\_{\bar{w}i} q\_i^w - \sum\_{w} \sum\_{k} \alpha\_{\bar{k}i}^{ij} \star (R\_{\bar{\epsilon}} - l\_{\bar{k}}^{ij}) \tag{15}$$

λ*wi*, μ*<sup>w</sup> ki*, <sup>ν</sup>*wi*, <sup>α</sup>*ij ke* are the Lagrange multipliers.

$$\begin{array}{ll}\text{If } l\_{\stackrel{kl}{k}}^{lj} \le R\_{\mathfrak{c}} & \alpha\_{\stackrel{kl}{k\stackrel{\nu}{\cdot}}}^{lj} = 0\\\text{If } l\_{\stackrel{kl}{k}}^{lj} \ge R\_{\mathfrak{c}} & \alpha\_{\stackrel{kl}{k\stackrel{\nu}{\cdot}}}^{lj} \ge 0 \end{array} \tag{16}$$

According to the Kuhn Tucker conditions [47], Equation (15) must satisfy the following conditions at the extreme point:

$$\frac{\partial L}{\partial f \stackrel{w}{k^w}} = 0, \; \mu^w\_{ki} f \stackrel{w}{k^i} = 0, \; \mu^w\_{ki} \ge 0 \quad \forall k, w, i \tag{17}$$

$$\frac{\partial L}{\partial q\_{\stackrel{w}{i}}^{\;\;w}} = 0, \; \nu\_{\text{w}i} q\_{\stackrel{w}{i}}^{\;\;w} = 0, \; \nu\_{\text{w}i} \ge 0 \quad \forall w, i \tag{18}$$

The partial derivative of *f <sup>w</sup> ki* for Equation (15) is derived as follows:

$$\frac{\partial L}{\partial f \frac{w}{ki}} = \sum\_{a} t\_a (\upsilon\_a) \delta\_{ak}^w + \frac{1}{\Theta\_i} \ln f \, \, \frac{w}{ki} - \lambda\_{zi} - \mu\_{ki}^w = 0 \quad \forall k, w, i \tag{19}$$

Note that if *f <sup>w</sup> ki* <sup>=</sup> 0, then <sup>∂</sup>*<sup>L</sup>* ∂ *f <sup>w</sup> ki* does not exist. The above formula is only valid when *f <sup>w</sup> ki* <sup>&</sup>gt; 0. So <sup>μ</sup>*<sup>w</sup> ki* = 0.

$$f^{\
u}\_{\ k\dot{i}} = \exp\left[-\theta\_i(\tilde{\mathbf{c}}^{\mu}\_{\ k\dot{i}} - \lambda\_{\text{wi}})\right] = \exp\left[-\theta\_i \tilde{\mathbf{c}}^{\mu\nu}\_{\ k\dot{i}}\right] \cdot \exp\left[-\theta\_i \lambda\_{\text{wi}}\right] \tag{20}$$

$$q\ \_i^w = \exp\left(-\theta\_i \lambda\_{\text{uv}}\right) \sum\_{k \in R\_{\text{uv}}} \exp\left(-\theta\_i \tilde{c}\_{ki}^w\right) \tag{21}$$

$$p\_{\
i
i}^{\
w} = \frac{\exp(-\theta\_i \overline{c}\_{ki}^{\rm w})}{\sum\_{k \in \, \mathbb{R}\_w} \exp(-\theta\_i \overline{c}\_{ki}^{\rm w})} \,\forall k, w, i \tag{22}$$

Equation (22) shows that type-*i* traveler follows the Logit model to select the travel path, which satisfies the SUE-ED condition described in Equation (8).

Calculate the partial derivative of *q <sup>w</sup> <sup>i</sup>* for Equation (15) as follows:

$$\frac{\partial L}{\partial q\_{\,\,i}^{\,w}} = -\frac{1}{\partial\_i} \ln q\_{\,\,i}^{\,w} - D\_{\rm wi}^{-1}(q\_{\,\,i}^{\,w}) + \lambda\_{\rm wi} \because \nu\_{\rm wi} = 0 \; \forall w, i \tag{23}$$

Because ln *q <sup>w</sup> <sup>i</sup>* should exist, so ν*wi* = 0. Then derive from Equation (21):

$$\theta\_i \lambda\_{wi} = \ln q\_i^w - \ln \sum\_{r \in R\_{\vartheta}} \exp \left[ -\theta\_i \tilde{c}\_{ri}^w \right] \tag{24}$$

Comparing Equation (23) and Equation (24), the inverse function of the flow demand function is as follows:

$$D\_{wi}^{-1}(q\,\,^{w}\_{i}) = -\frac{1}{\theta\_i} \ln \sum\_{r \in R\_{w}} \exp\left[-\theta\_i \tilde{c}\_{ri}^{w}\right] \forall w, i \tag{25}$$

Then *q <sup>w</sup> <sup>i</sup>* = *Dwi*(*Cwi*). This shows that Equations (9)–(14) can be used to represent a multi-user SUE problem under elastic demand. The BPR function *ta(va*) is a strictly monotonically increasing function of the link flow *va*. The objective function is a strict convex function about the link flow vector *v* and the path flow vector *f*. At the same time, the constraints of Equations (9)–(14) are linear equality constraints and non-negative constraints, so its solution space is a convex set. According to the optimization theory, the strict convex function defined on the convex set only has one optimal solution.
