- Email: [email protected]

Contents lists available at ScienceDirect

Automatica journal homepage: www.elsevier.com/locate/automatica

Optimal path planning and sensor placement for mobile target detection✩ Bomin Jiang a,b , Adrian N. Bishop b,c , Brian D.O. Anderson a,b , Samuel P. Drake d,e a

Research School of Engineering, Australian National University, Canberra, Australia

b

National ICT Australia (NICTA), Canberra, Australia

c

Centre for Health Technologies, University of Technology Sydney (UTS), Sydney, Australia

d

Defence Science and Technology Organisation, Edinburgh, Australia

e

School of Chemistry and Physics, The University of Adelaide, Adelaide, Australia

article

info

Article history: Received 14 June 2014 Received in revised form 29 March 2015 Accepted 18 June 2015

Keywords: Path planning Guidance navigation and control Perception and sensing Mobile robots Detectors

abstract For a flying military vehicle, avoiding detection can be a key objective. To achieve this, flying the leastprobability-of-detection path from A to B through a field of detectors is a fundamental strategy. While most of the previous optimization models aim to minimize the cumulative radar exposure, this paper derives a model that can directly minimize the probability of being detected. Furthermore, a variational dynamic programming method is applied to this model which allows one to find a precise locally optimal path with low computational complexity, even when Doppler effects increase the dimension of this problem. In addition, a homotopy method with exceptionally low computational complexity is derived to adjust the optimal path when the detection rate function changes due to the removal of detectors, the addition of detectors or other changes of detectors. Finally, the paper also shows how to apply a convex optimization method to find optimal positions of detectors when vehicles can do path planning. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction For flying military vehicles, avoiding detection is often a key objective. There has been much prior research leading to many different sophisticated strategies to avoid detection. Maithripala, Maithripala, and Jayasuriya (2008) explained how to deceive a ground radar network into seeing a spurious phantom track in its radar space by intercepting the pulse, introducing a time delay and

✩ B. Jiang, A.N. Bishop and B.D.O. Anderson are supported by NICTA, which is funded by the Australian Government as represented by the Department of Broadband, Communication and the Digital Economy and the Australian Research Council through the ICT Centre of Excellence program. A.N.Bishop is also supported by the Australian Research Council (ARC) via a Discovery Early Career Researcher Award (DE-120102873) and B.D.O. Anderson is also supported by Australian Research Council under grant DP110100538 and DP130103610. A.N. Bishop is also an Adjunct Fellow at the Australian National University (ANU). The material in this paper was partially presented at the 19th IFAC World Congress, August 24–29, 2014, Cape Town, South Africa. This paper was recommended for publication in revised form by Associate Editor Shun-ichi Azuma under the direction of Editor Toshiharu Sugie. E-mail addresses: [email protected] (B. Jiang), [email protected] (A.N. Bishop), [email protected] (B.D.O. Anderson), sam.d[email protected] (S.P. Drake).

http://dx.doi.org/10.1016/j.automatica.2015.07.007 0005-1098/© 2015 Elsevier Ltd. All rights reserved.

re-transmitting the radar’s pulses. Liu and Guo (2002) treat the use of radar stealth material to absorb radar pulses. For a detection minimizing problem, which is the focus of this paper, it is an obvious strategy to fly the least-probability-ofdetection path from the source to the destination through a field of detectors; this strategy has been studied with increasing intensity as vehicles have become increasingly intelligent. Much research emphasizes the modelling aspect of this problem. Examples of such work are Kim and Hespanha (2003), Murphey, Uryasev, and Zabarankin (2004) and Zabarankin, Uryasev, and Murphey (2006), in which risk functions are employed to represent the cumulative radar exposure when flying through the field of detectors. The objective of these models is to minimize the cumulative radar exposure. Other research focuses on path planning algorithms that generate optimal routes to minimize the cumulative radar exposure. Examples of such research are Royset, Carlyle, and Wood (2009) and Zabarankin et al. (2006), in which the problem is abstracted to be a weight constrained shortest path problem (WCSPP) and then solved. In the above papers, the overall objective is to minimize the cumulative radar exposure instead of the probability of being detected. These two objectives may intuitively seem to be the same but they are not precisely equivalent. For example, if a mobile

128

B. Jiang et al. / Automatica 60 (2015) 127–139

vehicle enters the non-escape zone of a radar, the model may still try to minimize the cumulative radar exposure which will ultimately make no difference in the probability of being detected (which always equals 1). Again, if a mobile vehicle is very far away from the detector and the signal to noise ratio is very low, minimizing the radar exposure is also of little relevance. As an example of a model whose objective is to minimize the probability of being detected, the paper (Royset et al., 2009) modelled radars using threat circles and solved the problem using Lagrangian relaxation plus (near-shortest path) enumeration (LRE). However, due to the idealized detection model and the high complexity of the algorithm, the method can only provide a very rough optimal path. This paper considers the problem where the objective is to minimize the probability of being detected rather than cumulative radar exposure. It is shown in literature that minimizing the probability of detection is closer to reality (Mahafza, 2013). While the majority of radars use a pulse-by-pulse mechanism, we assume that the radar repetition frequency is high enough so that the detection event can be considered as an event occurring in continuous time. As will be seen, the derived cost function has an identical form but different underlying meaning to those in cumulative exposure models. Consequently, most of the numerical algorithms are still applicable with a partial adjustment on the typical detection rate function. In part contrast to the high computational complexity of most previous numerical algorithms, this paper utilizes a variational dynamic programming approach (Barraquand & Ferbach, 1993) which can obtain accurate local optima. Furthermore, a homotopy method is derived to adjust the optimal path when the detection rate function changes due to the removal of detectors, the addition of detectors, the changes of understanding of detectors, and so on. The underlying principle of this homotopy method is to linearize the change of optimal path using the first order calculus of variations. The linearized model can be found with several efficient matrix operations; further, because the major coefficient matrix is a tridiagonal matrix, the computational complexity of this approach is exceptionally low. As a result, even miniature unmanned aerial vehicles (UAVs) such as quadrotors can deploy this method using on-board calculation. The novel contributions of the paper are as follows (a) using a detection rate function analogous to a Poisson process rate which is closer to reality, (b) applying a gradient based variational dynamic programming to the optimal path problem, which can obtain solutions when Doppler effect increases the dimension of the problem, and (c) deriving a homotopy method which can adjust the optimal path with extremely low computational complexity when the detection rate function changes. In comparison to the conference version of this paper (Jiang, Bishop, Anderson, & Drake, 2014), the optimization model in this paper takes into account the Doppler effect. Furthermore, the problem has been solved from the detectors’ perspective in addition to the vehicle perspective: a convex optimization method is applied to find optimal positions of detectors when vehicles can do path planning. 2. Problem statement using probability of detection measure Throughout this paper, we only discuss the first detection of radars or sensors; sensor tracking problem is not in the scope of this paper. We also assume that detection events in nonoverlapping time intervals are independent. Consider a closed bounded region C in R2 . Within this region, it is understood that there is a moving agent or target, and a set of detectors. The target wishes to remain undetected, and the detectors are not perfect. The target is assumed to be moving with constant (unit) speed between

two fixed points in C along a simple (i.e. non-self-intersecting), piecewise smooth curve. As the detectors are imperfect, associated with every point (x, y) ∈ C is a nonnegative value p(x, y), with the interpretation that the probability that a target located at the point (x, y) would be detected in a time interval of length ∆t is p(x, y)∆t + O (∆t ). The function p is assumed to be smooth. Detection events in nonoverlapping time intervals are independent. For a target stationary at (x, y), the assumptions made so far mean that the probability of first detection in a time interval of length T is 1 − exp[−p(x, y)T ]. The problem of interest is then to determine, given two fixed points in C and some function p(x, y), an optimal trajectory for the target between these two points such that the detection probability is minimized. 2.1. Computing the probability of detection for a path Consider a parameterized path (x(s), y(s)) of length L from (x0 , y0 ) to (xT , yT ). For a small positive ∆, determine the points (xi , yi ) = (x(i∆), y(i∆)) such that the unit speed assumption implies a sequence of ∆ spaced points (xi , yi ) along the path. The probability of non-detection for the target is the probability that there is no detection in each interval with end points (xi , yi ) and (xi+1 , yi+1 ), i = 0, 1, . . . , (L − 1)/∆. The probability that there is a detection in the ith interval is p(xi , yi )∆ + o(∆) and so the probability that there is no detection is simply [1 − p(xi , yi )∆ + o(∆)]. Since detection events are independent through time, the probability that there is no detection along the whole path is (L−1)/∆

Pno (x(·), y(·)) = Πi=1

[1 − p(xi , yi )∆ + o(∆)] .

(1)

Evidently, (L− 1)/∆

ln Pno (x(·), y(·)) =

ln[1 − p(xi , yi )∆ + o(∆)].

(2)

i =1

Letting ∆ go to zero, we see that the probability of detection, call it Pdet , obeys ln(1 − Pdet ) = ln Pno (x(·), y(·))

= − lim

∆→0

=− Q

(L− 1)/∆

p(xi , yi )∆

i=1

p(x(s), y(s))ds

(3)

where Q is a path starting from the point (x0 , y0 ) and ends at the point (xT , yT ). Note the integration is effectively a parameterized path integral. Equivalently, we thus have ñ ô Pdet = 1 − exp −

Q

p(x(s), y(s))ds .

(4)

Evidently, it is straightforward to compute the probability of detection for a prescribed trajectory and knowledge of p(x, y). However, the problem of interest requires the computation of a trajectory that minimizes the detection probability given two fixed points (x0 , y0 ), (xT , yT ) and knowledge of the detection probability function p(x, y). This is a calculus of variations problem which is equivalent to finding a given trajectory that minimizes the path integral p(x(s), y(s))ds. (5) Q

If we view (5) in isolation then the problem formulation resembles the widely considered risk function approach to path planning where p(x, y) could be regarded as a cumulative radar exposure per unit time; see Kim and Hespanha (2003), Murphey et al. (2004) and

B. Jiang et al. / Automatica 60 (2015) 127–139

Zabarankin et al. (2006). However, the interpretation as presented here is based on a derivation starting from an infinitesimal detection rate p(x, y) where the probability of detecting a target at (x, y) for one second is not p(x, y) but rather 1 − exp[−p(x, y)]. We believe the idea of defining an infinitesimal detection rate in this sense is more natural. Indeed, the function p(x, y) can be thought of as an intensity function for a Poisson process defined at each point (x, y) in the sense that for a stationary target at (x, y) the time of first detection could be thought of as a nonnegative random variable with distribution function 1 − exp[−p(x, y)T ] over T . In addition to detection events being independent across non-overlapping time intervals, it also follows that detection events for fixed targets at distinct points in R2 are also independent. Along any constant velocity path the intensity p(x, y) as a function of the path (and consequently time) essentially defines a non-homogeneous Poisson process (Franceschetti & Meester, 2007, p. 9). The typical detection rate function is given in Section 2.3 and that model will be used in the rest of this paper unless otherwise noted. Note that when deriving (5), the constant speed constraint (throughout taken to be unity) has been taken into account. A more convenient form of (5) for our purposes is given by

129

thresholding and depends naturally on the signal-to-noise ratio of the receiver. For a pulsed radar, let fr denote the pulse repetition frequency (alternatively 1/fr could correspond to some time period in a continuous-wave radar). Then the probability of detection during one pulse is 1 − exp[−p(x, y)/fr ]. Writing this in terms of the received signal statistics we get Ç å p(x, y) 1 − exp − fr Ç å ∞ 1 δ√ θ − a2 /δ 2 = I0 θ exp − dθ (8) a 2 vt 2

where τ ∈ [0, 1] is a given parameterization obeying (x(0), y(0)) = (x0 , y0 ) and (x(1), y(1)) = (xT , yT ) and ′ denotes ddτ . The objective is to find the function (x(·), y(·)) defined on [0, 1] that minimizes (6). Throughout we assume that the target trajectory is piecewise smooth as requiring total smoothness rules out paths with an isolated corner, and such paths may be optimal when compared to a totally smooth path.

where the term being integrated is just the probability density of the received signal strength using a common radar receiver model; see Mahafza (2013). In particular, vt is a chosen detection threshold, a is the signal amplitude at p(x, y), δ is the noise amplitude and I0 is the modified Bessel function of the first kind and zeroth order. For a typical radar, the signal-to-noise ratio a2 /δ 2 is proportional to k/ρ 4 (for active radars) or k/ρ 2 (for passive sensors) where ρ is the range between the radar and the target and k depends on the radar power, antenna gain etc. The choice of vt depends on the acceptable probability of false alarm denoted by Pfa ; i.e. the acceptable probability that the signal magnitude exceeds the threshold value when noise alone is present (Mahafza, 2013). Suppose SNR is the signal-to-noise ratio. The probability of detection in one can then be approximated by the Marcum √ pulse» q-function Q [ 2SNR, −2ln(Pfa )], which appears often in such detection problems; see Mahafza (2013) and Simon (1998). Examples of the probability of detection in one pulse are given in Fig. 1.

2.2. The connection to the Euler–Lagrange equations

2.4. The detection probability with multiple detectors

It is well known that optimization problems of the form of (6) can be treated classically by Euler–Lagrange equations, see for example (Goldstein, 1980, ch. 3.8). There is a clear parallel between the path of minimum detection as described by Eqs. (5) and (6) and Hamilton’s principle in classical mechanics. Hamilton’s principle is stated as (Goldstein, 1980, ch. 2). The motion of a system from time t1 to time t2 is such that the action functional t2 I = (7) L(qi , q˙i , t )dt

Consider a collection p1 (x, y), p2 (x, y), . . . , pi (x, y) of detection rate functions associated with, for example, multiple detectors. Then the relevant value of p(x, y) is simply the sum so long as the individual pi (x, y) allow one to define the probability of detection on an infinitesimal interval δ t and each detector operates independently (i.e. detection events occur independently). More formally, the probability that detector i does not detect the target is 1 − pi δ t and the probability that all detectors fail to detect the target is (1 − pi δ t ) owing to the independence assumption. Then the probability of detection is 1 − (1 − pi δ t ) = ( pi )δ t + o(δ t ).

0

1

p(x(τ ), y(τ )) x′τ 2 + y′τ 2 dτ

(6)

t1

has a stationary value for the correct path of motion, where qi s are the coordinates, the dot denotes the time derivative q˙i = dqi /dt, L is known as the Lagrangian and I is the action. Using Hamilton’s principle and the calculus of variation one can derive d dt

∂L ∂ q˙i

∂L ∂ qi

− = 0. By comparison the Euler–Lagrange equations with (6) we see that the problem of solving the optimal path reduces to solving the above Euler–Lagrange equations with L = »

p(x(τ ), y(τ )) + As a further connection with classical physics, the path in accordance with Fermat’s principle is also analogous to a light ray passing through an inhomogeneous medium with refractive index p(x, y). x′ 2

3. Incorporating the Doppler effect in radar probabilistic models 3.1. How the Doppler effect can make a difference in the models

y′ 2 .

2.3. A typical detection rate We briefly outline the relationship of our infinitesimal detection rate p(x, y) to a common signal detection probability function used in radar design. Target detection in practice is based on

When detecting a moving vehicle (for example, a car), other objects (for example, a road) will reflect some signal (also called clutter), which could significantly raise the background interference power level. However, as shown in Fig. 2, the reflection signal of the moving vehicle has some Doppler shift while the clutter signal does not. We assume that we use a stationary pulse Doppler radar in this scenario. Evidently, a Doppler shift can be used to potentially reduce the interference power level at the receiving frequency and thus increase the signal to interference power ratio (Skolnik, 2008, p. 4.2).

130

B. Jiang et al. / Automatica 60 (2015) 127–139

(a) Active radar.

(b) Passive sensor. Fig. 1. Examples of the probability of detection in one pulse.

3.3. Problem formulation when considering the Doppler effect When considering the Doppler effect, the probability of detection of a vehicle is dependent on the vehicle’s radial velocity with respect to each detector. In the first case, we suppose that a vehicle travels at a constant speed. Then the only factor affecting the Doppler shift is the direction of the velocity vector. Thus (6) becomes 1 p x(τ ), y(τ ), φ(x′τ , y′τ ) x′τ 2 + y′τ 2 dτ (9) 0

Fig. 2. Noise, clutter and target signal in a pulse Doppler radar.

3.2. How to incorporate the Doppler effect in radar probabilistic model Recall in a pulsed radar, the probability of detection is dependent on the signal to noise ratio SNR. When considering Doppler effect and clutter, the probability of detection is dependent on the signal to clutter and noise ratio SCNR (McMillan & Kohlberg, 2010). Suppose there is a stationary radar with known transmitting signal frequency in a environment with known clutter spectrum and white noise. Then the SCNR only depend on received signal strength and the Doppler shift caused by radial motion of the target with respect to the stationary detector. Thus there exists aSCNR = fc (aSNR , vr ). We introduce the function fc (aSNR , vr ) because the total interference power equals system noise plus clutter (Skolnik, 2008). With one can still use the Marcum q-function » √ this modification, Q [ 2SCNR, −2ln(Pfa )] in Section 2.3 to estimate the probability of being detected. The specific form of fc (aSNR , vr ) depends on how clutter signals are modelled; a reasonable setting is that the clutter is Gaussian distributed with mean at zero frequency shift, as shown in Mahafza (2013) and Skolnik (2008, p. 2.11). When the background noise dominates in interference power, for example, a ground radar looking into the sky, a Doppler shift will make no difference to the signal to clutter and noise ratio SCNR. That is because we assume that the background noise is white noise and uniformly distributed on frequency domain. On the other hand, when the clutter dominates in interference power, for example, a radar looking at the sea, a Doppler shift will reduce the SCNR Skolnik (2008, p. 4.48).

where φ(x′τ , y′τ ) denotes the direction of the vector (x′τ , y′τ ). In the second case, we suppose that the vehicle’ speed can change within a range. We consider this case because a vehicle does not necessarily need to travel at its highest speed to avoid detection. Although a higher speed can result in a shorter time exposure within a detection environment, it will also cause a higher probability of detection per unit time. As a result, a new variable vs (τ ), which is defined as the absolute value of the velocity of a vehicle at position x(τ ), y(τ ) , is introduced into the model. The integral becomes 0

1

1 · x′τ 2 + y′τ 2 dτ p x(τ ), y(τ ), φ(x′τ , y′τ ), vs ·

vs

(10)

and the optimization goal is to minimize the above integral. Note (9) can be regarded as a special case of (10). 3.4. Solution using variational dynamic programming Analytical solution to the stated optimization problem is not generally possible so we resort to numerical techniques. In this section, we will show how to use the gradient-based Variational Dynamic Programming method to obtain numerical solution taking into account the Doppler effect. First we consider the constant speed case as shown in (9). Consider a continuously parameterized target trajectory (x(τ ), y(τ )) from (x0 , y0 ) to (xT , yT ) in some compact set C ⊂ R2 . For a positive ∆, determine the points (xi , yi ) = (x(i∆), y(i∆)) along (x(τ ), y(τ )). Suppose ∆ is chosen such that there are N + 1 such points starting at (x0 , y0 ) and ending at (xN , yN ) = (xT , yT ). If we approximate (x(s), y(s)) by a straight line on each interval between (xi , yi ) and (xi+1 , yi+1 ) then the target path will be a piecewise smooth polygonal curve. Denote the cost to move from (xa , ya ) to (xb , yb ) by Ja,b . Then, under the discretization

B. Jiang et al. / Automatica 60 (2015) 127–139

131

assumptions thus far

1

0

N p x(τ ), y(τ ), φ(x′τ , y′τ ) · x′τ 2 + y′τ 2 dτ ≈ Ji−1,i

(11)

i=1

and the incremental cost to move from (xn−1 , yn−1 ) to (xn , yn ) is

Jn−1,n =

n N n−1 N

p(x, y, φ(x′τ , y′τ )) ·

x′τ 2 + y′τ 2 dτ .

(12)

Due to the straight line assumption, the above incremental cost is computable given (xn−1 , yn−1 ) and (xn , yn ). Suppose then that Ji∗,j is the optimal cost when the target moves from (xi , yi ) to (xj , yj ). Since the value of (x0 , y0 ) is prescribed it follows that J0∗,n is only a function of (xn , yn ). Then

J0∗,n =

min {J0∗,n−1 + Jn−1,n },

xn−1 ,yn−1

1 ≤ n ≤ N.

(13)

Define a grid on R2 as the embedding of a (countably infinite) square lattice graph G(V , E ) with vertices V and edges E. The spacing of the grid points (embedded vertices of G) defines the resolution of the grid. Now we have enough equations to state a variational dynamic programming solution to the given optimization problem. The broad idea of the algorithm is to guess a path, then consider variations to it based on local perturbations and when a reduction in the cost function is achieved, select the adjusted path. Let (x∗1 , y∗1 ), (x∗2 , y∗2 ), . . ., (x∗N −1 , y∗N −1 ) denote the optimal waypoints, then the algorithm is given by Algorithm 1. Because the probability of detection is always positive, the gradient algorithm must converge.

and if we relax the approximation that each interval is a straight line then the incremental cost to move from (xn−1 , yn−1 , vs,n−1 ) to (xn , yn , vs,n ) is

Algorithm 1 Variational Dynamic Programming

Jn−1,n =

Fig. 3. Comparison of the optimal path with or without Doppler effect.

1

0

1 p x(τ ), y(τ ), φ(x′τ , y′τ ), vs · · x′τ 2 + y′τ 2 dτ

≈

vs

N

1: 2: 3: 4: 5: 6: 7: 8: 9:

10: 11: 12: 13:

Select a check point number N Select a recursive loop number M Generate an initial path as near as possible to the optimal path represented by check points (xn , yn ) ∈ V where 0 ≤ n ≤ N. for i = 1; i ≤ M ; i = i + 1 do for n = 1, n ≤ N , n = n + 1 do Generate a finite set of points Cn whose elements include (xn , yn ) and a set of neighbour points of (xn , yn ) in G. end for for n = 2, n ≤ N , n = n + 1 do For each (xn , yn ) ∈ Cn , choose (xn−1 , yn−1 ) ∈ Cn−1 to minimize the function J0∗,n−1 + Jn−1,n . Then J0∗,n can be expressed as a function of (xn , yn ) only. end for Use the destination coordinates (xN , yN ) to find points (x∗N −1 , y∗N −1 ), (x∗N −2 , y∗N −2 ), . . . , (x∗1 , y∗1 ). Assign the value of (x∗1 , y∗1 ), (x∗2 , y∗2 ), . . . , (x∗N , y∗N ) to ordered pairs (x1 , y1 ), (x2 , y2 ), . . . , (xN , yN ) accordingly. end for

In the above algorithm, when generating the set Cn , a set of neighbour points was introduced. This set of neighbour points are points connected with the point (xn , yn ) in the graph G. In practice, G does not need to be a rectangular grid; a smarter way to choose G usually results in less computation. When considering the case where a vehicle can adjust its speed, the dimension of the problem will increase. Analogously to the setting in the above, we consider a continuously parameterized target trajectory (x(τ ), y(τ ), vs (τ )) from (x0 , y0 , vs0 ) to (xT , yT , vsT ) in some compact set C ⊂ R3 . Note vs (t ) denotes the speed at time t. Denote the cost to move from (xa , ya , vsa ) to (xb , yb , vsb ) by Ja,b . Then, under the discretization assumptions we have

(14)

Ji−1,i

i =1

n N n−1 N

p(x, y, φ(x′τ , y′τ ), vs ) ·

1

vs

·

x′τ 2 + y′τ 2 dτ .

(15)

The rest of the VDP algorithm remains the same except for the fact that now we are solving this problem in 3D rather than 2D. The extra dimension will not make the computation much more difficult if a gradient based method is used. However, it is clear that a global optimum will be difficult to obtain as the dimension of the problem increases. The gradient-based algorithm can only achieve local a minimum, therefore the generation of the initial path (i.e. the initial guess) has a critical impact on the final result. In the case of incorporating Doppler effect, the optimal path without considering Doppler effect can serve as a reasonable initial path. 3.5. Illustrative examples Suppose that a vehicle intends to fly from (−3, 3) to (5, 5) and it has a constant speed. We also assume that the Doppler shift of the clutter signal is normally distributed with mean of zero frequency shift (Mahafza, 2013). Fig. 3 shows the simulation result when there is a detector at the origin. The dash line shows the optimal path without considering the Doppler effect and the solid line shows the optimal path considering the Doppler effect. When the Doppler effect is taken into account, the vehicle intends to reduce the radial movement; although its path is closer to the detector, it can still hide within the clutter region in the frequency domain to reduce the probability of being detected. When considering the case where a vehicle can adjust its speed, the dimension of the problem will increase. Suppose the vehicle can choose its speed in [0 100], suppose further that initial and final velocities of the vehicle are both 0. Fig. 4(a) is the optimal path and Fig. 4(b) shows the speed of the vehicle along the path. When the vehicle speed is high, the vehicle tends to move along the tangential direction of the detector to reduce the Doppler shift.

132

B. Jiang et al. / Automatica 60 (2015) 127–139

(a) Optimal path.

(b) Velocity graph. Fig. 4. Optimal path and velocity considering Doppler effect and velocity change.

4. A homotopy method for fast trajectory computation One of the main contributions of this paper is a computationally fast homotopy method for deriving a modified optimal path given a change in the total detection probability rate p(x, y) and the previously derived optimal path. This method is practically important since the environment in which the target moves may be dynamic; e.g. the number of detectors may change as new detectors are added/discovered, detectors may be removed or switched off or the detector function at particular detectors may evolve, etc. In all such cases it could be very time consuming to re-run the optimal trajectory planning algorithm and obtain the complete path from the beginning. More formally, we suppose one already has an optimal solution to (6) (or (5)), for some p(x, y), expressed by x(τ ) and y(τ ). Given the addition/removal of detectors or change in detection rate at some detectors, we suppose the detection rate is changed to p(x, y) + q(x, y). The homotopy method introduced in this section allows one to modify the optimal path to the updated detection rate with rapid modification if the change is small. The computation of the updated path can be accomplished by iteratively performing certain matrix operations and the computational complexity is considerably less than Murphey et al. (2004).

where C is some constant. Now if the first detector is changed but remains unmoved, one can use a homotopy method to quickly modify its path. Suppose for a certain detection rate p(ρ) we know the optimal path ρ(θ ), θ ∈ [θ0 , θT ]. Then consider a small variation and let the detection rate function be p(ρ) + λq(ρ) where q(ρ) is a bounded function and λ is a sufficiently small number. Correspondingly, a small variation to the optimal path is expected and we model this as

ρ(θ ) + λσ (θ ) where σ (θ ) is also bounded. Furthermore, the constant on the right hand side of the Beltrami identity will become C + λγ where γ is a finite number. Now we have 2 p(ρ) + λq(ρ) · ρ(θ ) + λσ (θ ) … 2 = C + λγ . 2

ρ(θ ) + λσ (θ )

A special case will be considered first due to its analytic solvability. Furthermore, it can also serve as introduction to more complicated general cases. The special case is when the probability rate of detection is constant on any circle whose centre is the detector and there is only one such detector. Now set up a polar coordinate system whose origin is at the detector. Suppose ρ is the radial coordinate and θ is the angle. The probability of detection is only a function of ρ . If the desired initial and terminal points for the trajectory are on the same radial line then the optimal path is just a straight line. Otherwise, suppose θ0 and θT define the polar angles of the desired initial and terminal points. Then (5) becomes Ã å Ç θT dρ(θ ) 2 dθ (16) p(ρ(θ )) ρ 2 (θ ) + dθ θ0

ρ 2 (θ) +

dθ

(18)

dθ

where p(ρ), q(ρ) and ρ(θ ) are known functions and C is a known constant which can be obtained from (17). Now we have a linear differential equation in σ (θ ) a1 (θ )σ (θ ) + a2 (θ )

dσ dθ

+ a3 (θ )γ + a4 (θ ) = 0

(20) pρ 2

dρ

where a1 = ρ p(2 − 2 ρ dρ 2 ), a2 (θ ) = − 2 ddρθ 2 , a3 (θ ) = ρ +( dθ ) ρ +( dθ ) dρ 2 2 2 − ρ + ( dθ ) and a3 (θ ) = ρ q. The solution to this first order linear differential equation is θ ρ 2 + 2( dρ )2 dθ σ (θ ) = exp d θ dρ 2

θ0

ρ dθ

Ñ Ä ä2 é θ θ ρ 2 + 2 ddρθ · exp − dθ g (θ )dθ dρ θ0

from which, using the Beltrami identity, see Fox (1950), we find where (17)

) + λ dσd(θ θ

Therefore, on subtracting (17) from (18), rearranging, and eliminating the terms in higher order of λ we have Ã Ä äÄ ä Ç å ρσ + ddρθ ddσθ dρ 2 2 2 + + γ 2pρσ + ρ q = C … ρ (19) Ä ä2 dθ dρ

ρ2 +

4.1. The homotopy approach: A special case

p(ρ(θ ))ρ 2 (θ ) … Ä ä =C dρ(θ) 2

dρ(θ ) dθ

+

g (θ ) = −

ρ dθ

θ0

»

γ ρ 2 + ρ ′ 2 − ρ 2 q ρ 2 + ( ddρθ )2

pρ 2 ρ ′

(21)

.

(22)

B. Jiang et al. / Automatica 60 (2015) 127–139

The initial condition is σ (θ0 ) = 0 while σ (θT ) = 0 from which the constant γ can be obtained. It is not clear a priori whether (21) is bounded. Simulation result shows that (21) is bounded for many choices of q(θ ) but unbounded for others. Often (21) is unbounded when the detection rate falls off so fast that the optimal path is at infinity.

σx (τ ) and σy (τ ) can be obtained. The detail will be discussed in Section 4.4. As with the special case we have applied a linearization approximation and thus an error could occur in the form

(p + λq)(x′′ + σx′′ ) −

4.2. The homotopy method: The general case

+

Recall the equation 1

0

p(x(τ ), y(τ )) x′τ 2 + y′τ 2 dτ

(23)

∂p ′ 2 ∂p ′ ′ y + x y =0 ∂x τ ∂y τ τ

(24)

∂p ′ 2 ∂p ′ ′ x + x y = 0. ∂y τ ∂x τ τ

(25)

(p + λq)(y′′ + σy′′ ) −

where q(x, y) is a bounded function and λ is a sufficiently small number. Correspondingly, a small variation to the optimal path is expected x(τ ) + λσx (τ ) and y(τ ) + λσy (τ ) where σx (τ ) and σy (τ ) are also bounded. By following the same operations as with the special case we obtain the linear differential equations regarding σx (τ ) and σy (τ ) which are pσx ′′ + qx′′ − 2

(26)

and

−

∂p ′ ′ x σx ∂y

∂ q ′2 ∂ q ′ ′ ∂ p ′ ′ ∂ p ′ ′ x + xy + x σy + σx y = 0. ∂y ∂x ∂x ∂x

(27)

In the above equations, p(x, y), q(x, y), x(τ ) and y(τ ) are known functions. Consequently, we have a system of second order linear differential equations with two unknown functions σx (.) and σy (.) a1 (τ )σx ′′ (τ ) + b11 (τ )σx ′ (τ ) + b12 (τ )σy ′ (τ ) + c1 (τ ) = 0 a2 (τ )σy ′′ (τ ) + b21 (τ )σx ′ (τ ) + b22 (τ )σy ′ (τ ) + c2 (τ ) = 0.

∂(p + λq) ′ (x + σx′ )2 ∂y

∂(p + λq) ′ (x + σx′ )(y′ + σy′ ) = ϵ2 (τ ) (30) ∂x from which one can calculate the error ϵ1 (τ ) and ϵ2 (τ ). In order to adjust the obtained σx and σy and diminish ϵ1 (τ ) and ϵ2 (τ ), let +

x(τ ) + λσx (τ ) → x(τ ) + λσx (τ ) + ξx (τ )

Again following the same line of reasoning as with the special case we can derive now two differential equations regarding ξx (τ ) and

ξy (τ )

∂(p + λq) ′ (y + σy′ )ξy′ ∂x ∂(p + λq) ′ ∂(p + λq) ′ (x + σx′ )ξy′ + (y + σy′ )ξx′ (31) + ∂y ∂y

−ϵ1 (τ ) = (p + λq)ξx′′ − 2

p(x, y) + λq(x, y)

pσy ′′ + qy′′ − 2

(29)

y(τ ) + λσy (τ ) → y(τ ) + λσy (τ ) + ξy (τ ).

Now suppose for a certain detection rate p(x, y), we already know the optimal path expressed by x(τ ) and y(τ ) (e.g. from the variational method derived previously), and consider a small variation to the detection rate function given by

∂p ′ ′ y σy ∂x ∂q ∂q ∂p ∂p − y′2 + x′ y′ + x′ σy ′ + σx ′ y′ = 0 ∂x ∂y ∂y ∂y

∂(p + λq) ′ (x + σx′ )(y′ + σy′ ) = ϵ1 (τ ) ∂y

and

and p(x, y)y′′τ −

∂(p + λq) ′ (y + σy′ )2 ∂x

and

from (6). Recall we are seeking the path (x(τ ), y(τ )) such that (23) is minimized given an initial and terminal point for the trajectory and p(x, y). The Euler–Lagrange equations applied to (23) are p(x, y)x′′τ −

133

(28)

Because the above equations are second order differential equations, there are two constants of integration associated with each of σx (τ ) and σy (τ ). Correspondingly, σx (τ ) has two boundary conditions σx (0) = σx (1) = 0 while σy (τ ) also has two boundary conditions σy (0) = σy (1) = 0. Now the numerical solutions of

and

∂(p + λq) ′ (x + σx′ )ξx′ ∂y ∂(p + λq) ′ ∂(p + λq) ′ (yy + σy′ )ξx′ + (x + σx′ )ξy′ . + ∂x ∂x

− ϵ2 (τ ) = (p + λq)ξy′′ − 2

(32)

These equations have a similar form to (26) and (27) and one can then obtain ξx (τ ) and ξy (τ ) following the same operations used to find σx (τ ) and σy (τ ). 4.3. An iterative algorithm for trajectory updating The discussion thus far in this section allows us to compute the variation of the optimal path λσx (τ ) and λσy (τ ) corresponding to a small variation of the detection rate function λq(x, y). Suppose now there is a bounded variation q(x, y) occurring in the detection rate function p(x, y) and we want to find the optimal path when the detection rate function is p(x, y) + q(x, y). We construct a parameterized function H (x, y, Λ) = p(x, y) + Λq(x, y)

(33)

where Λ ∈ [0, 1]. Note the optimal path corresponding to H (x, y, 0) is known to us while the optimal path corresponding to H (x, y, 1) is desired. We then discretize the interval according to 0 = Λ0 < Λ1 < Λ2 < · · · < Λn = 1 where each increment Λi+1 − Λi (i = 0, 1, 2, . . . , n − 1) is a sufficient small number that can be thought of as corresponding to λ in the previous subsections. Consequently, the original problem is equivalent to solving this sequence of differential equations just previously derived. During each iteration, the objective is to find the optimal path with the detection rate H (x, y, Λi+1 ) and knowledge of the optimal path corresponding to the case H (x, y, Λi ). Note the homotopy method derived in this paper is not identical with the homotopy perturbation method in He (2003).

134

B. Jiang et al. / Automatica 60 (2015) 127–139

The homotopy perturbation method derived by He looked at series expansion around Λ = 0. It assumes that higher order terms of Λ are small and therefore in his method, convergence of the solution is not guaranteed. On the other hand, in our proposed homotopy method, we design λ in the manner of choosing the step size in numerical integration or Euler method. For detailed explanation see Section 4.5. ˆ as the Remark 1. Suppose a vehicle has already obtained curve AB optimal path from A to B, and then after the vehicle has travelled from A to an intermediate point C , it finds an additional detector. Now it can use the homotopy method to modify the optimal path. ı rather than AB ˆ as a whole because the Note that it should modify CB vehicle is not at A now; C should be considered as its new origin.

Fig. 5. Potential discontinuity in the homotopy method.

4.4. Numeric solution During each iteration the problem boils down to solving a system of second order linear differential equations with two unknown functions. When analytic solutions are not available, numeric methods can be used. The work in Keller (1992) shows a general approach for numerically solving second order ordinary linear differential equations with time complexity O(N ) where N is the dimension of the unknown vector. In our particular case, N is the number of waypoints. 4.5. Selection of λ When implementing the homotopy method, it is critical to choose the right λ in each step. Generally speaking, choosing λ is like choosing the step size in numerical integration or Euler method. Furthermore, since we know that the optimal path should fulfil the Euler–Lagrange equation (24) and (25), we can define an error indicator η ≥ 0 such that 1 Ç 2 ∂ p ′ 2 ∂ p ′ ′ 2 p(x, y)x′′ τ − y + x y η x(τ ), y(τ ) = ∂x τ ∂y τ τ 0 2 å ∂p 2 ∂p + p(x, y)y′′ τ − x′τ + x′τ y′τ dτ (34) ∂y ∂x as the two norm of the error of the Euler–Lagrange equation. This error indicator can be used to guide selection of λ given a desired level of accuracy and computational budget. If a large error indicator shows up, one should reduce the value of λ to achieve satisfactory accuracy; if a small error indicator shows up, one should increase the value of λ to save computation. An example appear below in Section 4.7. 4.6. Caution with discontinuities The change in the optimal path is unfortunately not always continuous in Λ which means a small variation in the detection rate may result in an abruptly large change to the optimal path. In such cases the homotopy method will fail to provide a good result. For example, suppose there are two detectors in R2 located symmetrically with respect to the line segment between the desired initial and terminal points for the trajectory; e.g. as shown by Fig. 5. At first, the target knows p(x, y) corresponding just to Detector 1 and computes an optimal path shown by Line 1. Then Detector 2 with detection rate q(x, y) is introduced. Further suppose that the detectors are of the same type but that Detector 2 has a much higher sensitivity than Detector 1. Now the optimal path corresponding to H (x, y, 0) is Line 1. When Λ increases from zero, the optimal path will be pushed down as shown by Line 2. However, because q(x, y) will ultimately dominate H (x, y, 1), the optimal path corresponding to H (x, y, 1) should be Line 3. This analysis leads to us to the conclusion that the optimal path will have to jump at some value of Λ.

Fig. 6. Simulation result of Section 4.1.

4.7. Illustrative examples and computational comparison We suppose there is a target in R2 that needs to travel from (−5, 5) to (5, 5) and there is a radar placed at the origin. The flying vehicle has knowledge about that detector and has previously computed an optimal path accordingly. The probability rate function in this section is as Section 2.3. Now a new active radar is added to the origin and the target then intends to quickly adjust its optimal path. As shown in Fig. 6, the dotted line is the original optimal path with only one detector; the solid line is the path found after adjusting the original path using the homotopy method and the dash line is the optimal path with both detectors found using variational dynamic programming. Now a new passive sensor is added to the point (6, 6) and the vehicle intends to quickly adjust its optimal path. As shown in Fig. 7, the dotted line is the optimal path before adjustment; the solid line is the path after adjustment using the homotopy method and the dash line is the optimal path generated using variational dynamic programming. In the above simulation, no adjustment using (31) and (32) is made in any steps of the homotopy method. Fig. 8 shows the result when adjustments are made. Now we give the computation comparison between firstly, the Network Flow optimization as proposed in Murphey et al. (2004), secondly, the Variational Dynamic Programming method as in Section 3.4 and thirdly, the homotopy method in this section. Suppose kn is the number of elements in Cn in Variational Dynamic programming (VDP), N and M are the checkpoint numbers and iteration numbers respectively in the VDP and homotopy methods. Similar with Murphey et al. (2004), we divide the computation time into two parts: the preprocessing time spent on calculating risk along each line segment, and the optimization time spent on carrying out the algorithm. Here are our observations:

B. Jiang et al. / Automatica 60 (2015) 127–139

135

5. Sensor positioning for detecting intelligent vehicles To this point of the paper, path planning to minimize detection has been viewed from the vehicle’s perspective rather than detector’s perspective. In this section, we are going to study the sensor placement problem where the objective now is to find sensor positions that maximize the probability of vehicle detection when vehicles can plan their path. In the literature generally, there are a number of papers discussing the problem of sensor placement for various applications (Barrett, 2007; Bishop, Fidan, Anderson, Dogancay, & Pathirana, 2010; Dhillon & Chakrabarty, 2003); e.g. optimal sensor placement for target localization. Nevertheless, to the best of our knowledge, the particular sensor placement problem in this section has not been previously studied. Suppose pi (x, y) is the detection rate function of the ith detector, which is defined in previous sections. Suppose further that one agent fly along a path Q from a certain source to a destination. The probability of being detected by multiple detectors is given by Ç å Pdet = 1 − exp − pi (x(t ), y(t ))dt . (35)

Fig. 7. Simulation result of Section 4.2.

Q

This expression holds when the positions of detectors are fixed. In the sensor positioning problem, the positions of detectors are variables, therefore the expression needs to be modified Ç å Pdet = 1 − exp − pi (x(t ), y(t ), Dxi , Dyi )dt (36) Q

where (Dxi , Dyi ) denotes the position of the ith detector. Our goal is to find optimal solutions (Dx∗i , Dy∗i ) of the optimization problem ´´ ® ® pi (x(t ), y(t ), Dxi , Dyi )dt . (37) max min Dxi ,Dyi

Fig. 8. Simulation result of Section 4.2. Table 1 Computation comparison with different λ. Value of λ

Error indicator η

Runtime/sec

0.2 0.1 0.02

1.4 × 10 3.3 × 10−16 5.3 × 10−19

0.025 0.047 0.206

−4

1. Network Flow optimization: preprocessing time is required to calculate the risk on each edge; numerical integration is required in this step. Optimization runtime depends on what algorithm is used to solve the Network Flow optimization problem. 2. Variational Dynamic Programming: preprocessing time is required to calculate the risk on each line segment; numerical integration is also required in this step. The operation count in variational dynamic programming is O(M · N0 kn ). 3. Homotopy method: During preprocessing, the risk on each waypoint is calculated, thus numerical integration is not required. The operation count is O(M · N ). Table 1 shows how the value of λ will affect the runtime and the value of the error indicator η. We use the same example as Fig. 7 and set the number of way points N as 100. In the above example, we assign λ = Λi − Λi−1 uniformly for comparison purposes. In this case only, the iteration number M = 1/λ.

x(t ),y(t )

Q

We assume that the vehicle will optimize its own path using the path planning model in Section 2. For each sensor placement setting, there is a corresponding optimal path, which gives a lower bound of the probability of detection for all admissible paths. The sensor placement problem is trying to maximize this lower bound. Note that the sensors do not have knowledge of the vehicle’s actual path. The vehicle can choose not to fly the optimal path, but this will only result in higher probability of detection, which is worse for the vehicle. The sensors can also choose not to place at optimal positions, but this will only result in lower probability of detection, which is worse for the sensors. So this is like a Nash equilibrium in game theory. 5.1. Single sensor positioning The solution of this problem without constraints is always extreme: if the detection rate is proportional to the inverse power of distance and there is only one detector, then the optimal position of the detector is either at the source or at the destination. That will ensure a probability of detection equal to 1. In a more realistic scenario, suppose an aircraft plans to fly across the river. It sets off from (−100, 0) and finally reaches (100, 0). As it has a path planning system and knows the position of the detector, an optimal path can be chosen to minimize the probability of being detected. A detector whose detecting rate is related to the inverse square of the distance from itself to the target is located on a ship; thus the detector can only be placed in the river. The goal is to find the optimal position of the detector. In Fig. 9, the circle and the cross are the source and the destination of the aircraft respectively. The dotted area represents the river. By doing the simulation, one can obtain the optimal detector positions (−60, 0) and (60, 0). The contour lines of the probability of being detected (under path planning for minimizing detection) can also

136

B. Jiang et al. / Automatica 60 (2015) 127–139

By discretizing the possible positions of sensors, the original problem can be converted to a sensor selection problem. Suppose there are n identical sensors and m possible positions to place those sensors. Note m > n and there can only be at most one sensor located at one possible position. Remark 3. We first consider the simplest case when all sensors are identical. Later we will consider sensor selection with nonidentical sensors. Suppose ξ = [ξ1 ξ2 · · · ξm ]T is a column vector and ξi represents the presence of a detector at the ith position. The optimization problem is maximize

f (ξ )

subject to

ξi ∈ {0, 1} 1T ξ = n

ξ

Fig. 9. Illustrative example.

where f (ξ ) = min

m

x(t ),y(t )

Fig. 10. Floating starting and ending points.

be obtained. This simulation result also justifies our statement that the solution of this problem without constraints is always extreme. The simulation result suggests that the optimal position of the detector is (1) collinear with the source and destination and (2) as close to either source or destination as possible. There is also a possibility that the start and end points of the aircraft can be optimized by the path planning system on the aircraft as well. As shown in Fig. 10, suppose the starting point of the aircraft can be anywhere on the line {(x, y)|x = −100, y ∈ [−100, 100]} and the end point can be anywhere on the line {(x, y)|x = 100, y ∈ [−100, 100]}. The aircraft can choose its optimal starting and end points by a path planning algorithm. Then the optimal position of the detector is (0, 0). In this simulation, the optimal position of the detector is right in the centre of the admissible area. 5.2. Positioning of multiple identical sensors via convex optimization When multiple sensors are involved, exhaustive search is no longer practical. In this paper, a convex optimization method is derived. Remark 2. The original sensor positioning problem is not a convex optimization problem, as demonstrated in Section 5.1, in which example there are more than one local optimum.

(38)

Q i =1

ξi · pi (x(t ), y(t ), Dxi , Dyi )dt

(39)

and Dxi and Dyi are the coordinates of the ith possible position of detectors. Notice that the detection rate function is not necessarily bounded in the above function, though the meaning of an unbounded function multiplied by a zero ξi is ambiguous. Let p ξ ( x, y ) = m i=1 ξi · pi (x(t ), y(t ), Dxi , Dyi ). Typically, when we use the detection rate function defined by the signal to noise ratio, the value of detection rate function will approach infinity if and only if the distance between the detector and vehicle approaches zero. Thus the isolated singularity should not cause a problem. We relax this problem so that ξi ∈ [0, 1]. We define the detection rate function when ξi ∈ [0, 1] to be ξi pi (x, y, Dxi , Dyi ). Note the quantity ξi does not have any physical meaning because a radar does not have state between being switched off and switched on in reality. Lemma 1. Suppose the function f (ξ ) is defined as in (39). ∀ β ∈ [0, 1], there holds f (βξ ) = β f (ξ ). Proof. f (βξ ) = min

x(t ),y(t )

m

Q i =1

β · ξi · pi (x(t ), y(t ), Dxi , Dyi )dt

Since scaling pi will not cause any change in Q , which is indicated by the system of Eqs. (24) and (25), one can write β outside of the minimization operator. Thus f (βξ ) = β f (ξ ). Lemma 2. Suppose the function f (ξ ) is defined by (39). Suppose further that a and b are two m dimensional vectors. There holds f (a) + f (b) < f (a + b). Proof. Let Qa represent the optimal path where sensors are selected according to a; let Qb represents the optimal path where sensors are selected according to b; let Qa+b represents the optimal path when sensors are selected according to a + b. It should be noted that Qa , Qb and Qa+b have the same start and end points. In addition, pa (x, y) is the overall detection rate function when sensors are selected according to a; pb (x, y) is the overall detection rate function when sensors are selected according to b. Thus we obtain pa (x, y)dt ≤ pa (x, y)dt . (40) Qa

Qa+b

B. Jiang et al. / Automatica 60 (2015) 127–139

This is because Qa is the optimal path of pa (x, y) while Qa+b is not. Similarly, we have Qb

pb (x, y)dt ≤

Qa+b

pb (x, y)dt .

(41)

As a result, f (a) + f (b) =

Qa

pa (x, y)dt +

≤ Qa+b

pb (x, y)dt

Qb

pa (x, y)dt +

Qa+b

pb (x, y)dt

= f (a + b).

(42)

Therefore, f (a) + f (b) ≤ f (a + b). Theorem 3. The function f (ξ ) is a concave function.

137

5.4. Discussion and interpretation of the results Now that we have shown that the relaxed problem is convex. It follows that there exist mature algorithms to solve such a problem. However, there is still work to do to obtain a solution to the original problem based on a solution to the relaxed problem. Suppose ξr is a solution of the relaxed problem, ξu is a solution to the unrelaxed problem and ξe is a estimated solution to the unrelaxed problem obtained from ξr . A simple strategy to obtain ξe is thresholding rr ; the threshold should be chosen such that all the numbers in ξe add up to nj . Furthermore, one can swap pairs in ξe to obtain more accurate local optimum; details of such sophisticated local optimization methods are given in Joshi and Boyd (2009). For a given solution ξe , one needs an accuracy measure to evaluate how well it approximates ξu . Although the real value of ξu is unknown to us, there holds f (ξe ) < f (ξu ) < f (ξr ), which sets an upper and lower bound on ξu . Thus a reasonable measure for accuracy could be

Proof. Suppose a and b are two vectors and β ∈ [0, 1]. According to Lemmas 1 and 2, there holds

κ=

∀β ∈ [0, 1], a, b

(43)

where a smaller κ implies higher accuracy.

f (β a + (1 − β)b) ≥ f (β a) + f ((1 − β)b) = β f (a) + (1 − β)f (b).

(44)

5.5. A fast modification following addition of sensors

Because we intend to maximize a concave function, the problem is indeed a convex optimization problem.

5.3. Sensor selection for different types of sensors via convex optimization In the prior section, the sensor selection algorithm assumed identical sensors. When there are different types of sensors to be deployed, a modification to the problem statement is required. Suppose there are k types of detectors to be deployed and the number of type j detector is nj ; note j = 1, 2, . . . , k. There are m possible positions to locate those detectors. Suppose further that

ξ11

ξ21 ξ = .. .

ξm1

ξ12 ξ22 .. .

ξm2

··· ··· .. . ···

ξ1k ξ2k .. .

f (ξ )

subject to

ξij ∈ {0, 1} 1T ξ = n

f (ξ ) = min

x(t ),y(t ) Q j=1 i=1

n Proof. If f n2 ξr (n1 , m) > f ξr (n2 , m) , then ξr (n2 , m) is not 1 a solution to the problem of selecting n2 from m, which is contradictory to the definition of ξr (n2 , m). On the other hand, if

ξ (n1 , m) < f ξr (n2 , m) , then according to Lemma 1, we

n2 n1 r

2

(45)

where n = [n1 n2 · · · nk ] and k m

Lemma 4. Suppose 0 < n1 < n2 ≤ m are three integers. If all n elements in the vector n2 ξr (n1 , m) are less than or equal to 1, then 1 Ç å n2 f ξr (n1 , m) = f ξr (n2 , m) . n1

n have f ξr (n1 , m) < f n1 ξr (n2 , m) . This equation means that

is a matrix and ξij represents the presence of a type j detector at the ith position. The optimization problem is ξ

(47)

f (ξr ) + f (ξe )

Define ξr (n, m) as a solution to the relaxed optimization problem for selecting n sensors from m possible positions. Suppose one has already obtained ξr (n1 , m); due to the addition of sensors, one wants to find ξr (n2 , m), where n2 > n1 . The problem is whether one can find a fast modification method and does not need to do the whole calculation again.

f

ξmk

maximize

f (ξr ) − f (ξe )

ξij · pi (x(t ), y(t ), Dxi , Dyi )dt

ξr (n1 , m) is not a solution to the problem of selecting n1 from m, which is contradictory to the definition of ξr (n1 , m).

Lemma 4 tells us that if n2 ξr (n1 , m) lies in the admissible area, 1 then it is a solution to the problem of selecting n2 sensors from m positions. On the other hand, if it is not in the admissible area, it gives a good initial value for a gradient-based convex optimization approach. As shown in simulation results, when n1 < n2 ≪ m, it n appears very likely that n2 ξr (n1 , m) lies in the admissible area. 1 In practice, we can run a convex optimization algorithm to obtain ξ r (1, m), and use this to compute solutions when n

n < floor 1/ max ξr (1, m) , where max ξr (1, m) denotes the

(46)

and Dxi and Dyi are the coordinates of the ith possible position of detectors. It is straightforward to see that after relaxation, the concavity of f (ξ ) still holds. Therefore, the problem remains a convex optimization problem. The dimension of the relaxed sensor selection problem is m · k, which is independent of the number of sensors.

maximum value of elements in ξr (1, m). 5.6. Simulation Here is a scenario. On a 2D plane, a vehicle intends to travel from (−5, 0) to (5, 0). There are 9 possible positions to put sensors, as shown in Fig. 11. One can select n positions from 9 and put sensors on them. The objective is to maximize the probability of detection given the fact that the vehicle has prior knowledge about

138

B. Jiang et al. / Automatica 60 (2015) 127–139

Fig. 11. Selecting 2 sensors from 9.

Fig. 13. Selecting 10 sensors from 81.

than cumulative radar exposure. Furthermore, the variational dynamic programming approach is utilized to obtain accurate local optima of a dedicated original path. In addition, the homotopy method is proposed to adjust the optimal path when the detection rate function changes. Extending (Jiang et al., 2014), here we add a detection rate function for radars with Doppler capability. In addition, we also study the use of convex optimization for selection of detectors when a target vehicle can plan its path to avoid detection. In the future, a further level of optimization could involve the introduction of jamming devices to the model. Such jamming devices would be exterior to the region where the sensing devices were located, and might result in a new optimal path of the vehicle and a relocation of the optimally positioned sensors. There are likely to be optimal positions also for the jamming devices. References Fig. 12. Selecting 4 sensors from 9.

the positions of these sensors and can optimize its path to avoid detection. The solution to the relaxed optimization problem when n = 2 and 4 are shown in Fig. 11 and Fig. 12 According to Fig. 11, (−3.5, 0) and (3.5, 0) are the best positions when choosing 2 sensors from 9 positions. Because choosing 2 sensors from 9 only has 36 possibilities, we can also do an exhaustive search on this, which gives exactly the same result. It is shown in the simulation that ξr (4, 9) = 2ξr (2, 9). This is also consistent with Lemma 4. Suppose {i, j} denotes the ith sensor in the jth row in Fig. 12. It is noticeable that the fourth and fifth largest number in ξr are the same in Fig. 12. Thus we can choose either {1, 2}, {2, 1}, {2, 2} and {2, 3} or {2, 1}, {2, 2}, {2, 3} and {3, 2}. These two choice are identical in theory because the sensor are located symmetrically. The optimal result obtained via exhaustive search is {1, 2}, {2, 1}, {2, 3} and {3, 2}. Solutions via exhaustive search are not typically feasibly obtained in larger scale problems. Nevertheless, methods such as swapping and branching may be helpful to find near optimal solutions. A larger scale example is given in the case that one intends to select 10 sensor from 81 possible positions. The result is shown in Fig. 13. Contour lines show the trend of the relaxed result and red triangles represent the chosen positions. 6. Conclusion In this paper, we proposed an optimization model where the objective is to minimize the probability of being detected rather

Barraquand, H., & Ferbach, P. (1993). Path planning through variational dynamic programming. Barrett, S. R. (2007). Optimizing sensor placement for intruder detection with genetic algorithms. In Intelligence and security informatics, 2007 IEEE (pp. 185–188). IEEE. Bishop, A. N., Fidan, B., Anderson, B. D. O., Dogancay, K., & Pathirana, P. N. (2010). Optimality analysis of sensor-target localization geometries. Automatica, 46(3), 479–492. Dhillon, S. S., & Chakrabarty, K. (2003). Sensor placement for effective coverage and surveillance in distributed sensor networks. In Wireless communications and networking, 2003 IEEE: Vol. 3 (pp. 1609–1614). IEEE. Fox, C. (1950). An introduction to the calculus of variations. Courier Dover Publications. Franceschetti, M., & Meester, R. (2007). Random networks for communication: from statistical physics to information systems: Vol. 24. Cambridge University Press. Goldstein, H. (1980). Classical mechanics. In Reading Massachusetts: Addison-Wesley series in physics (2nd ed.). Addison-Wesley Publishing Company. He, Ji-Huan (2003). Homotopy perturbation method: a new nonlinear analytical technique. Applied Mathematics and Computation, 135(1), 73–79. Jiang, B., Bishop, A. N., Anderson, B. D. O., & Drake, S. P. (2014). Path planning for minimizing detection. In IFAC World Congress. Joshi, S., & Boyd, S. (2009). Sensor selection via convex optimization. IEEE Transactions on Signal Processing, 57(2), 451–462. Keller, H. B. (1992). Numerical methods for two-point boundary-value problems. New York: Dover publications. Kim, J., & Hespanha, J. P. (2003). Discrete approximations to continuous shortestpath: Application to minimum-risk path planning for groups of uavs. In Decision and control, 2003. Proceedings. 42nd IEEE conference on: Vol. 2 (pp. 1734–1740). IEEE. Liu, S., & Guo, H. (2002). Electromagnetic interference shielding and waveabsorbing materials. Journal of Functional Materials and Devices, 8(3), 213–217. Mahafza, B. R. (2013). Radar systems analysis and design using MATLAB (3rd ed.). Chapman and Hall/CRC Press. Maithripala, D. H. A., Maithripala, D. H. S., & Jayasuriya, S. (2008). Unifying geometric approach to real-time formation control. In American control conference, 2008 (pp. 789–794). IEEE. McMillan, R. W., & Kohlberg, H. (2010). A probabilistic model of the radar signal-toclutter and noise ratio for weibull-distributed clutter. In Radar conference, 2010 (pp. 882–886). IEEE. Murphey, R., Uryasev, S., & Zabarankin, M. (2004). Optimal path planning in a threat environment. Recent Developments in Cooperative Control and Optimization, 349–406.

B. Jiang et al. / Automatica 60 (2015) 127–139 Royset, J. O., Carlyle, W. M., & Wood, R. K. (2009). Routing military aircraft with a constrained shortest-path algorithm. Military Operations Research, 14(3), 31–52. Simon, M. K. (1998). A new twist on the marcum q-function and its application. IEEE Communications Letters, 2(2), 39–41. Skolnik, M. I. (2008). Radar handbook (third ed.). Zabarankin, M., Uryasev, S., & Murphey, R. (2006). Aircraft routing under the risk of detection. Naval Research Logistics, 53(8), 728–747.

Bomin Jiang was born in Beijing, China. He received the B.E. degree with First Class Honours in Mechatronic Systems Engineering from the Australian National University and the B.E. degree in Industrial Engineering from Beijing Institute of Technology. Currently he is working as a research visitor in National ICT Australia. From September 2015, he will start his graduate study in Control, Instrumentation & Robotics at the Department of Mechanical Engineering, Massachusetts Institute of Technology. Bomin Jiang has been awarded with a University Medal upon completing his B.E. degree in July 2014. His current research interests include distributed control and optimization theory. Adrian N. Bishop has held academic positions at the Royal Institute of Technology (KTH) in Stockholm, at the Australian National University (ANU) in Canberra and at the University of Technology Sydney (UTS) in Sydney. He is also a Research Scientist at NICTA in Canberra and Sydney. He is funded by an ARC Discovery Early Career Research Award (DECRA) Fellowship, NICTA and the US Air Force among other funding bodies. His current research interests fall within the intersection of statistical control and estimation, statistical machine learning and distributed (or large-scale) applicability of such topics.

139

Brian D.O. Anderson was born in Sydney, Australia, and educated at Sydney University in Mathematics and Electrical Engineering, with Ph.D. in Electrical Engineering from Stanford University in 1966. He is a Distinguished Professor at the Australian National University and Distinguished Researcher in National ICT Australia (NICTA). His awards include the IFAC Quazza Medal in 1999, the IEEE Control Systems Award of 1997, the 2001 IEEE James H Mulligan, Jr Education Medal, and the Bode Prize of the IEEE Control System Society in 1992, and a number of other medals and best paper prizes. He is a Fellow of the Australian Academy of Science, the Australian Academy of Technological Sciences and Engineering, the Royal Society, Honorary Fellow of the Institution of Engineers, Australia, and a Foreign Associate of the US National Academy of Engineering. He holds honorary doctorates from a number of universities, including Université Catholique de Louvain, Belgium, and ETH, Zürich. He is a past president of the International Federation of Automatic Control and the Australian Academy of Science. He served as the first President of NICTA, and was a member of company boards, including Cochlear Ltd., the world’s major supplier of bionic ears, and a member of the Prime Minister’s Science Council under three prime ministers. His current research interests are in distributed control, sensor networks and econometric modelling.

Samuel P. Drake is a senior research scientist with the Defence Science and Technology Organisation (DSTO) Australia. He graduated with first class honours in physics from the University of Melbourne and completed his Ph.D. in general relativity at the University of Adelaide. He is on the editorial advisory board for the American Journal of Physics and lectures ‘‘Applications of Relativity’’ at the University of Adelaide.