The problem is autonomous coordinated search by an interconnected group of moving robots for the purpose of finding and localising a source of hazardous emissions (e.g., gas and particles). Dispersion of the emitted substance is assumed to be affected by turbulence, resulting in the absence of concentration gradients. The chapter proposes a search strategy that operates in a completely decentralised manner, as long as the communication network of the moving robots forms a connected graph. By decentralised operation, we mean that each moving robot is reasoning (i.e., estimating the source location and making decisions on robot motion) locally. Coordination of the group is achieved by consensus via communication with the neighbours only, in a manner which does not require global knowledge of the communication network topology.
- autonomous search
- machine intelligence
- sequential Monte Carlo estimation
Searching strategies for finding targets using appropriate sensing modalities are of great importance in many aspects of life. In the context of national security, there could be a need to find a source of hazardous emissions [1, 2, 3]. Similarly, rescue and recovery missions may be tasked with localising a lost piece of equipment that is emitting weak signals . Biological applications include, for example, protein searching for its specific target site on DNA , or foraging behaviour of animals in their search for food or a mate [6, 7]. The objective of search research  is to develop optimal strategies for localising a target in the shortest time (on average), for a given search volume and sensing characteristics.
The use of autonomous vehicles in dangerous missions, such as finding a source of hazardous emissions, has become widespread [9, 10, 11]. Existing approaches to the search and localisation in the context of atmospheric releases can be loosely divided into three categories: up-flow motion methods, concentration gradient-based methods and information gain-based methods, also known as
This chapter summarizes our recent results in development of an autonomous infotaxic coordinated search strategy for a group of robots, searching for an emitting hazardous source in open terrain under turbulent conditions. The assumption is that the search platforms can move and sense. Two types of sensor measurements are collected sequentially: (a) the concentration of the hazardous substance; (b) the platform location within the search domain. Due to the turbulent transport of the emitted substance, the concentration measurements are typically sporadic and fluctuating. The searching platforms form a moving sensor network, thus enabling the exchange of data and a cooperative behaviour. The multi-robot infotaxis have already been studied in [16, 17, 20, 24]. However, all mentioned references assumed
We develop an approach where the group of searching robots operate in a fully decentralised coordinated manner. Decentralised operation means that each searching robot performs the computations (i.e., source estimation and path planning) locally and independently of other platforms. Having a common task, however the robotic platforms must perform in a coordinated manner. This coordination is achieved by exchanging the data with immediate neighbours only, in a manner which does not require the global knowledge of the communication network topology. For this reason, the proposed approach is scalable in the sense that the complexities for sensing, communication, and computing per sensor platform are independent of the sensor network size. In addition, because all sensor platforms are treated equally (no leader-follower hierarchy), this approach is robust to the failure of any of the searching agents. The only requirement for avoiding the break-up of the searching formation is that the communication graph of the sensor network remains
2. Mathematical models
First, we describe the measurement model. The concentration measurements are modelled using a Lagrange encounters model developed in , based on an open field assumption and a two-dimensional geometry. Let th robotic vehicle position () at time be denoted by . Suppose that the emitting source is located at coordinates specified by the vector and its release rate, or strength, is . The goal of the search is to detect and estimate the source-parameter vector in the shortest possible time. The particles released from the source propagate with combined molecular and turbulent isotropic diffusivity , but can also be advected by wind. The released particles have an average lifetime before being absorbed. Let the
where , and are known environmental parameters, is the distance between the source and the th sensor platform, is the modified Bessel function of the second kind of order zero, and depends on environmental parameters only.
The probability that a sensor at location is hit by dispersed particles (where is a non-negative integer) during a time interval is Poisson distributed, i.e.,
Parameter in (2) is the mean number of particles expected to reach the sensor at location during interval . Eq. (2) expressed the likelihood function of a concentration measurement collected by th sensor, i.e., .
The motion model of a coordinated group of robots is described next. Let the pose vector of the th robot platform at time be denoted , where has already been introduced and is the vehicle heading. The group of searching vehicles moves in a formation. The centroid of the formation at time is specified by coordinates:
For each platform , the offset from the centroid is predefined and known to it (i.e., , ).
The measurements of concentration are taken at time instants , . Between two consecutive sensing instants, each platform is moving. Let the duration of this interval (referred to as the
Motion of the th platform during interval is controlled by linear velocity and angular velocity . Given that the motion control vector is applied to the th platform, its dynamics during a short integration time interval can be modelled by a Markov process whose transitional density is . The process noise covariance matrix captures the uncertainty in motion due to the unforeseen disturbances. The vehicle motion function is:
where vector is introduced to compensate for a distortion of the formation due to process noise with parameters:
Here, and are the estimates of the coordinates of the formation centroid at (that is of and , respectively) available to the th platform. Coordinates and refer to the
A robotic platform can communicate with another platform of the formation, if their mutual distance is smaller than a certain range . Because of process noise in motion, the distance between the vehicles in the formation will vary and consequently the topology of the communication network graph may also vary. For simplicity, we will assume that communication links (when established) are error free. Figure 1 illustrates the communication graphs of a formation consisting of searching platforms at two consecutive time instants.
3. Decentralised sequential estimation
Estimation and robot motion control are carried out using the measurement dissemination-based decentralised fusion architecture . Measurement locations2 and the corresponding measured concentration values, i.e., the triple , are exchanged via the communication network. The protocol is iterative. In the first iteration, platform broadcasts its triple to its neighbours and receives from them their measurement triples. In the second, third and all subsequent iterations, platform broadcasts its newly acquired triples to the neighbours, and accepts from them only the triples that this platform has not seen before (newly acquired). Providing that the communication graph is connected, after a sufficient number of iterations (which depends on the topology of the graph), a complete list of measurement triples from all platforms in the formation, denoted , will be available at each platform.
Suppose the posterior density function of the source at discrete-time and platform be denoted , where . Given and , the problem of sequential estimation is to compute the posterior at time , i.e., . Using the Bayes rule, the posterior is
where is the likelihood function. Assuming that individual platform measurements are conditionally independent, can be expressed as
is independent of . The posterior density is computed using the Rao-Blackwell dimension reduction scheme . Using the chain rule, the posterior can be expressed as:
where the posterior of source strength will be worked out analytically, while the posterior of source position will be computed using a particle filter. Following , we express the posterior with the Gamma distribution whose shape and scale parameters are and , respectively. That is
Since the conjugate prior of the Poisson distribution is the Gamma distribution , the posterior is also a Gamma distribution with updated parameters and , i.e., . The computation of and can be carried out analytically as a function of and the measurement set :
The parameters of the prior for source strength, are chosen so that this density covers a large span of possible values of .
Next, we turn our attention to the posterior of source position in the factorised form (8). Given , the update step of the particle filter using applies the Bayes rule:
where is a normalisation constant. The problem in using (11) is that the likelihood function is unknown; only of (6) is known. Fortunately, it is possible to derive an analytic expression for :
The Rao-Blackwellised particle filter (RBPF) fully describes the posterior by a particle system
Here, is the number of particles, is a (normalised) weight associated with the source position sample , while and are the parameters of the corresponding Gamma distribution for the source strength. Initially, at time , the weights are uniform (and equal to ), are the points on a regular grid covering a specified search area, while and . The sequential computation of the posterior using the RBPF is carried out by a recursive update of the particle system over time.
4. Decentralised formation control
In decentralised multi-robot search, each platform autonomously makes a decision at time about its next control vector , as described in Section 4.1. However, in order to maintain the geometric shape of the formation and thus avoid its break-up, there is a need to impose a form of coordination between the platforms. This will be explained in Section 4.2.
4.1 Selection of individual control vectors
A robot platform autonomously decides on the control vector using the infotaxis strategy , which can be formulated as a partially observed Markov decision process (POMDP) . The elements of POMDP are (i) the information state, (ii) the set of admissible actions and (iii) the reward function. The information state at time is the posterior density ; it accurately specifies the th platform current knowledge about the source position and its release rate. Admissible actions can be formed with one or multiple steps ahead. A decision in the context of search is the selection of a motion control vector which will maximise the reward function. According to Section 2, the space of admissible actions is continuous with dimensions: linear velocity , angular velocity and duration of motion . In order to reduce the computational complexity of numerical optimisation, is adopted as a discrete set with only myopic (one step ahead) controls. In addition, is time-invariant and identical for all platforms. If , and denote the sets of possible discrete-values of , and , respectively, then is the Cartesian product . The myopic selection of the control vector at time on platform is expressed as:
where is the reward function and is the future concentration measurement collected by the th platform if the platform moved under the control to position . In reality, this future measurement is not available (the decision has to be made at time ), and therefore the expectation operator with respect to the prior measurement PDF features in (13).
where is the current differential entropy, defined as
while is the future differential entropy (after a hypothetical control vector has been applied to collect ):
where . The expectation operator in (14) is with respect to the probability mass function , that is:
Given that is approximated by a particle system , one can approximately compute , which features in (14), as
In order to compute of (17), first note that , where is the predicted mean rate of chemical particle encounters at location (where the platform would move after applying a hypothetical control ), computed based on . According to Section 2,
where the product approximates the source release rate as the mean of the Gamma distribution with parameters . Next, we find the value of as the minimum value of such that the cumulative distribution function is greater than a certain threshold , where . The summation (17) is then computed only for . Computation of is carried out according to (18), except that is replaced with , where
Thus, (17) is approximated with
Pseudo-code of the routine for the computation of control vector on platform is given by Algorithm 1.
Algorithm 1 Computation of
2: Compute using (18)
3: Create admissible set
5: Compute the future platform location
6: Compute using (19)
7: Determine s.t.
8: Compute using (20)
9: Calculate the expected reward using (14)
11: Find using (13)
4.2 Cooperative control through consensus
So far, we have explained how platform would independent of the other platforms in the formation determine the best action for itself, i.e., . In general, individual platforms will disagree on the best action, and in the extreme . In order to maintain the shape of the formation during the motion period (from time to ), the platforms need to reach an agreement on the common action , to be applied to all platforms at the same time. But this is not sufficient; according to the motion model in Section 2, the platforms also need to agree on the formation centroid coordinates and the common heading angle to be applied in (4).
We apply decentralised cooperative control based on the average consensus [30, 31]. In a network of collaborating agents, consensus is an iterative protocol designed to reach an agreement regarding a certain quantity of interest. Suppose that every platform, as a node in the communication network, initially has an individual scalar value. The goal of average consensus is for every node in the network to compute the average of initial scalar values, in a completely decentralised manner: by communicating only with the neighbours in the communication graph (without knowing the topology of the communication graph).
In the problem we consider, there is not only a single individual scalar value, but six of them. They include three motion control parameters, i.e., for platform , , and , two formation centroid coordinates, i.e., , and the heading angle of each platform .
Let us denote the scalar value of interest by , that is
Ideally, we want every platform in the formation to compute the mean value . If all platforms in the formation were to use identical average values for motion control, centroid coordinates and heading, then their motion would be coordinated (except for process noise, which will be taken care of through vector in (4)) and the shape of the formation would be maintained (provided is adequate).
Average consensus is an iterative algorithm. At iteration , the node in the communication graph (the robot platform) will initialise its state using either a component of vector (if is a motion control parameter) or the platform pose (if is a formation centroid coordinate or heading angle). This value is locally available. The initial values of centroid coordinates are the actual th platform coordinates, i.e., and . At each following iteration , each platform updates its state with a linear combination of its own state and the states of its current neighbours. Let us denote the set of current neighbours of platform by . Then :
where is the number of neighbours of platform . This particular linear combination is based on the so-called
The search continues until the global stopping criterion is satisfied. The local stopping criterion is calculated on each platform independently based on the spread of the local positional particles , measured by the square-root of the trace of its sample covariance matrix . For example, if the spread of particles on platform is smaller than a certain threshold , then the local stopping criterion is satisfied and is given a value of one, otherwise it is zero. This local stopping criterion value (zero or one) becomes the initial state of the global stopping criterion on platform , denoted :
The global stopping criterion is computed on each platform using the average consensus algorithm, using (21), but with replaced by . After a sufficient number of iterations, , platform decides to stop the search if at least one of the platforms in the formation has reached the local stopping criterion, that is, if .
We point out that both estimation and control are based on the consensus algorithm. While the cooperative control is using the
5. Numerical results
The proposed search algorithm has been applied to an experimental dataset, collected by COANDA Research & Development Corporation using their large recirculating water channel. The emitting source was releasing fluorescent dye at a constant rate from a narrow tube. The dataset comprises a sequence of 340 frames of instantaneous concentration field measurements in the vertical plane and is sampled at every 10/23 s. The size of a frame is pixels, where a pixel corresponds to a square area of . As the size of the data is relatively small, we follow the approach used in : upscale each frame by a factor of 3 using bicubic interpolation and place the result in the top left corner of a search area. A measurement obtained by a platform is, thus, the integer value of the concentration of the dye taken from the closest spatial and temporal sample from the experimental data.
An example of the search algorithm running on the experimental data is shown in Figure 2. All physical quantities are in arbitrary units (a.u.). The following environmental/sensing parameters were used: , , , and . Algorithm parameters are selected as follows: , , number of particles , , degrees per unit of time and . The number of iterations, both for the exchange of measurement triples and in the consensus algorithm, was fixed to 30. The local search stopping threshold was .
Figure 2 displays the top-down view of the search progress at step indices . The formation consists of platforms, whose trajectories are shown in different colours. The search algorithm terminated at . Note that the plume size is much smaller than the search area. Panels
Using 200 Monte Carlo simulations, the mean search time for the algorithm was 2525 a.u., with a 5th and 95th quantile of 1840 and 3445 a.u., respectively. Note that in all simulations the formation started from the bottom right hand corner indicated in Figure 2(a).
The chapter presented a decentralised infotaxic search algorithm for a group of autonomous robotic platforms. The algorithm allows the platforms to search and locate a source of hazardous emissions in a coordinated manner without the need for a centralised fusion and control system. More precisely, this distributed coordination is achieved only by local exchange of measurement data between neighbouring platforms. Similarly, the movement decisions taken by the platforms were reached using a distributed average consensus algorithm over the whole formation. The key aspect is that individual platforms only require knowledge of their neighbours; the global knowledge of the communication network topology is unnecessary. An advantage of adopted distributed framework is that all platforms are treated equally, making the proposed search algorithm scalable and robust to the failure of a single platform. Numerical results using experimental data confirmed the robust performance of the algorithm. The main limitation of the algorithm is that the environmental parameters (such as diffusivity, the average direction and speed of the wind, particle lifetime), must be known. Future work will explore sensitivity to parametrisation and will aim to develop a team of “search and rescue” robots for further experimentation in realistic environments.
This research was supported in part by the Defence Science and Technology Group through its Strategic Research Initiative on Trusted Autonomous Systems.
- Robot locations are assumed to be non-coincidental with the source location r0.
- Because the measurement locations are assumed to be known exactly, they will not be treated as random variables.