- Email: [email protected]
Decentralized mobile sensor routing for parameter estimation of distributed systems
Proceedings of the First IFAC Workshop on Estimation and Control of Networked Systems September 24-26, 2009, Venice, Italy
Decentralized mobile sensor routing for parameter estimation of distributed systems Maciej Patan Institute of Control and Computation Engineering, University of Zielona G´ ora, Poland (e-mail: [email protected]). Abstract: The problem of determining optimal observation strategies for identification of unknown parameters in distributed systems is discussed. Particularly, a setting where the measurement process is performed by collecting spatial data from mobile nodes with sensing capacity forming an organized network is considered. The framework is based on the Fisher information matrix as a measure of the information content in the measurements. The approach is to convert the problem to a canonical optimal control one in Mayer form, in which the control forces of the sensors may be optimized. Then, through an adaptation of some pairwise communication algorithms a numerical scheme is developed, which decompose the resulting problem and distributes the computational burden between the network nodes. Keywords: Distributed parameter systems; mobile sensor networks; parameter estimation; sensor trajectories; distributed control. 1. INTRODUCTION Nowadays, due to the dynamically increasing complexity of modern measurement systems a strong necessity appears for more systematic approaches to the problem of data acquisition from various sources through optimization of observation strategies. This issue is especially related to one of the most general and important classes of systems, namely distributed parameter systems (DPSs), i.e., systems described by partial differential equations (PDEs), since they are widely encountered in numerous practical engineering domains. Although the importance of measurement system design for estimation of unknown coefficients in DPSs has been recognized for a long time, relatively few attempts have been made at solving this problem, cf. surveys in (Kubrusly and Malebranche, 1985; Uci´ nski, 2005; Patan, 2004). It has generated special interest in areas such as the design of air quality monitoring systems, groundwater-resources management, recovery of valuable minerals and hydrocarbon, model calibration in meteorology and oceanography, chemical engineering, hazardous environments and smart materials The operation and control of such systems usually requires precise information on the parameters which dictate the accuracy of the underlying mathematical model, but that information is only available through a limited number of possibly expensive sensors. The underlying idea is to select those sensor locations that lead to the best estimates of the process parameters. The optimality of the locations is judged by an appropriate measure of the estimate-error covariance matrix. In the framework of sensor location problem a decided majority of results communicated in the literature is limited to the selection of stationary sensor positions. A natural generalization which imposes itself is to apply sensors which are able to move in the spatial area providing the
978-3-902661-52-4/09/$20.00 © 2009 IFAC
best information about the parameters. The complexity of the resulting design problem is expected to be compensated by a number of benefits. In particular, sensors are not assigned to fixed positions which are optimal only on the average, but are capable of tracking points which provide at a given time the best information about the parameters to be identified. Consequently, by actively reconfiguring a sensor system we can expect to significantly improve the quality of the measurements or/and decrease the scale of the network. It should be pointed out that the optimal design of moving sensor trajectories is increasingly attracting attention in the context of sensor networks which play a role of importance in the research community, cf. (Zhao and Guibas, 2004; Cassandras and Li, 2005; Jain and Agrawal, 2005) Technological advances in communication systems and the growing ease in making small, low power and inexpensive mobile systems now make it feasible to deploy a group of networked vehicles in a number of ¨ environments, cf. Ogren et al. (2004). A cooperated and scalable network of vehicles, each of them equipped with a single sensor, has the potential to substantially enhance performance of the observation systems. The number of sensor placement techniques developed to manage the problems of a practical scale is very limited (cf. Uci´ nski (2005); Kubrusly and Malebranche (1985); Patan (2004)), nevertheless some effective approaches has been proposed to cover a number of different experimental settings related to parameter estimation, including stationary (Nehorai et al., 1995; Uci´ nski, 2000; Patan and Patan, 2005; Patan and Uci´ nski, 2008; Point et al., 1996), scanning (Patan and Uci´ nski, 2005; Patan, 2006, 2008) or moving observations (Jeremi´c and Nehorai, 2000; Uci´ nski and Chen, 2005; Patan et al., 2008). Furthermore, almost all developed solutions rely on centralized techniques, which assume an existence of some superior entity to
210
10.3182/20090924-3-IT-4005.0021
1st IFAC NECSYS (NECSYS'09) Venice, Italy, September 24-26, 2009
maintain the whole network and responsible for complex optimization of the observation strategy. The distributed nature of the design problem is taken into account very occasionally. Therefore, the purpose of the investigations undertaken here was to establish a practical approach to properly formulate and solve one of such problems, namely the problem of taking account of limited motion capabilities of mobile nodes while guiding them so as to observe the state of a DPS and then estimate its unknown parameters. Motivations come from technical limitations imposed on the vehicles conveying the measurement equipment. A decentralized formulation of the resulting optimal control problem is proposed based on the Fisher Information Matrix (FIM) as a measure of the information content in the measurements. 2. OPTIMAL SENSOR LOCATION PROBLEM Let Ω ⊂ Rn (n = 1, 2 or 3) be a bounded spatial domain with sufficiently smooth boundary Γ, and T = (t0 , tf ] a bounded time interval. Consider a distributed parameter system whose scalar state at a spatial point x ∈ Ω ⊂ Rn and time instant t ∈ T is denoted by y(x, t). System state is described by the partial differential equation ∂y = F x, t, y, θ in Ω × T , (1) ∂t where F is a well-posed, possibly nonlinear, differential operator which involves first- and second-order spatial derivatives and may include terms accounting for forcing inputs specified a priori. The PDE (1) is accompanied by the appropriate boundary and initial conditions B(x, t, y, θ) = 0 on Γ × T, (2) y = y0 in Ω × {t = 0}, (3) respectively, B being an operator acting on the boundary Γ and y0 = y0 (x) a given function. Conditions (2) and (3) complement (1) such that the existence of a sufficiently smooth and unique solution is guaranteed. We assume that the forms of F and B are given explicitly up to an m-dimensional vector of unknown constant parameters θ which must be estimated using observations of the system. The vector θ ∈ Rm is assumed to be estimated from measurements made by N moving sensors over the observation horizon T . We call xj : T → Ωad the trajectory of the j-th sensor, where Ωad ⊂ Ω ∪ Γ is a compact set representing the area where measurements can be made. The statistical model of observations is of the form z j (t) = y(xj (t), t) + ε(xj (t), t), t ∈ T, j = 1, . . . , N, (4) where ε constitutes the measurement noise which is assumed to be zero-mean, Gaussian, spatially uncorrelated and white. Although white noise is a physically impossible process, it constitutes a reasonable approximation to any disturbance whose adjacent samples are uncorrelated at all time instants for which the time increment exceeds some value which is small compared with the time constants of the DPS. The white-noise assumption is consistent with most of the literature on the subject. In the presented framework, the parameter identification problem is usually formulated as follows: Given the model (1)–(3) and the outcomes of the measurements z j along the trajectories xj , j = 1, . . . , N , determine b ∈ Θad (Θad being the set of admissible an estimate θ
parameters) which minimizes the output least-squares fitto-data functional given by N Z X j 2 b = arg min z (t) − y(xj (t), t; ϑ) dt (5) θ ϑ∈Θad
j=1
T
where y solves (1)–(3) when θ replaced by ϑ.
b From (5) it becomes clear that the parameter estimate θ depends on the trajectories xj s. This means that we have the possibility to improve the quality of estimates of the system parameters throughout proper selection of these design variables. To achieve such a goal, a quantitative measure Ψ for comparison of particular trajectories is required. A logical approach is to choose a criterion related to the expected accuracy of the parameter estimates to be obtained from the data collected. Such a measure is usually based on the concept of the Fisher Information Matrix (FIM), cf. (Uci´ nski, 2005), which is widely used in optimum experiment design (OED) theory for lumped systems, see (Walter and Pronzato, 1997; Atkinson et al., 2007). Under some not restrictive assumptions, the inverse of the FIM constitutes a good approximation of the covariance matrix for the estimate of θ (Walter and Pronzato, 1997; Atkinson et al., 2007). For simplicity of notations, let us write 1 x (t) q(t) = ... , ∀ t ∈ T.
(6)
N
x (t) Then, the FIM has the following representation, cf. (Uci´ nski, 2005; Quereshi et al., 1980): N Z X M (q) = g(xj (t), t)g T (xj (t), t) dt, (7) j=1
T
where g(x, t) = ∇ϑ y(x, t; ϑ) ϑ=θ0 denotes the vector of the so-called sensitivity coefficients, θ 0 being a prior estimate to the unknown parameter vector θ (Uci´ nski, 2005). As for a specific form of some scalar function Ψ of the FIM, various options exist, cf. (Walter and Pronzato, 1997; Atkinson et al., 2007), but the most popular criterion to be maximized, called the D-optimality criterion, is the log-determinant of FIM: Ψ[M ] = log det(M ). (8) The resulting D-optimum sensor trajectories leads to the minimum volume of the uncertainty ellipsoid for the estimates. 3. LIMITATIONS ON SENSOR MOVEMENTS We assume that the sensors are conveyed by vehicles whose motions are described by x˙ j (t) = h(xj (t), uj (t)) a.e. on T , n
ℓ
xj (0) = xj0 n
(9)
where the given function h : R × R → R is required to be continuously differentiable; xj0 ∈ Rn defines an initial sensor configuration, and uj : T → Rℓ is the control signal which must satisfy ul ≤ uj (t) ≤ uu a.e. on T (10) for given constant vectors ul and uu , j = 1, . . . , N .
211
1st IFAC NECSYS (NECSYS'09) Venice, Italy, September 24-26, 2009
For each j = 1, . . . , N , given any initial position xj0 and any control signal vector uj , there is a unique absolutely continuous function xj : T → Rn which satisfies (9) a.e. on T . In what follows, we will call it the sensor trajectory corresponding to xj0 and uj . For notational simplicity, in lieu of (9) we shall subsequently use one vector system of ODEs ˙ q(t) = d(q(t), u(t)) a.e. on T , q(0) = q 0 , (11) where 1 1 x0 u (t) (12) u(t) = ... , q 0 = ... ,
xN 0 f (x1 (t), u1 (t)) .. (13) d(q(t), u(t)) = . . N N f (x (t), u (t)) In real applications, apart from sensor dynamics some additional restrictions on the motions have to be imposed. For instance, all sensors should stay within the admissible region Ωad where measurements are allowed. It is assumed that Ωad is a compact set defined as follows: Ωad = {x ∈ Ω ∪ Γ | bi (x) ≤ 0, i = 1, . . . , I} (14) where bi ’s are given continuously differentiable functions. Accordingly, the conditions bi (xj (t)) ≤ 0, ∀ t ∈ T (15) must be fulfilled, where 1 ≤ i ≤ I and 1 ≤ j ≤ N . To shorten notation, after relabeling, we rewrite constraints (15) in the form γl (q(t)) ≤ 0, ∀ t ∈ T, (16) where γl , l = 1, . . . , ν tally with (15), ν = IN . In general, numerous additional constraints can be also considered, e.g. limited energy consumption, path lengths, node distances, etc. However there are usually case-specific and require individual formulation depending on application. uN (t)
4. OPTIMAL CONTROL FORMULATION The goal in the optimal measurement problem is to determine the controls applied to each vehicle with sensor node, which minimize a design criterion Ψ[M (q)] defined on the set of all real-valued information matrices of the form (7) under the constraints (10) on the magnitude of the controls and pathwise state constraints (16). In order to increase the degree of optimality, in our approach we will regard q 0 as the design parameter vector to be chosen in addition to the control signal vector u. In order to guarantee the correctness of such a formulation and further derivations, it is necessary to put some restrictions on the smoothness of sensitivity coefficients g. In what follows, we thus assume the continuity of g and ∂g/∂x. The above formulation can be interpreted as the following optimization problem Problem 1. Find the pair (q 0 , u) ∈ P which maximizes J1 (q 0 , u) = Ψ[M (q)] (17) over the set of feasible pairs P = (q 0 , u) | xj0 ∈ Ωad , uj : T → Rℓ , ul ≤ uj (t) ≤ uu a.e. on T , j = 1, . . . , N , (18) subject to the constraints (11) and (16).
Evidently, its high non-linearity excludes any possibility of finding closed-form formulas for its solution. Accordingly, we must resort to numerical techniques. A number of possibilities exist in this respect (Polak, 1997), but before exploiting them, observe that in spite of its apparently non-classical form, the resulting optimal-control problem can be easily cast as a classical Mayer problem where the performance index is defined only via terminal values of state variables. For notational convenience, define the function svec : Sm → Rm(m+1)/2 , where Sm denotes the subspace of all symmetric matrices in Rm×m , that takes the lower triangular part (the elements only on the main diagonal and below) of a symmetric matrix A and stacks them into a vector a: a = svec(A) = col[A11 , A21 , . . . , Am1 , A22 , A32 , . . . , Am2 , . . . , Amm ]. (19) Similarly, let A = Smat(a) be the symmetric matrix such that svec(Smat(a)) = a for any a ∈ Rm(m+1)/2 . Further, define the matrix-valued function N X Π(q(t), t) = g(xj (t), t)g T (xj (t), t).
(20)
j=1
Setting r : T → Rm(m+1)/2 as the solution of the differential equations ˙ r(t) = svec(Π(q(t), t)), r(0) = 0, (21) we have M (q) = Smat(r(tf )), (22) i.e., maximization of Ψ[M (r)] thus reduces to maximization of a function of the terminal value of the solution to (21). Introducing the augmented state vector q(t) w(t) = , r(t)
(23)
we obtain
q w0 = w(0) = 0 . (24) 0 Then, the equivalent canonical optimal control problem can be cast as Problem 2. (Mayer form). Find a pair (w0 , u) ∈ P¯ which minimizes the performance index J2 (w0 , u) = −ψ(w(tf )) (25) subject to ˙ w(t) = h(w(t), u(t)), w(0) = w0 , (26) γ¯l (w(t)) ≤ 0, ∀ t ∈ T, l = 1, . . . , ν, where P¯ = (w0 , u) | xj0 ∈ Ωad , uj : T → Rℓ , ul ≤ uj (t) ≤ uu a.e. on T , j = 1, . . . , N , (27) and d(q(t), u(t) h(w(t), u(t)) = , (28) svec(Π(q(t), t))
212
ψ(w(t)) = Ψ[Smat(q(t))], γ¯l (w(t)) = γl (q(t)).
(29) (30)
1st IFAC NECSYS (NECSYS'09) Venice, Italy, September 24-26, 2009
The problem formulated in this way is thoroughly treated in optimal control theory and can be solved in centralized manner using numerous existing packages or libraries for numerically solving dynamic optimization problems. On the other hand, it should be pointed out that the resulting problem involves dynamic constraints in the form of a large-scale system of ordinary differential equations of which some are functionally dependent on the solution of the sensitivity equations associated with the PDE modeling the distributed parameter system in question. Thus, even for small networks we encounter a high complexity of optimization process and potential severe difficulties in determining a solution via centralized algorithm. In order to deal with this issue a decentralized computation scheme is required and Problem 2 form a basis for such procedure. This constitutes the contents of the next section. 5. DISTRIBUTED SENSOR ROUTING A key idea underlaying a distributed policy for optimization of sensor trajectories consist in decomposition of the centralized optimal control problem into a set of simpler subproblems related to the separate design of trajectory for each particular sensor node. In the following we assume the asynchronous time model for the sensor configuration process. Let k = 0, 1, 2, . . . be the discrete time index, which partition the continuous configuration time axis into time slots Zk = [zk−1 , zk ). First, owing to the assumption of uncorrelated observations from distinct network nodes each sensor contributes to the FIM given by (7) independently. Furthermore, criterion in Problem 2 depends only on terminal values of the elements of FIM. Thus, the FIM entirely determines the sensor trajectories and is the only global quantity required for problem solution which cannot be calculated independently of other nodes. The combination of this two facts leads to the possibility of problem decomposition if the FIM could be calculated or estimated in a decentralized way. Fortunately, the D-optimality criterion is a homogenous function, therefore without loss of generality we can normalize FIM given by (7) with the number of nodes obtaining the so-called average FIM : N X ¯ (q) = 1 M (xj (t)) (31) M N j=1 where j
M (x (t)) =
Z
g(xj (t), t)g T (xj (t), t) dt.
(32)
T
It is clear, that average FIM is an arithmetic mean of the local information matrices given by (32). Thus, our task is closely related to the problem of distributed averaging on a network which appears in many applications and has been a subject of extensive studies (Bertsekas and Tsitsiklis, 1997; Boyd et al., 2006; Bullo et al., 2009). Distributed averaging can be achieved in many ways. One of straightforward techniques is a pairwise communication flooding, also known as a gossip scheme, which in its classic version assumes that at the k-th time slot the i-th sensor contacts some neighboring node j with probability Pij , i.e. a pair (i − j) is randomly and independently selected. At this time, both nodes set their values equal to the average ¯ ℓ(k) an estimate of of their current values. Denoting as M
global FIM maintained by ℓ-th sensor at time slot Zk we have the following updates for communicating pair (arrow denotes update operator) 1 ¯ (k−1) ¯ (k−1) ¯ (k) Mi , ℓ = i, j (33) +M M ← j ℓ 2 For each node, having updated an estimate of global FIM, the central idea is to determine its trajectory separately treating the contribution of other nodes to the global FIM as fixed variables which are not optimized, as we expect that such a procedure will improve the quality of current guidance policy. Mathematically, by analogy to the Problem 2, introducing the state vector j x (t) , (34) s(t) = p(t) where p : T → Rm(m+1)/2 is the solution to the Cauchy problem ˙ p(t) = svec g(xj (t), t)g T (xj (t), t) , p(0) = 0, (35) we have j (36) s0 = s(0) = x0 . 0 Hence, the control design problem for j-th node can be stated in the following form Problem 3. Find a pair (s0 , uj ) ∈ Pj which minimizes the performance index J3 (s0 , uj ) = −ψ(s(tf )) (37) subject to ˙ s(t) = h(s(t), uj (t)), (38) s(0) = s0 , γ˜l (s(t)) ≤ 0, ∀ t ∈ T, l = 1, . . . , I, where Pj = (s0 , uj ) | xj0 ∈ Ωad , uj : T → Rℓ , ul ≤ uj (t) ≤ uu a.e. on T , (39) and d(xj (t), uj (t) , (40) h(s(t), uj (t)) = svec g(xj (t), t)g T (xj (t), t) ψ(s(t)) = Ψ[ N1 Smat(p(t)) + j
γ˜l (s(t)) = γl (x (t)).
N −1 ¯ (k) N M j ],
(41)
(42)
After solving the Problem 3 and determining the controls for given network node, its estimate of global FIM can be updated according to the scheme 1 N − 1 ¯ (k) ¯ (k) Mj + Smat(p(tf )). (43) M ← j N N The first term is expected to ensure consensus among the nodes (represents the average information from the rest of network), while the second accounts for the increase in the total contribution of the considered node. Obviously, the solutions of the Problem 3 are dependent on the current estimates of global FIM stored in network nodes. Therefore, a very natural idea is to alternate design and estimation steps. At each configuration interval Zk , in turn, local sensor trajectories are determined based on the available estimates of FIM, then FIM is appropriately updated and propagated through network via pairwise communication process. In this general scheme embodied by Algorithm 1, it is supposed that each estimation phase
213
1st IFAC NECSYS (NECSYS'09) Venice, Italy, September 24-26, 2009
adaptively improves our knowledge about the FIM and this knowledge can then be used to improve the quality of the next design to be performed. Algorithm 1 Adaptive distributed sensor routing. Indexes j and i denote, respectively, data from local repository and obtained from neighbor. Function EXCHANGE is responsible for both sending and receiving data to/from connected neighbor (order depending on who initiated communication) and SolveOCP solves Problem 3 for given FIM estimate. 1: procedure Adaptive routing 2: Set k = 0. Set arbitrarily uj (t) and xj0 . ⊲ Initialization ¯ (k) = M(xj (t)). 3: Calculate xj (t) and M j 4: while Termination condition do ⊲ Repeat sequentially 5: Connect to the random neighbor i. ¯ (k) ,M ¯ (k) ) 6: Exchange(M ⊲ Sends and receives FIM j i (k) (k) (k) 1 ¯ ¯ ¯ 7: Update M ← M +M j
2
j
i
(k)
¯ 8: (xj0 , uj (t))=SolveOCP(M ⊲ Trajectory design j ) 9: Calculate xj (t) and M(xj (t)). ¯ (k) ← N −1 M ¯ (k) + 1 M(xj (t)) 10: Update M j j N N 11: k ←k+1 12: end while 13: end procedure
6. ILLUSTRATIVE EXAMPLE In order to demonstrate the developed method consider an example on transport-chemistry process of atmospheric pollutant over an urban area being normalized to a unit square. At the point x0 = (0.3, 0.6) an active source is located which emits the pollutant substance which is later spread over the area according to the combination of different processes such as turbulences, chemical decay and wind with velocity field v(x, t) = 2(y + 1/4), 2(t − x) varying in space and time (cf. Fig. 1). Mathematically, the process over the normalized observation interval T = [0, 1] can be described by the following advection-diffusionreaction PDE: ∂y(x, t) +∇ · v(x, t)y(x, t) (44) ∂t =∇ · κ(x)∇y(x, t) + η(x), x ∈ Ω subject to the following boundary and initial conditions: ∂y(x, t) = 0 on Γ × T, y(x, 0) = 0 in Ω, (45) ∂n 2 where the term η(x) = e−100kx−x0 k represents an intensity of pollution source and ∂y/∂n stands for the partial derivative of y with respect to the outward normal to the boundary Γ. In our illustrative simulations, the following form of the turbulent diffusivity coefficient was applied κ(x) = θ1 + θ2 x21 + θ3 x22 , (46) so parameters θ = (θ1 , θ2 , θ3 ) need to be estimated based on measurement data. Our goal is to determine the D-optimal trajectories for network of N = 4 movable sensor nodes. In order to verify the proposed decentralized approach, a Matlab program was written using a PC equipped with Centrino T9300 processor (2.5GHz, 3 GB RAM) running Windows Vista and Matlab 2007b. The nominal values of the systems parameters were assumed to be θ10 = 0.02 and θ20 = θ30 = 0.005. The system of PDEs was solved using COMSOL 3.5
environment on the spatial mesh composed of 682 triangles and an evenly partitioned time interval (30 subintervals). The sensitivity coefficients were then linearly interpolated and stored within the sensor nodes repositories. Finally, to adaptively solve the optimal control problems, the package RIOTS 95(Schwartz et al., 1997) was applied, which is designed as a Matlab toolbox suitable for solving a very broad class of control problems. In each configuration stage 10 pairwise averaging steps were performed before trajectory design and the uniform probability distribution for communications between the nodes was assumed. As regard to sensor dynamics a following simple model was adopted ˙ q(t) = u(t), q(0) = q 0 , and u was bounded according to |ui (t)| ≤ 0.3, ∀t ∈ T. The total simulation of k = 50 asynchronous configuration steps took about 4 minutes of computation time. The exemplary results are presented in Fig. 2 where we have trajectories calculated at consecutive distributed adaptation stages compared with D-optimal trajectories determined in a fully centralized way via direct solution of Problem 2. We can observe how sensors try to follow the complex pollutant concentration changes leading to quite sophisticated sensor motions. Also we can see that the Algorithm 1 after a relatively small configuration time k provides a solution (cf. Fig. 2(b)-(c)) which is very similar to the reference trajectory design shown in Fig. 2(d). This is because of small number of nodes, leading to very effective propagation of FIM estimates through the network. 7. CONCLUSION The problem of sensor routing was addressed for a monitoring network with mobile nodes providing informative observational data for the purpose of parameter estimation in DPSs. It was demonstrated that the relevant optimization task can be fit to the framework of an optimal control problem. Furthermore, a decentralized procedure for its solution is developed dedicated for scalable sensor networks. A proper reformulation of control problem in Mayer form makes it possible to directly apply numerous existing solvers for dynamic optimization. ACKNOWLEDGEMENTS The work was supported by the Polish Ministry of Science and Higher Education under Grant N N519 2971 33. REFERENCES Atkinson, A.C., Donev, A.N., and Tobias, R. (2007). Optimum Experimental Design, with SAS. Oxford University Press, Oxford. Bertsekas, D.P. and Tsitsiklis, J.N. (1997). Parallel and Distributed Computation: Numerical Methods. Athena Scientific, Belmont, MA. Boyd, S., Ghosh, A., Prabhakar, B., and Shah, D. (2006). Randomized gossip algorithms. IEEE Transactions on Information Theory, 52(6), 2508–2530. Bullo, F., Cort´ es, J., and Mart´ınez, S. (2009). Distributed Control of Robotic Networks. Applied Mathematics Series. Princeton University Press. To appear. Electronically available at http://coordinationbook.info. Cassandras, C.G. and Li, W. (2005). Sensor networks and cooperative control. European Journal of Control, 11(4–5), 436–463.
214
1st IFAC NECSYS (NECSYS'09) Venice, Italy, September 24-26, 2009
0.8
0.8
0.8
0.8
0.6
0.6
0.6
0.6
x2
1
x2
1
x2
1
x2
1
0.4
0.4
0.4
0.4
0.2
0.2
0.2
0.2
0 0
0.2
0.4
0.6
0.8
0 0
1
0.2
0.4
x
0.6
0.8
0 0
1
0.2
0.4
x
1
(a) t = 0.05
0.6
0.8
0 0
1
0.2
0.4
x
1
(b) t = 0.3
0.6
0.8
1
x
1
1
(c) t = 0.7
(d) t = 1
Fig. 1. Temporal changes in the wind velocity field and pollutant concentration (star indicates the location of source).
0.8
0.8
0.8
0.8
0.6
0.6
0.6
0.6
x2
1
x2
1
x2
1
x2
1
0.4
0.4
0.4
0.4
0.2
0.2
0.2
0.2
0 0
0.2
0.4
0.6
x
1
(a) k = 10
0.8
1
0 0
0.2
0.4
0.6
0.8
0 0
1
x
0.2
0.4
0.6
x
1
1
(b) k = 20
(c) k = 40
0.8
1
0 0
0.2
0.4
0.6
0.8
1
x
1
(d) Centralized approach
Fig. 2. Approximations of D-optimal sensor trajectories obtained in different stages of distributed configuration scheme (Algorithm 1) (a)-(c) versus trajectories determined in a centralized fashion (d). Open circles indicate the initial positions of mobile nodes. Jain, N. and Agrawal, D.P. (2005). Current trends in wireless sensor network design. International Journal of Distributed Sensor Networks, 1, 101–122. Jeremi´ c, A. and Nehorai, A. (2000). Landmine detection and localization using chemical sensor array processing. IEEE Transactions on Signal Processing, 48(5), 1295–1305. Kubrusly, C.S. and Malebranche, H. (1985). Sensors and controllers location in distributed systems — A survey. Automatica, 21(2), 117–128. Nehorai, A., Porat, B., and Paldi, E. (1995). Detection and localization of vapor-emitting sources. IEEE Transactions on Signal Processing, 43(1), 243–253. ¨ Ogren, P., Fiorelli, E., and Leonard, N.E. (2004). Cooperative control of mobile sensor networks: Adaptive gradient climbing in a distributed environment. IEEE Transactions on Automatic Control, 49(8), 1292–1302. Patan, M. (2004). Optimal Observation Strategies for Parameter Estimation of Distributed Systems. Zielona G´ ora: University Press. May be viewed at http://www.zbc.zgora.pl. Patan, M. (2006). Optimal activation policies for continous scanning observations in parameter estimation of distributed systems. International Journal of Systems Science, 37(11), 763–775. Patan, M. (2008). A parallel sensor scheduling technique for fault detection in distributed parameter systems. In Lecture Notes in Computer Science: Euro-Par 2008 : Parallel Procesing, volume 5168, 833–843. Springer-Verlag, Berlin Heidelberg. Patan, M., Chen, Y., and Tricaud, C. (2008). Resource-constrained sensor routing for parameter estimation of distributed systems. In Proc. 17th IFAC World Congress, Seoul, South Korea, 4–8 July 2008. Patan, M. and Patan, K. (2005). Optimal observation strategies for model-based fault detection in distributed systems. International Journal of Control, 78(18), 1497–1510. Patan, M. and Uci´ nski, D. (2005). Optimal activation strategy of discrete scanning sensors for fault detection in distributedparameter systems. In Proc. 16th IFAC World Congress, Prague,
Czech Republic, 4–8 July 2005. Published on CD-ROM. Patan, M. and Uci´ nski, D. (2008). Configuring a sensor network for fault detection in distributed parameter systems. International Journal of Applied Mathematics and Computer Science, 18(4), 513–524. Point, N., Vande Wouwer, A., and Remy, M. (1996). Practical issues in distributed parameter estimation: Gradient computation and optimal experiment design. Control Engineering Practice, 4(11), 1553–1562. Polak, E. (1997). Optimization. Algorithms and Consistent Approximations. Applied Mathematical Sciences. Springer-Verlag, New York. Quereshi, Z.H., Ng, T.S., and Goodwin, G.C. (1980). Optimum experimental design for identification of distributed parameter systems. International Journal of Control, 31(1), 21–29. Schwartz, A.L., Polak, E., and Chen, Y. (1997). A Matlab Toolbox for Solving Optimal Control Problems. Version 1.0 for Windows. URL http://www.schwartz-home.com/~adam/RIOTS/. Uci´ nski, D. (2000). Optimal selection of measurement locations for parameter estimation in distributed processes. International Journal of Applied Mathematics and Computer Science, 10(2), 357–379. Uci´ nski, D. (2005). Optimal Measurement Methods for DistributedParameter System Identification. CRC Press, Boca Raton, FL. Uci´ nski, D. and Chen, Y. (2005). Time-optimal path planning of moving sensors for parameter estimation of distributed systems. In Proc. 44th IEEE Conference on Decision and Control, and the European Control Conference 2005, Seville, Spain. Published on CD-ROM. ´ and Pronzato, L. (1997). Identification of Parametric Walter, E. Models from Experimental Data. Communications and Control Engineering. Springer-Verlag, Berlin. Zhao, F. and Guibas, L.J. (2004). Wireless Sensor Networks: An Information Processing Approach. Morgan Kaufmann Publishers, Amsterdam.
215
Decentralized mobile sensor routing for parameter estimation of distributed systems Maciej Patan Institute of Control and Computation Engineering, University of Zielona G´ ora, Poland (e-mail: [email protected]). Abstract: The problem of determining optimal observation strategies for identification of unknown parameters in distributed systems is discussed. Particularly, a setting where the measurement process is performed by collecting spatial data from mobile nodes with sensing capacity forming an organized network is considered. The framework is based on the Fisher information matrix as a measure of the information content in the measurements. The approach is to convert the problem to a canonical optimal control one in Mayer form, in which the control forces of the sensors may be optimized. Then, through an adaptation of some pairwise communication algorithms a numerical scheme is developed, which decompose the resulting problem and distributes the computational burden between the network nodes. Keywords: Distributed parameter systems; mobile sensor networks; parameter estimation; sensor trajectories; distributed control. 1. INTRODUCTION Nowadays, due to the dynamically increasing complexity of modern measurement systems a strong necessity appears for more systematic approaches to the problem of data acquisition from various sources through optimization of observation strategies. This issue is especially related to one of the most general and important classes of systems, namely distributed parameter systems (DPSs), i.e., systems described by partial differential equations (PDEs), since they are widely encountered in numerous practical engineering domains. Although the importance of measurement system design for estimation of unknown coefficients in DPSs has been recognized for a long time, relatively few attempts have been made at solving this problem, cf. surveys in (Kubrusly and Malebranche, 1985; Uci´ nski, 2005; Patan, 2004). It has generated special interest in areas such as the design of air quality monitoring systems, groundwater-resources management, recovery of valuable minerals and hydrocarbon, model calibration in meteorology and oceanography, chemical engineering, hazardous environments and smart materials The operation and control of such systems usually requires precise information on the parameters which dictate the accuracy of the underlying mathematical model, but that information is only available through a limited number of possibly expensive sensors. The underlying idea is to select those sensor locations that lead to the best estimates of the process parameters. The optimality of the locations is judged by an appropriate measure of the estimate-error covariance matrix. In the framework of sensor location problem a decided majority of results communicated in the literature is limited to the selection of stationary sensor positions. A natural generalization which imposes itself is to apply sensors which are able to move in the spatial area providing the
978-3-902661-52-4/09/$20.00 © 2009 IFAC
best information about the parameters. The complexity of the resulting design problem is expected to be compensated by a number of benefits. In particular, sensors are not assigned to fixed positions which are optimal only on the average, but are capable of tracking points which provide at a given time the best information about the parameters to be identified. Consequently, by actively reconfiguring a sensor system we can expect to significantly improve the quality of the measurements or/and decrease the scale of the network. It should be pointed out that the optimal design of moving sensor trajectories is increasingly attracting attention in the context of sensor networks which play a role of importance in the research community, cf. (Zhao and Guibas, 2004; Cassandras and Li, 2005; Jain and Agrawal, 2005) Technological advances in communication systems and the growing ease in making small, low power and inexpensive mobile systems now make it feasible to deploy a group of networked vehicles in a number of ¨ environments, cf. Ogren et al. (2004). A cooperated and scalable network of vehicles, each of them equipped with a single sensor, has the potential to substantially enhance performance of the observation systems. The number of sensor placement techniques developed to manage the problems of a practical scale is very limited (cf. Uci´ nski (2005); Kubrusly and Malebranche (1985); Patan (2004)), nevertheless some effective approaches has been proposed to cover a number of different experimental settings related to parameter estimation, including stationary (Nehorai et al., 1995; Uci´ nski, 2000; Patan and Patan, 2005; Patan and Uci´ nski, 2008; Point et al., 1996), scanning (Patan and Uci´ nski, 2005; Patan, 2006, 2008) or moving observations (Jeremi´c and Nehorai, 2000; Uci´ nski and Chen, 2005; Patan et al., 2008). Furthermore, almost all developed solutions rely on centralized techniques, which assume an existence of some superior entity to
210
10.3182/20090924-3-IT-4005.0021
1st IFAC NECSYS (NECSYS'09) Venice, Italy, September 24-26, 2009
maintain the whole network and responsible for complex optimization of the observation strategy. The distributed nature of the design problem is taken into account very occasionally. Therefore, the purpose of the investigations undertaken here was to establish a practical approach to properly formulate and solve one of such problems, namely the problem of taking account of limited motion capabilities of mobile nodes while guiding them so as to observe the state of a DPS and then estimate its unknown parameters. Motivations come from technical limitations imposed on the vehicles conveying the measurement equipment. A decentralized formulation of the resulting optimal control problem is proposed based on the Fisher Information Matrix (FIM) as a measure of the information content in the measurements. 2. OPTIMAL SENSOR LOCATION PROBLEM Let Ω ⊂ Rn (n = 1, 2 or 3) be a bounded spatial domain with sufficiently smooth boundary Γ, and T = (t0 , tf ] a bounded time interval. Consider a distributed parameter system whose scalar state at a spatial point x ∈ Ω ⊂ Rn and time instant t ∈ T is denoted by y(x, t). System state is described by the partial differential equation ∂y = F x, t, y, θ in Ω × T , (1) ∂t where F is a well-posed, possibly nonlinear, differential operator which involves first- and second-order spatial derivatives and may include terms accounting for forcing inputs specified a priori. The PDE (1) is accompanied by the appropriate boundary and initial conditions B(x, t, y, θ) = 0 on Γ × T, (2) y = y0 in Ω × {t = 0}, (3) respectively, B being an operator acting on the boundary Γ and y0 = y0 (x) a given function. Conditions (2) and (3) complement (1) such that the existence of a sufficiently smooth and unique solution is guaranteed. We assume that the forms of F and B are given explicitly up to an m-dimensional vector of unknown constant parameters θ which must be estimated using observations of the system. The vector θ ∈ Rm is assumed to be estimated from measurements made by N moving sensors over the observation horizon T . We call xj : T → Ωad the trajectory of the j-th sensor, where Ωad ⊂ Ω ∪ Γ is a compact set representing the area where measurements can be made. The statistical model of observations is of the form z j (t) = y(xj (t), t) + ε(xj (t), t), t ∈ T, j = 1, . . . , N, (4) where ε constitutes the measurement noise which is assumed to be zero-mean, Gaussian, spatially uncorrelated and white. Although white noise is a physically impossible process, it constitutes a reasonable approximation to any disturbance whose adjacent samples are uncorrelated at all time instants for which the time increment exceeds some value which is small compared with the time constants of the DPS. The white-noise assumption is consistent with most of the literature on the subject. In the presented framework, the parameter identification problem is usually formulated as follows: Given the model (1)–(3) and the outcomes of the measurements z j along the trajectories xj , j = 1, . . . , N , determine b ∈ Θad (Θad being the set of admissible an estimate θ
parameters) which minimizes the output least-squares fitto-data functional given by N Z X j 2 b = arg min z (t) − y(xj (t), t; ϑ) dt (5) θ ϑ∈Θad
j=1
T
where y solves (1)–(3) when θ replaced by ϑ.
b From (5) it becomes clear that the parameter estimate θ depends on the trajectories xj s. This means that we have the possibility to improve the quality of estimates of the system parameters throughout proper selection of these design variables. To achieve such a goal, a quantitative measure Ψ for comparison of particular trajectories is required. A logical approach is to choose a criterion related to the expected accuracy of the parameter estimates to be obtained from the data collected. Such a measure is usually based on the concept of the Fisher Information Matrix (FIM), cf. (Uci´ nski, 2005), which is widely used in optimum experiment design (OED) theory for lumped systems, see (Walter and Pronzato, 1997; Atkinson et al., 2007). Under some not restrictive assumptions, the inverse of the FIM constitutes a good approximation of the covariance matrix for the estimate of θ (Walter and Pronzato, 1997; Atkinson et al., 2007). For simplicity of notations, let us write 1 x (t) q(t) = ... , ∀ t ∈ T.
(6)
N
x (t) Then, the FIM has the following representation, cf. (Uci´ nski, 2005; Quereshi et al., 1980): N Z X M (q) = g(xj (t), t)g T (xj (t), t) dt, (7) j=1
T
where g(x, t) = ∇ϑ y(x, t; ϑ) ϑ=θ0 denotes the vector of the so-called sensitivity coefficients, θ 0 being a prior estimate to the unknown parameter vector θ (Uci´ nski, 2005). As for a specific form of some scalar function Ψ of the FIM, various options exist, cf. (Walter and Pronzato, 1997; Atkinson et al., 2007), but the most popular criterion to be maximized, called the D-optimality criterion, is the log-determinant of FIM: Ψ[M ] = log det(M ). (8) The resulting D-optimum sensor trajectories leads to the minimum volume of the uncertainty ellipsoid for the estimates. 3. LIMITATIONS ON SENSOR MOVEMENTS We assume that the sensors are conveyed by vehicles whose motions are described by x˙ j (t) = h(xj (t), uj (t)) a.e. on T , n
ℓ
xj (0) = xj0 n
(9)
where the given function h : R × R → R is required to be continuously differentiable; xj0 ∈ Rn defines an initial sensor configuration, and uj : T → Rℓ is the control signal which must satisfy ul ≤ uj (t) ≤ uu a.e. on T (10) for given constant vectors ul and uu , j = 1, . . . , N .
211
1st IFAC NECSYS (NECSYS'09) Venice, Italy, September 24-26, 2009
For each j = 1, . . . , N , given any initial position xj0 and any control signal vector uj , there is a unique absolutely continuous function xj : T → Rn which satisfies (9) a.e. on T . In what follows, we will call it the sensor trajectory corresponding to xj0 and uj . For notational simplicity, in lieu of (9) we shall subsequently use one vector system of ODEs ˙ q(t) = d(q(t), u(t)) a.e. on T , q(0) = q 0 , (11) where 1 1 x0 u (t) (12) u(t) = ... , q 0 = ... ,
xN 0 f (x1 (t), u1 (t)) .. (13) d(q(t), u(t)) = . . N N f (x (t), u (t)) In real applications, apart from sensor dynamics some additional restrictions on the motions have to be imposed. For instance, all sensors should stay within the admissible region Ωad where measurements are allowed. It is assumed that Ωad is a compact set defined as follows: Ωad = {x ∈ Ω ∪ Γ | bi (x) ≤ 0, i = 1, . . . , I} (14) where bi ’s are given continuously differentiable functions. Accordingly, the conditions bi (xj (t)) ≤ 0, ∀ t ∈ T (15) must be fulfilled, where 1 ≤ i ≤ I and 1 ≤ j ≤ N . To shorten notation, after relabeling, we rewrite constraints (15) in the form γl (q(t)) ≤ 0, ∀ t ∈ T, (16) where γl , l = 1, . . . , ν tally with (15), ν = IN . In general, numerous additional constraints can be also considered, e.g. limited energy consumption, path lengths, node distances, etc. However there are usually case-specific and require individual formulation depending on application. uN (t)
4. OPTIMAL CONTROL FORMULATION The goal in the optimal measurement problem is to determine the controls applied to each vehicle with sensor node, which minimize a design criterion Ψ[M (q)] defined on the set of all real-valued information matrices of the form (7) under the constraints (10) on the magnitude of the controls and pathwise state constraints (16). In order to increase the degree of optimality, in our approach we will regard q 0 as the design parameter vector to be chosen in addition to the control signal vector u. In order to guarantee the correctness of such a formulation and further derivations, it is necessary to put some restrictions on the smoothness of sensitivity coefficients g. In what follows, we thus assume the continuity of g and ∂g/∂x. The above formulation can be interpreted as the following optimization problem Problem 1. Find the pair (q 0 , u) ∈ P which maximizes J1 (q 0 , u) = Ψ[M (q)] (17) over the set of feasible pairs P = (q 0 , u) | xj0 ∈ Ωad , uj : T → Rℓ , ul ≤ uj (t) ≤ uu a.e. on T , j = 1, . . . , N , (18) subject to the constraints (11) and (16).
Evidently, its high non-linearity excludes any possibility of finding closed-form formulas for its solution. Accordingly, we must resort to numerical techniques. A number of possibilities exist in this respect (Polak, 1997), but before exploiting them, observe that in spite of its apparently non-classical form, the resulting optimal-control problem can be easily cast as a classical Mayer problem where the performance index is defined only via terminal values of state variables. For notational convenience, define the function svec : Sm → Rm(m+1)/2 , where Sm denotes the subspace of all symmetric matrices in Rm×m , that takes the lower triangular part (the elements only on the main diagonal and below) of a symmetric matrix A and stacks them into a vector a: a = svec(A) = col[A11 , A21 , . . . , Am1 , A22 , A32 , . . . , Am2 , . . . , Amm ]. (19) Similarly, let A = Smat(a) be the symmetric matrix such that svec(Smat(a)) = a for any a ∈ Rm(m+1)/2 . Further, define the matrix-valued function N X Π(q(t), t) = g(xj (t), t)g T (xj (t), t).
(20)
j=1
Setting r : T → Rm(m+1)/2 as the solution of the differential equations ˙ r(t) = svec(Π(q(t), t)), r(0) = 0, (21) we have M (q) = Smat(r(tf )), (22) i.e., maximization of Ψ[M (r)] thus reduces to maximization of a function of the terminal value of the solution to (21). Introducing the augmented state vector q(t) w(t) = , r(t)
(23)
we obtain
q w0 = w(0) = 0 . (24) 0 Then, the equivalent canonical optimal control problem can be cast as Problem 2. (Mayer form). Find a pair (w0 , u) ∈ P¯ which minimizes the performance index J2 (w0 , u) = −ψ(w(tf )) (25) subject to ˙ w(t) = h(w(t), u(t)), w(0) = w0 , (26) γ¯l (w(t)) ≤ 0, ∀ t ∈ T, l = 1, . . . , ν, where P¯ = (w0 , u) | xj0 ∈ Ωad , uj : T → Rℓ , ul ≤ uj (t) ≤ uu a.e. on T , j = 1, . . . , N , (27) and d(q(t), u(t) h(w(t), u(t)) = , (28) svec(Π(q(t), t))
212
ψ(w(t)) = Ψ[Smat(q(t))], γ¯l (w(t)) = γl (q(t)).
(29) (30)
1st IFAC NECSYS (NECSYS'09) Venice, Italy, September 24-26, 2009
The problem formulated in this way is thoroughly treated in optimal control theory and can be solved in centralized manner using numerous existing packages or libraries for numerically solving dynamic optimization problems. On the other hand, it should be pointed out that the resulting problem involves dynamic constraints in the form of a large-scale system of ordinary differential equations of which some are functionally dependent on the solution of the sensitivity equations associated with the PDE modeling the distributed parameter system in question. Thus, even for small networks we encounter a high complexity of optimization process and potential severe difficulties in determining a solution via centralized algorithm. In order to deal with this issue a decentralized computation scheme is required and Problem 2 form a basis for such procedure. This constitutes the contents of the next section. 5. DISTRIBUTED SENSOR ROUTING A key idea underlaying a distributed policy for optimization of sensor trajectories consist in decomposition of the centralized optimal control problem into a set of simpler subproblems related to the separate design of trajectory for each particular sensor node. In the following we assume the asynchronous time model for the sensor configuration process. Let k = 0, 1, 2, . . . be the discrete time index, which partition the continuous configuration time axis into time slots Zk = [zk−1 , zk ). First, owing to the assumption of uncorrelated observations from distinct network nodes each sensor contributes to the FIM given by (7) independently. Furthermore, criterion in Problem 2 depends only on terminal values of the elements of FIM. Thus, the FIM entirely determines the sensor trajectories and is the only global quantity required for problem solution which cannot be calculated independently of other nodes. The combination of this two facts leads to the possibility of problem decomposition if the FIM could be calculated or estimated in a decentralized way. Fortunately, the D-optimality criterion is a homogenous function, therefore without loss of generality we can normalize FIM given by (7) with the number of nodes obtaining the so-called average FIM : N X ¯ (q) = 1 M (xj (t)) (31) M N j=1 where j
M (x (t)) =
Z
g(xj (t), t)g T (xj (t), t) dt.
(32)
T
It is clear, that average FIM is an arithmetic mean of the local information matrices given by (32). Thus, our task is closely related to the problem of distributed averaging on a network which appears in many applications and has been a subject of extensive studies (Bertsekas and Tsitsiklis, 1997; Boyd et al., 2006; Bullo et al., 2009). Distributed averaging can be achieved in many ways. One of straightforward techniques is a pairwise communication flooding, also known as a gossip scheme, which in its classic version assumes that at the k-th time slot the i-th sensor contacts some neighboring node j with probability Pij , i.e. a pair (i − j) is randomly and independently selected. At this time, both nodes set their values equal to the average ¯ ℓ(k) an estimate of of their current values. Denoting as M
global FIM maintained by ℓ-th sensor at time slot Zk we have the following updates for communicating pair (arrow denotes update operator) 1 ¯ (k−1) ¯ (k−1) ¯ (k) Mi , ℓ = i, j (33) +M M ← j ℓ 2 For each node, having updated an estimate of global FIM, the central idea is to determine its trajectory separately treating the contribution of other nodes to the global FIM as fixed variables which are not optimized, as we expect that such a procedure will improve the quality of current guidance policy. Mathematically, by analogy to the Problem 2, introducing the state vector j x (t) , (34) s(t) = p(t) where p : T → Rm(m+1)/2 is the solution to the Cauchy problem ˙ p(t) = svec g(xj (t), t)g T (xj (t), t) , p(0) = 0, (35) we have j (36) s0 = s(0) = x0 . 0 Hence, the control design problem for j-th node can be stated in the following form Problem 3. Find a pair (s0 , uj ) ∈ Pj which minimizes the performance index J3 (s0 , uj ) = −ψ(s(tf )) (37) subject to ˙ s(t) = h(s(t), uj (t)), (38) s(0) = s0 , γ˜l (s(t)) ≤ 0, ∀ t ∈ T, l = 1, . . . , I, where Pj = (s0 , uj ) | xj0 ∈ Ωad , uj : T → Rℓ , ul ≤ uj (t) ≤ uu a.e. on T , (39) and d(xj (t), uj (t) , (40) h(s(t), uj (t)) = svec g(xj (t), t)g T (xj (t), t) ψ(s(t)) = Ψ[ N1 Smat(p(t)) + j
γ˜l (s(t)) = γl (x (t)).
N −1 ¯ (k) N M j ],
(41)
(42)
After solving the Problem 3 and determining the controls for given network node, its estimate of global FIM can be updated according to the scheme 1 N − 1 ¯ (k) ¯ (k) Mj + Smat(p(tf )). (43) M ← j N N The first term is expected to ensure consensus among the nodes (represents the average information from the rest of network), while the second accounts for the increase in the total contribution of the considered node. Obviously, the solutions of the Problem 3 are dependent on the current estimates of global FIM stored in network nodes. Therefore, a very natural idea is to alternate design and estimation steps. At each configuration interval Zk , in turn, local sensor trajectories are determined based on the available estimates of FIM, then FIM is appropriately updated and propagated through network via pairwise communication process. In this general scheme embodied by Algorithm 1, it is supposed that each estimation phase
213
1st IFAC NECSYS (NECSYS'09) Venice, Italy, September 24-26, 2009
adaptively improves our knowledge about the FIM and this knowledge can then be used to improve the quality of the next design to be performed. Algorithm 1 Adaptive distributed sensor routing. Indexes j and i denote, respectively, data from local repository and obtained from neighbor. Function EXCHANGE is responsible for both sending and receiving data to/from connected neighbor (order depending on who initiated communication) and SolveOCP solves Problem 3 for given FIM estimate. 1: procedure Adaptive routing 2: Set k = 0. Set arbitrarily uj (t) and xj0 . ⊲ Initialization ¯ (k) = M(xj (t)). 3: Calculate xj (t) and M j 4: while Termination condition do ⊲ Repeat sequentially 5: Connect to the random neighbor i. ¯ (k) ,M ¯ (k) ) 6: Exchange(M ⊲ Sends and receives FIM j i (k) (k) (k) 1 ¯ ¯ ¯ 7: Update M ← M +M j
2
j
i
(k)
¯ 8: (xj0 , uj (t))=SolveOCP(M ⊲ Trajectory design j ) 9: Calculate xj (t) and M(xj (t)). ¯ (k) ← N −1 M ¯ (k) + 1 M(xj (t)) 10: Update M j j N N 11: k ←k+1 12: end while 13: end procedure
6. ILLUSTRATIVE EXAMPLE In order to demonstrate the developed method consider an example on transport-chemistry process of atmospheric pollutant over an urban area being normalized to a unit square. At the point x0 = (0.3, 0.6) an active source is located which emits the pollutant substance which is later spread over the area according to the combination of different processes such as turbulences, chemical decay and wind with velocity field v(x, t) = 2(y + 1/4), 2(t − x) varying in space and time (cf. Fig. 1). Mathematically, the process over the normalized observation interval T = [0, 1] can be described by the following advection-diffusionreaction PDE: ∂y(x, t) +∇ · v(x, t)y(x, t) (44) ∂t =∇ · κ(x)∇y(x, t) + η(x), x ∈ Ω subject to the following boundary and initial conditions: ∂y(x, t) = 0 on Γ × T, y(x, 0) = 0 in Ω, (45) ∂n 2 where the term η(x) = e−100kx−x0 k represents an intensity of pollution source and ∂y/∂n stands for the partial derivative of y with respect to the outward normal to the boundary Γ. In our illustrative simulations, the following form of the turbulent diffusivity coefficient was applied κ(x) = θ1 + θ2 x21 + θ3 x22 , (46) so parameters θ = (θ1 , θ2 , θ3 ) need to be estimated based on measurement data. Our goal is to determine the D-optimal trajectories for network of N = 4 movable sensor nodes. In order to verify the proposed decentralized approach, a Matlab program was written using a PC equipped with Centrino T9300 processor (2.5GHz, 3 GB RAM) running Windows Vista and Matlab 2007b. The nominal values of the systems parameters were assumed to be θ10 = 0.02 and θ20 = θ30 = 0.005. The system of PDEs was solved using COMSOL 3.5
environment on the spatial mesh composed of 682 triangles and an evenly partitioned time interval (30 subintervals). The sensitivity coefficients were then linearly interpolated and stored within the sensor nodes repositories. Finally, to adaptively solve the optimal control problems, the package RIOTS 95(Schwartz et al., 1997) was applied, which is designed as a Matlab toolbox suitable for solving a very broad class of control problems. In each configuration stage 10 pairwise averaging steps were performed before trajectory design and the uniform probability distribution for communications between the nodes was assumed. As regard to sensor dynamics a following simple model was adopted ˙ q(t) = u(t), q(0) = q 0 , and u was bounded according to |ui (t)| ≤ 0.3, ∀t ∈ T. The total simulation of k = 50 asynchronous configuration steps took about 4 minutes of computation time. The exemplary results are presented in Fig. 2 where we have trajectories calculated at consecutive distributed adaptation stages compared with D-optimal trajectories determined in a fully centralized way via direct solution of Problem 2. We can observe how sensors try to follow the complex pollutant concentration changes leading to quite sophisticated sensor motions. Also we can see that the Algorithm 1 after a relatively small configuration time k provides a solution (cf. Fig. 2(b)-(c)) which is very similar to the reference trajectory design shown in Fig. 2(d). This is because of small number of nodes, leading to very effective propagation of FIM estimates through the network. 7. CONCLUSION The problem of sensor routing was addressed for a monitoring network with mobile nodes providing informative observational data for the purpose of parameter estimation in DPSs. It was demonstrated that the relevant optimization task can be fit to the framework of an optimal control problem. Furthermore, a decentralized procedure for its solution is developed dedicated for scalable sensor networks. A proper reformulation of control problem in Mayer form makes it possible to directly apply numerous existing solvers for dynamic optimization. ACKNOWLEDGEMENTS The work was supported by the Polish Ministry of Science and Higher Education under Grant N N519 2971 33. REFERENCES Atkinson, A.C., Donev, A.N., and Tobias, R. (2007). Optimum Experimental Design, with SAS. Oxford University Press, Oxford. Bertsekas, D.P. and Tsitsiklis, J.N. (1997). Parallel and Distributed Computation: Numerical Methods. Athena Scientific, Belmont, MA. Boyd, S., Ghosh, A., Prabhakar, B., and Shah, D. (2006). Randomized gossip algorithms. IEEE Transactions on Information Theory, 52(6), 2508–2530. Bullo, F., Cort´ es, J., and Mart´ınez, S. (2009). Distributed Control of Robotic Networks. Applied Mathematics Series. Princeton University Press. To appear. Electronically available at http://coordinationbook.info. Cassandras, C.G. and Li, W. (2005). Sensor networks and cooperative control. European Journal of Control, 11(4–5), 436–463.
214
1st IFAC NECSYS (NECSYS'09) Venice, Italy, September 24-26, 2009
0.8
0.8
0.8
0.8
0.6
0.6
0.6
0.6
x2
1
x2
1
x2
1
x2
1
0.4
0.4
0.4
0.4
0.2
0.2
0.2
0.2
0 0
0.2
0.4
0.6
0.8
0 0
1
0.2
0.4
x
0.6
0.8
0 0
1
0.2
0.4
x
1
(a) t = 0.05
0.6
0.8
0 0
1
0.2
0.4
x
1
(b) t = 0.3
0.6
0.8
1
x
1
1
(c) t = 0.7
(d) t = 1
Fig. 1. Temporal changes in the wind velocity field and pollutant concentration (star indicates the location of source).
0.8
0.8
0.8
0.8
0.6
0.6
0.6
0.6
x2
1
x2
1
x2
1
x2
1
0.4
0.4
0.4
0.4
0.2
0.2
0.2
0.2
0 0
0.2
0.4
0.6
x
1
(a) k = 10
0.8
1
0 0
0.2
0.4
0.6
0.8
0 0
1
x
0.2
0.4
0.6
x
1
1
(b) k = 20
(c) k = 40
0.8
1
0 0
0.2
0.4
0.6
0.8
1
x
1
(d) Centralized approach
Fig. 2. Approximations of D-optimal sensor trajectories obtained in different stages of distributed configuration scheme (Algorithm 1) (a)-(c) versus trajectories determined in a centralized fashion (d). Open circles indicate the initial positions of mobile nodes. Jain, N. and Agrawal, D.P. (2005). Current trends in wireless sensor network design. International Journal of Distributed Sensor Networks, 1, 101–122. Jeremi´ c, A. and Nehorai, A. (2000). Landmine detection and localization using chemical sensor array processing. IEEE Transactions on Signal Processing, 48(5), 1295–1305. Kubrusly, C.S. and Malebranche, H. (1985). Sensors and controllers location in distributed systems — A survey. Automatica, 21(2), 117–128. Nehorai, A., Porat, B., and Paldi, E. (1995). Detection and localization of vapor-emitting sources. IEEE Transactions on Signal Processing, 43(1), 243–253. ¨ Ogren, P., Fiorelli, E., and Leonard, N.E. (2004). Cooperative control of mobile sensor networks: Adaptive gradient climbing in a distributed environment. IEEE Transactions on Automatic Control, 49(8), 1292–1302. Patan, M. (2004). Optimal Observation Strategies for Parameter Estimation of Distributed Systems. Zielona G´ ora: University Press. May be viewed at http://www.zbc.zgora.pl. Patan, M. (2006). Optimal activation policies for continous scanning observations in parameter estimation of distributed systems. International Journal of Systems Science, 37(11), 763–775. Patan, M. (2008). A parallel sensor scheduling technique for fault detection in distributed parameter systems. In Lecture Notes in Computer Science: Euro-Par 2008 : Parallel Procesing, volume 5168, 833–843. Springer-Verlag, Berlin Heidelberg. Patan, M., Chen, Y., and Tricaud, C. (2008). Resource-constrained sensor routing for parameter estimation of distributed systems. In Proc. 17th IFAC World Congress, Seoul, South Korea, 4–8 July 2008. Patan, M. and Patan, K. (2005). Optimal observation strategies for model-based fault detection in distributed systems. International Journal of Control, 78(18), 1497–1510. Patan, M. and Uci´ nski, D. (2005). Optimal activation strategy of discrete scanning sensors for fault detection in distributedparameter systems. In Proc. 16th IFAC World Congress, Prague,
Czech Republic, 4–8 July 2005. Published on CD-ROM. Patan, M. and Uci´ nski, D. (2008). Configuring a sensor network for fault detection in distributed parameter systems. International Journal of Applied Mathematics and Computer Science, 18(4), 513–524. Point, N., Vande Wouwer, A., and Remy, M. (1996). Practical issues in distributed parameter estimation: Gradient computation and optimal experiment design. Control Engineering Practice, 4(11), 1553–1562. Polak, E. (1997). Optimization. Algorithms and Consistent Approximations. Applied Mathematical Sciences. Springer-Verlag, New York. Quereshi, Z.H., Ng, T.S., and Goodwin, G.C. (1980). Optimum experimental design for identification of distributed parameter systems. International Journal of Control, 31(1), 21–29. Schwartz, A.L., Polak, E., and Chen, Y. (1997). A Matlab Toolbox for Solving Optimal Control Problems. Version 1.0 for Windows. URL http://www.schwartz-home.com/~adam/RIOTS/. Uci´ nski, D. (2000). Optimal selection of measurement locations for parameter estimation in distributed processes. International Journal of Applied Mathematics and Computer Science, 10(2), 357–379. Uci´ nski, D. (2005). Optimal Measurement Methods for DistributedParameter System Identification. CRC Press, Boca Raton, FL. Uci´ nski, D. and Chen, Y. (2005). Time-optimal path planning of moving sensors for parameter estimation of distributed systems. In Proc. 44th IEEE Conference on Decision and Control, and the European Control Conference 2005, Seville, Spain. Published on CD-ROM. ´ and Pronzato, L. (1997). Identification of Parametric Walter, E. Models from Experimental Data. Communications and Control Engineering. Springer-Verlag, Berlin. Zhao, F. and Guibas, L.J. (2004). Wireless Sensor Networks: An Information Processing Approach. Morgan Kaufmann Publishers, Amsterdam.
215