{'fraction_non_alphanumeric': 0.045555843482767554, 'fraction_numerical': 0.012127494169473957, 'mean_word_length': 4.320099255583126, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Indirect Optimal Control of Advection-Diffusion Fields through Robotic Swarms

In this paper, we consider the problem of optimally guiding a large-scale swarm of underwater vehicles that is tasked with the indirect control of an advection-diffusion environmental field. The microscopic vehicle dynamics are governed by a stochastic differential equation with drift. The drift terms model the self-propelled velocity of the vehicle and the velocity field of the currents. In the mean-field setting, the macroscopic vehicle dynamics are governed by a Kolmogorov forward equation in the form of a linear parabolic advection-diffusion equation. The environmental field is governed by an advection-diffusion equation in which the advection term is defined by the fluid velocity field. The vehicles are equipped with on-board actuators that enable the swarm to act as a distributed source in the environmental field, modulated by a scalar control parameter that determines the local source intensity. In this setting, we formulate an optimal control problem to compute the vehicle velocity and actuator intensity fields that drive the environmental field to a desired distribution within a specified amount of time. In other words, we design optimal vector and scalar actuation fields to indirectly control the environmental field through a distributed source, produced by the swarm. After proving an existence result for the solution of the optimal control problem, we discretize and solve the problem using the Finite Element Method (FEM). The FEM discretization naturally provides an operator that represents the bilinear way in which the controls enter into the dynamics of the vehicle swarm and the environmental field. Finally, we show through numerical simulations the effectiveness of our control strategy in regulating the environmental field to zero or to a desired distribution in the presence of a double-gyre flow field.


INTRODUCTION
Given the exponential decrease in cost of electronic components over the last few decades, large collectives of robots, or robotic swarms, are becoming a viable option for a variety of missions of ever-increasing complexity, such as coverage, mapping, search-and-rescue, and surveillance (Dorigo et al., 2021). Controllers for robotic swarms should satisfy mission requirements while scaling gracefully as the number of robots N increases. Mean-field models of robotic swarms (see, e.g., Elamvazhuthi and Berman (2020)) describe a swarm as a set of probability densities over space and time; since these models are independent of N , they can be used to design controllers for arbitrarily large swarms, with the caveat that the distribution of the swarm is controlled rather than individual robots.
In this paper, we take advantage of the mean-field model's invariance to swarm size by using such a model to design scalable robotic swarm controllers that achieve indirect control of a distributed process that evolves according to a Partial Differential Equation (PDE), such as the concentration field of a contaminant in a fluid flow. The robots are underwater vehicles that are each equipped with an actuator that acts as a source for the process, and we aim to indirectly control the process through the coordinated motion of the swarm and the source actuation. A similar problem has previously been considered for finite teams of mobile robots, in which individual robot trajectories are controlled. For example, Cheng and Paley (2021) present an optimal control approach that uses an operator-valued Riccati equation to formulate the optimal actuation as Fig. 1. Illustration of the indirect control problem. The environmental field S diffuses from the segment of the boundary Γ d under the advection of a fluid velocity field F. The vehicle velocity field u and actuator intensity field k are designed to drive S to a target distribution. a function of the optimal guidance, and then recast the problem in terms of the latter alone to jointly optimize the guidance of the robots and their associated actuation. Demetriou (2021) describes a path-dependent reachability approach that takes into account constraints on the robots' motion and real-time implementation while regulating a spatially distributed process using local decentralized measurements only.
We define the indirect control problem for a distributed robotic swarm whose mean-field dynamics are modeled by a Kolmogorov forward equation in the form of a linear parabolic advection-diffusion PDE. We formulate an optimal control problem (OCP) for this mean-field model, which is coupled with the advection-diffusion dynamics of the environmental field and prove an existence theorem for the OCP using techniques for optimal control of PDEs. Then, we derive a set of first-order necessary optimality conditions and solve them numerically using a Finite Element Method (FEM) discretization. Finally, we evaluate the effectiveness of our control strategy through numerical simulations of regulation and target tracking problems in the presence of a double-gyre flow field.

PROBLEM FORMULATION
We consider a swarm of robots, labeled i = 1, ..., N , that move in a bounded fluid domain Ω ∈ R 2 . Robot i occupies position X i (t) ∈ Ω at time t and moves with velocity v(x, t) ∈ R 2 , which is the sum of its self-propelled velocity u(x, t) and the fluid velocity at its position, F(x). This motion is perturbed by a 2-dimensional Wiener process W(t), which models stochasticity arising from inherent sensor and actuator noise or intentionally programmed "diffusive" exploratory behaviors. The robot's position evolves according to the following Stochastic Differential Equation: where D q > 0 is a diffusion coefficient, n(x) is the unit normal to the domain boundary ∂Ω at x, and ψ(t) ∈ R is a reflecting function, which ensures that the swarm does not exit the domain. We assume that each robot i carries an on-board actuator that acts as a source of intensity k i (t) ∈ R in a scalar environmental field S(x, t), x ∈ Ω. We also assume that part of the boundary, Γ d ∈ ∂Ω, acts as a source with constant intensity S d for simplicity. The field S is an advection-diffusion process with diffusion coefficient D S > 0, modeled by the following PDE problem: where 1 Γ d is an indicator function. The right-hand side of the PDE consists of the cumulative effect of the point sources, which have the same dynamics as the robots since they are installed on-board. Note that the explicit dependence of S on space and time is omitted when clear from the context. In the limit as N → ∞, we obtain a mean-field model that describes the evolution of the probability density q(x, t) of a single robot occupying position x at time t, or alternatively, the swarm density at this position and time. For the robot dynamics we consider here, this model takes the form of a linear parabolic advection-diffusion problem with no-flux boundary conditions; see, e.g., Sinigaglia et al. (2021). The coupled system dynamics are therefore: on Ω × {0}. (1) Note that the right-hand side of the S dynamics is now the product of the intensity of the distributed actuation k and the local swarm density q. As a consequence of this Eulerian perspective, k = k(x, t) varies in both space and time. The robotic swarm thus aims at controlling the concentration S of an advection-diffusion field through a localized actuation field kq that depends on both the actuator intensity k and the swarm density q at each spatial location. See Figure 1 for an illustration of the problem.
Defining the state variables as q and S, the state dynamics consists of the one-way coupled system of PDEs (1). The objective of the control problem is then to find the optimal actuation for u(x, t) and k(x, t) to guide the space-time evolution of the field S to a target distribution z, which may be static or dynamic, at a given time T . This control objective can be easily encoded in a standard quadratic cost functional of the form: where α T , α, β, γ ≥ 0 are weighting constants. In the regulation problem, for example, we set z = 0 and seek the optimal trade-off between effectively regulating the field S and minimizing the overall energy expended by the swarm for propulsion motion and source actuation.

THE OPTIMAL CONTROL PROBLEM
In this section, we prove the existence of optimal controls, derive a system of first-order necessary optimality conditions using the Lagrangian method, and provide a consistent discretization of the OCP using the FEM.

Analysis
Both the q and S dynamics satisfy rather standard advection-diffusion equations of linear parabolic type; see, e.g., (Manzoni et al., 2021, Chapter 7). The natural functional space for the swarm density function q which is subjected to zero-flux Neumann boundary conditions is q ∈ L 2 (0, T, H 1 (Ω)), while the Dirichlet boundary suggests the choice ofS ∈ L 2 (0, T, H 1 0 (Ω)) for the "lifted" field variableS = S −S d , whereS d is a suitable extension of the boundary datum to the domain Ω -see, e.g., (Salsa, 2016, Chapter 8). It is also standard to select ∂q ∂t ∈ L 2 (0, T, H 1 (Ω) * ) and ∂S ∂t ∈ L 2 (0, T, H 1 0 (Ω) * ) so that the functional space for q is actu- As done in, e.g., Roy et al. (2018) and Sinigaglia et al. (2021) for similar problems involving the Kolmogorov forward equation alone, we consider L ∞ spaces for the control fields for which energy-like inequalities are readily available; that is, we select U = L 2 (0, T ; L ∞ (Ω) 2 ) × L 2 (0, T ; L ∞ (Ω)) as the control space, so that v = (u, k) ∈ U. Besides the choice of the functional spaces for states and controls, we make the following standard assumptions: The weak formulation of the PDE problem governing the swarm dynamics is: find q ∈ L 2 (0, T ; H 1 (Ω)) such that for a.e. t ∈ (0, T ), The weak formulation of the PDE problem for the "lifted" variableS is: findS ∈ L 2 (0, T ; H 1 0 (Ω)) such that for a.e. t ∈ (0, T ), We also define the linear functional F ∈ H −1 = H 1 0 (Ω) * as F φ = Ω k q φ dΩ − a S (S d , φ). In the following, we will need a bound on the operator norm of F , which we prove in the lemma below. Lemma 1. (Bound on F ). Let assumptions (A1), (A2), (A3), and (A4) hold. Then the following bound on the norm of F holds: Proof. From the definition of F and the Cauchy-Schwarz and Poincaré inequalities, we have: Regrouping the terms and using the definition of operator norm in H −1 (Ω), the result follows. 2 Existence and well-posedness of the state dynamics follow from the well-posedness of the q dynamics and basic energy estimates on the S dynamics. This is a consequence of the one-way coupling from q to S. Following the same arguments as in Sinigaglia et al. (2021), we have that . To prove the well-posedness of the S dynamics, we note that kq ∈ L 2 (0, T ; L 2 (Ω)) since kq 2 L 2 (0,T ;L 2 (Ω)) ≤ k 2 L 2 (0,T ;L ∞ (Ω)) q 2 L 2 (0,T ;L 2 (Ω)) ≤ C k 2 L 2 (0,T ;L ∞ (Ω)) u 2 L 2 (0,T ;L ∞ (Ω) 2 ) , and the latter quantity in the inequality is bounded by the definition of the control space U. Therefore, S satisfies an advection-diffusion equation with L 2 right-hand side for which existence and uniqueness results are available -see, e.g., (Salsa, 2016, Theorem 9.9).
We define the control-to-state operator as the map (S, q) = Ξ[u, k] which associates to each control function v ∈ U with a corresponding state y ∈ Y. The following result regarding the Fréchet differentiability of the control-tostate operator is also needed to prove existence of optimal controls. Lemma 2. (Differentiability of the control-to-state map). Let assumptions (A1), (A2), (A3), and (A4) hold. Then the control-to-state map (S, q) = Ξ[u, k] is Fréchet differentiable and the directional derivative (z S , z q ) = Ξ [u, k](h, l) at (u, k) ∈ U in the direction (h, l) ∈ U is the solution of the coupled PDE system: Proof. (Sketch) The derivation of the equations governing the sensitivity z q of the swarm dynamics, and bounds on the norm of z q , can be found e.g. in Roy et al. (2018) and Sinigaglia et al. (2021). On the other hand, the expression for the dynamics of z S can be obtained by formally computing the directional derivative, that is, the limit lim s→0 S(u+sh,k+sl)−S(u,k) s . Bounds on z S are obtained by adapting the above results on the norm of S, and noting that by the triangular inequality, kz q + lq L 2 (0,T ;L 2 (Ω)) ≤ kz q L 2 (0,T ;L 2 (Ω)) + lq L 2 (0,T ;L 2 (Ω)) , so that continuity of the sensitivity z q and z S with respect to variations of the control functions h, l can be easily obtained, thus proving the differentiability of the control-to-state map. 2 We are now ready to prove a result concerning the existence of optimal controls. Theorem 3. (Existence of optimal controls). Let assumptions (A1), (A2), (A3), and (A4) hold. Then there exists an optimal controlv = (ū,k) that minimizes the cost functional (2) subject to the dynamics (1), that is, y = (S,q) = Ξ[ū,k].
Therefore, the pair (ȳ,v) is an optimal pair for the considered optimal control problem. 2 We note that uniqueness of an optimal solution is not guaranteed, due to the bilinear way in which both controls u and k enter into the coupled system dynamics.s

Optimality Conditions
We now derive a system of first-order necessary optimality conditions using the Lagrangian multipliers method. For the problem at hand, the Lagrangian can be defined as Note that we have defined adjoint fields λ q and λ S that are related to the state dynamics of both q and S. The adjoint dynamics for λ q and λ S thus satisfy: in Ω × {T }, which are obtained by taking the first variation of the Lagrangian along variations in swarm density q and the environmental field S, respectively. Note that the coupling between the adjoint fields is dual with respect to the state dynamics. The coupling is from q to S in the state system, while it is from λ S to λ q in the adjoint system. The dual of the forcing term k q in the S dynamics is the forcing term kλ S in the λ q dynamics.
The reduced gradients of J with respect to u and k can therefore be expressed as by taking the first variations of the Lagrangian in the directions of u and k, respectively. Note that, despite k entering linearly into the S dynamics, the reduced gradient ∇J k depends on the dynamics of the swarm density q. This is a consequence of the multiplicative nature of the S forcing term k q.

Numerical Discretization
The OCP with coupled system dynamics (1) is discretized in the state variables S and q using the Finite Element Method (FEM). The discretized state dynamics are where B is the rank-3 transport coefficient tensor defined by B ijk = Ω ∂φj ∂x φ i φ k dΩ; Bu is the tensor vector product defined by (Bu) ij = Nu k=1 B ijk u k ; C is the reaction tensor defined by C ijk = Ω φ j φ i φ k dΩ; and M , B, and A are the usual FEM mass, transport, and stiffness matrices, respectively.
The discretization of the adjoint system is: Finally, the reduced gradient discretization is: We can apply the same reasoning as in Sinigaglia et al. (2021) to perform numerical gradient computation. Therefore, we use the Discretize-then-Optimize (DtO) approach (see e.g., (Manzoni et al., 2021, Chapter 8)) to numerically solve the problem while avoiding inconsistencies in the gradient computation. In order to do so, the discrete Lagrangian must be computed and differentiated. This computation, which is very similar to the one in our previous work (see Sinigaglia et al. (2021) for a more detailed treatment of the problem for the swarm dynamics alone), is not reported here due to space constraints.
Since the coupling is one-way, at each time step we advance the dynamics of the swarm density q and then solve the problem for S. Using similar reasoning, we first solve the discrete adjoint dynamics with respect to λ S , and then the adjoint problem for λ q . It can be checked that carrying out the optimization at the continuous level and then discretizing the optimality conditions, that is, adopting the Optimize-then-Discretize (OtD) approach, results in the same system of equations at the semi-discrete level; up to choosing a suitable time-discretization, the two approaches fully commute.

SIMULATION RESULTS
In this section, we present numerical simulation results that show the effectiveness of our control strategy. The computational domain Ω = [0, 1] 2 is discretized into a triangular mesh with N s = 2, 788 degrees of freedom, and the time interval [0, T ] is discretized into N t = 61 time steps, where the final time is T = 1.5(s). The resulting fully discrete optimization problem has 510,204 control variables and 170,068 state variables. A steady double-gyre flow field is chosen as the fluid velocity field F.
Using the DtO method, the reduced gradient is computed with respect to the control variables only. Computations are carried out in MATLAB using a modified version of the redbKit (Quarteroni et al. (2015)) library to assemble the FEM matrices and tensors and the Tensor Toolbox (Bader and Kolda (2006)) to perform efficient tensor computations. The nonlinear optimization software Ipopt (Wächter and Biegler (2006)) is then used to solve the resulting nonlinear optimization problem.
Two test cases are considered. In Test Case 1, we solve a regulation problem with target distribution z = 0 and the Dirichlet boundary condition illustrated in Figure 1, with S = S d = 10 along Γ d . Test Case 2 is a tracking problem with z = S 0 (x) and a homogeneous Dirichlet boundary condition. In Figure 2, we compare the controlled and uncontrolled dynamics of the total mass of the environmental field S, defined as m S (t) = Ω S(·, t)dΩ, for the two test cases. For Test Case 1, the uncontrolled steady-state value of m S (t) depends on the equilibrium balance between the mass generated along Γ d and the mass absorbed along the rest of the boundary. In the controlled case, however, the robotic swarm drives m S (t) to zero through coordinated motion, defined by their self-propelled velocity u, and the intensity k of their distributed actuation. In Test Case 2, the uncontrolled mass m S (t) exponentially converges to zero due to diffusion and the homogeneous Dirichlet boundary condition, while the controlled mass m S (t) is driven to Ω S 0 (·)dΩ due to the efforts of the swarm to maintain S at its initial condition S 0 (x). Figure 3 shows snapshots of the swarm density dynamics under the action of the controls and the fluid velocity field for Test Case 1. Figure 4 compares snapshots of the controlled and uncontrolled dynamics of S for Test Case 1, and Figure   5 presents snapshots of the optimal actuation k for this case. Finally, snapshots of the controlled dynamics of S for Test Case 2 are shown in Figure 6. Fig. 2. Controlled and uncontrolled dynamics of the total environmental field mass for both test cases. After a common transient phase during which the robotic swarm moves toward where its actuation is most effective, the controlled dynamics are driven to the total mass of the target distribution z. Fig. 3. Test Case 1. Top: Swarm density contours (black), optimal robot velocity control u(x, t) (blue), and fluid velocity field F (red) at times t = 0.25(s) and t = 1.5(s). Bottom: Euclidean norm ||u|| of the optimal control, with a maximum value of √ 2.