Cost and Risk Prediction in Road Transportation of Hazmat by ANFIS Trained with PSO, FA, HBBO and ICA

The risk of accident in road transportation of hazardous material can occur due to various mechanisms, and it is significantly dependent on human factors, environmental conditions, management factors and the equipment Failure. Different accident scenarios may result from a combination of events following road accident depending on the type of transported substance (explosion, fire and the release of toxic material). The damage resulting from road accidents gives rise to a number of consequences, which are classified in three main categories (environment, social and economic) see Figure 3.

In order to estimate the risk of accident aggregate with its cost this paper proposes a method based on ANFIS model, which is based on three main steps see Figure 1.

In the first phase the criteria (inputs) that influences the value of the CR are identified. The inputs and their descriptions are illustrated in Table 1. After that the road network on which the value of the input parameters and the output (CR) are given on each branch that separate the node (client, depot station ...) must be selected.

1.jpg

Figure 1. Main steps of ANFIS for determinng the CR value

2.jpg

Figure 2. ANFIS based PSO, FA, ICA and HBBO

Table 1. Criteria for defining the CR values on a transport network of city roads [3]

Inputs	Description of inputs
The carrier’s operating costs (X 1 ∈ [0,7500])	The cost of operator carrier is mainly computed according to the value of time, distance, cost of fuel etc. In this study carrier, operating costs is given based on the length of the road [13].
Fire Fighters response (X 2 ∈ [0,15])	Fire fighters response depend entirely on Traveled distance (Lij) and the costs of the transportation c (fuel consummation for example), which is given by: Lij* c. Fire fighters response is the necessary time taken by of emergency services to arrive at the location of accident. The average reaction time is expressed in this study in minutes, which is relative to the distance that separate emergency services center to the hot spot.
Losses related to the environment (X 3 ∈ [0,5])	This is the number of environment's critical point (Headwaters Rivers, Agricultural Urban, Forest fragment) that can be affected by the accident. These stragicals points that are within 800m from the accident location are taken into account [14].
Risk of an accident (X 4 ∈ [0,1])	Risk of an Accident is determined according to of the number of accident (f 1) recorded in database of the concerned road over the last 10 years. the number of branches (f 2), the percentage of heavy cars in the traffic (f 3) and the level of signaling (f 4). The total Risk of an Accident (X 4) is computed based on this equation X4 = 0.3 f1/ f1max + 0.2 f2/ f2max + 0.2 f3/ f3max + 0.3 f4/ f4max. the value of f4 is given by this following scale: 1 – traffic rules, 2 – regulation by means of traffic signs, 3 – regulation by means of light signals.
The gravity related to an accident (X 5 ∈ [0,9500])	Is given according the affected population who lives around the accident's zone. The control zone is determined as mentioned in the factor X3.
Losses related to the stratigical points (X 6 ∈ [0,10])	Risk related with stragicals point is the number of important infrastructures presented in the accident’s locations (railways, electric lines, industrial, commercial and transport facilities, schools, hospitals).
Risk of terror attack/hijack (X 7 ∈ [1,5])	Risks of terror attack is the consequences that caused by an eventual attack of terrorist in the control zone and its value is determent by adopting this following scale: 1- very low, 2 – low, 3 – moderate, 4 – height 5 – very height.

3.jpg

Figure 3. Hierarichal framwork of causes and consequances of raod accident

In the second step the type of the ANFIS and the membership function of each input and output must be selected then the first model is generated in order to map from the inputs to the output by the means of fuzification, product, normalization and defuzification layers (section 2 highlight how to map from the input parameters to the CR).

In order to optimize the performance of the proposed method and avoid the inclusion of expert for the designing parameters of the ANFIS four algorithm from meta-heuristic optimization (PSO, FA, ICA and HBBO) are used to tune the main parameters of ANFIS (the mean and variance of Gaussian membership function and the constant pi, qi and ri. See Figure 2). Finally after the termination criteria of the tuning process (the number of iterations) are satisfied then ANFIS model is ready to be used.

4.1 Particle swarm optimization (PSO)

PSO is nature-inspired methods, which is included in swarm-based techniques. It is based on a simple mathematical model, introduced by Kennedy and Eberhart in 1995 [17]. PSO imitate the social behavior of swarm when fling and seeking to reach a destination. It employs a number of particles (candidate solutions) which explore a given search space to find the best solution (i.e. the global best solution). Meanwhile, they update the best location (personal best solution) of each candidate in their paths. In other words, each particle moves through the decision space based on the individual and global best positions. In this algorithm the position (ANFIS parameters) of each particle are updated based on Eq. (9)-(10).

$\mathrm{V}_{\mathrm{i},}^{\mathrm{t}+1}=\mathrm{w} \mathrm{V}_{\mathrm{i}}^{\mathrm{t}}+\mathrm{c}_1 \mathrm{r}_{1 \mathrm{i}}^{\mathrm{t}}\left(\mathrm{p}_{\mathrm{i}, \text { best }}^{\mathrm{t}}-\mathrm{x}_{\mathrm{i}}^{\mathrm{t}}\right)+\mathrm{c}_2 \mathrm{r}_{2 \mathrm{i}}^{\mathrm{t}}\left(\mathrm{G}_{\text {best }}^{\mathrm{t}}-\mathrm{x}_{\mathrm{i}}^{\mathrm{t}}\right)$          (9)

$\mathrm{X}_{\mathrm{i},}^{\mathrm{t+1}}=\mathrm{V}_{\mathrm{i},}^{\mathrm{t}+1}+\mathrm{x}_{\mathrm{i}, \mathrm{j}}^{\mathrm{t}}$          (10)

where: $\mathrm{V}_{\mathrm{i},}^{\mathrm{t}+1}$, $\mathrm{V}_{\mathrm{i}}^{\mathrm{t}}$, indicate the velocity of each particles at the current iteration and previous iteration, $\mathrm{X}_{\mathrm{i}}^{\mathrm{t}+1}$ and $\mathrm{x}_{\mathrm{i}}^{\mathrm{t}}$ denotes the position of particles in current and previous iteration respectively, $\mathrm{p}_{\mathrm{i}, \mathrm{best}}^{\mathrm{t}}$ and $\mathrm{G}_{\text {best }}^{\mathrm{t}}$ represents the personal and the global best of each candidate and the hole swarm, w, $\mathrm{c}_{1,2}$ and r_1,2 are respectively the Inertia Weight, personal and global learning coefficient and random number.

4.2 Firefly algorithm (FA)

FA is bio-inspired meta-heuristic algorithm, which belong to swarm optimization technique. It was firstly developed by Yang in late 2007 and 2008 [18]. This algorithm is inspired by flashing pattern and behavior of fireflies at night. It engages a number of fireflies (candidate solutions) which fly around in the search space and move toward the attractive one based on its brightness regardless of their sexes, which constitute the best solution (i.e. the optimal position). It is worthy to notice that the dynamics of fireflies can be reduced to the dynamics of Particle swarm optimization, Harmony Search, simulating Annealing and differential evolution, which gains the ability of the automatic subdivision. The position of fireflies is computed according to the Eq. (11).

$\mathrm{x}(\mathrm{t}+1)_{\mathrm{i}}^{\mathrm{k}}=\mathrm{x}(\mathrm{t})_{\mathrm{i}}^{\mathrm{k}}+\beta_0 \exp \left(-\gamma \mathrm{r}_{\mathrm{ij}}^2\right)\left(\mathrm{x}(\mathrm{t})_{\mathrm{j}}^{\mathrm{k}}-\mathrm{x}(\mathrm{t})_{\mathrm{i}}^{\mathrm{k}}\right)+\alpha\left(\mathrm{r}-\frac{1}{2}\right)$          (11)

with

$\beta(\mathrm{r})=\beta_0 \exp \left(-\gamma \mathrm{r}_{\mathrm{ij}}^{\mathrm{m}}\right) ; \mathrm{m}<1$         (12)

$r_{i j}=\left\|x_j-x_i\right\|=\sqrt{\sum_{k=1}^d\left(x_j^k-x_i^k\right)^2}$          (13)

where: $\beta(\mathrm{r})$ indicates the attractiveness of firefly, $\beta_0$ represents the attractiveness of firefly at $\mathrm{r}=0$, Coefficient of light absorption, $\mathrm{x}_{\mathrm{i}}^{\mathrm{k}}$ is the position of i^th firefly at the k^th dimension, $\mathbf{x}_{\mathrm{j}}^{\mathrm{k}}$ denotes the position of the j^th firefly at the k^th dimension, $\mathrm{r}_{\mathrm{ij}}$ is the distance between $\mathrm{x}_{\mathrm{i}}$ and $\mathrm{x}_{\mathrm{j}}$, t iteration, r random number, α  randomization parameter.

4.3 Human based behavior optimization (HBBO)

HBBO is a new meta-heuristic optimization proposed by Ahmadi in the late of 2016 [19]. It imitates the behavior of humans in getting success in their life. Humans looks for the success in several fields (music, paint, engineering…), and after the domain of interest is opted based on their purposes and passion they moves to the expert to improve and excel their ability and hence gains more success.

HBBO in its first phase initialize set of population (candidate) and randomly spread them in deferent fields and the number of initial people in each field is obtained by:

N.Ind ${ }_i=$ round $\left\{\frac{N_{\text {pop }}}{N_{\text {field }}}\right\}$          (14)

where, $\mathrm{N}_{\text {pop }}$ and $\mathrm{N}_{\text {field }}$ are respectively the number candidates and fields.

In the second phase (education phase) this algorithm describe how the population move toward expert in their deferent fields. The experts are the candidates who possess the best solution. Each candidate will converge toward its relevant expert based on the spherical coordinate system which is given by:

$\mathrm{x}_{\mathrm{ind}}^{\mathrm{t}+1}=\mathrm{x}_{\mathrm{ind}}^{\mathrm{t}}+\mathrm{r} \cdot\left(\mathrm{k}_2-\mathrm{k}_1\right) \cdot\left(\mathrm{x}_{\mathrm{exp}}^{\mathrm{t}}-\mathrm{x}_{\mathrm{ind}}^{\mathrm{t}}\right)$          (15)

With $\mathbf{x}_{\text {ind }}^{\mathrm{t}+1}$ and $\mathrm{x}_{\mathrm{ind}}^{\mathrm{t}}$ denotes respectively the current and previous position of individuals, $\mathrm{x}_{\mathrm{exp}}^{\mathrm{t}}$ represent the position their expert, r is random variable, $\mathrm{k}_1$ and $\mathrm{k}_2$ are a constant.

After that in the consultation phase individuals may find within society an advisor who can change their minds and therefore they may changes their field. This consultation can be effective if we get a furthermore best solution consequently advisors changes the candidate of solution, otherwise nothing will be changed. The number of decision variable (parameters of ANFIS) that will be changed during this phase is computed in Eq. (16).

$\mathrm{N}_{\mathrm{c}}=\operatorname{round}\left(\sigma . \mathrm{N}_{\mathrm{var}}\right)$          (16)

where, $\mathrm{N}_{\mathrm{var}}$ is the total number of decision variable and $\sigma$ is a constant.

In the fourth phase (field changing probability) HBBO firstly ranks and sort the expert individuals based on their fitness score then calculate the changing probability by:

$\mathrm{p}_{\mathrm{i}}=\frac{\mathrm{O}_{\mathrm{i}}}{\mathrm{N}_{\text {field }}+1}$          (17)

where, $\mathrm{p}_{\mathrm{i}}$ and $\mathrm{O}_{\mathrm{i}}$ are ranking probability and the index of the i^th sorted field, respectively. After that a random variable (rand) between 0 and 1 is generated, this comparison is performed:

if rand $\leq \mathrm{p}_{\mathrm{i}} \rightarrow$ field changing occurs           (18)

In this phase, a selection probability for every person is computed by:

$\mathrm{P}_{.} \mathrm{S}_{\mathrm{j}}=\left|\frac{\mathrm{f}\left(\text { individual }_{\mathrm{j}}\right)}{\sum_{\mathrm{k}=1}^{\text {Nind }} \mathrm{f}\left(\text { individual }_{\mathrm{k}}\right)}\right|$           (19)

With $\mathrm{f}\left(\right.$individual$_{\mathrm{j}}$) is the fitness score of the j^th people and Nind is the number of people in the selected field. Finally, by adopting the roulette wheel selection technique [20], a candidate will be selected and moves to a random field.

Finally, like others algorithms HBBO perform the Finalization step (check the termination criteria).

4.4 Imperialist competitive algorithm (ICA)

ICA is a novel evolutionary algorithm originally developed by Atashpaz-Gargari and Lucas [21]. It mimics the imperialistic competition. The essence of this algorithm can be attained in 5 main steps (empire initialization, assimilation, revolution, imperialist competition and finalization).

In the first phase like other algorithms ICA generate randomly an initial population (countries), the countries with best score are called imperialists and the others are known as colonies, then spread them among other imperialist based on their relevant power, and by performing roulette wheel selection [20, 22] an initial empire will be formed. The number of imperialists is obtained by:

$\mathrm{N}_{\text {imp }}=\operatorname{round}\left\{\frac{\mathrm{N}_{\text {pop }}}{\mathrm{N}_{\mathrm{emp}}}\right\}$          (20)

With $\mathrm{N}_{\text {pop}}$ and $\mathrm{N}_{\mathrm{emp}}$ indicate the number of population and empires.

After that in the second phase (assimilation phase) all colonies moves to reach the position of their imperialist based on uniform distribution written in Eq. (21).

$\mathrm{x} \sim \mathrm{U}\left(\mathrm{x}_{\mathrm{col}}^{\mathrm{t}}, \beta \times \mathrm{d}\right)$          (21)

where, $\mathrm{x}_{\mathrm{col}}^{\mathrm{t}}$, $\beta$ and $\mathrm{d}$ denotes respectively the current position of colony, is assimilation coefficient and d is the distance between colony and its imperialist.

In the revolution phase (third phase) each country among empires seeks to increase their power by changing the position of their decision variable. During this phase the number of variable that will be changed are determined according to Eq. (16). Then the selected variables of imperialists are updated based on normal distribution described in Eq. (23), while for colonies a predefined probability (pc) is firstly considered after that a random variable (r) between 0 and 1 is generated then a comparison depicted in Eq. (22) is performed and if this comparison is verified they update their position based on Eq. (23).

if $\mathrm{r} \leq \mathrm{pc} \rightarrow$ rovolution process occurs          (22)

$\theta \sim \mathrm{N}\left(\mathrm{x}_{\mathrm{imp}, \mathrm{col} \quad}^{\mathrm{t}}, \gamma \cdot\left(\mathrm{x}_{\max }-\mathrm{x}_{\min }\right)\right)$          (23)

where, $\theta$ indicate a new position, $\mathbf{x}_{\mathrm{imp}, \mathrm{col}}^{\mathrm{t}}$ is the current position of imperialist and colony, $\gamma$ is revolution rate, $\mathrm{x}_{\max }$ and $\mathbf{x}_{\min }$ denotes the upper and lower bond of the decisions variable. Here it is worthy to notice that during this phase a colony may reach better position than its imperialist therefore it became a new leader of its relevant empire.

In order to perform the fourth step the total power of each empire is determined at each iteration by Eq. (24).

$\mathrm{TC}_{\mathrm{n}}=\operatorname{cost}\left(\right.$imperialist$\left._{\mathrm{n}}\right)+\xi \operatorname{mean}\left\{\operatorname{cost}\left(\right.\right.$colonies$\left.\left._{\mathrm{n}}\right)\right\}$          (24)

where, $\mathrm{TC}_n$ denote the total cost of an nth empire and ξ represent a positive number which is considered to be less than 1.

In the fourth phase (empire competition) each empire seeks to takes possession of the weakest colonies of other empire. The more colonies in empire the strongest it becomes the fewer colonies in empire the weaker it becomes. In this phase a weakest colony is firstly selected then the probability selection is calculated by:

$\mathrm{p}_{\mathrm{jemp}}=\exp \left(-\alpha \frac{\mathrm{TC}_{\mathrm{jemp}}}{\sum_{\mathrm{k}=1}^{\mathrm{Nemp}} \quad \mathrm{TC}_{\mathrm{k}}}\right)$          (25)

With $\mathrm{p}_{\text {jemp }}$ is the probability selection of the j^th empire, $\alpha$ is the selection pressure and $\mathrm{TC}_{\text {jemp }}$ is a total power of j^th empire. Finally based on the roulette wheel selection the winner empire is determined hence it possess the selected colony. The empires that are not able to gain more colonies collapses and therefore became a colony of other empires.

In the last phase of these algorithm (finalization) the above-mentioned process from phase 2 are repeated until the termination criteria are verified.