PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 5/08 K:IProj/Acc&Pres/ClRCIDateDuo‘indd OPTIMAL SAMPLING STRATEGIES USING SPATIALLY AVERAGED ADVECTION-DIFFUSION PARAMETER ESTIMATION By Andrew Philip White THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Mechanical Engineering 2008 ABSTRACT OPTIMAL SAMPLING STRATEGIES USING SPATIALLY AVERAGED ADVECTION-DIFFUSION PARAMETER ESTIMATION By Andrew Philip White In this paper, optimal sensing strategies are presented which use spatial averaging of noisy measurements of an advection-diffusion process for parameter estimation to trace the maximum of a field of interest. Each strategy utilizes nonlinear least- square (NLS) optimization to estimate the parameters of a locally modeled advection- diffusion process. In the first strategy, nonholonomic sensing agents are steered in the direction that maximizes an objective function. Using ideas from the first strategy, a second strategy is designed to find the maximum of a complicated advection—diffusion process in a non-convex surveillance region. Each time the NLS optimization is run, the target position for the nonholonomic sensing agents is updated with the estimated maximum concentration position of the locally modeled advection-diffusion process. To ensure the nonholonomic sensing agents do not collide with boundaries in the non-convex region, an object avoidance algorithm is utilized. ACKNOWLEDGMENTS I would like to thank my advisors Dr. Jongeun Choi and Dr. L. Guy Raguin, for their guidance on this project. I would also like to thank Dr. Subir Biswas and Dr. Seungik Baek for being members on my committee. iii TABLE OF CONTENTS LIST OF TABLES ............................... v LIST OF FIGURES .............................. vi 1 Introduction ................................. 1 2 Problem Statement ............................. 4 3 Averaged Process .............................. 6 4 Vehicle Formation Control for a Control Volume .......... 9 5 Feedback Linearization of Nonholonomic Agents .......... 13 5.1 Single-Integrator Compenstated System ................. 14 6 Control Strategies .............................. 16 6.1 Control strategy for Problem 1 ..................... 16 6.2 Control Strategy for Problem 2 ..................... 18 7 Simulation Results ............................. 24 7.1 Simulation Results from Problem 1 ................... 24 7.2 Simulation Results from Problem 2 ................... 26 8 Future Work ................................. 32 8.1 Numerical closed-form solution of the Local Process .......... 32 8.2 Function Approximation ......................... 34 9 Conclusion .................................. 37 APPENDICES ................................. 38 BIBLIOGRAPHY ............................... 41 iv LIST OF TABLES 7.1 Parameters in the simulation ........ 7.2 Estimated Parameters from each simulation oooooooooooooo 4.1 5.1 6.1 6.2 6.3 6.4 6.5 6.6 7.1 7.2 7.3 7.4 7.5 7.6 8.1 LIST OF FIGURES Images in this thesis are presented in color. Possible distribution of sensing agents in a control volume. ...... 10 Schematic of the ith Nonholonomic Agent ............... 14 Flow Chart of the Optimal Sampling Strategy ............. 17 Sampling Strategy Timeline ....................... 19 Contaminant Source with complex flow field .............. 19 Sample grid of a non-convex region ................... 22 Path found from (2,1) to (1,4) by A* search algorithm ........ 22 Search graph to get from node (2,1) to node (1,4) ............ 23 Nonholonomic Agents starting at x = 20 and y = 20 ......... 26 Nonholonomic Agents starting at x = 20 and y = 120 ........ 27 Nonholonomic Agents starting at x = 120 and y = 20 ........ 27 Comparison of the true plume and estimated plume after 150 iterations. 29 Comparison of performance with and without “anemotaxis”. ..... 30 Sensing Agents attemping to trace a more complex concentration field. 31 Local and Global Coordinate Relationship ............... 36 vi CHAPTER 1 Introduction Many physical transport phenomena exist that are governed by advection-difl'usion processes. Due to the harmful nature of some of these processes, researchers have been developing methods to track the source of such advection-diffusion processes [10, 5, 29, 14, 13, 8, 23, 28, 3]. The most common approach to this problem is bi- ologically inspired “chemotaxis” [1, 5], in which a mobile sensing vehicle is driven according to a local gradient of the chemical plume concentration. However, with this approach, the convergence rate can be slow and the mobile robot may get stuck in the local maxima of chemical plume concentration. If strong ventilation is present, an appropriate technique is “anemotaxis” in which a robot moves in the direction against the wind for locating the odor source. For example, a robot equipped with gas sensors and anemometric sensors moves in the direction of upwind [14]. A switch- ing mechanism between anemotaxis and chemotaxis according to the magnitude of measured odor concentration was designed in [13]. Before a sensing robot detects the plume concentration, popular biologically inspired exploratory movements of a robot include a biased random walk, a “run and tumble” flagellated bacteria like behavior, and “zigzag” and “spiral” trajectories. A methodology of identifying the source term for an instantaneous pollutant point source from measurements in different locations was presented in [15]. The parame- ters (such as chemical mass release, diffusivity, source height, location of the source and the time of the pollutant release) for the solution of an advection-diffusion model were estimated through nonlinear least-squares regression. This was possible since an analytical solution was obtained from a simple diffusion model. The advantage of parametric estimation is that it allows sensing robots to be steered in order to min- imize the uncertainty in the parameter estimation of a diffusion model. The control law of sensing robots for source localization was derived in [23, 28] by minimizing the error of the approximated maximum likelihood parameter estimation. The work [3] also described mobile sensing robots that estimate the parameters of the advection- diffusion equation. A major drawback of this approach is that the analytical solution of a chemical plume process is not available for most practical cases. For instance a gas release in an airport with a complex wind field (from ventilation and turbulence effects), it is impossible to pre-decide an analytical solution and its parameters that intelligent robots can learn about. With a wind field, dispersions of chemical plumes are more significant than diffusions. Given all boundary conditions and wind fields, the evolution of a plume can be simulated by the Large-Eddy Simulation (LES) code such as HIGRAD [4, 26] at the cost of high computational power. In addition, re— cently emerging distributed control strategies over mobile sensing networks have not been considered in the control algorithms for chemical plume tracing developed in [23, 28, 3]. In recent years, the cooperative control of mobile sensor networks has gained much attention from control engineers and scientists. The cooperative network of vehicles that performs adaptive gradient climbing in a distributed environment was presented in [11, 17]. The network can adapt its configuration in response to the sensed environment in order to optimize its gradient climb. Due to the complexity of realistc advection-diffusion processes that are seen in practice, it is clear that obtaining an analytical solution is not only unpractical but in many cases is also not possible. Realistic advection-diffusion processes can be nons- mooth, turbulent, and complex. However, if spatial averaging is used, measurements taken inside a control volume can be used to estimate the parameters of a locally averaged model. The locally averaged model can be smooth and simple with well defined parameters. These parameters can then be utilized to drive the robots to trace the field of interest that is undergoing a complex transport phenomena. Hense, the objective of this research is to synthesize a mobile sensor network that is capable of finding the source of contaminant. This paper is organized as follows. Chapter 2 describes the problems that are addressed in this research and a brief desciption of the methods used to solve them. Chapter 3 introduces the concept of an averaged process. Chapter 4 discusses the formation control censensus algorithms used to ensure that a consistent control vol- ume is maintained. Chapter 5 shows how the nonholonomic contraints of the sensing agents are handled. Chapter 6 details the control strategies used to solve problems 1 and 2. Chapter 7 shows the results that were obtained from each control strat- egy. Chapter 8 introduces additional concepts needed to proceed further with this research. Finally, in Chapter 9, conclusions are stated. CHAPTER 2 Problem Statement Problem 1: Consider an idealized advection-diffusion process in a convex surveil- lance region. Using a network of nonholonomic sensing agents, estimate the param— eters of the idealized advection-diffusion process. Moreover, coordinate the mobile sensor network to improve the quality of the estimates. Problem 2: Consider a realistic advection-diffusion process in a complicated non- convex surveillance region. Using a network of nonholonomic sensing agents, locate the maximum of the advection-diffusion process. The first problem is solved by considering a control volume over which spatial av- eraging takes place. To maintain a moving control volume, the robots are controlled such that they travel in a formation [27]. The mobile sensing agents take concentra- tion measurements inside the control volume. Nonlinear least-squares regression is used to estimate the parameters of the locally averaged process model. To determine the uncertainty in the estimated parameters, the Fisher Information Matrix is used to find the Cramer-Rao Lower Bound [16]. Then to improve the quality of the esti- mated parameters, the determinant of the Fisher Information Matrix is maximized (D-Optimality) [7]. By taking the partial derivative of the determinant of the Fisher Information Matrix with respect to x, control inputs for the formation control algo- rithm are generated which drive the formation in the direction that maximizes the determinant of the Fisher Information Matrix. The second problem is solved by using the same spatial averaging and control volume approach that the first problem utilizes. However, to solve the second problem, instead of using D-Optimality to drive the sensing agents, the initial position of the estimated diffusion impulse response is set as a target position. Due to the non-convex geometry under consideration, an obstacle avoidance algorithm is used to track a safe, collision-free path from the start position to the goal position. The A* graph search [20] is used to plan the collision-free path. CHAPTER 3 Averaged Process In this section, the averaged process that is used to model the local, spactially- averaged, advection-diffusion process is introduced. The transport of a solute of concentration C (x, t) in a medium is governed by the advection-diffusion equation BC where q is the flux of the solute and is decomposed into an advective and a diffusive term 11. is the velocity vector, and D3 the molecular diffusion coefficient of the solute in the medium. Assuming incompressible flow (V - u), with spatially averaged velocity field V, the governing equation can be rewritten in its volume-averaged form [22] BC“ + at where D is the dispersion tensor that contains both advective and diffusion terms, as v - v0 = DV2C', (3.3) well as effects of the geometry, and may be time-dependent [2]. D is symmetric and usually contains cross-terms. In three-dimensional space, V = (Vx, Vy, Vz)T and Data: Dry sz sz Dyz Dzz therefore the governing equation for C(x, y, z, t) can be written as 6C +V$ (Ci—0+ +Vy 6870+ Vz— 6C 075 82 8:0 0:0 a__20 6:6 a__:c: a_2_C*z (3.5) In an infinite medium, the response to an impulse function, C(x,y,z,t0) = CO 6(m0,y0, z0,t0), located at an unknown position X0 = ($0,310, 20) and unknown time to is then C0 \/(47r)3 (t — t0)3detD x e, [_(x—xo-v) ‘p 4(t—t0) CON) = (3.6) At this point, the scope of consideration is reduced to the two-dimensional process in the (x,y)-plane. The simplified response to the impulse function, C(x,y,t0) = Co 5($0,y0,t0), is -ny Dim: _ Co C(x, y, t) = \/<4vr>3 (t — t0)3(D:r:chy — 0.233,) 1 x ex — 3.7 p{ 4(Dxnyy — Dgy) (t " t0) ( ) T Dyy —ny X [(X — X0 — V(t - t0))l [ ] [(X - X0 — V(t — t0))l}- The goal is then to recover the following parameters: D333, Dyy, ny, CO, :30, yo, and t0, based on measurements of C (x, y, t) at different locations and times. The sensing agents are assumed to be able to measure wind velocity, such that V5,; and Vy are considered to be known. The concentration function is rewritten as C(96), (3.8) where 6 is the parameter defined by T 7 and the configuration Q is Q:=[x y t]TER3. (3.10) CHAPTER 4 Vehicle Formation Control for a Control Volume To describe the information flow between each of the sensing agents some graph notation is now introduced. A directed graph consists of a pair (N, 8) where N is a set of nodes and 8 is a set of ordered node pairs, also known as an edge (see [24]). For the formation control [27] used in this research, a graph known as a directed spanning tree is required. A directed spanning tree is formed by graph edges which connect all of the nodes of the graph. Examples of a directed spanning tree are illustrated in Fig. 4.1. To maintain a control volume for the averaged process estimation, the sensing agents are arranged into a controlled formation. The reference point of the control volume is refered to as the virtual center (xvcwvc). The position of all sensing agents is determined with respect to the virtual center based on the shape, size, and the number of agents in the control volume. Each agent has its own understanding of the location of the virtual center, which is defined as qz- = [1:30 yfc 6¥C]T for the ith agent. The path the formation follows is defined as qd = [avgC ygc 636]? According to [27], if the adjacency matrix GUC = v.0 [97; J ] has a directed spanning tree, then q,- —> qd asymptotically using the consensus (a) (b) . <9 Figure 4.1. Possible distribution of sensing agents in a control volume. algorithm 52(qd— «(xi—«1.1» +912 _19,-J-C[qJ- 7(qJ-qJ-H 4.: , (4.1 2 I8,- +Zj—1gijc ) where 6,- = 1 if the ith agent is a leader and zero otherwise, 9%)]? = 1 if agent j’s estimate is available to agent 2' and zero otherwise, and ’y > 0. To demonstrate that q,- —> qd asymptotically, the error dynamics are now investigated. First the denominator of equation (4.1) is multiplied over to the left hand side of the expression (‘3 + gvclqz' =3 [dd “((4 - 4(1)] l+ Z 9,- ”J-CldJ-V (qr- qJ)l (4-2) where n _ vc — Z gij j =1 Next, gchd is subtracted from each side of the expression (9 + gvclqz' gchd — €le 701-qu + j: g, ”JCIQJ-v - qu)l - $640). (43) After some more re—arranging, the expression becomes (‘3 +9116) (42"le —€7(qz -qd)+ZgJJ-cl(qJ'- qd)- 7am -qd)-(qJ--qd))l- (4.4) 10 Substituting (j,- = q,- — qd into the expression, the error dynamics are (8 + 9%” =4qu- + E g,- ”J-Cqu-v — (25)]. (4.5) which can be arranged into matrix form q‘: Bq" + Aq (4.6) where B and A are found to be '06 DC DC ‘ git 9% 9%? B: 1 921 922 92n (6+QUC) 5 3 ’ 9% 91%- A_ 1 7921 —7(€+Q ) 791,, — (6+Q'UC) ~ 3 .. 7933f -7(€+ch)- After some matrix operations, it can be shown that é=Fq with F = (I — B)—1A. For 7 > 0 [27], the eigenvalues of F have negative real parts, which shows that (j —> 0, implying that 92' —> qd asymptotically. It is assumed that the sensing agents have integrator dynamics such that ’izu'i) i=1,..-,n (4.7) where r,- = [1:23 31,-]T is the position and uzz ["32” um]T is the control input to the 3th vehicle. The agents are controlled by the following extended consensus algorithm as shown in [27]: km; uJ=r~§1—a,-( (“rz —'r§l )—Zgz-j ICU-[(rz r-—r )] (4.8) 11 where a,- > 0, kij > O, gij = 1 if information flows from vehicle j to vehicle 2' and 0 otherwise, and r? = [23%, yzd]T with rczd J _ [: 3:36 J + [ 003(636) —sin(6§’c) yg yyc sinwa) cos(0§’c) xiF yiF th vehicle’s position in the formation where [:132 F7 y,- F]T is the vector that gives the z' with respect to the virtual center. The extended concensus algorithm can be arranged into the following matrix form. a = Jf (4.9) where 1",- := r,- — 71?, 11:: u — 1'"? and P -01 + 231:1 gijkij 9121612 ginkln ' J _ 9211621 ~02 + 2321 Qijkij ... 9272.an 917.1an ~- ... ~an + $321 gzjkz'j . It is important to notice that J is diagonally dominant. From the Gershgorin disc theorem [12], it is easy to see that all of the eigenvalues of J have negative real parts. This guarantees that f ——> 0 exponentially, which shows that 7“,; —> r?, W. 12 CHAPTER 5 Feedback Linearization of Nonholonomic Agents The sensing agents to be used in the implementation of this algorithm have nonholo- nomic dynamics. Dynamics that are said to be “nonholonomic” have constaints which are not integrable. This essentially means that knowing an ending configuration \112 does not give information about the starting configuration \I/l. The nonholonomic dynamics simply constrain how the sensing agent can move from any starting con- figuration ‘Ill to any ending configuration \112. The kinematic equations that govern the dynamics of the 3th nonholonomic agent are 7252. = ’0,- cos 62-, f‘yi = ’07; sin 92', 9i = w, (5.1) h nonholonomic where (rxz-in) and 0,- give the position and orientation of the it agent and the inputs v,- and w, are, respectively, the forward and angular velocity of the vehicle (Fig. 5.1). There are many methods to feedback linearize equation (5.1) [21]. The controller discussed in this section will be used for trajectory tracking of a reference sensor position. Posture stabilization, a more difficult problem used in parking maneuvers, is not considered in this paper; since sensor position control is 13 YA vi (1172', yi) 6'; \‘7 “’2' ........ 7. ....... I.’ d2; (Til/'2’ ) Tyi) XV Figure 5.1. Schematic of the ith Nonholonomic Agent sufficient for the purposes of this research. Posture control can be found in [21]. 5.1 Single-Integrator Compenstated System To acheive single-integrator dynamics, the coordinate transformation of 11:,- : r557; + (1,; cos 92-, y,- = ryz- + dz- sin 02-, (5.2) is applied where (112,, 92) denotes the position of the ith sensor and (12- represents the distance from the wheel axis to the sensor (Fig. 5.1). This moves the point of interest from (Txi , r313) off the wheel axis by a distance of dz- to ($75,193)- Taking the derivative of equation (5.2) and entering the values for (mi, ryz.) from equation (5.1) produces 14 the following equations: 1152' = ’07; COS 02' — d7; sin(0i)wz-, 3],- : ’Ui sin 62' + di COS(67;)wz'. (5.3) The output of 3th nonholonomic agent is defined as 77,- : (xi,yz-)T. Taking the derivative of the output produces . :15- cosB- —d-sin6- v' w] .2] 2 [J ,2 -‘ll .[ (.4) 3,12 BID 2 (12 cos 6, to, which shows that both v,- and w,- are present in 77,-, meaning they can be used to h control the output of the it nonholonomic agent. Thus, the linearizing feedback control can be found by letting 1),- = (uxi,Uyz-)T such that use, = cos 0,- —d,- sin 6,- vz- . (5.5) “312' sm 6,- d,- cos 02- cu, By taking the inverse of the matrix above the inputs vi and w,- to the ith nonholonomic agent are found to be ”2' = cos 0,- sin 6,- “xi 5 6 _ 1 - 9. 1 0' ' ( - ) “’2' 3; sm 2 a; COS Z uyi By entering the values of v,- and w,- found in equation (5.6) into equation (5.3), it is lyzl=luyzl M It should be noted that a singularity occurs when (12' = 0. If the position of interest is easy to see that (mi, ryi), then double-integrator dynamics are necessary [25]. However, in this case, only control of the sensor position is necessary. Thus this will be sufficient since the sensor can be placed at an arbitrary distance 01,- from the wheel axis as shown in Fig. 5.1. 15 CHAPTER 6 Control Strategies 6.1 Control strategy for Problem 1 In this section, a navigation strategy that coordinates the mobile sensor network to improve the quality of the parameter estimates is provided. Let T; be the sensing h sampling time, the leader vehicle will collect the spatially sampling time. Every pt distributed measurements from the follower vehicles to update the parameter esti- mate 6. The estimate of the parameter is updated via the nonlinear least-squares optimization (for example, by the Levenberg-Marquardt method) Wk) = h($(tk), 905;, - T510». (6-1) It is then natural to coordinate the mobile vehicles a:(t) for t E (tk’tk + T 319] in order to maximize the determinant of 2 in equation (6.4) evaluated at the current parameter estimate and minimize the covariance of the estimated parameter 8. A flow chart of the optimal sampling strategy is given in Fig. 6.1. We introduce the Fisher Information Matrix (I), given by [16] . 2(6*) - h . n C, 0* C, 0* T 6 2 i=1 where C’(Q,,-,0) := [6%7C(Qz-,6)] . Since the true parameter 9* is unknown, the 3 j estimated parameter 6 is used to compute the Fisher Information Matrix in equation 16 Mollower Agenti ] Feedback Linearization (712': 6%) I Controller Ti— Formation Control Module 7‘ j KI-Q. 1.9.9.. ’06 qt Formation State Data (Cir (12') Nonholonomic _ Storage / Vehicle Estimator vc ea (02" (12') ‘1 j ‘12“ Y Communication Network Topology vc ’00 q 9- Data Optimal (qgc, 43¢) k ‘7 Storage Reference Formation State (G: (I) Generator Estimator .. vc Nonlinear Least 9(t + 1) qk rk _ rd Squares T‘k , k Formation _ d Control Module 7‘] — J “k (Ck, qk) Nonholonomic (12k, wk) Feed back Vehicle Linearization Controller [ LeaderAgentk 7 Figure 6.1. Flow Chart of the Optimal Sampling Strategy 17 (6.2) instead. The Cramer-Rao theorem states that the covariance matrix of any unbiased estimator d is lower bounded by the Cramer-Rao Lower Bound (CRLB), which is the inverse of the Fisher Information Matrix [7]. E {(é — 0*)(é — 6*)T} 2 I-1(a*) (6.3) In order to improve the quality of the estimated parameter 6, the determinant of the Fisher Information Matrix will be maximized (D-optimality [7]). 1 j = det :r = ——2—det2 (6.4) 0' The objective function J is maximized by steering the sensing vehicles so that they climb the gradient of the objective function with respect to x, which is given by a_.7 _ 2 6.7(2) 0_2 8X — 2,], 32 (9X 16 _ 32x) 2 §(det(2)(2) T) (73%)]c (6.5) The vehicle control is defined as 'v if llvll < vsat; -vc = . 6.6 qd { fl: [Usat 1f [[21]] Z ”sat; ( ) where v = 5%, ”salt is the saturation velocity of the vehicle and 5 > O is a gain. Also, to ensure that the sensing agents stay inside of the surveillance region Q, a projection algorithm is used, which updates the current position if it belongs to Q and keeps the previous position if it does not belong to Q [18]. 6.2 Control Strategy for Problem 2 The approach to solve this problem combines (A) the nonlinear least-squares es- timation approach, which was used to solve the simpler more ideal problem posed in the previous section, with (B) the wind tracking stategy known as “anemotaxis” 18 Path Planning NLS Path Planning (121—1 Sampling i ‘Id Figure 6.2. Sampling Strategy Timeline Figure 6.3. Contaminant Source with complex flow field 1+1 —"1d [14]. Anemotaxis, as stated in the introduction, is a technique used to send a robot in the direction against the wind for locating an odor source. A switching mechanism is used to determine whether the sensing agents will go to the goal position given by the estimated parameters or follow the wind velocity. The condition for choosing which strategr to use is based on the objective function (6.4). Only parameter estimates that produce a J greater than jmin are trusted. If a value of ,7 is computed to be lower than jmini then “anemotaxis” is implemented. Since the complicated non-convex region is considered (Fig. 6.3), an obstacle avoidance strategy is necessary to ensure that the sensing agents do not bump into the walls. In this work, it is assumed that the robots have knowledge of the nonconvex geometry. Each of the walls (as in Fig. 6.3) are placed inside a no-fiy zone of which the leader robot has full knowledge. The choice of the physical process that models the local behavior of the plume is important. Many options exist from which to choose. In the previous section, an advection-diffusion impulse response is used. For this problem, the following three options were considered: diffusion impulse response Eq. (3.6) (with V = O), advection- diffusion impulse response Eq. (3.6), and continuous source isotropic diffusion [19] I _ (loft —t0) . d C(X,t,x ,quo) — _2WKd erfc 2 K“ _ to) (67) where ‘10 is the continuous source rate, K is the isotropic diffusion coefficient, and d is the Euclidean distance from X to the source x’, d = ”X — x'llz. After each nonlinear least-square estimation is complete, an updated estimate of the target position is obtained. This new target position is then given to the obstacle avoidance algorithm, which creates a safe, collision free path from the starting postion to the target position. Before a path can be found, the surveillance region must be partitioned into a grid as illustrated by the sample region in Fig. 6.4. A node is placed at the center of each grid box. Also, the no—fly zone and boundary are indicated on the grid by giving any node in a grid box a terrain cost of 00, such that these nodes will not be considered during the path search. To find the path from the start node 20 to the target node, a search algorithm must be used to create a reasonable path. The search algorithm we used in this paper is the A* algorithm, which was introduced by Nilsson [20] and uses heuristic infomation to improve performance when compared to Dijkstra’s algorithm [6]. The A* graph search method creates a hierarchical tree search graph G (as displayed by the sample search in Fig. 6.6) to test many paths to find one with the lowest cost. The search graph G begines with the start node on the top of the tree. The possible nodes that can be travelled to from the start node are installed in G as children of the start node. Then, based on the heuristic cost [9] h: «56x min(|6g;|,]6y|)+ax ||6$| — lay”, (6.8) where a. is the axial distance between each path node, 6:1; is the distance between the current node and the target node in the a: coordinate, and 6y is the distance between the current node the goal node in the y coordinate, the next node to expand is chosen. This process is repeated until a path is found. The A* graph search is summerized below. 1. Create a search graph, G, which consists of only the start node, 3. Place 3 on a list called OPEN. 2. Create a list called CLOSED which initially has no nodes. 3. Start the search loop. If OPEN is empty, then exit with failure. 4. Select the first node on OPEN and name it 17.. Remove n from OPEN and place it on CLOSED. 5. If n is the goal node, then exit the loop with the path obtained by tracing the pointers from n to s. 6. Expand node n to find all possible children nodes which are not parent nodes of 71. Place these nodes in a set called M. 7. Create a pointer from the nodes in M, which are not already in OPEN or CLOSED, to 71. Add these members of M to OPEN. 21 8. Re-arrange OPEN according to the lowest heuristic cost. 9. Repeat the loop. @@@ 6D @@@ @@@ ®@O @ @ 6D Figure 6.4. Sample grid of a non- Figure 6.5. Path found from (2,1) to convex region (1,4) by A* search algorithm When this algorithm is carried out to find a path from a start node of (2,1) on the sample grid in Fig. 6.4 to a goal node of ( 1,4), the search graph is expanded as shown in Fig. 6.6 to find the path traced out in Fig. 6.5. Once a path, p = [110,111, . . . , pn], is found the reference trajectory for the forma- tion control, qfi(t), is computed by the leader vehicle as shown in Fig. 6.2. The reference trajectory is created by filtering a step function from p0 = [120,340] to p1 = [$1,111] through a low pass filter 11(3) = ”L“ (6.9) such that 43.10?) = (m — p0) (1 - exp 0%)) +100 (6-10) and (1210) = (1%) exr) (EL—TA). (6.11) 22 GD 6 0@ @® @® @° 3 bath Ends ] [Path Ends] [Path Ends] ® ® 0 We e 6 ea 96*) P e e @‘D e Figure 6.6. Search graph to get from node (2,1) to node (1,4). 23 CHAPTER 7 Simulation Results 7 .1 Simulation Results from Problem 1 The sampling strategy detailed Section 6.1 is tested using a Matlab simulation code. To simulate a physical process, a true set of parameters is defined. Table 7.1 displays the true parameters and the assumed apriori values used in the simulations. The true parameters are used to create the field of interest that the sensing agents take noisy measurements from. The concentration measurements are then used with the apriori parameter values to estimate the true parameters using nonlinear least-squares opti- mization. The apriori knowledge for the each of the parameters, with the exception of the initial position, was taken to be twenty percent away from the true value. The initial position of the plume was given a value far from the true value since the no apriori knowledge about the position of the plume is assumed except that it is inside of the surviellance region. In Table 7.2, the estimated values for each of the simulation trials are displayed. Figs. 7.1 - 7.3 display the sensing agents starting position in grey and their final position in black. The trajectory that each of the agents followed during the simulation is plotted in magenta and the reference trajectory is displayed in red. Figs. 7.1 - 7.3 also show that the sensing agents can be inserted into the field at any location without affecting their performance. However, Table 7.2 does show that 24 Table 7.1. Parameters in the simulation Parameters True Values apriori Values 2 Drum?) 0.07 0.084 2 D,,,,(-’%-) 0.07 0.084 2 D$J(IDS—) —0.01 —0.012 00(53):) 5.0x106 6.0x106 ($0, 90) (m) (80, 80) (1,1) to (s) 30.0 36.0 Table 7.2. Estimated Parameters from each simulation Estimated Values with starting position (31:3, ys) Parameters (20, 20) (120,20) 420, 120) 2 D21: (—mS—) 0.0699 0.0703 0.0702 2 pm, (my) 0.0697 0.0703 0.0702 2 my (—";—) -0.0099 —0.0101 —0.0101 00 ("52%) 5.0043 x 106 5.0149 x 106 5.0124 x 106 (80,110) (m) (7993,7989) (7993,8008) (8093,7995) t0 (.9) 24.2717 30.2394 29.4507 25 Figure 7.1. Nonholonomic Agents starting at x = 20 and y = 20 the parameters estimated by the robots does change slightly depending on the initial position of the agents. The estimated parameter which shows the highest amount of variation is the initial time, to. Fig. 7.4 shows the comparison of the plume computed using the true values and the plume computed using estimated values from the simulation starting at (4119.343) 2 (20, 20). 7.2 Simulation Results from Problem 2 The sampling strategy that was detailed in Section 6.2 was tested using a Matlab simulation code. The concentration field used in these simulations was created using Cornsol. The concentration field was created to model a constant contaminant release inside of a closed area. In Fig. 7.5 - 7.6, the nonholonomic agents are represented by a yellow plus sign at (1'3, ys) and by a white dot at their final location. Without utilizing the switching algorithm based on the objective function, the best results obtained are displayed in Fig. 7.5. The sensing agents clearly have not 26 Figure 7.2. Nonholonomic Agents starting at x = 20 and y = 120 Figure 7.3. Nonholonomic Agents starting at X : 120 and y = 20 27 found the maximum position in the field. This has happened because the physical process used to model the local plume does not match the local plume. To keep the sensing agents moving, the switching algorithm based on the objective function 6.4 was used. The simulation results obtained when using the switching algorithm are displayed in Fig. 7.5. In this simulation, the sensing agents do find the maximum of the field of interest. However, it is noticable that the parameter estimation algorithm does not get utilized until the sensing agents are almost at the maximum. It is also clear that with this field can be tracked with only the use of “anemotaxis”. To more completely test the strategy developed in Section 6.2, another concentra- tion field was created using Comsol. In this concentration field, a contaminate-free inlet was added below the contaminated inlet. The contaminate-free inlet is denoted in Fig. 7 .6 by the letter A, and the outlet is denoted by the letter B. The resulting concentration field cannot be traced with only the use of “anemotaxis”. However, as Fig. 7.6 illustrates, it appears as though the strategy developed cannot trace it either. This, again, is due to the mismatch between the modelled plume and the actual local plume. The physical process that most closely resembles the local plume behavior (of those considered) is the diffusion impulse response. The wind in the advection-diffusion impulse response moves the maximum concentration area of the plume away from the estimated initial point (320,310). Since the robots are moving towards the estimated initial point of the local plume, it is desireable that it does not move. Also, when investigating the continuous-source, isotropic-diffusion model, it was noticed that the source location has a very steep gradient that does not model the local behavior well at all. However, clearly, the diffusion impluse response does not model the local plume behavior well enough. For this reason, a numerical solution will be developed to better model the local plume behavior. The framework for creating such a numerical solution is laid out in the next section. 28 (a) True Plume 160 5 fl 2’ '2 160 (1)) Estimated Plume with starting position (1's, y...) : (20, 20) Figure 7.4. Comparison of the true plume and estimated plume after 150 iterations. 29 (a) Parameter Estimation Strategy Only 0 320 (1)) Switching between Parameter Estimation and “anemotaxis” f Figure 7.5. Comparison of performance with and without “anemotaxis’ . 30 Figure 7.6. Sensing Agents attemping to trace a more complex concentration field. 31 CHAPTER 8 Future Work 8.1 Numerical closed-form solution of the Local Process When considering different advection-diffusion processes, it is necessary to have a closed-form solution, like the one that provides the concentration field of the advection-diffusion process in equation (3.7). If an analytical closed-form solution for the local process of interest does not exist then a numerical closed-form solution must be found. The following steps can be done to find such a solution. Before the numerical analysis can be done, a dimensional analysis is carried out first so that non-dimensional parameters can be found. The steps used are outlined below. 1. First list all of the variables that are involved in the problem 0 = f(x. 9.61.90, D, V, V) (8-1) where C is the concentration, (23,31) is the position of the concentration C, d is the inlet size, ‘10 is the inlet rate, D is the isotropic diffusion coefficient, V is the inlet velocity, and V is the kinematic fluid viscosity. 32 2. Next the basic dimensions of each variable are expressed C 2 ML_3 :1: i L y i L d a L go a Mir—1 2 L2T_1 : LT—l V i L2T—1. 3. The number of pi terms is determined using the Buckingham pi theorem. The Buckingham pi theorem states that a relationship, with 19 variables, that is dimensionally homogeneous can be reduced to a relationship of 18—?" independent dimensionless variables, where r is the number of reference dimensions necessary to describe the variables. Since there are eight variables in (8.1) and three references dimensions are required to describe them, the number of pi terms is five k—r=8—3=5 (8.2) 4. Next, three repeating variables were selected to create the five non-dimensional pi terms. In this case, the repeated variables choosen are d, qO, and u. 5. The remaining variables are each multiplied by a combination of the repeating 33 variables to create the pi terms. 111 : _Cdy (8.3) ‘10 CE H2 = - (8.4) (1 H3 = % (8.5) 114 = V7“! (8.6) D 6. To make sure that each of the pi terms are dimensionless, the dimensions are double checked. (ML—3)(L>(L2T—1> r11 2 MOLOTO MT-1 112 _ éeMOLOTO L 113 .. ZéMOLOTO H w ; MOLOTO 4 L2T-1 ’— . L2T-1 . 0 0 0 7. Finally, the function is written as a relationship of the pi terms. Recognizing that inverse H5 is the Schmidt number, it’s inverse is used in the final relation- ship. CdV a: y Vd V E — f (J. J, 7, 5) (8.8) 8.2 Function Approximation To approximate the function given in (8.8), numerical simulations were completed for five different inlet velocities, five different inlet sizes, and five different isotropic diffusion coefficients. The data produced from the simulations was approximated with a combination of radial basis functions given by f9) := Z 4990,- : ¢T(A)6 (8.9) j=1 34 where ¢(A) and G are defined as 429) := [419) 4529) 4mm]. 9 := [9192 ---9ij The Gaussian radial basis functions {$00} are given by 2 1 H ’\’i - W H J i€{1,2,3,4} J where aj is the width of the Gaussian basis and F j is a normalizing constant. The centers of the basis functions {“ij I i E {1,2,3,4},j 6 {1,2, - -- ,m}} are uniformly distributed in the function space, where m = n1 x 712 x 723 x 714 and "2' is the number of centers needed to approximate A22 From the numerical simulations we have many observations of H1(’\k) = ¢T(/\k)0, where k is the measurement index. Based on the observations and refressors {(H1(/\k), ¢(/\k))}2=1, our objective is to find O which minimizes the least-square error aélH 15:1] |-r11(,\,,)-¢T(.\k)e|2. (8.11) For a set {(H1(Ak), q5(/\k))}2:=1, the optimal least-squares estimation solution is ecu.) = PTY, (8.12) where Y == [11101) H192) H1(»\n)lT€éR", = [491) 4592) 495226871”. P := ’9T9]_1 - —1 TL = kzl¢ (Ak) )¢T( (Ak) E iffmxm. Since the local coordinates of the process used to model the locally model the concentration of the plume is different than the global coordinates as in Fig. 8.1, a 35 ) E1 Figure 8.1. Local and Global Coordinate Relationship similarity transformation must be used. In this case, the similarity transformation will include a transformation and rotation. Any set of global coordinates (X,, Y,), can be written in local local coordinates (x,, y,) as follows: (Xi — X0)E1+(Yz' — Yo)§2 = 50221 + 9252- (8-13) The rotation between the local and global coordinate systems is handled by 51 _ cos6 sin 6 E1 62 - —sin6 cos6 E2 where 6 is the rotation angle. The following relationship is used to change the coor- , (8.14) dinates: { 11:, = (X,- — X0) cos6 + (Y, — YO) sin6 (8.15) y,- = —(XJ - X0) sin6 + (Y,- — Y0) cos6 ’ where 2:, and y,- are local coordinates and X, and Y,- are global coordinates. 36 CHAPTER 9 Conclusion This paper presents optimal sampling strategies for estimating the parameters of an advection—diffusion process using nonholonomic sensing agents. Two problems are stated, then the components needed to approach these problems are discussed. Simulation results are presented which show that the strategy proposed in Section 6.1 works. The simulation results for strategy proposed in Section 6.2 suggest that a more accurate model of the local plume is needed. Future work will be done to make a numerical closed-form solution to model the local plume behavior. Once the numerical closed-form solution is found, it will be used to trace the concentration field displayed in Fig. 7.6 and other more complicated concentration fields. Future work will also include using other strategies trace the source of realistic advection-diffusion processes. 37 Appendix Sensitivity Coefficients Computation In what follows, the sensitivity matrix is computed, which will require partial deriva- tives with respect to each of parameters. The analytical solution of the PDE equation is given by _ C C(z, y, t) = 0 \/(47T)3 (t — tol3fDxnyy — szy) X 1 ex — 1 p{ 4(Dxnyy — Dgy) (t _ t0) ( ) x [(x — X0 — V3 (t — t0)3(Dxnyy - 092.5)?! m where = 00 (A4703 9 - t0)3(Da:a:Dyy — 0.2g) 1 4(Dmeyy - Dgy) (t " t0) )2" if, if ] [(x — X0 — v9 — t0))l * Co 9::— X [(X -X() - V(t —t0) 38 and 6 * = Dyy >< — —Vt—t T an$( ) 4(Dxnyy _ D%y)2(t __ to) [(X X0 ( 0))l [ 33:, ‘01:? ] X M - X0 — V9 — t0))l X 1 T — x — — V t — t 4(Dxnyy — D,2,,,)(t — to) [0‘ X0 ( 0D] 00 01 8G —D:1::L'Co(t - t0)3 + 0* 6 = * ex *), (3) 8D” (2661703 (t - t0)3(Da::chy — D 1%?!)3 0 aDyy( )) Pf where X X [(X - X0 - V(t - t0))l M = X [(X — X0 - V(t - t0))l 0024 4(5..D,,_Dg,,2(._.,, Dyy ‘Drry x — —Vt—t x —D2y pm [(X X0 ( 0))l D 1 — _ V T _ X t-t 4( ICEDyy - Dgth —- to) [(X X0 ( 0))l X [33 x l>i 36 03900 = (*) ex (*), (4) 6013?! (\/(47r)3 (t _ t0)3(D$nyy _ Dgy)3+ CO BDyy ) p where 6 * : _ny _ _ _ T X [ 33:3, 31:31 x [(X "‘ X0 — V(t - t0))l 1 T - x — —Vt-t x [.01 91]x[(><—Xo-V(t-t0))l 39 ac? 1 = exp(*) (5) 800 \/(47r)3 (2 — t0)3(Dxnyy "' 63;) 370 — CO exp(*) (6) (x-X -V(t-t ))- Using chain—rule m _ £29. 8230 _ 6At 62:0 with At = t — to such that BAt _ ___ _1 Bio and 80 BC' '85— ‘ (‘54?) ‘7) where — = — + C ——(*) expe) O a“ ( 2\/(47r)3At5(Dxnyy — 02,) am ) and 8 —1 T —* = x — -—VAt amt) 4At2(omp,,y_pg,) (x X0 ( )) D -D W 23/ — — V At X 1 4At(Dxnyy — 052.3,) D —D 2V2" W “I (x - X0 — V> -ny Dame X 4O BIBLIOGRAPHY fl] [2] [4] [5] [6] [7] [8] Julius Adler. Chemotaxis in bacteria. Journal of Supramolecular Structure, 4(3)305—317,1966. M. H. G. Amin, S. J. Gibbs, R. J. Chorley, K. S. Richards, T. A. Carpenter, and L. D. Hall. Study of flow and hydrodynamics dispersion in a porous medium using pulsed field gradient magnetic resonance. Proc. R. Soc. Land. A, 453:489—513, 1997. Vassilios N. Chistopouios and Stergios Roumeliotis. Adaptive sensing for in- stantaneous gas release parameter estimation. In Proceeding of the 2005 IEEE International Conference on Robotics and Automation, 2005. David S. DeCroix. Large-eddy simulated of urban dispersion during the UR- BAN2000 field program IOP-lO, 25—26 October 2000. Amit Dhariwal, Gaurav S. Sukhatme, and Aristides A. G. Requicha. Bacterium- inspired robots for environmental monitoring. In Proceedings of the IEEE Inter- national Conference on Robotics and Automation, 2004. E. W. Dijkstra. A note on two problems in connexion with graphs. Numeriche Mathematik, pages 269—271, 1959. A F Emery and A V Nenarokomov. Meas. sci. technol. 9. pages 864-876, 1998. Jay A. Farrell. Chemical plume tracing via an autonomous underwater vehicle. IEEE Journal of Oceanic Engineering, 30(2):428—442, April 2005. 41 [9] [10] [11] [12] [13] [14] [15] [16] [17] Yongxing Hao and Sunil K. Agrawil. Planning and control of ugv formations in a dynamic environment] a practical framework with experiments. Robotics and Autononmous Systems, pages 101—110, 2005. Adam T. Hayes, Alcherio Martinoli, and Rodney M. Goodman. Swarm robotic odor localization. In Proceedings of the 2001 IEEE/RSJ International Conference on Intelligent Robotics and Systems, 2001. Petter Ogren, Edward Fiorelli, and Naomi Ehrich Leonard. Cooperative control of mobile sensor networks: Adaptive gradient climing in a distributed envriron- ment. IEEE Transaction on Automatic Control, 49(8):1292, August 2004. Roger A. Horn and Charles R. Johnson. Matrix Analysis. Cambridge University Press, 1985. H. Ishida, Y. Kagawa, T. Nakamoto, and T. Moriizumi. Odor-source localiza- tion in the clean room by an autonomous mobile sensing system. Sensors and Actuators B, 33:115—121, 1996. H. Ishida, Y. suetsugu, T. Nakamoto, and T. Moriizumi. Study of an au- tonomous mobile sensing system for localization of odor source using gas sensors and anemometric sensors. Sensors and Actuators A, 45:153—157, 1994. P. Kathirgamanathan, R. McKibbin, and RI. McLachlan. Source term estima- tion of pollution from an instantaneous point source. Research Letters in the Information and Mathmatical Sciences, 3(1):59—67, April 2002. Steven M. Kay. Fundamentals of Statistical Signal Processing: Estimation T he- ory. Prentice Hall, 1993. N. Leonardo and A. Robinson. Adaptive sampling and forecast- ing plan AOSN-II MB’03 project. Technical report, Available: http://www.princeton.edu/ dcsl / aosn / documents / ASFP.pdf, 2003. 42 [18] Lennart Ljung and Torsten S'o'derstrc'im. Theory and Practice of Recursive Iden- [19] [20] [21] [22] [23] [24] [25] tification. The MIT Press, Cambridge, Massachusetts, London, England, 1983. Jorg Matthas, Lutz Groll, and Hubert D. Keller. Source localization by spatially distributed electronic noses for advection and diffusion. IEEE Transactions on Signal Processing, pages 1711—1719, May 2005. Nils J. Nilsson. Principles of Artificial Intelligence. Tioga Publishing Co., 1980. Giuseppe Oriolo, Alessandro De Luca, and Marilena Vendittelli. er control via dynamic feedback linearization: Design, implementation, and experimental validation. IEEE Transactions on Control Systems Technology, pages 835—852, November 2002. J. C. Parker and M. T. van Genuchten. Flux-averaged and volume-averaged concentrations in continuum approaches to solute transport. Water Resour. Res, 20:866-872, 1984. Boaz Porat and Arye Nehorai. Localization vapor-emitting sources by moving sensors. IEEE Ti‘ansactions on Signal Processing, 44(4):1018—1021, April 1996. Wei Ren. Consensus seeking in multi-vehicle systems with a time-varying refer- ence state. Proceedings of the 2007 American Control Conference, pages 717—722, July 2007. Wei Ren. Second-order consensus algorithm with extensions to switching topolo- gies and reference models. Proceedings of the 2007 American Control Conference, pages 1431—1436, July 2007. [26] W. S. Smith, M. J. Brown, and D. S. DeCoix. Evaluation of CPD simulations using laboratory data and urban field experiments. 12th Joint Conference on the Applications of Air Pollution Meteorology with the Air and Waste Management Association, 2002. 43 [27] Nathan Sorensen and Wei Ren. A unified formation control scheme with a single [28] [29] or multiple leaders. Proceedings of the 2007 American Control Conference, pages 5412—5418, July 2007. Panos Tzanos, Milos Zefran, and Arye Nehorai. Information based distributed control for biochemical source detection and localization. In Proceedings of the 2005 IEEE International Conference on Robotics and Automation. Dimitri Zarzhitsky, Diana F. Spears, and William M. Spears. Distributed robotics approach to chemical plume tracing. In Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, 2005. 44