A Bi-Level Programming Model of Liquefied Petroleum Gas Transportation Operation for Urban Road Network by Period-Security

As a clean energy, Liquefied Petroleum Gas (LPG) is consistent with the coordinated and sustainable development of both the economy and environment. However, LPG is a hazardous material (hazmat) and is thus always transported in cylinders by vehicles on urban road networks to meet varying demand. This transport can threaten the surrounding citizens, vehicles, and even the whole urban area. Hence, LPG transportation should be focused on maintaining its security while simultaneously minimizing shipping costs. When LPG is moved through an urban area, its threat level fluctuates with the network congestion level, which continually varies by different time periods. So, variation in the magnitude of the threat posed by LPG transportation causes additional changes in the safe-related cost as well as the shipping cost. This study aims to solve the problem of an LPG transportation operation on an urban road network according to congested periods; the solution is based on cutting its two types of cost. In general, we should choose an LPG transport period that results in a lower safety cost, however optimization of an LPG transportation operation must minimize both the safety cost and shipping cost. This paper presents the problem of LPG flow distribution and vehicle dispatch scheme by “period-security” to rationalize the LPG transport risk level. Firstly, the impedance function of LPG flow distribution was constructed with a focus on the safety cost in different periods. Meanwhile, a bi-level programming model was built, in which the upper mixed binary integer programming model aims to minimize the shipping cost and the lower model is a user-equilibrium model that is aimed at calculating the distribution of the LPG demands on the given lines and in feasible periods. Then, we designed a heuristic algorithm based on the Genetic Algorithm to solve the upper model and embedded the Frank-Wolfe Algorithm to get the optimal LPG flow distribution solution. Numerical examples are presented which validate that the LPG optimal operation can realize a minimal safety cost and the minimum shipping cost for three LPG demand values by considering the congestion situation.


Introduction
Liquefied Petroleum Gas (LPG) is a clean energy that, by nature, is less polluting, generates high heat, and burns cleaner with less carbon build-up and oil contamination [1].Hence, it has become a convenient, environmental, and sustainable fuel that is used in many fields and for many applications, such as non-ferrous metal industry, heating appliances, and cooking equipment.Nevertheless, LPG has dangerous features which mainly include its inflammability, explosiveness, and poisonousness [2][3][4][5], all of which may be high risk or may even cause accidents during storage and transportation.Although an LPG pipe network can transfer large amounts of gas and has been built in some cities, transporting LPG in cylinders by vehicles routed through an urban road network is still quite necessary for sporadic demands.According to the Chinese Provision CJ/T 35-2004 [6], LPG transportation should avoid the routes and periods with a high density of the population and vehicles.During the transportation process, potential LPG related risks and accident damage are augmented if there are many surrounding passengers and vehicles.So, choosing the optimal period to decrease the overall threat of the above occurrences is significant.An LPG cylinder supply station (LCSS) is in charge of gas storage, cylinder filling, and distribution and transportation to users (including transfer centers and households).Due to its hazardous features, LPG cylinder transportation must be secure and safe, which could lead to a much higher safety cost related fee that is covered by users, while the LCSS hopes to cut down its transportation cost.These two kinds of costs conflict with each other because, for example, the shortest line with the lowest transport cost may not be the safest one.Theoretically, all roads crossing the city should be available to these LPG cylinder carrying vehicles throughout the whole day, and the LCSS should thus be able to choose any route in any time period.However, in practice, some routes in given periods (for instance, the peak hours with congestion) are not suitable for LPG transport and generate higher safety cost, which could effectively make some routes in given periods unavailable.This is measured by period-security in this paper.Therefore, a LCSS must balance between transportation feasibility, depending on period-security, and its capacity for all periods in order to meet the fluctuating user demand.All of these considerations cause an LPG vehicle periodic and uneven operation of LPG vehicles.This paper examines the coupling relationship between safety cost and transportation cost for different period-security conditions and optimizes the LPG vehicle distribution scheme by periods.
In terms of LPG, there are a few studies in the literature that introduce its features and transportation organization.LPG's potential risks and accident damages have urgently compelled the creation of a series of directives in countries worldwide [7][8][9] which clarify the dangers, the methods for avoiding accidents, and even the rescue protocols for such accidents.Although some researchers (Godwin Glivin et al. [10], Ramchandra Bhandari et al. [11]) are seeking the substitution of LPG energy, such as biowaste and electricity, there are still many scientists that regard LPG as one of the most sustainable green energies.Laurencas Raslavičius et al. [12] evaluated LPG as a clean, relatively low cost, and abundant energy source that provides affordable fuel-efficient transportation, and they discussed different types of existing LPG fuel delivery systems and requirements for gas equipment.To enhance the safety of rail or truck transport of LPG, Swoveland C et al. [13] built a model to compute the expected net benefits by estimating costs, potential effects on LPG release rates and consequences, and the unit costs of incorporating the results into new and existing rail tank cars and tank trucks.To address LPG transportation accidents, Lã 3 Pez-Atamoros et al. [14] studied the main causes of accidents in Mexico and created a classification guide in order to harmonize criteria for the assessment of the risk index parameters.Bariha et al. [15] introduced two software programs-Area Locations of Hazardous Atmospheres and Process Hazard Analysis Software Tool-to model LPG accident scenarios and simulate the results of a fireball, jet flame radiation, and explosion overpressure.
Transporting LPG in the urban road network must comply with relevant safety policies and optimizing regulations.There are many literatures of city logistics and its sustainability.Meiling He et al. [16] and Nathanail E et al. [17] studied urban logistics and freight transportation and promoted the methods of keeping sustainable and livable urban development.Being a special type of freight, LPG transportation must comply with those basic principles of urban sustainability.Furthermore, LPG transport should integrate into the urban logistics system by obeying concrete directives when it moves across the urban road network.Special requirements of LPG transportation in urban network cover periods limitation, speed upper bound, non-LPG zones, sunshade in summer, and limitation of height and speed when passing channels, culverts, and interchange bridges [6].
Zhang [18] analyzed LPG transport condition and its rescuing methods to ensure sustainable and safe social development.
Although studies on LPG are limited, hazmat transport has been actively researched for a long time and can be extrapolated to LPG quite easily.Qiao Y et al. [19] estimated the accident frequency from hazardous materials transportation by utilizing publicly available databases and expert knowledge.Cozzani V et al. [20], Chakrabarti U K et al. [21], and Zhao et al. [22] processed a comprehensive risk assessment for hazmat routing and transporting.Bianco L et al. [23] studied a novel toll setting policy, and the tolls deterred carriers from using links that had a high total risk.Bagheri M et al. [24] compared the characteristics of hazmat transportation by rail and by road and demonstrated that rail transport reduced the risk, irrespective of the risk measure and the transport corridor.
The hazmat routing problem has been studied by many researchers for a long time.Caramia M et al. [25] selected k efficient paths with respect to the minimization of length, time (cost), and risk; in particular, the selection was made by choosing k representative paths with high spatial dissimilarity.Huang R L C B. [26] evaluated the risk of hazmat transportation by integrating Geographic Information Systems (GISs) and Genetic Algorithms (GAs) using the criteria of safety, costs, and security.Akgün V et al. [27] considered weather systems for routing for hazmat.Tarantilis C et al. [28] determined the routes that were used by a fleet of trucks serving a set of customers with a focus on population exposure risk mitigation via the production of truck routes by solving a variant of the Vehicle Routing Problem.Some researchers have focused on the vehicle routing problem with the time window aspect of the hazmat transportation problem [29][30][31][32][33]. Ren C X et al. [34] discussed the route-independent and route-dependent factors, including intrinsic road characteristics, meteorological conditions, traffic characteristics, and population density.Pradhananga et al. [35] studied the bi-objective optimization of the hazardous materials vehicle routing and scheduling problem with time windows, with objectives focused on minimizing the total scheduled travel time and the total risk for the transportation process.Androutsopoulos et al. [36] presented a formulation of the hazardous material distribution problem as a bi-objective time-dependent vehicle routing problem that focused on both the travel time and cost of the routes.Bronfman A et al. [37] introduced the maxi-sum hazmat routing problem, which maximizes the sum of the population-weighted distances from vulnerable centers to their closest point on the routes using multiple origin-to-destination (OD) pairs.Bronfman A et al. [38] proposed an approach that maximized the weighted distance between the route and its closest vulnerable centre in order to minimize the potential consequences at the most exposed center.Nembhard D A et al. [39] determined a hazmat path that maximizes a multi-attribute, non-order-preserving value function based on transportation cost and risk to the population.As is continuously emphasized, the value of the risk is an important factor affecting hazmat transport.Hosseini et al. [40] proposed a special value-at-risk (VaR) approach to denote rail hazmat shipments' features over a given railroad network with a limited number of train services with the transport risk measured by the minimal VaR.Ma et al. [41] built a multi-objective robust optimization model with the objective of minimizing not only transportation risk, however also transportation time, and a robust counterpart of the model was established by applying the Bertsimas-Sim robust optimization theory.Ma et al. [42] built a multi-objective optimization model to minimize total transportation energy consumption and transportation risk while studying the multi-depot vehicle routing problem for hazardous materials.Their result was to fix the set of customer points served by each depot and to determine their service order.
User Equilibrium (UE) was developed by Wardrop [43] and many researchers use it to address traffic distribution problems.Friesz T L et al. [44] proposed a new class of UE models with a variation inequality formulation of dynamic user equilibrium while simultaneously making route choice and departure time decisions.Huang H J et al. [45] considered a simultaneous route and departure time choice equilibrium assignment problem in a network with queues based on discrete-time and formulated it as an equivalent "zero-extreme value" minimization problem.Bell et al. [46] proposed a deterministic UE traffic assignment as an equivalent to the mixed-strategy Nash equilibrium of a n-player, non-cooperative game, where n network users seek their best routes and OD-specific demons penalize the network users maximally by failing links.Cantarella G E [47] presented a fixed-point formulation of a multi-mode multi-user equilibrium assignment with elastic demand.Prakash et al. [48] formulated a consistent, reliability-based user equilibrium problem with risk-averse travelers which incorporates endogenous link flow variances and covariances and flow-dependent and endogenous link travel time correlations.Carotenuto P et al. [49] studied the selection of paths that minimize the total risk of hazmat shipments while spreading the risk to the population in an equitable way.
The bi-level programming model with the UE principle is a hot topic worldwide and is used to solve some traffic problems as well.Maher M J et al. [50] combined trip matrix estimation with traffic signal optimization on congested road networks and built a bi-level programming model with stochastic user equilibrium assignment to describe a transport network.Gao Z et al. [51] designed an integrated method to maximize the reserve capacity of a road network with both the maximum possible increase in traffic demand and maximum road reserve capacity.Migdalas A [52], with the goal of transport planning, considered the interests and rights of both the network user and the public sector.Taslimi M et al. [53] discussed a bi-level hazmat transportation network design problem in a road network between specified OD points.Ceylan H et al. [54] addressed a congested area traffic control problem as an upper-level problem, while the user-equilibrium traffic assignment was treated as a lower-level problem.
According to the above literature review, we find that most researchers study the hazmat routing problem by combining risk equity with various planning models.However, they seldom consider how different traffic situations affect hazmat transportation at distinct time periods of the day, which should have a greater impact on the hazmat transportation safety level.To clarify the influence of different time periods on hazmat transport security, we initiated the "period-security" factor to measure the LPG transportation risk level and to maintain similar safety costs for all periods that were examined in this study.Additionally, the UE principle is used and the impedance function is constructed initiatively designed for distributing LPG flow.
This paper is organized as follows.Section 2 presents the LPG transport network (LTN) composition with period-security and a description of LTN congestion with the UE principle.Also, the impedance function is proposed here.Section 3 formulates the LPG flow distribution problem as a bi-level programming model.Section 4 proposes a genetic algorithm to solve the upper model and the Frank-Wolfe algorithm for the lower model.In Section 5, numerical examples are given to evaluate the bi-level model and the algorithm.Section 6 draws conclusions.

LPG Hazardous Features
An LPG is a colorless volatile liquid that is obtained by pressurizing cooling liquefaction of natural gas or oil in a refinery, and is usually supplied in pressurized steel cylinders [55].It has the following hazards: (i) explosiveness: if a cylinder is subjected to a fire of sufficient duration and intensity, it can undergo a boiling liquid expanding vapor explosion; (ii) inflammability: an LPG can flame easily due to its low flash point and combustion point; (iii) poisonousness: according to the National Institute for Occupational Safety and Health [56], at levels of 2000 ppm (3600 mg/m 3 ), an LPG is considered immediately dangerous to life and health by breathing it in, skin contact, and eye contact; (iv) expansibility: once there is any thermal condition (for example, high temperature during the noon in summer), an LPG can expand fast in cylinders; (v) volatility: as its boiling point is below room temperature, an LPG will evaporate quickly at normal temperatures and pressures [18].

Impacting Factors to LPG Accidents
There are a variety of LPG accidents including leakage, fire, explosion, and poisoning.The LPG's safe transportation condition can be easily damaged and many factors can cause an LPG accident when storing, filling, and transporting, such as excessive filling, thermal, mechanical collision, cylinder and vehicle defect, and operation failure [18].In order to decrease the overall LPG accident possibility, selecting the reasonable line and the suitable time period, which can avoid the above factors occurrence, makes significant sense when transporting an LPG across the urban area.

Threat of LPG Accident
Any LPG accident can possibly be out of manual control and can even lead to social harm.It threatens passenger and citizen health, causes damage of LPG and non-LPG vehicles, interferes with urban road networks, and, in particular, disturbs the urban environment security (water, air, and earth).These can seriously hinder environmental and urban sustainable development [56].

LPG Transport Network Description
An LCSS serves as an LPG storing compound, with vehicles delivering the LPG to certain users in its market area.The LCSS, its users, and all lines constitute an LTN.We denote S = {1, 2, • • • , m} as the set of users and K = {1, 2, • • • , n} as the set of lines.Figure 1 is an example of an LTN composed of one LCSS, 10 users, and n lines presented by different layers (on the left side of the figure).The directed lines in the LTN mean that vehicles send LPG cylinders to users while taking the empty ones back sequentially.Clearly, line 1 covers users of {6, 5, 4, 3, 2, 1, 7}, while line n serves {2, 3, 6, 5, 10, 9, 8}.The remainder of Figure 1 shows that the six concrete lines are fixed for the numerical examples that are given later in this paper.Notably, in the physical LTN, a given line k does not necessarily cover all users in S, however all chosen lines must be selected in such a way that all users' demands can be met.It is clear that each user can be served by a set of k lines, however their costs change with their distance.Therefore, choosing a line for a user is a combinational optimization problem.
x k p represents the LPG volume being transported by line k in period p.All 10 users need LPG from the LCSS; seven users (s = 1, 2, 3, 4, 5, 6, 7) can get it when line 1 is available during a certain period p.Clearly, three users (s = 8, 9, 10) cannot get service when x 1 p > 0, and their demands must be met by other lines and in other periods.Therefore, the LPG transportation operation of the LCSS is determined by both the line and period parameters simultaneously.

Period Definition
In an LTN, there are obvious fluctuations in LPG flow along a line at different time periods throughout a whole day.This may lead to the safety cost remaining high along certain lines in a certain period.
In this work, we divided 24 h into eight periods from 00:00 to 24:00.To reflect periods of the peak commuting hours in urban traffic and to simplify further calculation, each 3 h interval is regarded as one period represented by p , and the period set is fixed as the following: {1:00-4:00, 4:00-7:00, …, 22:00-1:00}.Clearly, the third and sixth periods are the congested peak periods, and their safety cost should be higher due to longer duration of exposure by non-LPG vehicles and a higher population joining the traffic network.Additionally, the temperature during noon periods is normally higher owing to strong sunlight in summer, therefore the relevant safety cost remains high.Hence, LPG transporting risk level of each line and each period is determined by the temperature, the population and traffic volume involved in the line, vehicles' running time, and so on.Safety cost is introduced and described in this paper.

LTN Congestion by Periods
A congested LTN is the subject of LPG transportation in this paper and differs from traditional traffic congestion.Although there are K lines available to both the LCSS and users throughout the whole day, transportation disutility still exists on these lines in a given time period for many reasons.These include the following main scenarios that are derived from the investigation of the LCSS and relevant management agencies: (i) non-LPG traffic flow interrupts LPG transport in the LTN; (ii) the transporting capacity of the LCSS in one period is fixed and could be less than the total demands of all users simultaneously; and (iii) line k may be unfavorable for the LCSS if the number of people

Period Definition
In an LTN, there are obvious fluctuations in LPG flow along a line at different time periods throughout a whole day.This may lead to the safety cost remaining high along certain lines in a certain period.
In this work, we divided 24 h into eight periods from 00:00 to 24:00.To reflect periods of the peak commuting hours in urban traffic and to simplify further calculation, each 3 h interval is regarded as one period represented by p, and the period set P = {1, 2, • • • , 8} is fixed as the following: {1:00-4:00, 4:00-7:00, . . ., 22:00-1:00}.Clearly, the third and sixth periods are the congested peak periods, and their safety cost should be higher due to longer duration of exposure by non-LPG vehicles and a higher population joining the traffic network.Additionally, the temperature during noon periods is normally higher owing to strong sunlight in summer, therefore the relevant safety cost remains high.Hence, LPG transporting risk level of each line and each period is determined by the temperature, the population and traffic volume involved in the line, vehicles' running time, and so on.Safety cost is introduced and described in this paper.

LTN Congestion by Periods
A congested LTN is the subject of LPG transportation in this paper and differs from traditional traffic congestion.Although there are K lines available to both the LCSS and users throughout the whole day, transportation disutility still exists on these lines in a given time period for many reasons.These include the following main scenarios that are derived from the investigation of the LCSS and relevant management agencies: (i) non-LPG traffic flow interrupts LPG transport in the LTN; (ii) the transporting capacity of the LCSS in one period is fixed and could be less than the total demands of all users simultaneously; and (iii) line k may be unfavorable for the LCSS if the number of people involved with this line and the public areas surrounding it are large, meaning that other lines would then be congested.The final situation is that the LPG flow distribution to each line in period p is similar and minimum, which obeys the UE principle.These three points may cause LPG demand to aggregate to line k in the exact transport period selected, and we refer to this phenomenon as "LTN congestion".
Therefore, the choice of an LPG transport line in period p can be reasonably assumed to obey the UE principle.The disutility of available transport lines can be represented by the safety cost of LPG transport corresponding to the chosen period.

LPG Line Choice: UE Principle
The safety cost of an LPG transported by line k is affected by three factors: (i) the population exposed to the LPG along line k in period p; (ii) the number of public areas (for instance, schools, hospitals, and military buildings) along line k; (iii) the non-LPG transport flow sharing line k in period p.In this paper, the safety cost of the LPG being transported on line k in period p is expressed by the parameter v k p , and its value is determined by methods of statistical data analysis and experts investigation.
For the congested LTN, we define the impedance function of the LPG flow distribution along line k in period p as follows: where v k p (0) is the comprehensive safety cost of line k in period p with zero flow; n k p is a given constant and denotes the maximum LPG flow volume along line k from the LCSS in period p; α is a parameter to be calibrated.It is clear that the above impedance value of LPG transportation reflects the overall safety cost, which is generated by both individual risk and social risk [57].This impedance value increases as the LPG flow climbs and complies with the common characteristics of the traditional impedance function.
When transporting LPG, the line and the period with the minimum impedance value are always chosen by the LCSS, which may cause "LTN congestion" and, naturally, this impedance level would increase synchronously.Then, this could cause the LCSS to abandon this choice and try others.After this cycling process, the final result is the UE status, that is, the impedance in each period on each line is similar and remains at the minimum value.This UE status, with optimal impedance value of each line and each period, can ensure the minimal "LTN congestion" level and, thus, can benefit the environmental and social justice in the urban area.Therefore, our first step is to describe the LPG flow distribution problem with the UE principle.

Problem Formulation and Hypotheses
In different periods of one day, the safety cost of LPG transportation fluctuates with the period-security level of the urban road network.The optimal LPG flow distribution should assure the minimum total safety cost.Meanwhile, the LCSS searches for the LPG vehicle distribution schemes with the lowest transport cost.To deal with these two conflicting costs, a bi-level programming model was built in this study, with the lower level model ensuring the minimal safety cost and the upper level model seeking the minimum transport cost.
In order to simplify the problem, this study assumed the following: (1) Each line in the urban road network is a closed cycle, and vehicles transfer LPG cylinders sequentially along the line after leaving the LCSS and ultimately return back to the LCSS.(2) The third and sixth periods are unavailable for LPG transport because of the higher risk level and lower period-security during these peak hours.

Decision Variables and Parameters
The decision variables and parameters in this paper are listed below.
x k p : the lower decision; variable presents the LPG flow volume being transported by line k in period p, ∀p ∈ P, k ∈ K; g p,k s : the intermediate variable; is the total LPG flow volume moving to s user by line k in period p, ∀p ∈ P, k ∈ K, s ∈ S; y k p : the upper decision variable; presents the number of vehicles carrying LPG cylinders from the LCSS to users by line k in period p, ∀p ∈ P, k ∈ K; q s : the LPG demand of user s, ∀s ∈ S; τ k p : the LPG transportation cost per vehicle by line k in period p, ∀p ∈ P, k ∈ K; β: the maximum number of LPG cylinders each vehicle can carry; a p : the maximum number of vehicles departing from the LTN in period p, ∀p ∈ P; 1, if the LPG vehicle connecting user s by link k in period p 0, otherwise.∀p ∈ P, k ∈ K, s ∈ S.

Lower Programming Model
The LCSS always chooses the period with the minimum impedance level, and this may cause the level to increase; a higher impedance then leads to the period being abandoned as feedback.This kind of impedance fluctuation with dynamic selection would finally result in UE of the LPG flow distribution.For example, all k lines have impedance values that are similar to each other as long as there is any LPG flow in period p.Let O k be the minimum impedance of line k, then, the UE principle can be expressed by the following formula: It is obvious that the UE principle (2) is equivalent to the following optimized model of the even LPG flow distribution function: Constraint (4) meets the LPG flow conservation requirement, that is, all the LPG transport demands are distributed to all p periods on k lines.Constraint (5) ensures that the LPG flow value from the LCSS to s user is positive.Constraint (6) shows the flow composition on line k in period p.

Upper Programming Model
There are different available choices of LPG transportation in period p. Different choices may cause diverse running costs.The upper model objective is to minimize the total LPG transportation cost for the LCSS, and its function is expressed as follows: min subject to y k p ≥ 0 , and is integer (10) Constraint (8) means that the total carrying capacity of vehicles on line k in period p is not less than the distributed flow volume that is obtained from the lower model.Constraint (9) represents the capacity of the vehicles dispatched by the LCSS.Constraint (10) ensures that the number of vehicles on line k in period p is a positive integer.

Algorithm Design
The bi-level programming model is NP-hard (Bard J et al. [58]) and several generally used solution methods, such as KKT-based (Kara B Y et al. [59], Bianco et al. [60]) and duality-based optimization approaches (Marcotte et al. [61]), could possibly lead to locally optimal solution (Vicente L et al. [62]).The Genetic Algorithm (GA) is characterized by global superiority and dramatic convergence, and its essence is the simulation of nature's survival of the fittest (Changxi Ma et al. [63]).Additionally, the lower problem is a UE model, and the Frank-Wolfe algorithm is widely regarded as an efficient method.So, in this study, a heuristic genetic algorithm embedded with the Frank-Wolfe algorithm was designed to solve this bi-level programming model.

Chromosome Coding
Chromosome integer coding is used in GA, which includes eight gene segments expressing eight transporting periods of the LTN during the whole day.Besides, each gene segment has six gene positions presenting the vehicle number on six LPG transporting lines.Therefore, each gene segment owns 48 (8×6) positions and can be described by y k p |k ∈ K, p ∈ P , say, one chromosome is a possible LPG transport operation with the line choice and the vehicle number.Every gene position is assigned an integer value that it belongs to 0, a p that represents the vehicle number in the corresponding LPG transport operation.The upper limit a p is the maximum vehicle number that the LCSS can dispatch in period p.The chromosome structure is shown in Figure 2.

Upper Programming Model
There are different available choices of LPG transportation in period p. Different choices may cause diverse running costs.The upper model objective is to minimize the total LPG transportation cost for the LCSS, and its function is expressed as follows: , and is integer (10) Constraint ( 8) means that the total carrying capacity of vehicles on line k in period p is not less than the distributed flow volume that is obtained from the lower model.Constraint (9) represents the capacity of the vehicles dispatched by the LCSS.Constraint (10) ensures that the number of vehicles on line k in period p is a positive integer.

Algorithm Design
The bi-level programming model is NP-hard (Bard J et al. [58]) and several generally used solution methods, such as KKT-based (Kara B Y et al. [59], Bianco et al. [60]) and duality-based optimization approaches (Marcotte et al. [61]), could possibly lead to locally optimal solution (Vicente L et al. [62]).The Genetic Algorithm (GA) is characterized by global superiority and dramatic convergence, and its essence is the simulation of nature's survival of the fittest (Changxi Ma et al. [63]).Additionally, the lower problem is a UE model, and the Frank-Wolfe algorithm is widely regarded as an efficient method.So, in this study, a heuristic genetic algorithm embedded with the Frank-Wolfe algorithm was designed to solve this bi-level programming model.

Fitness Function
The fitness function is used to evaluate each chromosome's quality and is the basis to get the optimum result.The fitness function is calculated by the following formula.

Fitness Function
The fitness function is used to evaluate each chromosome's quality and is the basis to get the optimum result.The fitness function is calculated by the following formula.
Z is the objective value of any chromosome in the present generation, while Z max and Z min are the maximum and minimum ones, respectively.When Z = Z max , we obtain the lowest fitness 0 and the highest fitness 1 when Z = Z min .This can describe quite well the objective function in the upper model.

Genetic Manipulation
Conventionally, genetic manipulation includes the process of selection, crossover, and mutation, all of which create the new generation.

Algorithm Implementation
Step 1: Initialization Set the parameters including the maximum generation times Ge, crossover probability p c , and the mutation rate p m .According to the above chromosome coding rule and the LCSS capacity constraint (9), the feasible initial population is generated randomly with scale of popsize.Let the number of iteration ϕ = 0.

Step 2: UE flow distribution
Frank-Wolfe algorithm is introduced to calculate the lower programming model for every individual in the present generation.Consequently, we work out the accumulative LPG volume x k p being transported by line k in period p.This process is listed as below: ( (5) Update the impedance value of all periods by the following recursion formula: If the current chromosome does not meet constraint (9), then randomly choose a gene position from all non-zero gene positions in the corresponding chromosome segment.Set The algorithm procedure is shown in Figure 3.

Numerical Examples
Figure 1 shows an example of an LTN with one LCSS, 10 users, and six different lines for assessing the LPG transportation operation bi-level optimized model and the genetic algorithm embedded with the Frank-Wolfe algorithm.

Parameters of the Impedance Function
According to statistics, 15 min is needed to fill the cylinders that can be carried by one vehicle according to investigations of the LCSSs.Then, the number of maximum vehicles for each period is a p = 180 ÷ 15 = 12, ∀p ∈ P, which means 12 vehicles can depart from the LCSS to users, one by one, during period p. Suppose that the same types of vehicles are used by the LCSS to send LPG cylinders and that each vehicle has a carrying capacity of 100 cylinders.In addition, each cylinder can be filled by 15 kg of LPG according to "Rules for Package and Transportation of LPG(CG/T )", developed and issued by the Ministry of Housing and Urban-Rural Development of P. R. China.Hence, β = 100 × 15 = 1500(kg/vehicle).
According to the statistics of the LCSSs, the frequency of sending vehicles to line k in one period is determined by the number of users on the line and their demands, which are denoted by the corresponding percentage Ω k .So, the maximum volume of sending dispatched LPG along line k from the LCSS in period p can be calculated by the formula n k p = a p • β • Ω k .Also, according to the principles mentioned earlier, transporting LPG can cause the safety cost to change with different period choices, so the value of v k p can be fixed as well.Table 1 shows all parameter values for the impedance function calculation.According to formula (1), the value of α influences the impedance value of the LPG flow distribution directly, for example, compared with the safety cost, the bigger the value of α is, the more LPG congestion affects the final impedance.In this paper, α = 1.5, which better reflects the equal status of both components of the impedance function.

Parameters of the Lower Programming Model
Table 2 explains the parameter values for the lower programming model with the UE principle.According to the earlier analysis, q s represents the LPG demands of user s.To test the effectiveness of both the LPG transportation operation bi-level programming model and the algorithm, examples with three different values of the total LPG demand of 10 users are provided in this paper: 60,000 kg, 72,000 kg and 90,000 kg.To shorten the paper, here we just list the value of q s for the 60,000 kg demand.Furthermore, to begin the iteration of the algorithm, the initial value matrix ε p,k s is listed below, which shows the linking status of s users on k lines in p periods.

Parameters of the Upper Programming Model
The length of line k determines its transportation cost.Table 3 shows the value of τ k p which represents the vehicle unit cost of transporting LPG by line k in period p.The LTN in Figure 1 includes 15 edges and their lengths are specified in Figure 4. From this, we get the total length of k links and their transport cost which are listed in Table 4. Line 1

Parameters of the Upper Programming Model
The length of line k determines its transportation cost.Table 3 shows the value of k p τ which represents the vehicle unit cost of transporting LPG by line k in period p .The LTN in Figure 1 includes 15 edges and their lengths are specified in Figure 4. From this, we get the total length of k links and their transport cost which are listed in Table 4.   3 0 1 1 1 1 1 0 0 0 0  Line 4 1 1 1 0 1 1 1 0 1 1  Line 5 1 1 1 1 1 1 1 1 1 0  Line 6 0 1 1 0 1 1 0 1 1 1    (1 CNY = 0.1439 USD) 600 620 500 740 640 580

Comprehensive Analysis
The algorithm was realized on a computer configured with an Intel(R) Core(TM) 2 Duo CPU P7450 with 2.93 GHz and 4 GB, and the numerical example was solved by a program developed with Matlab2014a software and took 3.03 min.The population was set to a size of 200, the number of iterations was 100, the crossover ratio was 0.7, and the mutation ratio was 0.07.The objective function value of the upper programming model tended toward convergence after 38 generations and the optimized transportation cost appeared at the 42nd iteration.Figure 5 shows the changes in minimum transportation cost in the upper model with respect to the number of generations.

Comprehensive Analysis
The algorithm was realized on a computer configured with an Intel(R) Core(TM) 2 Duo CPU P7450 with 2.93 GHz and 4 GB, and the numerical example was solved by a program developed with Matlab2014a software and took 3.03 min.The population was set to a size of 200, the number of iterations was 100, the crossover ratio was 0.7, and the mutation ratio was 0.07.The objective function value of the upper programming model tended toward convergence after 38 generations and the optimized transportation cost appeared at the 42nd iteration.Figure 5 shows the changes in minimum transportation cost in the upper model with respect to the number of generations.
For the optimized LPG distribution scheme of the lower model, with the optimal 42nd generation in the upper model, the solution accuracy of the Frank-Wolfe algorithm was 1asn -4 after 15 iterations.The iteration steps changed from 0.7064 to 9.01×10 −5 being shown in Figure 6.

Analysis of the Optimal LPG Transportation Operation Results with 60,000 kg Demand
According to the optimized LPG transportation operation results with a 60,000 kg LPG demand, there are 48 vehicles dispatched to transport the LPG cylinders from the LCSS to all 10 users along six lines in eight distinct periods.This transportation operation can meet all users' LPG demands and its details are shown in Figure 7.For the optimized LPG distribution scheme of the lower model, with the optimal 42nd generation in the upper model, the solution accuracy of the Frank-Wolfe algorithm was 1asn -4 after 15 iterations.The iteration steps changed from 0.7064 to 9.01×10 −5 being shown in Figure 6.For the optimized LPG distribution scheme of the lower model, with the optimal 42nd generation in the upper model, the solution accuracy of the Frank-Wolfe algorithm was 1asn -4 after 15 iterations.The iteration steps changed from 0.7064 to 9.01×10 −5 being shown in Figure 6.   5).In the above formula, due to the peak hour congestion, only six periods are clearly available for LPG transport out of all eight periods of the whole day according to the analysis in Section 2.2.Therefore, when the total LPG demand is 90,000 kg and the transporting capacity of the LCSS is close to the maximum, 15 vehicles are needed to fulfill all transporting tasks, however this number exceeds the maximum vehicle capacity (12 vehicles) of the LCSS in period p .This illustrates that the UE solution for a 90,000 kg demand can be achieved only if we change other relevant parameters: for example, by enhancing the maximum vehicle capacity p a of the LCSS in each period or increasing the carrying capacity of each vehicle β .

$
For the LPG demand of 90,000 kg, 66 vehicles are needed (seen in Figure 9).It is quite necessary to clarify that the maximum transporting capacity of the LCSS one day can be calculated by the following formula: In the above formula, due to the peak hour congestion, only six periods are clearly available for LPG transport out of all eight periods of the whole day according to the analysis in Section 2.2.Therefore, when the total LPG demand is 90,000 kg and the transporting capacity of the LCSS is close to the maximum, 15 vehicles are needed to fulfill all transporting tasks, however this number exceeds the maximum vehicle capacity (12 vehicles) of the LCSS in period p.This illustrates that the UE solution for a 90,000 kg demand can be achieved only if we change other relevant parameters: for example, by enhancing the maximum vehicle capacity a p of the LCSS in each period or increasing the carrying capacity of each vehicle β.
For the LPG demand of 90,000 kg, 66 vehicles are needed (seen in Figure 9).In the above formula, due to the peak hour congestion, only six periods are clearly available for LPG transport out of all eight periods of the whole day according to the analysis in Section 2.2.Therefore, when the total LPG demand is 90,000 kg and the transporting capacity of the LCSS is close to the maximum, 15 vehicles are needed to fulfill all transporting tasks, however this number exceeds the maximum vehicle capacity (12 vehicles) of the LCSS in period p .This illustrates that the UE solution for a 90,000 kg demand can be achieved only if we change other relevant parameters: for example, by enhancing the maximum vehicle capacity p a of the LCSS in each period or increasing the carrying capacity of each vehicle β .
For the LPG demand of 90,000 kg, 66 vehicles are needed (seen in Figure 9).7.
In addition, the more LPG transporting demand is, the more congested each period will be.This can be expressed by the comparison of k p x value with 60,000 kg, 72,000 kg, and 90,000 kg. Figure 10 shows the difference of k p x of these three LPG flow distribution situations.7.In addition, the more LPG transporting demand is, the more congested each period will be.This can be expressed by the comparison of x k p value with 60,000 kg, 72,000 kg, and 90,000 kg. Figure 10 shows the difference of x k p of these three LPG flow distribution situations.To meet transportation demand and to maintain the minimum period impedance level, the LCSS distributes a suitable number of vehicles to users in each period and by each line.Clearly, the optimal LPG distribution scheme for a 90,000 kg demand cannot be realized with the UE principle until the relevant parameters or constraints are modified, which will be subject to further studies.

Conclusions
In this study, we built a bi-level programming model of an LPG transportation operation in an urban road network using the period-security factor.In the lower level model, we propose a novel LPG flow distribution problem that considers both the lines and period-security in the context of the UE principle.This problem is different from the conventional one because it incorporates UE to assure a similar impedance value for each line and each period.Therefore, the optimal LPG flow distribution with the minimum comprehensive safety cost can be worked out.In the upper level model, the optimal LPG vehicle distribution scheme assures the minimal transportation cost for the LCSS.
(1) For the LPG flow distribution with the UE principle, the following initiatives and conclusions are made: • To formulate the problem, we first clarify LPG hazards and threat.With the analysis of impacting factors to LPG accidents, we conclude that selecting the reasonable transportation lines and periods can benefit urban sustainability.

•
We represent a given LTN for which the LCSS makes line choices according to periodsecurity.

•
"LTN congestion" by period is defined, and the UE principle is introduced to describe LTN congestion characteristics.

•
An impedance function consisting of the game relationship between safety cost and the impacts of LTN congestion is proposed to reflect the LCSS's preference of lines and periodsecurity choices, and to describe the spirit of environmental, urban, and social justice.To meet transportation demand and to maintain the minimum period impedance level, the LCSS distributes a suitable number of vehicles to users in each period and by each line.Clearly, the optimal LPG distribution scheme for a 90,000 kg demand cannot be realized with the UE principle until the relevant parameters or constraints are modified, which will be subject to further studies.

Conclusions
In this study, we built a bi-level programming model of an LPG transportation operation in an urban road network using the period-security factor.In the lower level model, we propose a novel LPG flow distribution problem that considers both the lines and period-security in the context of the UE principle.This problem is different from the conventional one because it incorporates UE to assure a similar impedance value for each line and each period.Therefore, the optimal LPG flow distribution with the minimum comprehensive safety cost can be worked out.In the upper level model, the optimal LPG vehicle distribution scheme assures the minimal transportation cost for the LCSS.
(1) For the LPG flow distribution with the UE principle, the following initiatives and conclusions are made: • To formulate the problem, we first clarify LPG hazards and threat.With the analysis of impacting factors to LPG accidents, we conclude that selecting the reasonable transportation lines and periods can benefit urban sustainability.

•
We represent a given LTN for which the LCSS makes line choices according to period-security.

•
"LTN congestion" by period is defined, and the UE principle is introduced to describe LTN congestion characteristics.

•
An impedance function consisting of the game relationship between safety cost and the impacts of LTN congestion is proposed to reflect the LCSS's preference of lines and period-security choices, and to describe the spirit of environmental, urban, and social justice.
(2) The conclusions for the bi-level programming model and the algorithm can be summarized by the following:

•
The LPG transportation operation bi-level programming model was developed with the objectives of arriving at the UE solution in the lower model and minimizing the LCSS transportation cost in the upper model.

•
To solve this model, a genetic algorithm embedded with the Frank-Wolfe algorithm for UE-based LPG flow assignment was proposed.This algorithm was found to have reasonably

4. 1
.1.Chromosome Coding Chromosome integer coding is used in GA, which includes eight gene segments expressing eight transporting periods of the LTN during the whole day.Besides, each gene segment has six gene positions presenting the vehicle number on six LPG transporting lines.Therefore, each gene segment owns 48 (8×6) positions and can be described by { } , k p y k K p P ∈ ∈ , say, one chromosome is a possible LPG transport operation with the line choice and the vehicle number.Every gene position is assigned an integer value that it belongs to 0, p a     that represents the vehicle number in the corresponding LPG transport operation.The upper limit p a is the maximum vehicle number that the LCSS can dispatch in period p.The chromosome structure is shown in Figure 2.

1 )( 2 )( 4 )
Initialize.Let impedance function c k(0) p (x k p ) = c k p (0), ∀p, k and distribute the LPG flow x k p from the LCSS to the users in all p periods using the All-or-Nothing assignment method, which is to obtain the flow value of each user g p,k(1) s .Then, let iteration time λ = 1.Calculate the impedance value of all periods c k(λ) p (x k p ) = c k p (x k p ), ∀p, k after iterating λ times.(3) Seek the feasible LPG flow distribution.Based on c k(λ) p (x k p ) , distribute the LPG flow volume of x k(λ) p to each line and each period by the All-or-Nothing assignment method so we get an assistant flow volume b k(λ) p .Calculate the iteration steps α (λ) by solving the following two functions: min

( 6 )
Stop the iteration if meeting the convergence principle and calculate the fitness of each individual in the upper model and turn to Step 3; otherwise, λ = λ + 1 and turn to Step 1.

5 :Step 6 :
this new chromosome.If it meets constraint(9), then turn to Step 3; otherwise, repeat the above process.Likewise, if the current chromosome does not meet constraint (4), then randomly choose a gene position from all non-zero gene positions in the corresponding chromosome segment.Set 1 new chromosome.If it meets constraint (4), then turn to Step 3; otherwise, repeat the above process.Step Termination of testing If the iteration time Ge ϕ > , then turn to Step 6; otherwise, turn to Step 2. Stopping and outputting Output the decision variable k p y with the maximum fitness in the current generation and the LPG flow distribution result of k p x in the lower programming model.

Figure 1
Figure 1 shows an example of an LTN with one LCSS, 10 users, and six different lines for assessing the LPG transportation operation bi-level optimized model and the genetic algorithm embedded with the Frank-Wolfe algorithm.

Figure 4 .
Figure 4. Length of each edge in LTN.

Table 4 .
Parameter values of length of line k and its cost of

Figure 4 .
Figure 4. Length of each edge in LTN.

Table 4 . 2 .
Parameter values of length of line k and its cost of τ k p .Analysis of Calculation Result

Figure 5 .
Figure 5.The minimum objective value in each generation of the upper model.

Figure 5 .
Figure 5.The minimum objective value in each generation of the upper model.

Figure 5 .
Figure 5.The minimum objective value in each generation of the upper model.

Figure 7 .
Figure 7.The number of vehicles departing from the LCSS (LPG cylinder supply station) of a 60,000 kg LPG demand.

Figure 8 .
Figure 8.The number of vehicles departing from the LCSS of a 72,000 kg LPG demand.

Figure 8 .
Figure 8.The number of vehicles departing from the LCSS of a 72,000 kg LPG demand.

Figure 8 .
Figure 8.The number of vehicles departing from the LCSS of a 72,000 kg LPG demand.

Figure 9 .
Figure 9.The number of vehicles departing from the LCSS of a 90,000 kg LPG demand.

Figure 9 .
Figure 9.The number of vehicles departing from the LCSS of a 90,000 kg LPG demand.5.2.4.Comparison of LPG Flow Distribution Results of Three Different LPG DemandsBased on the above analysis, the optimized LPG transporting scheme must ensure similar and minimum impedance value c k p for each period.These values of each LPG demand are shown clearly in Table7.

Figure 10 .
Figure 10.Comparison of k p x of three different LPG demands.

Figure 10 .
Figure 10.Comparison of x k p of three different LPG demands.

Table 1 .
Parameter values for the impedance function.

Table 2 .
Parameter values of user s LPG demand q s .

Table 3 .
Value of ε

Table 3 .
Value of

Table 7 .
Impedance value of c k p for different LPG demands.

Table 7 .
Impedance value of