Differential Evolution Based Numerical Variable Speed Limit Control Method with a Non-Equilibrium Traffic Model

Abstract

This paper introduces a numerical variable speed limit (VSL) control method on a motorway, modeled by the system of partial differential equations (PDEs) of a non- equilibrium continuum traffic model. The method consists of a macroscopic simulation (i.e., numerical solution of the system of PDEs of the continuum model), introduction of the solution-based cost function and numerical optimization with a differential evolution algorithm (DE). Due to the numerical solution scheme, the method enables application of a wide range of continuum traffic models without prior discretization of PDEs. In this way, the method overcomes the limitations of the basic continuum models and represents a step towards more accurate traffic modelling in control strategies. In this paper, we determine optimal variable speed limits with the DE algorithm on a motorway section modeled by the modified switching curve model, which is a non-equilibrium continuum model consistent with the three-phase traffic flow theory. The effectiveness of the determined variable speed limits is validated using microsimulations of the test section, which show promising reductions of queue lengths and number of stops.

1. Introduction

Extensive research has been conducted on the traffic breakdown mechanism ever since the first continuum traffic flow model LWR was developed [1,2]. The development of Intelligent Transportation Systems with their wide potential to control traffic emphasizes the significance of appropriate control policies, aiming at traffic breakdown prevention. Emerging technologies, such as intelligent speed adaptation, promise higher variable speed limit (VSL) observance and thus greater potential for influencing traffic flow. Using VSLs as a control measure has been the subject of scientific interest in recent years. Various control strategies have been proposed aiming to establish a control strategy for preventing traffic breakdown based on model-prediction of traffic flow dynamics, such as model predictive control to minimize travel time [3,4,5] and proactive variable speed limit control [6,7,8]. Another type of control strategy, which is concerned with maximizing freeway throughput by reducing the effect of capacity drop at bottlenecks, is feedback-based mainstream control [9,10]. The results of all these strategies are continuous values for VSLs, which are later discretized. As a response to significant loss of performance due to subsequent discretization of VSLs, a hybrid model predictive control with discrete speed limit signals using genetic algorithms for optimization was proposed in [11]. Genetic algorithms have been identified as a powerful optimization tool that has been shown to outperform sequential quadratic programming for the VSL control problem [12].

For the description of traffic dynamics in all of the aforementioned control strategies, METANET [13] is used, which is based on Payne’s [14] continuum traffic flow model. However, this model is known to have some severe drawbacks and inconsistencies, the most important being the violation of the anisotropy principle [15]. This led Aw and Rascle [16] and Zhang [17] to develop a new model, also called the ARZ model, which has also been applied in an optimal control strategy [18]. While the ARZ model omits the drawbacks of Payne-type models, it has been shown that it loses the ability of the Payne model in stimulating unstable traffic and vehicle clusters [19,20]. Another model that has been used in optimal control strategies is the cell transmission model CTM [21,22,23,24], which is a numerical scheme for solving problems associated with the LWR (or modified LWR) model. A more comprehensive review of VSL control strategies can be found in [25,26,27].

Several alternative control mechanisms to omit modelling traffic dynamics have been proposed such as rule-based algorithms for prescribing VSLs based on pre-specified threshold values of traffic characteristics [28], fuzzy logic [29,30] and reinforcement learning [31,32,33,34,35]. In [36], a decision support model for VSL control in recurrent congestion based on simulation of 1024 traffic environments was developed, emphasizing the large variety of possible traffic situations. However, in reality it is nearly impossible for an exact traffic situation to emerge (and evolve) twice. Moreover, accurate prediction of traffic flow dynamics enables timely control, i.e., controlling traffic flow before the actual onset of a shockwave.

This highlights the importance of accurate traffic flow modelling to capture the phenomena of traffic breakdown. According to recent extensive research [37,38], traffic breakdown has a more sophisticated nature than first thought. The findings of the above mentioned research suggest that the true nature of traffic flow is better described by Kerner’s three phase traffic flow theory [39,40] with three distinct steady phases, namely free flow phase, synchronized flow and wide moving jams. The observed breakdown characteristics are conceptually different from the breakdown description of the previous continuum macroscopic traffic flow models, such as LWR [1,2], Payne [14] and Aw-Rascle [16], which represent basically all of the models used for model predictive control policies [3,4,5,6,7,8,9,10,11,12,18,19,20,21,22,23,24]. The authors have emphasized that empirical features of traffic breakdown should not be neglected in freeway control strategies [28,40]. While non-equilibrium continuum traffic flow models, that are able to describe (at least most of) important traffic flow characteristics (i.e., capacity drop, stop-and-go waves, multiple-regime, three-phase traffic theory, anisotropy) have been developed [41,42,43,44,45], to our knowledge no control policy has been proposed for such a model.

The gap between existing continuum traffic flow models and continuum models used in VSL control strategies is evident from the above literature review. In the usual approach, the PDEs of the continuum traffic flow model are first discretized to avoid differential formulation and to enable the formulation of the control strategy. This PDE discretization can be a difficult mathematical task and can greatly reduce the available set of continuum traffic flow models and/or lose some model properties. However, as the main objective of the VSL control is to prevent or at least minimize traffic breakdown, the need for more accurate modelling, consistent with empirical features of traffic breakdown, is obvious.

The present research is motivated by the described modelling gap. Therefore, the focus of the research is the modelling part of the control strategy. While the purpose of control strategies is to prevent traffic breakdown, the objective of this research is to consider the empirical features of traffic breakdown when developing a control strategy, i.e., to use a non-equilibrium continuum model that captures empirical features of the traffic breakdown. The main contributions of the paper are summarized as follows:

• We introduce a novel numerical VSL control method that opens control strategies for a wide range of physically and empirically proven continuum models in their original differential formulation.
• We consider three phase traffic flow theory in model predictive control and provide an insight into the effects of the synchronized flow phenomena on control strategies.
• We define a cost function based on the numerical solution of the continuum model with the aim of improving traffic flow in two ways: improving network performance and improving traffic safety.

The paper is organized as follows: Section 2 with its subsections provides the detailed description of the methodology, Section 3 presents a case study, Section 4 discusses the results and Section 5 summarizes the findings of this paper and considers future research directions.

2. Methodology

In the proposed method, we solve the system of PDEs of the continuum traffic model with a numerical scheme. This transfers the whole problem formulation into a more sophisticated computational environment, so special consideration for the numerical form of the cost functional along with a powerful optimization tool is required. The proposed VSL control strategy therefore consists of two major parts: macroscopic modeling and optimization. The first macroscopic modeling part is explained in detail in Section 2.1 and Section 2.2 and the optimization part is described in Section 2.3 and Section 2.4. The conceptual and methodological framework is presented in the flowchart in Figure 1.

In the macroscopic modeling part, the systems of PDEs of a continuum traffic model are solved using a numerical solution scheme. In the optimization process, we are concerned with the optimal control problem. First, we define the cost function, based on a numerical solution of the traffic model. Secondly, we use a differential evolution (DE) algorithm [46,47] to perform the optimization and determine the optimal set of VSLs.

In the proposed method, VSLs are determined in such a way that the value of the cost function, obtained from the numerical solution of the macroscopic continuum model, is optimized. Consequently, traffic optimization with the determined VSLs in the macroscopic model is self-explanatory. For this reason, we have additionally validated the effectiveness of the determined VSLs with microsimulations of the test section to capture the effect of speed limits on driver behavior.

2.1. Numerical Solution of the Continuum Traffic Model

Continuum traffic flow models are mathematical models that use macroscopic traffic flow description, i.e., traffic flow is described with characteristics at an aggregate level. They consist of nonlinear systems of hyperbolic PDEs of type:

$$\frac{\partial }{\partial x}u+\frac{\partial }{\partial t}f\left(u\right)=g\left(u\right),$$

(1)

where $u$ is a m-dimensional vector of macroscopic quantities and $f\left(u\right)$ and $g\left(u\right)$ are m-dimensional vector functions. To find the solution to the system of PDEs of the continuum macroscopic model in the conservative form of type (1), we employ the MUSCL-type finite volume method in combination with a two-stage Runge-Kutta scheme and dynamic timestep adaptation, using the CFL condition [48,49], which is a second order accurate numerical method for solving systems of PDEs.

The systems of PDEs of the continuum model is considered with appropriate boundary conditions. At the downstream boundary we introduce Neumann boundary condition, which assumes unchanged traffic state at the boundary:

$$\frac{\partial u}{\partial x}\left(L,t\right)=0,$$

(2)

where $L$ is the length of the road section under consideration. We can simulate real traffic by imposing a Dirichlet boundary condition at the upstream boundary, i.e., the inputs for traffic simulations are the measured traffic flow characteristics arriving from the previous road section:

$$u\left(0,t\right)=u\left(t\right).$$

(3)

However, this needs to be done carefully, as the Dirichlet boundary condition is not able to handle queue spillback. This can be managed by dynamically switching between Dirichlet and Neumann boundary condition based on the current traffic conditions at the upstream boundary [50].

Furthermore, the effect of VSLs on traffic behavior must be considered in the solution of the continuum traffic flow model. The usual approach, where the VSL is considered as the new free flow speed parameter of the fundamental diagram, was argued to be inaccurate [5], since the VSL should not affect the entire speed–density relationship, i.e., in the fundamental part of the diagram, where density dependent speed is lower than the VSL, the VSL should have no effect. In our study, the imposed VSL is considered as the upper bound for velocity, whereas the speed-density relationship still follows the same fundamental diagram [51]. This is set by equation:

$$v^e\left(x,t\right)=\min \left(VSL,v^e\left(\rho (x,t\right)\right)),$$

(4)

where $VSL$ stands for the imposed VSL and $v^e\left(\rho (x,t\right))$ is the equilibrium velocity determined for the no-VSL case. In case of several equilibrium speeds (i.e., a two or more regime fundamental diagram), the above equation applies to each, i.e., for a two regime fundamental diagram we have:

$$v_{1,VSL}^e\left(\rho ,VSL\right)=\min \left(VSL,v_1^e\left(\rho \right)\right)$$

(5)

and

$$v_{2,VSL}^e\left(\rho ,VSL\right)=\min \left(VSL,v_2^e\left(\rho \right)\right).$$

(6)

2.2. Modified Switching Curve Model

In this research, we choose to apply the “switching curve” model [42], with modified relaxation term [44] so that the model is capable of reproducing the features of three-phase traffic theory and is able to describe other empirically observed features (i.e., capacity drop, stop-and-go waves, multiple-regime, three-phase traffic theory, anisotropy).

The governing equations of the model are:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho v\right)}{\partial x}=\frac{\rho v}{l}\frac{\partial l}{\partial x}+{\nu }_{rmp}\left(x\right),$$

(7)

$$\frac{\partial v}{\partial t}+\left(v+\rho \frac{v_1^e\left(\rho \right)}{d\rho }\right)\frac{\partial v}{\partial x}=\left\{\begin{array}{l}\frac{v_1^e\left(\rho \right)-v}{\tau },\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{if} \rho <{\rho }_{min}^{syn}\mathrm{or}v>v_1^e\left(\rho \right) \mathrm{and} \rho <{\rho }_{max}^{free},\\ \frac{1}{\tau }{v}^{*}\left(\rho \right),\hspace{1em}\hspace{1em}\mathrm{if} v_2^e\left(\rho \right)<v<v_1^e\left(\rho \right) \mathrm{and} {\rho }_{min}^{syn}<\rho <{\rho }_{max}^{free},\\ \frac{v_2^e\left(\rho \right)-v}{\tau },\hspace{1em}\mathrm{if} \rho >{\rho }_{max}^{free}\mathrm{or}v<v_2^e\left(\rho \right) \mathrm{and} {\rho }_{min}^{syn}<\rho <{\rho }_{max}^{free}.\end{array}\right.$$

(8)

In Equations (7) and (8), $\rho \left(x,t\right)$ stands for traffic density and $v\left(x,t\right)$ stands for velocity, which represent two basic traffic flow characteristics. Equation (7) resembles the basic LWR model, with right side modified to account for on-ramps, off-ramps and changes in number of lanes [50]—it is only nonzero within the merging zone of on- and off ramps or within the transition zones, where some of the lanes drop or new lanes begin. The first term $\frac{\rho v}{l}\frac{\partial l}{\partial x}$, accounts for changes in the number of lanes (notation $l$). The second term ${\nu }_{rmp}$ describes on- and off- ramps and is equal to $\frac{Q_{rmp}}{l{L}_{rmp}}$ within the merging zones, where $Q_{rmp}$ is the on-ramp inflow ($Q_{rmp}>0$) or the off-ramp outflow ($Q_{rmp}<0)$, $l$ is the number of lanes (not counting the ramp itself as a lane) and $L_{rmp}$ is the ramp length (in the case of a numerical method, $L_{rmp}$ must be a multiple of distance between mesh points). The right side of Equation (8) is the relaxation term with relaxation time $\tau$, defined based on the traffic density region, and $v_{e1}\left(\rho \right)$ and $v_{e2}\left(\rho \right)$ represent two distinct equilibrium velocities. The latter represent the basic feature of the switching curve model, accounting for overtaking difficulties at higher densities [42]. The basis of the model modification in [44] is $v^{*}$, which stands for linear interpolation between speeds $v_1^e\left(\rho \right)-v$ and $v_2^e\left(\rho \right)-v$. The boundaries for relaxation term definitions ${\rho }_{min}^{syn}$ and ${\rho }_{max}^{free}$ correspond to regions in the fundamental diagram, i.e., ${\rho }_{min}^{syn}$ represents the minimum density below which synchronized flow cannot occur and ${\rho }_{max}^{free}$ represents the limit density for free flow existence.

We rewrite the model using a new variable $y=\rho \alpha$ and put it into conservative form of type (1), so we can solve it with the aforementioned numerical method.

$$\mathit{u}=\left[\begin{array}{c}\rho \\ y\end{array}\right], \mathit{f}\left(\mathit{u}\right)=\left[\begin{array}{c}y+\rho v_1^e\left(\rho \right)\\ y\left(\frac{y}{\rho }+v_1^e\left(\rho \right)\right)\end{array}\right], \mathit{g}\left(\mathit{u}\right)=\left[\begin{array}{c}\frac{y+\rho v_1^e\left(\rho \right)}{l}\frac{\partial l}{\partial x}+{\nu }_{rmp}\left(x\right)\\ K_{\tau }\left(\rho , v\right)\end{array}\right],$$

(9)

where $K_{\tau }\left(\rho , v\right)$ is the relaxation term in conservative form:

$$K_{\tau }(\rho ,v)=\left\{\begin{array}{c}\frac{-y}{\tau },\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{if} \rho <{\rho }_{min}^{syn}\mathrm{or}v>v_{e1}\left(\rho \right) \mathrm{and} \rho <{\rho }_{max}^{free},\\ \frac{\rho v^{*}\left(\rho \right)}{\tau },\hspace{1em}\hspace{1em}\mathrm{if} v_2^e\left(\rho \right)<v<v_1^e\left(\rho \right) \mathrm{and} {\rho }_{min}^{syn}<\rho <{\rho }_{max}^{free},\\ \frac{\rho v_2^e\left(\rho \right)-{\rho v}_1^e\left(\rho \right)-y}{\tau },\hspace{1em}\mathrm{if} \rho >{\rho }_{max}^{free}\mathrm{or}v<v_2^e\left(\rho \right) \mathrm{and} {\rho }_{min}^{syn}<\rho <{\rho }_{max}^{free}.\end{array}\right.$$

(10)

The function forms of the equilibrium velocity curves $v_1^e\left(\rho \right)$ and $v_2^e\left(\rho \right)$ and its parameters, the boundaries for relaxation term definitions ${\rho }_{min}^{syn}$ and ${\rho }_{max}^{free}$, and the relaxation time $\tau$ are subject to calibration of empirical data.

2.3. Optimal Control Problem

Optimal control theory is a rich mathematical theory concerned with finding a control law for a dynamic system such that a certain optimality criterion is achieved. A dynamical system is in general described with partial differential equations

$$\dot{x}\left(t\right)=f\left(t, x\left(t\right), u\left(t\right)\right)$$

(11)

where $t$ is time, $x\left(t\right)$ is the state function and $u\left(t\right)$ is the control function. In our case the dynamic system is described with a system of PDEs for the continuum traffic flow model. The cost functional is generally written as

$$J=\Phi \left(T,x\left(T\right)\right)+{{\displaystyle \int }}_0^Tℒ\left(t, x\left(t\right), u\left(t\right)\right)dt,$$

(12)

where $T$ is the terminal time, $x\left(0\right)=x_0$ is the initial state, $\Phi$ is the terminal cost and $ℒ$ is the running cost. Even proving the existence of the optimal solution to such problems is a difficult mathematical task, let alone finding the solution. For further mathematical considerations we refer to [52].

Colombo and Groli [53] studied optimal control problems related to continuum traffic flow models by minimizing the total variation of the traffic flow density to optimize traffic flow. They proved analytically the existence of an optimal solution. However, due to the complex form and specificity of the proposed cost functional, we cannot directly apply their results. Moreover, the total variation of the density in case of VSL control does not behave well enough for stable computations. In our research, we decided to use a cost functional that behaves more stably, yet is still concerned directly with traffic density. We consider the differences between the actual traffic flow density $\rho \left(x,t\right)$ and some optimal density ${\rho }_{opt}$. Herein, the optimal density is considered as some acceptable value of density, that does not yet have negative impact on traffic flow and is not desired to be exceeded. This means that the cost functional, to be minimized, should consist of densities higher than this parameter value, so we finally introduce the cost functional:

$$J=\{\begin{array}{c}{{\displaystyle \int }}_0^L{{\displaystyle \int }}_0^T\left(\rho \left(x,t\right)-{\rho }_{opt}\right)dt dx, \mathrm{if} \rho \left(x,t\right)>{\rho }_{opt}, \\ 0, \mathrm{otherwise},\end{array}$$

(13)

where $L$ is the total length of the road section and $T$ is the length of the analyzed time period.

The solution of systems of PDEs is obtained numerically, i.e., the density function is given in the form of solutions in $n$ discrete points $x_i$, where $n$ is a positive integer. Let $m$ denote the total number of time steps $t_j$ and $\rho \left(x_i,t_j\right)$ the solution of the continuum macroscopic traffic flow model, where $i=1,\dots ,n$ and $j=1,\dots ,m$. The corresponding numerical form of the proposed cost functional $J$ used in numerical computations can be written as

$$J=\{\begin{array}{c}{\displaystyle \sum _1^n}{\displaystyle \sum _1^m}\left(\rho \left(x_i,t_j\right)-{\rho }_{opt}\right)∆x_i∆t_j, \mathrm{if} \rho \left(x_i,t_j\right)>{\rho }_{opt,} \\ 0, \mathrm{otherwise}.\end{array}$$

(14)

The term $∆x_i∆t_j$ denotes the size of the cell in the numerical method for solving partial differential equations, which does not affect the optimization process. The final form of our cost functional $J$ used in numerical computations is therefore

$$J=\{\begin{array}{c}{\displaystyle \sum _1^n}{\displaystyle \sum _1^m}\left(\rho \left(x_i,t_j\right)-{\rho }_{opt}\right), \mathrm{if} \rho \left(x_i,t_j\right)>{\rho }_{opt,}\\ 0, \mathrm{otherwise}.\end{array}$$

(15)

2.4. Differential Evolution Optimization Algorithm

Differential evolution (DE) [46,47] is a powerful optimization tool that is well suited for numerical optimization problems and a discrete environment. The (DE) algorithm does not require analytical calculation of gradient, which means we can use the above defined cost function based on the numerical solution of the system of PDEs. Based on convergence behavior [54,55] and theoretical results [56,57], we selected variant rand/1/bin, with the following configuration parameters [58,59]: crossover rate $CR=0.7$, mutation factor $F=0.8$ and population size $\mathrm{NP}=20$.

The DE rand/1/bin algorithm is:

• Initialization: generate initial population of vectors $x_i^t$ In our case, the components of vectors $x_i^t$ are candidate values for VSLs at successive discrete points.

• Mutation: for each individual $x_i^t,$ DE creates a donor vector $v_i^t$ the $t-$ th generation:
  $$v_i^t=x_{r_1}^t+F\left(x_{r_2}^t-x_{r_3}^t\right),$$

  (16)

  where the indices ${\mathrm{r}}_1$, ${\mathrm{r}}_2$ and ${\mathrm{r}}_3$ are uniformly random integers mutually different and distinct from the running index $i$.

• Crossover: In crossover, the trial vector $u_i^t$ is produced using the binomial crossover;
  $$u_{i,j}^t=\{\begin{array}{c}{v}_{i,j}^t, if rand\left(0,1\right)\le CR or j=j_{rand},\\ x_{i,j}^t, otherwise, \end{array}$$

  (17)

  where $rand\left(0,1\right)$ is a uniformly distributed random number within $\left[0,1\right]$ and $j_{rand}$ is a random integer in $\left[1,\mathrm{n}\right]$.

• Selection: a new population is generated with regard to the value of the cost function;

$$x_i^{t+1}=\{\begin{array}{c}{u}_i^t, if J\left(u_i^t\right)<J \left(x_i^t\right),\\ x_i^t, otherwise. \end{array}$$

(18)

If the termination condition is satisfied, DE rand/1/bin stops;, or else a new iteration from step 2 starts.

The framework of the proposed VSL control method is the described DE rand/1/bin algorithm with the cost functional, described with Equation (13), with ${\rho }_{opt}={\rho }_{min}^{syn}$. The numerical solution of the modified switching curve model for different VSL candidates is calculated in each DE iteration with appropriate boundary conditions and on-ramp/of-ramp inflow to account for the measured traffic volumes. The code is written in MATLAB version 9.3 (R2017b, The Mathworks, Inc., Natick, MA, USA). For purposes of real time computation, the code enables the use of multiple cores (if available) for optimization.

It is important to note, that in general the described DE optimization algorithm is independent of the model equations and depends solely on its solution, which has an important consequence; the proposed optimal control can be used for any continuum model (with incorporated speed limits), that can be solved with any numerical scheme.

3. Case Study

With the above-described method, we determined VSLs in two discrete points on a motorway section in the western Ljubljana ring in Slovenia that represents one of the busiest motorway sections in the country. The analyzed section is fully equipped with video detection, microwave sensors and inductive loops within the traffic management and control system (TMCS) and traffic monitoring of the section is constantly performed. Two full-graphic variable message portals are installed on the section and these enable application of VSLs. All the on-ramps, as well as the motorway diverge represent bottlenecks, can become active at any time. The aim of the case study is to determine optimal VSLs on the two portals locations and analyze the effectiveness of such VSL control. The scheme of the analyzed motorway section is shown on Figure 2.

For generality purposes, VSL control is determined and tested for five different traffic volume scenarios. Traffic volumes in the first scenario were selected from the TMCS database in such a way that one hour traffic volumes in locations A (main direction) and B (on-ramp from Gorenjska motorway) represent one of the highest traffic volume scenarios—highly busy AM peak hour. We should note that these volumes are at the capacity of the motorway interchanges prior to the analyzed section (i.e., this is the maximum throughput of these interchanges and higher demands result in congestion outside of the analyzed section). By reducing traffic volumes at the two locations we obtained four different scenarios:

• Scenario 2: Traffic volumes at location B are reduced by 5%;
• Scenario 3: Traffic volumes at location A are reduced by 5%;
• Scenario 4: Traffic volumes at location B are reduced by 10%;
• Scenario 5: Traffic volumes at location A are reduced by 10%.

All the other volumes remained the same in all four scenarios. In this way, we varied both size and distribution of traffic volumes.

3.1. Model Parameters Calibration

The modified switching curve model enables calibration of the field data by determination the function forms of the equilibrium velocity curves $v_1^e\left(\rho \right)$ and $v_2^e\left(\rho \right)$ and its parameters, the boundaries for relaxation term definitions ${\rho }_{min}^{syn}$ and ${\rho }_{max}^{free}$ and the relaxation time $\tau$. Note that these parameters define the fundamental diagram and their calibration enables adjustment to local traffic flow properties/situation, such as external conditions (e.g., weather conditions [60]), traffic regulations, properties of the road (width of the lanes, grade), the composition of the flow, etc.) In the calibration process, the data collected under prevailing weather conditions were considered.

The extensive data from the TMCS system could theoretically be used as calibration data. However, the usual microscopic calibration schemes may be unsuccessful in a complex macroscopic situation, due to the objective function containing many secondary minima and fluctuations; therefore, it is usually better to calibrate to global properties of traffic flow, such as travel time or jam front propagation [50].

Therefore, we decided to use an indirect approach to field data calibration. Initial parameter values and function shapes were determined to approximately fit the TMCS data. Then we prepared microsimulations for five traffic volume scenarios with 30 random seeds in PTV VISSIM (PTV Group) software, calibrated it to TMCS data and collected the average travel times and queue lengths. Then we prepared macrosimulations in which we varied all the parameters to best fit the microsimulation results. For higher equilibrium velocity $v_1^e\left(\rho \right),$ we choose the Northwestern model [61].

$$v_1^e\left(\rho \right)=v_0{e}^{{\left(-0.5\left(1-\frac{\rho }{ {\rho }_{ref}}\right)\right)}^2},$$

(19)

where $v_0$ is free flow speed and ${\rho }_{ref}$ reference density. For $v_2^e\left(\rho \right)$ we choose logarithmic relationship:

$$v_2^e\left(\rho \right)=a_1{\left(\ln \left(\frac{{\rho }_{max}}{\rho }\right)\right)}^{a_2},$$

(20)

where ${\rho }_{max}$ is maximum density, with $a_1$ and $a_2$ as additional calibration parameters.

Parameters after calibration for the test section are shown in

Table 1. Model parameters after calibration.

Parameter	Section 1	Section 2
$v_0 \left[\mathrm{km}/\mathrm{h}\right]$	115	110
${\rho }_{ref}\left[\mathrm{eov}/\mathrm{km}\right]$	50	50
${\rho }_{max}\left[\mathrm{eov}/\mathrm{km}\right]$	140	140
$a_1\left[\mathrm{km}/\mathrm{h}\right]$	35	28
$a_2$	1.40	1.68
${\rho }_{min}^{syn}\left[\mathrm{eov}/\mathrm{km}\right]$	24	24
${\rho }_{max}^{free}\left[\mathrm{eov}/\mathrm{km}\right]$	45	45
$\tau \left[\mathrm{s}\right]$	15	15

(see Figure 2 for subsections).

3.2. Macroscopic Simulations

The calibrated parameters were then used for macroscopic simulations, i.e., numerical solutions of the modified switching curve model for five different traffic volume scenarios. TMCS 1-min data in scenario 1 and the reduced traffic volumes in all the other scenarios are used as the boundary conditions and on-ramp/off-ramp flows. Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7 show density in time in space in different traffic volume scenarios. Stop and go waves can clearly be seen from the density fluctuations. As expected, severity of the traffic jams differs among scenarios, and in scenario 5 no traffic jam appears, while synchronized flow is present. This will enable a more illustrative effect of the VSL control in different conditions.

4. Results and Discussion

Finally, we used the proposed VSL control method to determine optimal VSLs in all five traffic volume scenarios. The VSL candidates’ solution set consisted of VSLs from 60 km/h to the stationary speed limit of the section with increments of 10 km/h. While VSLs lower than 60 km/h could theoretically be beneficiary, they can be inappropriate outside of congested areas from a drivers’ perspective. Moreover, experiments show that speed limits lower than 60 km/h are not found optimal even if admissible, since they often result in a new shockwave upstream of the speed limit location. Obviously, the new shockwave would directly result in high cost function value due to high densities. Therefore, such low VSLs are eliminated by the optimization process. It should be noted that prescribing low speed limits within a congestion area for safety reasons is of a different nature. As frequent and abrupt VSL changes can be unsafe [62], the change in two successive VSLs (both in time and space) was limited to 20 km/h.

Experiments show fast optimization for the selected test section—in all cases at most 15 s due to multi-core optimization, which enables real time application.

Figure 8 shows macroscopic simulation of the traffic volume scenario 1 with determined VSL control. One can observe (probably a too) promising theoretical result, as the stop and go waves are completely eliminated (however, density is in the synchronized flow region, therefore, the traffic flow is not fluent). Since this is the result for our highest traffic volumes scenario, the macroscopic simulations for other scenarios are obvious—no stop-and-go waves—and are not presented.

From the results of the macroscopic simulation we can also calculate the total travel time $TTT$ as [50]

$$TTT={{\displaystyle \int }}_0^L{{\displaystyle \int }}_0^T\rho \left(x,t\right)dt dx.$$

(21)

In the numerical form, the $TTT$ is calculated as

$$TTT=\sum _1^n\sum _1^m\rho \left(x_i,t_j\right)∆x_i∆t_j.$$

(22)

Note that by comparing these two equations with Equations (13) and (14), we can obtain physical explanation for our cost function—the cost function to be minimized is the total travel time spent in a traffic state above optimal density.

Table 2. Macroscopic simulation results.

Scenario	$\mathbf{Lane}\mathbf{Averaged}\mathit{T}\mathit{T}\mathit{T}\mathbf{with}\mathbf{VSL}$	$\mathbf{Lane}\mathbf{Averaged}\mathit{T}\mathit{T}\mathit{T}\mathbf{w}/\mathbf{o}\mathbf{VSL}$	Difference
Scenario 1	90.4 h	102.4 h	−13.0%
Scenario 2	86.2 h	91.0 h	−5.3%
Scenario 3	74.0 h	80.6 h	−8.2%
Scenario 4	78.7 h	78.6 h	+0.1%
Scenario 5	67.8 h	64.5 h	+5.1%

shows the calculated $TTT$ for all the five traffic volume scenarios for the VSL controlled and uncontrolled case.

Table 2 shows that total travel time is decreased by VSL control in the first three scenarios, remains basically unchanged in the fourth scenario and increases in the last scenario. As mentioned in the previous section, a wide moving jam in scenario 5 does not occur, while there are several instances of synchronized flow phase (see Figure 7). According to Kerner’s three-phase traffic flow theory, only small local delays occur in the synchronized flow phase due to short stops. Prescribing VSL control oriented at achieving more fluent traffic with fewer velocity oscillations in this scenario apparently leads to a slight extension of the travel time.

Since macroscopic models represent traffic flow in an aggregate manner, they are basically the aggregation of individual vehicles. Therefore, these models assume that all vehicles and drivers are essentially the same or at most there exists a limited (very small) number of vehicle/driver types. Moreover, the speed limit observance rate is 100%, which could be theoretically set to a lower percentage but without distinguishing individual driver behavior. Although macroscopic modelling is a useful tool for determining appropriate VSL control, the effectiveness of such VSL control should be estimated at a more detailed level. Even more, we believe that the stochastic nature of traffic flow, which is apparent from both real-life observation and microsimulations with different random seeds, should not be overlooked. Therefore, we tested the proposed VSL control algorithm on a sample of 30 microsimulations of the test section (30 random seeds) for all the five traffic volume scenarios in PTV VISSIM (PTV Group). Results are presented in

Table 3. Microsimulation results and comparison.

Scenario	Performance Parameter	Mean with VSL	Mean w/o VSL	T Value	Mean Difference
Scenario 1	Total travel time	240.5 h	290.4 h	−4.810	−17.2%
	Total stops	2006.3	6350.2	−6.879	−68.4%
	Queue length	1222.1 m	2257.0 m	−5.724	−45.9%
Scenario 2	Total travel time	216.6 h	227.4 h	−1.533	−4.7%
	Total stops	789.5	2479.2	−3.997	−68.2%
	Queue length	773.7 m	1321.6 m	−4.006	−41.5%
Scenario 3	Total travel time	187.0 h	201.5 h	−1.959	−7.2%
	Total stops	538.5	2738.7	−4.096	−80.3%
	Queue length	447.0 m	1039.3 m	−4.689	−57.0%
Scenario 4	Total travel time	196.0 h	192.0 h	1.040	+2.1%
	Total stops	193.6	811.9	−2.956	−76.2%
	Queue length	523.2 m	819.1 m	−3.909	−36.1%
Scenario 5	Total travel time	160.8 h	153.9 h	1.814	+4.5%
	Total stops	160.0	666.1	−2.158	−76.0%
	Queue length	196.2 m	443.9 m	−2.819	−55.8%

. Statistically significant differences are underlined. Note that travel time values from macroscopic and microscopic simulations cannot be directly compared, since macroscopic simulations only provide lane averaged total travel time.

Table 3 shows significant reduction in total number of stops and queue lengths in VSL control case for all scenarios, which are the most obvious indicators of shock waves. With VSL control, the total number of stops is reduced by 68.2% to 80.3% and queue length by 36.1% to 55.8%. This indicates significant reduction in stop-and-go waves, which are less frequent and/or less intense and cause fewer oscillations in velocity.

The effect of VSL control on total travel time in microscopic simulations is similar to the effect in macroscopic simulations. While VSL control in Scenario 1 contributes to a significant reduction in travel time by 17.2%, the total travel time in VSL controlled case in scenario 5 is extended by 4.5%. Scenario 5 is the scenario with no wide moving jam but with several instances of synchronized flow. VSL control in this scenario leads to more fluent traffic but with slight travel time extension, therefore we can conclude that synchronized flow phenomena has a substantial impact and should be accounted for in VSL control strategies.

At this point, it is reasonable to discuss the selection of parameter optimal density ${\rho }_{opt}={\rho }_{min}^{syn}$. Our selection means that we practically do not tolerate any density above the lowest synchronized flow threshold. This means zero tolerance for high velocity oscillations but can in general lead to total travel time extension. Selection of higher value would allow for more oscillations (which represent dangerous traffic situations), but could potentially mean lower travel time, since the actual onset of traffic breakdown is of a probabilistic nature in the synchronized flow region.

5. Conclusions

In this paper, we proposed a numerical VSL control method that aims at reducing stop-and-go waves caused by one or more stationary bottlenecks. The method presented in this paper consists of a macroscopic simulation, introduction of a cost functional and numerical optimization with DE algorithm.

The method represents a novel approach to VSL control using a numerical scheme for solving systems of PDEs of a continuum traffic flow model without previous discretization of PDEs, which is usually used to avoid differential formulation. This is enabled by a DE algorithm that does not require analytical calculation of gradient and enables using a cost functional, that is not directly dependent on the model formulation but solely on its solution. Moreover, the presented DE optimization algorithm and cost functional can in general be used for any continuum model (with incorporated speed limits) that can be solved with any numerical scheme. This means reducing the gap between the wide set of existent continuum models and the limited set of models applied to VSL control strategies. While the basis of model predictive control is a deterministic traffic flow model, we believe it is of great importance to be aware of the synchronized flow region, where the breakdown phenomena appear. We employ a modified switching curve model, that is in compliance to three phase traffic flow theory, which enables integration of more empirical features of traffic breakdown into model predictive control.

The application of the proposed VSL method is presented on a motorway section in Slovenia with multiple bottlenecks. The efficiency of the proposed VSL control is tested with microsimulations of the test sections for five traffic volume scenarios, each with 30 random seeds. The results provide insight into the power of influencing traffic flow with VSLs and raise some interesting consideration regarding control policies with respect to the actual nature of traffic breakdown in the synchronized flow region. While appropriate VSL control results in greatly reducing oscillations in velocity in the synchronized flow region, the total travel time does not necessarily act in the same manner. This result should certainly be considered when defining any other cost function in future control strategies. The presented research provides some insight on the trade-off between travel time and speed oscillations by VSL control, which will be investigated in detail in future research.

The case study with two VSL locations shows fast computation, which enables real-time application. The method is not limited to two VSL locations, since DE algorithm operates with vectors of optimization variables (the vector can have n components that represent VSLs in n locations). However, the number of VSL locations affects the computation time, since greater number of VSL locations means a greater number of possible solution candidates with more components. Since the code enables multi-core optimization, the computation speed is greatly dependent on the type of the processor and number of cores. The scalability of the solution method for larger networks with several VSL locations is a subject for future research. In this context, we should also mention the case where a queue spills back past the location of VSL. The method works by setting VSLs only until the queue reaches the VSL location. To achieve meaningful and comprehensive control after this point, this case should be additionally controlled by a VSL outside of the section, i.e., this case also requires a study of a larger network. It should be noted that since the proposed method presents a completely novel approach, there are several difficulties and drawbacks. First is a higher computational burden, which can be managed with the appropriate optimization method, i.e., VSLs for the test sections can be computed in real-time. Appropriate modelling of the section and its possible bottlenecks can take same additional mathematical effort and requires more attention than usual strategies in order to be consistent with the system of PDEs. While we successfully managed the two mentioned difficulties, there exists a major challenge. In this study, we used model predictions for the test section by prescribing that initial conditions of a numerical interval are equal to the numerical solution of the previous interval. However, the initial condition of the PDEs system in each step enables prescribing actual field data, i.e., a feedback loop could be integrated. However, this traffic flow estimation should be performed with careful mathematic consideration and represents a subject for future research.