An orienteering model for the search and rescue problem

  • Published on
    14-Dec-2016

  • View
    216

  • Download
    3

Transcript

Comput Manag SciDOI 10.1007/s10287-013-0179-1ORIGINAL PAPERAn orienteering model for the search and rescueproblemAdel Guitouni Hatem MasriReceived: 6 February 2013 / Accepted: 3 July 2013 Springer-Verlag Berlin Heidelberg 2013Abstract In this paper, we propose a new model for the search and rescue problem.We focus on the case of a single airborne search asset through a connected spaceand continuous time with a maximum travel time T . The intent is to maximize thedetection of a cooperative target (search and rescue). The proposed model is basedon the assumption of existing a priori information (e.g., result of information fusionprocess) to establish a spatial distribution of probability of containment in possiblegeographic locations. The possibility area is defined using a cut threshold on theprobability of containment and the search path as well as the allocation of the level ofeffort to each region in the search space is obtained based on an orienteering model.We illustrate the application of the proposed model on an empirical example.Keywords Search and rescue problem Orienteering problem1 IntroductionIn this paper, we consider the search path problem for a single airborne search assetin a continuous space and time for a maximum mission time T , through a connectedregion. The objective is to maximize the detection of a cooperative target (search andrescue). A priori information (e.g., sensors data, experts opinion) is available to buildA. GuitouniCommand and Control Decision Support Systems Section, Defence R&D Canada,2459 Pie XI North, Quebec, QC G3J 1X5, Canadae-mail: adel.guitouni@drdc-rddc.gc.caH. Masri (B)College of Business Administration, University of Bahrain,P. O. Box 32038, Sakhir, Kingdom of Bahraine-mail: hmasri@uob.edu.bh123A. Guitouni, H. Masria spatial distribution of the possible location of the target. Such distribution might bea probability, possibility, fuzzy or any other type of uncertainty modelling functions.A seeker or a searcher is an airborne platform equipped with sensors (e.g., scanningradar, imaging radar, electro-optic, infrared). The airborne platform has a maximumfly time that might be seen as time to refuel or any other system or human constraint.The problem is therefore to sequence a set of search activities in time and space inorder to maximize the chances to find the target. Search activities are in this case thetransit between two regions and active search in each region. The search plan shouldthen provide a path plan as well as the level of effort to be spent in each location.The search-path problem has been closely related to search effort allocation prob-lem (Kierstead and Balzo 2003). Search theory supposes that targets and seekers aremodelled probabilistically (Stone 1975). Search theory was developed for the first timeby Koopman (1980) and the US Navys operations research group in order to developbetter strategies for anti-submarine warfare during the Second World War. Since then,search theory was mainly applied in application fields such as: surveillance, SAR mis-sions and exploration (e.g. finding the US Nuclear Submarine Scorpion lost in 1968and Clearing unexploded ordnance in the Suez Canal).The search path planning problem has been subject to many ongoing researchactivities. In this paper, we propose a novel modelling for the search problem froma mathematical programming perspective. This paper is organized as follows. Wepresent an overview of the search and rescue problem and we review related work inSect. 2. In Sect. 3, we describe a variant of the search and rescue problem and theway that we are going to proceed to solve that problem. First, we start by buildingthe possibility area in Sect. 4 and then adapting in Sect. 5 an orienteering model tobuild the search path and to allocate the search effort to each region in the possibilityarea. We present an empirical validation of the model in Sect. 6. The conclusion ofthis paper is given in Sect. 7.2 Search and rescue problemResource management for conducting search and rescue, surveillance and reconnais-sance mission is a very important problem for many organizations. This problem ischaracterized by the employment of mobile (e.g. maritime patrol aircraft, helicopters,UAVs, ships) and fixed surveillance assets (e.g. land radar) to a large geographic areain order to identify, assess and track the maximum number of moving, stopped or drift-ing objects. The observed objects are not necessarily aware of being observed and arecooperative or non-cooperative, and friendly or hostile. The scarce surveillance (e.g.electro-optical, infra-red, and synthetic aperture radar sensors), tracking capabilities(normal radar modes), scarce platform, time and space constraints makes it a verydifficult problem.A search mission is usually planned for a given search area while minimizing costin term of time, distance and exposure to threats, and maximizing the chances ofsuccess subject to a set of constraints imposed by the amount of effort available forthe search. The search path planning problem is intractable and complex because it is123Orienteering model for the search and rescue problemmultifaceted and NP-hard even for stationary targets (Trummel and Weisinger 1986).Many elements should be considered including: the platform configurations and properties which generally related to the physicaland functional properties of the platform (aircraft), the number of platforms and their configurations, the sensors, or the instruments used to carry out the seeker such as the eye, radar,sonar, television, and cameras, the object or the target being searched which may be stationary or moving, coop-erative or non cooperative, friendly or hostile, the geographic or physical search area that may be continuous (in Euclideann-space) or discrete (a set of cells), closed or open, and the available search effort may be continuous (i.e. measured by time or track length)or discrete (i.e. measured by a finite number of scans or looks).Abi-Zeid and Frost (2005) developed a decision support system (SARPlan) based onsearch theory and optimization to enhance the SAR mission effectiveness. SARPlanaims at maximizing the efficient use of resources and increasing the chances of findingsurvivors in a shorter time. Based on search theory, SARPlan suggests an optimal over-land air plan for a SAR mission. In SARPlan, Constraint Satisfaction Programming(CSP) and traditional optimization techniques are used to derive the optimal plan.Note that an optimal plan is a plan that maximizes the probability of finding the searchobject and the efficient use of resources while minimizing search time, search area,and associated costs. Dor et al. (2009) proposed an extension of the search theoryfor search and rescue. They proposed the theory of belief functions for informationcombination and update for the optimal search planning context.The discrete path problem is based on the assumption that the searcher or seekerand the target are roaming in a discrete time and space dimensions. Moreover, discrete-time moving-target search model usually assumes the Markovian property of thetargets movement (Hong et al. 2009). Dambreville and Cadre (2002) proposed analgorithm for optimal search for a target following a Markovian movement or a con-ditionally deterministic motion when search resources are renewed with generalizedlinear constraints.An optimal dynamic programming method is proposed by Eagle (1984). Stewart(1979) has proposed a search efforts relocation heuristic where the effort is restrictedto paths. Eagle and Yee (1990) proposed an extension of the work of Stewart (1979).They proposed an optimal branch-and-bound algorithm limited 7 7 search grids.Other extensions of this work might be found in Dell et al. (1996) and Hohzaki andIida (1997).Hong et al. (2009) proposed a pseudo-polynomial heuristic for the single-searcherpath-constrained discrete-time Markovian-target search. The heuristic is based on anapproximate of non-detection probability computed from the conditional probabil-ity that reflects the search history over a fixed time windows. Jacobson and McLay(2006) proposed simultaneous generalized hill climbing algorithms to determine opti-mal search strategies over multiple search platforms. Kierstead and Balzo (2003)proposed a genetic algorithm (GA) for planning search paths against a moving target.In his paper, Janez (2007) proposed to model the sensors planning problem as vehicle123A. Guitouni, H. MasriDestress CallOrIntelligencePrioritization ProcessFusion ProcessSensor 2Sensor 3Sensor 1Spatial Density (Possibility) DistributionFiltering ProcessResource Allocation Path PlanningTasking OrdersPlatform Flying the missionAction/EffectsAndAssessment Fig. 1 Simplified processrouting problem (VRP). The formulation is based on the assumption that the coordi-nates of the targets to be visited is known a priori. Unlike previously mentioned works,targets are stationary and their location is well defined. The problem addressed is thenthe allocation of different surveillance assets to visiting each target.3 The problem descriptionLet a seeker s be a combination of a given platform with appropriate sensors for agiven search situation. Multiple seekers might be assigned to different search areas.Search path in each area might follow a static or dynamic pattern. This paper addressesoptimal path planning given a particular seeker. We focus on airborne seekers.We assume the process shown in Fig. 1. A distress call might be received by asearch and rescue centre, which will activate its prioritization process. Informationis received from different information sources (e.g., sensors). An information fusionprocess produces a spatial distribution of possible location of the target. Through afiltering process, one can generate a synthesized distribution function over the searcharea.123Orienteering model for the search and rescue problemWe assume that a priori information is available through information fusionprocesses about a plausible distribution that represent the possibility that the target iswithin a particular containment area (Kao et al. 2001). The intensity of the distributionmight refer to a higher degree of confidence (e.g., probability) that the object mightbe in a particular more confined region of that containment area. At this point, thedistribution function is not necessarily a probability distribution function. Differentclustering techniques might be applied in order to create a more intelligible distributionmap.A good search plan should maximize the probability of success (POS), which is theprobability of finding the search object. To compute the POS, two other probabilitiesare used. The probability of containment (POC), which is the probability that the searchobject is in a particular cell of the searched area, and the probability of detection (POD),which is the probability that the sensor detects the search object given that it is in thearea searched. Note that the probability POD is function of the effort provided by thesensor on a specific area. For simplicity of the model, we assume that the probabilityof detecting the target is function of the time spent searching in a particular geographicregion. We simplify the conditional probability to be expressed in term of the lateralrange function of the sensor, its sweep width W , and the total time spent in the region(Frost 1998).We assume that the probability of detection is given by the following expression:POD = 1 e(W j v jA j)t j (1)where W j is the effective sweep width in region j, A j is the surface of region j, v j isthe platform speed in region j and t j is the level of effort allocated to region j . In thispaper, we propose a path planning algorithm to maximize the probability of findingthe cooperative target within a given time window. The effort may be measured bytime, area searched, track length or any other appropriate measure. The relationshipbetween these three probabilities is as follows:POS = POC POD (2)Now, the problem is to define the possibility area and to fly an airborne seeker withthe intent to find the target while optimizing a set of objectives under constraints.4 The possibility areaThe first step of this process consists in identifying the boundaries of the possibilityarea (PA) that may contain the search object. Subjective (e.g., rules, experts judg-ments) and objective (e.g., historical data, fused sensors data) knowledge are generallyused to define the limits of this area. For instance, the International Aeronautical andMaritime Search and Rescue (IAMSAR, 1999) Manual contains some guidelines toestablish the possibility area. Other approaches like possibility distribution functionor fuzzy functions might be used. Then, the possibility area is divided into J cells123A. Guitouni, H. MasriFig. 2 2D view of a distribution of possible locations of the targetFig. 3 Cut( j = 1, . . ., J ) that represent the smallest regions over which search effort can beallocated. A search within each cell might follow a predefined or dynamic pattern.Lets assume the search problem shown in Fig. 2. An airborne platform shouldfly over a geographic area following a given path to maximize the chances to findthe target. The aircraft has a limited flying time and the target might have a limitedsurvival time. The different contours represent distribution function of the probabilityof containment. The concentration of the contours represents higher probabilities.We propose to apply a cut threshold and then we obtain the distribution of Fig. 3.The value of threshold may have an impact on the search area. A higher value of thethreshold can allow the airborne platform to fly over different regions and spend moretime where the contours are denser and less elsewhere. Such flying path might berepresented by Fig. 4.123Orienteering model for the search and rescue problemFig. 4 Flight plan12i j(tij)i j(tij)345768(ti)(t3)(t4)(t5)(t6)(t7)(t8)(t2)Fig. 5 Example of flying pathTherefore, the search regions are the vertices of a graph G = (V, E) (see Fig. 5)where the search platform starts at the vertex i = 1 to visit the other regions (vertices).To travel between two vertices i and j , the search platform takes ti j time (tii = 0 forall i). Note that, we consider only vertices that verify {i | t1i + ti1 TM } where TMis the time window for the search and rescue mission.5 A orienteering modelThe search and rescue problem using only one search platform (SARP1) presentssome similarity compared to the orienteering problem with variable profits (OPVP).The OPVP is to find a set of vehicle tours that begins and ends at the start vertex,traverses a subset of the set of vertices and must have a total duration not exceeding123A. Guitouni, H. Masria given limit (maximum mission time) where at every vertex the vehicle collects apercentage of the profit that depends on the time spent on that vertex (Erdogan andLaporte 2013). The difference between the OPVP and the following proposed modelfor the SARP1 is that here we consider the MillerTuckerZemlin (MTZ) formulationto deal with the subtour elimination problem.5.1 Variables xi j is a binary variable. If xi j = 1 then the search platform, after visiting the regioni , will over fly the region j . Otherwise, such a travel between region i and j willnot be considered in the optimal flight plan. ti is the time of visit that the search platform will take to over fly the region i (levelof effort allocated to the region i)5.2 The constraintsThe first constraints restricts the time of search and rescue to the given time of missionTMni=1nj=1j =iti j xi j +ni=2ti TM (3)t j TM ni=1i = jj = 2, ..., n (4)The second set of constraints ensures the flow conservation in vertices:ni=1i = jxi j ns=1s = jx js = 0, j = 2, ..., n (5)The third constraints concerns that any tour must include the vertex 1ni=2xi1 = 1 (6)The last set of constraints has to ensure the non-existence of subtours. In the literature,many formulations were proposed for the sub tours elimination problem (Miller et al.1960). Among them the MTZ formulation (Kara et al. 2004):u j ui 1 + p(xi j 1), i = j, i = 2, . . . , n, j = 2, . . . , n1 u j n 1, j = 2, . . . , n (7)123Orienteering model for the search and rescue problemwhere ui , i = 2, . . . , n, are positive variables and p is the number of visited verticesin the tour. For the case p = n, the MTZ formulation for the travel salesman problemmay be used. Many extensions and lifting of the MTZ formulation were proposed inthe literature essentially for the travel salesman problem (Sherali and Driscoll 2002).One of these extensions of the MTZ formulation was provided by Desrochers andLaporte (1991) to the vehicle routing problem:u j ui ci j + M(xi j 1), i = j, i = 2, . . . , n, j = 2, . . . , ns j + (c1 j s j )x1 j u j L t j (c j1 r j )x j1, j = 2, . . . , n (8)where ci j is the distance associated to the edge (i, j), s j the length of a shortest pathfrom vertex 1 to vertex j, r j is the length of a shortest path from vertex j to vertex 1,L is the maximum tour length and M maxi, j{L s j ri + ci j }.Under the following correspondences: L = TM , ci j = ti j + t j , s j = c1 j andr j = c j1, the Desrochers and Laporte formulation (8) became as followsu j ui t j + ti j + M(xi j 1), i = j, i = 2, . . . , n, j = 2, . . . , nt1 j + t j u j TM t j1, j = 2, . . . , n (9)where M maxi, j,i = j{TM t1i + ti j + t j t j1}. The set of constraints (9) is notvalid for the SARP1 as we suppose that u j t1 j + t j for all j 2, which will obligethe optimal tours to include all vertices. This lower bound is only valid if the vertex j isvisited by optimal tours (i.e. i=1ni = j xi j = 1). Therefore, we propose to reformulatethe Desrochers and Laporte formulation as follows:u j ui t j + ti j + M(xi j 1), i = j, i = 2, . . . , n, j = 2, . . . , n(t1 j + t j )ni=1i = jxi j u j , j = 2, . . . , nu j TM t j1, j = 2, . . . , n(10)where M maxi, j,i = j{TM ti1 + ti j + t j}.Proposition 1 The formulation (10) is valid for the SARP1.Proof To prove that the formulation (10) is valid for SARP1, we have to verify thatconstraints in (10) eliminate subtours and that the length of feasible subtour does notexceed mission time TM .Consider a subtour that does not include the vertex 1: l1, . . . , l p where lk = 1 fork = 1, . . . , p and l p = l1. Therefore, xlklk+1 = 1, k = 1, ..., p1 and summing from1 to p1 the constraints ulk+1 ulk tlk+1 +tlk lk+1 yields top1k=1 (tlk+1 + tlklk+1) 0which is a contradiction and then all tours must include the vertex 1.123A. Guitouni, H. MasriWe are going to prove that any subtour length does not exceed the time of missionTM . Consider the following subtour: 1, l1, . . . , l p,1 where lk = 1 for k = 1, . . . , pand Ls the length of this subtour. Summing up the constraints ulk+1 ulk tlk+1 +tlklk+1, k = 1, . . . , p 1 givesp1k=1(tlk+1 + tlklk+1) ul p ul1 (11)Since t1lk tlk ulk TM tlk 1 for k = 1, . . . , p thenLs = t1l1 + tl1 + tl p1 +p1k=1(tlk+1 + tlklk+1) t1l1 + tl1 + tl p1 + ul p ul1 TM (12)unionsq5.3 The objective functionThe objective function is to maximize the sum of the probability of success POSi ofthe visited (vertices) regionni=2POCi(1 e(Wi viAi)ti)(13)where POCi is the probability of containment for each vertex i , which represents theprobability that the object we are looking for is in the region i 2(POC1 = 0).5.4 The nonlinear mixed integer programming modelThe previous set of equations provides the following nonlinear mixed integer program:Maxni=2POCi(1 e(Wi viAi)ti)s.t.ni=1nj=1j =iti j xi j +ni=2ti TMt j TM ni=1i = jxi j , j = 2, ..., n (14)123Orienteering model for the search and rescue problemni=1i = jxi j ns=1s = jx js = 0, j = 1, ..., nni=2xi1 = 1u j ui t j + ti j + M(xi j 1), i = j, i = 2, ..., n, j = 2, ..., nu j (t1 j + t j ) ni=1i = jxi j , j = 2, ..., nu j (TM t j1) ni=1i = jxi j , j = 2, ..., nxi j {0, 1} , i = 1, ..., n, j = 1, ..., n, i = ju j , j = 2, ..., nt j , j = 1, ..., nwhere M maxi, j,i = j{TM ti1 + ti j + t j }.6 Empirical resultsIn this paper, we consider the problem shown in Fig. 6. A fixed wings surveillanceaircraft is based in location 0.The information fusion process has established the possibility area given in theFig. 6. After we apply a cut threshold, we obtain 10 regions to be the potential locationscontaining the target with the highest probability of containment.In this application, the contours of the probability of containment in each regionare circular. The aircraft transit speed between two regions is constant at 650Km/hspeed. Table 1 shows computed travel time between two pair of regions centres. Thetravel time is computed using the following formula:ti j = vR cos1[cos(longi long j ) cos(lati ) cos(lat j ) + sin(lati ) sin(lat j )](15)where R is the earth radius, longi and lati are respectively the radius longitude andlatitude of the centre of region i . This formula is based on distance calculation betweentwo points in a spherical earth, which has error typically up to 0.3 %.The probability of containment (POC) distribution for each region is obtained froman information fusion process. Table 2 provides the cumulative probability of contain-ment for each region. For simplicity, the sweeping width W j of the airborne sensor issupposed to be constant and the same for all regions at 300 m (0.3 km). Nonetheless,123A. Guitouni, H. MasriFig. 6 Surveillance problemTable 1 Travel time between each pair regions centresti j 0 1 2 3 4 5 6 7 8 9 100 0 0.296 0.6 0.352 0.399 0.809 0.997 1.074 1.453 1.158 1.3251 0.296 0 0.377 0.313 0.443 0.776 0.844 1.133 1.428 1.112 1.1972 0.6 0.377 0 0.324 0.438 0.545 0.483 0.965 1.149 0.827 0.8423 0.352 0.313 0.324 0 0.141 0.481 0.648 0.82 1.136 0.827 0.9744 0.399 0.443 0.438 0.141 0 0.413 0.664 0.701 1.055 0.761 0.9585 0.809 0.776 0.545 0.481 0.413 0 0.39 11.3 0.655 0.349 0.586 0.997 0.844 0.483 0.648 0.664 0.39 0 0.756 0.742 0.448 0.3587 1.074 1.133 0.965 0.82 0.701 11.3 0.756 0 0.491 0.491 0.48 1.453 1.428 1.149 1.136 1.055 0.655 0.742 0.491 0 0.322 0.529 1.158 1.112 0.827 0.827 0.761 0.349 0.448 0.491 0.322 0 0.36810 1.325 1.197 0.842 0.974 0.958 0.58 0.358 0.4 0.52 0.368 0the proposed model handles well a variable sweeping width. The aircraft search speedv j is constant at 350 km/h. Here again, it is possible to consider variable search speedsfor each region. Table 2 provides the value of KA(KA = W j v jA j ) and POC valuesnecessary to compute the objective function (see Eq. 13). The overall mission time is20 hours (the aircraft should land at the base airport no later than t0 + 20h).The proposed surveillance problem is then formulated according to the model givenin Eq. (14). We obtain a non-linear mixed integer program with: Decision variables: 144 including 110 integers Number of constraints: 156 with 113 nonlinear123Orienteering model for the search and rescue problemTable 2 Parameters of thesurveillance problem for eachregionRegion POC Wj(Km) vj(Km/h) KA e(KA)0 0.3 3501 0.031 0.3 300 0.974 0.3772 0.091 0.3 350 1.243 0.2883 0.084 0.3 350 1.888 0.1514 0.135 0.3 350 0.804 0.4475 0.142 0.3 350 1.079 0.3406 0.011 0.3 350 9.512 0.0007 0.152 0.3 350 0.731 0.4818 0.196 0.3 350 0.608 0.5449 0.017 0.3 350 1.787 0.16710 0.140 0.3 350 0.618 0.539Table 3 Optimal solutionRegion 0 1 2 3 4 5 6 7 8 9 10Search time 0.000 1.576 1.465 2.320 2.479 2.636 3.429 2.164Arrival time 0.000 0.600 18.182 15.721 2.721 12.384 8.464 5.780Departure time 0.000 2.176 19.647 18.041 5.200 15.020 11.893 7.944Fig. 7 Mapping of the optimal search pathWe build up a computer program using the Lingo 12.0 solver. This experiment is testedwith a Core2 duo 2 GHZ laptop. The solution is obtained after 41 min 29 s runtime.The estimated POS (value of the objective function) is 0.606167. The solution is givenin Table 3.The optimal path is 0-2-5-10-8-9-4-3-0 as shown in Fig. 7 (time is given in hours).123A. Guitouni, H. MasriThe search starts at t0 = 0. The aircraft transit from the airport and start searchingin region 2. The overall effort allocated to region 2 is 1.576 h. The search pattern inregion 1 is another optimization problem solved by the aircrew. If the target has notbeen located in region 2, then the aircraft will move to region 5 and spend 2.479 hsearch time. If the target is found, the aircraft might resume its routine operations orreturn to base (airport 0). The transit between regions is not always passive. The aircraftcontinues deploying its sensors and collecting data (the probability of containment isusually higher than 0 along the transit path).The model proposed optimizes both the path as well as the allocation of the level ofeffort to each region in the search space. The runtime to solve this problem might beseen as a limitation. The solver took most of the time to improve the optimal solutionfrom 98 to 99 %. If we accept a 98 % solution, the solver generates a potential optimalsolution in Orienteering model for the search and rescue problemFrost J (1998) The theory of search: a simplified explanation. Report by Soza & Company Ltd and Officeof Search and Rescue US Coast GuardHohzaki R, Iida K (1997) Optimal strategy of route and look for the path constrained search problem withreward criterion. Eur J Oper Res 100(1):236249Hong SP, Cho SJ, Park MJ (2009) A pseudo-polynomial heuristic for path-constrained discrete-timeMarkovian-target search. Eur J Oper Res 193(2):351364IAMSAR (1999) International Aeronautical and Maritime Search And Rescue Manual. ICAO/IMOpublicationsJacobson SH, McLay LA (2006) Optimal search strategies using simultaneous generalized hill climbingalgorithms. Math Comput Model 43(910):10611073Janez F (2007) Optimization method for sensor planning. Aerosp Sci Technol 11(4):310316Kao D, Dungan J, Pang A (2001) Visualizing 2D probability distributions from EOS satellite image-deriveddata sets: a case study. In: Proceedings of the conference on visualization 01 San Diego. California SanDiego, USA, pp 457460Kara I, Laporte G, Tolga B (2004) A note on the lifted MillerTuckerZemlin subtour elimination constraintsfor the capacitated vehicle routing problem. Eur J Oper Res 158(3):793795Kierstead DP, Del Balzo DR (2003) A genetic algorithm applied to planning search paths in complicatedenvironments. Mil Oper Res 8(2):4559Koopman BO (1980) Search and screening: general principles with historical applications. Pergamon Press,NYMiller C, Tucker A, Zemlin RA (1960) Integer programming formulation of traveling salesman problems.J ACM 7(4):326329Sherali HD, Driscoll PJ (2002) On tightening the relaxations of MillerTuckerZemlin formulations forasymmetric traveling salesman problem. Oper Res 50(4):656669Stewart TJ (1979) Search for a moving target when searcher motion is restricted. Comput Oper Res 6(3):129140Stone LD (1975) Theory of optimal search. Academic Press, LondonTrummel JR, Weisinger JR (1986) The complexity of the optimal searcher path problem. Oper Res34(2):324327123An orienteering model for the search and rescue problemAbstract1 Introduction2 Search and rescue problem3 The problem description4 The possibility area5 A orienteering model5.1 Variables5.2 The constraints5.3 The objective function5.4 The nonlinear mixed integer programming model6 Empirical results7 ConclusionReferences