{'fraction_non_alphanumeric': 0.04658969764329621, 'fraction_numerical': 0.01601379091342119, 'mean_word_length': 4.237970773434715, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
PDE-Based Optimization for Stochastic Mapping and Coverage Strategies using Robotic Ensembles

This paper presents a novel partial differential equation (PDE)-based framework for controlling an ensemble of robots, which have limited sensing and actuation capabilities and exhibit stochastic behaviors, to perform mapping and coverage tasks. We model the ensemble population dynamics as an advection-diffusion-reaction PDE model and formulate the mapping and coverage tasks as identification and control problems for this model. In the mapping task, robots are deployed over a closed domain to gather data, which is unlocalized and independent of robot identities, for reconstructing the unknown spatial distribution of a region of interest. We frame this task as a convex optimization problem whose solution represents the region as a spatially-dependent coefficient in the PDE model. We then consider a coverage problem in which the robots must perform a desired activity at a programmable probability rate to achieve a target spatial distribution of activity over the reconstructed region of interest. We formulate this task as an optimal control problem in which the PDE model is expressed as a bilinear control system, with the robots' coverage activity rate and velocity field defined as the control inputs. We validate our approach with simulations of a combined mapping and coverage scenario in two environments with three target coverage distributions.


Introduction
Over the past few decades, partial differential equation (PDE) models of multiagent systems have been used extensively in mathematical biology to analyze collective behaviors such as chemotaxis, flocking, schooling, predator-prey interactions, and pattern formation [33]. Many of these models are linear or nonlinear advection-diffusion type PDEs, which describe the spatiotemporal evolution of probability densities of agents. Mathematical tools such as bifurcation analysis, optimization, and control theory can be applied to these continuum macroscopic models to make qualitative and quantitative predictions about the system behavior. Typically, each PDE model corresponds to a discrete microscopic model that captures the stochastic and deterministic actions of individual agents. While these microscopic models are more accurate descriptions of the agents' behavior, the macroscopic models enable tractable analysis for large agent numbers.
Recently, this work has motivated the use of similar types of PDEs to model and control the spatiotemporal dynamics of very large collectives, or swarms [15], of small, resource-constrained robots (e.g., [24,39]) that are currently being developed for applications such as environmental monitoring, exploration, surveillance, disaster response, and biomedical procedures. PDEs have been used to characterize the distributions of chemotactic robots in a diffusive fluid environment [17], miniature robots inspecting a model of jet turbine blades [37], and honeybee-inspired agents that aggregate at the optimal value of a scalar field [7]. The parameters of these PDE models can be mapped to control inputs that drive the robots' motion and probability rates of switching between states or tasks, and the collective behavior of the robots follows the PDE model prediction in expectation. Several works have exploited this correspondence to control the spatial distribution of a ensemble [31,14]. These control approaches can be viewed as extensions of stochastic task allocation schemes based on nonspatial rate equation models [5,8,28,27]. Other applications of continuum population dynamical models to multi-agent control include optimized confinement strategies [22], consensus using the theory of mean field controls [32], controlled flocking [35] that includes non-parallel motions [21], and pattern generation in the presence of obstacles [36]. There has also been some recent work on using PDEs to model Laplacian network dynamics of agents for formation control; see [16,29,11] and references therein.
We apply this PDE-based modeling framework to develop a control approach for allocating tasks among an ensemble of robots. In our scenarios, a task is defined as a desired activity that a robot performs in a certain spatial region of the environment. The tasks can be performed in parallel, and multiple robots can be simultaneously allocated to each task. While various deterministic approaches have been developed for multi-robot task allocation, including centralized and decentralized market-based techniques [9,6] and centralized methods for optimal task assignment and trajectory planning [1,41], their computation and/or communication requirements do not scale well to very large numbers of robots and tasks. In contrast to these works, we develop a stochastic approach in which tasks are performed at random times by unidentified robots with limited computing capabilities and no communication or global localization. Such limitations will be common in swarm robotic platforms, e.g. micro aerial vehicles [24] and microrobots [39], and in scenarios where the robots operate in GPS-denied environments where communication is impractical or unreliable. In our proposed approach, a task allocation emerges from the collective ensemble activity.
We first consider a mapping problem in which the objective is to estimate a scalar spatial field from unlocalized data obtained by the robots. We then define a coverage problem in which the ensemble must produce a target spatial density of activity over a region of interest, which may be estimated in the mapping problem. For this problem, we express the PDE model as a bilinear control system [4] and formulate an optimal control problem that computes the control inputs. Since we do not assume that agents are capable of global localization or estimation of the local agent population density, we frame the coverage problem as an open-loop control problem that does not require feedback on agent positions or densities. We follow the variational approach described in [40] for optimal control of the PDE model. While there has been some prior work on bilinear optimal control of systems of PDEs [34,3], these works do not address the types of PDEs that we consider. An optimal control problem for a bilinear parabolic PDE was formulated in [34] with the diffusion coefficient as the control. In [3], bilinear control of a class of advection-reaction systems was considered; unlike our PDE models, these systems did not include diffusion.
The mapping and controller synthesis approaches described in this paper require a central supervisor with the computational capabilities necessary to solve the associated optimization problems. Despite this centralized component, the approaches are scalable with the number of agents in the ensemble since each agent executes the same controllers with the same control variables, which are preprogrammed or broadcast by the supervisor. In our coverage strategy, there are only three control variables to be computed; in contrast, the most naive approach to controlling an ensemble of N agents moving in d dimensions would require the computation of N d control inputs.
We first presented our coverage approach in [12], where we introduced a similar optimal control problem, derived the gradient of the objective functional with respect to the control parameters, and used a gradient descent algorithm to compute the optimal control. This paper provides a complete analysis of our approach in [12] by investigating the well-posedness of the PDE model and the optimal control problem. The theory of weak solutions that we use to establish the well-posedness of the PDE model is classical [13]. However, to the best of our knowledge, there have been no prior results on well-posedness that can be directly applied to our model, which is a system of PDEs in which diffusion is present only in one of the species, the control variables are time-dependent, and a zero-flux boundary condition is imposed on the boundary of a Lipschitz domain. In this paper, we prove the existence and uniqueness of solutions of our PDE model by deriving suitable energy estimates for the solutions. We also use these derived energy estimates to ensure that the computation of the gradient, performed using the adjoint equation approach, is well-posed. Moreover, we prove the existence of an optimal control for the problem using standard compactness arguments adapted to the PDE control setting [40]. In addition to this analysis, our formulation of the mapping problem in the same framework is a novel contribution of this paper; in [12] it was assumed that the environment is known beforehand.
The paper is organized as follows. Section 2 describes the robot capabilities and their programmed behaviors during the mapping and coverage assignments, and Section 3 defines the microscopic and macroscopic models of the ensemble and its activity during each assignment. Section 4 defines key mathematical terminology that is used in Sections 5 and 6 to formulate and analyze the mapping and coverage objectives, respectively, as optimization problems that incorporate the macroscopic models. We validate our approach in Section 7 with simulations in which a region of interest must first be mapped and then covered with a target distribution of robot activity, and we conclude in Section 8.

Task Descriptions and Assumptions
We consider a scenario in which (1) a small number of agents must map a region of interest in an unknown, bounded environment, which we refer to as the mapping assignment, and then (2) a larger ensemble of agents must produce a target spatial distribution of activity within the mapped regions, which we call the coverage assignment. For instance, this activity could consist of sensor measurements, or as in our previous work [12], contacts with flowers to effect crop pollination. The mapping and coverage assignments will be formulated in a decoupled manner by posing them as two separate optimization problems in terms of their associated mean-field PDE models. We will then demonstrate through numerical simulations that these two problems can be solved sequentially in order to achieve the desired coverage objective.

Robot capabilities
We assume that the agents lack global localization, inter-agent communication, and prior data about the environment. Each agent is equipped with local sensing capabilities, allowing it to detect and distinguish between different types of regions within its sensing range, and a compass, enabling it to move with a specified heading. Additionally, the agents have sufficient memory to store the times at which they record observations of regions of interest. Similarly, it is assumed that agents have sufficient memory to store time-dependent velocity and task-switching parameters for the coverage assignment.
Both the mapping and coverage assignments involve numerical optimization computations that are performed offline by an external supervisor. After the mapping assignment, the supervisor collects recorded information from the robots and uses this information to reconstruct the map of the environment. We note that the supervisor does not have information on the individual identities of the robots. Based on this map, the supervisor calculates the parameters of the agent behaviors for a specified coverage objective and broadcasts these parameters to the agents before they are deployed for the coverage assignment. This broadcast architecture has been previously proposed for controlling large numbers of robots [30]. Although the external supervisor constitutes a potential single point of failure, it enables the control of large ensembles of agents that are subject to the constraints that we consider.

Robot controllers
The agents traverse the environment, a bounded open connected set Ω ⊂ R 2 with Lipschitz continuous boundary ∂Ω, with a combination of deterministic and stochastic motion during a deployment. Each agent moves with a timedependent velocity field v(t) ∈ R 2 , which may be broadcast to the agents or imposed on them using external stimuli, such as magnetic fields in microrobotic applications. Concurrently, each agent exhibits random motion that may be programmed, for instance to perform probabilistic search and tracking, or that arises from inherent sensor and actuator noise. We assume that this random movement can be modeled as a Brownian motion that drives diffusion with an associated diffusion coefficient D. This approach to modeling the motion of members of a robotic swarm has, for example, been used previously in [20].
The agents switch stochastically between behavioral states at constant or time-dependent transition rates, which define the probability per unit time for an agent to switch from one state to another. We define a region of interest, which may be disconnected and is assumed to be Lebesgue measurable, as Γ ⊂ Ω. During mapping, a continually moving agent records observations of Γ at a constant rate k o . During coverage, an agent that is moving over Γ pauses to perform an action (such as a sensor measurement) at a time-dependent rate k(t), and it resumes moving at a constant rate k f , which determines the time taken to complete the action. This coverage strategy can be extended to scenarios where there are different types of regions, and the agents perform actions over each region at a different rate, as in our prior work [12].
We specify that the agents' velocity field and transition rates are controllable parameters. In the mapping assignment, v(t) is designed to ensure thorough exploration of the environment (for instance, following a lawnmower pattern), and k o is assigned a high value to yield a large number of observations and thus produce an accurate map. In the coverage assignment, k f is fixed while v(t) and k(t) are computed prior to the agents' deployment according to the optimal control approach in Section 6.

Microscopic Models
The microscopic model is used to simulate the individual robots' motion and probabilistic decisions that are produced by the robot controllers in Section 2.2.

Remark 1
Here we describe the microscopic model at a formal level as a discretetime stochastic process. A rigorous correspondence between the microscopic model and the macroscopic model for the case Ω = R 2 was shown recently by the authors and collaborators in [43].
We model a robot's changes in state and performance of desired actions as a Chemical Reaction Network (CRN) in which the possible species are M , a moving robot; S, a stationary robot; and A, an instance of a desired robot activity. These reactions can only occur when a robot is located in the region of interest Γ. For example, in an artificial pollination scenario [12], Γ could represent the subset of the domain in which flowers are distributed. A robot's mapping activity consists of a single irreversible reaction, where A is the robot's observation of the region of interest. A robot's coverage activity consists of two irreversible reactions, where A is a desired action that the robot performs.
In addition, we model the displacement of a robot i over a timestep ∆t using the standard-form Langevin equation [18], where X i (t) ∈ R 2 is the position of robot i ∈ {1, .., N } at time t and Z i (t) ∈ R 2 is a vector of independent standard normal random variables that are generated at time t. Here, Y i (t) is a binary variable associated with the discrete state of robot i at time t. When the robot is performing the mapping assignment, Y i (t) = 1 for all t. When it is performing the coverage assignment, where H Γ is the indicator function of the set Γ and P is the probability measure induced by the process (X, Y ) on the set of sample paths. The robot avoids collisions with the domain boundary by performing a specular reflection when it encounters this boundary. We emphasize that the velocity v(t) and diffusion coefficient D are the same for each robot i. This assumption enables the coverage problem to be posed as a PDE control problem, as shown in Sections 5 and 6.

Macroscopic Models
The macroscopic model consists of a set of advection-diffusion-reaction (ADR) PDEs that describe the evolution of probability densities of agents that follow the microscopic model, conditioned on the initial distribution of the agents' states. Since we assume that there are no interactions between the agents, the random variables associated with the agents are independent and identically distributed, and hence we can drop the subscript i from X i (t), Y i (t), and Z i (t) and represent the macroscopic model as a single system of PDEs rather than N systems of PDEs. The state variables of this model are the population density fields of moving robots, y 1 (x, t); stationary robots, denoted by y 2 (x, t) in the coverage model; and instances of the desired robot activity, denoted by y 2 (x, t) in the mapping model and by y 3 (x, t) in the coverage model. We define Q = Ω×(0, T ) and Σ = ∂Ω×(0, T ) for some fixed final time T . The vector n(x) ∈ R 2 is the outward normal to ∂Ω, defined for almost every x ∈ ∂Ω. We also define an initial density of moving robots y 0 (x) over Ω that is normalized such that The macroscopic model of mapping is given by: with the no-flux boundary condition and initial conditions y 1 (·, 0) = y 0 , y 2 (·, 0) = 0 on Ω, where 0 is the function on Ω that takes the value 0 almost everywhere. The microscopic model of mapping defined in Section 3.1 is related to the above macroscopic model by the approximation P( is the expected number of observations recorded by the robots in the region Λ up until time t.
The macroscopic model of coverage is defined as: with the no-flux boundary condition (6) and initial conditions In this macroscopic model, the density fields y 1 (x, t) and y 2 (x, t) are related to the microscopic model of coverage, given by the stochastic process (X(t), Y (t)) defined in Section 3.1, as follows: where Λ is a measurable subset of Ω. Additionally, Λ y 3 (x, t)dx represents the expected number of instances of robot coverage activity performed in the region Λ up until time t.

Mathematical Preliminaries
We now define some mathematical terms that are used in later sections. Given two Banach spaces P and Q, L(P, Q) is the space of bounded operators from P to Q. For a finite collection of Banach spaces, {Z 1 , ..., Z M }, we equip the Cartesian product of the spaces, We define L 2 (Ω) as the Hilbert space of square integrable functions over Ω, where Ω ⊂ R n is a bounded open subset of a Euclidean domain of dimension n. The Hilbertian structure of L 2 (Ω) is induced by the standard inner product, . For a positive integer n ∈ Z + , P n refers to the Cartesian product of n copies of the space P ; that is, P n = n i=1 P . We use P * to denote the space of linear continuous functionals on the Banach space P . The bilinear form that induces the duality between P and P * is given by ·, · P * ,P : P * × P → R, where x, y P * ,P = x(y) for each x ∈ P * and each y ∈ P .
We define the Sobolev space H 1 (Ω) = z ∈ L 2 (Ω) : ∂z ∂x1 , ∂z ∂x2 ∈ L 2 (Ω) . Here, the spatial derivatives are to be understood as weak derivatives defined in the distributional sense. See [13] for this notion of a derivative of a function that is not necessarily differentiable in the classical sense. We equip the space We say that a sequence f n is weakly converging to f in X, consists of all continuous functions u : We will also use the trace operator, τ : H 1 (Ω) → L 2 (∂Ω), defined as the unique linear bounded operator that satisfies (τ u)(x) = u(x) for every x ∈ ∂Ω whenever u is in C 1 (Ω) (that is, u is at least once differentiable in the classical sense). Informally speaking, the trace operator gives meaning to the evaluation of an element of H 1 (Ω) on ∂Ω, which is a set of measure zero.

Optimization Problem for the Mapping Assignment
This section formulates the mapping assignment as an optimization problem that identifies the unknown coefficient H Γ in the mapping macroscopic model (5). This coefficient is reconstructed from the cumulative robot sensor dataĝ(t), the total number of observations made over the entire domain by all the robots up until time t. We recall that the robots record these observations stochastically according to the chemical reaction (1) and that the agents' spatial states evolve according to equation (4). Here, the velocity input v(t) is predetermined by the user and should be chosen such that the agents approximately cover the entire domain. An example of such a choice for v(t) is one that would drive the robots along a lawnmower path, as shown in Fig. 1 in Section 7.
In the mapping macroscopic model (5), the cumulative robot sensor data is interpreted as a continuous quantity given by an integral over the domain, g(t) = Ω y 2 (x, t)dx. We address the existence and uniqueness of the solutions of the macroscopic model (5) in Section 6, since this model is a special case of the coverage macroscopic model (8) for which k f = 0, k(t) = k o , and the stationary robot state is absent.
We defineĤ as an estimate of H Γ . Due to the one-sided coupling between the state variables y 1 and y 2 in the PDE model (5), the variable y 1 is not affected by the coefficient H Γ . Hence, assuming that the PDE is well-posed and that a unique solution exists in H 1 (Ω) × L 2 (Ω), we can pose the mapping problem follows: Problem 2 (Mapping Problem) Given an initial condition y 0 (x) of PDE model (5) and the corresponding solution y 1 (x, t), find an estimateĤ ∈ S ad , where that satisfies the following equation: Ideally, the estimateĤ would be constrained to be a function that takes values of 0 or 1 over its domain. However, this constraint would make the resulting inverse problem, i.e. solving equation (11) forĤ in Problem 2, non-convex when considered as an optimization problem. To make the optimization problem tractable, we relax the range ofĤ and consider it as an essentially bounded element of L 2 (Ω). Generally, equations of the form (11) need not have unique solutions unless some special conditions on k o y 1 (x, t), the kernel of the integral operator K : L 2 (Ω) → L 2 (0, T ), can be guaranteed. To resolve the ill-posedness of the inverse problem, it can be alternatively posed as an optimization problem with a functional J(Ĥ) that is convex, but not necessarily strictly convex: To ensure a unique solution forĤ, the optimization problem can be recast with a strictly convex functional: where λ > 0 is the regularization parameter that is often used in the so-called "Tikohnov regularization" of inverse problems [25]. Since S ad is convex in L 2 (Ω), the existence and uniqueness of regularized problems of the form (13) can be guaranteed and has been been well-studied in the theory of inverse problems [25]. Hence, we have the following result in Theorem 3. The assumption in the theorem that y 1 ∈ C([0, T ]; L 2 (Ω)) is justified in Section 6.2.
We use a gradient descent method to solve the optimization problem (13). To apply this method, we need to characterize the derivative of the objective functional J λ (Ĥ). This functional is differentiable in the Fréchet sense. Since K ∈ L(L 2 (Ω), L 2 (0, T )), the derivative of K is itself. Then by the chain rule of differentiation, the Fréchet derivative of J λ (Ĥ), denoted by J λ (Ĥ), is defined as Using Riesz representation [13][Appendix D, Theorem 2], we can obtain an explicit representation of the gradient of J λ (Ĥ), where K * ∈ L(L 2 (0, T ), L 2 (Ω)) is defined as To verify that this characterization of K * is correct, it is straightforward to check that ∀Ĥ ∈ L 2 (Ω), ∀G ∈ L 2 (0, T ).

Optimal Control Problem for the Coverage Assignment
In Section 6.1, we formulate the coverage assignment as an optimal control problem that is used to compute the parameters v(t) and k(t) in the coverage macroscopic model (8). We analyze the existence and uniqueness of solutions to the PDE model (8) in Section 6.2, the well-posedness of the optimal control in Section 6.3, and the differentiability of the objective functional in Section 6.4.

Formulation of the Optimal Control Problem
We first express the PDE model (8) as a bilinear control system with control inputs v(t) and k(t). Toward this end, we supply the following definitions.
Recall from Section 4 that V = H 1 (Ω), X = V ×L 2 (Ω) 2 , and X * := V * ×L 2 (Ω) 2 , and that L(X, L 2 (Ω) 3 ) denotes the space of linear bounded operators from X to L 2 (Ω) 3 . We define the operators A, B i ∈ L(X, L 2 (Ω) 3 ), i = 1, 2, 3, as follows: We also define the spaces F = L 2 (0, T ; L 2 (Ω) 3 ), G = L 2 (0, T ; L 2 (∂Ω)), Y = L 2 (0, T ; X), and Y * = L 2 (0, T ; X * ). These definitions will be used to establish the well-posedeness of the PDE (8) using the classical theory of weak solutions of linear PDEs [13]. For f ∈ F , g ∈ G, and y 0 ∈ L 2 (Ω) 3 , let y ∈ Y be a function with time derivative ∂y/∂t ∈ Y * . Denoting the components of the robots' velocity field by v(t) = [v x (t) v y (t)] T , we define the control inputs u 1 (t) = v x (t), u 2 (t) = v y (t), and u 3 (t) = k(t). We can now write the PDE (8) with boundary condition (6) and initial condition (9) as a bilinear control system in the following form: where Q and Σ are as defined in Section 3.2. We note that the solution of the PDE model in equations (8), (6), (9) is the solution of the system (18) with f = 0 and g = 0. We consider the more general form (18) for the purpose of analyzing the differentiability properties of the control-to-state map (defined in Section 6.4) and the objective functional in the optimal control problem. A function y is a weak solution of system (18) if for all φ ∈ L 2 (0, T ; X). Here, A g (t) : X → X * is the variational form of the operator A for each t ∈ (0, T ). The boundary condition (6) is equipped with A g in the variational formulation using Green's theorem as, where M g (t) : V → V * for each t ∈ (0, T ) is the Laplacian in the variational formulation and is defined as for all y ∈ V and φ ∈ V * . The coverage problem can now be framed as follows.
Problem 4 (Coverage Problem) Define a target spatial distribution of robot coverage activity, y Ω ∈ L 2 (Ω) 3 , to be achieved by time T , and a set of admissible control inputs, , and define a function W ∈ L(L 2 (Ω) 3 ) that weights the relative significance of minimizing the distances between different state variables and their target distributions. Then, given an initial condition y 0 (x), solve the optimal control problem min (y,u)∈Ȳ ×U ad subject to equation (18) with f = 0 and g = 0.
Note that, due to the essential bounds on u, we have that u ∈ L ∞ (0, T ) 3 . Here, we only consider essentially bounded control inputs since we only prove existence of solutions for control variables in L ∞ (0, T ) 3 . The existence and uniqueness of solutions for control variables in L p (0, T ) 3 , p < ∞, is a nontrivial issue and is beyond the scope of this work.

Energy Estimates
Energy estimates refer to bounds on the solutions of a PDE system with respect to certain parameters of interest, such as the initial condition, coefficients, or boundary parameters. Whereas in the theory of weak solutions of PDEs [13], these energy estimates are used to show the existence of solutions, in optimal control analysis they are used to study the differentiability properties of the control-to-state map. In this section, we derive energy estimates for the solutions of the PDE model (19) and use them to establish existence and uniqueness of solutions to the model, a result (Lemma 7) that we previously stated without proof in our work [12]. More importantly, the estimates are then used to prove existence of solutions to the optimal control problem. The differentiability of the control-to-state map also follows from these energy estimates.
Then we have the following energy estimate for all y ∈ V : which holds for someα,β > 0 that depend only on Ω, max Proof. Setting φ = y in definition (21), we obtain the inequalities from which we can conclude that where τ is the trace operator defined in Section 4. Then it follows from a modified form of the Trace Theorem for domains with Lipschitz boundaries [19][Theorem 1.5.1.10] that there exists a constant K, independent of ∈ (0, 1) and y ∈ V , such that the following inequality holds for every ∈ (0, 1): Thus, for each ∈ (0, 1), the inequality (25) becomes Note that we can choose an > 0 arbitrarily small, independent of y ∈ V , such that D >b . We fix such an > 0. Then we setβ = D −b andα = max{ 1 2 , C }. This concludes the proof.

Lemma 7
Given f ∈ L 2 (0, T ; L 2 (Ω) 3 ), g ∈ L 2 (0, T ; L 2 (∂Ω)), and an initial condition y 0 ∈ L 2 (Ω) 3 , there exists a unique solution y ∈ C([0, T ]; L 2 (Ω) 3 ) to system (18). We have the following energy estimate for this solution: where K > 0 depends only on Ω, max Proof. We first determine a bound on the sum y C([0,T ];L 2 (Ω) 3 ) + y L 2 (0,T ;X) in the energy estimate (28). Let y be a weak solution of system (19), and set φ = y in (19). Then Combining this equality with the inequality (27), and observing that ∂y Using Cauchy's inequality and Young's inequality [13], we derive the following inequality: where r is the maximum of the operator norms of the operators B i , i = 1, 2, 3. Then we have a bound on the first term on the right side of inequality (29) as follows: Applying Young's inequality, we can bound the term on the right side of inequality (30): Combining these upper bounds with the inequality (29), we find that: and therefore d dt for C = E 2 β + 1 + 2α. This inequality implies that Setting , and applying Gronwall's lemma [13] with these functions, we get: Then integrating both sides of inequality (32) over the time interval (0, T ), and using the inequality above, we obtain a bound on the sum of the first two terms in the energy estimate (28): where C = max{e CT , C β }. Next, we derive a bound on the third term in estimate (28), ∂y/∂t L 2 (0,T ;X * ) . We know that B i is a bounded operator from X to L 2 (Ω) 3 for each i ∈ {1, 2, 3}. Therefore, B i q L 2 (Ω) 2 ≤ c i q X for all q ∈ X and for some positive constants c i . Similarly, A g ∈ L(X, X * ) when g = 0. It follows that A g q 2 X ≤ c A ( q 2 X + g 2 L 2 (∂Ω) ) for all q ∈ X and for a fixed g ∈ L 2 (∂Ω) and some positive constant c A . Let v ∈ X such that v X ≤ 1. Then, using inequality (33) and the constants c A and c i , i = 1, 2, 3, we find that any weak solution y of system (18) satisfies the following inequality, for almost every t ∈ (0, T ) and for the constant Here, we are bounding the left-hand side of equation (19) by setting φ(t) = v for almost every t ∈ (0, T ).
The inequality above implies that The bound on the term ∂y/∂t L 2 (0,T ;X * ) follows from the above inequality. Hence, the bounds in the energy estimate (28) hold for K = max{2C , 2C }.
To confirm the existence and uniqueness of the solution to the PDE (8), we note that the weak form of the PDE (19) is in the abstract form dy dt =Ã(t)y +f (t), Here, the term ∂Ω g(x, t)ψ 1 (x)dx defines a bounded functional on X due to the bounds on the trace map τ

Existence of Optimal Control
In this section, we prove the existence of a solution to the optimal control problem (22). Toward this end, we define the control-to-state mapping, Ξ: U ad →Ȳ , which maps a control, u, to y, the corresponding solution to system (18) for f = 0 and g = 0. We introduce the following reduced optimization problem, min u∈U adĴ where J is defined in problem (22). This problem incorporates the PDE system (18) into the objective functional, rather than defining it as a set of constraints as in the original problem formulation (22). We shall henceforth analyze the reduced problem (35).
Theorem 8 An optimal control u * exists that minimizes the objective func-tionalĴ in problem (35).
Proof. The functionalĴ(u) is bounded from below, which implies that the infimum is a finite real number. Therefore, q = inf u∈U adĴ (u) exists. We now determine an optimal pair (y * , u * ), for which J(y * , u * ) = q. Let {u n } ∞ n=1 be a minimizing sequence such thatĴ(u n ) → q as n → ∞. U ad is a bounded and closed convex set, and thus is weak sequentially compact. Hence, there exists a subsequence {u n } ∞ n=1 such that with u * ∈ U ad . Recall from Section 4 that X = V × L 2 (Ω) 2 . Then, the uniform boundedness of the solution, y, in the energy estimate (28) allows us to extract a subsequence {y n } ∞ n=1 , where y n = Ξ(u n ), such that y n y * in L 2 (0, T ; X).
It is necessary to confirm that Ξ(u * ) = y * , since we do not know whether Ξ is weakly continuous. Using the energy estimate (28), we have uniform bounds on y n L 2 (0,T ;X) + ∂y n /∂t L 2 (0,T ;X * ) . Therefore, by the Aubin-Lions lemma [38], we have that there exists a subsequence such that y n 1 → y * 1 in L 2 (0, T ; L 2 (Ω)).
The energy estimate (28) implies that the terms ∇y n 1 , ∂y n ∂t , and k f y n 2 are uniformly bounded. Hence, we can also conclude that there exist subsequences that satisfy: in L 2 (0, T ; L 2 (Ω)).
The first two components of u n are denoted by v n and the third component by k n . From the strong convergence of y n 1 in L 2 (0, T ; L 2 (Ω)) and the weak convergence of u n in L 2 (0, T ) 3 , we can further deduce that: To prove convergence of the terms in the weak form of the PDE (19) that arise from the boundary condition (6), we apply Green's theorem v n · ∇y n 1 , φ L 2 (0,T ;L 2 (Ω)) for all φ ∈ L 2 (0, T ; V ). Due to the strong convergence of y n 1 in L 2 (0, T ; L 2 (Ω)) and the weak convergence of v n in L 2 (0, T ) 2 , we have that v n · y n 1 v * · y * 1 in L 2 (0, T ; L 2 (Ω)).
The convergence of all the sequences defined thus far implies that the sequence of solutions y n = Ξ(u n ), which satisfies with g = 0, converges to the solution y * = Ξ(u * ), which satisfies with g = 0.
It remains to be shown thatĴ(u * ) = q. Because J is weakly lower semicontinuous, we can state that q = lim n→∞ J(y n , u n ) ≤ J(y * , u * ).
Since q was defined earlier as the infimum ofĴ(u), we conclude that which completes the proof.

Differentiability of the Control-to-State Map
In this section, we summarize several results from our prior work [12] that we apply to derive the gradient ofĴ(u). We use the gradient in the method of projected gradient descent to numerically compute a locally optimal control. The numerical computation of this gradient is made possible by deriving an expression for the gradient in terms of the adjoint equation, which is a system of PDEs. This method for computing the control variables is explained in [10][Section 5.2]. See also [40][Section 3.7.1], where the method is discussed in the more general context of optimal control of parabolic PDEs.

Proposition 9
The mapping Ξ is directionally differentiable at every u ∈ U ad along each direction h ∈ L ∞ (0, T ) 3 . Its directional derivative Ξ (u) : U ad →Ȳ in the direction h ∈ L ∞ (0, T ) 3 , denoted by Ξ (u)h, is given by the solution w to the following equation, Theorem 10 The reduced objective functionalĴ is directionally differentiable at every u ∈ U ad along all h ∈ L ∞ (0, T ) 3 . Its directional derivative along h ∈ L ∞ (0, T ) 3 has the form where p is the solution of the backward-in-time adjoint equation, with the Neumann boundary conditions and final time condition where W is defined in Problem 4.

Remark 11
The adjoint equation (47)-(49) should not be confused with the backward heat equation, which is known to be ill-posed in H 1 (Ω). The wellposedness of system (47)-(49) can be established using a standard variable transformation of the form p * (t) = p(T − t). The resulting equation, with the Neumann boundary conditions and initial condition can be shown to have a unique solution p * using arguments similar to those in the proof of Lemma 7.

Simulation Results and Discussion
In this section, we validate our approaches to the mapping assignment (Section 5) and coverage assignment (Section 6) with numerical simulations. The PDE models of ensemble mapping and coverage, presented in Section 3.2, are numerically solved using the method of lines. In this method, the operators associated with the PDE are initially discretized in space. For this spatial discretization, we use a finite-difference and flux limiter based scheme on a uniform mesh. Time is viewed as a continuous variable, and the resulting system of ordinary differential equations (ODEs) is solved using built-in ODE solvers in MATLAB. See [23] for further details on this approach to numerically solving advectiondiffusion-reaction type PDEs.

Mapping Assignment
We validate our mapping approach in two test cases. In both scenarios, the domain is Ω = [0, 100] 2 m 2 , the ensemble has N = 30 agents, the diffusion coefficient is D = 10 −4 m 2 /s, and the agents' rate of recording observations over a region of interest is k o = 100 s −1 , yielding a high probability rate of registering observations. The agents' velocity field v(t) is defined to drive the agents along a lawnmower path, as illustrated by the two agent trajectories in Figure 1. This path is perturbed by stochastic fluctuations arising from the agents' diffusive motion. We define two environments, Case 1 and Case 2, with different regions of interest, labeled Actual H Γ (x) in Figure 2. The results of our mapping approach are shown in the plots labeled Estimated H Γ (x). To use the estimated maps in simulations of our coverage approach (Section 7.2), we thresholded the estimated H Γ (x) in both environments by defining the Thresholded H Γ (x), also shown in Figure 2, as H T (x) = 1 for each x ∈ Ω such that H Γ (x) ≥ 0.5 and H T (x) = 0 otherwise. In both cases, our optimization method is able to reconstruct the spatial coefficient H Γ (x) with considerable accuracy, even though a relatively small number of robots was used (N = 30). The largest error in the estimates occurs in the top half of the Case 2 environment, which can be attributed to the increased dispersion of the agents, due to their diffusive motion, as they reach the upper portion of the domain. We would expect a larger ensemble of agents to generate a more accurate map, since the microscopic model converges to the macroscopic model as N → ∞ [43].

Coverage Assignment
We validate our coverage approach in three test cases, Cases 1, 2, and 3. The region of interest in Case 1 is defined as Thresholded H Γ (x) in Figure 2(a), and the region of interest in Cases 2 and 3 is defined as Thresholded H Γ (x) in Figure 2(b). In all scenarios, the domain is Ω = [0, 100] 2 m 2 , and the diffusion coefficient is D = 5 × 10 −4 m 2 /s. To illustrate the scalability of our control methodology to large numbers of agents, we simulate an ensemble with N = 1000 agents. All agents are initialized in the moving state, with an initial density given by a Gaussian distribution centered at [x o1 x o2 ] = [ 10 10]. Thus, the initial conditions (9) are: , y 2 (x, 0) = 0, with σ x1 = σ x2 = 0.02 and A defined such that the Gaussian distribution integrates to 1 over the domain. The final time is set to T = 800 s for Cases 1 and 2 and T = 300 s for Case 3.
We define y Ω , the target spatial distribution in Problem 4, as follows. We partition the domain Ω into P = 20 cells, denoted by {Ω nm }. The cell Ω 00 Actual H Γ (x) occupies the region [0, 100 P ] × [0, 100 P ], and all other cells are defined as Ω nm = 100 m P , 100 m+1 P × (100 n P , 100 n+1 P ], where n, m ∈ {0, 1, ..., P − 1} and n = 0 or m = 0. For Case 1 and Case 2, we set the target number of instances of desired robot activity in each cell Ω mn to be z * mn = C × µ(Ωnm ∩ Γ) µ(Ωmn) , where C is a positive constant and µ is the Lebesgue measure on Ω. The target distribution of robot coverage activity, which is the third component of y Ω , is defined as y * 3 (x, T ) = C 50 for all x ∈ Γ. Since we do not require the first two components of y Ω (the densities of moving and stationary robots) to reach target distributions at time T , we set the function W in Problem 4 to be a diagonal operator matrix of the form diag([0 0 I]), where I is the identity operator on L 2 (Ω). We specify the following sub-cases: C = 450 in Case 1a and Case 2a, and C = 3600 in Case 1b and Case 2b. For Case 3, coverage activity is desired only in the upper half of the domain; we set y * 3 (x, T ) = 36 for x ∈ Γ with x 2 ≥ 60, and y * 3 (x, T ) = 0 otherwise. Figure 3 plots the time evolution of the objective function J for each case. For Cases 1a,b and 2a,b, the low values of J at the final time T = 800 s indicate that the target coverage density is nearly achieved. This is in part due to the accuracy of the thresholded estimate of H Γ (x) obtained from the mapping task (Section 7.1). The lower error in the estimated H Γ (x) for Case 1 compared to this mapping error for Case 2 contributes to the better coverage performance (i.e., lower values of J) in Case 1a versus Case 2a and in Case 1b versus Case 2b. For Case 3, the value of J is higher at the final time T = 300 s than the values of J for Cases 1a,b and 2a,b at T = 800 s, indicating that the final coverage density is relatively farther from the target density. The poorer coverage performance in Case 3, in which coverage activity is limited to a subset of the region of interest, is a consequence of the limited controllability of the system, which can be attributed to three factors. First, only three control variables are used to control the PDE model, which is an infinite-dimensional dynamical system. Second, the system is constrained to achieve the target coverage density by time T . Third, the time-dependent diffusion, reaction, and advection operators in the PDE model commute [26], and this commutativity of the control and drift vector fields degrades the controllability properties of the system [2]. Thus, the ensemble would achieve better coverage performance in assignments with less stringent requirements on the system controllability properties. Such assignments include those in which the target coverage distribution is proportional to an environmental parameter, or in which the objective is to achieve a minimum coverage density in each region of interest rather than an exact coverage density. Figure 4 plots the target coverage density z * mn for Cases 1b and 2a alongside the corresponding expected coverage density, y 3 (x, T ) from the macroscopic PDE model (18), and the achieved coverage density, z Ωmn of agents (N = 1000) in the ensemble. As demonstrated in [43], the discrepancy between these two density fields will tend to zero as N → ∞.

Conclusion
In this paper, we presented stochastic approaches to mapping and coverage assignments for an ensemble of autonomous robots in a PDE control framework. This framework enables the modeling and control of a robotic ensemble in a rigorous way that is scalable with the number of robots. We demonstrated that temporal data obtained by a small ensemble of diffusive agents can provide rich information about the spatial distribution of a region of interest, despite severe restrictions on the agents' sensing, localization, tracking, and computational capabilities. We also showed that we can pose the coverage task as an optimal control problem that computes the agents' control inputs to achieve a target distribution of coverage activity over the previously mapped regions of interest.
In future work, we will investigate the controllability properties of the PDE models that we have presented in this paper. Additionally, we plan to incorporate pairwise interactions between agents, such as those defined by attractionrepulsion potentials, in order to increase the cohesiveness of the ensemble and improve the reachability properties of the system. The component of the robots' velocity field that is induced by pairwise interactions would be included in the advection term, resulting in a nonlinear PDE as the macroscopic model. This type of model would require more advanced analytical tools than the ones used in this paper.