Comparison of stability regions for a line distribution network with stochastic load demands

Abstract

We compare stability regions for different power flow models in the process of charging electric vehicles (EVs) by considering their random arrivals, their stochastic demand for energy at charging stations, and the characteristics of the electricity distribution network. We assume the distribution network is a line with charging stations located on it. We consider the Distflow and the Linearized Distflow models, and we assume that EVs have an exponential charging requirement, that voltage drops on the distribution network stay under control, and that the number of charging stations N goes to infinity. We investigate the stability of utility-optimizing power allocations in large distribution networks for both power flow models by controlling the arrival rate of EVs to charging stations. For both power flow models, we show that, to obtain stability, the maximum feasible arrival rate, i.e., stability region of vehicles, is decaying as $1/N^2$, and the difference between those arrival rates is up to constants, which we compare explicitly.

1 Introduction

The popularity of electric vehicles (EVs) has been growing due to their increased range, lower cost of batteries, and governmental subsidies. However, EVs need to be charged and the current infrastructure cannot support their increasing demand for energy. This can therefore cause capacity problems in electricity networks in the (very near) future [14]. Modifying the existing infrastructure is costly and constrained by the limits of the electricity network. Therefore, we need to explore the impact of the increasing charging demand on the quality of service to EV drivers, subject to network constraints. More specifically, in this paper, we characterize the set of feasible arrival rates of EVs to charging stations, guaranteeing the quality of service and satisfying network constraints.

It has been shown that for the penetration of EVs to 20–40% of the households, existing electricity networks should be able to deal with an increasing demand for energy from EVs if the majority is restricted to low charging rates at off-peak times [1]. However, controlling the rate at which EVs charge can lead to better utilization of existing networks. For example, in [25], an optimal charging rate for each EV is determined. The goal is to maximize the total power that can be delivered to the vehicles while operating within network limits. In [13], individual preferences regarding their charging rates are also taken into account. Those users who want a higher rate can receive it, but must pay a higher price.

Irrespective of the allocation of charging capacity to EVs, it is important to respect network limits. In an electricity network, an important constraint is the requirement of keeping voltage losses, or in other words, the voltage drop, on a cable in the network under control. This voltage drop comes from the transportation of electricity due to physical properties of cables in the network. Keeping the voltage losses under control ensures that every customer receives safe and reliable energy at a voltage that is within some standard range, which varies from one country to another [16]. For example, according to the Dutch law, the voltage loss in distribution networks is not allowed to be more than 4.5% [30].

Motivated by this, we consider EVs charging in a neighborhood such that voltage drops on the distribution network stay under control. EVs arrive randomly (at charging stations) to get charged and have different stochastic demands for energy. Moreover, the charging rates allocated to EVs depend on the number of cars charging simultaneously in the neighborhood. We model this process as a queue, with EVs representing jobs, charging stations classified as servers, and the service being delivered as the power supplied to EVs, constrained by the physical limitations of the distribution network. These physical limitations are modeled by different power flow models.

The goal of this paper was to compare different power flow models while guaranteeing stability of the queuing model under distribution network constraints. Whenever we write stability, we mean stability of the queuing model, unless stated otherwise. In our setting, we assume that cars arrive at the same rate at each charging station and we want to find the maximum feasible arrival rate in terms of the number of charging stations N such that voltage drops are within constraints. This is very challenging due to the uncertainty of the stochastic electricity demand of cars. Another difficult aspect is the modeling of the power flow dynamics in the distribution network. In this paper, we approximate the power flow dynamics in the distribution network by the commonly used Distflow and the Linearized Distflow models, because they are good approximations and of their appealing mathematical properties. Moreover, studying these power flow models provides insight into what happens by linearizing physical models. The main stability results (Theorems 3.1, 3.2) characterize the maximal feasible arrival rates under distribution network constraints for different power flow models. More specifically, they show, on the one hand, that for large distribution networks (to obtain stability) the maximal feasible arrival rates under both power flow models decay at the same rate, and on the other hand that the difference between these arrival rates is up to constants, which we compare explicitly.

1.1 Distribution networks

An electric grid is a connected network that transfers electricity from producers to consumers. It consists of generating stations that produce electric power, high voltage transmission lines that carry power from distant sources to demand centers, and distribution lines that connect individual customers, e.g., houses, charging stations, etc. We focus on a network that connects a generator to individual customers with only distribution lines. Such a network is called a distribution network. Such a network is not necessarily limited to EV drivers that want to charge their cars. All individual consumers of energy with the same characteristics as EVs, such as arrival and demand patters, fit this framework.

In the rest of the paper, we assume that the distribution network, consisting of one generator, several charging stations, and distribution lines, has a line topology. The generator that produces electricity is called the root node. Charging stations consume power and are called the load nodes. Thus, we represent the distribution network by a graph (here, a line) with a root node, load nodes, and edges representing the distribution lines.

1.2 Power flow models

In order to model the power flow in the network, we use approximations of the alternating current (AC) power flow equations [24]. These power flow equations characterize the steady-state relationship between power injections at each node, the voltage magnitudes, and phase angles that are necessary to transmit power from generators to load nodes. First, we study a load flow model known as the branch flow model or the Distflow model [4], and second, a linearized version of the Distflow model called the Linearized Distflow model. Both power flow models focus on quantities such as complex current and power flowing on the distribution lines. Due to the specific choice for the distribution network as a line, both power flow models have a recursive structure, which we can exploit.

The accuracy and effectiveness of the Linearized Distflow model has been numerically justified in the literature [4]. Its use is justified by the fact that the nonlinear terms in the equations of the Distflow model represent the losses that appear if electric power is transferred over the lines. These losses should be, in practice, much smaller than the active and reactive power terms that show up in the equations. In this paper, we show analytically that the impact of neglecting these losses, in terms of stability, is insignificant if the network is large.

1.3 Stability conditions for EV charging

Customers can charge their EVs at charging stations. In practice, EVs are served simultaneously, because they require concurrent usage of all upstream distribution lines between the location of an EV and the generator of the distribution network. This feature of sharing service capacity under concurrent users can be captured adequately in queuing theory by a resource-sharing network and specifically a bandwidth-sharing network. Bandwidth-sharing networks have been successfully applied in many areas, for example, in computer systems and communication networks [21, 22, 34].

There is vast literature on stability conditions in bandwidth-sharing networks. The question of stability was first posed under Markovian assumptions in [10, 32], and cover the stability of such networks for weighted proportionally fair policies. In [5], the latter results are generalized to a more general family of policies called weighted $\alpha$-fair policies [6]. Stability results for more general arrival and service distributions have also been derived [8, 12, 23]. In particular, [20] has shown the stability of the proportionally fair policy under phase-type service distributions. In a continuous-time Markovian setting, [26, 27] present stability conditions for a more general class of utility-maximizing allocations, which include weighted $\alpha$-fair allocations. For a more detailed overview, see [33] or a recent survey paper by Williams [34].

The literature on stability results for EV charging is limited to numerical experiments. An early paper on stability analysis in EV charging is Huang et al. [15]. The authors present a new quasi-Monte Carlo stability analysis method to assess the dynamic effects of plug-in electric vehicles in power systems. They conclude that improvements for stability control are worth further study since the number of EVs is growing.

Simulation studies are conducted to obtain stability conditions in Carvalho et al. [7]. Here, the authors find that there is a threshold value $\lambda _c$ on the arrival rates, such that if the actual arrival rate $\lambda$ is greater than this threshold, i.e., $\lambda >\lambda _c$, some vehicles have to wait for increasingly long times to fully charge. Similar findings are obtained in [9, 11, 29, 35]. Up to a certain point, the distribution network is able to serve all EVs properly, but if the number of EVs in the system is too high, it cannot be guaranteed, e.g., that voltage drops stay under control, that distribution transformers, distribution lines, and cables are not overloaded, that power losses are not too big, or that peak demands are not too high.

The present paper builds upon Aveklouris et al. [2]. In this study, the authors consider a distribution network used to charge EVs such that voltage drops stay under control, taking into account randomness of future arriving EVs and power demands. The work focuses on a fluid approximation for the number of uncharged EVs, while we focus on explicit conditions to ensure stability of the queuing model under both power flow models. We present maximal arrival rates to obtain stability for both power flow models. We use differential and integral calculus to compute these maximal feasible arrival rates under the Distflow model and to study the difference between maximal feasible arrival rates obtained under the Linearized Distflow and the Distflow models. The main conclusion that we draw is that these rates are the same as the size of the distribution network increases, up to constants. The existence of conditions guaranteeing stability has been shown in [26, 27], where the authors establish stability for networks of interacting queues governed by utility-maximizing service rate allocations using direct applications of Lyapunov–Foster-type criteria.

The structure of the paper is as follows. In Sect. 2, we provide a detailed model description. In particular, we introduce the queuing model, the distribution network model, and the power flow models. In Sect. 3, we find stability regions for both power flow models. The stability results are presented in Sects. 3.1 and 3.2. The comparison between the stability regions of both power flow models is made in Sect. 3.3. The aforementioned stability results under the Linearized Distflow and the Distflow models are proven in Sects. 4 and 5, respectively. From an engineering point of view, the computation of the stability region for the Distflow model, as in Sect. 3.2, can also be done numerically, via an iterative approach. Therefore, we show how to use Newton’s method to compute the maximal feasible arrival rate under the Distflow model in Sect. 6. In Sect. 7, we state most important outcomes of this paper. The appendix focuses on proofs and a theorem that are used in Sects. 2.3.1, 4, and 5. A more detailed description is found in Appendix.

In the remainder of this section, we list the notation we use throughout the paper. An overview of notations is also given in Appendix F.

1.4 Notation

All vectors and matrices are denoted by bold letters. Vector inequalities hold coordinate-wise; namely, ${\textbf{x}}>{\textbf{y}}$ implies that $x_i>y_i$ for all i. The imaginary unit is denoted by $\textrm{i}$, and the absolute value of a complex number $z=x+\textrm{i} y$ is $\left| z\right| =\sqrt{x^2+y^2}$. Also, the complex conjugate $x-\textrm{i}y$ of z is denoted by ${\overline{z}}$.

The following operations are defined on $x,y\in {\mathbb {R}}$:

$$\begin{aligned} \lfloor x\rfloor&:=\max \{m\in {\mathbb {Z}}:m\le x\}, \\ \lceil x\rceil&:=\min \{n\in {\mathbb {Z}}:n\ge x\}. \end{aligned}$$

Denote the space of functions $f:[0,T]\rightarrow {\mathbb {R}}$ that are right continuous with left limits, i.e., càdlàg functions, by

$$\begin{aligned} {\textbf{D}}[0,T]. \end{aligned}$$

Furthermore, we define the space ${\textbf{D}}_{\ge 1}[0,T]$ as

$$\begin{aligned} {\textbf{D}}_{\ge 1}[0,T] := \{f\in {\textbf{D}}[0,T]: \inf _{t\in [0,T]} f(t)\ge 1.\}. \end{aligned}$$

For a function $f(\cdot )$ defined on (a subset of) ${\mathbb {R}}$, $f'(t)$ denotes its derivative at t and $f''(t)$ its second derivative at t, when these exist.

2 Model description

This section describes the three main components of the EV-charging model. In Sect. 2.1, we describe the characteristics of the queuing model, i.e., the evolution of the number of cars charging at each charging station. In Sect. 2.2, we specify the distribution network model, and in Sect. 2.3, we introduce the Distflow and the Linearized Distflow power flow models.

2.1 Queuing model of EV charging

We use a queuing model to study the process of charging EVs in a distribution network. Recall that in this setting, EVs referred to as jobs require service. This service is delivered by charging stations, referred to as servers, and the service being delivered is the power supplied to EVs. Since we are interested in the maximal feasible arrival rate, we assume that we can simultaneously charge an infinite number of EVs so that we are not limited by space, i.e., a finite number of charging stations or limited numbers of parking lots. At the same time, this is not out of the physical realm if we consider, for example, wireless charging.

Thus, in the queuing system, we consider N single-server queues, each having its own arrival stream of jobs. Denote by ${\textbf{X}}(t) = (X_1(t),\ldots ,X_N(t))$ the vector giving the number of jobs at each queue at time t. We make the following assumption on the arrival rates and service requirements of all EVs.

Assumption 2.1

At all charging stations, all EVs arrive independently according to Poisson’s processes with the same rate $\lambda$ and have independent service requirements which are Exp(1) random variables.

We make this assumption, since our goal is to compare different power flow models. Therefore, we exclude other sources of variability from our model, such as different arrival and demand patterns.

At each queue, all jobs are served simultaneously and start service immediately (there is no queuing). Furthermore, each job receives an equal fraction of the service capacity allocated to a queue. Denote by ${\tilde{\textbf{p}}}(t) = ({\tilde{p}}_1(t),\ldots ,{\tilde{p}}_N(t))$ the vector of service capacities allocated to each queue at time t. This represents the active power that is allocated to each node and is dependent on the number of EVs that are charging at each charging station. This means that at each queue j, it takes 1/${\tilde{p}}_j(t)$ time units to serve one job. However, service capacities are state-dependent and subject to changes, and the dynamics are more complicated. See below for details.

We can then represent the number of electric vehicles charging at every station as an N-dimensional continuous-time Markov process. The evolution of the queue at node j is given by

$$\begin{aligned} X_j(t) \rightarrow X_j(t)+1~\text {at rate}~\lambda \end{aligned}$$

and

$$\begin{aligned} X_j(t) \rightarrow X_j(t)-1~\text {at rate}~{\tilde{p}}_j(t). \end{aligned}$$

From now on, for simplicity, we drop the dependence on time t from the notation. For example, we write $X_j$ and ${\tilde{p}}_j$ instead of $X_j(t)$ and ${\tilde{p}}_j(t)$. We assume that the rates ${\tilde{\textbf{p}}}$ may be allocated according to the current vector ${\textbf{X}}=(X_1,\ldots ,X_N)$ of number of jobs.

As discussed in introduction, we use the existence of stability conditions for a general class of utility-maximizing allocations, which include the popular and well-known $\alpha$-fair algorithms [6]. Therefore, we discuss the use of an $\alpha$-fair algorithm to allocate the charging capacity to EVs in more detail.

In state ${\textbf{X}}$, an $\alpha$-fair algorithm allocates ${\tilde{p}}_j/X_j$ to each EV at charging station j, with ${\tilde{\textbf{p}}}=({\tilde{p}}_1,\ldots ,{\tilde{p}}_N)$ the solution of the utility optimization problem

$$\begin{aligned} {\tilde{\textbf{p}}} \in \,{\mathrm{arg~max}}\, \sum _{j=1}^N X_j U_j^{(\alpha )}\left( \frac{{\tilde{p}}_j}{X_j}\right) , \end{aligned}$$

(2.1)

subject to physical constraints on the vector ${\tilde{\textbf{p}}}$ of allocated power and where

$$\begin{aligned} U_j^{(\alpha )}(x_i) = {\left\{ \begin{array}{ll} \log x_i &{}\ \text {if}\ \alpha =1, \\ x_i^{1-\alpha }/(1-\alpha ) &{} \ \text {if}\ \alpha \in (0,\infty )\backslash \{1\}, \end{array}\right. },\ x_i\ge 0. \end{aligned}$$

A fair allocation is one that does not penalize some users arbitrarily. In our setting, that could be an allocation where the charging rate allocated to an EV does not depend on the specific charging station it is connected to.

The parameter $\alpha$ measures the degree of fairness of the allocation. Popular choices are $\alpha \rightarrow 0$, where the total allocated power tends to be maximized, but the allocation is very unfair. The choice $\alpha =1$, where we end up with proportional fairness, that tends to maximize the utility of the total power allocation or $\alpha =2$, that corresponds to the minimum potential delay allocation, where the total charging time tends to be minimized. Last, we have the limiting case $\alpha \rightarrow \infty$, where the minimum of the power allocated to any charging station tends to be maximized, namely max–min fairness [21].

2.2 Distribution network model

The distribution network is modeled as a directed graph ${\mathcal {G}}=({\mathcal {I}},{\mathcal {E}})$, where we denote by ${\mathcal {I}} = \{0,\ldots ,N\}$ the set of nodes and by ${\mathcal {E}}$ its set of directed edges, assuming that node 0 is the root node. We assume that ${\mathcal {G}}$ has a line topology. Each edge $\epsilon _{j-1,j}\in {\mathcal {E}}$ represents a line connecting nodes $j-1$ and j, where node j is further away from the root node than node $j-1$. Each edge $\epsilon _{j-1,j}\in {\mathcal {E}}$ is characterized by the impedance $z=r+\textrm{i} x$, where $r,x\ge 0$ denote the resistance and reactance along the lines, respectively. We make the following assumption. This is a natural assumption given that $r>>x$ in distribution networks [17, 28].

Assumption 2.2

All edges have the same resistance value $r>0$ and reactance value $x=0$.

Furthermore, let ${\tilde{s}}_j = {\tilde{p}}_j + \textrm{i} {\tilde{q}}_j$ be the complex power consumption at node j. Here, ${\tilde{p}}_j$ and ${\tilde{q}}_j$ denote the active and reactive power consumption at node j, respectively (cf. (2.1)). By convention, a positive active (reactive) power term corresponds to consuming active (reactive) power. Since EVs can only consume active power [7], it is natural to make the following assumption.

Assumption 2.3

The active power ${\tilde{p}}_j$ is non-negative and the reactive power ${\tilde{q}}_j$ is zero at all charging stations $j\in {\mathcal {I}}$.

Let ${\tilde{V}}_j$ denote the voltage at node j. Given Assumptions 2.2 and 2.3, observe [e.g., from Equations (2.6)–(2.9)] that the voltages at each node j, ${\tilde{V}}_j$, can be chosen to have zero imaginary components [2, 7]. For each $\epsilon _{j-1,j}\in {\mathcal {E}}$, let $I_{j-1,j}$ be the complex current and ${\tilde{S}}_{j-1,j}={\tilde{P}}_{j-1,j}+\textrm{i}{\tilde{Q}}_{j-1,j}$ be the complex power flowing from node $j-1$ to j. Here, ${\tilde{P}}_{j-1,j}$ and ${\tilde{Q}}_{j-1,j}$ denote the active and reactive power flowing from node $j-1$ to j. The model is illustrated in Fig. 1.

Fig. 1

Line network with N charging stations and arriving vehicles at rate $\lambda$

Voltage drop constraint The distribution network constraints, that is in our case only the voltage drop constraint, are described by a set ${\mathcal {C}}$. The set ${\mathcal {C}}$ is contained in an N-dimensional vector space and represents feasible power allocations. In our setting, a power allocation is feasible if the maximal voltage drop; i.e., the relative difference between the base voltage ${\tilde{V}}_0$ and the minimal voltage in all buses between the root node and any other node is bounded by some value $\Delta \in (0,\frac{1}{2}]$. Thus, the distribution network constraints can be described as

$$\begin{aligned} {\mathcal {C}}:= \left\{ {\textbf{p}}:\ \frac{{\tilde{V}}_0-\min _{1\le j\le N} {\tilde{V}}_j}{{\tilde{V}}_0}\le \Delta ,\quad 0<\Delta \le \frac{1}{2} \right\} . \end{aligned}$$

(2.2)

In Sect. 2.3.1, we show there is a technical reason why we use the explicit interval $(0,\frac{1}{2}]$ for $\Delta$. However, this is not out of the physical realm, since in practice, we want the maximal voltage drop to be no more than a small percentage of the base voltage, e.g., $\Delta = 0.05$ or $\Delta =0.1$. In Sect. 2.3, we give more concrete definitions of the constraint set ${\mathcal {C}}$.

2.3 Power flow models

We introduce two commonly used models to represent the power flow that are valid for radial systems, i.e., systems where all charging stations have only one (and the same) source of supply. They are called the Distflow and Linearized Distflow model [4]. Both models are valid when the underlying network topology is a tree, which is indeed the case in this paper (as we consider a line topology). Moreover, we show that for a line topology, the power flow model equations can be rewritten recursively for the voltages. For an overview of other representations of power flow, we refer the reader to Molzahn and Hiskens [24, Chapter 2].

Given the impedance r, the voltage at the root node ${\tilde{V}}_0$, and the power consumptions ${\tilde{p}}_j, j=1,\ldots ,N$, both power flow models satisfy three relations. First, we have power balance at each node: For all $j\in {\mathcal {I}}\backslash \{0\}$,

$$\begin{aligned} {\tilde{S}}_{j-1,j}-r\left| I_{j-1,j}\right| ^2 = {\tilde{s}}_j + {\tilde{S}}_{j,j+1}. \end{aligned}$$

(2.3)

Here, on the one hand, the quantity $r|I_{j-1,j}|^2$ represents line loss so that ${\tilde{S}}_{j-1,j}-r|I_{j-1,j}|^2$ is the receiving-end complex power at node j from node $j-1$. On the other hand, the delivering-end complex power is the sum of the consumed power at node j and the complex power flowing from node j to node $j+1$. Second, by Ohm’s law, we have for each edge $\epsilon _{j-1,j} \in {\mathcal {E}}$,

$$\begin{aligned} {\tilde{V}}_{j-1} - {\tilde{V}}_j = rI_{j-1,j} \end{aligned}$$

(2.4)

and third, due to the definition of complex power, we have for each edge $\epsilon _{j-1,j} \in {\mathcal {E}}$,

$$\begin{aligned} {\tilde{S}}_{j-1,j} = {\tilde{V}}_{j-1}{\overline{I}}_{j-1,j}. \end{aligned}$$

(2.5)

Decomposing Eq. (2.3) in real and imaginary parts leads to Eqs. (2.6) and (2.7). Isolating the voltage ${\tilde{V}}_j$ from (2.4) and taking the squared magnitude on both sides, lead to (2.8). Similarly, taking the squared magnitude on both sides of (2.5) leads to (2.9). Thus, we have that an equivalent formulation of (2.3)–(2.5) is given by:

$$\begin{aligned}&{\tilde{p}}_j = {\tilde{P}}_{j-1,j}-r\left| I_{j-1,j}\right| ^2-{\tilde{P}}_{j,j+1},&\quad j\in \{1,\ldots ,N\}, \end{aligned}$$

(2.6)

$$\begin{aligned}&0 = {\tilde{Q}}_{j-1,j}-{\tilde{Q}}_{j,j+1},&\quad j\in \{1,\ldots ,N\}, \end{aligned}$$

(2.7)

$$\begin{aligned}&{\tilde{V}}_j^2 = {\tilde{V}}_{j-1}^2-2r{\tilde{P}}_{j-1,j}+\left| r\right| ^2\left| I_{j-1,j}\right| ^2,&\quad \epsilon _{j-1,j}\in {\mathcal {E}}, \end{aligned}$$

(2.8)

$$\begin{aligned}&\left| I_{j-1,j}\right| ^2 = \frac{\left| {\tilde{S}}_{j-1,j} \right| ^2}{{\tilde{V}}_{j-1}^2},&\quad \epsilon _{j-1,j}\in {\mathcal {E}}, \end{aligned}$$

(2.9)

with the understanding that when j is a leaf node ${\tilde{P}}_{j,j+1}=0$ and ${\tilde{Q}}_{j,j+1}=0$ in (2.6) and (2.7), respectively. Equations (2.6)–(2.9) are the Distflow model equations. Given Assumptions 2.2 and 2.3, these equations yield a first-order recurrence relation in two real variables: the active power flow between two subsequent buses and the square of the voltages at each node.

In what follows, in the case of the Distflow model, we rewrite Equations (2.6)–(2.9) to obtain a second-order recurrence relation in a single variable: the voltage at each node. In the case of the Linearized Distflow model, we make an additional assumption to simplify Equations (2.6)–(2.9) first and then find recursive equations for the voltages.

Using these recursive expressions for the voltages under both power flow models, we determine whether a given power allocation $\tilde{{\textbf{p}}}$ is feasible or not under both power flow models (given properties of the distribution network, such as the impedance r and base voltage ${\tilde{V}}_0$). In other words, in Sects. 2.3.1 and 2.3.2, we specify the constraint set ${\mathcal {C}}$ in (2.2) for both power flow models.

Furthermore, to avoid a notational issue, we denote the expressions for the different models using the superscripts D (Distflow model equations) and L (Linearized Distflow equations). If there is no superscript specified, the expressions hold for both models.

2.3.1 Distflow

The first power flow model we discuss is called the Distflow model. This model was first proposed in [3] and [4] and derived under the assumption that the underlying network topology is a tree. Further simplifying the network topology into a line and assuming that load nodes consume active power only, direct manipulation and rewriting of Eq. (2.8) with the help of Eqs. (2.6), (2.7), and (2.9) yield a second-order recursive relation between the voltages,

$$\begin{aligned} {\tilde{V}}_{N-1}^D&= {\tilde{V}}_N^D + \frac{r{\tilde{p}}_N}{V_N^D}, \end{aligned}$$

(2.10)

$$\begin{aligned} {\tilde{V}}_{j-1}^D&= 2{\tilde{V}}_j^D-{\tilde{V}}_{j+1}^D+\frac{r{\tilde{p}}_j}{{\tilde{V}}_j^D}, \quad j = 1,\ldots ,N-1 \end{aligned}$$

(2.11)

where the voltage at the root node, ${\tilde{V}}_0^D$, and the charging rates ${\tilde{p}}_j$ for $j=1,\ldots ,N$ are given. Relabeling Eqs. (2.10) and (2.11) by ${V_j}^D:= {\tilde{V}}_{N-j}^D$ and $p_j:={\tilde{p}}_{N-j}$ results in a new recursive relation such that each further term of the sequence is defined as a function of the preceding terms. Notice that in terms of the new variables ${V_N}^D$ and for $j=0,\ldots ,N-1$, $p_j$ are given. Thus, we have the following second-order recursive relation between the voltages,

$$\begin{aligned}&{V_{1}}^D = {V_0}^D + \frac{rp_0}{{V_0}^D}, \end{aligned}$$

(2.12)

$$\begin{aligned}&{V_{j+1}}^D = 2{V_j}^D-{V_{j-1}}^D+\frac{rp_j}{{V_j}^D}, \quad j = 1,\ldots ,N-1 \end{aligned}$$

(2.13)

where ${V_N}^D$ is known.

Given the recurrence relation in (2.12) and (2.13), it is not straightforward to find out if a given power allocation ${\tilde{\textbf{p}}}$ is feasible or not. Computing if a given power allocation is feasible, requires finding a voltage ${V_0}^D$ such that the known voltage ${V_N}^D$ satisfies the maximal voltage drop constraint as in Eq. (2.2).

However, in what follows, we find an equivalence of the recursive relations in (2.12) and (2.13) between the situation where the voltage ${V_N}^D$ is known and the situation where the voltage ${V_0}^D$ is known, such that the voltage drop constraint is satisfied. In the latter situation, one can find out easily if a given power allocation ${\textbf{p}}$ (with distribution network parameter r and voltage ${V_N}^D$) is feasible: Initialize the recursion with ${V_0}^D>0$ and iterate through (2.12) and (2.13) and check, using the voltage ${V_N}^D$, if the maximal voltage drop constraint is satisfied. To this end, we let the voltage ${V_0}^D$ be positive and define $v_j: {\mathbb {R}}_+ \rightarrow {\mathbb {R}}_+$ be the function defined by the equation $v_j({V_0}^D) = {V_j}^D$. The equivalence is found in Proposition 2.1.

Proposition 2.1

Let the resistance r and power consumptions $p_j, j=0,\ldots ,N-1$ be given. Then, the following equivalence, for $0< \Delta \le \frac{1}{2}$, for the voltages (under the Distflow model) at the root node and the node at the end of the line holds:

$$\begin{aligned}{} & {} \Big (\exists \,x \ge (1-\Delta )c~\mathrm{s.t.}~{V_0}^D = x~\text {and}~v_N({V_0}^D) = c\Big )~\mathrm{if\, and \,only\, if}\nonumber \\{} & {} \Big ({V_0}^D = (1-\Delta )c~\textrm{and}~v_N({V_0}^D) \le c\Big ), \end{aligned}$$

(2.14)

where $c>0$.

To prove the desired equivalence, we follow the approach in [31, Chapter 5]. We first derive an alternative expression for the voltages under the Distflow model as a double summation, instead of the voltages given by Eqs. (2.12) and (2.13). Then, we show that the voltage at node j, for $j=0,\ldots ,N$, as a function of the voltage ${V_0}^D$ is an increasing function. Here, we rely on the fact that the sequence of voltages $V_0^D,\ldots ,V_N^D$ is increasing. This is proven in Lemma 2.3. The proofs of Proposition 2.1 and Lemmas 2.1–2.3 are found in Appendix A.

In Lemma 2.1, we derive an alternative expression for the voltages under the Distflow model.

Lemma 2.1

Let ${V_j}^D, j=0,1,\ldots ,N$, be as in (2.12) and (2.13). Then,

$$\begin{aligned}&{V_{j+1}}^D-{V_j}^D=\sum _{i=0}^j \frac{rp_i}{{V_i}^D},\quad j = 0,\ldots ,N-1, \end{aligned}$$

(2.15)

$$\begin{aligned}&{V_{j}}^D = {V_0}^D + \sum _{n=0}^{j-1} \sum _{i=0}^n \frac{rp_i}{{V_i}^D},\quad j = 0,\ldots ,N. \end{aligned}$$

(2.16)

Now that we have an alternative expression for the voltages under the Distflow model, we show, in Lemma 2.2 that the voltages under the Distflow as a function of the voltage ${V_0}^D$ are increasing functions.

Lemma 2.2

Let the voltage ${V_0}^D$ be positive $({V_0}^D>0)$. If ${V_N}^D\le 2{V_0}^D$, then

$$\begin{aligned} 0\le \frac{dv_j}{d{V_0}^D}\le 1, \end{aligned}$$

for all $j=0,\ldots ,N$.

In Lemma 2.2, we rely on the fact that the sequence of voltages $V_0^D,\ldots ,V_N^D$ is increasing. This is established by Lemma 2.3.

Lemma 2.3

The voltages $V_0^D,\ldots ,V_N^D$ form an increasing sequence.

Moreover, in Lemma 2.2, we have restricted ourselves to a network of nodes where the voltage ${V_N}^D$ can maximally be two times the voltage ${V_0}^D$. However, this is not an additional restriction since the inequality $0<\Delta \le \frac{1}{2}$ immediately implies that ${V_N}^D\le 2{V_0}^D$.

Due to Proposition 2.1, we can, given the voltage $V_0^D$, recursively determine the voltage $V_N^D$. The value for the voltage ${V_0}^D$ is arbitrary, so for simplicity, we set ${V_0}^D = 1$ (and thus need to determine whether ${V_N}^D\le \frac{1}{1-\Delta }$). For the rest of the paper, we work with the convention that ${V_0}^D=1$, i.e.,

$$\begin{aligned}&{V_{1}}^D = 1 + rp_0, \end{aligned}$$

(2.17)

$$\begin{aligned}&{V_{j+1}}^D = 2{V_j}^D-{V_{j-1}}^D+\frac{rp_j}{{V_j}^D}, \quad j = 1,\ldots ,N-1 . \end{aligned}$$

(2.18)

Thus, following the expression for the voltage drop constraint in (2.2), the constraint set for the Distflow model reduces to

$$\begin{aligned} {\mathcal {C}}&= \left\{ {\textbf{p}}: {V_N^D}\le \frac{1}{1-\Delta } ,\quad 0<\Delta \le \frac{1}{2}\right\} . \end{aligned}$$

(2.19)

2.3.2 Linearized Distflow

The second power flow model we consider is the Linearized Distflow model. While the first power flow model eventually results in the nonlinear Eq. (5.4) when wishing to compute the voltage at each node, the Linearized Distflow model has a simpler representation. In the Linearized Distflow model, we assume that the active and reactive power losses $r|I_{j-1,j}|^2$ and $x|I_{j-1,j}|^2$, respectively, are much smaller than the active and reactive power flows ${\tilde{P}}_{j-1,j}$ and ${\tilde{Q}}_{j-1,j}$; i.e., the Linearized Distflow approximation neglects the loss terms associated with the squared current magnitudes $|I_{j-1,j}|^2$. Rewriting Eq. (2.8) using Eqs. (2.6), (2.7), and (2.9) under the assumption that $|I_{j-1,j}|^2=0$ for all $\epsilon _{j-1,j}\in {\mathcal {E}}$, and relabeling the resulting equation such that each further term of the sequence is defined as a function of the preceding terms, yields a linear relationship for the squared voltages at each node,

$$\begin{aligned} ({{V'_{j+1}}^L})^2 - ({{V'_{j}}^L})^2&= 2r\sum _{m=0}^{N-1-j} p_m,\quad j=1,\ldots ,N-1. \end{aligned}$$

(2.20)

where $({{V'_N}^L})^2$ is given. Notice that we do not need an equivalence as in the case of the voltages under the Distflow model since we have an explicit expression for ${V'_0}^L$. In this case, the voltage $({V'_0}^L)^2$ can be computed directly; iteratively using (2.20), we find a closed-form solution for the squared voltage ${({{V'_{0}}^L})}^2$, i.e.,

$$\begin{aligned} {({{V'_{0}}^L})}^2 = {({V'_N}^L)}^2-2r\sum _{j=0}^{N-1}\sum _{m=0}^{N-1-j} p_m. \end{aligned}$$

(2.21)

In this case, the value for the voltage $({V'_N}^L)^2$ is arbitrary. To get an equivalent expression for the voltage drop constraint in (2.2) in terms of the squared voltages, we set $({V'_N}^L)^2 = 1/(1-\Delta )^2$ and need to determine whether $({V'_0}^L)^2 \ge 1$. Thus, the constraint set for the Linearized Distflow model reduces to

$$\begin{aligned} {\mathcal {C}}&= \left\{ {\textbf{p}}: \left( {{V_0'}^L}\right) ^2\ge 1 ,\quad 0<\Delta \le \frac{1}{2}\right\} . \end{aligned}$$

(2.22)

Remark 2.1

Note that $V_j, j=1,\ldots ,N$ are dependent on the vector ${\textbf{p}}$ and resistance r. We write $V_j({\textbf{p}},z)$ as a continuous function of the power allocation ${\textbf{p}}$ and distribution network parameter r, respectively, when we wish to emphasize this. Furthermore, the constraint set ${\mathcal {C}}$ is dependent on the specific power flow model that is used. However, we write ${\mathcal {C}}$ as the constraint set for both power flow models as it is clear from the constraints in the set ${\mathcal {C}}$ or from the context which power flow model is used.

3 Main results

In this section, we present our main results, which concern the comparison of the stability of the queuing models under distribution network constraints. The stability of the system is defined as positive recurrence of the Markov process ${\textbf{X}}$.

We compute, in Theorems 3.1 and 3.2, the maximal feasible arrival rates when the number of charging stations, denoted by N, goes to infinity, such that the queuing model is stable given distribution network constraints. In other words, we compute the arrival rates such that the maximal voltage drop is attained. Then, in Sect. 3.3, we compare these arrival rates explicitly.

3.1 Stability conditions under the linearized Distflow model

In this section, we show the computation of the specific arrival rate under the Linearized Distflow model such that the process is stable. That is, we show there exists an explicit threshold $\lambda _N^L$ that only depends on the distribution network parameter r, the maximal voltage drop parameter $\Delta$, and the number of charging stations N in the network such that for all arrival rates below this threshold, the queuing model is stable.

Theorem 3.1

Consider the queuing model in Sect. 2.1 and the power flow model in Sect. 2.3.2. The Markov process ${\textbf{X}}$ is positive recurrent, if

$$\begin{aligned} \lambda < \lambda _N^L := \frac{\frac{1}{r}\left( \left( \frac{1}{1-\Delta }\right) ^2-1 \right) }{N(N+1)}. \end{aligned}$$

(3.1)

Moreover, $N^2\lambda _N^L\rightarrow \lambda _c^L$ as $N\rightarrow \infty$ with

$$\begin{aligned} \lambda _c^L = \frac{1}{r}\left( \left( \frac{1}{1-\Delta } \right) ^2-1 \right) . \end{aligned}$$

(3.2)

The proof is given in Sect. 4. It exploits the explicit expression of the squared voltage ${W_{0,0}^L}$ and uses [26, Theorem 11]. We show for the arrival rate $\lambda _N^L$ that the maximal allowed voltage drop is attained, so that $\varvec{\lambda }_N^L$ is still in the constraint set ${\mathcal {C}}$.

Remark 3.1

Under the Linearized Distflow power flow model, defined in (4.1), the maximum feasible arrival rate is decaying as the inverse of the square of the number of nodes.

We now proceed with the stability result under the Distflow model.

3.2 Stability conditions under the Distflow model

In this section, we show there exists a specific arrival rate $\lambda _N^D$ that only depends on the resistance r, the maximal voltage drop parameter $\Delta$, and the number of charging stations N in the network such that for every arrival rate below this threshold the system is stable. Furthermore, as the number of charging stations N goes to infinity, we show the convergence of a suitably scaled version of $\lambda _N^D$ to a critical arrival rate $\lambda _c^D$.

Theorem 3.2

Consider the queuing model in Sect. 2.1 and the power flow model in Sect. 2.3.1. There exists $\lambda _N^D$ such that the Markov process ${\textbf{X}}$ is positive recurrent, if

$$\begin{aligned} \lambda < \lambda _N^D. \end{aligned}$$

Moreover, $N^2\lambda _N^D \rightarrow \lambda _c^D$ as $N\rightarrow \infty$ with

$$\begin{aligned} \lambda _c^D = \frac{\pi }{2r}\text {erfi}^2\left( \sqrt{\ln \left( \frac{1}{1-\Delta }\right) } \right) . \end{aligned}$$

(3.3)

The imaginary error function $\text {erfi}(z)$ is the function defined by

$$\begin{aligned} \text {erfi}(z) = \frac{2}{\sqrt{\pi }}\int _0^z \exp (v^2)dv = \frac{1}{\textrm{i}}\text {erf}(\textrm{i}z), \end{aligned}$$

where $\text {erf}(x) = \frac{2}{\sqrt{\pi }}\int _0^x \exp (-t^2)dt$ is the standard error function.

The proof is given in Sect. 5. Unlike the proof of Theorem 3.1, the key difference is that we do not have an explicit expression for the voltage $V_N^D$. Thus, we require a different approach to find the maximal feasible arrival rate $\lambda _N^D$. Instead of solving for the value $\lambda _N^D$ such that the maximal allowed voltage drop is attained, we approximate the voltage $V_N^D$ by a continuous function. Then, we show that this continuous function converges to the voltage $V_N^D$ as N goes to infinity and compute the arrival rate such that the maximal voltage drop is attained under this new approximation.

From Theorems 3.1 and 3.2, we see that the arrival rates for both power flow models decay at the same rate if the number of charging stations grows to infinity. However, this does not provide insights into the actual difference between the critical arrival rates if the number of charging stations grows to infinity. In the next section, we quantify this difference explicitly.

3.3 Comparison of stability regions

In this section, we compare the critical arrival rates under both power flow models. The expression for the critical arrival rates under the Linearized Distflow and Distflow model is given in (3.2) and (3.3), respectively. Then, we have for the ratio between $\lambda _c^D$ and $\lambda _c^L$,

$$\begin{aligned} \frac{\lambda _c^D}{\lambda _c^L}&= \frac{2(1-\Delta )^2\left( \int _0^{\sqrt{\ln \left( \frac{1}{1-\Delta }\right) }}\exp (u^2)du \right) ^2}{\Delta (2-\Delta )} \nonumber \\&=: P(\Delta ). \end{aligned}$$

(3.4)

We study the behavior of the function $P(\Delta )$ in Theorem 3.3.

Theorem 3.3

Let $P(\Delta )$ be defined as in (3.4) for $0<\Delta \le \frac{1}{2}$. Then, $P(\Delta )$ is a strictly decreasing function from 1 at $\Delta =0$ to $\frac{\pi }{6}\text {erfi}\left( \sqrt{\ln (2)} \right) ^2\approx 0.77$ at $\Delta =\frac{1}{2}$.

The proof of Theorem 3.3 is found in Appendix B. The importance of Theorem 3.3 lies in the fact that it implies that the critical arrival rate under the Distflow is always smaller than the critical arrival rate under the Linearized Distflow model and that we are able to quantify the difference between the critical arrival rates.

It is shown that the Linearized Distflow model overestimates the voltages in comparison with the voltages given by the Distflow model [18, Lemma 12]. Therefore, any feasible power allocation under the Linearized Distflow model is also a feasible power allocation under the Distflow model. Hence, given service requirements, any stable arrival rate under the Linearized Distflow model is also a stable arrival rate under the Distflow model.

For several values of $\Delta$, Table 1 shows the ratio between $\lambda _c^D$ and $\lambda _c^L$.

Table 1 The fraction $\lambda _c^D/\lambda _c^L$ for specific values of $\Delta$

Fig. 2

The fraction $\frac{\lambda _c^D}{\lambda _c^L}$ via (3.4) for $0<\Delta \le \frac{1}{2}$

Similar to Table 1, one can see from Fig. 2 that the ratio between $\lambda _c^D$ and $\lambda _c^L$ decreases as $\Delta$ increases. When we allow a voltage drop of $50\%$, i.e., $\Delta = \frac{1}{2}$, the maximal feasible arrival rate under the Linearized Distflow model is $\approx 20\%$ higher than under the Distflow model. In more realistic situations, such as in Dutch distribution networks, when the relative drop between voltages cannot be more than $4.5\%$, or in other words, when $\Delta = 0.045$, the difference between these critical arrival rates is even smaller.

4 Proof of Theorem 3.1

Our work is based on [2] and [26]. Both papers provide a way to prove stability. So, before we continue with the proofs on stability in Sects. 4 and 5, we discuss these two ways to prove stability. First, we start with a brief discussion of the approach in [2]—the method of convex relaxation. Then, we study an extension of Shneer and Stolyar [26, Theorem 11] that we use to prove Theorems 3.1 and 3.2 and provides stability of the Markov process ${\textbf{X}}$ as well.

The approach of convex relaxation in [2] goes as follows. It is known that the utility optimization problem in (2.1) subject to the constraint set defined in (2.19) is in general non-convex. For non-convex sets, the choice of the utility function is crucial to ensure stability [6]. However, since our underlying network topology is a line, we obtain that the convex relaxation of (2.1) under the constraint defined in (2.19) is exact [2]. When the constraint set is convex, the stability of utility-maximizing allocations is known to be independent of the considered utility function [6]. Then, also for all $\alpha$-fair allocations [c.f. (2.1)], an allocation can be given such that the Markov process ${\textbf{X}}$ is stable [6].

The approach in [26] is more general in the sense that it does not assume $\alpha$-fair allocations or convexity to prove stability. It holds for a larger class of functions and for non-convex constraint sets. If the constraint set ${\mathcal {C}}$ is compact and coordinate-convex, then the Markov process ${\textbf{X}}$, which represents the number of EVs charging at every station, is stable if there exists a vector $\varvec{\nu }\in {\mathcal {C}}$ such that $\varvec{\lambda }<\varvec{\nu }$ according to [26]. In fact, in Appendix C, we show for general arrival and service rates that we do not need coordinate-convexity, but only compactness of the constraint set ${\mathcal {C}}$ to prove stability. The fact that we do not need coordinate-convexity is important, since we can construct examples where this property does not hold for the constraint set of the Distflow model [cf. (2.19)].

Therefore, the proof consists of three steps. First, we show that our queuing model description, as in Sect. 2.1, is contained in the general model framework that is given by [26]. Then, in order to apply their stability result, we show that the constraint set ${\mathcal {C}}$ is compact. Last, we show that the maximal feasible arrival rate $\lambda _N^L$ is given by Eq. (3.1).

4.1 Queuing model framework

We consider the general model framework in the continuous-time setting in [26, Section 5.1] and compare it with the continuous-time model of Sect. 2. In [26], the authors consider a finite network of queues such that individual instantaneous service rates may depend on the state of the entire system. This corresponds to our network of charging stations such that the power allocation depends on the total number of cars charging at the same time. The network state represents a set ${\textbf{X}}(t)=(X_i(t),i\in N)$ of the queue lengths $X_i$ at the network nodes $i\in N$ at time t. In our case, the network state is represented by the number of cars charging at each charging station. Arrivals into queue i occur according to a Poisson process with a constant rate $\lambda$, independent of all other processes. The instantaneous service rate $\mu _i$ of a node represents the intensity of the Poisson process modeling departures (service completions) of the node. In [26], the authors consider a service rate allocation algorithm $\varvec{\psi }({\textbf{X}}(t))$, mapping a network state ${\textbf{X}}(t)$ into a set of instantaneous service rates $\varvec{\mu }$, is all that is needed to specify the service allocation algorithm. In our case, the power allocation is given by a constrained optimization problem, mapping the network state, the number of EVs at each charging station, into a power allocation.

In this general model framework, if there exists $\varvec{\nu }\in {\mathcal {C}}$ such that $\varvec{\lambda }<\varvec{\nu }$ and ${\mathcal {C}}$ is compact and coordinate-convex (i.e., if a vector $\mu$ belongs to ${\mathcal {C}}$ and $\mu '\le \mu$ coordinate-wise, then $\mu '\in {\mathcal {C}}$), then the Markov process ${\textbf{X}}$ is positive recurrent according to Shneer and Stolyar [26, Theorem 11]. In Appendix C, we show that the set ${\mathcal {C}}$ does not need to be coordinate-convex, since we can construct a coordinate-convex set ${\mathcal {C}}_1$ such that the following are equivalent:

$$\begin{aligned} \text {there exists}\ \varvec{\nu }\in {\mathcal {C}}\ \text {such that}\ \varvec{\lambda }<\varvec{\nu }, \end{aligned}$$

and

$$\begin{aligned} \text {there exists}\ \varvec{\nu }\in {\mathcal {C}}_1\ \text {such that}\ \varvec{\lambda }<\varvec{\nu }. \end{aligned}$$

We are thus left to show that there exists $\varvec{\nu }\in {\mathcal {C}}$ such that $\varvec{\lambda }<\varvec{\nu }$ and that ${\mathcal {C}}$ is a compact set. First, we show the compactness of the set ${\mathcal {C}}$.

4.2 Compactness of the constraint set ${\mathcal {C}}$

To simplify notation, it is easier to work with the squared voltages to prove compactness of the constraint set ${\mathcal {C}}$, defined in (2.22). Thus, we may rewrite (2.21) by defining ${W_{j,j}}^L:=({{V'_j}^L})^2$, for $j=0,\ldots ,N$. Then, (2.21) has the linear expression,

$$\begin{aligned} {W_{0,0}}^L = {W_{N,N}}^L - 2r\sum _{j=0}^{N-1}\sum _{m=0}^{N-1-j} p_m. \end{aligned}$$

(4.1)

In this section, we show that the constraint set ${\mathcal {C}}$ is compact, i.e., we show that the set is closed and bounded.

Lemma 4.1

The constraint set ${\mathcal {C}}$, defined in (2.22), is compact for the Linearized Distflow model.

Proof

Consider N fixed. Recall from the relabeling for the Linearized Distflow model that $W_{N,N}^L$ is set to a fixed value. Define ${\textbf{C}} = -2r\begin{pmatrix} N&\ldots&1\end{pmatrix}\in {\mathbb {R}}^{1\times N}$ and ${\textbf{p}}=\begin{pmatrix} p_0&\ldots&p_{N-1}\end{pmatrix} \in {\mathbb {R}}^{N\times 1}$. Then, one can derive from (4.1) that $W_{0,0}^L$ is given by an affine transformation $l:{\mathbb {R}}^N \rightarrow {\mathbb {R}}$ where $l({\textbf{p}})={\textbf{C}}{\textbf{p}}+W_{N,N}^L$. Furthermore, let $\epsilon >0$ and choose $\delta = \epsilon /\Vert {\textbf{C}}\Vert$. Then, for ${\textbf{p}}^*\in {\mathbb {R}}^{N\times 1}$ such that $\Vert {\textbf{p}}-{\textbf{p}}^*\Vert <\delta$, we have

$$\begin{aligned} \Vert l({\textbf{p}})-l({\textbf{p}}^*)\Vert&= \Vert ({\textbf{C}}{\textbf{p}}+d)-({\textbf{C}}{\textbf{p}}^*+d)\Vert = \Vert {\textbf{C}}({\textbf{p}}-{\textbf{p}}^*)\Vert \\&\le \Vert {\textbf{C}}\Vert \Vert {\textbf{p}}-{\textbf{p}}^*\Vert \\&< \Vert {\textbf{C}}\Vert \delta = \epsilon . \end{aligned}$$

This implies that $W_{0,0}^L$ as a function of ${\textbf{p}}$ is continuous. Therefore, the constraint set ${\mathcal {C}}$ is closed. Furthermore, by definition of the constraint set in (2.22), we have

$$\begin{aligned} {W_{0,0}^L}&= {W_{N,N}^L}-2r\sum _{j=0}^{N-1}\sum _{m=0}^{N-1-j} p_m\\&= \left( \frac{1}{1-\Delta } \right) ^2-2r\sum _{j=0}^{N-1}\sum _{m=0}^{N-1-j} p_m \ge 1. \end{aligned}$$

Thus,

$$\begin{aligned} \sum _{j=0}^{N-1}\sum _{m=0}^{N-1-j} p_m \le \frac{1}{2r}\left( \left( \frac{1}{1-\Delta } \right) ^2-1\right) . \end{aligned}$$

(4.2)

Therefore, ${\textbf{p}}$ is bounded. Hence, the set of constraints ${\mathcal {C}}$ is bounded and closed, which implies compactness. $\square$

Next, we show there exists $\varvec{\nu }\in {\mathcal {C}}$ such that $\varvec{\lambda }<\varvec{\nu }$.

4.3 Computation of $\lambda _N^L$ such that maximal voltage drop is attained

In this section, we compute the maximal feasible arrival rate $\lambda _N^L$ such that $\varvec{\lambda }_N^L = (\lambda _N^L,\ldots ,\lambda _N^L)\in {\mathcal {C}}$ (where we define $\lambda _N^L$ below). Then, by Shneer and Stolyar [26, Theorem 11], for all $\lambda <\lambda _N^L$, the Markov process ${\textbf{X}}$ is positive recurrent. First, we make the following remark.

Remark 4.1

We do not necessarily need that all arrival rates for all nodes are the same. It is sufficient to find a vector $\varvec{\lambda }^*$ such that (4.2) is tight because for this vector $\varvec{\lambda }^*$ the maximal voltage drop is attained. In this case, the Markov process ${\textbf{X}}$ is positive recurrent if $\varvec{\lambda } < \varvec{\lambda }^*$ for some arrival rate vector $\varvec{\lambda }$. However, to make a comparison between the two power flow models, we eliminate all sources of noise by letting all arrival rates for all nodes be the same, i.e., $\varvec{\lambda }_N^L = (\lambda _N^L,\ldots ,\lambda _N^L)$.

Now, we establish that $\varvec{\lambda }_N^L$ is contained in the constraint set ${\mathcal {C}}$. By definition, we have

$$\begin{aligned} \lambda _N^L = \frac{\left( \frac{1}{1-\Delta } \right) ^2-1}{rN(N+1)}. \end{aligned}$$

Using this arrival rate in the expression for the squared voltage $W_{0,0}^L(\varvec{\lambda }_N^L)$ [cf. Equation (4.1)] yields

$$\begin{aligned} W_{0,0}^L(\varvec{\lambda }_N^L)&= W_{N,N}^L-2r\sum _{j=0}^{N-1}\sum _{m=0}^{N-1-j} \lambda _N^L \\&= \left( \frac{1}{1-\Delta } \right) ^2-2r\lambda _N^L \frac{1}{2}N(N+1)\\&= 1. \end{aligned}$$

Hence, the vector $\varvec{\lambda }_N^L = (\lambda _N^L,\ldots ,\lambda _N^L)$ is contained in the constraint set ${\mathcal {C}}$ and for $\lambda <\lambda _N^L$ the Markov process ${\textbf{X}}$ is positive recurrent.

5 Proof of Theorem 3.2

The proof consists of four steps. Similar to the proof of Theorem 3.1, we show that our queuing model description, as in Sect. 2.3.1, is contained in the general model framework that is given by Shneer and Stolyar [26] (see Appendix C) and that the constraint set ${\mathcal {C}}$ is compact for the Distflow model. Then, we need an intermediate step before we can compute the arrival rate $\lambda _N^D$ such that the maximal voltage drop is attained. Unlike the proof of Theorem 3.1, we do not have an explicit expression for the voltage $V_N^D$. Therefore, our next step is to approximate the voltage $V_N^D$ by a continuous function, as $N\rightarrow \infty$. In Appendix D, we show the convergence of the approximation of the voltage to the actual value $V_N^D$, as $N\rightarrow \infty$. This approximation allows us to compute the arrival rate $\lambda _N^D$ such that the maximal voltage drop is attained, which concludes the proof.

5.1 Queuing model framework

Again, we consider the continuous-time model as in Sect. 2 and we follow the exact same reasoning as in Sect. 4.1. The only difference in this description is given by the constrained optimization problem. In this case, the voltages are computed according to the Distflow model instead of the Linearized Distflow model.

5.2 Compactness of the constraint set ${\mathcal {C}}$

It is again easier to work with the squared voltages to prove the compactness of the constraint set ${\mathcal {C}}$ for the Distflow model. Therefore, for $\epsilon _{j-1,j}\in {\mathcal {E}}$, we apply the transformation,

$$\begin{aligned} {\textbf{W}}(\epsilon _{j-1,j}) = \begin{pmatrix} {{V_{j-1}}}^2 &{} V_{j-1}V_j \\ V_jV_{j-1} &{} {{V_j}}^2 \end{pmatrix}:= \begin{pmatrix} W_{j-1,j-1} &{} W_{j-1,j} \\ W_{j,j-1} &{} W_{jj} \end{pmatrix} \end{aligned}$$

(5.1)

to (2.17) and (2.18). This transformation, for $j=1,\ldots ,N-1$, leads to the linear Equations (5.2) and (5.3) (in terms of ${\textbf{W}}(\epsilon _{ij})$), i.e.,

$$\begin{aligned} {W_{j,j+1}}^D&= 2{W_{j,j}}^D-{W_{j-1,j}}^D+rp_j, \end{aligned}$$

(5.2)

where

$$\begin{aligned} {W_{0,0}}^D=1\quad \text {and}\quad {W_{0,1}}^D = 1+rp_0 \end{aligned}$$

(5.3)

are the initial values. However, if we want to compute the recursion variables ${W_{j,j+1}}^D, j=1,\ldots ,N-1$ explicitly, we need the relationship defined in (5.4), which is not linear:

$$\begin{aligned} W_{j,j}^D&= \frac{{{W_{j-1,j}}^D}^2}{{W_{j-1,j-1}}^D}. \end{aligned}$$

(5.4)

We want to keep the maximal voltage drop under control. Using the transformation, the voltage drop constraint reduces to

$$\begin{aligned} W_{N,N}^D\le \left( \frac{1}{1-\Delta } \right) ^2, \end{aligned}$$

(5.5)

where $0<\Delta \le \frac{1}{2}$.

To establish compactness of ${\mathcal {C}}$ in the case of the Distflow model, we use the fact that the voltages $V_0^D,\ldots ,V_N^D$ form an increasing sequence. This is proven in Lemma 2.3. Now, we present an analogous result in Corollary 5.1.

Corollary 5.1

The squared voltages $W_{0,0}^D, W_{0,1}^D, W_{1,1}^D,\ldots ,W_{N-1,N}^D,W_{N,N}^D$ are an increasing sequence.

Proof

This is an immediate consequence of Lemma 2.3. Since $V_j^D\ge 1$, for all $j=0,\ldots ,N$, we have

$$\begin{aligned} W_{j,j}^D = V_{j}^DV_{j}^D \ge V_{j-1}^DV_{j}^D = W_{j-1,j}^D \end{aligned}$$

for all $j=0,\ldots ,N$. $\square$

Now, we are in position to prove Lemma 5.1.

Lemma 5.1

The constraint set ${\mathcal {C}}$, defined in (2.19), is compact for the Distflow model.

Proof

Since $W_{N,N}^D({\textbf{p}},z)$ is continuous, the set ${\mathcal {C}}$ is closed. By the transformation in (5.1), we have

$$\begin{aligned} \left( \frac{1}{1-\Delta } \right) ^2 \ge W_{N,N}^D({\textbf{p}},z). \end{aligned}$$

Then, we show a lower bound on the squared voltage $W_{N,N}^D$. By Corollary 5.1 and Eq. (5.2), we have

$$\begin{aligned} {W_{N,N}^D}({{\textbf{p}}},z)&\ge {W_{N-1,N}^D}({{\textbf{p}}},z) \nonumber \\&= 2{W_{N-1,N-1}^D}({{\textbf{p}}},z)-{W_{N-2,N-1}^D}({{\textbf{p}}},z)+rp_{N-1} \nonumber \\&\ge 2{W_{N-2,N-1}^D}({{\textbf{p}}},z)-{W_{N-2,N-1}^D}({{\textbf{p}}},z)+rp_{N-1}\nonumber \\&= {W_{N-2,N-1}^D}({{\textbf{p}}},z) +rp_{N-1}. \end{aligned}$$

(5.6)

Now, apply the definition of $W_{N-2,N-1}^D({\textbf{p}},z)$ in Eq. (5.2) and Corollary 5.1 to Eq. (5.6) to find the inequality

$$\begin{aligned} {W_{N-2,N-1}^D}({{\textbf{p}}},z) \ge {W_{N-3,N-2}^D}({{\textbf{p}}},z)+rp_{N-2}. \end{aligned}$$

(5.7)

Inserting (5.7) in (5.6) gives

$$\begin{aligned} {W_{N,N}^D}({{\textbf{p}}},z)&\ge {W_{N-3,N-2}^D}({{\textbf{p}}},z)+rp_{N-2}+rp_{N-1}. \end{aligned}$$

Applying the definition of ${W_{j-1,j}^D}({{\textbf{p}}},z)$ and Corollary 5.1 over and over again, for $j=1,\ldots ,N-3$ in reversed order, results in the inequality

$$\begin{aligned} {W_{N,N}^D}({{\textbf{p}}},z)&\ge 1+r\left( p_0+\cdots +p_{N-1} \right) . \end{aligned}$$

Thus,

$$\begin{aligned} \left( p_0+\cdots +p_{N-1} \right) \le \frac{1}{r}\left( \frac{1}{1-\Delta } \right) ^2. \end{aligned}$$

Therefore, ${\textbf{p}}$ is bounded. Hence, the set of constraints ${\mathcal {C}}$ is bounded and closed, which implies compactness. $\square$

5.3 Approximation of $V_N^D$ by continuous counterpart

Unlike the proof of Theorem 3.1, we do not have an explicit expression for the voltage $V_N^D$. Therefore, we define a continuous extension of the voltages $V_j^D,\ j=0,1,\ldots ,N$ to a function $V_N^D(t): [0,N] \rightarrow {\mathbb {R}}_+$. That is, for all $j\in \{0,\ldots ,N\}$, we show $V_{j}^D = V_N^D(j)$. A suitable and natural extension, cf. (2.16), is given by

$$\begin{aligned} V_N^D(t)&:= 1+k_N \int _0^{\lfloor t\rfloor }\int _0^{\lceil s\rceil }\frac{1}{V_N^D({u})} du\ ds,\ t\in [0,N]; \end{aligned}$$

(5.8)

where $k_N = r\lambda _N^D$. Here, the arrival rate $\lambda _N^D$ is the same for all charging stations. We include N in the notation of the function $V_N^D(t)$ and $k_N$ to stress its dependence on the number of charging stations N.

Then, we scale the continuous extension $V_N^D(t)$ in time by N and in space by scaling the parameter $k_N$ by $N^2$, such that $k_N = \frac{a}{N^2}$. Therefore, we use the following definition of the in time-scaled continuous extension of $V_N^D(t)$.

Definition 1

Denote the voltages under the Distflow model in continuous time by $V_N^D(t)$ as in (5.8), and scale them in time by N,

$$\begin{aligned} {\overline{V}}_N^D(t) := V_N^D(Nt), \ t\in [0,1]. \end{aligned}$$

(5.9)

This scaling allows us to show (see Sect. D.1.1) that the representation of the scaled version of the voltages ${\overline{V}}_N^D(t)$ can be written as

$$\begin{aligned} {\overline{V}}_N^D(t)&= 1+a\int _0^{t}\int _0^{s}\frac{1}{{\overline{V}}_N^D(u)}du\ ds +{\overline{R}}_N(t), \end{aligned}$$

(5.10)

where the remainder term ${\overline{R}}_N(t)$ vanishes as $N\rightarrow \infty$. Hence, (5.10) can be expected to converge to the integral equation given by

$$\begin{aligned} V(t) = 1+a\int _0^t\int _0^x \frac{1}{V(y)}dy\ dx. \end{aligned}$$

(5.11)

This convergence is exactly shown in Proposition 5.1 and proven in Sect. D.1.

Proposition 5.1

Let ${\overline{V}}_N^D(t)$, for $t\in [0,1]$, as defined in (5.10) and V(t) as in (5.11). Then,

$$\begin{aligned} \sup _{0\le t\le 1} |{\overline{V}}_N^D(t)-V(t)|\rightarrow 0\ \text {as}\ N\rightarrow \infty . \end{aligned}$$

Thus, the appropriately scaled voltages under the Distflow model ${\overline{V}}_N^D(t)$ can be approximated by V(t) as the number of charging stations goes to infinity. The solution to the integral equation V(t) is found in Appendix E. Using this approximation, we are able to compute the arrival rate $\lambda _N^D$ such that the maximal allowed voltage drop is attained.

5.4 Computation of $\lambda _N^D$ such that maximal voltage drop is attained

In this section, we compute the arrival rate $\lambda _N^D$ such that the maximal allowed voltage drop is attained. First, we show that we can approximate $V_N^D$ by V(1) [cf. (5.11)]. Then, using this approximation, we compute the value a, which is related to the arrival rate $\lambda _N^D$ [cf. (5.14)], that solves the equation $V(1)=\frac{1}{1-\Delta }$. In other words, we compute the value a such that the maximal voltage drop is attained. However, this solution a is computed using the approximation V(1). So, we conclude by showing that this solution a also maximizes the voltage drop using $V_N^D$.

The recursion for the voltages under the Distflow model [cf. Eqs. (2.17) and (2.18)] with arrival rate $\lambda _N^D$, is given by,

$$\begin{aligned} V_0^D&= 1\ \text {and}\ V_1^D = 1 + k_N, \end{aligned}$$

(5.12)

$$\begin{aligned} V_{j+1}^D&= 2V_j^D-V_{j-1}^D+\frac{k_N}{V_j^D}, \quad j = 1,\ldots ,N-1. \end{aligned}$$

(5.13)

where

$$\begin{aligned} k_N=r\lambda _N^D = \frac{a}{N^2}\ge 0. \end{aligned}$$

(5.14)

Then, according to the voltage drop constraint in (2.19), the maximum feasible arrival $\lambda _N^D$ is such that the voltage $V_N^D$ is equal to $1/(1-\Delta )$. In Sect. 5.3, we found an equivalent expression for the voltage $V_N^D$, namely ${\overline{V}}_N^D(1)$. Furthermore, using Proposition 5.1, we found an approximation of the voltage $V_N^D$, namely V(1), as in (5.11).

Differentiating (5.11) twice gives the differential equation form of the integral equation (5.11), i.e.,

$$\begin{aligned} V''(t) = \frac{a}{V(t)}, \end{aligned}$$

(5.15)

for $a\in {\mathbb {R}}_+$, with initial conditions

$$\begin{aligned} V(0)=1\ \text {and}\ V'(0)=0. \end{aligned}$$

(5.16)

The solution to the differential equation (5.15) and (5.16) is given in terms of the function $f_0(x)=\exp (U^2(x))$ where U(x) is related to the inverse of the imaginary error function. We prove this result in Lemma E.1. Moreover, it follows that

$$\begin{aligned} V(t) = f_0(t\sqrt{a}). \end{aligned}$$

Therefore, an approximation of the voltage $V_N^D$, namely V(1), is given by,

$$\begin{aligned} V(1)=f_0(\sqrt{a}). \end{aligned}$$

(5.17)

Although it is possible to compute the arrival rate $\lambda _N^D$ such that $V_N^D = \frac{1}{1-\Delta }$ using an iterative approach, we can obtain a convenient approximation by using (5.17). In Sect. 6, we show how to compute the arrival rate $\lambda _N^D$ using Newton’s method and that approximation (5.17) yields already a good approximation for the arrival rate $\lambda _N^D$.

With the approximation V(1) (as $N\rightarrow \infty$) at hand, we compute the value of a that solves the equation $V(1)=\frac{1}{1-\Delta }$, i.e., $f_0(\sqrt{a}) = \frac{1}{1-\Delta }$ and use that $r\lambda _N^D = \frac{a}{N^2}$ [cf. (5.14)] to relate it to a critical arrival rate. Solving for a such that $f_0(\sqrt{a}) = \frac{1}{1-\Delta }$ yields,

$$\begin{aligned} a = \left( f_0^{-1}\left( \frac{1}{1-\Delta }\right) \right) ^2. \end{aligned}$$

(5.18)

Notice that the inverse function $f_0^{-1}$ is given by

$$\begin{aligned} f_0^{-1}(x) = \sqrt{\frac{\pi }{2}}\text {erfi}\left( \sqrt{\ln (x)} \right) . \end{aligned}$$

Finally, we use [19, Lemma 3.1] to conclude that a in (5.18) also maximizes the voltage drop using $V_N^D$. In [19, Lemma 3.1], the authors consider a sequence $R_n, n=1,2,\ldots$ on a Banach space X, such that $R_n(x)\rightarrow R(x)$ uniformly for $x\in U$, where U is a normed space. Then, the following holds:

• 1. $\inf \left\{ R_n(x)\ \text {s.t.}\ x\in U\right\} \rightarrow \inf \left\{ R(x)\ \text {s.t.}\ x\in U\right\}$ as $n\rightarrow \infty$,

and

• 2. If $x_n$ is a solution to the optimization problem

  $$\begin{aligned} \min \left\{ R_n(x)\ \text {s.t.}\ x\in U\right\} , \end{aligned}$$

  then $x_n$ is also a solution to the optimization problem

  $$\begin{aligned} \min \left\{ R(x)\ \text {s.t.}\ x\in U\right\} , \end{aligned}$$

  (5.19)

  and for a continuous function R is every limit point in U of $x_n$ a solution to (5.19).

Applying this result, we have, for any fixed a, by Theorem 5.1 that $|{\overline{V}}_N^D(1)-V(1)|\rightarrow 0$ uniformly as $N\rightarrow \infty$ and by construction that $V(1)=\frac{1}{1-\Delta }$. Now, let

$$\begin{aligned} \phi _N := \max \Big \{{\overline{V}}_N^D(1)\ \text {s.t.}\ {\overline{V}}_N^D(1)\le \frac{1}{1-\Delta }\Big \}, \end{aligned}$$

(5.20)

and

$$\begin{aligned} \phi := \max \Big \{V(1)\ \text {s.t.}\ V(1)\le \frac{1}{1-\Delta }\Big \}, \end{aligned}$$

denote the maximal voltages and the approximation of the maximal voltages under the voltage drop constraint, respectively. Then, by [19, Lemma 3.1], $\phi _N$ converges to $\phi$ as $N\rightarrow \infty$ and since $a=N^2r\lambda _N^D$ maximizes V(1), a maximizes ${\overline{V}}_N^D(1)$ as well and the limit point $r\lambda _c^D$ of $N^2r\lambda _N^D$ is a solution to (5.20), i.e., $N^2r\lambda _N^D \rightarrow r\lambda _c^D = \frac{\pi }{2}\text {erfi}^2\left( \sqrt{\ln \left( \frac{1}{1-\Delta } \right) }\right)$, as $N\rightarrow \infty$.

6 Newton’s iterations for the Distflow model

In Sects. 5.3 and 5.4, we approximated the voltage $V_N^D$ to make the computation of the arrival rate $\lambda _N^D$ such that $V_N^D=\frac{1}{1-\Delta }$ easier. However, it is possible to compute the arrival rate $\lambda _N^D$ such that $V_N^D=\frac{1}{1-\Delta }$, iteratively, by using the recursive Eqs. (5.12) and (5.13). In other words, for a fixed value of the number of charging stations N and a maximum allowed voltage drop $\Delta >0$, we compute the number a such that the voltage drop between the root node and node N is exactly equal to $\Delta$. Furthermore, we present some numerical tests showing the convergence to the number a.

In what follows, we use Newton’s method to find the solution a that yields $V_N^D=\frac{1}{1-\Delta }$. We need a particular $a_0$ to initialize Newton’s method. Here, we are led by Theorem 5.1. We look for an initial guess $a_0$ such that $V(1) = f_0\left( \sqrt{a}\right) = \frac{1}{1-\Delta }$, we choose

$$\begin{aligned} a_0&= \left( f_0^{-1}\left( \frac{1}{1-\Delta } \right) \right) ^2 \nonumber \\&= \frac{\pi }{2}\text {erfi}^2\left( \sqrt{\ln \left( \frac{1}{1-\Delta }\right) } \right) , \end{aligned}$$

and iterate according to

$$\begin{aligned} a_{j+1}=a_j-\frac{V_N^D-(\frac{1}{1-\Delta })}{Y_N^D/N^2},\quad j=0,1,\ldots \end{aligned}$$

(6.1)

where $Y_N^D$ denotes the derivative of $V_N^D$ with respect to $k_N$. Note that $Y_N^D/N^2$ denotes the derivative of $V_N^D$ with respect to a. We need to compute $Y_N^D$, and this can be done by differentiating in (5.12) and (5.13) with respect to $k_N$, and so that

$$\begin{aligned} Y_0^D=0, Y_1^D=1 \end{aligned}$$

(6.2)

and

$$\begin{aligned} Y_{n+1}^D - 2Y_n^D+Y_{n-1}^D = \frac{1}{V_n^D}-\frac{k_NY_n^D}{(V_n^D)^2}. \end{aligned}$$

(6.3)

Observe that the right-hand side of (6.3) also involves $V_n^D$. Thus, one should compute $V_n^D$ and $Y_n^D$ simultaneously and recursively according to the initialization in (5.12) and (6.2), and the recursion step for $n=1,\ldots ,N-1$ in (5.13) and (6.3). Furthermore, the iteration step in (6.1) is well defined if $Y_N^D$ is positive. This is, for reasonable values of a, shown in Lemma 6.1. Before we move to Lemma 6.1, we comment on the range of reasonable values for a. Notice, from the condition that $f_0(\sqrt{a}) = \frac{1}{1-\Delta }$, it follows that

$$\begin{aligned} \int _0^{\sqrt{\ln (\frac{1}{1-\Delta })}}\exp (u^2)du = \sqrt{\frac{a}{2}}. \end{aligned}$$

For small $\Delta$, which is usually the case (cf. Sect. 2.2), it follows from a Taylor expansion about $\Delta =0$ that $\Delta \approx \frac{a}{2}$, which is also seen in Tables 2, 3, and 4. So as long as $0<\Delta <1$, which is the case, we have $0<a<2$. Therefore, if we consider $0<a<2$, $Y_N^D$ is positive and the iteration step is well defined.

Lemma 6.1

Consider the recursion in (6.2) and (6.3) with $0<a<2$. Then, the sequence $Y_0,Y_1,\ldots ,Y_N$ is positive, increasing, and convex.

Proof

The right-hand side of (6.3) is positive for all $n=1,2,\ldots ,N-1$ such that $Y_n^D < V_n^D/k_N$. Therefore, since $Y_0^D = 0 < 1 = Y_1$, the sequence $Y_n^D, n=1,2,\ldots$ is increasing and convex (and thus positive) as long as $Y_n^D < V_n^D/k_N$. For such n, we have certainly

$$\begin{aligned} Y_{n+1}^D-2Y_n^D+Y_{n-1}^D = \frac{1}{V_n^D}\left( 1-\frac{k_NY_n^D}{V_n^D} \right) \le 1, \end{aligned}$$

since

$$\begin{aligned} \frac{k_NY_n^D}{(V_n^D)^2}<\frac{1}{V_n^D}\le 1. \end{aligned}$$

Here, we used Lemma 2.3. So, using $Y_0^D=0,Y_1^D=1$ and induction, we get

$$\begin{aligned} Y_{n+1}^D-Y_n^D \le (n+1), \end{aligned}$$

and consequently for these n, by summing from $m=0,\ldots ,n$,

$$\begin{aligned} Y_{n+1}^D = \sum _{m=0}^n Y_{m+1}^D-Y_m^D \le \sum _{m=0}^n m+1 = \frac{1}{2}(n+1)(n+2). \end{aligned}$$

Therefore, the condition $Y_n^D<V_n^D/k_N$ continues to be satisfied certainly when

$$\begin{aligned} \frac{1}{2}n(n+1)<V_n^D/k_N. \end{aligned}$$

(6.4)

The right-hand side of (6.4) is at least equal to $N^2/a$, and the left-hand side of (6.4) is at most equal to $\frac{1}{2}(N-1)N$, since $n=0,\ldots ,N-1$. Thus, (6.4) holds for all $n=1,2,\ldots ,N-1$ when $\frac{1}{2}(N-1)< N/a$, i.e., when

$$\begin{aligned} a < \frac{2N}{N-1}. \square \end{aligned}$$

Table 2 Illustration of Newton scheme in the computation of the solution ${\bar{a}}$ such that $V_N^D=\frac{1}{1-\Delta }$, for $a=0.1$

Table 3 Illustration of Newton scheme in the computation of the solution ${\bar{a}}$ such that $V_N^D=\frac{1}{1-\Delta }$, for $a=0.2$

Table 4 Illustration of Newton scheme in the computation of the solution ${\bar{a}}$ such that $V_N^D=\frac{1}{1-\Delta }$, for $a=1$

6.1 Numerical validation of Newton’s method

Below, we summarize some numerical tests. For fixed values of a, ranging from 0.01 to 1 and several values of $N=10,\ldots ,10^5$, we compute $V_N^D$ and we put

$$\begin{aligned} \frac{1}{1-\Delta } = V_N^D, \end{aligned}$$

and using the Newton scheme above, we compute the solution ${\bar{a}}$ such that $V_N^D = \frac{1}{1-\Delta }$ with stopping criteria $\frac{|a_{j+1}-a_j|}{a_j}<10^{-10}$.

It appears that the initialization $a_0$ is already a good approximation for a as the differences between $a_0$ and a for all situations are small. For example, the relative error of the initialization $a_0$ with respect to the true value $a=0.01$ is $10\%$ in the case of $N=10$ and decreases even further when N increases. Moreover, the number of iterations needed to converge to a final estimate ${\bar{a}}$ is small.

7 Conclusion

In this paper, we have compared the stability regions of utility-optimizing power allocations in large distribution networks for the Distflow and the Linearized Distflow models. We have expressed the stability regions in terms of their maximum feasible arrival rates. As expected, the Linearized Distflow power flow model allows for too optimistic arrival rates. However, we are able to quantify the difference between the maximal feasible arrival rates under the Distflow model and the Linearized Distflow model. This allows us to use the computationally favorable Linearized Distflow model over the Distflow model and adjust according to the quantified difference. Even in the case of non-equal arrival rates, this could serve as a conservative approach to get insights in maximal feasible arrival rates under the Distflow model.

Appendices

Appendix A Equivalence for the Distflow model

In Sect. 2.3, we derive ways to establish if a given power allocation is feasible under both power flow models. For the Linearized Distflow model, we can directly use a recursion to compute the voltages at all nodes. However, for the Distflow model this is not possible. Therefore, in Appendix A, we find equivalence (2.14) between the voltages under the Distflow model at the root node and the node at the end of the line. To prove the equivalence in Proposition 2.1, we need Lemmas 2.1 and 2.2.

1.1 Proof of Proposition 2.1

Proof

In order to prove the implication from left to right, we take the negation of the right-hand side of (2.14) and show that this implies the negation of the left-hand side of (2.14). Suppose $v_N({V'_0}^D)>c$. By Lemma 2.2, $v_N$ is an increasing function in ${V'_0}^D$, so $c<v_N((1-\Delta )c)\le v_N(x)$ for $x\ge (1-\Delta )c$. Hence, there exists no $x\ge (1-\Delta )c$ such that $v_N(x)=c$.

For the other implication, we first observe that $v_N$ is a continuous function, because it is a composition of continuous functions itself. To prove the implication from the right-hand side of (2.14) to the left-hand side of (2.14), we assume that at ${V'_0}^D=(1-\Delta )c$, we have $v_N({V'_0}^D)\le c$. By Lemma 2.2, we know $v_N((1-\Delta )c)\le c$. Then, by the intermediate value theorem, we know that there exists $(1-\Delta )c\le x\le c$ such that $v_N(x)=c$. $\square$

1.2 Proof of Lemma 2.1

Proof

Both identities follow from the initialization in (2.12) and the definition of the sequence in (2.13).

• 1. We have from (2.13) for $j=1,2,\ldots ,N-1$,

  $$\begin{aligned} {V_{j+1}}^D-{V_j}^D = {V_j}^D-{V_{j-1}}^D+\frac{rp_j}{{V_j}^D}. \end{aligned}$$

  (A.1)

  We have from the identity in (A.1) by summation that

  $$\begin{aligned} {V_{i+1}}^D-{V_i}^D&= ({V_1}^D-{V_0}^D) + \sum _{j=1}^i \frac{rp_j}{{V_j}^D}\nonumber \\&= \frac{rp_0}{{V_0}^D} + \sum _{j=1}^i \frac{rp_j}{{V_j}^D}\nonumber \\&= \sum _{j=0}^i \frac{rp_j}{{V_j}^D},\quad i = 0,\ldots ,N-1. \end{aligned}$$

  (A.2)

• 2. We have from the identity in (A.2) by summation that

  $$\begin{aligned} {V_{n}}^D-{V_0}^D = \sum _{i=0}^{n-1}\sum _{j=0}^i \frac{rp_j}{{V_j}^D}, \quad n = 0,\ldots ,N. \end{aligned}$$

  Hence, with the initialization in (2.12) statement (2) follows immediately.

$\square$

1.3 Proof of Lemma 2.2

Proof

This is true for $j=0$, since $\frac{dv_j}{d{V_0}^D} = \frac{d{V_0}^D}{d{V_0}^D}=1.$ Suppose that there is some $0\le J\le n-1$ such that $0\le \frac{dv_j}{d{V_0}^D}\le 1$ for all $0\le j\le J$, then by Lemma 2.1 and the definition of $v_j$,

$$\begin{aligned} \frac{dv_{J+1}}{d{V_0}^D}&= \frac{d}{d{V_0}^D}\left( {V_0}^D + \sum _{i=0}^J\sum _{j=0}^i \frac{rp_j}{v_j({V_0}^D)} \right) \nonumber \\&= 1-\sum _{i=0}^J\sum _{j=0}^i \frac{rp_j}{v_j^2({V_0}^D)}\frac{dv_j}{d{V_0}^D}\nonumber \\&= 1-\sum _{i=0}^J\sum _{j=0}^i \frac{rp_j}{{({V_j}^D})^2}\frac{dv_j}{d{V_0}^D} . \end{aligned}$$

(A.3)

First, we show that $\frac{dv_{J+1}}{d{V_0}^D}\ge 0$. Recall that we establish in Lemma 2.3 that the sequence of voltages ${V_j}^D, j=0,\ldots ,N$, is an increasing sequence. This implies that $-\frac{1}{V_0^D}\le -\frac{1}{V_j^D}$ for all $j=0,\ldots ,N$. Using the fact that $-\frac{1}{V_0^D}\le -\frac{1}{V_j^D}$ for all $j=0,\ldots ,N$ and $\frac{dv_j}{d{V_0}^D}\le 1$ yields,

$$\begin{aligned} \frac{dv_{J+1}}{d{V_0}^D}&\ge 1-\frac{1}{{V_0}^D}\sum _{i=0}^J\sum _{j=0}^i\frac{rp_j}{{V_j}^D}. \end{aligned}$$

By (2.16) in Lemma 2.1, we have consequently

$$\begin{aligned} \frac{dv_{J+1}}{d{V_0}^D}&\ge 1-\frac{1}{{V_0}^D}\left( {V_{J+1}}^D-{V_0}^D\right) \\&= 2-\frac{{V_{J+1}}^D}{{V_0}^D}. \end{aligned}$$

Combining the fact that the ${V_j}^D$ is an increasing sequence and the assumption that ${V_N}^D\le 2{V_0}^D$, gives the inequality $-{V_{J+1}}^D \ge -2{V_0}^D$. Hence,

$$\begin{aligned} \frac{dv_{J+1}}{d{V_0}^D}&\ge 0. \end{aligned}$$

Now, we show that $\frac{dv_{J+1}}{d{V_0}^D}\le 1$. By the induction hypothesis, we have $-\frac{dv_j}{d{V_0}^D}\le 0$ for all $j=0,\ldots ,N$. Starting from Eq. (A.3), we immediately get,

$$\begin{aligned} \frac{dv_{J+1}}{d{V_0}^D}&= 1-\sum _{i=0}^J\sum _{j=0}^i \frac{rp_j}{{({V_j}^D})^2}\frac{dv_j}{d{V_0}^D} \le 1. \end{aligned}$$

$\square$

1.4 Proof of Lemma 2.3

Proof

We give a proof by induction on j. Recall that $V_0^D=1$. Then, from (2.17), we have $V_1^D = 1+rp_0\ge 1 = V_0^D$, and for $j=1$, we have from (2.18)

$$\begin{aligned} V_2^D&= 2V_1^D-V_0^D+\frac{rp_1}{V_1^D} = (1+rp_0)+rp_0+\frac{rp_1}{1+rp_0} \ge 1+rp_0 = V_1^D. \end{aligned}$$

Suppose that $V_j^D \ge V_{j-1}^D$ for some $j\in \{2,\ldots ,n-1\}$. Then, by (2.18) and the induction hypothesis we have

$$\begin{aligned} V_{j+1}^D&\ge V_j^D + V_{j-1}^D-V_{j-1}^D + \frac{rp_j}{V_j^D} = V_j^D + \frac{rp_j}{V_j^D} = V_j^D\left( 1+\frac{rp_j}{{V_j^D}^2}\right) \ge V_j^D \end{aligned}$$

since $r\ge 0,p_j\ge 0$ and $V_j^D\ge V_0^D=1$. $\square$

Appendix B Comparison of stability regions

In this appendix, we prove the result to compare the stability regions under both power flow models.

The critical arrival rates for the Distflow and the Linearized Distflow models are given in (3.3) and (3.2), respectively. To compare these critical rates, we use

Proof

We let

$$\begin{aligned} P(\Delta ) = 2Q^2(\Delta ), \end{aligned}$$

where

$$\begin{aligned} Q(\Delta ) := \frac{1-\Delta }{\sqrt{1-(1-\Delta )^2}}\int _0^{\sqrt{\ln \left( \frac{1}{1-\Delta } \right) }}\\ \exp (u^2)du,\ 0< \Delta \le \frac{1}{2}. \end{aligned}$$

Set $y=\sqrt{\ln \left( \frac{1}{1-\Delta } \right) }\in [0,\infty )$. Then,

$$\begin{aligned} Q(\Delta ) = \frac{\exp (-y^2)}{\sqrt{1-\exp (-2y^2)}}\int _0^y \exp (u^2)du = \frac{1}{\sqrt{\exp (2y^2)-1}}\int _0^y \exp (u^2)du =: H(y). \end{aligned}$$

Since y is a strictly increasing function of $\Delta \in [0,\frac{1}{2}]$, it is sufficient to show that H(y) is a strictly decreasing function of y. We compute

$$\begin{aligned} H'(y)&= (\exp (2y^2)-1)^{-\frac{1}{2}}\left( -\frac{2y\exp (2y^2)}{\exp (2y^2)-1}\int _0^y \exp (u^2) du+\exp (y^2) \right) . \end{aligned}$$

Therefore,

$$\begin{aligned} H'(y)<0 \iff \int _0^y \exp (u^2) > \frac{\exp (y^2)}{\frac{2y\exp (2y^2)}{\exp (2y^2)-1}} = \frac{\sinh (y^2)}{y}. \end{aligned}$$

We have,

$$\begin{aligned} \int _0^y \exp (u^2)du\Big |_{y=0} = 0 = \frac{\sinh (y^2)}{y}\Big |_{y=0}, \end{aligned}$$

so it is enough to show that

$$\begin{aligned} \frac{d}{dy}\left( \int _0^y \exp (u^2)du \right) = \exp (y^2)> \frac{d}{dy}\left( \frac{\sinh (y^2)}{y} \right) ,\ y>0. \end{aligned}$$

We compute

$$\begin{aligned} \frac{d}{dy}\left( \frac{\sinh (y^2)}{y} \right)&= \frac{2y\cosh (y^2)y-\sinh (y^2)}{y^2} \\&= 2\cosh (y^2)-\frac{\sinh (y^2)}{y^2} = \exp (y^2)+\exp (-y^2)-\frac{\sinh (y^2)}{y^2}, \end{aligned}$$

, and so, it is sufficient to show

$$\begin{aligned} \exp (-y^2)-\frac{\sinh (y^2)}{y^2}<0, y>0. \end{aligned}$$

We have,

$$\begin{aligned} \exp (-y^2)-\frac{\sinh (y^2)}{y^2}&= -\exp (-y^2)\frac{\exp (2y^2)-1-2y^2}{2y^2} < 0, \end{aligned}$$

since $\exp (t)-1-t>0$ for $t>0$. Furthermore, $P(0) = 1$ since $H(y)\rightarrow \frac{1}{\sqrt{2}}$ as $y\rightarrow 0$ and $P(\frac{1}{2})=\frac{\pi }{6}\text {erfi}\left( \sqrt{\ln (2)} \right) ^2\approx 0.77$. $\square$

Appendix C Stability result

Our paper uses the stability results derived in [26, 27]. In Appendix C, we consider the general model framework in [26] and show that we can drop one of their assumptions and still apply their theorem indicating stability, which we use in Sects. 2.3.1 and 2.3.2.

Assume that we have an arbitrary network of queues such that individual instantaneous service rates may depend on the state of the entire system. The network state is represented as a set ${\textbf{X}}(t)=(X_i(t),i\in N)$ of the queue lengths $X_i$ at the network nodes $i\in N$ at time t. Arrivals into queue i occur according to a Poisson process with a constant rate $\lambda _i$, independent of all the other processes. For this section, denote the vector of arrival rates by $\varvec{\lambda } = (\lambda _1,\ldots ,\lambda _N)$. The instantaneous service rate $\mu _i$ of a node represents the intensity of the Poisson process modeling departures (service completions) of the node. In this case, the service rate allocation algorithm $\varvec{\psi }({\textbf{X}}(t))$, mapping a network state ${\textbf{X}}(t)$ into a set of instantaneous service rates $\varvec{\mu } = (\mu _1,\ldots ,\mu _N)$, is all that is needed to specify the service allocation algorithm.

For simplicity, in this section, we drop the dependence on time t from the notation and we assume that

$$\begin{aligned} \varvec{\psi }({\varvec{x}})\in \mathop {\mathrm{arg~max}}\limits _{\varvec{\mu }\in {\mathcal {C}}}\sum _{i}g(x_i)h(\mu _i), \end{aligned}$$

(C.1)

where the set ${\mathcal {C}}$ is compact. We would like to stress that this set ${\mathcal {C}}$ can be any constraint set (as long as it is compact), so not necessarily the one that we have used in our work. We impose, in addition, the following conditions:

• Condition (H): the function $h:[0,\infty )\rightarrow {\mathbb {R}}$ is strictly increasing, differentiable, and concave (both the cases $\lim _{y\downarrow 0}h(y)=h(0)>-\infty$ and $\lim _{y\downarrow 0}h(y)=-\infty$ are allowed); and

• Condition (G): the function $g:{\mathbb {Z}}_{+}\rightarrow [0,\infty )$ is strictly increasing and such that

  $$\begin{aligned} \frac{g(y)}{\delta (y)}\rightarrow \infty \ \text {as}\ y\rightarrow \infty , \end{aligned}$$

  (C.2)

  where $\delta (y)=g(y+1)-g(y)$. Note that (C.2) is equivalent to

  $$\begin{aligned} \frac{g(y+1)}{g(y)}\rightarrow 1,\ \text {as}\ y\rightarrow \infty . \end{aligned}$$

We are going to assume that

$$\begin{aligned} \text {there exists}\ \varvec{\nu }\in {\mathcal {C}}\ \text {such that}\ \varvec{\lambda }<\varvec{\nu }, \end{aligned}$$

(C.3)

where ${\mathcal {C}}$ can be any constraint set. Note that stability condition (C.3) is equivalent to

$$\begin{aligned} \text {there exists}\ \varvec{\nu }\in {\mathcal {C}}_1\ \text {such that}\ \varvec{\lambda }<\varvec{\nu }, \end{aligned}$$

(C.4)

where

$$\begin{aligned} {\mathcal {C}}_1 = \{\varvec{\nu }: \exists \varvec{\mu }\in {\mathcal {C}}\ \text {s.t.}\ \varvec{\nu }\le \varvec{\mu }\}. \end{aligned}$$

Indeed, for the implication from (C.4) to (C.3), if there exists $\varvec{\nu }\in {\mathcal {C}}_1$ such that $\varvec{\lambda }<\varvec{\nu }$, then $\varvec{\lambda }<\varvec{\nu }\le \varvec{\mu }$ for $\varvec{\mu }\in {\mathcal {C}}$, so there exists $\varvec{\mu }\in {\mathcal {C}}$ such that $\varvec{\lambda }<\varvec{\mu }$. For the implication from (C.3) to (C.4), assume there exists $\varvec{\nu }\in {\mathcal {C}}$ such that $\varvec{\lambda }<\varvec{\nu }=\varvec{\nu }$; hence, $\varvec{\nu }\in {\mathcal {C}}_1$ such that $\varvec{\lambda }<\varvec{\nu }$.

Note that ${\mathcal {C}}_1$ is coordinate-convex by construction. Note also that if

$$\begin{aligned} \sum _i g(x_i)h(\psi _i({\varvec{x}}))\ge \sum _i g(x_i)h(\mu _i) \end{aligned}$$

for all $\varvec{\mu }\in {\mathcal {C}}$, then

$$\begin{aligned} \sum _i g(x_i)h(\psi _i({\varvec{x}}))\ge \sum _i g(x_i)h(\nu _i) \end{aligned}$$

for all $\varvec{\nu }\in {\mathcal {C}}_1$, since the function h is strictly increasing and g is non-negative. Hence,

$$\begin{aligned} \varvec{\psi }({\varvec{x}})\in \mathop {\mathrm{arg~max}}\limits _{\varvec{\nu }\in {\mathcal {C}}_1}\sum _i g(x_i)h(\nu _i). \end{aligned}$$

Then, the rest of the proof of [26, Theorem 11] follows.

Appendix D Convergence of ${\overline{V}}_N^D(t)$ to V(t) as $N\rightarrow \infty$

An important step in the proof of Theorem 3.2 is the approximation of the voltage $V_N^D$ by a continuous counterpart. In this section, we prove the convergence of a scaled version of the voltage $V_N^D$, i.e., ${\overline{V}}_N^D(t)$, as introduced in Equation (5.9), to the continuous function V(t), as introduced in Equation (5.11), as the network size N grows to infinity (cf. Proposition 5.1) and present its numerical validation.

1.1 Proof of Proposition 5.1

The proof consists of several steps that we discuss globally in Sect. 5.3 already, and for which we now present the details.

In Sect. D.1.1, we rewrite the scaled version of the voltage ${\overline{V}}_N^D(t)$, in the same form as the continuous function V(t). We show that the continuous extension $V_N^D(t)$ [as in Eq. (5.8)] of the discrete voltages $V_n^D$ is such that $V_n^D = V_N^D(n)$ for all $n=0,1,\ldots ,N$. Then, we scale the continuous extension $V_N^D(t)$ by N, such that ${\overline{V}}_N^D(t) = V_N^D(Nt)$ [as in Eq. (5.9)] and show that this scaled version ${\overline{V}}_N^D(t)$ can be written as ${\overline{V}}_N^D(t) = 1+a\int _0^t\int _0^s \frac{1}{{\overline{V}}_N^D(u)}du\ ds + {\overline{R}}_N(t)$ [as in Eq. (5.10)].

Then, for convenience, we introduce the map H (cf. Definition 2) that allows us to write,

$$\begin{aligned} {\overline{V}}_N^D(t)-V(t) = {\overline{R}}_N(t) + (H{\overline{V}}_N^D)(t)-(HV)(t). \end{aligned}$$

(D.1)

The goal is to show that the supremum over the interval [0, 1] of the absolute difference between ${\overline{V}}_N^D(t)$ and V(t) goes to zero as the size of the network N goes to infinity. Having (D.1) in mind, we consider the remainder term ${\overline{R}}_N(t)$ and the term $(H{\overline{V}}_N^D)(t)-(HV)(t)$ of the right-hand side of (D.1) separately.

In Sect. D.1.2, we show that this remainder term ${\overline{R}}_N(t)$ vanishes as $N\rightarrow \infty$ for every t in the interval [0, 1] and in Sect. D.1.3, we show that the map H is a contractive mapping on the interval [0, h], where $h<1$. This means,

$$\begin{aligned} \sup _{0\le t\le h} |(H{\overline{V}}_n^D)(t)-(HV)(t)|\le \kappa (h)\sup _{0\le t\le h} |V(t)-{\overline{V}}_N^D(t)| \end{aligned}$$

for $\kappa (h)<1$. Combining the results in Sects. D.1.2 and D.1.3, we then find that

$$\begin{aligned} \sup _{0\le t\le h} |{\overline{V}}_N^D(t)-V(t)| \rightarrow 0\ \text {as}\ N\rightarrow \infty . \end{aligned}$$

In Sect. D.1.4, we show by induction that the convergence of the supremum of the absolute difference between ${\overline{V}}_N^D(t)$ and V(t) on the interval [0, h] can be extended to the interval [0, 1], as desired in Proposition 5.1.

1.1.1 Derivation of ${\overline{V}}_N^D(t)$ as in Eq. (5.10)

In this section, we first derive an expression for the voltages under the Distflow model as a double summation, instead of the voltages given by Eqs. (5.12) and (5.13). Then, we show, in Lemma D.1, that the function $V_N^D(t)$ (as in Equation (5.9)) is a continuous extension of the discrete voltages $V_n^D$, using the alternative expression for the voltages from Lemma 2.1.

For convenience, we state the alternative expression for the voltages under the Distflow model from Lemma 2.1; i.e., for $n=0,\ldots ,N$,

$$\begin{aligned} V_{n}^D = 1+k_N\sum _{i=0}^{n-1}\sum _{j=0}^i \frac{1}{V_j^D}, \end{aligned}$$

(D.2)

with the convention $\sum _{i=0}^{-1}\cdot = 0$ and $k_N = r\lambda _N^D$. Again, observe that the alternative expression in (D.2) is only defined for integer values. In the following lemma, we define a continuous extension of the alternative expression for the voltages under the Distflow model.

Lemma D.1

(Continuous extension) The sequence $V_n^D, n=0,1,\ldots ,N$ defined in (D.2) can be represented by

$$\begin{aligned} V_N^D(t)&:= 1+k_N \int _0^{\lfloor t\rfloor }\int _0^{\lceil s\rceil }\frac{1}{V_N^D(u)}du\ ds,\ t\in [0,N]; \end{aligned}$$

(D.3)

i.e., for all integer values $n=\lfloor t\rfloor$, $V_N^D(n) = V_n^D$.

Proof

Let $t\in [0,N]$. Due to the definition of $V_N^D(t)$ in (D.3), we have

$$\begin{aligned} V_N^D(t)&= 1+k_N \int _0^{\lfloor t\rfloor }\int _0^{\lceil s\rceil }\frac{1}{V_N^D(u)}du\ ds = V_N^D(\lfloor t\rfloor ). \end{aligned}$$

(D.4)

Thus, for all integer values $n = \lfloor t\rfloor$, we have the equality

$$\begin{aligned} V_N^D(n)&= 1+k_N \int _0^{n}\int _0^{\lceil s\rceil }\frac{1}{V_N^D(u)}du\ ds. \end{aligned}$$

We show $V_N^D(n)=V_n^D$ for all $n=0,\ldots ,N$ by strong induction. Note that, for $n=0$, we have

$$\begin{aligned} V_N^D(0) = 1 = V_0^D \end{aligned}$$

by definition. Now, assume that,

$$\begin{aligned} V_N^D(n)=1+_N\int _0^n\int _0^{\lceil s\rceil }\frac{1}{V_N^D(u)}du\ ds = V_n^D, \end{aligned}$$

(D.5)

for $0\le n\le m$. Then, we need to show that $V_N^D(m+1)=V_{m+1}^D$. By definition and linearity of the integral, we have

$$\begin{aligned} V_N^D(m+1)&= 1+k_N\int _0^{m+1}\int _0^{\lceil s\rceil }\frac{1}{V_N^D(u)}du\ ds \\&= 1+k_N\int _0^{m}\int _0^{\lceil s\rceil }\frac{1}{V_N^D(u)}du\ ds + k_N\int _m^{m+1}\int _0^{\lceil s\rceil }\frac{1}{V_N^D(u)}du\ ds. \end{aligned}$$

By the induction hypothesis in (D.5), we have

$$\begin{aligned} V_N^D(m+1) = V_m^D + k_N\int _m^{m+1}\int _0^{\lceil s\rceil }\frac{1}{V_N^D(u)}du\ ds. \end{aligned}$$

(D.6)

For $s\in [m,m+1]$, the inner integral in (D.6) simplifies to

$$\begin{aligned} \int _0^{m+1}\frac{1}{V_N^D(u)}du&= \int _0^1 \frac{1}{V_N^D(u)}du + \int _1^2 \frac{1}{V_N^D(u)}du + \cdots + \int _{m}^{m+1} \frac{1}{V_N^D(u)}du. \end{aligned}$$

Now, again by the induction hypothesis in (D.5) and (D.4), we have

$$\begin{aligned} \int _0^{m+1}\frac{1}{V_N^D(u)}du = \frac{1}{V_0^D} +\frac{1}{V_1^D}+\cdots +\frac{1}{V_m^D} = \sum _{j=0}^m \frac{1}{V_j^D}. \end{aligned}$$

(D.7)

Continuing from (D.6), and using the result in (D.7), leads to

$$\begin{aligned} V_N^D(m+1)&= V_m^D + k_N\int _m^{m+1} \sum _{j=0}^m \frac{1}{V_j^D} ds \\&= V_m^D + k_N\sum _{j=0}^m \frac{1}{V_j^D}. \end{aligned}$$

Here, we recognize Eq. (2.15) from Lemma 2.1. Hence,

$$\begin{aligned} V_N^D(m+1) = V_{m+1}^D. \end{aligned}$$

This concludes the proof. $\square$

So far, we have shown that the function ${V}_N^D(t)$ is a continuous extension of the discrete voltages $V_n^D$. Now, we scale the continuous extension $V_N^D(t)$, by ${\overline{V}}_N^D(t) = V_N^D(Nt)$ (cf. Definition 1), into an integral equation that, as we show later, differs from the integral equation in (5.11) by terms that vanish as $N\rightarrow \infty$. We have, for $t\in [0,1]$ and $k_N = \frac{a}{N^2}$, that

$$\begin{aligned} {\overline{V}}_N^D(t+\frac{1}{N})&= 1+\frac{a}{N^2}\int _0^{\lfloor Nt\rfloor +1}\int _0^{\lceil s\rceil }\frac{1}{V_N^D(u)}du\ ds\nonumber \\&= 1+\frac{a}{N^2}\int _0^{\lfloor Nt\rfloor +1}\int _0^{\lceil s\rceil }\frac{1}{{\overline{V}}_N^D\left( \frac{u}{N}\right) }du\ ds. \end{aligned}$$

(D.8)

Substituting $s=Nx$ and $u=Ny$ in the integral in (D.8), we get

$$\begin{aligned} {\overline{V}}_N^D(t+\frac{1}{N})&=1+\frac{a}{N^2} N^2\int _{0}^{\frac{\lfloor Nt\rfloor +1}{N}}\int _0^{\frac{\lceil Nx\rceil }{N}}\frac{1}{{\overline{V}}_N^D(y)}dy\ dx\nonumber \\&= 1+a\int _{0}^{t}\int _0^{\frac{\lceil Nx\rceil }{N}}\frac{1}{{\overline{V}}_N^D(y)}dy\ dx + a\int _{t}^{\frac{\lfloor Nt\rfloor +1}{N}}\int _0^{\frac{\lceil Nx\rceil }{N}}\frac{1}{{\overline{V}}_N^D(y)}dy\ dx. \end{aligned}$$

(D.9)

Finally, we rewrite (D.9) as

$$\begin{aligned} {\overline{V}}_N^D(t+\frac{1}{N})&= 1+a\int _{0}^{t}\left( \int _0^{x}\frac{1}{{\overline{V}}_N^D(y)}dy + \int _x^{\frac{\lceil Nx\rceil }{N}}\frac{1}{{\overline{V}}_N^D(y)}dy \right) dx +\nonumber \\&\quad +\, a\int _{t}^{\frac{\lfloor Nt\rfloor +1}{N}}\int _0^{\frac{\lceil Nx\rceil }{N}}\frac{1}{{\overline{V}}_N^D(y)}dy\ dx \nonumber \\&= 1+a\int _{0}^{t}\int _0^{x}\frac{1}{{\overline{V}}_N^D(y)}dy\ dx + \int _0^t\int _x^{\frac{\lceil Nx\rceil }{N}}\frac{1}{{\overline{V}}_N^D(y)}dy\ dx + \nonumber \\&\quad +\, a\int _{t}^{\frac{\lfloor Nt\rfloor +1}{N}}\int _0^{\frac{\lceil Nx\rceil }{N}}\frac{1}{{\overline{V}}_N^D(y)}dy\ dx \nonumber \\&= 1+a\int _{0}^{t}\int _0^{x}\frac{1}{{\overline{V}}_N^D(y)}dy\ dx + R_N(t), \end{aligned}$$

(D.10)

where

$$\begin{aligned} R_N(t):= \int _0^t\int _x^{\frac{\lceil Nx\rceil }{N}}\frac{1}{{\overline{V}}_N^D(y)}dy\ dx + a\int _{t}^{\frac{\lfloor Nt\rfloor +1}{N}}\int _0^{\frac{\lceil Nx\rceil }{N}}\frac{1}{{\overline{V}}_N^D(y)}dy\ dx. \end{aligned}$$

(D.11)

By adding ${\overline{V}}_N^D(t)$ and subtracting ${\overline{V}}_N^D(t+\frac{1}{N})$ on both sides of (D.10), we get

$$\begin{aligned} {\overline{V}}_N^D(t)&= 1+a\int _0^t\int _0^x \frac{1}{{\overline{V}}_N^D(y)}dy\ dx + {\overline{V}}_N^D(t)-{\overline{V}}_N^D(t+\frac{1}{N})+ R_N(t). \end{aligned}$$

(D.12)

Observe that the representation in (D.12) is equal to the one in (5.10), as desired, if we define

$$\begin{aligned} {\overline{R}}_N(t) := {\overline{V}}_N^D(t)-{\overline{V}}_N^D(t+\frac{1}{N})+R_N(t). \end{aligned}$$

(D.13)

For convenience, we introduce an operator H. This is useful for the proofs in Sect. D.1.3.

Definition 2

Let $H:{\textbf{D}}_{\ge 1}[0,1]\rightarrow {\textbf{D}}_{\ge 1}[0,1]$ such that

$$\begin{aligned} (Hf)(t) = 1+a\int _0^t\int _0^s\frac{1}{f(u)}du\ ds,\ f\in {\textbf{D}}_{\ge 1}[0,1]. \end{aligned}$$

This allows us to write the difference between ${\overline{V}}_N^D(t)$ and V(t) as,

$$\begin{aligned} {\overline{V}}_N^D(t)-V(t) = {\overline{R}}_N(t) + (H{\overline{V}}_N^D)(t)-(HV)(t). \end{aligned}$$

So far, we reached Eq. (D.1) and made all the necessary preparation work to consider the remainder term ${\overline{R}}_N(t)$ and the term $(H{\overline{V}}_N^D)(t)-(HV)(t)$ of the right-hand side of (D.1) in the next two sections, separately.

1.1.2 Convergence of supremum of remainder term to zero

In this section, we show that the remainder term ${\overline{R}}_N(t)$ in (D.13) vanishes as $N\rightarrow \infty$, i.e., we show convergence of to zero as $N\rightarrow \infty$. Therefore, by (D.13), we show the convergence of as $N\rightarrow \infty$ in Lemma D.2 and as $N\rightarrow \infty$ in Lemma D.3. Then, as an immediate consequence of Lemmas D.2 and D.3, we show as $N\rightarrow \infty$.

We begin with the convergence of as $N\rightarrow \infty$ in Lemma D.2.

Lemma D.2

Let ${\overline{V}}_N^D(t)$ as defined in (5.10), then

$$\begin{aligned} \sup _{t\in [0,1]} \left| {\overline{V}}^D_N(t)-{\overline{V}}_N^D(t+\frac{1}{N})\right| \rightarrow 0\ \text {as}\ N\rightarrow \infty . \end{aligned}$$

Proof

By definition of (5.10), we have

$$\begin{aligned} {\overline{V}}_N^D(t)-{\overline{V}}_N^D(t+\frac{1}{N}) = a\bigg (\int _0^{\frac{\lfloor Nt-1\rfloor +1}{N}}\int _0^{\frac{\lceil Nx\rceil }{N}}\frac{1}{{\overline{V}}_N^D(y)}dy\ dx-\nonumber \\ -\int _0^{\frac{\lfloor Nt\rfloor +1}{N}}\int _0^{\frac{\lceil Nx\rceil }{N}}\frac{1}{{\overline{V}}_N^D(y)}dy\ dx \bigg ). \end{aligned}$$

(D.14)

However, (D.14) simplifies to,

$$\begin{aligned} {\overline{V}}_N^D(t)-{\overline{V}}_N^D(t+\frac{1}{N}) = a\left( -\int _{\frac{\lfloor Nt\rfloor }{N}}^{\frac{\lfloor Nt\rfloor +1}{N}}\int _0^{\frac{\lceil Nx\rceil }{N}} \frac{1}{{\overline{V}}_N^D(y)}dy\ dx \right) . \end{aligned}$$

(D.15)

In what follows, we bound the absolute value of the right-hand side of (D.15) in terms of N. Since ${\overline{V}}_N^D(t)\in {\textbf{D}}_{\ge 1}[0,1]$, we have, for all $t\in [0,1]$ that ${\overline{V}}_N^D(t)\ge 1$. Applying this bound yields,

$$\begin{aligned} |{\overline{V}}_N^D(t)-{\overline{V}}_N^D(t+\frac{1}{N})|&\le a\left( \int _{\frac{\lfloor Nt\rfloor }{N}}^{\frac{\lfloor Nt\rfloor +1}{N}}\int _0^{\frac{\lceil Nx\rceil }{N}} 1 dy\ dx \right) \\&= a\left( \int _{\frac{\lfloor Nt\rfloor }{N}}^{\frac{\lfloor Nt\rfloor +1}{N}} \frac{\lfloor Nt\rfloor +1}{N}dx \right) . \end{aligned}$$

Observe that $\frac{\lfloor Nt\rfloor +1}{N}-\frac{\lfloor Nt\rfloor }{N}=\frac{1}{N}$. Thus, for $t\in [0,1]$, we get

$$\begin{aligned} |{\overline{V}}_N^D(t)-{\overline{V}}_N^D(t+\frac{1}{N})|&\le a\left( \frac{N+1}{N^2} \right) . \end{aligned}$$

As a consequence, we get the desired result,

$$\begin{aligned} \sup _{t\in [0,1]} |{\overline{V}}_N^D(t)-{\overline{V}}_N^D(t+\frac{1}{N})| \rightarrow 0\ \text {as}\ N\rightarrow \infty . \end{aligned}$$

$\square$

Now, it remains to be shown that as $N\rightarrow \infty$. This is done in Lemma D.3.

Lemma D.3

Let $R_N(t)$ be defined as in (D.11), then

$$\begin{aligned} \sup _{t\in [0,1]} |R_N(t)| \rightarrow 0\ \text {as}\ N\rightarrow \infty . \end{aligned}$$

Proof

Since ${\overline{V}}_N^D(t)\in {\textbf{D}}_{\ge 1}[0,1]$, we have, for all $t\in [0,1]$ that ${\overline{V}}_N^D(t)\ge 1$, we get from (D.11), that

$$\begin{aligned} |R_N(t)|&\le a\int _0^t\int _x^{\frac{\lceil Nx\rceil }{N}}1 dy\ dx + a\int _t^{\frac{\lfloor Nt\rfloor +1}{N}}\int _0^{\frac{\lceil Nx\rceil }{N}}1 dy\ dx \nonumber \\&= a\int _0^t\int _x^{\frac{\lceil Nx\rceil }{N}}1 dy\ dx + a\int _t^{\frac{\lfloor Nt\rfloor +1}{N}} \frac{\lceil Nx\rceil }{N}\ dx. \end{aligned}$$

(D.16)

Observe that, for the first integral in (D.16), that $\frac{\lceil Nx\rceil }{N}-x \le \frac{Nx+1}{N}-x = \frac{1}{N}$, and for the second integral, that $\frac{\lfloor Nt\rfloor +1}{N}-t\le \frac{Nt+1}{N}-t = \frac{1}{N}$. Thus, for $t\in [0,1]$,

$$\begin{aligned} |R_N(t)|&\le a\int _0^t \frac{1}{N}dx + a\int _t^{\frac{\lfloor Nt\rfloor +1}{n}}\frac{\lceil Nt\rceil }{N}dx \nonumber \\&= \frac{a}{N}+\frac{a}{N} = \frac{2a}{N}. \end{aligned}$$

Hence,

$$\begin{aligned} \sup _{t\in [0,1]}|R_N(t)| \rightarrow 0\ \text {as}\ N\rightarrow \infty . \end{aligned}$$

$\square$

The convergence of the supremum of the remainder term $|{\overline{R}}_N(t)|$ follows immediately from Lemmas D.2 and D.3 and is summarized in the following corollary.

Corollary D.1

Let ${\overline{R}}_N(t)$ be defined as in (D.13), then

$$\begin{aligned} \sup _{t\in [0,1]}|{\overline{R}}_N(t)|\rightarrow 0\ \text {as}\ N\rightarrow \infty . \end{aligned}$$

1.1.3 Contractive property of map H

In this section, we show that $H:{\textbf{D}}_{\ge 1}[0,h]\rightarrow {\textbf{D}}_{\ge 1}[0,h]$ is a contractive mapping for a suitably chosen $h<1$. This allows us to show that as $N\rightarrow \infty$.

In Lemma D.4, we show that H is a contraction on ${\textbf{D}}_{\ge 1}[0,h]$ for a suitably chosen $h<1$.

Lemma D.4

Let H as in Definition 2. There exists $0\le h<1$, such that H is a contractive mapping on ${\textbf{D}}_{\ge 1}[0,h]$, i.e.,

$$\begin{aligned} \sup _{0\le t\le h} |(H{\overline{V}}_n^D)(t)-(HV)(t)|\le \kappa (h)\sup _{0\le t\le h} |{\overline{V}}_N^D(t)-V(t)| \end{aligned}$$

with $\kappa (h)<1$.

Proof

For $0\le t\le h$, we have

$$\begin{aligned} (H{\overline{V}}_N^D)(t)-(HV)(t)&= a \int _0^t\int _0^x \frac{1}{{\overline{V}}_N^D(y)}-\frac{1}{V(y)}dy\ dx \\&= a \int _0^t\int _0^x \frac{V(y)-{\overline{V}}_N^D(y)}{{\overline{V}}_N^D(y)V(y)}dy\ dx \end{aligned}$$

Since ${\overline{V}}_N^D(t)\in {\textbf{D}}_{\ge 1}[0,1]$, we have, for all $t\in [0,1]$ that ${\overline{V}}_N^D(t)\ge 1$. Therefore, we get,

$$\begin{aligned} |(H{\overline{V}}_n^D)(t)-(HV)(t)|&\le a \int _0^t\int _0^x |V(y)-{\overline{V}}_N^D(y)|dy\ dx. \end{aligned}$$

Furthermore, by definition of the supremum, we get

$$\begin{aligned} |(H{\overline{V}}_n^D)(t)-(HV)(t)|&\le a \int _0^t\int _0^x \sup _{0\le y\le h} |V(y)-{\overline{V}}_N^D(y)| dy\ dx\\&= a\sup _{0\le y\le h} |V(y)-{\overline{V}}_N^D(y)|\int _0^t\int _0^x 1 dy\ dx. \end{aligned}$$

Thus, for $0\le t\le h$,

$$\begin{aligned} |(H{\overline{V}}_n^D)(t)-(HV)(t)|&\le \frac{ah^2}{2}\sup _{0\le y\le h} |V(y)-{\overline{V}}_N^D(y)|. \end{aligned}$$

Hence,

$$\begin{aligned} \sup _{0\le t\le h} |(H{\overline{V}}_n^D)(t)-(HV)(t)|\le \kappa (h)\sup _{0\le t\le h} |V(t)-{\overline{V}}_N^D(t)| \end{aligned}$$

where $\kappa (h)=\frac{ah^2}{2}$. If $a>2$, take $h<\sqrt{\frac{2}{a}}$. Then $\kappa (h) = \frac{ah^2}{2}<1$. Otherwise, if $a\le 2$, then $\kappa (h) = \frac{ah^2}{2} \le h^2 < 1$, since $h<1$. $\square$

From the proof of Lemma D.4, it follows that the contraction property of the map H on the interval [0, h] depends on the size of a. We choose

$$\begin{aligned} h<{\left\{ \begin{array}{ll} \sqrt{\frac{2}{a}} &{} \text {if}\ a>2, \\ 1 &{} \text {if}\ a\le 2. \end{array}\right. } \end{aligned}$$

(D.17)

Recall from (5.14) that $a = {\hat{r}}\lambda _N^DN^2\ge 0$. Thus, for every $a\ge 0$, we choose a suitable h, according to (D.17), such that H is a contraction on the interval [0, h].

1.1.4 Convergence of supremum on different intervals of ${\overline{V}}_N^D(t)$ to V(t)

In Sects. D.1.2 and D.1.3, we have shown that the remainder term ${\overline{R}}_N(t)$ vanishes as $N\rightarrow \infty$ for every t in the interval [0, 1] and that the map H is a contractive mapping on the interval [0, h], respectively. Using both results yields the convergence of as $N\rightarrow \infty$ in Lemma D.5. Then, we show, by induction, that as $N\rightarrow \infty$ for every $m=0,\ldots ,M$ such that $Mh=1$ in Lemma D.6. This allows us to show that as $N\rightarrow \infty$, as desired in Proposition 5.1.

Using the contraction property as shown in Lemma D.4, we show the convergence of ${\overline{V}}_N^D(t)$ to V(t) in the supremum norm.

Lemma D.5

Let h as in (D.17) and ${\overline{V}}_N^D(t)$ as defined in (5.10) and V(t) as in (5.11). Then, we have,

$$\begin{aligned} \sup _{0\le t\le h} |{\overline{V}}_N^D(t)-V(t)|\rightarrow 0,\ \text {as}\ N\rightarrow \infty . \end{aligned}$$

Proof

By Eq. (D.1) and the triangle inequality, we have

$$\begin{aligned} \sup _{0\le t\le h} |{\overline{V}}_N^D(t)-V(t)|&= \sup _{0\le t\le h} |{\overline{R}}_N(t)+(H{\overline{V}}_N^D)(t)-(HV)(t)| \nonumber \\&\le \sup _{0\le t\le h}|{\overline{R}}_N(t)|+\sup _{0\le t\le h}|(H{\overline{V}}_N^D)(t)-(HV)(t)|. \end{aligned}$$

(D.18)

Using Lemma D.4 in the right-hand side of (D.18), we get

$$\begin{aligned} \sup _{0\le t\le h} |{\overline{V}}_N^D(t)-V(t)|&\le \sup _{0\le t\le h}|{\overline{R}}_N(t)|+\kappa (h)\sup _{0\le t\le h}|{\overline{V}}_N^D(t)-V(t)|. \end{aligned}$$

Rearranging terms yields,

$$\begin{aligned} \sup _{0\le t\le h}|{\overline{V}}_N^D(t)-V(t)|\le \frac{1}{1-\kappa (h)}\sup _{0\le t\le h}|{\overline{R}}_N(t)|. \end{aligned}$$

Remark that $\sup _{0\le t\le h}|{\overline{R}}_N(t)| \le \sup _{0\le t\le 1}|{\overline{R}}_N(t)|$. Hence, by Corollary D.1,

$$\begin{aligned} \sup _{0\le t\le h}|{\overline{V}}_N^D(t)-V(t)|\rightarrow 0\ \text {as}\ N\rightarrow \infty . \end{aligned}$$

$\square$

Now, we are in position to prove that on every small interval $[mh,(m+1)h]$ (where $m=0,\ldots ,M$ such that $Mh=1$) the supremum over $t\in [mh,(m+1)h]$ of the absolute difference between ${\overline{V}}_N^D(t)$ and V(t) goes to zero as $N\rightarrow \infty$.

Lemma D.6

Let h as in (D.17). For every $0\le m\le M$, let $t\in [mh,(m+1)h]$. Then,

$$\begin{aligned} \sup _{mh\le t\le (m+1)h} |{\overline{V}}_N^D(t)-V(t)|\rightarrow 0,\ \text {as}\ N\rightarrow \infty . \end{aligned}$$

(D.19)

Proof

We prove the statement in (D.19) by strong induction. By Lemma D.5, the statement holds for $m=0$, i.e.,

$$\begin{aligned} \sup _{0\le t\le h}|{\overline{V}}_N^D(t)-V(t)|\rightarrow 0\ \text {as}\ N\rightarrow \infty . \end{aligned}$$

Suppose that the statement holds for $0\le t\le mh$:

$$\begin{aligned} \sup _{0\le t\le mh} |{\overline{V}}_N^D(t)-V(t)| \rightarrow 0\ \text {as}\ N\rightarrow \infty . \end{aligned}$$

(D.20)

Then, we need to show that for $mh\le t\le (m+1)h$, the statement in (D.19) holds. Before we move on with expressions of the form $\sup _{mh\le t\le (m+1)h}|{\overline{V}}_N^D(t)-V(t)|$, we first try to bound the expression $|{\overline{V}}_N^D(t)-V(t)|$ itself. By Eqs. (5.11) and (5.8), we have

$$\begin{aligned} {\overline{V}}_N^D(t)-V(t)&= a\int _0^t\int _0^x \frac{1}{{\overline{V}}_N^D(y)}-\frac{1}{V(y)}dy\ dx + {\overline{R}}_N(t) \nonumber \\&= a\int _0^{mh}\int _0^x \frac{1}{{\overline{V}}_N^D(y)}-\frac{1}{V(y)}dy\ dx +\nonumber \\&\quad + \,a\int _{mh}^{t}\int _0^x \frac{1}{{\overline{V}}_N^D(y)}-\frac{1}{V(y)}dy\ dx + {\overline{R}}_N(t), \end{aligned}$$

and since ${\overline{V}}_N^D(t)\in {\textbf{D}}_{\ge 1}[0,1]$, we have, for all $t\in [0,1]$ that ${\overline{V}}_N^D(t)\ge 1$. Therefore, we get,

$$\begin{aligned} |{\overline{V}}_N^D(t)-V(t)|&= \le a\int _0^{mh}\int _0^x |V(y)-{\overline{V}}_N^D(y)| dy\ dx + \nonumber \\&\quad +\, a\int _{mh}^{t}\int _0^x |V(y)-{\overline{V}}_N^D(y)|dy\ dx + |{\overline{R}}_N(t)|. \end{aligned}$$

(D.21)

In what follows, we find for each term in (D.21) an upper bound in either terms of $\sup _{0\le t\le mh}|{\overline{V}}_N^D-V(t)|$ or $\sup _{mh\le t\le (m+1)h}|{\overline{V}}_N^D(t)-V(t)|$ or terms that vanish as $N\rightarrow \infty$. Then, by rearranging terms and using the induction step in (D.20), we find that the statement in (D.19) holds.

Consider the first term on the right-hand side of (D.21). By definition of the supremum, we get

$$\begin{aligned} a\int _0^{mh}\int _{0}^x |V(y)-{\overline{V}}_N^D(y)|dy\ dx&\le a\int _0^{mh}\int _0^x \sup _{0\le y\le mh}|V(y)-{\overline{V}}_N^D(y)| dy\ dx \nonumber \\&= \frac{a(mh)^2}{2}\sup _{0\le y\le mh}|V(y)-{\overline{V}}_N^D(y)|. \end{aligned}$$

Now, consider the second term on the right-hand side of (D.21). This integral can be bounded as follows:

$$\begin{aligned} a\int _{mh}^{t}\int _0^x |V(y)-{\overline{V}}_N^D(y)|dy\ dx&= a\int _{mh}^t\int _0^{mh} |V(y)-{\overline{V}}_N^D(y)|dy\ dx \\&\quad + a\int _{mh}^t\int _{mh}^x |V(y)-{\overline{V}}_N^D(y)|dy\ dx \\&\le ah(mh)\sup _{0\le y\le mh}|V(y)-{\overline{V}}_N^D(y)| \\&\quad +\frac{ah^2}{2}\sup _{mh\le y\le (m+1)h}|V(y)-{\overline{V}}_N^D(y)|. \end{aligned}$$

Finally, the third term on the right-hand side of (D.21) is the remainder term $|{\overline{R}}_N(t)|$. Continuing from (D.21) and taking the supremum over the interval $m\le t\le (m+1)h$ of $|{\overline{V}}_N^D(t)-V(t)|$ on the left-hand side yield,

$$\begin{aligned} \sup _{mh\le t\le (m+1)h} |{\overline{V}}_N^D(t)-V(t)|\le & {} \frac{a(mh)^2}{2}\sup _{0\le y\le mh}|V(y)-{\overline{V}}_N^D(y)|\\{} & {} +\,amh^2\sup _{0\le y\le mh}|V(y)-{\overline{V}}_N^D(y)|\\{} & {} +\, \frac{ah^2}{2}\sup _{mh\le y\le (m+1)h}|V(y)-{\overline{V}}_N^D(y)| + |{\overline{R}}_N(t)| \end{aligned}$$

Simplifying and rearranging terms yields,

$$\begin{aligned}{} & {} \sup _{mh\le t\le (m+1)h} |{\overline{V}}_N^D(t)-V(t)|\\\le & {} \frac{1}{1-\frac{1}{2}ah^2} \bigg (amh^2\left( \frac{m}{2}+1\right) \sup _{0\le y\le mh}|V(y)-{\overline{V}}_N^D(y)| +|{\overline{R}}_N(t)| \bigg ). \end{aligned}$$

Remark that $\sup _{mh\le t\le (m+1)h}|{\overline{R}}_N(t)| \le \sup _{0\le t\le 1}|{\overline{R}}_N(t)|$. Hence, by Corollary D.1 and from the induction hypothesis that , it follows that

$$\begin{aligned} \sup _{mh\le t\le (m+1)h} |{\overline{V}}_N^D(t)-V(t)|\rightarrow 0\ \text {as}\ N\rightarrow \infty . \end{aligned}$$

$\square$

Now, we have all the ingredients to prove Theorem 5.1.

Proof of Proposition 5.1

Let $t\in [0,1]$ and ${\overline{V}}_N^D(t)$ as defined in (5.10) and V(t) as in (5.11). Furthermore, let h as in (D.17) and consider the partition $\{[0,h],[h,2h],\ldots ,[(M-1)h,Mh]\}$ of the interval [0, 1]. Then,

$$\begin{aligned} \sup _{0\le t\le 1}|{\overline{V}}_N^D(t)-V(t)|\le & {} \sup _{0\le t\le h} |{\overline{V}}_N^D(t)-V(t)| \\{} & {} +\, \sup _{h\le t\le 2h} |{\overline{V}}_N^D(t)-V(t)| + \cdots + \sup _{(M-1)h\le t\le Mh} |{\overline{V}}_N^D(t)-V(t)|. \end{aligned}$$

By Lemma D.6, we have for the supremum over each interval $[mh,(m+1)h]$, for $m=0,\ldots ,M-1$, of the absolute difference of ${\overline{V}}_N^D(t)$ and V(t) that it converges to zero as $N\rightarrow \infty$. Hence,

$$\begin{aligned} \sup _{0\le t\le 1}|{\overline{V}}_N^D(t)-V(t)|\rightarrow 0\ \text {as}\ N\rightarrow \infty . \end{aligned}$$

$\square$

1.2 Numerical validation of convergence in Proposition 5.1

We validate the convergence established in Theorem 5.1. We run the recursion $V_n^D, n=0,\ldots ,N$ until node N and evaluate the approximation V(t) in the point $t=1$ for different values of the parameter a. The results obtained via the approximation are compared with the recursion values, i.e., we compare V(1) [computed via (5.17) where the function $f_0$ is defined in Lemma E.1] with $V_N^D$ for different values of a. For each value of a, we compute V(1) and $V_N^D$ for several values of N, ranging from 10 to $10^4$. One example of V(1) and $V_N^D$, where $a=0.05$ is shown in Fig. 3 and Table 5. In the case that N is small, e.g., $N=10$, the absolute and relative errors are only $\approx 0.2\%$.

Fig. 3

Recursion values $V_N^D$ (blue solid circles) and approximation V(1) (red solid line), for $a=0.05$ (Color figure online)

Table 5 Absolute and relative error of different values of N for $a=0.05$

One can see from Fig. 3 that, as N gets bigger, $V_N^D$ decreases to V(1). Thus, the approximation V(1) offers a close approximation of the voltage $V_N^D$ under the Distflow model, especially for large N. The errors are bigger for small N.

Appendix E Integral equation

The solution to the integral equation in (5.15) and (5.16) is given in this appendix.

Lemma E.1

For $t\ge 0,k>0,y>0,w\ge 0$, the nonlinear differential equation

$$\begin{aligned} f''(t) = \frac{k}{f(t)} \end{aligned}$$

with initial conditions $f(0)=y$ and $f'(0)=w$ has the exact solution $f(t) = \gamma f_0(\alpha +\beta t)$. Here, $f_0$ is given by

$$\begin{aligned} f_0(x) = \exp (U^2(x)),\quad \textrm{for}~x\ge 0, \end{aligned}$$

where U(x), for $x\ge 0$, is given by

$$\begin{aligned} \int _0^{U(x)}\exp (u^2)~du = \frac{x}{\sqrt{2}}, \end{aligned}$$

(E.1)

and

$$\begin{aligned} \alpha&= \sqrt{2}\int _0^\frac{w}{\sqrt{2k}} \exp (u^2)~du,\\ \beta&= \frac{\sqrt{k}}{y}\exp \left( \frac{w^2}{2k}\right) ,\\ \gamma&= y\exp \left( \frac{-w^2}{2k} \right) . \end{aligned}$$

Proof of Lemma E.1

From the differential equation

$$\begin{aligned} f''(t) = 1/f(t),\ t\ge 0 \end{aligned}$$

(E.2)

we get

$$\begin{aligned} f'(u)f''(u) = f'(u)/f(u),\ 0\le u\le t. \end{aligned}$$

(E.3)

Integrating Eq. (E.3) over u from 0 to t using $f(0)=1, f'(0)=0$, we get

$$\begin{aligned} \int _0^t f'(u)f''(u)du = \frac{1}{2}(f'(t))^2 = \int _0^t \frac{kf'(u)}{f(u)} = \ln (f(t)). \end{aligned}$$

Hence, for $t>0$,

$$\begin{aligned} \frac{f'(t)}{\left( 2\ln (f(t))\right) ^{\frac{1}{2}}} = 1. \end{aligned}$$

(E.4)

Integrating $f'(u)/(2\ln (f(t)))^{\frac{1}{2}})=1$ from $u=0$ to $u=t$, while substituting $s=f(u)\in [1,f(t)]$, we get

$$\begin{aligned} \int _1^{f(t)} \frac{1}{(2\ln (s))^{\frac{1}{2}}}ds = t. \end{aligned}$$

(E.5)

Substituting $v=(\ln (s))^{\frac{1}{2}}, s=\exp (v^2),ds=2v\exp (v^2)dv$ in the integral (E.5), we get

$$\begin{aligned} \int _0^{(\ln (f(t))^{\frac{1}{2}}}\frac{2v\exp (v^2)}{\sqrt{2}v} dv = t, \end{aligned}$$

i.e.,

$$\begin{aligned} \int _0^{(\ln (f(t))^{\frac{1}{2}}}\exp (v^2)dv = t/\sqrt{2}. \end{aligned}$$

Hence, when we define U(t) by

$$\begin{aligned} \int _0^{U(t)}\exp (v^2)dv = t/\sqrt{2}, \end{aligned}$$

(E.6)

we have

$$\begin{aligned} (\ln (f(t))^{\frac{1}{2}}=U(t), \text {i.e.}, f(t)=\exp (U^2(t)). \end{aligned}$$

(E.7)

Observe that from (E.4) and (E.7) we get,

$$\begin{aligned} f'(t) = \sqrt{2}U(t), t\ge 0. \end{aligned}$$

(E.8)

We denote in the sequel the solution of $f''(t)=1/f(t), t\ge 0$, with $f(0)=1,f'(0)=0$ by $f_0(t)$. Next, for given $w\ge 0$ and $y,k>0$ we consider the following initial value problem:

$$\begin{aligned} f''(t) = \frac{k}{f(t)}~\text {for}~ t\ge 0; f(0)=y\quad \text {and}\quad f'(0)=w. \end{aligned}$$

(E.9)

We find unique, explicit and positive values $\gamma ,\alpha$ and $\beta$ such that

$$\begin{aligned} f(t) = \gamma f_0(\alpha +\beta t)~\text {for}~t\ge 0. \end{aligned}$$

(E.10)

The function f(t) in (E.10) fulfills the conditions in (E.9),

$$\begin{aligned} {\left\{ \begin{array}{ll} f''(t) = {\beta }^2{\gamma }^2f_0^{''}(\alpha +\beta t)f_0(\alpha +\beta t) \overset{(E.2)}{=} {\beta }^2{\gamma }^2 = k,\\ f(0) = \gamma f_0(\alpha ) = y,\\ f'(0) = \beta \gamma f_0^{'}(\alpha ) = w, \end{array}\right. } \end{aligned}$$

i.e., if and only if

$$\begin{aligned} {\left\{ \begin{array}{ll} \gamma \beta = \sqrt{k},\\ \gamma f_0(\alpha ) = y,\\ f_0^{'}(\alpha ) = \frac{w}{\sqrt{k}}. \end{array}\right. } \end{aligned}$$

From $f_0'(\alpha )=\frac{w}{\sqrt{2k}}$, we get $U(\alpha )=\frac{w}{\sqrt{2k}}$ by (E.8), and w, from (E.6), we find

$$\begin{aligned} \int _0^\frac{w}{\sqrt{2k}} \exp (u^2)du = \frac{\alpha }{\sqrt{2}}. \end{aligned}$$

Hence, subsequently we get

$$\begin{aligned} \alpha&= \sqrt{2}\int _0^\frac{w}{\sqrt{2k}} \exp (u^2)du, \\ \gamma&= \frac{y}{f_0(\alpha )} = \frac{y}{\exp (U(\alpha )^2)} = \frac{y}{\exp \left( \frac{w^2}{2k}\right) }, \\ \beta&= \frac{\sqrt{k}}{\frac{y}{\exp \left( \frac{w^2}{2k}\right) }} = \frac{\sqrt{k}}{y}\exp \left( \frac{w^2}{2k}\right) . \end{aligned}$$

$\square$

Notice that we do not find an elementary closed-form solution of the function f, since f is given in terms of U(x). The function U(x), for $x\ge 0$, is given by Eq. (E.1). The left-hand side of (E.1) is equal to $\frac{1}{2}\sqrt{\pi } \text {erfi}(U(x))$, where $\text {erfi}(z)$ is the imaginary error function.

Appendix F Notation

• N: number of charging stations.

• $\lambda$: arrival rate of EVs at each charging station.

• $\varvec{\lambda } = (\lambda ,\ldots ,\lambda )$: the arrival rate at each charging station.

• $\lambda _c$: critical arrival rate for EVs.

• $\lambda _c^D$: critical arrival rate under the Distflow model.

• $\lambda _c^L$: critical arrival rate under the Linearized Distflow model.

• ${\textbf{X}}(t)=(X_1(t),\ldots ,X_N(t))$: the number of EVs at each charging station at time t.

• ${\textbf{p}}=(p_1(t),\ldots ,p_N(t))$: the power allocated to each charging station at time t.

• $g(\cdot ),h(\cdot )$: auxiliary functions to define utility-maximizing mechanism.

• ${\mathcal {C}}$: N-dimensional vector space that contains the distribution network constraints.

• ${\mathcal {G}}=({\mathcal {I}},{\mathcal {E}})$: directed graph.

• ${\mathcal {I}}$: set of nodes.

• ${\mathcal {E}}$: set of edges.

• $\epsilon _{ij}$: edge.

• z: impedance on edge.

• r: resistance on edge.

• x: reactance on edge.

• ${\tilde{V}}_j$: real voltage at node $j\in {\mathcal {I}}$.

• ${\tilde{V}}_j^L$: real voltage at node $j\in {\mathcal {I}}$ under the Linearized Distflow model.

• ${\tilde{V}}_j^D$: real voltage at node $j\in {\mathcal {I}}$ under the Distflow model.

• $W_{ij}$: transformed voltage; product of real voltage at node $i,j\in {\mathcal {I}}$ (after relabeling).

• $W_{ij}^L$: transformed voltage; product of real voltage at node $i,j\in {\mathcal {I}}$ under the Linearized Distflow model (after relabeling).

• $W_{ij}^D$: transformed voltage; product of real voltage at node $i,j\in {\mathcal {I}}$ under the Distflow model (after relabeling).

• ${\tilde{s}}_i$: complex power consumption at node $i\in {\mathcal {I}}$.

• ${\tilde{p}}_i$: active power consumption at node $i\in {\mathcal {I}}$.

• ${\tilde{q}}_i$: reactive power consumption at node $i\in {\mathcal {I}}$.

• $I_{ij}$: complex current on edge $\epsilon _{ij}\in {\mathcal {E}}$.

• ${\tilde{S}}_{ij}$: complex power flowing over edge $\epsilon _{ij}\in {\mathcal {E}}$.

• ${\tilde{P}}_{ij}$: active power flowing over edge $\epsilon _{ij}\in {\mathcal {E}}$.

• ${\tilde{Q}}_{ij}$: reactive power flowing over edge $\epsilon _{ij}\in {\mathcal {E}}$.

• $\Delta$: bound on the maximal voltage drop.

• $k_N$: product of r and arrival rate $\lambda _c^D$.