Design of control policies for spatially inhomogeneous robot swarms with application to commercial pollination

We present an approach to designing scalable, decentralized control policies that produce a desired collective behavior in a spatially inhomogeneous robotic swarm that emulates a system of chemically reacting molecules. Our approach is based on abstracting the swarm to an advection-diffusion-reaction partial differential equation model, which we solve numerically using smoothed particle hydrodynamics (SPH), a meshfree technique that is suitable for advection-dominated systems. The parameters of the macroscopic model are mapped onto the deterministic and random components of individual robot motion and the probabilities that determine stochastic robot task transitions. For very large swarms that are prohibitively expensive to simulate, the macroscopic model, which is independent of the population size, is a useful tool for synthesizing robot control policies with guarantees on performance in a top-down fashion. We illustrate our methodology by formulating a model of rabbiteye blueberry pollination by a swarm of robotic bees and using the macroscopic model to select control policies for efficient pollination.


I. INTRODUCTION
Tasks that require parallelism, redundancy, and adaptation to dynamic environments can potentially be performed very robustly and efficiently by a swarm robotic system, which would consist of hundreds or thousands of autonomous robots with limited sensing, communication, and computation capabilities.Key challenges in the development of such systems include the accurate prediction of swarm behavior and the design of robot controllers that can be proven to produce a target macroscopic outcome.The controllers should be scalable, meaning that they maintain the successful operation of the system regardless of the number of robots.
To ensure scalability, we propose a control approach that is decentralized, in which robots make decisions using only local information obtained via sensing and/or communication without knowledge of the global system state.The robots that we consider are unidentified and each programmed with the same set of control algorithms.We assume a broadcast architecture [26] in which a supervisory agent computes the parameters that govern the robots' motion and task transitions and transmits them to the swarm, without requiring knowledge of the population size or individual S. Berman and R. Nagpal are with the School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138 (e-mail: {sberman,rad}@eecs.harvard.edu) V. Kumar is with the GRASP Laboratory, University of Pennsylvania, Philadelphia, PA 19104 (e-mail: kumar@seas.upenn.edu) We gratefully acknowledge support from NSF Award CCF-0926148, ARO Grant W911NF-05-1-0219, and ONR Grants N00014-07-1-0829 and N00014-08-1-0696.S. Berman was supported by the National Science Foundation Graduate Research Fellowship and the National Defense Science and Engineering Graduate Fellowship when she was affiliated with the University of Pennsylvania.robot actions.This paradigm facilitates the communication of information to all the robots while avoiding the lossiness and inefficiencies of a multi-hop ad-hoc network.
In previous work [5], [24], we described a methodology for producing target swarm behaviors that can be implemented using this control paradigm.Our approach consists of developing an abstraction of a swarm and using this model to analyze system performance and synthesize robot control policies that cause the populations of different swarm elements, defined as robots at various tasks and objects of different types with which they interact, to evolve in a desired way.The methodology was developed for swarms in which robots execute task transitions stochastically, which allows us to produce any of a range of population distributions by designing an appropriate set of transition probability rates.For systems in which there are task transitions that are triggered by encounters between elements, we specified that swarm elements are always uniformly randomly distributed throughout the environment, which can be realized when the robots execute random walks.This allowed us to model the swarm as a well-mixed chemical reaction system, which could be abstracted to a set of ordinary differential equations.
In this paper, we extend our methodology to spatially inhomogeneous swarms, in which elements are arbitrarily distributed throughout space.We consider systems in which robots follow a deterministic velocity field in addition to random motion.We generate the continuous dynamics of individual elements as they move and interact in a physical environment using a micro-continuous model.In this model, we define the robots' motion using the technique of random walk particle tracking [10], [33], a method from statistical physics, and implement the robot task-switching using the formulation given in [4].In the limit of infinite swarm populations, the system abstracts to a set of advection-diffusionreaction (ADR) partial differential equations (PDE's) [8], the macro-continuous model, which governs the time evolution of concentrations of different elements over a spatial domain.Similar to this approach is the work [15] on developing rigorous abstractions of a swarm to a PDE model based on the Fokker-Planck equation, which describes the time evolution of the probability density function of the position of a robot.By designing the parameters of the macro-continuous model and using them to define the robot motion controllers and stochastic policies for task transitions, a supervisory agent can produce a target collective behavior.This constitutes a novel application of a "top-down" controller synthesis approach to spatially inhomogeneous swarms.
The general ADR equations that define the macro-continuous model cannot be solved analytically, and instead must be solved using numerical methods [17].A widely applied technique is a grid or mesh-based method such as finite difference, finite volume, and finite element methods [6].Grid-based methods are unsuitable for systems that consist of discrete physical entities (such as robots) rather than a continuum [21], and they can suffer from numerical dispersion and artificial oscillations when simulating advection-dominated systems [16], [33], which we may want to design.A more appropriate technique to use for our type of system is a meshfree method [20], which solves integral equations and PDE's with all kinds of boundary conditions using a set of arbitrarily distributed nodes or particles without imposing connectivity among these elements with a mesh.We numerically solve the ADR equations using smoothed particle hydrodynamics (SPH), a meshfree particle method that was originally developed to simulate astrophysical phenomena and has been extended to a variety of problems in computational fluid dynamics and solid mechanics [21], [22].We apply our methodology to synthesize control policies for a swarm of insect-inspired micro air vehicles [37] to pollinate an orchard of rabbiteye blueberries (Vaccinium ashei), a scenario of interest to the Robobees project [1], whose objective is to develop a colony of robotic bees.Insect pollination is necessary for adequate commercial yields of this crop; almost all rabbiteye cultivars are self-sterile and require cross-pollination with a compatible cultivar, of which at least two should be interplanted [9], [25].Our investigation of effective swarm strategies and their required robot capabilities will inform the design of the robotic bees.We describe micro-continuous and macro-continuous models of the scenario, demonstrate a close correspondence between the models, and illustrate how the macro-continuous model can be used to select motion controllers and stochastic policies for the robots that produce efficient pollination.

A. Robot Motion and Concentration Distribution
We consider systems in which the position q i ∈ Y ⊂ R d , d ∈ {1, 2, 3}, of each robot i (and any other swarm element) evolves according to a deterministic velocity field v(q i ) and a Brownian motion that drives diffusion, where D is the associated diffusion coefficient.The velocity field can be designed, and it may have a component that is fixed by the bulk motion of the medium, such as wind or water, in which the robots operate.Here diffusion models randomness in robot movement that either arises from inherent noise due to sensor and actuator errors or is actively added by the robot's motion controllers, or both.D is the sum of D inh , a constant determined by the first source of randomness, and D act , a tunable parameter that produces the second source.A robot i, which can be represented as a point kinematic agent, updates its displacement at each (very small) timestep ∆t according to the Itō-Taylor integration scheme [12] where I ∈ R d×d is the identity matrix and Z ∈ R d is a vector of independent, normally distributed random variables with zero mean and unit variance.A swarm of N robots that move according to equation (1) can be viewed as a random walk particle tracking (RWPT) scheme, a Lagrangian approach to solving advection-diffusion problems.In the limit N → ∞, the Fokker-Planck equation becomes equivalent to the advectiondiffusion PDE (defined in Section III), which governs the time evolution of the concentration of robots over the domain, x(q, t) [10], [33].Using the particle tracking analogy, we consider the robots to be a discrete distribution of the mass of the swarm and associate each robot i with a mass m i = m.The robot masses can be mapped to concentration values as described in [2], [36].The concentration x(q, t) at a point q ∈ Y can be approximated by the smoothed integral interpolation x(q, t) where δ(q) is the Dirac delta function and ζ(q) is a projection function with finite support that satisfies Y ζ(q)dq = 1 and is ideally invariant with respect to coordinate transformations.A particle approximation for x(q, t) is given by Most particle tracking implementations use ζ(q) = 1/V for points within a small cube of volume V around q and ζ(q) = 0 otherwise [16].This corresponds to dividing the domain Y into a grid of G cubic cells with volume V c , finding the mass of all particles inside each cell c, M c , and defining the concentration at all points q inside cell c as x(q, t) = M c /V c .

B. Robot Task Switching
The robots' transitions between tasks can be modeled as chemical reactions: the robots switch stochastically, either spontaneously or upon encountering certain objects or other robots, at rates that are determined by constants that are analogous to reaction rate constants.Borrowing chemical reaction terminology, a species X i symbolizes a robot that is performing task i or an object of type i with which the robots interact, and a complex is a combination of species that appears before or after a reaction arrow.The rate constant k ij associated with a reaction that converts complex i into complex j is a function of a parameter c ij , where c ij dt is the average probability that a particular combination of reactant elements will transition to product elements in the next infinitesimal time interval dt.
A robot's decision to switch tasks spontaneously is modeled as a unimolecular reaction, for which k ij = c ij [13].The robot control policy that implements this transition has no dependence on spatial properties of the system: a robot that is performing task i computes a uniformly distributed random number in [0, 1], u ∼ U (0, 1), at each timestep ∆t and executes the transition if u < c ij ∆t.
A robot's decision to switch tasks upon encountering another element is modeled as a bimolecular reaction.Let r = q o (t) − q p (t) be the initial relative displacement of two reactants o and p and let ∆q = q p (t + ∆t) − q o (t + ∆t) + r.The probability of the reactants occupying the same position at time t+∆t is P (∆q = r).Defining p(r) as the probability density of ∆q and specifying that all reactants are associated with the same arbitrary mass m, c ij is given by [4] We consider reactions in which one reactant is a robot that moves according to equation ( 1) and the other is stationary.The function p(r) is then the probability density of the moving reactant's displacement ∆q i (t) = v(q i (t))∆t + √ 2D∆tIZ(t), which is a bivariate normal distribution with mean v(q i (t))∆t and covariance matrix 2D∆tI.The standard deviation of the robot' As in the case of a well-mixed system [14], the robot control policy that implements the transition can be defined by decomposing the probability c ij ∆t into the product of c e ij , the probability that the robot will encounter a particular reactant within its sensing range in the next ∆t, and c d ij ∆t, the probability that the robot will decide to follow through with the transition given an encounter.We assume that a robot's sensing range is a sphere in R d whose radius r is on the order of the robot's typical random displacement per timestep, reflecting the robot's very limited sensing capacity.The probability that a robot will encounter a reactant at a position r relative to the robot in the next timestep is c e ij = Ω p(r )dr , where Ω ⊂ R d is a sphere of radius r centered at the relative position r.For small r, Ω p(r )dr ≈ p(r)V Ω , where V Ω is the volume of Ω.Using equation ( 5), we can decompose c ij as At each timestep, a robot moves according to equation ( 1) and, if it senses a reactant, computes u ∼ U (0, 1) and executes the transition if u < c d ij ∆t.The robot then becomes a species in the product of the reaction and may follow a different displacement equation of the form (1).
In simulation, we may implement an encounter-based transition in an alternate way.At each timestep ∆t, robot i is moved a deterministic displacement v(q i )∆t from equation (1).We define a sphere of radius r sim , centered at the robot's new position, which contains the vast majority of reactants that the robot is likely to reach in ∆t due to its diffusive motion.For instance, we may set r sim = 2σ.The value of c ij is computed from equation (5) for each potential reactant in this region, along with u ∼ U (0, 1), and the reaction is executed with the reactant for which u < c ij ∆t, if any.The robot is moved to the position of the chosen reactant.If no reaction occurs, the robot's diffusive motion is simulated as √ 2D∆tIZ(t) from equation (1).

III. MACRO-CONTINUOUS MODEL
The macro-continuous model of the time evolution of the species concentrations, x s (q, t), s = 1, ..., S, is given by a set of advection-diffusion-reaction (ADR) equations, which describe the conservation of chemical species in a fluid.Define v s (q) as a velocity field that specifies the deterministic motion of species s, D s as the diffusion coefficient of species s, K s as the set of rate constants k ij that are associated with the reactions involving species s, and R s (K s , x) as the sum of the corresponding reaction rates evaluated at local concentrations x(q).The ADR equations are: The ADR equation for x s is known as the advection-diffusion equation when R s (K s , x) = 0 and the reaction-diffusion equation when v s = 0. To synthesize a target macroscopic behavior in a swarm modeled by equations ( 7), a supervisory agent computes the k ij and the tunable components of v s (q) and D s , s = 1, ..., S.
The technique of smoothed particle hydrodynamics (SPH) can be used to numerically solve the ADR equations (7).In this method, a function f (q, t) over a domain is represented in terms of its values at a set of arbitrarily distributed particles.This converts the governing PDE equations of a system into a set of ODE's, each of which describes the time evolution of a variable associated with a particle.The value of a variable at a particle is influenced by values at particles within a local neighborhood only.
The SPH formulation is derived by first approximating f (q) with a smoothed integral interpolation, similarly to equation (2): Here w(q − q , h) is a differentiable kernel function with smoothing length h and compact support.The kernel w satisfies the properties where κ is a constant that defines the support domain of the kernel of point q.The kernel function is frequently defined as a Gaussian or a spline that approximates a Gaussian [21].
The integral (8) can be approximated using a Monte Carlo integration scheme as follows.The spatial distribution of the mass in the system is represented by a set of Q particles.Particle i has position r i , mass m i , and density ρ i , where m i /ρ i is the volume associated with the particle.The number of particles per unit volume, n i = ρ i /m i , is calculated as [35] Replacing dq in equation ( 8) by the volume m i /ρ i = 1/n i at the particle locations, the particle approximation for f (q) is given by The error in approximating equation ( 8) by equation ( 13) is O(h 2 ) [27].Computing the gradient of f (q) using the particle approximation entails differentiating the kernel w rather than the function itself.
In recent years, the SPH technique has been used to define decentralized motion control laws for robot swarms at the micro-continuous level to achieve pattern generation in obstacle-filled environments [30], [31] as well as deployment, sensor coverage, patrolling, flocking, and formation control behaviors [28], [29].These works model a swarm as a fluid, which may be subjected to an external force, and represent each robot as an SPH particle that has an associated position, velocity, mass, density, energy, and pressure.The robots use the SPH formulations of the governing equations for the conservation of mass, momentum, and energy, along with an equation of state for pressure, to update their associated quantities and compute their control law.The approach is scalable because these updates only require information from robots within a local neighborhood that corresponds to the support domain of w.
In contrast, our application of the SPH method at the macro-continuous level does not associate particles with robots, but rather employs the particles as computational elements that track concentration values at points in space.We use an SPH scheme similar to those that have been recently applied to model solute transport in porous media [16], [35].We define a set of particles for each of the S ≤ S distinct velocity fields v s in equations (7).The set associated with v s , s ∈ {1, ..., S }, tracks the concentrations of all species s for which v s = v s .The velocity of each particle i ∈ {1, ..., Q s } in this set is The diffusion term in equations (7), D s ∇ 2 x s , could be evaluated by differentiating the interpolated concentration twice.However, the resulting expression contains the second derivative of w, which can be noisy and sensitive to particle disorder, so instead we use an SPH discretization of the Laplace operator that involves only first order derivatives of w [18]: where ∇ rij w(r ij , h) is the directional derivative of w along r ij .For a set of irregularly spaced particles, n i = n j in general, which produces an asymmetry in the magnitude of the contribution of concentration from particle i to particle j and vice versa.To rectify this, n j is replaced by either the arithmetic or harmonic average of n i and n j [16].Choosing the harmonic average and multiplying equation (15) by D s , the particle approximation of the diffusion term becomes In the reaction rate terms of equations ( 7), the concentration at r i of a species that is not tracked by particle i is evaluated using the particle approximation (13) with the set of particles that do track the species.
The SPH formulation of the ADR equation for each species s that is tracked by particle i ∈ {1, ..., Q s } is The SPH method is implemented by initializing the particle positions and concentrations and then numerically integrating equations ( 14) and ( 18) using standard techniques such as Runge-Kutta methods and the Velocity Verlet scheme [21].

IV. APPLICATION: COMMERCIAL POLLINATION A. Micro-Continuous Model
We simulate a 50 ft × 18 ft section of a rabbiteye blueberry orchard that consists of alternating rows of two cultivars, as is commonly recommended in the literature [23].The section contains four rows of three plants each, which are spaced 6 ft apart within a row and 12 ft apart between rows, the industry standard [32].Each plant has 10 4 flowers, the quantity for a mature rabbiteye plant [34], which are uniformly randomly distributed throughout the plant canopy.Initially, a colony of robotic bees occupies a 2 ft × 3 ft hive to the left of the rows; the robots are uniformly randomly distributed throughout the hive.Fig. 1 illustrates the orchard layout and the initial placement of the robots.
The supervisory agent in our architecture can be a computer at the hive [3] with the computational resources to calculate the robot motion and task transition parameters for a specific macroscopic objective, which may involve solving the ADR equations to predict the system performance.The robots would have sufficient power to undertake brief flights away from the hive [19].They would return to the hive to recharge, upload data, and receive parameters from the hive computer.Here we simulate a portion of one flight.
We assume that each robotic bee is equipped with a compass and thus can fly with a constant heading.The robots search for flowers by flying with a deterministic component to the right superimposed with a random walk according to equation ( 1), with v(q i ) = v = [v 0] T .We specify the robot flight speed as v + 2σ/∆t = 16.4 ft/s, which is believed to be a realistic speed for robotic bees as it falls within the upper range of speeds attainable by wild orchid bees [7].
The robotic bees in this scenario are assumed to be capable of recognizing a flower that is very close by, potentially using an ultraviolet light sensor; flying to the flower; and hovering briefly while vibrating the flower to release its pollen as is done by efficient blueberry pollinators [9].The flower visits are modeled as being instantaneous.We assume that a rabbiteye flower, each of which produces 8000-17000 pollen grains [11], transfers pollen to each robot that visits it and that sufficient pollen remains on the robot to pollinate an arbitrary number of flowers on subsequent visits.
We define a set of reactions that represents the actions of pollen retrieval, pollination, and unproductive flower visits.The robot species are defined as follows: B 0 is a robot without pollen, B i is a robot carrying pollen from cultivar i ∈ {1, 2}, and B 3 is a robot carrying pollen from both cultivars.U i and P i symbolize an unpollinated flower and a pollinated flower, respectively, of cultivar i ∈ {1, 2}.W is a "waste product" indicating that a flower visit has not resulted in pollination or retrieval of a pollen type that is not already on the robot.Each reaction is associated with the same rate constant k, since it is unrealistic for the robots to be able to distinguish between flowers of different cultivars or flowers that have been visited by other robots.The reactions are: The robot behavior of encountering a flower in very close proximity with its limited sensing range and deciding whether to visit it is simulated according to the procedure described in the last paragraph of Section II-B.
We compute the concentration fields of the species as described in Section II-A.The species initially present in the simulation are B 0 , U 1 , and U 2 .The initial nonzero concentrations of B 0 were set to 1 at the center q c of each grid cell c within the hive region.The mass associated with each of the N B0 robots was computed as m = 1 The U 1 , U 2 flowers were assigned the same mass m as the robots by setting their initial nonzero concentrations to at the center of each of the G U k grid cells in the rows of cultivar k ∈ {1, 2}.To initialize the robot and flower positions, we computed their populations per cell as N c s = x s (q c , 0)V c /m, s ∈ {B 0 , U 1 , U 2 }, rounded this number to the nearest integer, and distributed the N c s robots uniformly randomly inside the cell.

B. Macro-Continuous Model
The time evolution of the concentrations of the species in the set of reactions ( 19) is governed by the following ADR equations, where Our SPH formulation of the system uses a set B of Q B particles with positions r i that move with velocity v and track the concentrations of robot species, as well as a set F of Q F stationary particles with positions s j that track the concentrations of flower species and unproductive flower visits.Each set is arranged on a lattice with particle spacing ∆r = h/h 0 , where h 0 is a constant.The SPH formulation of model ( 20), ( 21) consists of equation ( 14) with v s (r i ) = v and v s (s j ) = 0 and the set of equations for dx B k /dt| ri , k = 0, 1, 2, 3; dx U l /dt| sj , dx P l /dt| sj , l = 1, 2; and dx W /dt| sj that are defined by equation (18).In the reaction rate terms of equation (18), the concentration of a species s that is tracked by one set of particles at a particle in the other set is approximated according to equation (13): where l = 1, 2, m = 0, 1, 2, 3, and n P k is the number of particles in set P per unit volume evaluated at particle k.We define the kernel w as a cubic spline with compact support and a shape similar to a Gaussian: where R = ||r ij ||/h and γ = 15 7πh 2 in two dimensions.This function has been widely used in the SPH literature [21].
Because the lattice spacings do not change over time, we can precompute the identities of a particle's neighbors in its lattice (i.e., those within the support domain of the kernel w), as well as all quantities derived from its distance to its neighbors.In general, however, these variables must be updated at every iteration.The initial nonzero species concentrations are x B0 (r i ) = 1 at B particles within the hive region and x U k (s j ) = x 0 U k at F particles in the rows of cultivar k ∈ {1, 2}.At each timestep, the particle positions and concentrations are computed by numerically integrating the SPH equations using the Euler method.

C. Results
For the analysis in this section, we ran the microcontinuous model and the SPH formulation of the macrocontinuous model with ∆t = 0.002 s.The micro-continuous model used a grid of 50×18 cells over the domain [0 50] f t× [0 18] f t.In the SPH method, set B was comprised of a lattice of 53 × 35 particles that initially covered the domain [− 11 15] f t×[0 17] f t, and set F was a lattice of 101 × 37 particles that covered the domain [0 50] f t×[0 18] f t.Using the suggested values of h 0 = 1.2 and κ = 2 [21] as a guideline, we defined h 0 = 2 and κ = 2.
1) Accuracy of the Macro-Continuous Model: To verify that the macro-continuous model is an accurate abstraction of the system, we compared concentrations computed by the micro-and macro-continuous models for the four parameter sets in Table I.The quantity σ/(v∆t) is a ratio of diffusive to deterministic robot motion per timestep, where v and σ are subject to the flight speed constraint in Section IV-A.The micro-continuous model simulated 10 4 robotic bees.Fig. 2 and 3 show that there is a close match between the concentrations of B 3 , P 1 , P 2 , and W computed by the two models at specified times for parameter sets 1 and 2. Note that raising σ/(v∆t) results in a broader pollinated region and higher maximum concentrations of both pollinated flowers and unproductive flower visits.The first row of cultivar 1 is not pollinated because the robots visiting that row have not yet acquired pollen from cultivar 2.
We quantify the degree of closeness between the models with the error metric , where x micro s is the vector of concentrations of species s computed at each of the G cell centers in the micro-continuous model and x macro s is the corresponding vector computed at SPH particles in set F that overlap the cell centers.Fig. 4 shows that the concentration fields of P 1 and P 2 computed with parameter sets 3 and 4 have similar degrees of closeness to those that are visually compared in Fig. 2 and 3, respectively.
2) Macro-Continuous Model as a Tool for Top-Down Design: The correspondence between the micro-and macrocontinuous models demonstrates that the macro-continuous model can be used to design the parameters v, D, and k that will produce desired collective behaviors in the physical system.The advantage of using the macro-continuous model for this purpose is that, unlike the micro-continuous model, its computational time is invariant to the sizes of the species populations.Fig. 5 demonstrates this advantage with the ratio of τ micro , the runtime on a standard 2.53 GHz laptop per iteration of the micro-continuous model averaged over 201 iterations in the interval t = 1.6 − 2.0 s, to τ macro = 1.42 s, the same quantity in the macro-continuous model.As the number of robots N increases above ∼ 3000, it is faster to run the macro-continuous model.Hence, for very large populations, it is more suitable to use the macro-continuous model in an optimization technique such as a Monte Carlo method to compute the parameters that maximize some metric of success.We illustrate this concept by using the macro-continuous model to find the value of k that maximizes a metric of pollination efficiency for fixed v and D. Computing the mass of species s We define η poll as the mass fraction of pollinated flowers, (M P1 + M P2 )/M 0 , and η waste as the ratio of the "mass" of unproductive flower visits to the total mass of flowers, M W /M 0 .Thus, η poll /η waste is a measure of the number of pollination events to the number of unproductive visits.As Fig. 6 this metric is a maximum at k = 3 and then decreases below 1 for k > 9, since η waste increases faster than η poll .

V. CONCLUSIONS AND FUTURE WORK
We have described an approach to synthesizing robot control policies that produce a target collective behavior in a spatially inhomogeneous swarm that emulates an advectiondiffusion-reaction chemical system.Our approach relies on a rigorous correspondence between two models of a swarm.The micro-continuous model offers a realistic representation of individual robot activities and captures variability in macroscopic outcomes that an actual system would produce over multiple trials.However, it becomes more computationally intensive to run as the swarm population is increased, and studying the effect of parameter variation requires simulations under many different conditions.These factors motivate us to abstract the system to a macro-continuous model, which is independent of the population size and amenable to techniques for analysis, control, and optimization that yield theoretical guarantees on performance.This model becomes a more accurate description as the swarm population increases, and it does not capture variability in performance.Using the application of commercial pollination with robotic insects, we demonstrate the utility of this model in designing the parameters that govern the deterministic and random components of the robot motion and the robot probabilistic task transitions.
Future work includes the investigation of analysis and controller synthesis methodologies for the PDE macrocontinuous model.These may include nonlinear stability analysis, nonlinear and robust control methods for parabolic PDE's, and stochastic optimization techniques.Parameter optimization in higher dimensions can be facilitated by speeding the simulation execution time with parallelization of the SPH code [21] and the optimization technique.
We are also interested in applying our controller synthesis methodology to expanded models of the pollination scenario.We would like to model different hive placements, more complex velocity fields, time delays associated with pollination, and inter-robot communication.We can add unimolecular reactions that represent a robotic bee failing in the field or returning to the hive at a time-dependent rate that is related to its energy consumption.The flower visit rate constants can be modified to vary with time, depend on the presence of pollen on a robot, and incorporate the probability of fertilization given pollination [34].Finally, we would like to include feedback from the robots about their flower visits upon their return to the hive.The hive computer would use this data to recalculate the parameters for the next round of flights to produce robot behavior that compensates for pollination deficiencies due to environmental disturbances and robot failures.

Fig. 1 .
Fig. 1.An initial distribution of robotic bees (black) and unpollinated flowers of cultivars 1 (yellow) and 2 (cyan) in the micro-continuous model.

Fig. 3 .
Fig. 3. Snapshots of the concentration distributions of B 3 at time t = 7.5 s and P 1 , P 2 , and W at t = 15 s over subsets of the domain [0 50] f t × [0 18] f t.Concentrations were computed for parameter set 2 using the micro-continuous model (left column) with 10 4 robots and the SPH method (right column).

Fig. 4 .Fig. 5 .
Fig.4.Normalized errors µ P 1 , µ P 2 over time for all four parameter sets and 10 4 robots.Each plot ends at a time t when the robotic bees have covered an average x distance slightly larger than vt = 50 f t, the width of the field.

Fig. 6 .
Fig.6.Values of pollination metrics η poll , ηwaste, and η poll /ηwaste computed at t = 6.2 s using the macro-continuous model over a range of k with v = 8.20 f t/s, D = 0.01681 f t 2 /s.
Snapshots of the concentration distributions of B 3 at time t = 3 s and P 1 , P 2 , and W at t = 6 s over subsets of the domain [0 50] f t × [0 18] f t.Concentrations were computed for parameter set 1 using the micro-continuous model (left column) with 10 4 robots and the SPH method (right column).