Route Planning for Military Ground Vehicles in Road Networks under Uncertain Battlefield Environment

Route planning for military ground vehicles in the uncertain battlefield is a special kind of route planning problem, as the military vehicles face a great of uncertain and unpredicted attacks. This paper models these uncertainties in the road network by a set of discrete scenarios. A 𝑘 th shortest-path method is introduced to find intact routes from the origin to the destination for each vehicle. A binary integer programming is presented to formulate the problem. As the combination of the uncertainties results in a huge number of scenarios, we employed the sample average approximation method to obtain a robust solution for the problem. The solution approach is illustrated and tested through three road networks with different scales. The computational results show that, for networks of small scale, our method can provide a good solution with a sample of small size, while, for the large network, with sample of small size, this method usually leads to a suboptimal solution, but a good solution can still be obtained as the sample size grows bigger. In addition, variation trend of the deviation with different sample size indicates that a sample of larger size can bring more stability to the results.


Introduction
During wartime, Military Ground Vehicles (MGVs) usually have to maneuver out of their rear bases or concealed places (origins), pass through a complex road network, and arrive at the predetermined areas (destinations) to complete their operational missions, as in Figure 1. For example, mobile missile launchers depart from the military bases and move to the missile launch sites to complete the launch mission. Generally, origins are well fortified and disguised, which can ensure the safety of the MGVs. However, when the MGVs depart from the origins and move on the road network, they may be detected by the enemy's reconnaissance sensors. Because of the lack of fortification and camouflage resources, MGVs are exposed to man-made attacks during their movement in the road network. Such attacks, mainly arising from ground assault and air raid, are likely to occur frequently in battlefield. In fact, according to South Asia Terrorism Portal (SATP) [1], there were 309 attacks to the North Atlantic Treaty Organization (NATO) supply transport in Pakistan from 2008 to 2014. The number of deaths due to such attacks attains to 143. The problem of maneuvering in networks prone to attacks is prominently important due to the numerous zones that are currently at war around the world. According to the data given by the Global Security Organization (GSO) [2] there were 35 major ongoing wars and more than 25 minor conflicts around the world by the end of 2015. All these regions are highly prone to attacks where enemies perform attacks to military staff and vehicles travelling in the road networks. Thus, properly planning routes for MGVs in uncertain battlefield road networks, so as to reduce the risk of being attacked, is an important issue worth studying in the military field.
A classical problem related to MGVs' route planning is the Convoy Movement Problem (CMP) [3][4][5]. In the CMP, there is a set of convoys with a specific origin-and-destination (OD) pair. The objective is to find the routes for every convoy that minimize the total overall travel time while adhering to some strategic constraints [6]. Time is the most important factor in this problem since present-day military engagements emphasize the need for mobility more than ever [3]. However, if we consider route planning associated with a mission, MGVs need not accomplish the maneuver mission as early as possible but must reach destinations before the mission deadline. Thus, in this paper we will consider the routing time of each MGV as a strong constraint rather than the objective for costs. Given the deadline of the mission, all MGVs involved in the maneuvering are requested to arrive at their destinations in a time no longer then this timeline.
Another important feature of the studied MGV routing problem is that the destination for each MGV is not deterministic and there are a number of candidate destinations for the MGVs to choose, which is quite different from the oneto-one OD pairs in CMP. Destination location is a critical decision for MGVs, which can be coordinately optimized with the routes. This is a practical requirement in military operations; for example, there usually are a set of potential launch areas for the mobile missile launchers to complete their mission and the final launch destination has to be determined according to the situation of the battlefield and the road network. Thus, the MGV routing problem includes two main parts: for each vehicle, select a destination among the candidate locations and a route among the th-shortest paths between its origin and destination.
The emergency evacuation problem is concerned with shelter selection and route planning for evacuates, which displaces people from evacuation sources (areas suffered natural or human-made disasters) and transports them to evacuation destinations (safe shelters) [7]. The evacuation problem is related to our research in the location and routing decisions, but most of the work in this field needs not to consider the uncertain attacks during the movement in the road network. A comprehensive review of this topic can be seen in the work of Altay and Green III [8] and Galindo and Batta [9].
In the evacuation problem, there are usually a couple of candidate shelters and a given road network with deterministic delivery capacity in each road arc. The decision maker has to determine the suitable shelters and choose paths for all evacuees to optimize the evacuation cost and efficiency. Yamada [10] modeled the city road network as an undirected graph, and by solving a shortest-path problem on this graph, they obtained the optimal evacuation plan for assigning each resident in the city to one of the places for refuge. Cova and Johnson [11] focused on traffic delays occurring at intersections. To reduce these delays, they presented a network flow model for identifying optimal lane-based evacuation routing plans in a complex road network. Besides network-flow-based approach, the th-shortest-path method is another method widely employed in the evacuation problem. Campos et al. [12] first presented this method for vehicle flow allocation in emergency transportation planning. Their goal is to send a greater number of vehicles in a minimum time outside a region under menace of some catastrophic event. Stepanov and Smith [13] applied the / / / state dependent queuing models to cope with congestion and time delays on arcs. In their model, the set of feasible and potential egress routes is defined with the th-shortest-path algorithm. He et al. [7] presented a -shortest-path-based evacuation routing method with consideration of police resource allocated in transportation network. They proposed a nonlinear mixed integer programming to minimize overall evacuation clearance time. In general, the main purpose of the evacuation is to minimize the overall travelling time, and the uncertainty in the evacuation problem usually arises from the congestions and time delay occurring at intersections or arcs, which is different from the uncertain attacks considered in this paper.
To the best of our knowledge, existing studies rarely considered the arcs' uncertainty that arcs may be disrupted or travelling vehicles may be attacked during moving on the arc. In wartime, it is important for MGVs, such as missile launchers and movable radars, to quickly and safely arrive at their destination for completing the mission. In this work, we model the road uncertainty by a set of discrete scenarios. Each scenario, occurring with a specific probability, specifies a subset of arcs that become potentially insecure, on which MGVs travelling may be destroyed by intentional attacks. The optimization objective is to maximize the expectation that all MGVs accomplish the maneuver mission securely.
The main contributions of this paper are as follows.
We employ the -shortest-path method to describe the scenario-based route planning problem. A binary integer programming (BIP) model is developed for the MGVs' route planning in battlefield road network, considering arcs' insecurity and the uncertainty in travelling time.
In each scenario, the instance of the problem is deterministic and can be solved by CPLEX solver [14]. To cope with the huge number of scenarios caused by uncertainties in the road network, we adapt the sample average approximation (SAA) method to address the computational complexity and propose a solution approach based on SAA.
Three road networks with different scales are presented to test the performance of the SAA-based algorithm. Sensitivity analysis on the sample size is conducted to search for a balance between efficiency and accuracy.
The rest of this paper is organized as follows. First we present the formulation of the problem in Section 2; next we introduce the SAA background and describe the solving process based on SAA in Section 3; then the computational experiments are presented and results are discussed in Section 4; the conclusion is given in the final section.

Model Formulation
In this section, we present the problem description in detail, the underlying assumptions, and the model formulation.
We define the problem on an undirected road network graph = ( , ), where is the node set and is the set of arcs. There are a number of origins for the MGVs, and a set of potential destinations. The total number of MGVs in all origins is , which satisfies ⩽ ⩽ . The security of arcs is assumed to be independent. Considering a route , containing arcs 1 , 2 , . . . , , with secure probability (i.e., the probability a vehicle travels through the arc safely) 1 , 2 , . . . , , respectively, the secure probability of route can be formulated as (1) We are supposed to select a route for each MGV so that the sum of all routes' is maximized. We employ the th-shortest-path method to help find the possible routes for each OD pair in the road network, based on which the routing problem of MGV is transformed into MGV allocation problem. Through allocation of MGVs to different th-shortest paths of OD pairs, their route and destinations can be optimized to maximize the overall security.

Assumptions.
In order to facilitate the formulation of the problem, several underlying assumptions are proposed.
There are two statuses for each arc in the road network: secure and insecure. For an arc in the secure status, MGVs can safely travel through the arc, while, for an arc in the insecure status, there is a probability that the MGV is attacked.
Each road arc is associated with a capacity limiting the maximum number of MGVs that can travel through. We consider the capacity of each arc is limited in the whole maneuvering process. Unlike the hypothesis in CMP [3] that no pair of convoys is allowed to occupy the same part of the network at the same time, we employ the vehicle hiding description presented by Yang et al. [15] as the hypothesis. In their work, they set V ( ) as standardized value of the flow between node and in a given transport program , and the information entropy between node and is expressed as It is observed that the information entropy (V ( )) increases along with the V ( ) decreasing. If all vehicles are distributed to different roads, we can obtain higher network information entropy, which also provides a better performance in elusiveness. Therefore, the capacity of each arc is assumed to be artificially limited in the entire maneuver process in our paper.
The MGVs considered are identical and are able to execute missions independently. To ensure the safety of vehicles and reduce losses during the moving process, sending vehicles to separate terminals is a bright choice. In this paper, we suppose that each destination can accommodate one vehicle at most.
As destinations are also assumed to be homogenous in this paper, the main issue we need to focus on is to select a route to an alternative destination for each vehicle.

Notations.
To describe the problem more clearly, the notations used in the model are presented as follows.

Indices
: index of origins, and = 1, 2, . . . , , where is the total number of origins; : index of alternative destination, and = 1, 2, . . . , , where is the total number of destinations; : index of th-shortest path, and for each pair of OD = 1, 2, . . . , , where is the total number of shortest paths considered between origin and destination ; : index of arc , and = 1, 2, . . . , , where is the total number of arcs involved in the road network; : index of scenarios, and = 1, 2, . . . , , where is the total number of scenarios considered.

Parameters
: the number of MGVs dispatched from origin ; : the probability that scenario occurs; : the probability that arc is in insecure state; : the probability that a vehicle travels through arc safely in scenario , where the arc can be either secure or insecure; : the probability a vehicle travels through the insecure arc safely; : 1 if arc is in secure state in scenario , 0 otherwise; : the secure probability of the th-shortest path from origin to destination in scenario ; : a set of arcs, composed of all arcs in the th-shortest path from origin to destination ; : time consumed for travelling through arc in scenario , which is obtained considering distance of the arc and the vehicle velocity in the specific scenario; : time length of the mission, that is, the time allowed to finish the movement; : 1 if arc is included in the th-shortest path from origin to destination , 0 for otherwise; : the capacity of arc .

Decision Variable
: 0-1 variable. 1 for dispatching a vehicle from origin to destination through the th-shortest path, 0 for otherwise.
As each scenario specifies a subset of arcs which become insecure, given scenario (in which , ∀ are known), 4 Journal of Advanced Transportation the probability that it occurs can be obtained through the following function: (1 − ) . (3)

Formulation.
We formulated the problem as follows: The objective function of the problem is to maximum the sum of expectations of all vehicles travelling securely to their destinations. Constraints (5) ensure that all vehicles at origin are dispatched. Constraints (6) restrict that the number of vehicles in each destination should be no more than one. Constraints (7) present that in any scenario, each vehicle's moving time should be less than the time length restricted by the mission. In our paper, the time consumed in each arc, , is assumed to different under different scenario. Constraints (8) show that if an arc is in secure status, the MGV can pass through it safely with probability 1, and the safely travelling probability turns to be if the arc is insecure. Constraints (9) are the variants of formula (1) in an appropriate form. Constraints (10) show that the flow in each arc is limited. Constraints (11) state that decision variables, , are all 0-1 variables.
It is observed that in a specified scenario, the problem is a deterministic BIP, which can be solved by CPLEX solver. If there are arcs involved in a problem and each arc can either be in secure status or insecure status, 2 scenarios need to be considered. Thus, the number of scenarios grows exponentially with the number of arcs. It would be a tough work to get the optimal solution with all scenarios considered. To cope with this problem, we introduce the SAA method to solve the problem in the next section.

Solving Methodology
In this section, we present the solution approach based on the SAA method.

Problem Description Based on SAA.
The SAA method is a Monte Carlo simulation based approach to stochastic discrete optimization problems [16]. In this method, the expected value of the objective function is approximated by solving the problem for a sample of scenarios. The sample size of the SAA problem is much smaller than the number of scenarios in the true problem and as the sample size increases, SAA converges exponentially fast to the true problem [16]. The SAA method has been applied in many stochastic programs lacking exact algorithms, for example, the network design [17][18][19] and empty container repositioning [20,21].
The SAA method repeats the sampling and solving process for several times to select a best solution. In each sample, scenarios are randomly generated, and thus the objective function (4) of the true problem is approximated witĥ( For each scenario in a sample, both the time, , and the secure probability, , of each arc are fixed. The true problem (4)- (11) in each scenario can be rewritten as ST (5) = ∏ ∈ , ∀ , , .
This problem is a deterministic BIP, which can be solved by CPLEX solver. However, since the time, , varies from scenario to scenario, a solution suiting a scenario well may be an infeasible solution to another scenario. The solution we selected should be universally applicable. A solving process to this matter is given in Procedure 1.
First, we solve the problem of each scenario and get the optimal solutions separately. The solutions recurring for many times are thought to be of better performances than these hardly recurring, which are optimal in just some extreme scenarios. With all scenarios considered, we keep several solutions with the most recurring times. (In fact, in the experiment it is observed that the first three or four solutions consist of best solutions to more than a half of all scenarios.) From these solutions, the solution which is feasible in all scenarios and has the largest value of formula (12) is selected as the solution to the sample, denoted by .
The process of sampling scenarios and solving to obtain and is repeated times. Calculate the average of the It is noteworthy that the sample size is a tradeoff between computational complexity and accuracy. A larger value of brings a more accurate estimate of the true problem, but an increase in computing complexity as well. The size can be selected dynamically. The SAA solving process starts with a smaller size, and then we can regulate the sample size based on the acquired accuracy and the computational complexity. The selected solution, given a candidate solution , for example, a feasible solution from the optimal solutions of samples, needs to be evaluated. The evaluation of objective function is conducted by taking sample size of ( ≫ ). That is, a number of scenarios are randomly generated for the evaluation. Let Since is much larger than , the objective value obtained ( ) is considered as an estimator of the true objective function (3) [19]. In a common SAA problem, the quality of the solution is evaluated by calculating the optimality gap and the variance of optimality gap. In our problem, as the involved expectations are all numbers around 1, to highlight the dispersion of the result, we replace the variance of optimality gap with the deviation of optimality gap. Here, we define the optimality gap as gap fl −̂( ) .
And we define the deviation of optimality gap aŝ +̂fl̂+̂̂, wherêfl 3.2. SAA Algorithm. Based on description above, the SAA algorithm process is presented in Algorithm 1.

Experiment Design.
In this part, we examine the performance of the algorithm on different dataset. We selected three road networks: a small-scale designed network with 13 nodes    The th-shortest-path datasets were derived by the methods in [22][23][24]. Other key parameters are given as follows: is a constant and is set to be 0.2 in all experiments; the value of is set to be the road length divided by the velocity, which varies in a range around ±10% of a given velocity as the scenario varies; is randomly generated.
The data for all the experiments can be found in the supplementary material noted as dataset. Table 1 shows a brief look of the problems considered in 3 networks, respectively. For the fulfillment of SAA, the value of is set to be 10 and keeps constant for all sample sizes. The sample size ranges from 100 through 1000, with an interval of 100. And the value of is set to be 10,000. The algorithms were coded in MATLAB by invoking CPLEX v.12.6. And all the programs were executed on a 64-bit computer with a 3.2 GHz CPU and 8.0 GB of physical RAM.

Computational Results and Analysis.
The computational results for the three networks are presented in Tables 2,  3, and 4, respectively. We give the minimum, maximum, average, and corresponding deviation of each sample size, as the results of the SAA process. The calculating time of each SAA process is also listed as the CPU time. It can be seen that the calculating time increases as the sample size grows, where the growth rate keeps almost constant for each network separately. This is an acceptable result, since that we viewed the working time of each module of the program and found the scenario generating and CPLEX solving parts, both of   which require repeating (the value of sample size) times and take the most of the time (more than 95% of the whole). The evaluation part shows the results of employing the selected solution of the SAA process to a much bigger sample of size . And the gaps between these results are given in the optimality gap part. Several discussions about the optimality gap and the deviations are given as in Tables 2, 3, and 4. Figure 5 shows the trends of gap (%) for the 3 networks while the sample size increases. Here the index, gap (%), means the closeness between the SAA results and evaluation results, which can be viewed as the estimator of the true problem. For Network 1, the gaps (%) are all below 1% and fluctuate irregularly, while, for Network 2 and Network 3, the values start from 4.46% and 5.69%, respectively, as the sample size is 100 and decrease distinctly as the sample size increases. Although there are several fluctuations in both of them, the downtrends are obvious. These differences between Network 1 and Network 2/3 arise from the solutions selected in each network. In fact, the solution selected for Network 1 stays constant from sample size 100 through 1000.
The fluctuation is caused by the uncertainty of sampling. However, the solutions selected for Network 2 and Network 3 vary as the sample size increases. Better solution brings smaller optimality gap, which indicates that better solution to the true problem is obtained. Figure 6 presents comparison for the trends of SAA results and evaluation results. We can see the evaluation results of Network 1 stay relative constant, while the SAA results fluctuate around the evaluation results in a quite small range. For Network 2, the SAA results have a distinct downtrend as the sample size increases, while the trend of the evaluation results is opposite. The gap between the two values is first below 1% when sample size is 500, after which the SAA results fluctuate around the evaluation results. The smallest gap between the two values reaches 0.013% when sample size is 900. The discussion for Network 2 can be applied to Network 3 as well, in which these trends are more obvious. The SAA results rise and evaluation results fall distinctly as the sample size increases from 100 through 600. From 600 through 1000, the SAA results fluctuate. And the   smallest gap between the two is 0.08%, when the sample size is 700.
The results of Network 1 indicate that, for a network of small scale, a good solution can be obtained with a considerably small sample size. But for networks of bigger scales, such as Network 2 and Network 3, as the number of arcs considered in the problem increases, the number of possible scenarios increases. A sample of small size cannot reflect the panorama of the problem. In our algorithm, 10 samples are generated and solved, and the sample solution with the biggest expectation value is selected, while, for a sample of small size, a suboptimal solution may have a better performance than the "good" solution (as we cannot tell whether it is the optimal one) does. This explains why for networks of bigger scales, suboptimal solutions are usually selected when the sample size is small. But if the sample size is big enough, the good solutions can still be obtained constantly. Figure 7 presents the deviation trends with the sample size increasing. All the deviations of the three networks have downward trends, which are especially evident in Network 2  and Network 3. The deviation consists of the sampling deviation and the evaluation deviation. The evaluation deviation reflects the quality of a given solution. As all solutions selected in sample sizes are identical in Network 1, the evaluation deviations vary a little. The variation of the deviation is mainly caused by the sampling deviations, which manifests volatility. For networks 2 and 3, of each sample size, the evaluation deviation is much bigger than the sampling deviation. And as the selected solutions vary, the evaluation deviations have evident downtrends. In fact, the sampling deviations of Network 2 and 3 are fluctuating, of which the overall trends are drops as well. As the longest time spent on sampling and calculating is just 288.78 s (for Network 3, of sample size 1000), the time spent on calculating is relatively short and grows linearly with the growth of the sample size, we can solve the problem with bigger sample size to get a good solution of higher accuracy to reflect the true problem better.

Conclusion
In this paper, we addressed the problem MGVs maneuver in unsecure road networks in battlefield environment. We describe the uncertainty of the network based on discrete scenarios and introduce the th-shortest-path method to calculate the possible paths between origins and destinations. The problem is formulated into a binary integer programming model. Due to the huge number of scenarios generated by the combination of uncertain arcs, we proposed a SAA method to solve the model. Experiments based on three networks with different scales were conducted to test the efficiency of the solution approach. Computational results show that bigger sample brings relevant progress of the approximation and an acceptable increasing in computational complexity.
In our paper, we did not consider the fortification or camouflage resource for the maneuver, and we only need to decide routes for all vehicles. In some cases, there are protective resources available in the MGV moving process, which can be allocated to some road arcs and reduce the risk for being attacked. Therefore, it is a valuable extension for integrating the protective resources in the MGV routing problem in future works.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.