- Email: [email protected]
Constrained Mobile Sensor Routing for Parameter Estimation of Spatiotemporal Processes
16th IFAC Symposium on System Identification The International Federation of Automatic Control Brussels, Belgium. July 11-13, 2012
Constrained Mobile Sensor Routing for Parameter Estimation of Spatiotemporal Processes Dariusz Uciński, Maciej Patan ∗ ∗
Institute of Control and Computation Engineering University of Zielona Góra, ul. Podgórna 50, 65–246 Zielona Góra, Poland (e-mail: {D.Ucinski,M.Patan}@issi.uz.zgora.pl) Abstract: An approach is proposed to joint activation of stationary sensor nodes and design of mobile sensor trajectories in a hybrid sensor network collecting measurements for parameter estimation of a process described by a partial differential equation. The ultimate objective is maximization of the log-determinant of the information matrix associated with the estimated parameters. The search for the optimal solution is performed using the branch-and-bound method in which a block coordinate ascent method is employed to produce an upper bound to the maximum objective function. It alternates between solving a relaxed combinatorial problem for the selection of active stationary nodes and a relaxed optimal control problem for the design of sensor trajectories. The former is achieved using a simplicial decomposition algorithm in which the restricted master problem is solved using a multiplicative algorithm for optimal design. In turn, the latter is solved by another algorithm for optimal design, namely the Wynn-Fedorov algorithm, which is capable of finding an optimal element in the convex hull of the set of attainable information matrices and can be easily implemented by using a standard optimal control solver as its component. 1. INTRODUCTION Severe requirements imposed by modern process control are frequently associated with using very accurate models in which non-negligible spatial dynamics has to be included in addition to the temporal one. The processes in question are usually termed distributed parameter systems (DPSs) and they are described by partial differential equations (PDEs). A major difficulty in calibration of such models is the impossibility to measure process variables over the entire spatial domain. What is more, the observations are inexact by virtue of inherent measurement errors. This leads to the question of where to locate sensors so that the information content of the resulting signals be as high as possible. Over the past years, laborious research has been observed on the development of strategies for efficient sensor placement (for reviews, see papers by Kubrusly and Malebranche [1985], van de Wal and de Jager [2001] and El Jai and Hamzaoui [2009] or the comprehensive monograph by Uciński [2005]). The adopted optimization criteria are essentially the same, i.e., various scalar performance measures based on the Fisher Information Matrix (FIM) associated with the parameters to be identified are maximized, see Rafajłowicz [1981, 1983]. In recent years, this methodology has been substantially refined to extend its applicability. A comprehensive treatment of both theoretical and algorithmic aspects of the resulting sensor location strategies is contained in the monographs by Uciński [2005] and Sun [1994]. A revived interest in optimal sensor location is correlated with advances in Sensor Networks (SNs) [Zhao and 978-3-902823-06-9/12/$20.00 © 2012 IFAC
Guibas, 2004, Hirsch et al., 2008, Cassandras and Li, 2005, Bauer, 2008]. In this kind of network, each sensor node has a sensing capability, as well as limited energy supply, computing power, memory and communication ability. These inexpensive, low-power communication devices can be placed throughout the physical space, providing dense sensing close to physical phenomena. On one hand, SNs have recently come into prominence because they hold the potential to revolutionize observation systems. On the other hand, however, completely new design problems have to be solved. In this context, an appealing alternative to stationary sensor nodes is to apply mobile ones. The complexity of the resulting optimization problem is compensated by the clear advantage that sensors are not assigned to fixed positions which are optimal only on the average, but are capable of tracking points which provide best information about the parameters to be identified. Consequently, by actively reconfiguring a sensor system we can expect the minimal value of an adopted design criterion to be lower than the one for the stationary case. In his seminal article on experimental design for mobile observations, Rafajłowicz [1986] considers the D-optimality criterion and seeks an optimal time-dependent measure, rather than the trajectories themselves. On the other hand, Uciński [2005, 2000], apart from generalizations of Rafajłowicz’s results, develops some computational algorithms based on the FIM. He reduces the problem to a stateconstrained optimal-control one for which solutions are obtained via the methods of successive linearizations which is capable of handling various constraints imposed on sensor motions. Quite similar techniques were used in the
1317
10.3182/20120711-3-BE-2027.00363
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
work by Zhao and Nehorai [2006] to detect and estimate the chemical dispersion of possibly moving sources using model-based integrated sensor array processing. In the same setting, the Ds -optimality criterion was employed by Uciński and Patan [2010] to determine mobile sensor trajectories to detect and localize moving contamination sources. A receding-horizon scheme was then used to make this technique work on-line. In turn, the work by Uciński and Chen [2005] was intended as an attempt to properly formulate and solve the timeoptimal problem for moving sensors which observe the state of a DPS so as to estimate some of its parameters. Possibilities of exploiting parallel computations on clusters of PCs to reduce the time spent on determining optimal sensor node trajectories were discussed by Zięba and Uciński [2008]. More applications in the area of SNs can be found in the recent monograph by Song et al. [2009]. In this paper, we wish to demonstrate how the existing sensor placement techniques can be combined to solve the node activation problem for a hybrid SN which is supposed to collect data for parameter estimation of the observed DPS. The considered network consists of a number of mobile nodes which can move in a given spatial domain and, therefore, we would like their trajectories to be optimal in a sense. In addition to that, the data from mobile sensors are to be complemented by the ones gathered by a given number of nodes selected from among a greater number of nodes whose locations in space are fixed. Therefore, a decision must be made about which subset of stationary sensors is to be activated. The logdeterminant of the FIM serves as the optimality measure of this hybrid observation system. Mathematically, the problem is a mixed discrete optimal control one and, due to its potential high dimensionality, naive solutions are deemed to failure. In this contribution we deal with this problem by applying the branch-and-bound (BB) method to drasticallylly reduce the search space. The device to prune the BB tree is an original procedure producing upper bounds to the D-optimal design criterion. The key idea behind it is alternation between two relaxed problems, namely a discrete optimization one related to stationary sensors and an optimal control one associated with moving sensors.
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 `-dimensional vector of unknown constant parameters θ which must be estimated using observations of the system. The implicit dependence of the state y on the parameter vector θ will be be reflected by the notation y(x, t; θ). The parameter vector θ is to be restored based on measurements from both stationary and mobile sensors. 2.2 Sensor activation model A subnetwork of our SN is supposed to consist of n nodes activated from among a total of N available ones located at positions xi , i = 1, . . . , N . We assume that their observations are modelled as zis (t) = vi y(xi , t) + ε(xi , t) , t ∈ T, (4) i = 1, . . . , N, where ε(x, t) denotes the measurement disturbance which is customarily assumed to be a zeromean Gaussian white noise in both time and space, and the flag vi takes the value 1 or 0 depending on whether a sensor located at xi is or is not active, respectively. Obviously, the vector v = (v1 , . . . , vN ) must belong to the feasible set N X o n vi = n, vi = 0 or 1, ∀i ∈ I , (5) V = (v1 , . . . , vN ) i=1
where the symbol I denotes the index set {1, . . . , N } of possible sensor locations. The FIM which associated with the resulting measurements can then be written as N X M s (v) = vi Mis , (6) i=1
where Mis =
Z
g(xi , t)g T (xi , t) dt.
(7)
T
The quantity g(x, t) = ∇ϑ y(x, t; ϑ) ϑ=θ0
(8)
0
is the so-called sensitivity vector, θ being a prior estimate to the unknown parameter vector θ, see Uciński [2000, 2005].
2. SENSOR LOCATION PROBLEM 2.1 System description
2.3 Mobile sensor model d
Let Ω ⊂ R be a bounded spatial domain with sufficiently smooth boundary Γ, and let T = (0, tf ] be a bounded time interval. Consider a DPS whose scalar state at a spatial ¯ ⊂ Rd and time instant t ∈ T¯ is denoted by point x ∈ Ω y(x, t). Mathematically, the system state is governed by the PDE ∂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
Node dynamics. In addition to n stationary sensors, we exploit m mobile nodes. We assume that all sensors are conveyed by identical vehicles whose motions are described by x˙ j (t) = f (xj (t), uj (t)) a.e. on T , xj (0) = xj0 , (9) d r d where a given function f : R × R → R is required to be continuously differentiable, xj0 ∈ Rd defines an initial sensor configuration, and uj : T → Rr is a measurable control function which satisfies ul ≤ uj (t) ≤ uu a.e. on T (10) for some constant vectors ul and uu , j = 1, . . . , m.
1318
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
For each j = 1, . . . , m, given any initial position xj0 and any control function, there is a unique absolutely continuous function xj : T → Rd which satisfies (9) a.e. on T . In what follows, we will call it the state trajectory corresponding to xj0 and uj . The measurements along these trajectories are of the form zjm (t) = y(xj (t), t) + ε(xj (t), t), t ∈ T, (11) for j = 1, . . . , m. Pathwise state constraints. In practice, some restrictions on the motions are inevitably induced. First of all, all sensors should stay within an admissible region Ωad where measurements are allowed. We assume that it is a compact set defined as follows: Ωad = {x ∈ Ω ∪ Γ | bl (x) ≤ 0, l = 1, . . . , L} (12) where bl ’s are given continuously differentiable functions. Accordingly, the conditions bl (xj (t)) ≤ 0, ∀ t ∈ T (13) must be fulfilled, where 1 ≤ l ≤ L and 1 ≤ j ≤ m. Parameterization of vehicle controls. From now on we make the assumption that the controls of the available vehicles can be represented in the parametric form uj (t) = η(t, aj ), t ∈ T, (14) where η denotes a given function such that η( · , aj ) is continuous for each fixed aj and η(t, · ) is continuous for each fixed t, the constant parameter vector aj ranging over a compact set A ⊂ Rq . An exemplary parameterization can rely on using B-splines as employed in numerous optimal control solvers, e.g., RIOTS_95 Schwartz et al. [1997]. Based on this parameterization, we can define the mapping χ which assigns every cj = (xj0 , aj ) ∈ Ωad × A a trajectory xj = χ(cj ) through solving (9) for the initial position xj0 and control defined by (14). Since only the controls and trajectories satisfying the imposed constraints are interesting, we introduce the set m C = (c1 , . . . , cm ) ∈ A × Ωad η( · , aj ) satisifes (10), χ(cj ) satisfies (13), j = 1, . . . , m (15) and assume that it is nonempty. A trivial verification shows that C is also compact.
of linear and non-linear constraints, where the decision variables are both continuous (c) and discrete (v). At first sight, one might think of solving it by determining γ(u) = arg maxc∈C P(v, c) for each v ∈ V using optimal control techniques as described, e.g., in [Uciński, 2005, Ch. 4], and then setting v ? = arg maxv∈V P(v, γ(v)) and c? = γ(v ? ). Unfortunately, the inherent combinatorial nature here implies an extremely high dimensionality of the set V even for moderate values of N and n. Consequently, even if we have at our disposal an extremely fast solver for state-constrained optimal control problems, the exhaustive search may readily consume appreciable computer time and space. Therefore, our objective is to employ an algorithm which would use an implicit enumeration of V to determine an optimal design in an efficient manner. In what follows, we propose to solve it using the branch-andbound (BB) method. 3. BRANCH-AND-BOUND SOLUTION The BB constitutes a general algorithmic technique for finding optimal solutions of various optimization problems, especially discrete or combinatorial, see Floudas [2001], Bertsekas [1999]. If applied carefully, it can lead to algorithms that run reasonably fast on average. Our implementation of BB for Problem P involves the partition of the feasible set V into subsets. It is customary to select subsets of the form V (I0 , I1 ) = v ∈ V | vi = 0, ∀i ∈ I0 , vi = 1, ∀i ∈ I1 , (19) where I0 and I1 are disjoint subsets of I. Consequently, V (I0 , I1 ) is the subset of V such that a sensor is activated at the locations with indices in I1 , no sensor is activated at the locations with indices in I0 , and a sensor may or may not be activated at the remaining locations. Each subset V (I0 , I1 ) is identified with a node in the BB tree. The key assumption in the BB method is that for every nonterminal node V (I0 , I1 ), i.e., the node for which I0 ∪ I1 6= I, there is an algorithm that determines an ¯ 0 , I1 ) to the maximum design criterion upper bound P(I over V (I0 , I1 ), i.e., ¯ 0 , I1 ) ≥ P(I max P(v, c), (20) (v,c)∈V (I0 ,I1 )×C
Given m sensors, we thus obtain trajectories xj ( · ) corresponding to vectors c ∈ Rm(d+q) . The FIM corresponding to mobile nodesZcan then be written down as M m (c) = g(x(t), t)g T (x(t), t) dt. (16) j
and a feasible solution (v, c) ∈ V × C for which P(v, c) can serve as a lower bound to the maximum design criterion ¯ 0 , I1 ) by solving the over V . We propose to compute P(I following relaxed problem:
2.4 Ultimate optimal design problem
e where Problem R(I0 , I1 ): Find a pair (¯ v , c¯) ∈ Ve (I0 , I1 )× C, Ve (I0 , I1 )
T
x=χ(c )
The D-optimum experimental design problem for the overall SN takes the following form:
= co(V (I0 , I1 )) n = v ∈ RN vi = 0, ∀i ∈ I0 , vi = 1, ∀i ∈ I1 ,
Problem P: Find a pair (v ? , c? ) ∈ V × C to maximize P(v, c) = log det(M (v, c)), (17) where the total FIM is M (v, c) = M s (v) + M m (c). (18) This constitutes a mixed discrete non-linear programming problem which seeks a global optimum subject to a set
0 ≤ vi ≤ 1, ∀i ∈ I \ (I0 ∪ I1 ),
N X
(21)
o vi = n ,
i=1
e = co(C), C
(22)
e being the convex hulls of V (I0 , I1 ) and C, Ve (I0 , I1 ) and C respectively.
1319
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
In Problem R(I0 , I1 ) all 0–1 constraints on the variables vi are relaxed by allowing them to take any value in the interval [0, 1], except that the variables vi , i ∈ I0 ∪ I1 are fixed at either 0 or 1. Let us also note that c is now sought e and not in C. A simple and efficient in the convex hull C, method for its solution is given in the sequel. As a result ¯ 0 , I1 ) = P(¯ of its application, we set P(I v ). As for (v, c), we can specify it as the best feasible solution (i.e., an element of V × C) found so far. If no solution has been found yet, we can either set the lower bound to −∞, or use an initial guess about the optimal solution (experience provides evidence that the latter choice leads to much more rapid convergence). The result of solving Problem R(I0 , I1 ) can serve as a basis to construct a branching rule for the binary BB tree. We adopt here the approach in which the node/subset V (I0 , I1 ) is expanded (i.e., partitioned) by first picking out all fractional values from among the values of the relaxed variables, and then rounding to 0 and 1 a value which is the most distant from both 0 and 1. Specifically, we apply the following steps: (i) Determine i? = arg
min
|vi − 0.5|.
i∈I\(I0 ∪I1 )
(23)
(ii) Partition V (I0 , I1 ) into V (I0 ∪ {i? }, I1 ) and V (I0 , I1 ∪ {i? }) whereby two descendants of the node in question are defined. A recursive application of the branching rule starts from the root of the BB tree, which corresponds to the trivial subset V (∅, ∅) = V and the fully relaxed problem. Each node of the BB tree corresponds to a continuous relaxed problem, R(I0 , I1 ), while each edge corresponds to fixing one relaxed variable at 0 or 1. The above scheme has to be complemented with a search strategy to incrementally explore all the nodes of the BB tree. Here we use a common depth-first technique which always expands the deepest node in the current fringe of the search tree. The search then proceeds immediately to the deepest level of the search tree, where the nodes have no successors. In this way, lower bounds on the optimal solution can be found or improved as fast as possible. A recursive version of the resulting depth-first branchand-bound is implemented in Algorithm 1. The operators involved in this implementation are as follows: • Singularity-Test(I0 , I1 ) returns true only if expansion of the current node will result in a singular FIM. • Relaxed-Solution(I0 , I1 ) returns a solution to Problem R(I0 , I1 ). • Det-FIM(v) returns the log-determinant of the FIM corresponding to v. • Integral-Test(v) returns true only if the current solution v is integral. • Index-Branch(v) returns the index defined by (23). 4. BLOCK COORDINATE ASCENT TO DETERMINE RELAXED SOLUTIONS The main contribution of this article consists in setting forth an efficient scheme to implement the procedure
Algorithm 1 A recursive version of the depth-first branch-and-bound method. It uses two global variables, LOWER and v_best, which are respectively the maximal value of the FIM determinant over feasible solutions found so far and the solution at which it is attained. 1: procedure Recursive-DFBB(I0 , I1 ) 2: if |I0 | > N − n or |I1 | > n then 3: return 4: end if 5: if Singularity-Test(I0 , I1 ) then 6: return 7: end if 8: (v_relaxed, c_relaxed) 9: ← Relaxed-Solution(I0 , I1 ) 10: det_relaxed ← Det-FIM(v_relaxed, c_relaxed ) 11: if det_relaxed ≤ LOWER then 12: return 13: else if Integral-Test(v_relaxed ) then 14: v_best ← v_relaxed 15: LOWER ← det_relaxed 16: return 17: else 18: i? ← Index-Branch(v_relaxed ) 19: Recursive-DFBB(I0 ∪ i? , I 1 ) 20: Recursive-DFBB(I0 , I1 ∪ i? ) 21: end if 22: end procedure Relaxed-Solution. It is based on the block coordinate ascent method [Bertsekas, 1999, Ch. 2.7]. Specifically, the method proceeds iteratively, and at each iteration, P is maximized with respect to one variable (either v or c) while the other variable is held fixed. In detail, its scheme is embodied by Algorithm 2. Algorithm 2 Algorithm model for block coordinate ascent to obtain a relaxed solution. Step 0: (Initialization) Select v (0) ∈ Ve (I0 , I1 ) with 1 (0) vi = , i ∈ I \ (I0 ∪ I1 ), r where r = n − |I1 | and set k ← 0. Choose > 0 as the accuracy in finding the solution. Step 1: (Block ascent with respect to c) Find c(k+1) = arg max P(v (k) , c). e c∈C Step 2: (Block ascent with respect to v) Find v (k+1) = arg max P(v, c(k+1) ). e (I0 ,I1 ) v∈V Step 3: (Termination check) If |P(v (k+1) , c(k+1) ) − P(v (k) , c(k) )| < then STOP. Otherwise, set k ← k + 1 and go to Step 1.
The convergence of the block coordinate ascent method requires that P be strictly concave and continuously differentiable, and the feasible coordinate sets be convex. In our setting this is clearly satisfied, cf. [Uciński, 2005, Ch. 3].
1320
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
Due to space limitation, we are not able to characterize all technical details of Steps 2 and 3, but the following remarks make them rather clear:
1
• In order to implement Step 2, the well-known WynnFedorov algorithm is a very convenient tool due to its simplicity. It constitutes the most common method in optimum experimental design and its detailed description can be found in [Uciński, 2005, p. 51] or in Fedorov and Hackl [1997]. Its major computational effort is related to finding a candidate for a new support point in the current design. But here this can be achieved using any available numerical solver specialized in optimal optimal control, e.g., RIOTS_95, see Schwartz et al. [1997]. • Implementation of Step 3 proceeds basically on exactly the same lines as the implementation of the relaxed problem in Uciński and Patan [2007], Uciński [2009].
0.8
x2
0.6
0.4
0.2
0 0
0.2
0.4
0.6
0.8
1
x1
5. ILLUSTRATIVE EXAMPLE Consider a transport-chemistry process of an atmospheric pollutant over an urban area which is normalized to a unit square. At the point x0 = (0.3, 0.6) an active source is located which emits the pollutant being later spread over the area as a result of a combination of different processes such as turbulence, chemical decay and wind with the 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 observation interval T = [0, 1] can be described by the following advection-diffusion-reaction PDE: ∂y(x, t) +∇ · v(x, t)y(x, t) (24) ∂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 Ω, (25) ∂n 2 where the term η(x) = e−100kx−x0 k models the pollution source and ∂y/∂n stands for the derivative of y in the direction of the outward normal to the boundary Γ. The following form of the turbulent diffusivity coefficient is applied: κ(x) = θ1 + θ2 x21 + θ3 x22 , (26) and parameters θ = (θ1 , θ2 , θ3 ) are to be estimated based on measurement data. Our goal is to determine a D-optimal sensor configuration for the hybrid SN consisting of m = 1 mobile and n = 5 stationary sensor nodes selected from among a total of 121 ones residing at nodes of the even 11 × 11 spatial grid of Fig. 2. A Matlab program implementing the recursive version of the BB procedure was written using a PC equipped with Centrino T9300 processor (2.5GHz, 3 GB RAM) running Windows Vista and Matlab 2009b. The system of PDEs was solved using efficient solvers of COMSOL 3.5 based on the finite element method (FEM). The nominal values of the systems parameters were assumed to be θ10 = 0.02 and θ20 = θ30 = 0.005. Calculations were performed for a FEM mesh composed of 1800 triangles and 961 nodes and an evenly partitioned time interval (30 subintervals). The sensitivity coefficients
Fig. 2. D-optimal sensor configuration. Discs and open circles indicate the potential locations and activated nodes, respectively. A square and the solid line represent the initial position and path of mobile node, respectively were then linearly interpolated and stored in memory. In order to solve the optimal control problems related to the design of mobile sensor path, the package RIOTS_95 was applied. As regard sensor dynamics, the following simple model was adopted: x˙ 1 (t) = u1 (t), x1 (0) = x10 , and the bounds for u1 were taken as |u1k (t)| ≤ 0.5, ∀t ∈ T, k = 1, 2. (27) After 27 recursive steps of BB search, which took about 150 minutes of computational time, the algorithm produced the D-optimal solution presented presented in Fig. 2. From the point of view of the estimation accuracy, the mobile sensor is expected to follow the path following greatest changes in the pollutant concentration, but due to the limited velocity, it is not able to reach the boundary of Ω. Therefore, the stationary sensors are located in those positions which are not available to mobile sensor and allow to gather additional data which complement the information about the distributed process. 6. CONCLUSIONS This work was intended as an attempt to solve a D-optimal design problem for a hybrid sensor network consisting of two subnetworks, namely those including stationary and mobile nodes. At the current stage, the proposed BB technique is rather suited for situations when the number of mobile sensors is rather low, but the main reason is the high cost incurred by numerically solving optimal control problems. This challenge leaves room for further refinements and improvements. Note, however, that the BB technique can be easily parallelized due to numerous decompositions applied here. This task constitutes a primary topic for our further research.
1321
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
0.8
0.8
0.8
0.8
0.6
0.6
0.6
0.6
x
x
2
1
2
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
1
0 0
x1
0.2
0.4
0.6
0.8
0 0
1
x1
(a) t = 0.05
0.2
0.4
0.6
x1
(b) t = 0.3
(c) t = 0.7
0.8
1
0 0
0.2
0.4
0.6
0.8
1
x1
(d) t = 1
Fig. 1. Temporal changes in the wind velocity field and pollutant concentration (a star indicates the location of the source). REFERENCES Peter H. Bauer. New challenges in dynamical systems: The networked case. International Journal of Applied Mathematics and Computer Science, 18(3):271–277, 2008. Dimitri P. Bertsekas. Nonlinear Programming. Optimization and Computation Series. Athena Scientific, Belmont, MA, 2nd edition, 1999. Christos G. Cassandras and Wei Li. Sensor networks and cooperative control. European Journal of Control, 11 (4–5):436–463, 2005. Abdelhaq El Jai and Houria Hamzaoui. Regional observation and sensors. International Journal of Applied Mathematics and Computer Science, 19(1):5–14, 2009. Valerii V. Fedorov and Peter Hackl. Model-Oriented Design of Experiments. Lecture Notes in Statistics. Springer-Verlag, New York, 1997. Christodoulos A. Floudas. Mixed integer nonlinear programming, MINLP. In C. A. Floudas and P. M. Pardalos, editors, Encyclopedia of Optimization, volume 3, pages 401–414. Kluwer Academic Publishers, Dordrecht, The Netherlands, 2001. Michael J. Hirsch, Panos M. Pardalos, Robert Murphey, and Don Grundel, editors. Advances in Cooperative Control and Optimization. Proceedings of the 7th International Conference on Cooperative Control and Optimization. Springer-Verlag, Berlin, 2008. C. S. Kubrusly and H. Malebranche. Sensors and controllers location in distributed systems — A survey. Automatica, 21(2):117–128, 1985. Ewaryst Rafajłowicz. Design of experiments for eigenvalue identification in distributed-parameter systems. International Journal of Control, 34(6):1079–1094, 1981. Ewaryst Rafajłowicz. Optimal experiment design for identification of linear distributed-parameter systems: Frequency domain approach. IEEE Transactions on Automatic Control, 28(7):806–808, July 1983. Ewaryst Rafajłowicz. Optimum choice of moving sensor trajectories for distributed parameter system identification. International Journal of Control, 43(5):1441–1451, 1986. Adam Lowell Schwartz, Elijah Polak, and YangQuan Chen. A Matlab Toolbox for Solving Optimal Control Problems. Version 1.0 for Windows, May 1997. URL http://www.schwartz-home.com/˜adam/RIOTS/. Zhen Song, YangQuan Chen, Chellury Ram Sastry, and Nazif Cihan Tas. Optimal Observation for Cyberphysical Systems: A Fisher-information-matrix-based
Approach. Springer-Verlag, London, 2009. Ne-Zheng Sun. Inverse Problems in Groundwater Modeling. Theory and Applications of Transport in Porous Media. Kluwer Academic Publishers, Dordrecht, The Netherlands, 1994. Dariusz Uciński. Optimal sensor location for parameter estimation of distributed processes. International Journal of Control, 73(13):1235–1248, 2000. Dariusz Uciński. D-optimum sensor activity scheduling for distributed parameter systems. In Preprints of the 15th IFAC Symposium on System Identification, Saint-Malo, France, July 6-8,, 2009. Published on CD-ROM. Dariusz Uciński. Optimal Measurement Methods for Distributed-Parameter System Identification. CRC Press, Boca Raton, FL, 2005. Dariusz Uciński and YangQuan Chen. 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, 2005. Published on CDROM. Dariusz Uciński and Maciej Patan. D-optimal design of a monitoring network for parameter estimation of distributed systems. Journal of Global Optimization, 39:291–322, 2007. Dariusz Uciński and Maciej Patan. Sensor network design for the estimation of spatially distributed processes. International Journal of Applied Mathematics and Computer Science, 20(3):459–481, 2010. Marc van de Wal and Bram de Jager. A review of methods for input/output selection. Automatica, 37:487–510, 2001. Feng Zhao and Leonidas J. Guibas. Wireless Sensor Networks: An Information Processing Approach. Morgan Kaufmann Publishers, Amsterdam, 2004. Tong Zhao and Arye Nehorai. Detecting and estimating biochemical dispersion of a moving source in a semiinfinite medium. IEEE Transactions on Signal Processing, 54(6):2213–2225, June 2006. Tomasz Zięba and Dariusz Uciński. Mobile sensor routing for parameter estimation of distributed systems using the parallel tunneling method. International Journal of Applied Mathematics and Computer Science, 18(3): 307–318, 2008.
1322
Constrained Mobile Sensor Routing for Parameter Estimation of Spatiotemporal Processes Dariusz Uciński, Maciej Patan ∗ ∗
Institute of Control and Computation Engineering University of Zielona Góra, ul. Podgórna 50, 65–246 Zielona Góra, Poland (e-mail: {D.Ucinski,M.Patan}@issi.uz.zgora.pl) Abstract: An approach is proposed to joint activation of stationary sensor nodes and design of mobile sensor trajectories in a hybrid sensor network collecting measurements for parameter estimation of a process described by a partial differential equation. The ultimate objective is maximization of the log-determinant of the information matrix associated with the estimated parameters. The search for the optimal solution is performed using the branch-and-bound method in which a block coordinate ascent method is employed to produce an upper bound to the maximum objective function. It alternates between solving a relaxed combinatorial problem for the selection of active stationary nodes and a relaxed optimal control problem for the design of sensor trajectories. The former is achieved using a simplicial decomposition algorithm in which the restricted master problem is solved using a multiplicative algorithm for optimal design. In turn, the latter is solved by another algorithm for optimal design, namely the Wynn-Fedorov algorithm, which is capable of finding an optimal element in the convex hull of the set of attainable information matrices and can be easily implemented by using a standard optimal control solver as its component. 1. INTRODUCTION Severe requirements imposed by modern process control are frequently associated with using very accurate models in which non-negligible spatial dynamics has to be included in addition to the temporal one. The processes in question are usually termed distributed parameter systems (DPSs) and they are described by partial differential equations (PDEs). A major difficulty in calibration of such models is the impossibility to measure process variables over the entire spatial domain. What is more, the observations are inexact by virtue of inherent measurement errors. This leads to the question of where to locate sensors so that the information content of the resulting signals be as high as possible. Over the past years, laborious research has been observed on the development of strategies for efficient sensor placement (for reviews, see papers by Kubrusly and Malebranche [1985], van de Wal and de Jager [2001] and El Jai and Hamzaoui [2009] or the comprehensive monograph by Uciński [2005]). The adopted optimization criteria are essentially the same, i.e., various scalar performance measures based on the Fisher Information Matrix (FIM) associated with the parameters to be identified are maximized, see Rafajłowicz [1981, 1983]. In recent years, this methodology has been substantially refined to extend its applicability. A comprehensive treatment of both theoretical and algorithmic aspects of the resulting sensor location strategies is contained in the monographs by Uciński [2005] and Sun [1994]. A revived interest in optimal sensor location is correlated with advances in Sensor Networks (SNs) [Zhao and 978-3-902823-06-9/12/$20.00 © 2012 IFAC
Guibas, 2004, Hirsch et al., 2008, Cassandras and Li, 2005, Bauer, 2008]. In this kind of network, each sensor node has a sensing capability, as well as limited energy supply, computing power, memory and communication ability. These inexpensive, low-power communication devices can be placed throughout the physical space, providing dense sensing close to physical phenomena. On one hand, SNs have recently come into prominence because they hold the potential to revolutionize observation systems. On the other hand, however, completely new design problems have to be solved. In this context, an appealing alternative to stationary sensor nodes is to apply mobile ones. The complexity of the resulting optimization problem is compensated by the clear advantage that sensors are not assigned to fixed positions which are optimal only on the average, but are capable of tracking points which provide best information about the parameters to be identified. Consequently, by actively reconfiguring a sensor system we can expect the minimal value of an adopted design criterion to be lower than the one for the stationary case. In his seminal article on experimental design for mobile observations, Rafajłowicz [1986] considers the D-optimality criterion and seeks an optimal time-dependent measure, rather than the trajectories themselves. On the other hand, Uciński [2005, 2000], apart from generalizations of Rafajłowicz’s results, develops some computational algorithms based on the FIM. He reduces the problem to a stateconstrained optimal-control one for which solutions are obtained via the methods of successive linearizations which is capable of handling various constraints imposed on sensor motions. Quite similar techniques were used in the
1317
10.3182/20120711-3-BE-2027.00363
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
work by Zhao and Nehorai [2006] to detect and estimate the chemical dispersion of possibly moving sources using model-based integrated sensor array processing. In the same setting, the Ds -optimality criterion was employed by Uciński and Patan [2010] to determine mobile sensor trajectories to detect and localize moving contamination sources. A receding-horizon scheme was then used to make this technique work on-line. In turn, the work by Uciński and Chen [2005] was intended as an attempt to properly formulate and solve the timeoptimal problem for moving sensors which observe the state of a DPS so as to estimate some of its parameters. Possibilities of exploiting parallel computations on clusters of PCs to reduce the time spent on determining optimal sensor node trajectories were discussed by Zięba and Uciński [2008]. More applications in the area of SNs can be found in the recent monograph by Song et al. [2009]. In this paper, we wish to demonstrate how the existing sensor placement techniques can be combined to solve the node activation problem for a hybrid SN which is supposed to collect data for parameter estimation of the observed DPS. The considered network consists of a number of mobile nodes which can move in a given spatial domain and, therefore, we would like their trajectories to be optimal in a sense. In addition to that, the data from mobile sensors are to be complemented by the ones gathered by a given number of nodes selected from among a greater number of nodes whose locations in space are fixed. Therefore, a decision must be made about which subset of stationary sensors is to be activated. The logdeterminant of the FIM serves as the optimality measure of this hybrid observation system. Mathematically, the problem is a mixed discrete optimal control one and, due to its potential high dimensionality, naive solutions are deemed to failure. In this contribution we deal with this problem by applying the branch-and-bound (BB) method to drasticallylly reduce the search space. The device to prune the BB tree is an original procedure producing upper bounds to the D-optimal design criterion. The key idea behind it is alternation between two relaxed problems, namely a discrete optimization one related to stationary sensors and an optimal control one associated with moving sensors.
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 `-dimensional vector of unknown constant parameters θ which must be estimated using observations of the system. The implicit dependence of the state y on the parameter vector θ will be be reflected by the notation y(x, t; θ). The parameter vector θ is to be restored based on measurements from both stationary and mobile sensors. 2.2 Sensor activation model A subnetwork of our SN is supposed to consist of n nodes activated from among a total of N available ones located at positions xi , i = 1, . . . , N . We assume that their observations are modelled as zis (t) = vi y(xi , t) + ε(xi , t) , t ∈ T, (4) i = 1, . . . , N, where ε(x, t) denotes the measurement disturbance which is customarily assumed to be a zeromean Gaussian white noise in both time and space, and the flag vi takes the value 1 or 0 depending on whether a sensor located at xi is or is not active, respectively. Obviously, the vector v = (v1 , . . . , vN ) must belong to the feasible set N X o n vi = n, vi = 0 or 1, ∀i ∈ I , (5) V = (v1 , . . . , vN ) i=1
where the symbol I denotes the index set {1, . . . , N } of possible sensor locations. The FIM which associated with the resulting measurements can then be written as N X M s (v) = vi Mis , (6) i=1
where Mis =
Z
g(xi , t)g T (xi , t) dt.
(7)
T
The quantity g(x, t) = ∇ϑ y(x, t; ϑ) ϑ=θ0
(8)
0
is the so-called sensitivity vector, θ being a prior estimate to the unknown parameter vector θ, see Uciński [2000, 2005].
2. SENSOR LOCATION PROBLEM 2.1 System description
2.3 Mobile sensor model d
Let Ω ⊂ R be a bounded spatial domain with sufficiently smooth boundary Γ, and let T = (0, tf ] be a bounded time interval. Consider a DPS whose scalar state at a spatial ¯ ⊂ Rd and time instant t ∈ T¯ is denoted by point x ∈ Ω y(x, t). Mathematically, the system state is governed by the PDE ∂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
Node dynamics. In addition to n stationary sensors, we exploit m mobile nodes. We assume that all sensors are conveyed by identical vehicles whose motions are described by x˙ j (t) = f (xj (t), uj (t)) a.e. on T , xj (0) = xj0 , (9) d r d where a given function f : R × R → R is required to be continuously differentiable, xj0 ∈ Rd defines an initial sensor configuration, and uj : T → Rr is a measurable control function which satisfies ul ≤ uj (t) ≤ uu a.e. on T (10) for some constant vectors ul and uu , j = 1, . . . , m.
1318
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
For each j = 1, . . . , m, given any initial position xj0 and any control function, there is a unique absolutely continuous function xj : T → Rd which satisfies (9) a.e. on T . In what follows, we will call it the state trajectory corresponding to xj0 and uj . The measurements along these trajectories are of the form zjm (t) = y(xj (t), t) + ε(xj (t), t), t ∈ T, (11) for j = 1, . . . , m. Pathwise state constraints. In practice, some restrictions on the motions are inevitably induced. First of all, all sensors should stay within an admissible region Ωad where measurements are allowed. We assume that it is a compact set defined as follows: Ωad = {x ∈ Ω ∪ Γ | bl (x) ≤ 0, l = 1, . . . , L} (12) where bl ’s are given continuously differentiable functions. Accordingly, the conditions bl (xj (t)) ≤ 0, ∀ t ∈ T (13) must be fulfilled, where 1 ≤ l ≤ L and 1 ≤ j ≤ m. Parameterization of vehicle controls. From now on we make the assumption that the controls of the available vehicles can be represented in the parametric form uj (t) = η(t, aj ), t ∈ T, (14) where η denotes a given function such that η( · , aj ) is continuous for each fixed aj and η(t, · ) is continuous for each fixed t, the constant parameter vector aj ranging over a compact set A ⊂ Rq . An exemplary parameterization can rely on using B-splines as employed in numerous optimal control solvers, e.g., RIOTS_95 Schwartz et al. [1997]. Based on this parameterization, we can define the mapping χ which assigns every cj = (xj0 , aj ) ∈ Ωad × A a trajectory xj = χ(cj ) through solving (9) for the initial position xj0 and control defined by (14). Since only the controls and trajectories satisfying the imposed constraints are interesting, we introduce the set m C = (c1 , . . . , cm ) ∈ A × Ωad η( · , aj ) satisifes (10), χ(cj ) satisfies (13), j = 1, . . . , m (15) and assume that it is nonempty. A trivial verification shows that C is also compact.
of linear and non-linear constraints, where the decision variables are both continuous (c) and discrete (v). At first sight, one might think of solving it by determining γ(u) = arg maxc∈C P(v, c) for each v ∈ V using optimal control techniques as described, e.g., in [Uciński, 2005, Ch. 4], and then setting v ? = arg maxv∈V P(v, γ(v)) and c? = γ(v ? ). Unfortunately, the inherent combinatorial nature here implies an extremely high dimensionality of the set V even for moderate values of N and n. Consequently, even if we have at our disposal an extremely fast solver for state-constrained optimal control problems, the exhaustive search may readily consume appreciable computer time and space. Therefore, our objective is to employ an algorithm which would use an implicit enumeration of V to determine an optimal design in an efficient manner. In what follows, we propose to solve it using the branch-andbound (BB) method. 3. BRANCH-AND-BOUND SOLUTION The BB constitutes a general algorithmic technique for finding optimal solutions of various optimization problems, especially discrete or combinatorial, see Floudas [2001], Bertsekas [1999]. If applied carefully, it can lead to algorithms that run reasonably fast on average. Our implementation of BB for Problem P involves the partition of the feasible set V into subsets. It is customary to select subsets of the form V (I0 , I1 ) = v ∈ V | vi = 0, ∀i ∈ I0 , vi = 1, ∀i ∈ I1 , (19) where I0 and I1 are disjoint subsets of I. Consequently, V (I0 , I1 ) is the subset of V such that a sensor is activated at the locations with indices in I1 , no sensor is activated at the locations with indices in I0 , and a sensor may or may not be activated at the remaining locations. Each subset V (I0 , I1 ) is identified with a node in the BB tree. The key assumption in the BB method is that for every nonterminal node V (I0 , I1 ), i.e., the node for which I0 ∪ I1 6= I, there is an algorithm that determines an ¯ 0 , I1 ) to the maximum design criterion upper bound P(I over V (I0 , I1 ), i.e., ¯ 0 , I1 ) ≥ P(I max P(v, c), (20) (v,c)∈V (I0 ,I1 )×C
Given m sensors, we thus obtain trajectories xj ( · ) corresponding to vectors c ∈ Rm(d+q) . The FIM corresponding to mobile nodesZcan then be written down as M m (c) = g(x(t), t)g T (x(t), t) dt. (16) j
and a feasible solution (v, c) ∈ V × C for which P(v, c) can serve as a lower bound to the maximum design criterion ¯ 0 , I1 ) by solving the over V . We propose to compute P(I following relaxed problem:
2.4 Ultimate optimal design problem
e where Problem R(I0 , I1 ): Find a pair (¯ v , c¯) ∈ Ve (I0 , I1 )× C, Ve (I0 , I1 )
T
x=χ(c )
The D-optimum experimental design problem for the overall SN takes the following form:
= co(V (I0 , I1 )) n = v ∈ RN vi = 0, ∀i ∈ I0 , vi = 1, ∀i ∈ I1 ,
Problem P: Find a pair (v ? , c? ) ∈ V × C to maximize P(v, c) = log det(M (v, c)), (17) where the total FIM is M (v, c) = M s (v) + M m (c). (18) This constitutes a mixed discrete non-linear programming problem which seeks a global optimum subject to a set
0 ≤ vi ≤ 1, ∀i ∈ I \ (I0 ∪ I1 ),
N X
(21)
o vi = n ,
i=1
e = co(C), C
(22)
e being the convex hulls of V (I0 , I1 ) and C, Ve (I0 , I1 ) and C respectively.
1319
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
In Problem R(I0 , I1 ) all 0–1 constraints on the variables vi are relaxed by allowing them to take any value in the interval [0, 1], except that the variables vi , i ∈ I0 ∪ I1 are fixed at either 0 or 1. Let us also note that c is now sought e and not in C. A simple and efficient in the convex hull C, method for its solution is given in the sequel. As a result ¯ 0 , I1 ) = P(¯ of its application, we set P(I v ). As for (v, c), we can specify it as the best feasible solution (i.e., an element of V × C) found so far. If no solution has been found yet, we can either set the lower bound to −∞, or use an initial guess about the optimal solution (experience provides evidence that the latter choice leads to much more rapid convergence). The result of solving Problem R(I0 , I1 ) can serve as a basis to construct a branching rule for the binary BB tree. We adopt here the approach in which the node/subset V (I0 , I1 ) is expanded (i.e., partitioned) by first picking out all fractional values from among the values of the relaxed variables, and then rounding to 0 and 1 a value which is the most distant from both 0 and 1. Specifically, we apply the following steps: (i) Determine i? = arg
min
|vi − 0.5|.
i∈I\(I0 ∪I1 )
(23)
(ii) Partition V (I0 , I1 ) into V (I0 ∪ {i? }, I1 ) and V (I0 , I1 ∪ {i? }) whereby two descendants of the node in question are defined. A recursive application of the branching rule starts from the root of the BB tree, which corresponds to the trivial subset V (∅, ∅) = V and the fully relaxed problem. Each node of the BB tree corresponds to a continuous relaxed problem, R(I0 , I1 ), while each edge corresponds to fixing one relaxed variable at 0 or 1. The above scheme has to be complemented with a search strategy to incrementally explore all the nodes of the BB tree. Here we use a common depth-first technique which always expands the deepest node in the current fringe of the search tree. The search then proceeds immediately to the deepest level of the search tree, where the nodes have no successors. In this way, lower bounds on the optimal solution can be found or improved as fast as possible. A recursive version of the resulting depth-first branchand-bound is implemented in Algorithm 1. The operators involved in this implementation are as follows: • Singularity-Test(I0 , I1 ) returns true only if expansion of the current node will result in a singular FIM. • Relaxed-Solution(I0 , I1 ) returns a solution to Problem R(I0 , I1 ). • Det-FIM(v) returns the log-determinant of the FIM corresponding to v. • Integral-Test(v) returns true only if the current solution v is integral. • Index-Branch(v) returns the index defined by (23). 4. BLOCK COORDINATE ASCENT TO DETERMINE RELAXED SOLUTIONS The main contribution of this article consists in setting forth an efficient scheme to implement the procedure
Algorithm 1 A recursive version of the depth-first branch-and-bound method. It uses two global variables, LOWER and v_best, which are respectively the maximal value of the FIM determinant over feasible solutions found so far and the solution at which it is attained. 1: procedure Recursive-DFBB(I0 , I1 ) 2: if |I0 | > N − n or |I1 | > n then 3: return 4: end if 5: if Singularity-Test(I0 , I1 ) then 6: return 7: end if 8: (v_relaxed, c_relaxed) 9: ← Relaxed-Solution(I0 , I1 ) 10: det_relaxed ← Det-FIM(v_relaxed, c_relaxed ) 11: if det_relaxed ≤ LOWER then 12: return 13: else if Integral-Test(v_relaxed ) then 14: v_best ← v_relaxed 15: LOWER ← det_relaxed 16: return 17: else 18: i? ← Index-Branch(v_relaxed ) 19: Recursive-DFBB(I0 ∪ i? , I 1 ) 20: Recursive-DFBB(I0 , I1 ∪ i? ) 21: end if 22: end procedure Relaxed-Solution. It is based on the block coordinate ascent method [Bertsekas, 1999, Ch. 2.7]. Specifically, the method proceeds iteratively, and at each iteration, P is maximized with respect to one variable (either v or c) while the other variable is held fixed. In detail, its scheme is embodied by Algorithm 2. Algorithm 2 Algorithm model for block coordinate ascent to obtain a relaxed solution. Step 0: (Initialization) Select v (0) ∈ Ve (I0 , I1 ) with 1 (0) vi = , i ∈ I \ (I0 ∪ I1 ), r where r = n − |I1 | and set k ← 0. Choose > 0 as the accuracy in finding the solution. Step 1: (Block ascent with respect to c) Find c(k+1) = arg max P(v (k) , c). e c∈C Step 2: (Block ascent with respect to v) Find v (k+1) = arg max P(v, c(k+1) ). e (I0 ,I1 ) v∈V Step 3: (Termination check) If |P(v (k+1) , c(k+1) ) − P(v (k) , c(k) )| < then STOP. Otherwise, set k ← k + 1 and go to Step 1.
The convergence of the block coordinate ascent method requires that P be strictly concave and continuously differentiable, and the feasible coordinate sets be convex. In our setting this is clearly satisfied, cf. [Uciński, 2005, Ch. 3].
1320
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
Due to space limitation, we are not able to characterize all technical details of Steps 2 and 3, but the following remarks make them rather clear:
1
• In order to implement Step 2, the well-known WynnFedorov algorithm is a very convenient tool due to its simplicity. It constitutes the most common method in optimum experimental design and its detailed description can be found in [Uciński, 2005, p. 51] or in Fedorov and Hackl [1997]. Its major computational effort is related to finding a candidate for a new support point in the current design. But here this can be achieved using any available numerical solver specialized in optimal optimal control, e.g., RIOTS_95, see Schwartz et al. [1997]. • Implementation of Step 3 proceeds basically on exactly the same lines as the implementation of the relaxed problem in Uciński and Patan [2007], Uciński [2009].
0.8
x2
0.6
0.4
0.2
0 0
0.2
0.4
0.6
0.8
1
x1
5. ILLUSTRATIVE EXAMPLE Consider a transport-chemistry process of an atmospheric pollutant over an urban area which is normalized to a unit square. At the point x0 = (0.3, 0.6) an active source is located which emits the pollutant being later spread over the area as a result of a combination of different processes such as turbulence, chemical decay and wind with the 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 observation interval T = [0, 1] can be described by the following advection-diffusion-reaction PDE: ∂y(x, t) +∇ · v(x, t)y(x, t) (24) ∂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 Ω, (25) ∂n 2 where the term η(x) = e−100kx−x0 k models the pollution source and ∂y/∂n stands for the derivative of y in the direction of the outward normal to the boundary Γ. The following form of the turbulent diffusivity coefficient is applied: κ(x) = θ1 + θ2 x21 + θ3 x22 , (26) and parameters θ = (θ1 , θ2 , θ3 ) are to be estimated based on measurement data. Our goal is to determine a D-optimal sensor configuration for the hybrid SN consisting of m = 1 mobile and n = 5 stationary sensor nodes selected from among a total of 121 ones residing at nodes of the even 11 × 11 spatial grid of Fig. 2. A Matlab program implementing the recursive version of the BB procedure was written using a PC equipped with Centrino T9300 processor (2.5GHz, 3 GB RAM) running Windows Vista and Matlab 2009b. The system of PDEs was solved using efficient solvers of COMSOL 3.5 based on the finite element method (FEM). The nominal values of the systems parameters were assumed to be θ10 = 0.02 and θ20 = θ30 = 0.005. Calculations were performed for a FEM mesh composed of 1800 triangles and 961 nodes and an evenly partitioned time interval (30 subintervals). The sensitivity coefficients
Fig. 2. D-optimal sensor configuration. Discs and open circles indicate the potential locations and activated nodes, respectively. A square and the solid line represent the initial position and path of mobile node, respectively were then linearly interpolated and stored in memory. In order to solve the optimal control problems related to the design of mobile sensor path, the package RIOTS_95 was applied. As regard sensor dynamics, the following simple model was adopted: x˙ 1 (t) = u1 (t), x1 (0) = x10 , and the bounds for u1 were taken as |u1k (t)| ≤ 0.5, ∀t ∈ T, k = 1, 2. (27) After 27 recursive steps of BB search, which took about 150 minutes of computational time, the algorithm produced the D-optimal solution presented presented in Fig. 2. From the point of view of the estimation accuracy, the mobile sensor is expected to follow the path following greatest changes in the pollutant concentration, but due to the limited velocity, it is not able to reach the boundary of Ω. Therefore, the stationary sensors are located in those positions which are not available to mobile sensor and allow to gather additional data which complement the information about the distributed process. 6. CONCLUSIONS This work was intended as an attempt to solve a D-optimal design problem for a hybrid sensor network consisting of two subnetworks, namely those including stationary and mobile nodes. At the current stage, the proposed BB technique is rather suited for situations when the number of mobile sensors is rather low, but the main reason is the high cost incurred by numerically solving optimal control problems. This challenge leaves room for further refinements and improvements. Note, however, that the BB technique can be easily parallelized due to numerous decompositions applied here. This task constitutes a primary topic for our further research.
1321
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
0.8
0.8
0.8
0.8
0.6
0.6
0.6
0.6
x
x
2
1
2
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
1
0 0
x1
0.2
0.4
0.6
0.8
0 0
1
x1
(a) t = 0.05
0.2
0.4
0.6
x1
(b) t = 0.3
(c) t = 0.7
0.8
1
0 0
0.2
0.4
0.6
0.8
1
x1
(d) t = 1
Fig. 1. Temporal changes in the wind velocity field and pollutant concentration (a star indicates the location of the source). REFERENCES Peter H. Bauer. New challenges in dynamical systems: The networked case. International Journal of Applied Mathematics and Computer Science, 18(3):271–277, 2008. Dimitri P. Bertsekas. Nonlinear Programming. Optimization and Computation Series. Athena Scientific, Belmont, MA, 2nd edition, 1999. Christos G. Cassandras and Wei Li. Sensor networks and cooperative control. European Journal of Control, 11 (4–5):436–463, 2005. Abdelhaq El Jai and Houria Hamzaoui. Regional observation and sensors. International Journal of Applied Mathematics and Computer Science, 19(1):5–14, 2009. Valerii V. Fedorov and Peter Hackl. Model-Oriented Design of Experiments. Lecture Notes in Statistics. Springer-Verlag, New York, 1997. Christodoulos A. Floudas. Mixed integer nonlinear programming, MINLP. In C. A. Floudas and P. M. Pardalos, editors, Encyclopedia of Optimization, volume 3, pages 401–414. Kluwer Academic Publishers, Dordrecht, The Netherlands, 2001. Michael J. Hirsch, Panos M. Pardalos, Robert Murphey, and Don Grundel, editors. Advances in Cooperative Control and Optimization. Proceedings of the 7th International Conference on Cooperative Control and Optimization. Springer-Verlag, Berlin, 2008. C. S. Kubrusly and H. Malebranche. Sensors and controllers location in distributed systems — A survey. Automatica, 21(2):117–128, 1985. Ewaryst Rafajłowicz. Design of experiments for eigenvalue identification in distributed-parameter systems. International Journal of Control, 34(6):1079–1094, 1981. Ewaryst Rafajłowicz. Optimal experiment design for identification of linear distributed-parameter systems: Frequency domain approach. IEEE Transactions on Automatic Control, 28(7):806–808, July 1983. Ewaryst Rafajłowicz. Optimum choice of moving sensor trajectories for distributed parameter system identification. International Journal of Control, 43(5):1441–1451, 1986. Adam Lowell Schwartz, Elijah Polak, and YangQuan Chen. A Matlab Toolbox for Solving Optimal Control Problems. Version 1.0 for Windows, May 1997. URL http://www.schwartz-home.com/˜adam/RIOTS/. Zhen Song, YangQuan Chen, Chellury Ram Sastry, and Nazif Cihan Tas. Optimal Observation for Cyberphysical Systems: A Fisher-information-matrix-based
Approach. Springer-Verlag, London, 2009. Ne-Zheng Sun. Inverse Problems in Groundwater Modeling. Theory and Applications of Transport in Porous Media. Kluwer Academic Publishers, Dordrecht, The Netherlands, 1994. Dariusz Uciński. Optimal sensor location for parameter estimation of distributed processes. International Journal of Control, 73(13):1235–1248, 2000. Dariusz Uciński. D-optimum sensor activity scheduling for distributed parameter systems. In Preprints of the 15th IFAC Symposium on System Identification, Saint-Malo, France, July 6-8,, 2009. Published on CD-ROM. Dariusz Uciński. Optimal Measurement Methods for Distributed-Parameter System Identification. CRC Press, Boca Raton, FL, 2005. Dariusz Uciński and YangQuan Chen. 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, 2005. Published on CDROM. Dariusz Uciński and Maciej Patan. D-optimal design of a monitoring network for parameter estimation of distributed systems. Journal of Global Optimization, 39:291–322, 2007. Dariusz Uciński and Maciej Patan. Sensor network design for the estimation of spatially distributed processes. International Journal of Applied Mathematics and Computer Science, 20(3):459–481, 2010. Marc van de Wal and Bram de Jager. A review of methods for input/output selection. Automatica, 37:487–510, 2001. Feng Zhao and Leonidas J. Guibas. Wireless Sensor Networks: An Information Processing Approach. Morgan Kaufmann Publishers, Amsterdam, 2004. Tong Zhao and Arye Nehorai. Detecting and estimating biochemical dispersion of a moving source in a semiinfinite medium. IEEE Transactions on Signal Processing, 54(6):2213–2225, June 2006. Tomasz Zięba and Dariusz Uciński. Mobile sensor routing for parameter estimation of distributed systems using the parallel tunneling method. International Journal of Applied Mathematics and Computer Science, 18(3): 307–318, 2008.
1322