A Physarum-inspired model for the probit-based stochastic user equilibrium problem

TL;DR

This paper introduces a Physarum-inspired model tailored for the probit-based stochastic user equilibrium problem in traffic networks, offering a novel approach to find equilibrium flows considering two-way traffic and normal distribution of perceived costs.

Contribution

The work modifies the Physarum model to handle two-way traffic networks and probabilistic cost perception, providing a new method for stochastic user equilibrium computation.

Findings

The proposed model effectively finds equilibrium flows in traffic networks.

It performs favorably compared to the Method of Successive Averages.

Numerical examples demonstrate its applicability and efficiency.

Abstract

Stochastic user equilibrium is an important issue in the traffic assignment problems, tradition models for the stochastic user equilibrium problem are designed as mathematical programming problems. In this article, a Physarum-inspired model for the probit-based stochastic user equilibrium problem is proposed. There are two main contributions of our work. On the one hand, the origin Physarum model is modified to find the shortest path in traffic direction networks with the properties of two-way traffic characteristic. On the other hand, the modified Physarum-inspired model could get the equilibrium flows when traveller's perceived transportation cost complies with normal distribution. The proposed method is constituted with a two-step procedure. First, the modified Physarum model is applied to get the auxiliary flows. Second, the auxiliary flows are averaged to obtain the equilibrium flows. Numerical examples are conducted to illustrate the performance of the proposed method, which is compared with the Method of Successive Average method.

Tables

Table 1. Table 1: Parameters of the link cost functions for each link, t a ​ ( x a ) = α a + β a ​ x a 4 subscript t a subscript 𝑥 𝑎 subscript 𝛼 𝑎 subscript 𝛽 𝑎 superscript subscript 𝑥 𝑎 4 {{\rm{t}}_{\rm{a}}}({x_{a}})={\alpha_{a}}+{\beta_{a}}{x_{a}}^{4}

From node	TO node	${\alpha }_a$	${\beta }_a$
1	2	20	0.0056
1	5	18	0.0078
2	1	20	0.0071
2	6	19	0.0033
2	3	23	0.0086
3	2	23	0.0108
3	7	16	0.0101
3	4	17	0.0063
4	3	17	0.0116
4	8	22	0.0138
5	1	18	0.0131
5	6	14	0.0093
5	9	24	0.0026
6	2	19	0.0048
6	5	14	0.0041
6	7	17	0.0123
6	10	20	0.0056
7	3	16	0.0078
7	6	17	0.0071
7	8	13	0.0033
7	11	26	0.0086
8	4	22	0.0108
8	7	13	0.0101
8	12	19	0.0063
9	5	24	0.0016
9	10	7	0.0138
10	9	7	0.0131
10	6	20	0.0093
10	11	18	0.0026
11	10	18	0.0048
11	7	26	0.0141
11	12	17	0.0123
12	8	19	0.0056
12	11	17	0.0078

Table 2. Table 2: The link traffic flow calculated by both algorithms in 4.1

From node	TO node	The MSA algorithm	The proposed algorithm
outer iteration $n$		12233	236
computer time(s)		2.317002	0.069317
1	2	10.3639	10.2070
1	5	9.6361	9.5445
2	1	0	0
2	6	4.4459	4.4894
2	3	5.9180	5.7079
3	2	0	0
3	7	2.7803	2.5665
3	4	3.1377	3.1324
4	3	0	0
4	8	3.1377	3.1328
5	1	0	0
5	6	4.9213	4.7524
5	9	4.7148	4.7896
6	2	0	0
6	5	0	0
6	7	5.6918	5.4874
6	10	3.6754	3.7607
7	3	0	0
7	6	0	0
7	8	7.6230	7.5404
7	11	0.8492	0.5210
8	4	0	0
8	7	0	0
8	12	10.7607	10.6752
9	5	0	0
9	10	4.7148	4.7948
10	9	0	0
10	6	0	0
10	11	8.3902	8.5612
11	10	0	0
11	7	0	0
11	12	9.2393	9.0669
12	8	0	0
12	11	0	0

Table 3. Table 3: The link traffic flow calculated by both algorithms in 4.2

From node	TO node	The MSA algorithm	The proposed algorithm
outer iteration $n$		686	179
computer time(s)		2.162422	0.068314
1	2	10.3988	10.1945
1	5	9.6058	9.4830
2	1	0	0
2	6	3.6292	3.5431
2	3	6.7686	6.6450
3	2	0	0
3	7	2.2849	2.1953
3	4	4.4803	4.4454
4	3	0	0
4	8	4.4803	4.4424
5	1	0	0
5	6	5.1153	4.7598
5	9	4.4905	4.7273
6	2	0	0
6	5	0	0
6	7	6.4263	6.3017
6	10	2.3182	1.9937
7	3	0	0
7	6	0	0
7	8	8.7109	8.4797
7	11	0	0.0325
8	4	0	0
8	7	0	0
8	12	3.2044	3.0647
9	5	0	0
9	10	4.4905	4.7273
10	9	0	0
10	6	0	0
10	11	6.8088	6.7691
11	10	0	0
11	7	0	0
11	12	6.8088	6.7691
12	8	0	0
12	11	0	0

Equations

$$c_k^{rs}=t_a\cdot {\delta }_{a,k}^{rs}$$

$$\mathbf{c}={△}^{\mathbf{T}}\cdot \mathbf{t}$$

$$T_a=t_a(x_a)+{\epsilon }_a,\forall a$$

$${\epsilon }_a\sim N(0,\gamma t_{\mathrm{a}}^0),\forall a$$

$$T_a\sim N({\mathrm{t}}_{\mathrm{a}},\gamma t_{\mathrm{a}}^0),\forall a$$

$$C_k^{rs}=\sum _a{T}_a.{\delta }_{a,k}^{rs},\forall r,s,k\in K_{rs}$$

$$c_k^{rs}=\sum _a{t}_a.{\delta }_{a,k}^{rs},\forall r,s,k\in K_{rs}$$

$$f_k^{rs}=q_{rs}.P_k^{rs},\forall r,s,k\in K_{rs}$$

$$P_k^{rs}=P(C_k^{rs}\le C_l^{rs},\forall l\ne k),\forall r,s,k\in K_{rs}$$

$$\min Z(x)=-\sum _{\mathrm{rs}}{q}_{rs}.\mathrm{E}[\underset{k\in {\mathrm{K}}_{rs}}{\min }{\{\mathrm{C}}_k^{rs}\}|{\mathbf{c}}^{rs}(\mathbf{x})]+\sum _a\{x_a.t_a(x_a)-\underset{0}{\overset{a}{\int }}{t}_a(w)dw\}$$

$$\frac{\partial }{\partial {\mathrm{c}}_{\mathrm{k}}^{rs}}E[\underset{k\in {\mathrm{K}}_{rs}}{\min }{\{\mathrm{C}}_k^{rs}\}|{\mathbf{c}}^{rs}(\mathbf{x}){]=\mathrm{P}}_k^{rs}$$

$${\mathrm{x}}_a^{(n+1)}=x_a^n+\frac{1}{n}({\hat{y}}_a^n-x_a^n),\forall a$$

$$\sqrt{\sum _{a\in A}({\mathrm{x}}_a^{(n+1)}-{\hat{y}}_a^n)}\le {\epsilon }_0$$

$$Q_{ij}=\frac{\pi {\mathrm{r}}_{ij}^4}{8\eta L_{ij}}(p_i-p_j)=\frac{D_{ij}}{L_{ij}}(p_i-p_j)$$

$$\sum Q_{ij}=0,(\mathrm{j}\ne 1,2)$$

$$\sum _i{Q}_{i1}+I_0=0$$

$$\sum _i{Q}_{i2}-I_0=0$$

$$\sum\limits_{i}{\frac{{{D_{ij}}}}{{{L_{ij}}}}({p_{i}}-{p_{j}})}=\left\{\begin{array}[]{l}-1~{}~{}~{}~{}~{}{\rm{}}for~{}{\rm{}}j=1,\\ +1~{}~{}~{}~{}~{}{\rm{}}for~{}{\rm{}}j=2,\\ 0~{}~{}~{}~{}~{}~{}~{}{\rm{}}otherwise{\rm{}}{}{}{}\\ \end{array}\right.$$

$$\frac{d}{dt}{D}_{ij}=f(|Q_{ij}|)-\alpha D_{ij}$$

$$\sum\limits_{i}{(\frac{{{D_{ij}}}}{{{L_{ij}}}}+\frac{{{D_{ji}}}}{{{L_{ji}}}})({p_{i}}-{p_{j}})}=\left\{\begin{array}[]{l}-1~{}~{}~{}~{}~{}{\rm{}}for~{}{\rm{}}j=1,\\ +1~{}~{}~{}~{}~{}{\rm{}}for~{}{\rm{}}j=2,\\ 0~{}~{}~{}~{}~{}~{}~{}{\rm{}}otherwise{\rm{}}{}{}{}\\ \end{array}\right.$$

$${Q_{ij}}=\left\{\begin{array}[]{l}\frac{{{D_{ij}}}}{{{L_{ij}}}}({p_{i}}-{p_{j}}),~{}~{}~{}~{}~{}{\rm{}}\frac{{{D_{ij}}}}{{{L_{ij}}}}({p_{i}}-{p_{j}})>0\\ 0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm{}}otherwise{\rm{}}{}{}{}\\ \end{array}\right.$$

$$\sum _i{Q}_{io}+I_o=0,o\in O$$

$$\sum _i{Q}_{id}-I_d=0,d\in D$$

$$\sum\limits_{i}{(\frac{{{D_{ij}}}}{{{L_{ij}}}}+\frac{{{D_{ji}}}}{{{L_{ji}}}})({p_{i}}-{p_{j}})}=\left\{\begin{array}[]{l}-{I_{o}},~{}~{}~{}~{}~{}{\rm{}}\forall o\in O,\\ +{I_{d}},~{}~{}~{}~{}~{}{\rm{}}\forall d\in D,\\ 0,~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm{}}otherwise\\ \end{array}\right.$$

$$C_a^n=\frac{C_a^{n-1}+t_a(x_a)}{2}$$

$${\mathrm{t}}_{\mathrm{a}}(x_a)={\alpha }_a+{\beta }_a{x_a}^4$$

Full text

1 Introduction

The traffic assignment problem (TAP) refers to assign traffic trip of each origin-destination (OD) pair to the links in the transportation networks and give an OD trip matrix [1, 2, 3, 4, 5, 6]. Traditionally, the TAP problem falls into two major classes, known as user equilibrium(UE) and stochastic user equilibrium(SUE) [7, 8]. Considering the negative effect of road traffic congestion upon travel time, the user equilibrium was conceptualised by Wardrop [9]. Assuming that the travellers know the precise route cost and choose the route with minimum cost, the UE principle is reached when no traveller can reduce transportation cost by changing routes. To overcome the unrealistic assumption of precise perception of route travel time across travellers, the stochastic user equilibrium was firstly defined by Daganzo and Sheffi [10]. The SUE principle is obtained when no traveller’s perceived transportation cost can be reduced by unilaterally changing routes.

In the exist literatures, the SUE problem was classed into two types: the logit-based SUE and the probit-based SUE, according to those random costs following Gumble or normal distribution [11]. Due to explicit form and calculation, the logit-based SUE model has paid great attention [12, 13, 14, 15, 16, 17] . However, the probit-based SUE model behaves more appealing, attributing to the fact that it can take no account of overlapping, or correlated routes [18]. Sheffi and Powell proposed the well-known Method of Successive Average (MSA) with predetermined step sizes for solving the probit-based SUE problem[19]. Maher modified the model to decrease the computation complexity by choosing the optimal step length along the search direction[18]. Though distributed computing approaches were executed to reduce the computing time for the probit-based SUE problem [20], the development of probit-based SUE models still has some limitations, especially on large networks, which can be explained by the difficulty of completing path enumeration or Monte Carlo simulation.

Considering that performing the SUE principle in conventional probit-based models is difficult, we propose solving the stochastic traffic assignment problem by a Physarum-inspired model. The plasmodium of Physarum polycephalum is a large amoeboid organism, which contains a great number of nuclei and tubular structures [21]. These tubular structures will distribute protoplasm as a transportation network. Recently, it is shown that Physarum has the capacity of finding the short path between two points in a given labyrinth [22]. Tero et al. [23] inspired an mathematical model that can capture the basic dynamics of network adaptability through iteration of local rules and produces solutions with properties comparable to or better than those of real-world infrastructure networks. Bonifaci [24] has proved that the mathematical model can convergence to the shortest path. Later, the Physarum model was used to design and simulate transport network [25, 26, 27, 28], find the short path [29, 30, 31, 32]. To handle the uncertainty in the real application [33], the Physarum model can also solve shortest path under uncertain environment[34, 35].

Considering the continuity of the flow and protoplasmic network adaptivity, now the UE problem can also be solved by the Physarum model [36]. Note that UE problem is just a subsection of the SUE problem, here we present a Physarum-inspired model for the probit-based SUE problem. In the proposed model, the origin Physarum model is modified to adapt the directed network with multiple sources and directions and the link travel time is regard as the length of Physarum tubular structures.

This paper is organized as follows. In Section 2, the SUE assignment problem in traffic networks is reviewed and the Physarum polycephalum model is briefly introduced. In Section 3, a Physarum-inspired model for the SUE problem is presented. In Section 4, numerical examples are given to prove the rationality and convergence properties of the proposed model. Finally, the paper ends with conclusions in Section 5.

2 PRELIMINARIES

In this section, the basic theories, including the probit-based stochastic user equilibrium problem and Physarum polycephalum model, are briefly introduced.

2.1 Probit-based user problem

2.1.1 *Notations, assumptions and definitions *

Given a strongly connected transportation network $G=(N,A)$, where $N$ and $A$ denote the sets of nodes and links,respectively. Network attributes are denoted by notations as follow:

{eqlist*}

[\eqliststarinit]

Set of origin nodes, $R\subseteq N$

Set of destination nodes, $S\subseteq N$

An origin node, ${r\in R}$

An destination node, ${s\in S}$

Set of all the paths between OD pair rs.

Travel demand between OD pair rs, and all the OD travel demands are grouped into column vector, ${\bf{q}}={(\cdots,{q_{rs}},\cdots)^{T}}$, ${r\in R}$, ${s\in S}$

Traffic flow on path $k$ between OD pair rs, ${k\in{K}_{rs}}$.

Column vector of traffic flows on the paths between OD pair rs, ${{\bf{f}}^{rs}}={(\cdots,f_{k}^{rs},\cdots)^{T}}$, ${k\in{K}_{rs}}$.

Column vector of traffic flows on the all paths,${\bf{f}}={(\cdots,{{\bf{f}}^{rs}},\cdots)^{T}}$, ${r\in R}$, ${s\in S}$.

Travel time on path $k$ between OD pair rs, ${k\in{K}_{rs}}$

Column vector of traffic time on the paths between OD pair rs, ${{\bf{c}}^{rs}}={(\cdots,c_{k}^{rs},\cdots)^{T}}$, ${k\in{K}_{rs}}$

Column vector of traffic time on the all paths,${\bf{c}}={(\cdots,{{\bf{c}}^{rs}},\cdots)^{T}}$, ${r\in R}$, ${s\in S}$

Traffic flow on link a, $a\in A$

Column vector of all link flows, ${\bf{x}}{\rm{=}}{(\cdots,{x_{a}},\cdots)^{T}}$, ${a\in A}$.

Asymmetric travel time on link a, $a\in A$.

Column vector of all the link-travel-time functions, ${\bf{t}}{\rm{=}}{(\cdots,{t_{a}},\cdots)^{T}}$, ${a\in A}$

$\delta_{a,k}^{rs}=1$ if ${k\in K_{rs}}$ between OD pair rs traverses link ${a\in A}$, $\delta_{a,k}^{rs}=0$, otherwise.

link/path incidence matrix associated with OD pair rs, ${\bf\triangle^{rs}}=(\delta_{a,k}^{rs},a\in A,k\in{K}_{rs})$

link/path incidence matrix for the entire network, ${\bf\triangle}={(\cdots,{\bf\triangle}^{rs},\cdots)}$

According to the cost flow superposition principle, the path travel time can be valued as the summation of link travel time [19], which can be expressed as:

[TABLE]

Compactly, the relation between path travel time and link travel time can be expressed in vector form, namely:

[TABLE]

Assuming that the network users perceived link travel time is consist of the determined link travel time and random error term. The perceived link travel time is thus expressed as:

[TABLE]

where link travel time function ${t_{a}}({x_{a}})$ is positive, continuously differentiable and strictly monotone increasing. The error term $\varepsilon_{a}$ associated with link $a$ is a normally distributed random variable with zero mean for the probit-based SUE problem [19], which can be expressed as following:

[TABLE]

where $\gamma$ is a proportionality constant parameter and ${t_{a}^{0}}$ is a constant which usually equals free-flow link. Similarly, the perceived link travel time is also a normally distributed random variable, namely:

[TABLE]

Due to the linearity of the incidence relationships, the perceived path travel time also follows a multivariate normal distribution leading to the probit model for the path choice, which can be expressed as:

[TABLE]

Obviously, according to the accumulation of random variables, $C_{k}^{rs}$ is also a normally distributed random variable with $c_{k}^{rs}$ mean :

[TABLE]

Let $f_{k}^{rs}$ denote the traffic flow on path $k$ between OD pair rs, it can be expressed as the following equation:

[TABLE]

where $P_{k}^{rs}$ denotes the path choice probability for path $k$ between OD pair rs. On the basis of economics principles, $P_{k}^{rs}$ denotes the probability of path $k$ being the shortest one for given path travel time between OD pair rs [37], namely:

[TABLE]

The stochastic user equilibrium is reached when no user can reduce his perceived travel time by unilaterally changing routes [10]. The objective function of SUE problem was first proposed and proved by Sheffi and Power [19]:

[TABLE]

where we have used the results of Williams that

[TABLE]

2.1.2 The Method of Successive Average

The MSA algorithm developed by Sheffi and Powell [19] was the first algorithm applied to solve the SUE problem. In the MSA process, the link costs are calculated by the current link flows. An auxiliary link flow pattern is produced through a stochastic network loading procedure. And the search direction is obtained by the difference between the auxiliary link flow and the current link flow. The step size is predetermined by a descent sequence with respect to the iterations. The procedures of MSA method are summarized as following:

1. Step 1.1:

Choose initial link travel costs $\{t_{a}^{0},\forall a\}$, usually free-flow costs. Find an initial feasible flow pattern $\{x_{a}^{1},\forall a\}$ by carrying out, for example, a pure stochastic loading using mean costs. Set the iteration count ${n}$ to $1$. 2. Step 1.2:

According to the current flow pattern $\{x_{a}^{n},\forall a\}$, calculate the current travel costs $\{t_{a}^{n}(x_{a}^{n}),\forall a\}$. 3. Step 1.3:

Given the mean travel costs $\{t_{a}^{n},\forall a\}$ and the demands of OD pairs, find the auxiliary flow pattern $\{\hat{y}_{a}^{n},\forall a\}$ by carrying out a pure stochastic loading. 4. Step 1.4:

Calculate the new current solution according to the equation:

[TABLE] 5. Step 1.5:

Convergence test. If the following condition is fulfilled, then stop and output . Otherwise, ${n=n+1}$, go to step 2.

[TABLE]

The search direction is found by using the auxiliary flow pattern $\{\hat{y}_{a}^{n},\forall a\}$, which is computed through Monte Carlo simulation methods:

1. Step 2.1:

Initialize counter ${i=1}$. 2. Step 2.2:

Sample one realization from each link, using $f({T_{a}}|t_{a}^{n}))$. 3. Step 2.3:

Assign ”all or nothing” from each origin to each destination. This results in the auxiliary flow pattern $y_{a}^{(i)}$. 4. Step 2.4:

Average the flow for each link, $\bar{y}_{a}^{(i)}=[(i-1)\bar{y}_{a}^{(i-1)}+y_{a}^{(i)}]/i$. 5. Step 2.5:

If the stopping criterion is met, set ${\rm{\hat{y}}}_{a}^{n}{\rm{=}}\bar{y}_{a}^{(i)}\forall a$; If not,set ${i=i+1}$ and go to step 2.2.

where $f({T_{a}}|t_{a}^{n})$ is the probability density of ${T_{a}}$, and the travel time of link ${a}$ can be calculated according to Eq.(5) The stopping criterion referred to at step 2.5 may be based on the reduction of the variance of $\bar{y}_{a}^{(i)}$ as i grows [7, 19], such as a fixed number of drawings $I_{0}$.

2.2 Physarum polycephalum model

Physarum polycephalum is a single-celled amoeboid organism, which is also called as plasmodium in the vegetative phase. It is able to solve the shortest path selection, basing on its special foraging mechanism: the transformations of tubular structures and a positive feedback from flow rates. The high rates of the flow motivate tubes to thicken, and the diameter of the tube minishes at a low flow rate. A total introduction for the physarum polycephalum model is given below.

Supposing the shape of the network formed by the Physarum represented by a graph, plasmodial tube refers to an edge of the graph and a junction between tubes refers to a node. Assuming a set of nodes ${N}$, ${N_{1}}$ and ${N_{2}}$ are signed as the source and destination nodes, any others are labeled as ${N_{3}}$, ${N_{4}}$, ${N_{6}}$, ${N_{7}}$, etc. The edge connecting nodes ${N_{i}}$ and ${N_{j}}$ is remarked as ${M_{ij}}$. The flux from node $N_{i}$ to node $N_{j}$ through edge ${M_{ij}}$ is remarked as $Q_{ij}$, which we can expressed as [23] :

[TABLE]

where $\eta$ is the viscosity of the fluid and ${D_{ij}}{\rm{=}}\pi{\rm{r}}_{ij}^{4}/8\eta$ is measure of the conductivity of the edge $M_{ij}$ tube. $p_{i}$ is the measure of the pressure at the node $N_{ij}$ and $L_{ij}$ is the length of the edge of $M_{ij}$. According to the conservation law of flow, the inflow and outflow must be balanced, namely:

[TABLE]

For the source nodes $N_{1}$ and $N_{2}$, the flux equations can be denoted as:

[TABLE]

[TABLE]

where $I_{0}$ is the flux from the source node to the destination node, which is assumed as a constant in the model. According to the Eqs (14)-(17),the network Poisson equation for the pressure is derived as following:

[TABLE]

by further setting $p_{2}=0$ as the basic pressure level, the pressure of all nodes can be determined according to Eq.(18) and all $Q_{ij}$ can also be determined by solving Eq.(14).

To accommodate the adaptive behavior of the plasmodium, the conductivity $D_{ij}$ is assumed to change when adapting to the flux $Q_{ij}$. And tubes with zero conductivity will die out. The conductivity of each tube is described as the following equation [23]:

[TABLE]

where $\alpha$ is the decay rate of the tube and $f$ is monotonically increasing continuous function which satisfies $f(0)=0$. The Physarum can converge to the shortest path when $f(|{Q_{ij}}|)=\left|Q\right|$ and $\alpha=1$ [24]. Obviously, the positive feedback exists in the model.

3 PROPOSED METHOD

In this section, we employ the propoesd Physarum model to solve the stochastic user equilibrium problem. Generally speaking, there are three problems to be addressed:

The original Physarum model is used to solve the shortest path problem in undirected graphs (${L_{ij}=L_{ji}}$) while most network is directed graphs (${L_{ij}\neq L_{ji}}$) in real traffic assignment problem. 2. 2.

There is only one source node in the shortest path finding mode, but we should solve the traffic assignment problem with multiple sources and sinks. 3. 3.

The modified Physarum model should approach the optimal flow distribution in the traffic assignment problem.

3.1 Physarum-based model for the shortest path in the directed network

In the original Physarum model, each arc shown in Figure 1.a is bidirectional, which means the distance from node $i$ to node $j$ is same as that from node $j$ to node $i$. To solve the constrained shortest path problem, Wang et al. [30] proposed a modified Physarum model shown in Figure 1.b where each edge is regarded as two tubes with opposite direction and equal weight. And there is only one direction between two nodes, which means that the flux can flow from node $s$ to node $t$. The modified Physarum model has the ability to find the shortest path in the directed network. While most roads in the city have the properties of two-way traffic characteristic and opposite directions are separated with each other, the Physarum model modified by Wang et al. couldn’t work in these networks. Here we proposed a new modified Physarum model shown in Figure 1.c. There are two opposite directions between node $s$ and node $t$, and the length of two opposite directions is denoted by ${L_{st}}$ and ${L_{ts}}$. Basing on the feature of foraging behavior, the conductivity matrix ${D}$ implies not only the conductivity but also the direction of each tube, namely ${D_{ij}\neq D_{ji}}$ during iterations. In order to implement this idea into the original Physarum model, Eq. (18) is modified as following:

[TABLE]

To keep the validity of conductivity, the conductivity equation defined in Eq.(14) should be improved as following:

[TABLE]

Particularlly when ${L_{st}}=inf$ or ${L_{ts}}=inf$, which means that the flux can only flow from node $t$ to node $s$ or from node $s$ to node $d$, our modified model is the same as that of Wang et al. Exactly, the model modified by Wang et al. is a section of our modified model.

3.2 Physarum-based model for multiple sources and directions

In the original Physarum model, there is only one source node and one direction node. While in the stochastic user equilibrium problem, there are always multiple OD pairs. Assuming ${O}$ denoting the set of origin nodes, $O\subseteq N$, and ${D}$ denoting the set of destination nodes, $D\subseteq N$, we can modify Eqs. (16) and (17) as following:

[TABLE]

[TABLE]

where ${I_{o}}$ is the inflow at the origin node ${o}$,${I_{d}}$ is the outflow at the destination node ${d}$. To ensure the flow is distributed in a optimal way, here we use the modified model proposed by Zhang [28] to replace Eq.(18):

[TABLE]

3.3 Physarum-inspired model for the probit-based SUE problem

Now, we study how to solve the probit-based stochastic user equilibrium problem basing the Physarum-inspired model. Due to the feature of foraging behavior, the flow and the conductivity along each link are continuous in the process of Phyasrun approaching the shortest path. While in other classical shortest path algorithms, such as Dijkstra algorithm [38], Floyd algorithm [39], algorithms approach the shortest path by traversing all the nodes until the destination node is visited, which is totally uncontinuous.

Considering the continuity and dynamic reconfiguration of Physarum model, we can update the link travel time within each iteration. The flux will be redistributed by the modified Physarum model when the link travel time is updated during iterations. Here we adopt the following equation to update the length of link $a$:

[TABLE]

where $x_{a}$ denotes the traffic flow on link $a$ at the $(n)_{th}$ iteration, $C_{a}^{n}$ and $C_{a}^{n-1}$ represent the length of link $a$ at the $n_{th}$ and $(n-1)_{th}$ iteration. And the search direction of link length $C_{a}$ is guided by ${t_{a}}(x_{a})$. Note that in equilibrium, there will be $C_{a}=t_{a}(x_{a})$, which means the length of link $a$ equals the travel time along link $a$.

The main steps of the proposed method for the probit-based stochastic user equilibrium problem is presented in Algorithm 1. In the process of Monte Carlo simulation, we can use Physarum model to replace ”All or nothing” method to calculate auxiliary flow pattern. In the iteration, $C^{n}$ is the current link travel cost matrix at the $n$ iteration, $C_{ij}^{n}$ represents the current travel cost from node $i$ to node $j$, at the $n$ iteration. Particularly, $C^{0}$ is the free-flow link travel cost matrix, $C_{ij}^{0}$ represents the free-flow travel cost from node $i$ to node $j$.

Different from the MSA algorithm, the current travel cost($C_{ij}^{n}$) is calculated by the modified auxiliary flow in the proposed algorithm, attributed to the continuity of Physarum model. During the the process of Monte Carlo simulation the flow and the conductivity along each link are continuous.

More importantly, the current travel cost $C_{ij}^{n}$ doesn’t equal $C_{ji}^{n}$ when $\hat{Q}_{ij}^{(I_{0})}\neq\hat{Q}_{ji}^{(I_{0})}$, which means the same edge have different travel costs in two opposite directions. This peculiarity is rather important in nowadays traffic network. Because most roads in the city don’t interfere in two opposite directions, opposite directions are separated with each other. Indeed, the flow in edge $L_{ij}$ does’t influence the travel cost in edge $L_{ji}$.

4 NUMERICAL EXAMPLES

In this section, two examples are designed to prove the rationality and convergence properties of the proposed algorithm, a one source and sink node network and a multiple sources and sinks network. The inner iteration $I_{0}$ and the outer iteration $n$ are compared with those in the MSA algorithm.

Both tests are investigated using a simple network shown in Figure 2, which is introduced by Sheffi and Powell [19]. Link costs are calculated by the US Bureau of Public Roads (BPR) function, which is expressed as following:

[TABLE]

where parameters of the link cost functions for each link, $\alpha_{a}$ and$\beta_{a}$ are shown in Table 1.

According to Eq.(5), the value of $\gamma$ determines the variance of the perceived link travel cost ,which has a great effect on the convergence properties of the algorithm. In both tests, $\gamma$ is kept constant at 0.3. In Eq.(13), $\varepsilon_{0}$ is the condition of stopping iteration, if the value of $\varepsilon_{0}$ is to small, the process of both algorithms cost much calculating time. However, we can’t get the final travel flux if its value is too large. Hence, in order to compare the speed of convergence, the value of $\varepsilon_{0}$ is kept constant at 0.1 in the proposed algorithm and the MSA algorithm. And all computational experiments are executed using Matlab on Intel(R) Core(TM) i5-5200U processor (2.2Ghz) with 8.00 GB of RAM under Win 8.

4.1 Example 1

In this example, there is only one source-direction from node 1 to node 12 with travel rate of $q_{1,12}=20$ vehicles per unit time. To study the effect of the inner iteration $I_{0}$, we examined the flow on particular link corresponding the number of inner iterations. Here, we choose the traffic flow of link $L_{6,7}$ in both algorithms. The effect of inner iterations is illustrated in Figure 3.

It’s clearly that the convergence per equilibrium iteration improves when the inner iteration $I_{0}$ augments. However, in the proposed algorithm, the deviation between lower iteration and equilibrium flow is obviously smaller than that in the MSA algorithm when the counter of outer iterations is small.

What’s importantly, both algorithms can get the same equilibrium flow when the inner iterations are different. Sheffi and Powell has proved that very few inner simulation iterations (possibly just one) may be sufficient to achieve a reasonable convergence rate of the equilibrium iterations by the SMA algorithm [19]. Note that the proposed algorithm can also get the similar result when the inner iteration equals 1 in Figure 3, we speculate that one iteration can also achieve a reasonable convergence rate of the equilibrium iterations. The numerical example is presented as below.

To prove the rationality and convergence properties of the proposed algorithm, the inner iteration $L_{0}$ should be kept same in both algorithms. So we should keep the inner iteration $I_{0}$ equal 1 in both algorithm. The results of link traffic flow calculated by both algorithms are shown in Table 2.

It’s obviously there are same traffic paths in the network. And traffic rates in each link are almost similar calculated by both algorithms, which differ by no more than $0.32$ vehicles per unit time. Besides, the convergence rate of the proposed algorithm is faster than that of the SMA algorithm. This peculiarity will become much more important especially when the network is quite large. Considering that responsiveness of the traffic flow assignment is much more significant in nowadays traffic network, the proposed algorithm contributes a positive idea to reduce computing time.

4.2 Example 2

In this example, we also used the traffic network shown in Figure 2 and Table 1. Different from Example 1, the origin-destination demands are assumed as $q_{1,12}=10$ and $q_{1,8}=10$, which denote the rate of vehicles per unit time. The inner iteration was set as $10$ in this example. The results of link traffic flow calculated by both algorithms are shown in Table 3. Clearly, traffic rates calculated by both algorithms are also similar and the proposed method obviously get the equilibrium flow faster than the MSA algorithm. The maximum error of the link flows is no more than $0.35$ vehicles per unit time. The computing time don’t increase compared with that in Example 1.

5 CONCLUSIONS

Considering of deviation between traveller’s perceived transportation cost and actual cost, the stochastic user equilibrium is much more significant than user equilibrium. This paper presents a Physarum-inspired model for the probit-based stochastic user equilibrium problem. The Physarum model is modified to solve the SUE problem in the first time. To satisfy the characteristic of the real traffic networks, the origin Physarum model is modified to find the shortest path in direction networks with multiple sources and directions. Considering the foraging behavior of Physarum, the Physarum could find the shortest travel time path between each OD pair. The equilibrium flows could be obtained whenPhysarum couldn’t find a shorter travel time path.

We compared the proposed algorithm with the MSA algorithm. And Numerical results showed that the proposed algorithm can effectively achieve the SUE solution in practice. If the inner iteration properly assigned, the proposed algorithm is faster and more efficient than the MSA algorithm. Note that many investigations about paralleled Physarum model have been achieved [40], the time consumption of the proposed algorithm will obviously reduced in concurrent computation. Besides, the proposed method is easy to combine with other algorithms [41].

ACKNOWLEDGMENTS

The work is partially supported by National Natural Science Foundation of China (Grant No. 61671384), Natural Science Basic Research Plan in Shaanxi Province of China (Program No. 2016JM6018), the Fund of SAST (Program No. SAST2016083), the Seed Foundation of Innovation and Creation for Graduate Students in Northwestern Polytechnical University (Program No. Z2016122).