Controlled transport of fluid particles by microrotors in a Stokes flow using linear transfer operators

The manipulation of a collection of fluid particles in a low Reynolds number environment has several important applications. As we demonstrate in this paper, this manipulation problem is related to the scientific question of how fluid flow structures direct Lagrangian transport. We investigate this problem of directing the transport by manipulating the flow, specifically in the Stokes flow context, by controlling the strengths of two rotors fixed in space. We demonstrate a novel dynamical systems approach for this problem and apply this method to several scenarios of Stokes flow in unbounded and bounded domains. Further, we show that the time-varying flow field produced by the optimal control can be understood in terms of dynamical structures such as coherent sets that define Lagrangian transport. We model the time evolution of the fluid particle density using finite dimensional approximations of the Liouville operators for the micro-rotor flow fields. Using these operators, the particle transport problem is framed as an optimal control problem, which we solve numerically. This framework is then applied to the problem of transporting a blob of fluid particles in domains with different boundary conditions: free space, near to a plane wall, in a circular confinement, and the transport of two distributions of particles to a common target. These examples demonstrate the effectiveness of the proposed framework and also shed light on the effects of boundaries on the ability to achieve a desired fluid transport using a rotor-driven flow.

pair of micro-rotors by modulating the torques of these rotors.
In recent decades, significant attention has been given to the application of dynamical systems theory to problems of fluid transport at low Reynolds number, with much of this attention focusing on mixing by chaotic advection [18][19][20] .Such research has been largely inspired by applications including industrial mixing, design of microfluidic lab-on-a-chip devices, and biomedical applications such as cell sorting and targeted drug delivery.While mixing is important to many such applications, many also require an ability to transport packets of fluid or a concentrated, passive scalar in a controlled way to a target destination while minimally mixing, stretching, or distributing the blob.Despite its growing practical importance, this area has received considerably less research attention.In particular, the application considered in this work is motivated by recent works which have shown that biological and artificial microswimmers can generate significant transport of passive fluid particles [21][22][23][24] and finite-sized particles in the fluid [25][26][27] , and the potential usefulness of this observation for the design of micro-scale robots for controlled fluid or cargo transport 17,[28][29][30][31] .
In this work, we study the problem of steering an ensemble of fluid particles in a Stokes flow from an initial particle distribution to a final distribution, where the particles are advected by the flow field generated by a pair of fixed rotors.In our formulation, the distribution of fluid particles is described by a density function, and a data-driven method based on a finite dimensional approximation of the Liouville operator associated with the rotor-driven flow is developed to approximate the density transport dynamics.With this model, we show that the problem of controlled density transport can be posed as an optimal control problem which we solve using differential dynamic programming, an iterative trajectory optimization scheme.To apply this framework to the problem of steering a density using fixed rotors, we model the rotors as rotlets, the singularity solution of the Stokes equations associated with a point torque 32,33 .This work is an extension of the authors' recent conference paper 34 .In this work, we seek to further highlight the fluid mechanical applications of the proposed method and use it to study the effects of boundaries on this fluid transport problem.Further, we show that the Perron-Frobenius operator associated with the time-varying flow field produced by the optimal control can be used to analyze the coherent sets and coherent structures of this flow.This enables us to make a useful qualititative connection between the optimal control and the dynamical structures which determine the Lagrangian transport, an area which has received increased interest in recent years 35,36 .
The rotlet singularity model has become commonly used as an approximation for flows generated by rotating bodies at small length and velocity scales.Meleshko and Aref 37 studied the flow generated by the so-called blinking rotlet model consisting of two rotlets at fixed positions in a circular domain which run for a fixed time period in an alternating pattern, itself a Stokes flow alternative to the blinking vortex model introduced by Aref 18 .These models have received interest as minimalistic examples of the concept of chaotic advection, the notion that fluid particle trajectories in time dependent laminar flows can exhibit chaotic motions, even in two spatial dimensions.Van der Woude et al. 38 considered a similar blinking rotlet problem in a rectangular cavity and considered mixing by sinusoidal stirring patterns as well as the typical blinking pattern.While these works introduce a time dependence by explicitly varying the strengths of fixed rotors in time, more recent works have studied effect of a time dependence in the fluid flow due non-stationary rotors, typically where each rotor is advected by the flow field generated by all other rotors [39][40][41] .Recently, Piedra et al. 42 analyzed the flow produced by an electromagnetically driven rotor both experimentally and theoretically, also linking the mixing and transport properties of the flow to the associated coherent structures.In this work, while we only consider the case of fixed rotors, we develop methods to stir in a controlled way to steer a distribution to a desired location while minimizing the spread of the particle distribution.
Several works have framed fluid mechanical transport problems as optimization or optimal control problems [43][44][45][46][47][48] , with most of these focusing on optimizing mixing performance.Mathew et al 43 studied the problem of optimally modulating (in time) a finite set of spatially varying force fields to optimize mixing over a fixed timespan and for a fixed action integral, using a conjugate gradient descent method to numerically approximate the optimal control.Zhang and Balasuriya 48 develop a method to determine an optimal spatiotemporally varying addi-tive control velocity field for two problems: Lagrangian mixing and to drive trajectories to desired end states in a finite time.In this work, we present a numerical method to optimally modulate two flow fields (corresponding to rotors in fixed positions) in order to drive an initial distribution of fluid particles to a desired final distribution, as specified by the moments of the density function of the distribution.Further, we examine the structure of the optimal flow field by calculating the coherent sets and the associated flow structures produced by this flow field.These results show that the optimal control typically produces a flow field which generates a transport barrier dividing the coherent sets which passes through the blob location at the initial time and connects to the target location at the final time, effectively directing the particle distribution toward the target.
Our method relies on a finite-dimensional approximation of the Liouville operator, the infinitesimal generator of the semi-group of Perron-Frobenius operators 49 , which describe the density transport dynamics for a given flow map.The use of data-driven approximations of transfer operators in modelling fluid flows and in problems with actuation has been an active area of research in recent years 50,51 , with many of the most common methods having their origin in the analysis of fluid flows 52,53 .Refs. 54 and 55 develop a convex optimization formulation based on transfer operators to determine optimal local perturbations of a flow field to enhance mixing of a fluid.Klünker et al 56 recently studied mixing in open flows in terms of spectral properties of a finite-rank approximation of the Perron Frobenius operator.Sinha et al 57 use the Perron Frobenius and Koopman generators associated with a given velocity field to choose an optimal location of release of a dispersant in the flow field.Brockett 58,59 proposed the optimal control of the Liouville equation with applications in ensemble control, but assumes a control input that can be varied arbitrarily in space and time.Relatedly, Grover and Elamvazuthi 60,61 use transfer operators and their generators in a graphbased approach to solving the optimal transport problem, motivated control problems for multi-agent and swarm systems, in which the control is also taken to vary spatiotemporally.The problem considered in this work can be viewed as a variation of those in 60,61 with one significant distinction: in this work the control input u does not vary with the spatial location of the particle.The flow field is restricted to those that can be generated as linear combinations of the flow fields of two fixed micro-rotors and the strengths of the micro rotors in turn influence the flow field.
The remainder of the paper is structured as follows.In Sec.II, we review methods from the operator theoretic view of dynamical systems for modelling the transport of density functions through a dynamical system and present a numerical method for the computation of a finite dimensional approximation of the Liouville operator.In Sec.II C we demonstrate that this method can be naturally extended to account for the effects of actuation on a dynamical system, allowing the use of this framework to express the density transport problem as an optimal reference tracking problem.In Sec.II D we discuss how the operator theoretic methods relate to the computation of finite-time coherent sets for a time-varying flow field.In Sec.III we briefly review the method of differential dynamic programming, an iterative trajectory optimization scheme which we implement to numerically solve this optimal control problem.In Sec.IV, we implement these methods on the problem of steering a density of fluid particles using a pair of fixed microrotors.In Secs.
V and VI, we study the effects of plane wall and circular boundaries on this transport problem in comparison to the case of an unbounded flow.In Sec.VII, we consider the ability to manipulate multiple density functions simultaneously in this system.

II. DENSITY TRANSPORT
In order to formulate the problem of controlling the motion of ensembles of fluid particles, we will first specify the distribution of such an ensemble by a density function.In this section, we will review the methods used to study the evolution of such a density function over time, given that the individual particle motion is specified by a known dynamical system.

Consider a dynamical system
on a measure space (X, A, µ) where x ∈ X is the state, X ⊂ R n is the state space, A is the Borel σ-algebra on X, and µ is a measure on X. Denote the time-t flow map from an initial state x 0 by Φ t (x 0 ).We will further assume that the measure µ is absolutely continuous with respect to the Lebesgue measure, so that µ can be expressed in terms of a density, ρ ∈ L 1 (X), such that dµ(x) = µ(dx) = ρ(x)dx.With this, the Perron Frobenius operator, P t : L 1 (X) → L 1 (X) corresponding to the flow Φ t can be defined as the unique operator 49 such that for any A ∈ A, t ≥ 0. The family of these operators, parameterized by time, t, have been shown to satisfy the properties of a semigroup 49 .The infinitesimal generator of this semigroup, denoted here by L, is known as the Liouville operator or the Perron-Frobenius generator, and defined as where I is the identity operator.Alternatively, as this operator expresses the deformation of a density function under an infinitesimal action of the operator P t , the Liouville operator can be thought of as expressing a continuity equation for the number of particles in the state space 49,62 ; that is, From this definition, we can immediately derive the following important property of the Liouville operator.
Lemma 1 Suppose the Liouville operator associated with a vector field f 1 : X → R n is denoted by L 1 and the Liouville operator associated with the vector field f 2 : X → R n by L 2 , then the Liouville operator associated with the vector field The numerical method used for the computation of the Perron-Frobenius operator and Liouville operator is derived from the relationship between the Perron-Frobenius operator and the Koopman operator.The Koopman operator K t : L ∞ (X) → L ∞ (X) is the operator which propagates observable functions h ∈ L ∞ (X) forward in time along trajectories of the system and is defined as The Koopman and Perron-Frobenius operators are adjoint to one another, with the adjoint relationship given by

B. Numerical approximation
One of the most common methods of approximating the Perron-Frobenius operator is a set-oriented approach known as Ulam's method 63 , in which a domain of interest is discretized into cells, a large number of short-time trajectories are simulated, and then the operator is computed as the matrix containing the approximate transition probabilities between the cells 64 .It has been shown that this method can be viewed as a Galerkin projection of the Perron Frobenius operator onto the function space spanned by indicator functions corresponding to the discrete cells 65 .In recent works involving numerical approximation of the Koopman operator, one of the most common approaches is that of extended dynamic mode decomposition (EDMD) 66 , in which the operator is computed by solving a least squares problem, which can also be viewed as a Galerkin projection of the operator onto a function space spanned by a predefined set of basis functions 65,66 .By exploiting the adjoint relationship between the Perron Frobenius and Koopman operators, it has been shown that methods typically used for one operator can be used to compute the other.Based on this idea, recent works have developed variations of EDMD for the computation of the Perron-Frobenius operator 65,67,68 .In this work, we also implement EDMD for the computation of the Perron-Frobenius operator, which we outline below, largely following Klus et al. 65 .
The method requires a predefined dictionary D of k scalar-valued basis functions, D = {ψ 1 , ψ 2 , . . ., ψ k }, where ψ i : X → R for i = 1, . . ., k and trajectory data collected from the dynamical system with fixed timestep, ∆t, arranged into snapshot matrices as where the subscript i = 1, . . ., m is a measurement index and x + i = Φ ∆t (x i ).Then, given an the observable function h and density ρ, these functions are approximated by their projections onto the space spanned by elements of D as where ĥ, ρ ∈ R k are column vectors containing the projection coefficients and Ψ : X → R k is a columnvector valued function where the elements are given by [Ψ(x)] i = ψ i (x).Substituting these expansions into Eq.( 6) yields Then noting that [K ∆t Ψ](x) = Ψ(x + ) and assuming that P ∆t can be approximated by a matrix P operating on the coordinates ρ, it is clear that in the limit of a large dataset m → ∞, the above expression becomes where e is a residual error arising due to the matrix approximation of P ∆t by P .This can be posed as a leastsquares problem for the matrix where Ψ X ,Ψ Y ∈ R k×m are matrices with columns containing Ψ evaluated on the columns of X and Y respectively.The analytical solution of this least squares problem is where (•) † is the Moore-Penrose pseudoinverse.Given this matrix approximation of the operator, P , if the timestep ∆t chosen in the data collection is sufficiently small, the corresponding matrix approximation L of the Liouville operator can be approximated based on the limit definition of the generator in Eq. 3. as where I k is the k ×k identity matrix.The matrix approximation P of the operator P ∆t approximates the propagation of a density function ρ by advancing the projection coordinates ρ forward for a finite time, ∆t.Similarly, the matrix approximation L of the generator L approximates the infinitesimal action of the operator P t by approximating the time derivative of the projection coordinates C. Extension to controlled systems In the field of control theory, much attention has been given in recent years to applications the Koopman operator to control systems 50,51,69 , including several recent works which have noted the usefulness of formulating the problem in terms of the Koopman generator, rather than the Koopman operator [70][71][72][73][74] .Such a formulation in terms of the Koopman generator typically results in a lifted system that is bilinear in the control and lifted state, as the effect of the control vector fields is expressed in a way that is also dependent on the lifted state.This approach allows for a better approximation of the effects of control as compared to other common approaches 73 , especially for systems in control-affine form where the u i are control inputs and n c is the number of control inputs affecting the system.Here we apply a similar approach to the density transport problem, expressed in terms of the Perron-Frobenius generator.As shown by Peitz et al. 72 for the Koopman generator, by the property of the Perron-Frobenius generator given in Lemma 1, if the dynamics are control-affine, then the generators are also control affine, as can be seen by application of Eq. 4. This leads to density transport dynamics of the following form where L 0 is the Perron Frobenius generator associated with the vector field f (x) and similarly, the B i are the Perron Frobenius generators associated with the control vector fields g i (x).Therefore, given the finite dimensional approximation of these generators, we can approximate the density transport dynamics as where the matrices L 0 and B i are the matrix approximations of the operators in Eq. 18.These matrix approximations can be computed using the method outlined in Sec.II B for uncontrolled systems.This is done by first computing L 0 by Eq. 15 using trajectory data from the system with all control inputs set to zero.Once L 0 is found, each of the B i can be computed similarly by first computing a matrix L i by Eq. 15 using trajectory data from the system collected with u i = 1 and u j = 0 for j ̸ = i.This matrix L i approximates the Liouville operator corresponding to a vector field f + g i .The matrix approximation B i of the operator corresponding to the vector field g i alone is then found using Lemma 1 as For the systems of microrotors considered in this work, the control inputs are taken to be the strengths γ i of a pair of micro-rotors and the states are taken to be the position coordinates of a fluid particle.In this application, it will be shown (see, e.g.Eq. 34), that the control system is drift-free.That is, the vector field f = 0 in Eq. 17, and therefore, the corresponding Liouville operator L 0 = 0, as well for these systems.This is due to the typical quasistationary assumption of Stokes flows, which indicates that any change in the flow field is established instantaneously, without transience 75 .

Propagation of moments
In what follows, the problem of driving an initial density to a desired final density will be posed as an optimal control problem.The control inputs for this problem are the strengths of a finite number of micro-rotors, meaning that this problem involves steering a function using only a finite number of control inputs.To make this problem more tractable, we instead consider the problem of steering the moments of the density function to match the moments of a desired final density function.In the remainder of this section, an approximation of the moments of a density function ρ(x) are derived in terms of the projection of ρ onto the space spanned by the elements of D.
Given a projection of ρ onto D, as in Eq. 10, the first moment (mean), m 1 is written as where we use the superscript i in the moment to indicate the coordinate index and the subscript indicates the order of the moment being considered.Therefore, the first moment of ρ can be approximated as a linear combination of the means of the dictionary functions in Ψ, weighted by the projection coefficients ρ.This is also true for higher order raw moments, whereas higher order central moments become polynomial in ρ due to their dependence on the mean.Since ρ will be treated as the 'lifted state' in the control formulation, it is desirable to consider moments which are linear in ρ, so for this reason we will work with raw moments in what follows.
Here, for the dictionary functions, we use Gaussian radial basis functions of the form where c l is the center of the l th basis function, and s is a scaling parameter affecting the spread.Computing the integral in Eq. 20, in terms of this dictionary, the mean is approximated as where c i l is the i th coordinate of the l th basis function center.Similarly, the second raw moment can be written as where the last integral reduces to for a given basis function ψ l (x) where superscripts i and j are coordinate indices.

D. Finite-time coherent set detection
For autonomous dynamical systems, methods based on the Perron-Frobenius operator have been used to compute invariant or almost invariant sets of the system 76 .This is typically done by studying eigenfunctions of the Perron-Frobenius operator with eigenvalues, λ ≈ 1.Such eigenfunctions correspond to invariant or almost invariant densities, which describe groups of states which are left nearly unchanged by the flow of the system.These methods have also been extended to time-varying systems, in which the goal is to identify finite-time coherent sets [77][78][79][80] .Such sets are defined as sets in the state-space which are maximally coherent, or minimally dispersive, over a certain finite time interval.That is, they describe sets of states which may be transported as a whole by the flow, but with minimal transport outside of the coherent set or between coherent sets.These methods are also closely related to the Perron-Frobenius operator and are commonly seen as a probabilistic alternative to geometric methods related to the identification of invariant manifolds, dominant material lines, or Lagrangian coherent structures (see Refs. 79 and 81 for a review).
Here, we will apply the methods of Ref. 78 to the time-varying flow field generated by the solution to the optimal control problem to elucidate the flow structures associated with the optimal control.In this section, we will briefly summarize the method for the detection of coherent structures used here and its relation to the finitedimensional operator approximation defined in the previous section.
We assume that a dataset is given of m points, {(x i , y i )} m i=1 , where x i is the position of the i th particle at the initial time, t 0 and y i is the position of the particle at a later time t f .That is, t0 is the flow map associated with the non-autonomous system from time t 0 to t f .Given that the data lies in a set X at time t 0 and a set Y at time t f , our goal is to partition this dataset into two sets, X 1 and X 2 at time t 0 and Y 1 and Y 2 at time t f , such that points in X 1 are mapped into Y 1 by the flow and points in X 2 are mapped into Y 2 .This partition is designed by constructing partition functions f X and f Y which partition the space based on their sign.For example, we can define Then the problem of identifying coherent sets can be framed as choosing the functions f X and f Y to maximize the objective which can be thought of as an approximation of the an inner product where K .Note that this objective is only reasonable if an overall scale is imposed on the magnitude of the functions f X and f Y .If we approximate the partition functions f X and f Y by their projection onto the space spanned by the dictionary D, then the objective is approximated as where If we impose a scale by requiring that a T a = ãT ã = 1, then this maximization can be solved by singular value decomposition, with the optimal a and ã given by left and right singular vectors, respectively, as shown in Refs.77 and 78.This problem can be solved trivially by choosing f X to be uniform over X and choosing f Y to be uniform over Y -this solution typically corresponds to the singular vector associated with the largest singular value.Therefore, the singular vectors associated with the 2nd largest singular value give the optimal non-trivial solution, which divides the domain into partitions of roughly equal size 78 .

III. CONTROL FORMULATION
In Sec.II C, it was shown that the problem of steering a density ρ to a desired final density can be expressed as an output tracking problem on a lifted, bilinear system given by Eq. 19, where the projection coefficients ρ can be interpreted as a lifted state.Then, if the first and second raw moments are taken to be the relevant output, this can be expressed linearly in the lifted state, y = C ρ, where the elements of the output matrix C are given by rewriting Eqs.22, 23 in matrix form.
For the optimal output tracking problem, we consider a discrete time optimal control problem min u1,u2,...,u H−1 where H is the number of timesteps in the time horizon and Eq.28b represents the discrete time version of Eq. 19.
In particular, for output tracking, we consider inhorizon and terminal cost functions l and l H of the following quadratic forms where S, R, and S H are weighting matrices which define the penalty weight on tracking error, control effort, and error in the terminal state, respectively.Since the output y is linear in the lifted state ρt , this cost can be rewritten as a quadratic cost in terms of ρt , with an added linear term.
It is well known that for optimal control problems on bilinear systems with quadratic cost, an effective way of solving the problem is by iteratively linearizing and solving a finite time linear quadratic regulator (LQR) problem about a nominal trajectory, utilizing the Ricatti formulation of that problem 82 .For this reason, we solve the optimal control problem using differential dynamic programming (DDP) 83,84 , which is closely related to the method of iterative LQR.We briefly recount the primary steps of this algorithm below.
DDP computes a locally optimal control around a nominal trajectory by minimizing a quadratic approximation of the value function along this trajectory, and then doing this iteratively about the new trajectories obtained by applying the locally optimal control.First define the value function V (ρ t , t) at time t as, which expresses the optimal cost-to-go from ρt , where Denote by Q(δ ρ, δu) the change in the value function due to applying change in control input δu about the nominal trajectory and consider its quadratic approximation where these derivatives are given by where the notation (•) ′ indicates the next time step.The algorithm proceeds by computing these derivatives by recursing backward in time along the nominal trajectory from the end of the horizon.At each iteration, the control policy is improved by optimizing this quadratic expansion with respect to δu This can be seen as providing a descent direction in the space of control policies.An updated nominal control is then computed by a line search over a stepsize parameter α to update the policy, that is and this new control is applied to obtain a new nominal trajectory, and this procedure is iterated until the relative change in cost falls to less than a specified tolerance.For full details of the algorithm, the reader should refer to Refs. 83,84.

IV. TRANSPORT BY ROTORS IN FREE SPACE
To describe the fluid flow produced by a microscale rotor, we employ a model of a point torque in a two dimensional Stokes flow.Mathematically, this flow is described by a rotlet 33 , whose stream function is given by where x = (x, y) is the position a point in the fluid, x r = (x r , y r ) is the position of the rotlet, and γ is the strength of the rotlet.Physically, γ describes the magnitude of the point torque or the angular velocity of the rotor.The linearity of Stokes flows allows for the velocity fields produced by multiple rotlets to be determined by superposition of the velocity field produced by each rotlet individually.Therefore, for n r rotors, the resulting fluid flow results in the following fluid velocity field: where x i is the location of the i-th rotlet and Clearly, this results in a flow with a singularity at x r , circular streamlines around the singularity with counterclockwise flow for positive γ, and a fluid velocity that decays as r −2 going away from the rotor.
Here, we consider the case of rotors fixed in place on the x-axis at (−1, 0) and (1, 0), respectively, denoting their strengths by γ L and γ R for left and right.We consider a problem of manipulating a collection of fluid particles initially distributed at time t = 0 according to a normal distribution, with mean m 1 (0) = (1, 1) and covariance Σ(0) = 0.025I 2 , where I 2 is the 2 × 2 identity matrix.That is, ρ(x(0)) = N (1, 1), 0.025 I 2 .From this initial fluid particle distribution, we seek a sequence of rotor strengths over a timespan of 5 time units to drive the fluid particles to a final distribution with a mean of m 1 (5) = (−1, −1), while minimizing the variance.For this, we use the relationship between the second raw moment and the variance where the σ ij are the elements of the covariance matrix Σ, to convert the desired final variance to a desired second moment.With this, a rotor control sequence is found by solving the optimization problem as in Eq. 28a using the DDP scheme described in Sec.III.In solving this, the cost function weights are chosen to be S = 0.1I 5 , R = I 2 , and S H = 10 3 I 5 .That is, the error in the moments is penalized very little from t = 0 until t = 5, with a large penalty placed on the moments at t = 5.This choice allows the optimizer the flexibility to steer the distribution in a way that may temporarily increase the error if it results in a lower error in the moments at t = 5.
For the computation of the Liouville operators for this case, data is collected by simulating a grid of 2500 initial conditions, evenly spaced over [−2, 2] 2 forward for time interval ∆t = 0.005.The Perron-Frobenius operators are computed using Eq.14 with this trajectory data and a 25 × 25 grid of Gaussian radial basis functions with centers evenly spaced over the same domain, excluding small radii around the rotors.The Liouville operators are obtained from this using Eq. 15 as described in Sec.II C. Fig. 1 shows the effect of the rotor control on the motion of a distribution of 10 4 fluid particles sampled according to the initial density and displayed as a histogram approximation of the density.The rotor positions are indicated by the circle-cross and the position of the target mean is shown by the green circle.The white-filled circle indicates the mean of this sample, while the black-filled circle indicates the mean as predicted using the Liouville operator.The streamlines in the figure indicate direction of the fluid velocity field produced by the rotors at the indicated time instant.Fig. 2 shows the rotor strengths γ L and γ R selected by the DDP algorithm.For the first 1.25 seconds, the rightmost rotor has a positive strength of near 1 to generate a counterclockwise flow, pulling the distribution of particles toward the origin, while the leftmost rotor has a low strength near zero.As the distribution nears the origin, the strength of the right rotor decreases, while the magnitude of the strength of the left rotor increases to generate a clockwise flow, which pulls the distribution toward the target mean.
Also shown in Fig. 2 are plots of the elements of the first and second moment over time as computed from the sample shown in Fig. 1 (labelled 'True') and as predicted using the finite approximation of the Liouville operators (labelled 'Predicted').

A. Finite time coherent sets
With the optimal control determined, Eq. 34 gives a nonautonomous dynamical system.We can then apply the methods outlined in Sec.II D to this system to identify coherent sets to better understand the underlying structure of the flow field produced by the optimal control.For this computation, we use a dataset of 10 4 data pairs, initially spaced on a uniform grid over [−2, 2] 2 .For the basis, we use a set of 2501 basis functions consisting of Gaussian radial basis functions uniformly spaced on a 50 × 50 grid over [−2, 2] 2 and the constant function, ψ = 1.Fig. 3 shows the time evolution of the data set, with the points colored according to the partition function f X as approximated by the 2nd left singular vector of A = 1 m Ψ X Ψ T Y .Also shown is the evolution of the f X = 0 contour, which approximates the barrier between the coherent sets, as depicted by the black line.Finally, level sets and mean of the density function, ρ(x(t)), as approximated from a sample of 10 4 points from the initial density, are shown by the purple contours and purple markers, respectively.The level sets shown correspond to values of the initial density at one and two standard deviations from the mean, respectively.
To quantify the coherence of the sets identified by the partition functions f X and f Y , we use a modification of the objective in Eq. 24, which only considers the sign of f X and f Y , which effectively gives the fraction of the data points which are classified correctly by the partition functions (for which the partition functions do not change sign from initial to final time).For the case shown in Fig. 3, we have ḡ = 0.9904.This computation of the coherent sets shows that the flow field generated by the optimal control is such that a transport barrier is formed over the 5s time interval, with the barrier passing through the particle distribution at the initial time and connecting to the target location at the final time.Many previous works 35,[85][86][87] have studied the relationship between the optimal control problem of steering a particle efficiently in an unsteady flow and the coherent structures associated with that flow.Typically in these studies, the problem being considered is motivated by the efficient navigation of an underwater vehicle to a target in an unsteady ocean flow.For this reason, the control input is usually taken to be a propulsive velocity which is added to the unsteady flow field, as could be generated by a thruster onboard an underwater vehi- cle, and coherent structures associated with the unsteady flow are used to identify efficient routes.Our work takes a different perspective, where instead of controlling individual particles in a given unsteady flow field, we solve an optimal control problem to determine the optimal timevarying flow field to steer the initial particle distribution to the target, where the unsteady flow field is constrained to be a superposition of flow fields produced by the two rotors at each time instant.In previous works 35,85 , it was seen that the optimal routes of an underwater vehicle tend to follow the coherent structures which guide the particle towards the target for energy-optimal navigation.Here we see that the optimal flow field produces a flow structure which guides the distribution of particles from the initial condition to the target, as shown in Fig.

3.
This sort of flow structure seems to be typical of the optimal control solutions in this setting.To verify this, we solve the control problem with the same parameters but with an initial density centered about (−0.5, 1).That is, the initial density is ρ(x(0)) = N ((−0.5, 1), 0.025I 2 ).Fig. 4 shows the 10 4 data points colored according to the 3rd left singular vector of A. With both the initial distribution and the target in the left half of the domain, the second left singular vector simply divides the domain roughly into its left and right halves.However, for this case the third singular vector shows a partition which indicates a coherent structure that extends from the initial blob location at the initial time (see Fig. 4 (a)) to the target at the final time (see Fig. 4 (b)).Evaluating the objective in Eq. 36 for this case, we have that Y refer to the partition functions given by the second and third singular vectors, respectively.

V. TRANSPORT BY ROTORS NEAR AN INFINITE PLANE WALL
For the case of a rotlet located at a point x r above and infinite plane wall at y = w, the fluid flow must satisfy the additional boundary conditions of no slip and no penetration at the plane wall.The stream function associated with this flow is given by 33,88 where x im = (x im , y im ) = (x r , 2w − y r ) is the location of an image singularity which has the effect of making the fluid velocity vanish at the plane wall.Similarly, r im = |x − x im |.Therefore, the flow produced for this case is u w = (u w , v w ) where and the velocity field for multiple rotlets above a plane wall can be found by summing the individual velocity fields as in Eq. 34.
With these governing equations for the flow produced by microrotors in the presence of a plane wall, we consider a similar transport problem to the one considered in Sec.IV in order to examine the boundary effects of the plane wall on the transport problem.As in Sec.IV, the same rotor positions of (−1, 0) and (1, 0), initial density of ρ(x(0)) = N (1, 1), 0.025 I 2 , timespan of 5 units, target moments, and cost function are considered.The Liouville operators are computed using trajectory data from the same grid of initial conditions and basis functions positioned on the same grid as in Sec.IV, but with any points in these grids lying outside of the fluid domain (below the plane wall) neglected.
Fig. 5 shows the resulting flow field and its effect on the motion of the particle distribution for the case of a plane wall located at w = −1.25 from the control computed using the DDP algorithm.From this figure, it is clear that the effect of the wall is to stretch the particle distribution along the wall due to the vanishing fluid velocity at the wall.Due to this effect, the control tends to pull the distribution to the left in the early stages of the trajectory using a larger positive (counterclockwise) strength of the left rotor than in the free space case.Related to this, in the middle stages of the trajectory, a larger positive strength of the right rotor is needed to supplement effects of the left rotor, as compared to the free space case.These effects can also be clearly seen in Fig. 6 (b), which shows a time sequence of the rotor strengths for this problem for varying wall locations.Fig. 6 (a) shows an overlay of the final particle distribution at t = 5 for the same wall locations.From this figure, it can be seen that effect of the wall is to elongate the distribution more for cases where the wall is closer to the target mean location.Fig. 6 (c) shows a comparison of the optimal cost found from the DDP algorithm at varying wall locations, which demonstrates that the cost increases significantly as the wall nears the target mean position.This is due to both to the increased control effort (rotor strength) needed to steer the distribution as well as well as greater error in the moments due to the stretching effect of the wall.Fig. 7 shows the coherent sets for the case shown in Fig. 5 at the initial and final times.As in the free space case, the optimal control forms a coherent structure which passes near to the initial blob location at the initial time and extends toward the target at the final time.

VI. TRANSPORT BY ROTORS WITHIN A CIRCULAR BOUNDARY
For the case of a rotlet positioned at a point x r inside of a circular boundary of radius, a, centered about the origin, again the no-slip and no penetration boundary conditions must be satisfied by the flow at the boundary, and again, these can be satisfied by modifying the stream function to include image terms to cancel out the flow at the wall.The stream function satisfying these conditions can be shown to be 33,37,88 where R = |x| and R r = |x r | are the radial distances from the center of the circle to the evalutation point and to the rotlet, respectively, and x r is the location of the image system.That is, the image is located outside of the circular boundary at a point along the line between the center of the circle and the rotlet at a radial distance of a 2 /R r from the center of the circle.Then the flow field  for this case is given by u c = (u c , v c ) where With these governing equations for the flow produced by microrotors within a circular boundary, we consider the same transport problem considered in previous cases in order to examine the boundary effects of the circu-lar boundary on the transport problem.As before, the rotor positions of (−1, 0) and (1, 0), initial density of ρ(x(0)) = N (1, 1), 0.025 I 2 , timespan of 5 units, the same target moments and cost function are considered.The Liouville operators are computed using trajectory data from the same grid of initial conditions and basis functions positioned on the same grid as in Sec.IV, but with any points in these grids lying outside of the fluid domain (beyond the circular boundary) neglected.
Fig. 8 shows the resulting flow field from the control and its effect on the motion of the particle distribution for the case of the two rotors within a circular boundary of radius a = 2.5.Similarly to the case next to a plane wall, the reduced fluid velocity near the circular boundary leads to a stretching effect on the distribution, especially when a significant part of the particle distribution lies in regions near to the boundary.Since this is encountered at the initial condition, significantly more particles remain in the upper, trailing 'tail' of the distribution due to the drag effects of the boundary in the upper right quadrant.This effect becomes more apparent for smaller boundary radius.Due to this effect, more control effort must be exerted by the rotors in the early stages of the trajectory to overcome this drag.A secondary effect of this is that the leading tail of the distribution, which consists of particles closer to the interior of the circle and further from the boundary, tends to stretch more, leading it to wrap around the rightmost rotor in the later stages of the trajectory in a way that was not seen in the previous cases.These qualitative differences are highlighted in Fig. 9 (a), which shows an overlay of the final particle distribution at t = 5 for the same boundary radius.Fig. 9 (b) shows a time sequence of the rotor strengths for this problem for varying wall locations.Fig. 10 shows the coherent sets for the case shown in Fig. 8 at the initial and final times.As in the previous cases, the optimal control produces a flow field a coherent structure which passes near to the initial blob location at the initial time and extends toward the target at the final time.

VII. TRANSPORT OF TWO DENSITIES
We now return to the case of two micro-rotors in free space to consider the problem of manipulating two distinct distributions of fluid particles to a common target mean and second moment.This requires a reformulation of the optimal control problem as posed in Eq. 28a.In that formulation, the state of the control problem was taken to be the vector of projection coefficients ρ.Here we consider an augmented state containing the projection coefficients of the two density functions.Denoting these two density functions as ρ A and ρ B , and their corresponding projection coefficients by ρA and ρB , the augmented state for this case is [(ρ A ) T , (ρ B ) T ] T .Similarly, we consider an output vector which concatenates the first and second moments for the two density functions y = [(m A ) T , (m B ) T ] T , where m A and m B are vectors containing the moments of the densities ρ A and ρ B respectively, as in Eq. 27.The same Liouville operators are used to propagate each of these densities forward in time.From this point, an appropriate cost function can be specified and the rotor control can be optimized using the DDP scheme as before.
With this formulation, we consider the problem of manipulating two densities using two rotors fixed at the same locations as before, (−1, 0) and (1, 0).We take the initial density for one of the distributions to be the same as the previous examples, ρ A (x(0)) = N (1, 1), 0.025 I 2 , and consider a second distribution starting from an initial density of ρ B (x(0)) = N (b, 1), 0.025 I 2 , where b is a parameter to be varied.This formulation allows us to examine the ability to steer two distributions starting from varying initial distances apart.We consider the problem of choosing the rotor strengths to steer both of these distributions to a final distribution with a mean of m A 1 (5) = m B 1 (5) = (−1, −1), while minimizing the variances.For this problem, the cost function is taken to be of the same form as Eq. 29 with the weights chosen to be S = 0.1I 10 , R = I 2 , and S H = 500I 10 .That is, the terminal cost is chosen to be half that of the previous cases since it is being applied to the error in the moments of two density functions and summed.distribution of ρ B starts farther from the initial distribution of ρ A , a higher negative spin is applied by the left rotor in the early stages of the trajectory, producing a flow that is more symmetric as the blobs are pulled toward the middle, but with a similar flow near the end of the trajectory as the distributions near the target.In the last case shown, where b = −1.0, it appears that a transition has occurred and a qualitatively different optimal trajectory is found in which the leftmost distribution ρ B is stirred counterclockwise around the left rotor rather than through the region between the rotors.This is done by a positive torque applied from the left rotor, which also results in the rightmost blob ρ A being pulled to a position above the left rotor.As a result of this, at the end of the sequence, the rightmost rotor generates a counterclockwise flow which pushes the two distributions down toward the target.This is in contrast to the other cases, where the right rotor generates a clockwise flow near the end in order to steer the particles from right to left toward the target.These effects can also be seen by examining the rotor strengths directly, as shown in Fig. 11.
Fig. 13 shows the coherent sets at the initial and final time for the cases shown in Fig. 12 where two distributions are to be steered to the common target.In the first three cases considered, the coherent structure which divides the coherent sets at the initial time passes through the regions of high concentration of both initial distributions.At the final time, this structure moves toward the target, effectively pulling both distributions toward the goal.

VIII. CONCLUSION
A promising new approach has been developed and demonstrated for computing the optimal control to transport a distribution of states whose dynamics are governed by a control affine system to a desired final state distribution in a fixed, finite time.We demonstrate the usefulness of this method by highlighting a fluid mechanical application, in which the relevant state is the position of a fluid particle, the distribution describes a blob of fluid particles, and the controls are the torques applied by a pair of fixed rotors, which stir the flow in circular patterns.In this setting, we used the proposed approach to analyze the effects of fixed boundaries on the transport problem.We believe that such control strategies will be very useful in applications, such as targeted drug delivery, particle manipulation, and cell sorting in which the relevant transport problem is not to mix the fluid, but to transport a concentrated distribution of particles in a controlled way to a desired location.In future works, we plan to study similar transport problems in which the flow is generated by non-stationary stirrers, such as a moving rotors or microswimming robots 89 , or by boundary controls.Other interesting use case of the work presented here could be using this algorithm to optimize rotor placement for a given task.This application could be especially relevant for the design of microfluidic devices where fluid transport is critical.
While it was demonstrated on and motivated by problems in the fluids setting, we believe that the proposed approach can have much broader application in control systems, where the density of states can be taken to represent an uncertainty distribution 90 .Other exciting extensions of this work could include understanding the relationship between this method and the formation, motion, and manipulation of transport barriers in a flow field.
are the Koopman and Perron-Frobenius operators associated with this time-varying flow and Φ t0 t f = Φ t f t0 −1

FIG. 1 .
FIG. 1. Controlled transport of a distribution of fluid particles by two micro-rotors fixed at (−1, 0) and (1, 0) from an initial density ρ(x(0)) = N ((1, 1), 0.025I2) to a target mean at (−1, −1) (green circle).White filled circle indicates the sample mean and black filled circle indicates the mean predicted by the proposed method.Streamlines depict the flow field produced by the rotor control at the instant shown.

2 FIG. 2 .
FIG.2.Rotor controls and fluid particle distribution moments for transport by rotors in free space as shown in Fig.1.Left, Top: Rotor strengths for the left rotor (γL) and right rotor (γR).Left, Bottom: Mean of particle distribution from sample shown in Fig.1 ('True') and as predicted using the Liouville operator ('Predicted').Right: Second raw moment.Markers correspond to instants shown in Fig.1.

FIG. 3 .
FIG. 3. Finite-time coherent sets produced by the optimal rotor-driven flow field corresponding to the case shown in Figs. 1 and 2. The sequence shows the time evolution of 10 4 particles initially placed on a uniform grid on [−2, 2] 2 and colored according to the partition function fX , as approximated by the left singular vector of A = 1 m ΨX Ψ T Y .The black line in the first image depicts the fX = 0 level set, and the sequence shows the evolution of this line as the dataset is advected by the flow.The purple marker and lines show the mean value and level sets of the density function corresponding to values of the initial density at one and two standard deviations from the mean, respectively.

FIG. 4 .
FIG. 4. Finite-time coherent sets produced by the optimal rotor-driven flow field with initial density ρ(x(0)) = N ((−0.5, 1), 0.025I2) .The subfigures show the 10 4 data points at the initial time t = 0 and final time t = 5, colored according to the 3rd left singular vector of A, along with the fX = 0 contour (black line).Purple contours depict the mean value and level sets of the density function corresponding to values of the initial density at one and two standard deviations from the mean, respectively.

FIG. 5 .
FIG. 5. Controlled transport of a distribution of fluid particles near an infinite plane wall at y = −1.25 (black line) by two micro-rotors fixed at (−1, 0) and (1, 0) from an initial density ρ(x(0)) = N ((1, 1), 0.025I2) to a target mean at (−1, −1) (green circle).White filled circle indicates the sample mean and black filled circle indicates the mean predicted by the proposed method.Streamlines depict the flow field produced by the rotor control at the instant shown.

FIG. 6 .FIG. 7 .
FIG. 6.Comparison of boundary effects due to an infinite plane wall at various locations, y = w.(a) Overlaid final particle distributions at t = 5 for varying wall locations.(b) Left and right rotor strengths as selected through the DDP algorithm for varying wall locations.(c) Total cost from the optimal control problem near a wall, Jw normalized by the cost for the free space case, J0.

FIG. 8 .
FIG. 8. Controlled transport of a distribution of fluid particles within a circular boundary of radius a = 2.5 (black circle) by two micro-rotors fixed at (−1, 0) and (1, 0) from an initial density ρ(x(0)) = N ((1, 1), 0.025I2) to a target mean at (−1, −1) (green circle).White filled circle indicates the sample mean and black filled circle indicates the mean predicted by the proposed method.Streamlines depict the flow field produced by the rotor control at the instant shown.

FIG. 9 .
FIG. 9. Comparison of boundary effects due to a circular boundary of varying radius, a.(a) Overlaid final particle distributions at t = 5 for varying boundary radius.(b) Left and right rotor strengths as selected through the DDP algorithm for varying boundary radius.(c) Total cost from the optimal control problem in a circular boundary, Jc normalized by the cost for the free space case, J0.

FIG. 10 .
FIG. 10.Finite-time coherent sets for the case shown in Fig. 8 with a circular boundary of radius a = 2.5.Left: the initially uniform data set colored according to the partition function fX , as approximated by the left singular vector of A = 1 m ΨX Ψ T Y .Right: The same points at the final time.The black line shows the fX = 0 level set.The purple marker and lines show the mean value and level sets of the density at values of the initial density at one and two standard deviations from the mean.

FIG. 13 .
FIG. 13.Finite-time coherent sets for the the transport of two distributions, as shown in Fig. 12. Left column: the initially uniform data set colored according to the partition function fX , as approximated by the left singular vector of A = 1 m ΨX Ψ T Y .Right: The same points at the final time.The black line shows the fX = 0 level set.The purple marker and lines show the mean value and level sets of the density at values of the initial density at one and two standard deviations from the mean.Each row shows the results for a different value of b: first row: b = 0.5, second row: b = 0, third row: b = −0.5, fourth row: b = −1.