A DC programming approach for planning a multisensor multizone search for a target

Computers & Operations Research 41 (2014) 231–239

Contents lists available at ScienceDirect

Computers & Operations Research journal homepage: www.elsevier.com/locate/caor

A DC programming approach for planning a multisensor multizone search for a target Hoai An Le Thi a,n, Duc Manh Nguyen b, Tao Pham Dinh b a b

Laboratory of Theoretical and Applied Computer Science, UFR MIM, University of Lorraine, Ile du Saulcy, 57045 Metz, France Laboratory of Mathematics, National Institute for Applied Sciences - Rouen, 76801 Saint-Etienne-du-Rouvray, France

a r t i c l e i n f o

abstract

Available online 20 July 2012

In this paper, we consider a well-known problem in the general area of search theory: planning a multisensor in multizone search so as to maximize the probability of detection of a target under a given resource effort to be shared. We propose a new optimization model that is a nonlinear mixed 0–1 programming problem. This problem is then reformulated as a DC (Difference of Convex) functions program via an exact penalty technique. DC programming and DCA (DC algorithm) have been investigated for solving the resulting DC program. Numerical experiments demonstrate the efficiency and the superiority of the proposed algorithm in comparison with the existing method. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Search theory Hierarchical optimization Combinatorial optimization DC programming and DCA Nonlinear mixed 0–1 programming Exact penalty

1. Introduction

 the local detection probability that a particular amount of local

Search theory is defined by Cadre and Souris [3] as a discipline that treats the problem of how a missing object can be searched optimally, when the amount of searching time is limited and only probabilities of the possible position of the missing object are given. The theory of how to search for missing objects has been a subject of serious scientific research for more than 50 years. It is a branch of the broader applied science known as operations research [6]. In fact, Search theory was first established during World War II by the work of Koopman and his colleagues [12] in the Antisubmarine Warfare Operations Research Group (ASWORG). The applications of search theory were first made on military operations [22]. Koopman [11] stated that the principles of search theory could be applied effectively to any situation where the objective is to find a person or object contained in some restricted geographic area. After military applications, it was also applied to different problems such as surveillance, explorations, medicine, industry and search and rescue operations [8]. The aim of searching in the context of Aeronautical Search and Rescue (ASAR), for instance, is to find the missing aircraft effectively and as quickly as possible with the available resources [21]. A search theory problem, in general, is characterized by three pieces of data [3]:

 the total amount of searching effort available.

search effort could detect the target;

 the probabilities of the searched object (the ‘‘target’’) being in various possible locations;

n

Corresponding author. Tel.: þ33 387315441; fax: þ33 387315309. E-mail address: [email protected] (H.A. Le Thi).

0305-0548/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cor.2012.07.006

The problem is to find the optimal distribution of this total effort that maximizes the probability of detection. Stone [21] summarized major steps in the development of search theory, for instance, stationary target problems, moving target problems, optimal searcher path algorithms, and dynamic search games. The last item (search games) is the primary focus of recent research, including numerous sub-domains such as mobile evaders, avoiding target, ambush games, inspection games, and tactical games. For moving target problems, decisive progress has been made in developing search strategies that maximize the probability of detecting the (moving) target within a fixed amount of time. In particular, Brown [2] and Washburn [24] have proposed an iterative algorithm in which the motion space and the time frame have been discretized, and the search effort available in each period is infinitely divisible between the grid cells of the target motion space. In this approach, the search effort available in each period is bounded above by a constant that does not depend on the allocations made during any other periods. In the ‘‘classical’’ search theory, the target is said to be detected if a detection occurs any time during the measurement epoch. Here, the target track will be said to be detected if (multiple) elementary detections occur at various times, and this is the fundamental difference. Thus, there is a test for acceptance (or detection) of a target track [5,23]. Track detection is also associated with a spatio-temporal modeling of the target track. A quite common hypothesis in classical search theory is that the objective function, which is generally defined as the non-detection probability in terms of the search efforts (variables), is a

232

H.A. Le Thi et al. / Computers & Operations Research 41 (2014) 231–239

convex function. This is quite a reasonable assumption in this context, which ensures convergence of the iterative algorithms. In this paper we consider an instance of search theory problems that can be described as follows: suppose that a space of search is partitioned into zones of reasonable size. A unique sensor must be able to explore efficiently a whole zone. Each zone is itself partitioned into cells. A cell is an area in which all points have the same properties, according to the difficulty of detection (altitude, vegetation, etc.). Each sensor has its own coefficient of visibility over a cell. The visibility coefficients depend also on the kind of the target that is searched. Here, there is a unique target to detect. The objective is to assign sensors to search zones and to find the best (minimum) resource shared by all the sensors so that the probability of detection of the target is maximum. The optimization model of this problem is hierarchical:

 At upper level: finding the best allotment of sensors to search zones (a sensor is allotted to a unique zone).

 At lower level: determining the best resource shared by all sensor, in order to have an optimal surveillance over the allotted zone. At the upper level, the objective function can be non-convex or implicitly defined via an algorithm applied to the lower level. This makes the problem very hard. In [20], Simonin et al. have proposed a hierarchical approach for solving this problem where a cross-entropy (CE) algorithm [1,4,19] has been developed for the upper level while an optimization method based on the algorithm of de Guenin [7] for detecting a stationary target has been used in the lower level. Besides this paper, we do not find in the literature any other work considering this problem. We can say that, up to now, there are no deterministic methods for solving this problem. In this work we develop a deterministic continuous approach based on DC (Difference of Convex functions) programming and DCA (DC optimization Algorithms). The contribution of the paper is 3-fold: Firstly, we propose in this paper an original mathematical model for the problem. By introducing the binary variables (the assignment variables) we formulate this hierarchical optimization problem as a mixed 0–1 nonlinear program. Due to the very large dimension of this problem in practice and the nonlinearity of the objective function, the standard methods in combinatorial optimization such as branch and bound, branch and cut, cutting plane cannot be applied. Attempting to develop robust numerical solution approaches, we try to reformulate the problem as a DC program. Secondly, we investigate an exact penalty technique to reformulate the mixed 0–1 nonlinear program into a continuous optimization problem that is in fact a DC program. We prove also that a penalty parameter can be estimated explicitly. Thirdly, we investigate DC programming and DCA for solving the related DC program. DC programming and DCA ([15,17,18] and references therein) aim to solve a general DC program of the form

a ¼ inff f ðxÞ :¼ gðxÞhðxÞ : x A Rp g ðPdc Þ,

ð1Þ

where g, h are lower semicontinuous proper convex functions on Rp . Such a function f is called a DC function, and gh, a DC decomposition of f while g and h are the DC components of f. The construction of DCA involves convex DC components g and h but not the function f itself: each iteration k of DCA consists of computing yk A @hðxk Þ,

xk þ 1 A arg minfgðxÞhðxk Þ/xxk ,yk S : x A Rp g

ðPk Þ:

Hence, for a DC program, each DC decomposition corresponds to a different version of DCA. Since a DC function f has an infinitely many DC decompositions which have crucial impacts on the qualities (speed of convergence, robustness, efficiency, globality of computed solutions, etc.) of DCA, the search for a ‘‘relevant’’ DC decomposition is important from an algorithmic point of view. Moreover, despite its local character, in practice DCA, with appropriately chosen initial points, converges quite often to global solutions. Finding a ‘‘good’’ initial point is then also an important stage of DCA. How to develop an efficient algorithm based on the generic DCA scheme for a practical problem is thus a judicious question to be studied, and the answer depends on the specific structure of the problem being considered. In the current work, using an appropriate DC decomposition we propose a DCA scheme that is very inexpensive in terms of CPU time thanks to the rapidity of the algorithm for solving the subproblem ðP k Þ. The paper is organized as follows. The problem statement is presented in Section 2. Section 3 introduces a new mathematical formulation of the considered problem, that is a mixed 0–1 nonlinear program. The solution method by DC programming and DCA via an exact penalty technique is investigated in Section 4. The computational results are reported in Section 5 while some conclusions and perspectives are discussed in Section 6.

2. Problem statement First, let us introduce the notation employed in the paper. E

space of search, Z: number of zones, S: number of sensors zone index cell index sensor index prior on the initial location of the target quantity of resource available for sensor s to search the space coefficient that characterizes the acuity of sensor s over cell i of the zone z (visibility coefficient).

z i s

a Fs wi,z,s

We state below the problem as described in [20] (see [20] for more details): The space of search: The search space, named E, is a large space with spatially variable search characteristics. The search space E is divided into Z search zones, denoted Ez ,z ¼ 1,2, . . . ,Z, each of them z is partitioned into Cz cells, denoted fci,z gCi ¼ so that: 1 Z [

E¼

Ez ,

z¼1 Cz [

Ez ¼

ci,z ,

Ez \ Ez0 ¼ |, 8z a z0 , ci,z \ cj,z ¼ |, 8ia j:

i¼1

A cell ci,z represents the smallest search area in which the search parameters are constant. For example, it can be a part of land with constant characteristics (latitude, landscape). Each zone must have a reasonable size in order to be explored by a sensor within a fixed time interval. The target: The target is hidden in one unit of the search space. Its location is characterized by a prior ai,z . Thus, we have Cz Z X X

ai,z ¼ 1:

z¼1i¼1

The means of search: Means of search can be passive (e.g. IRST, ESM) or active sensors (radars). We will consider that searching the target will be carried out by S sensors. Due to operational constraints, each sensor s A S must be allotted to a unique search

H.A. Le Thi et al. / Computers & Operations Research 41 (2014) 231–239

zone. For example, it could be the exploration time to share between units of a zone. At the lower level, for the sensor s allotted to the zone Ez , the amount of search resource allocated to the cell ci,z is denoted xi,z,s . It can represent the time spent on searching the cell ci,z (passive sensor), the intensity of emissions or the number of pulses (active sensors), etc. Furthermore, each sensor s has a search amount Fs , it means that if the sensor s is allotted to the zone Ez then Cz X

P The constraints Zz ¼ 1 uz,s ¼ 1, 8s ¼ 1,2, . . . ,S mean that each sensor s A S must be allotted to a unique search zone. With the usual model Ps [10]: Ps ðxi,z,s Þ ¼ expðwi,z,s xi,z,s Þ, the objective function f becomes f ðx,uÞ ¼

Cz Z X X

ai,z exp 

z¼1i¼1

xi,z,s r Fs :

To characterize the effectiveness of the search at the cell level, we consider the conditional non-detection probability Ps ðxi,z,s Þ which represents the probability of not detecting the target given that the target is hidden in ci,z and that we apply an elementary search effort xi,z,s on ci,z . Some assumptions are made to model Ps ðxi,z,s Þ. For all sensors, the function xi,z,s /P s ðxi,z,s Þ is convex and nonincreasing (law of diminishing return). Assuming independence of elementary detections, a usual model is P s ðxi,z,s Þ ¼ expðwi,z,s xi,z,s Þ, where wi,z,s is a (visibility) coefficient which characterizes the reward for the search effort put in ci,z by the sensor s. An additional assumption is that sensors act independently at the cell level. This means that if S sensors are allotted to ci,z the probability of not detecting a target hidden in ci,z is simply QS s ¼ 1 P s ðxi,z,s Þ. Let m : f1,2, . . . ,Sg-f1,2, . . . ,Zg be a mapping allotting sensors to search zones. Our aim is to find both the optimal mapping m (i.e. the best sensor-to-zone allotment) and the optimal local distributions x in order to minimize the non-detection probability. The objective function is then Cz Z X X z¼1i¼1

S X

! wi,z,s xi,z,s uz,s :

s¼1

Set

i¼1

f ðm,xÞ ¼

233

ai,z

Y

P s ðxi,z,s Þ,

n ¼ Z  S:

Let D and M be the sets defined by ( ) Cz X d xi,z,s r Fs ,z ¼ 1, . . . ,Z, s ¼ 1, . . . ,S , D ¼ x ¼ ðxi,z,s Þ A R þ : i¼1

( M¼

u ¼ ðuz,s Þ A f0,1g

ZS

:

Z X

) uz,s ¼ 1, s ¼ 1, . . . ,S :

z¼1

Then we can write the problem (3) in the form 8 ! Cz Z X S > X X > < min ai,z exp  wi,z,s xi,z,s uz,s x,u ðPÞ z ¼ 1i ¼ 1 s¼1 > > : s:t: x A D,u A M, which is a nonlinear mixed 0–1 programming problem. It is easy to see that the objective function of (P), say f, is convex in x for each fixed u, and similarly, it is convex in u for each fixed x. Moreover, f is infinitely differentiable.

4. Solution method

s A m1 ðzÞ

which leads to the following problem [20]: 8 Cz Z X Y X > > > f ðm,xÞ ¼ ai,z P s ðxi,z,s Þ min > > m,x > > z ¼ 1i ¼ 1 s A m1 ðzÞ > > > < Cz X xi,z,s r Fs , 8z, 8s A m1 ðzÞ, s:t: > > > i¼1 > > > > xi,z,s Z0, 8i A C z ,8s ¼ 1, . . . ,S, 8z ¼ 1, . . . ,Z, > > > : m mapping : f1,2, . . . ,Sg-f1,2, . . . ,Zg:

d ¼ S  ðC 1 þ C 2 þ    þ C z Þ,

We first reformulate (P) as a DC program by using an exact penalty technique. 4.1. DC reformulation ð2Þ

Consider the function p and the bounded polyhedral convex set K defined, respectively, by pðuÞ ¼

Z X S X

uz,s ð1uz,s Þ

z¼1s¼1

and (

3. A new mathematical formulation of problem (2)

K¼

Z:S

u ¼ ðuz,s Þ A ½0,1

:

Z X

) uz,s ¼ 1, s ¼ 1, . . . ,S :

z¼1

Let us introduce the allocation variable uz,s defined by ( 1 if the sensor s is allotted to the zone z, uz,s ¼ 0 otherwise:

We notice that p is finite and concave on RZS , non-negative on K and M ¼ fu A K : pðuÞ r 0g:

Let variable xi,z,s be the quantity of resource of the sensor s allotted to the cell ci,z in the zone z. We can rewrite (2) in the following form: 8 Cz Z X S Y X > > > > min ðf ðx,uÞ ¼ ai,z Ps ðxi,z,s uz,s ÞÞ > > x,u > > z ¼ 1i ¼ 1 s¼1 > > > > Cz X > > < s:t: xi,z,s r Fs , z ¼ 1, . . . ,Z,s ¼ 1, . . . ,S, i¼1 > > > Z > X > > > uz,s ¼ 1, s ¼ 1, . . . ,S, > > > > z¼1 > > > : uz,s A f0,1g, xi,z,s Z 0, s ¼ 1, . . . ,S, z ¼ 1, . . . ,Z, i ¼ 1, . . . ,C z : ð3Þ

Hence Problem (P) can be rewritten as

a ¼ minff ðx,uÞ : x A D, u A K, pðuÞ r 0g: The exact penalty result is given by the following theorem. Theorem 1. (i) Let 2

t 0 ¼ maxfJru ðf ðx,uÞÞJ : u A ½0,1n , x A Dg, then 8t 4t 0 the problem (P) minf f ðx,uÞ : x A D,u A Mg

ðPÞ

is equivalent to the problem minf f ðx,uÞ þtpðuÞ : x A D,u A Kg

ðPt Þ

234

H.A. Le Thi et al. / Computers & Operations Research 41 (2014) 231–239

in the following sense: they have the same optimal value and the same optimal solution set. (ii) If ðxn ,un Þ is a local minimizer to problem ðPt Þ then ðxn ,un Þ is a feasible solution to the problem (P).

Proposition 3. We have

Proof. Let jðuÞ :¼ minx A D ðf ðx,uÞ þ tpðuÞÞ. By definition of t0, for each fixed x in D, the function f ðx,uÞ þ t 0 pðuÞ is concave in u, thus the function minx A D ðf ðx,uÞ þt 0 pðuÞÞ is concave in u, and consequently, for any t 4t 0 , j is strongly concave in u. If ðxn ,un Þ is an optimal solution to the problem (P), then pðun Þ ¼ 0 and

where

f ðx,uÞ Z f ðxn ,un Þ,

8x A D, 8u A M:

ð4Þ

Hence f ðx,uÞ þ tpðuÞ Zf ðx,uÞ Z f ðxn ,un Þ þ tpðun Þ,

8x A D, 8u A M,

Since extrðKÞ ¼ M (see Lemma 2 below) and j is concave, we have minfjðuÞ : u A Kg ¼ minfjðuÞ : u A Mg:

ð6Þ

Combining (5) and (6) we obtain 8x A D, 8u A K:

ð7Þ

So ðx ,u Þ is also an optimal solution to the problem ðPt Þ. Conversely, suppose that ðxn ,un Þ is an optimal solution of the problem ðPt Þ, i.e. (7) holds. Thus jðuÞ Z jðun Þ,8u A K. Since j is strongly concave, un is an extreme point of K and therefore un A M. It means that ðxn ,un Þ is a feasible solution of the problem (P). This and (7) imply (4), i.e. ðxn ,un Þ is an optimal solution to (P). (ii) If ðxn ,un Þ is a local minimizer of problem ðPt Þ, then there exists a neighborhood V of xn , U of un , with U is convex and compact such that f ðx,uÞ þ tpðuÞ Z f ðxn ,un Þ þ tpðun Þ,8x A D \ V,8 u A K \ U. Let c be the function defined by n

cðuÞ ¼ minf f ðx,uÞ þ tpðuÞ : x A D \ Vg: It is clear that c is strongly concave in u A K \ U, and cðuÞ Z cðun Þ,8u A K \ U. So un is an extreme point of the convex set K \ U, thus un is an extreme point of the convex set K, say un A M. Therefore ðxn ,un Þ is a feasible solution to the problem (P). & Note that the part (ii) in this theorem is very useful for DCA applied to ðP t Þ because DCA usually produces a local minimizer. Lemma 2. We have extrðKÞ ¼ M, where extr(K) is the set of extreme points of K. Hence coðMÞ ¼ K. Proof. The set K can be expressed as the product of S simplex D, say K ¼D  D    D ( with D :¼

s

Fs ,

s¼1

a ¼ maxfai,z : z ¼ 1, . . . ,Z,i ¼ 1, . . . ,C z g, F ¼ maxfFs : 1, . . . ,Sg, w ¼ maxfwi,z,s : z ¼ 1, . . . ,Z,i ¼ 1, . . . ,C z ,s ¼ 1, . . . ,Sg: Proof. See Appendix. In the following sections we will investigate DCA for solving problem ðP t Þ with a sufficiently large value of t 0 r Faw2 PS s ¼ 1 Fs . 4.2. DC algorithm (DCA)

jðuÞ Z jðun Þ, 8u A M:

n

S X

ð5Þ

which implies that

f ðx,uÞ þ tpðuÞ Zf ðxn ,un Þ þ tpðun Þ,

t 0 :¼ maxfJr2u ðf ðx,uÞÞJ : u A ½0,1n ,x A Dg r Faw2

Z

u ¼ ðuz,s Þz A ½0,1 :

Z X

) uz,s ¼ 1 :

z¼1

It is easy to see that the set of extreme points of the simplex D is given by ( ) Z X s Z extrðDÞ ¼ u ¼ ðuz,s Þz A f0,1g : uz,s ¼ 1 : z¼1

Hence extrðKÞ ¼ extrðDÞ  extrðDÞ      extrðDÞ ¼ M:

&

Since D is compact and f A C 1 , maxfJr2u ðf ðx,uÞÞJ : u A ½0,1n , x A Dg exists. The following proposition gives an estimation of t 0 :

4.2.1. Outline of DC programming and DCA Pham Dinh Tao in 1985 introduced DC programming and DCA, which constitute the backbone of (smooth or nonsmooth) nonconvex programming and global optimization. These theoretical and algorithmic tools have been extensively developed, since 1994, by Le Thi Hoai An and Pham Dinh Tao to become now classic and increasingly popular (see [13–15,17,18], and references therein). DC programming and DCA address the problem of minimizing a function f which is a difference of convex functions on the whole space Rp or on a convex set C  Rp . Generally speaking, a DC program takes the form

a ¼ infff ðxÞ :¼ gðxÞhðxÞ : x A Rp g ðPdc Þ,

ð8Þ

where g, h are lower semicontinuous proper convex functions on

Rp . Such a function f is called DC function, and gh, DC decomposition of f while g and h are DC components of f. The convex constraint x A C can be incorporated in the objective function of ðP dc Þ by using the indicator function of C denoted wC which is defined by ( 0 if x A C, wC ðxÞ ¼ þ1 otherwise: If either g or h is a polyhedral convex function then ðP dc Þ is called a polyhedral DC program. Let g n ðyÞ :¼ supf/x,ySgðxÞ : x A Rp g be the conjugate function of g. Then, the following program is called the dual program of ðP dc Þ:

aD ¼ inffhn ðyÞg n ðyÞ : y A Rp g ðDdc Þ:

ð9Þ

One can prove that a ¼ aD (see e.g. [15,17]) and there is the perfect symmetry between primal and dual DC programs: the dual to ðDdc Þ is exactly ðP dc Þ. DCA is based on the local optimality conditions of ðP dc Þ, namely @hðxn Þ \ @gðxn Þ a|

ð10Þ

and |a @hðxn Þ  @gðxn Þ:

ð11Þ

A point x satisfying (10) is called a critical point of gh, or (10) is the generalized KKT condition for the DC program ðP dc Þ. The condition (11) is necessary for local optimality of ðP dc Þ. It is also sufficient for many classes of DC programs quite often encountered in practice (see e.g. [15,17]). The idea of DCA is simple: each iteration of DCA approximates the concave part h by its affine majorization (that corresponds to taking yk A @hðxk Þ) and minimizes the resulting convex function (that is equivalent to determining xk þ 1 A @g n ðyk ÞÞ. n

H.A. Le Thi et al. / Computers & Operations Research 41 (2014) 231–239

Generic DCA scheme Initialization: Let x0 A Rp be a best guess, 0’k. Repeat Calculate yk A @hðxk Þ Calculate xk þ 1 Aarg minfgðxÞhðxk Þ/xxk ,yk S : x A Rp g ðPk Þ kþ 1’k Until convergence of xk .

235

and nr o min JuJ2 /u,vk S : u A K : 2

ð14Þ

We are now in a position to summarize the DCA for solving Problem ðPt Þ Step1. Initialization: let ðx0 ,u0 Þ satisfy the constraints of the problem. Choose E1 4 0, E2 40 and k¼0. Step2. Compute ðyk ,vk Þ ¼ rhðxk ,uk Þ, with

It is important to mention the following main convergence properties of DCA:

yk ¼ rxk rx f ðxk ,uk Þ,

 DCA is a descent method (the sequences fgðxk Þhðxk Þg and

Step3. Compute ðxk þ 1 ,uk þ 1 Þ by solving the two convex quadratic problems (13) and (14). Step4. Iterate Steps 2 and 3 until

fh ðyk Þg n ðyk Þg are decreasing) without linesearch. If the optimal value a of the problem ðP dc Þ is finite and the infinite sequences fxk g and fyk g are bounded then every limit point xn (resp. yn ) of the sequence fxk g (resp. fyk gÞ is a critical n point of gh (resp. h g n Þ. DCA has a linear convergence for general DC programs. DCA has a finite convergence for polyhedral DC programs. n



 

It is worth noting that the general DCA scheme for solving general DC programs is rather a philosophy than an algorithm. In fact, there is not only one DCA but infinitely many DCAs for a considered DC program. DCA’s distinctive feature relies upon the fact that DCA deals with the convex DC components g and h but not with the DC function f itself. This fact is crucial for nonconvex nonsmooth programs for which DCA is one of the rare effective algorithms. The solution of a practical nonconvex program by DCA must have two stages: the search of an appropriate DC decomposition and the search of a good initial point. An appropriate DC decomposition, in our sense, is the one that corresponds to a DCA which is not expensive and has interesting convergence properties. In the next subsection we will develop an instance of DCA to solve the problem ðPt Þ, and study the convergence properties of the proposed algorithm.

s¼1

where N ¼ maxfC z : z ¼ 1, . . . ,Zg: PS

s¼1

Fs þ Nwag,

minfgðx,uÞhðx,uÞ : ðx,uÞ A D  Kg,

9ðghÞðxk þ 1 ,uk þ 1 ÞðghÞðxk ,uk Þ9 r E1 ð1þ 9ðghÞðxk þ 1 ,uk þ 1 Þ9Þ or Jðxk þ 1 ,uk þ 1 Þðxk ,uk ÞJ1 r E2 ð1 þJðxk þ 1 ,uk þ 1 ÞJ1 Þ:

Theorem 4 (Convergence properties of Algorithm DCA. For simplicity’s sake, we omit the dual part of these properties). (i) DCA generates the sequence fðxk ,uk Þg such that the sequence fðghÞðxk ,uk Þg is decreasing convergent. (ii) The sequence fðxk ,uk Þg converges to a KKT point for the problem (12). Proof. Since the sequence fðxk ,uk Þg is bounded, the DCA’s convergence for general DC programs [15,17] states that the sequence ff ðxk ,uk Þ þtpðuk ÞÞg is decreasingly convergent and every limit point of fðxk ,uk Þg is a generalized KKT point for problem (12). Property (ii) concerning the convergence of the whole sequence fðxk ,uk Þg to a KKT point for problem (12) is a stronger result for DC programs with subanalytic data [16]. & 4.3. Efficient ways for solving problems (13) and (14)

4.2.2. Description of the DCA applied to ðP t Þ From the computations in the Appendix we have ( ) S X JHðf ÞJ1 r max Saw2 þSaw2 F þ aw,2Faw2 Fs þNwa ,

Hence, let r :¼ maxfSaw2 þ Saw2 F þ aw,2Faw2 a DC formulation of (Pt) can be

vk ¼ ruk ru f ðxk ,uk Þ þ tð2uk eÞ:

ð12Þ

where

In fact, we can continuously divide the problems (13) and (14) into some smaller problems with unique constraints (except positive constraints). By denoting yk ¼ ðyi,z,s Þ A Rd ,vk ¼ ðvz,s Þ A RZS , solving the problem (13) is replaced by solving S  Z subproblems with unique inequality constraints 8 ! Cz X > 2yki,z,s > 2 > > min xi,z,s  x > > x r i,z,s > > i,z,s i ¼ 1 < Cz X ð15Þ > xi,z,s r Fs , s:t: > > > > i¼1 > > > : xi,z,s Z 0, 8i ¼ 1,2, . . . ,C z ,

DCA applied to DC program (12) consists of computing, at each iteration k, the two sequences fðyk ,vk Þg and fðxk ,uk Þg such that ðyk ,vk Þ A @hðxk ,uk Þ and ðxk þ 1 ,yk þ 1 Þ is an optimal solution of the following convex quadratic program: nr o Jðx,uÞJ2 /ðx,uÞ,ðyk ,vk ÞS : ðx,uÞ A D  K , min 2

and solving the problem (14) is replaced by solving S subproblems with unique equality constraints ! 8 Z X 2vkz,s > 2 > > min u  u z,s > z,s > uz,s r > > z¼1 < Z X ð16Þ > uz,s ¼ 1, > > s:t: > > z¼1 > > : uz,s Z 0, 8z ¼ 1,2, . . . ,Z:

which can be decomposed into two convex quadratic ones nr o JxJ2 /x,yk S : x A D min 2

The problem (16) can be solved efficiently by the Block Pivotal Principal Pivoting Algorithm [9]. We also adapt this algorithm for solving the problem (15) as follows.

gðx,uÞ :¼

r 2

Jðx,uÞJ2

and

hðx,uÞ :¼

r 2

Jðx,uÞJ2 f ðx,uÞtpðuÞ:

ð13Þ

236

H.A. Le Thi et al. / Computers & Operations Research 41 (2014) 231–239

Consider the convex quadratic problem of the form (15): 8  n  X 1 2 > > > xi ai xi > min > x 2 > > i¼1 < n X ð17Þ > xi r F, > > s:t: > > i¼1 > > : x Z 0, 8i ¼ 1,2, . . . ,n: i

Set I ¼ fi A f1,2, . . . ,ng : ai 4 0g: P We firstly observe that if i A I ai r F then xn ¼ ðxn1 , . . . ,xnn Þ, where n xi ¼ ai if iA I, 0 otherwise, is an optimal solution to Problem (17). In contrast, i.e. X ai 4 F, ð18Þ iAI

then it is not difficult to show that the optimal solution of (17) belongs to the following set: n X

xi ¼ F:

Table 1 Parameters for dataset 1.

i¼1

Therefore, we can use the Block Pivotal Principal Pivoting Algorithm for solving (17).

5. Numerical experimentation The search space is the lake of Laouzas in France [20] (Fig. 1). We consider two sets of data:

 The first dataset: the search space is divided into Z¼4 zones



Fig. 2. Partition 1 of the lake Laouzas.

and there are n¼30 cells in each zone (see Fig. 2). We have six sensors. All sensors have the same amount of resource F. The coefficients are given in Table 1. The second dataset: the search space is divided into Z ¼6 zones and there are n ¼80 cells in each zone (see Fig. 3). We have 10 sensors. All sensors have the same amount of resource F. The coefficients are given in Table 2.

Type of cell Prior of target

Forest Water Plat land Rough land Very rough land Town

Visibility of sensors Sensor 1

Sensor 2

Sensor 3

Sensor 4

Sensor 5

Sensor 6

0.0085 0.001 0.0115 0.013 0.014

0.4 0.9 0.3 0.2 0.1

0.5 0.1 0.1 0.7 0.6

0.6 0.1 0.4 0.8 0.7

0.8 0.1 0.6 0.2 0.1

0.5 0.3 0.5 0.4 0.3

0.1 0.5 0.2 0.6 0.5

0

0.8

0.9

0.1

0.7

0.6

0.2

Fig. 3. Partition 2 of the lake Laouzas.

Fig. 1. An aerial photograph of the lake of Laouzas.

The program is written in C using Microsoft Visual Cþþ 2008, and implemented on a notebook with chipset Intel(R) Core(TM) Duo CPU 2.0 GHz, RAM 3 GB.

H.A. Le Thi et al. / Computers & Operations Research 41 (2014) 231–239

Table 2 Parameters for dataset 2.

Forest Water Plat land Rough land Very rough land Town

Prior of target

Table 3 DCA and CE: dataset 1. Visibility of sensors

Amount of resource

S1

S2

S3

S4

S5

S6

S7

S8

S9

S10

0.002125 0.000250 0.002875 0.003250 0.003500

0.4 0.9 0.3 0.2 0.1

0.5 0.1 0.1 0.7 0.6

0.6 0.1 0.4 0.8 0.7

0.8 0.1 0.6 0.2 0.1

0.5 0.3 0.5 0.4 0.3

0.1 0.5 0.2 0.6 0.5

0.4 0.3 0.7 0.2 0.2

0.5 0.5 0.5 0.5 0.3

0.6 0.4 0.2 0.6 0.6

0.9 0.7 0.4 0.1 0.4

0.000000

0.8 0.9 0.1 0.7 0.6 0.2 0.9 0.4 0.6 0.3

We compare our DCA with the existing algorithm for this problem, the CE algorithm developed in [20]. The main idea of this algorithm is to generate particular allotments of sensors to search zones that will be evaluated and then selected, in order to obtain a drawing law which will converge toward the optimal allotment. The steps of this CE algorithm can be described as follows: Step1: Initialize M ¼ M 0 ¼ ðpM0 ðz9sÞÞ a uniform distribution, i.e. pM0 ðz9sÞ ¼

1 , Z

s ¼ 1, . . . ,S, z ¼ 1, . . . ,Z,

and choose y A ð0,1Þ. (Note that pðz9sÞ represents the probability to assign sensor s to zone z.) Step2: Generate N allotment mappings m1 ,m2 , . . . ,mN according to M, and compute f ðmk Þ :¼ minx f ðmk ,xÞ,k ¼ 1,2, . . . ,N, where f ðm,xÞ is objective function of (2). (In other words, we obtain f ðmk Þ by solving Problem (2) when mk is known. It is a convex programming problem.) Step3: Sort the sequence ff ðmk ÞgN k ¼ 1 in the increasing orders. Let f ðmsð1Þ Þ rf ðmsð2Þ Þ r    rf ðmsðNÞ Þ, where s is a permutation of the set f1,2, . . . ,Ng. Set T ¼ byNc, then choose T best draws X sð1Þ ,X sð2Þ , . . . ,X sðTÞ . Step4: Update M by the formula sðkÞ

pM ðz9sÞ :¼

DCA method

CE method

F

cardfk A fsð1Þ, sð2Þ, . . . , sðTÞg : xs T

¼ zg

:

Step5: Iterate steps 2–4 until convergence. The starting point ðx0 ,u0 Þ of DCA, where x0 ¼ ðx0i,z,s Þ, u0 ¼ ðu0z,s Þ are chosen as follows: x0i,z,s ¼ F=n, u0z,s ¼ 1=Z,s ¼ 1, . . . ,S,z ¼ 1, . . . , Z,i ¼ 1, . . . ,C z . The parameters r and t are dynamically adjusted during DCA’s iterations. From an algorithmic point of view, the smaller the r and t are, the more efficient the DCA are (because the smaller the concave part of the objective function of Problem (Pt ) is, the better its affine approximation is). Hence it would be efficient to choose the values t and r as small as possible. In our experiments, we start DCA with a quite small value of r and t (depending on the parameters F), and then increase both parameters r and t by 1% at each iteration of DCA. The choice of the first value of r and t must ensure that during the algorithm all the values r and t are smaller than their overestimation given in Section 4. We take E1 ¼ E2 ¼ 107 . Since the CE method is heuristic, for each parameter F we run this algorithm 10 times and then take the average for the 10 results. The parameters for CE are as follows: for each iteration, y ¼ 0:3, the number of samples is N ¼50 for the first dataset, and is N ¼100 for the second dataset. The number of iterations is limited to 15. We observe from Tables 3 and 4 that the DCA produced better solutions than the CE in almost all cases while DCA’s CPU time is shorter. Especially, in the case of large dimension (the second dataset), the DCA is much faster than the CE: the ratio of CPU time

F¼5 F ¼ 10 F ¼ 15 F ¼ 20 F ¼ 25 F ¼ 30 F ¼ 35 F ¼ 40

Objective function

Time (s)

Objective function

Time (s)

0.818115 0.684438 0.579880 0.496526 0.422119 0.363833 0.323616 0.277450

0.328 0.266 0.334 0.672 0.484 0.703 0.797 0.828

0.818577 0.684701 0.580673 0.500550 0.428900 0.372219 0.327574 0.292742

1.192 1.427 1.925 2.284 2.689 2.620 2.285 3.688

Table 4 DCA and CE: dataset 2. Amount of resource

DCA method

CE method

F

F¼5 F ¼ 10 F ¼ 15 F ¼ 20 F ¼ 25 F ¼ 30 F ¼ 35 F ¼ 40

Objective function

Time (s)

Objective function

Time (s)

0.906248 0.831474 0.766634 0.706112 0.652908 0.610872 0.564653 0.524267

2.227 1.072 1.311 1.508 1.816 1.917 2.114 2.284

0.905917 0.831506 0.766327 0.707495 0.655548 0.608630 0.565490 0.528006

9.587 22.532 25.420 31.006 35.894 45.529 36.396 56.692

1 0.9 0.8 Non−detection probability

Type of cell

237

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

10

20

30

40

50

60

70

80

90

100 110 120

Amount of resource Fig. 4. Relation of non-detection probability given by DCA and amount of resource.

between the two algorithms varies from 4.30 to 24.82. That is because the CE method requires a large computation at lower level. Fig. 4 shows the results (non-detection probability) obtained by DCA on the first dataset when the amount of resource, F, varies. We see that there is a negative relationship between F and the non-detection probability. When F ¼ 110, the non-detection probability is 0.060488.

238

H.A. Le Thi et al. / Computers & Operations Research 41 (2014) 231–239

6. Conclusion

þ9ai1 ,z1 wi1 ,z1 ,s1 9Þexp 

We have presented a new approach for solving the problem of planning a multisensor multizone search for a target. This is the first time a deterministic optimization approach is investigated in the literature for solving this problem. This constitutes an interesting contribution of the paper. For solving the combinatorial optimization problem by DCA, an innovative continuous approach in nonconvex programming framework, we first reformulated the original problem in the form of a DC program. We prove an exact penalty result in which the penalty parameter can be estimated. That is our second important contribution. The third contribution deals with the development of an efficient DCA scheme for solving the resulting problem. In this scheme, the subconvex programs is decomposed on some simpler convex programs that can be solved by an efficient algorithm. Numerical results show that our approach can solve this problem more efficiently than the CE method, an efficient algorithm for the problem. Whenever there is a good starting point, we can get a better local solution or even a global solution. In a future work, we will investigate DCA for solving the moving target case of this problem.

S X

! wi1 ,z1 ,s xi1 ,z1 ,s uz1 ,s

s¼1 S X

¼

9ai1 ,z1 wi1 ,z1 ,s1 uz1 ,s1 wi1 ,z1 ,s2 uz1 ,s2 9

s2 ¼ 1

þ

S X

9ai1 ,z1 wi1 ,z1 ,s1 uz1 ,s1 wi1 ,z1 ,s3 xi1 ,z1 ,s3 9

s3 ¼ 1

! S X  wi1 ,z1 ,s xi1 ,z1 ,s uz1 ,s : 9ai1 ,z1 wi1 ,z1 ,s1 9 exp  s¼1

Let

a ¼ maxfai,z : z ¼ 1,2, . . . ,Z,i ¼ 1,2, . . . ,C z g, F ¼ maxfFs : 1,2, . . . ,Sg, w ¼ maxfwi,z,s : z ¼ 1,2, . . . ,Z,i ¼ 1,2, . . . ,C z ,s ¼ 1,2, . . . ,Sg, N ¼ maxfC z : z ¼ 1,2, . . . ,Zg: Because 0 r uz,s r 1,8s A S,z A Z, uz1 ,s Þ r 1, we have

and

expð

PS

s¼1

wi1 ,z1 ,s xi1 ,z1 ,s

Si1 ,z1 ,s1 rSaw2 þ Saw2 F þ aw ¼ Saw2 ðF þ 1Þ þ aw, 8z1 ¼ 1,2, . . . ,Z,i1 ¼ 1,2, . . . ,C z1 ,s1 ¼ 1,2, . . . ,S: We have ! C z1 S X X @f ðx,uÞ ¼  ai,z1 wi,z1 ,s1 xi,z1 ,s1 exp  wi,z1 ,s xi,z1 ,s uz1 ,s , @uz1 ,s1 s¼1 i¼1

Appendix A. Estimate the norm of Hessian matrix of function f ðx,uÞ and maxfJr2u ðf ðx,uÞÞJ1 : u A ½0,1n ,x A Dg x ¼ ðx1,1,1 , . . . ,xC Z ,Z,S Þ A Rp ,

u ¼ ðu1,1 , . . . ,uZ,S Þ A RZS :

We have @f @xi1 ,z1 ,s1

ðx,uÞ ¼ ai1 ,z1 wi1 ,z1 ,s1 uz1 ,s1 exp 

S X

! wi1 ,z1 ,s xi1 ,z1 ,s uz1 ,s ,

s¼1

if z2 a z1 ,

if z2 ¼ z1 :

s¼1

2

@ f ðx,uÞ @xi1 ,z1 ,s1 @xi2 ,z2 ,s2 8 0 > > > > < ai1 ,z1 wi1 ,z1 ,s1 uz1 ,s1 wi1 ,z1 ,s2 uz1 ,s2 ! ¼ S X > > >  exp  w x u > i1 ,z1 ,s i1 ,z1 ,s z1 ,s :

@2 f ðx,uÞ @uz1 ,s1 @uz2 ,s2 8 0 > > > C > z1 > > X > > < ai,z1 wi,z1 ,s1 xi,z1 ,s1 wi,z1 ,s2 xi,z1 ,s2 ¼ i¼1 ! > > S > X > > >  exp  w x u > z ,s i,z ,s i,z ,s 1 1 1 > :

if ði2 ,z2 Þ aði1 ,z1 Þ,

@2 f ðx,uÞ @uz1 ,s1 @xi2 ,z2 ,s2 8 0 > > > > > a w x w u > > > i2 ,z1 i2 ,z1 ,s1 i2 ,z1 ,s1 i2 ,z1 ,s2 z1 ,s2 ! > > S > X > > <  exp  wi2 ,z1 ,s xi2 ,z1 ,s uz1 ,s ¼ s¼1 > > > ðwi2 ,z1 ,s1 ai2 ,z1 þ ai2 ,z1 wi2 ,z1 ,s1 xi2 ,z1 ,s1 wi2 ,z1 ,s1 uz1 ,s1 Þ > > ! > > S > X > > >  exp  wi2 ,z1 ,s xi2 ,z1 ,s uz1 ,s > > : s¼1

if ði2 ,z2 Þ ¼ ði1 ,z1 Þ,

s¼1

@2 f ðx,uÞ @xi1 ,z1 ,s1 @uz3 ,s3 8 0 > > > > > ai ,z wi ,z ,s uz1 ,s1 wi1 ,z1 ,s3 xi1 ,z1 ,s3 > > ! > 1 1 1 1 1 > S > X > > > wi1 ,z1 ,s xi1 ,z1 ,s uz1 ,s <  exp  s¼1 ¼ > > > ðai1 ,z1 wi1 ,z1 ,s1 þ ai1 ,z1 w2i ,z ,s uz1 ,s1 xi1 ,z1 ,s1 Þ > 1 1 1 > > ! > > S > X > > > w x u  exp  z ,s i ,z ,s i ,z ,s > 1 1 1 1 1 :

if z2 a z1 ,

if z2 ¼ z1 ,s2 a s1 , 8i2 A z1 ,

if z2 ¼ z1 ,s2 ¼ s1 , 8i2 Az1 :

if z3 a z1 , On the line ðz1 ,s1 Þ of the Hessian matrix of the function f, we have the sum of absolute of all elements:

if z3 ¼ z1 ,s3 as1 ,

Sz1 ,s1 ¼

C z1 S X X

9ai,z1 wi,z1 ,s1 xi,z1 ,s1 wi,z1 ,s2 xi,z1 ,s2 9

s2 ¼ 1 i ¼ 1

if z3 ¼ z1 ,s3 ¼ s1 :

s¼1

exp 

S X

! wi2 ,z1 ,s xi2 ,z1 ,s uz1 ,s

s¼1

On the line ði1 ,z1 ,s1 Þ of the matrix Hessian of the function f, we have the sum of absolute of all elements: Si1 ,z1 ,s1 ¼

S X

9ai1 ,z1 wi1 ,z1 ,s1 uz1 ,s1 wi1 ,z1 ,s2 uz1 ,s2 9exp 

s2 ¼ 1

þ

S X s3 ¼ 1

S X

þ

C z1 S X X

! wi1 ,z1 ,s xi1 ,z1 ,s uz1 ,s

exp 

s¼1

ð9ai1 ,z1 wi1 ,z1 ,s1 uz1 ,s1 wi1 ,z1 ,s3 xi1 ,z1 ,s3 9

9ai2 ,z1 wi2 ,z1 ,s1 xi2 ,z1 ,s1 wi2 ,z1 ,s2 xi2 ,z1 ,s2 9

s 2 ¼ 1 i2 ¼ 1 S X

! wi2 ,z1 ,s xi2 ,z1 ,s uz1 ,s

s¼1

þ

C z1 X i2 ¼ 1

9wi2 ,z1 ,s1 ai2 ,z1 9exp 

S X s¼1

! wi2 ,z1 ,s xi2 ,z1 ,s uz1 ,s

H.A. Le Thi et al. / Computers & Operations Research 41 (2014) 231–239

r

C z1 S X X

Faw2 xi,z1 ,s2 þ

s2 ¼ 1 i ¼ 1

r 2Faw2

C z1 S X X

Faw2 xi2 ,z1 ,s2 þNwa

s2 ¼ 1 i2 ¼ 1

S X

Fs þ Nwa:

s¼1

The infinity norm of Hessian matrix of function f satisfies ( ) S X 2 2 2 Fs þNwa ¼ r: JHðf ÞJ1 r max Saw þSaw F þ aw,2Faw s¼1

On the line ðz1 ,s1 Þ of the Hessian matrix ru f ðx,uÞ of the function f ðx,Þ, we have the sum of absolute of all elements T z1 ,s1 ¼

C z1 S X X

9ai,z1 wi,z1 ,s1 xi,z1 ,s1 wi,z1 ,s2 xi,z1 ,s2 9

s2 ¼ 1 i ¼ 1

exp 

S X

! wi2 ,z1 ,s xi2 ,z1 ,s uz1 ,s

s¼1 C z1 S X X

r

s2 ¼ 1 i ¼ 1

Faw2 xi,z1 ,s2 r Faw2

S X

Fs , 8x A D, 8u A ½0,1n :

s¼1

Thus, we have supfJr2u ðf ðx,uÞÞJ : u A ½0,1n ,x A Dg r Faw2

S X

Fs :

s¼1

References [1] deBoer PT, Kroese DP, Mannor S, Rubinstein RY. A tutorial on the crossentropy method. Annals of Operations Research 2005;134:19–67. [2] Brown SS. Optimal search for a moving target in discrete time and space. Operations Research 1979;32(5):1107–15. [3] Cadre JP, Souris G. Searching tracks. IEEE Transactions on Aerospace and Electronic Systems 2000;36(4):1149–66. [4] Costa A, Jones OD, Kroese D. Convergence properties of the cross-entropy method for discrete optimization. Operations Research Letters 2007;35(5): 573–80.

239

[5] Dobbie JM. Transfer of detection contacts to tracking contacts in surveillance. Operations Research 1966;14:791–800. [6] Frost JR. Principles of search theory, Part III: probability density distributions. Response 1999;17(3):1–10. [7] deGuenin J. Optimum distribution of effort: an extension of the Koopman theory. Operations Research 1961;9(1):1–7. [8] Haley KB, Stone LD, editors. Search theory and applications. New York: Plenum Press; 1980. [9] Ju´dice J, Raydan M, et al. On the solution of the symmetric eigenvalue complementarity problem by the spectral projected gradient algorithm. Numerical Algorithms 2008;47(4):391–407. [10] Koopman BO. Search and its optimization. Mathematical Monthly 1979;7: 527–40. [11] Koopman BO. Search and screening: general principles with historical applications. New York: Pergamon Press; 1980. [12] Koopman BO. Search and screening: general principle with historical applications. Alexandria, VA: MORS Heritage Series; 1999. [13] Le Thi HA. DC programming and DCA. Available on /http://lita.sciences. univ-metz.fr/lethi/DCA.htmlS. [14] Le Thi HA, Pham Dinh T. Large scale molecular optimization from distance matrices by a DC optimization approach. SIAM Journal on Optimization 2003; 14(1):77–116. [15] Le Thi HA, Pham Dinh T. The DC (difference of convex functions) programming and DCA revisited with DC models of real world nonconvex optimization problems. Annals of Operations Research 2005;133:23–46. [16] Le Thi HA, Huynh VN, Pham Dinh T. Convergence analysis of DCA for DC programming with subanalytic data. Research Report, LITA, University of Paul Verlaine Metz; 2010. [17] Pham Dinh T, Le Thi HA. Convex analysis approach to DC programming: theory, algorithms and applications (dedicated to Professor Hoang Tuy on the occasion of his 70th birthday). Acta Mathematica Vietnamica 1997;22: 289–355. [18] Pham Dinh T, Le Thi HA. DC optimization algorithms for solving the trust region subproblem. SIAM Journal on Optimization 1998;8:476–505. [19] Rubinstein RY, Kroese D. The cross-entropy method: a unified approach to combinatorial optimization, Monte´ Carlo simulation, and machine learning. Berlin: Springer; 2004. [20] Simonin C, Le Cadre JP, Dambreville F. A hierarchical approach for planning a multisensor multizone search for a moving target. Computers & Operations Research 2009;36(7):2179–92. [21] Stone LD. What’s happened in search theory since the 1975 Lanchester prize? Operations Research 1989;37(3):501–6. [22] Stone LD. Theory of optimal search, 2nd ed. Arlington, VA: Operations Research Society of America. ORSA Books; 1989. [23] Van Keuk G. Sequential track extraction. IEEE Transactions on Aerospace and Electronic Systems 1998;34(4):1135–48. [24] Washburn AR. Search for a moving target, The FAB algorithm. Operations Research 1983;31(4):739–51.