# Controlled density transport by microrotors in a Stokes flow using linear
transfer operators

Jake Buzhardt  Phanindra Tallapragada<EMAIL_ADDRESS>,
<EMAIL_ADDRESS>Department of Mechanical Engineering, Clemson University

(20 October 2023)

###### Abstract

We study the problem of transporting a distribution of fluid particles in a
Stokes flow to a desired final distribution in a fixed, finite time by
controlling the torques of a pair of microrotors at fixed positions in the
flow. Our approach is based on a finite dimensional approximation of the
Liouville operator, the infinitesimal generator of the semi-group of Perron
Frobenius operators, which describes the density transport dynamics associated
with the microrotor flow fields. Using this operator, we express the transport
problem as an optimal control problem in terms of the moments of the density
function of the particle distribution. The finite time optimal control problem
is then solved using differential dynamic programming, an iterative trajectory
optimization method. We apply this framework to the microrotor driven flow on
four related problems: transport using rotors in an unbounded flow, transport
near to an infinite plane wall, transport within a circular domain, and the
simultaneous transport of two particle distributions to a common final
distribution in an unbounded flow. These examples demonstrate the
effectiveness of the proposed framework and also allow us to better understand
the effects of boundaries on the ability to achieve a desired fluid transport
using a rotor-driven flow.

††preprint: AIP/123-QED

> The problem of transporting a blob of fluid particles in a Stokes flow to a
> desired distribution in a fixed, finite time period has several important
> applications. This 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.
> Manipulating the flow allows control of dynamical structures such as almost
> invariant 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 this operator,
> 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 free space, near to a plane wall,
> in a circular confinement, and the transport of two blobs to a common
> target. These examples demonstrate the effectiveness of the proposed method
> and help us to understand the abilities and limitations of fluid transport
> in a rotor-driven flow.

##  I Introduction

Understanding and controlling the motion of fluid particles in the low
Reynolds number regime has become increasingly significant in recent years,
particularly in the realm of microrobotics and microfluidics. Microrotors and
micropumps propelled by various mechanisms have been proposed as a useful
means of transporting fluid particles or other submersed cargo in a
microfluidic solution Iverson and Garimella (2008); Zhou  _et al._ (2016); Xu
_et al._ (2017); Wang and Fu (2018). In this paper, we develop a method based
on recent advances in data-driven dynamical systems to model and control the
transport of distributions of fluid particles using microrotors in a Stokes
flow.

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
Aref (1984); Aref  _et al._ (2017); Ottino (1989). 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 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 Chwang and Wu (1975); Pozrikidis
(1992). This work is an extension of the authors’ recently submitted
paperBuzhardt and Tallapragada (2023). 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.

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 Meleshko and Aref (1996) 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 Aref (1984). 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 dimensions. Van der Woude et al. Van Der Woude
_et al._ (2007) 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 Lushi and Vlahovska (2015); Tallapragada and Sudarsanam
(2019); Delmotte (2019). 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 Mathew  _et al._ (2007); D’Alessandro, Dahleh, and
Mezic (1999); Cortelezzi, Adrover, and Giona (2008); Lin, Thiffeault, and
Doering (2011); Hassanzadeh, Chini, and Doering (2014); Zhang and Balasuriya
(2020), with most of these focusing on optimizing mixing performance. Mathew
et al Mathew  _et al._ (2007) 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
BalasuriyaZhang and Balasuriya (2020) develop a method to determine an optimal
spatiotemporally varying additive 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 Lasota and Mackey (1994), 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 Otto and Rowley (2021);
Kaiser, Kutz, and Brunton (2020), with many of the most common methods having
their origin in the analysis of fluid flows Rowley  _et al._ (2009); Mezić
(2013). Refs. Froyland, González-Tokman, and Watson, 2016; Froyland and
Santitissadeekorn, 2017 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 Klünker, Padberg-Gehle, and
Thiffeault (2022) recently studied mixing in open flows in terms of spectral
properties of a finite-rank approximation of the Perron Frobenius operator.
Sinha et alSinha, Vaidya, and Rajaram (2016) 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. BrockettBrockett (2007,
2012) 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 Grover and
Elamvazhuthi (2018); Elamvazhuthi and Grover (2018) use transfer operators and
their generators in a graph-based 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 Grover and
Elamvazhuthi (2018); Elamvazhuthi and Grover (2018) 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.3 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.4 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.

###  II.1 Perron Frobenius operator and generator

Consider a dynamical system

$\frac{dx}{dt}=f(x)$ (1)

on a measure space $(\mathbf{X},\mathcal{A},\mu)$ where $x\in\mathbf{X}$ is
the state, $\mathbf{X}\subset\mathbb{R}^{n}$ is the state space, $\mathcal{A}$
is the Borel $\sigma$-algebra on $\mathbf{X}$, and $\mu$ is a measure on
$\mathbf{X}$. Denote the time-$t$ flow map from an initial state $x_{0}$ by
$\Phi^{t}(x_{0})$. We will further assume that the measure $\mu$ is absolutely
continuous with respect to the Lebesgue measure, so that $\mu$ can be
expressed in terms of a density, $\rho\in L_{1}(\mathbf{X})$, such that
$d\mu(x)=\mu(dx)=\rho(x)dx$. With this, the Perron Frobenius operator,
$\mathcal{P}^{t}:L_{1}(\mathbf{X})\mapsto L_{1}(\mathbf{X})$ corresponding to
the flow $\Phi^{t}$ can be defined as the unique operator Lasota and Mackey
(1994) such that

$\int_{A}\mathcal{P}^{t}\rho(x)dx=\int_{(\Phi^{t})^{-1}(A)}\rho(x)dx$ (2)

for any $A\in\mathcal{A}$, $t\geq 0$. The family of these operators,
parameterized by time, $t$, have been shown to satisfy the properties of a
semigroup Lasota and Mackey (1994). The infinitesimal generator of this
semigroup, denoted here by $\mathcal{L}$, is known as the Liouville operator
or the Perron-Frobenius generator, and defined as

$\mathcal{L}\rho=\lim_{t\to 0}\frac{\mathcal{P}^{t}\rho-\rho}{t}=\lim_{t\to
0}\left(\frac{\mathcal{P}^{t}-\mathcal{I}}{t}\right)\rho$ (3)

where $\mathcal{I}$ is the identity operator. Alternatively, as this operator
expresses the deformation of a density function under an infinitesimal action
of the operator $\mathcal{P}^{t}$, the Liouville operator can be thought of as
expressing a continuity equation for the number of particles in the state
space Lasota and Mackey (1994); Cvitanovic  _et al._ (2005); that is,

$\frac{\partial\rho}{\partial t}=\mathcal{L}\rho=-\nabla_{x}\cdot(\rho f)~{}.$
(4)

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}:\mathbf{X}\mapsto\mathbb{R}^{n}$ is denoted by $\mathcal{L}_{1}$ and
the Liouville operator associated with the vector field
$f_{2}:\mathbf{X}\mapsto\mathbb{R}^{n}$ by $\mathcal{L}_{2}$, then the
Liouville operator associated with the vector field $f(x)=f_{1}(x)+f_{2}(x)$,
is $\mathcal{L}=\mathcal{L}_{1}+\mathcal{L}_{2}$.

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
$\mathcal{K}^{t}:L^{\infty}(\mathbf{X})\mapsto L^{\infty}(\mathbf{X})$ is the
operator which propagates observable functions $h\in L^{\infty}(\mathbf{X})$
forward in time along trajectories of the system and is defined as

$\mathcal{K}^{t}h=h\circ\Phi^{t}~{}.$ (5)

The Koopman and Perron-Frobenius operators are adjoint to one another, with
the adjoint relationship given by

$\int_{\mathbf{X}}[\mathcal{K}^{t}h](x)\rho(x)dx=\int_{\mathbf{X}}h(x)[\mathcal{P}^{t}\rho](x)dx\,.$
(6)

###  II.2 Numerical approximation

One of the most common methods of approximating the Perron-Frobenius operator
is a set-oriented approach known as Ulam’s methodUlam (1960), 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 cellsDellnitz,
Froyland, and Junge (2001). 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 Klus,
Koltai, and Schütte (2016). In recent works involving numerical approximation
of the Koopman operator, one of the most common approaches is that of extended
dynamic mode decomposition (EDMD) Williams, Kevrekidis, and Rowley (2015), 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 Williams, Kevrekidis, and
Rowley (2015); Klus, Koltai, and Schütte (2016). 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 Huang and Vaidya (2018);
Goswami, Thackray, and Paley (2018); Klus, Koltai, and Schütte (2016). In this
work, we also implement EDMD for the computation of the Perron-Frobenius
operator, which we outline below, largely following Klus _et al_. Klus,
Koltai, and Schütte (2016).

The method requires a predefined dictionary $\mathbb{D}$ of $k$ scalar-valued
basis functions, $\mathbb{D}=\\{\psi_{1},\psi_{2},\dots,\psi_{k}\\}$, where
$\psi_{i}:\mathbf{X}\mapsto\mathbb{R}$ for $i=1,\dots,k$ and trajectory data
collected from the dynamical system with fixed timestep, $\Delta t$, arranged
into snapshot matrices as

$\displaystyle X$
$\displaystyle=\begin{bmatrix}x_{1}&,~{}\cdots~{},&x_{m}\end{bmatrix}$ (7)
$\displaystyle Y$
$\displaystyle=\begin{bmatrix}x_{1}^{+}&,~{}\cdots~{},&x_{m}^{+}\end{bmatrix}$
(8)

where the subscript $i=1,\dots,m$ is a measurement index and
$x_{i}^{+}=\Phi^{\Delta t}(x_{i})$.

Then, given an the observable function $h$ and density $\rho$, these functions
are approximated by their projections onto the space spanned by elements of
$\mathbb{D}$ as

$\displaystyle h(x)$ $\displaystyle\approx\hat{h}^{T}\Psi(x)$ (9)
$\displaystyle\rho(x)$ $\displaystyle\approx\Psi^{T}(x)\hat{\rho}$ (10)

where $\hat{h},\,\hat{\rho}\in\mathbb{R}^{k}$ are column vectors containing
the projection coefficients and $\Psi:\mathbf{X}\mapsto\mathbb{R}^{k}$ is a
column-vector valued function where the elements are given by
$[\Psi(x)]_{i}=\psi_{i}(x)$. Substituting these expansions into Eq. (6) yields

$\int_{\mathbf{X}}\mathcal{K}^{\Delta
t}[\hat{h}^{T}\Psi]\Psi^{T}\hat{\rho}\,dx=\int_{\mathbf{X}}\hat{h}^{T}\Psi\mathcal{P}^{\Delta
t}[\Psi^{T}\hat{\rho}]\,dx\,.$ (11)

Then noting that $[\mathcal{K}^{\Delta t}\Psi](x)=\Psi(x^{+})$ and assuming
that $\mathcal{P}^{\Delta t}$ can be approximated by a matrix $P$ operating on
the coordinates $\hat{\rho}$, it is clear that in the limit of a large dataset
$m\to\infty$, the above expression becomes

$\Psi_{Y}\Psi_{X}^{T}=\Psi_{X}\Psi_{X}^{T}P+e$ (12)

where $e$ is a residual error arising due to the matrix approximation of
$\mathcal{P}^{\Delta t}$ by $P$. This can be posed as a least-squares problem
for the matrix $P$

$\min_{P}\|\Psi_{Y}\Psi_{X}^{T}-\Psi_{X}\Psi_{X}^{T}P\|_{2}^{2}$ (13)

where $\Psi_{X}$,$\Psi_{Y}\in\mathbb{R}^{k\times m}$ are matrices with columns
containing $\Psi$ evaluated on the columns of $X$ and $Y$ respectively. The
analytical solution of this least squares problem is

$P=\left(\Psi_{X}\Psi_{X}^{T}\right)^{\dagger}\Psi_{Y}\Psi_{X}^{T}$ (14)

where $(\cdot)^{\dagger}$ is the Moore-Penrose pseudoinverse.

Given this matrix approximation of the operator, $P$, if the timestep $\Delta
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

$L\approx\frac{P-I_{k}}{\Delta t}$ (15)

where $I_{k}$ is the $k\times k$ identity matrix. The matrix approximation $P$
of the operator $\mathcal{P}^{\Delta t}$ approximates the propagation of a
density function $\rho$ by advancing the projection coordinates $\hat{\rho}$
forward for a finite time, $\Delta t$. Similarly, the matrix approximation $L$
of the generator $\mathcal{L}$ approximates the infinitesimal action of the
operator $\mathcal{P}^{t}$ by approximating the time derivative of the
projection coordinates

$\frac{d\hat{\rho}}{dt}=L\hat{\rho}\,.$ (16)

###  II.3 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 Korda and Mezić
(2018); Kaiser, Kutz, and Brunton (2020); Otto and Rowley (2021), including
several recent works which have noted the usefulness of formulating the
problem in terms of the Koopman generator, rather than the Koopman operator
Goswami and Paley (2017); Klus  _et al._ (2020); Peitz, Otto, and Rowley
(2020); Bruder, Fu, and Vasudevan (2021); Folkestad and Burdick (2021). 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 approachesBruder, Fu, and Vasudevan
(2021), especially for systems in control-affine form

$\frac{dx}{dt}=f(x)+\sum_{i=1}^{n_{c}}g_{i}(x)u_{i}$ (17)

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. Peitz, Otto, and Rowley (2020) 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

$\frac{\partial\rho}{\partial
t}=\mathcal{L}_{0}\rho+\sum_{i=1}^{n_{c}}u_{i}\mathcal{B}_{i}\rho$ (18)

where $\mathcal{L}_{0}$ is the Perron Frobenius generator associated with the
vector field $f(x)$ and similarly, the $\mathcal{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

$\frac{d\hat{\rho}}{dt}=L_{0}\hat{\rho}+\sum_{i=1}^{n_{c}}u_{i}B_{i}\hat{\rho}$
(19)

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.2 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\neq 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 $B_{i}=L_{i}-L_{0}$.

For the systems of microrotors considered in this work, the control inputs are
taken to be the strengths $\gamma_{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 $\mathcal{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 Happel and Brenner (1983).

####  II.3.1 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 $\rho(x)$ are derived in
terms of the projection of $\rho$ onto the space spanned by the elements of
$\mathbb{D}$.

Given a projection of $\rho$ onto $\mathbb{D}$, as in Eq. 10, the first moment
(mean), $m_{1}$ is written as

$m_{1}^{i}=\int x^{i}\rho(x)dx=\hat{\rho}^{T}\int x^{i}\Psi(x)dx$ (20)

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 $\rho$ can be approximated as a linear
combination of the means of the dictionary functions in $\Psi$, weighted by
the projection coefficients $\hat{\rho}$. This is also true for higher order
raw moments, whereas higher order central moments become polynomial in
$\hat{\rho}$ due to their dependence on the mean. Since $\hat{\rho}$ will be
treated as the ‘lifted state’ in the control formulation, it is desirable to
consider moments which are linear in $\hat{\rho}$, 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

$\psi_{l}(x)=\exp\left(-\frac{(x-c_{l})^{T}(x-c_{l})}{2s^{2}}\right)$ (21)

where $c_{l}$ is the center of the $l^{\text{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

$m_{1}^{i}=2\pi s^{2}\sum_{l=1}^{k}\hat{\rho}_{l}c^{i}_{l}$ (22)

where $c_{l}^{i}$ is the $i^{\text{th}}$ coordinate of the $l^{\text{th}}$
basis function center. Similarly, the second raw moment can be written as

$m_{2}^{ij}=\int x^{i}x^{j}\rho(x)dx=\hat{\rho}^{T}\int x^{i}x^{j}\Psi(x)dx$
(23)

where the last integral reduces to

$\int x^{i}x^{j}\psi_{l}(x)dx=\begin{cases}2\pi
s^{2}(s^{2}+(c_{l}^{i})^{2})&i=j\\\\[4.30554pt] 2\pi
s^{2}c_{l}^{i}c_{l}^{j}&i\neq j\end{cases}$

for a given basis function $\psi_{l}(x)$ where superscripts $i$ and $j$ are
coordinate indices.

###  II.4 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 Dellnitz and Junge (1999). This is typically done by studying
eigenfunctions of the Perron-Frobenius operator with eigenvalues,
$\lambda\approx 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
setsFroyland, Santitissadeekorn, and Monahan (2010); Williams, Rypina, and
Rowley (2015); Allshouse and Peacock (2015). 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. Allshouse and Peacock, 2015; Hadjighasem  _et al._ ,
2017 for a review).

Here, we will apply the methods of Ref. Williams, Rypina, and Rowley, 2015 to
the time-varying flow field generated by the solution to the optimal control
problem to illucidate 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 finite-dimensional
operator approximation defined in the previous section.

We assume that dataset is given of $m$ points,
$\\{(x_{i},y_{i})\\}_{i=1}^{m}$, where $x_{i}$ is the position of the
$i^{\text{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,
$y_{i}=\Phi_{t_{0}}^{t_{f}}(x_{i})$, where $\Phi_{t_{0}}^{t_{f}}$ 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 $X_{1}=\\{x\in X|f_{X}(x)>0\\}$. Then the problem
of identifying coherent sets can be framed as choosing the functions $f_{X}$
and $f_{Y}$ to maximize the objective

$g(f_{X},f_{Y})=\frac{1}{m}\sum_{i=1}^{m}f_{X}(x_{i})f_{Y}(y_{i})$ (24)

which can be thought of as an approximation of the an inner product

$\displaystyle g(f_{X},f_{Y})$ $\displaystyle\approx\langle
f_{X},\mathcal{K}_{t_{0}}^{t_{f}}f_{Y}\rangle=\int_{X}f_{X}(x)f_{Y}\left(\Phi_{t_{0}}^{t_{f}}(x)\right)dx$
(25a)
$\displaystyle\approx\langle\mathcal{P}_{t_{0}}^{t_{f}}f_{X},f_{Y}\rangle=\int_{Y}f_{X}\left(\Phi_{t_{f}}^{t_{0}}(y)\right)f_{Y}(y)dy$
(25b)

where $\mathcal{K}_{t_{0}}^{t_{f}}$ and $\mathcal{P}_{t_{0}}^{t_{f}}$ are the
Koopman and Perron-Frobenius operators associated with this time-varying flow
and $\Phi_{t_{f}}^{t_{0}}=\left(\Phi_{t_{0}}^{t_{f}}\right)^{-1}$. 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 $\mathbb{D}$,

$f_{X}(x)\approx\Psi^{T}(x)a,\qquad\qquad f_{Y}(y)\approx\Psi^{T}(y)\tilde{a}$

then the objective is approximated as

$g(f_{X},f_{Y})\approx\frac{1}{m}\sum_{i=1}^{m}a^{T}\Psi(x_{i})\Psi^{T}(y_{i})\tilde{a}=a^{T}A\tilde{a}$
(26)

where $A=\frac{1}{m}\Psi_{X}\Psi_{Y}^{T}$. If we impose a scale by requiring
that $a^{T}a=\tilde{a}^{T}\tilde{a}=1$, then this maximization can be solved
by singular value decomposition, with the optimal $a$ and $\tilde{a}$ given by
left and right singular vectors, respectively, as shown in Refs. Froyland,
Santitissadeekorn, and Monahan, 2010; Williams, Rypina, and Rowley, 2015. 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 Williams, Rypina, and Rowley (2015).

##  III Control formulation

In Sec. II.3, it was shown that the problem of steering a density $\rho$ 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
$\hat{\rho}$ can be interpreted as a lifted state. Then, if the first and
second raw moments are taken to be the relevant output,

$y=\begin{bmatrix}m_{1}^{1}&m_{1}^{2}&m_{2}^{11}&m_{2}^{22}&m_{2}^{12}\end{bmatrix}^{T}$
(27)

this can be expressed linearly in the lifted state, $y=C\hat{\rho}$, 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

$\displaystyle\min_{u_{1},u_{2},\dots,u_{H-1}}$
$\displaystyle\sum_{t=1}^{H-1}l(\hat{\rho}_{t},u_{t})+l_{H}(\hat{\rho}_{H})$
(28a) $\displaystyle\mathrm{s.t.}\qquad$
$\displaystyle\hat{\rho}_{t+1}=F(\hat{\rho}_{t},u_{t})$ (28b) $\displaystyle
y_{t}=C\hat{\rho}_{t}$ (28c)

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 in-horizon and terminal cost
functions $l$ and $l_{H}$ of the following quadratic forms

$\displaystyle l(\hat{\rho}_{t},u_{t})$
$\displaystyle=(y_{t}-y^{\mathrm{ref}}_{t})^{T}S(y_{t}-y^{\mathrm{ref}}_{t})+u_{t}^{T}Ru_{t}$
(29a) $\displaystyle l_{H}(\hat{\rho}_{H})$
$\displaystyle=(y_{H}-y^{\mathrm{ref}}_{H})^{T}S_{H}(y_{H}-y^{\mathrm{ref}}_{H})$
(29b)

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
$\hat{\rho}_{t}$, this cost can be rewritten as a quadratic cost in terms of
$\hat{\rho}_{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
Hofer and Tibken (1988). For this reason, we solve the optimal control problem
using differential dynamic programming (DDP) Tassa, Erez, and Todorov (2012);
Yakowitz and Rutherford (1984), which is closely related to the method of
iterative LQR. We briefly recount the primary steps of this algorithm below.

Figure 1: Controlled transport of a distribution of fluid particles by two
micro-rotors fixed at $(-1,0)$ and $(1,0)$ from an initial density
$\rho(x(0))=\mathcal{N}((1,1),0.025I_{2})$ 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.

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(\hat{\rho}_{t},t)$ at time $t$ as,

$V(\hat{\rho}_{t},t)=\min_{u_{t}}[l(\hat{\rho}_{t},u_{t})+V(\hat{\rho}_{t+1},t+1)]$
(30)

which expresses the optimal cost-to-go from $\hat{\rho}_{t}$, where
$V(\hat{\rho}_{H},H)=l_{f}(\hat{\rho}_{H})$. Denote by
$Q(\delta\hat{\rho},\delta u)$ the change in the value function due to
applying change in control input $\delta u$ about the nominal trajectory and
consider its quadratic approximation

$\begin{split}Q(\delta\hat{\rho},\delta
u)\approx&\>Q_{\hat{\rho}}\delta\hat{\rho}+Q_{u}^{T}\delta
u+\delta\hat{\rho}^{T}Q_{\hat{\rho}u}\delta u\\\\[4.30554pt]
&+\frac{1}{2}\delta\hat{\rho}^{T}Q_{\hat{\rho}\hat{\rho}}\delta\hat{\rho}+\frac{1}{2}\delta
u^{T}Q_{uu}\delta u\end{split}$ (31)

where these derivatives are given by

$\begin{split}Q_{\hat{\rho}}&=l_{\hat{\rho}}+F_{\hat{\rho}}^{T}V_{\hat{\rho}}^{\prime}\\\
Q_{u}&=l_{u}+F_{u}^{T}V_{\hat{\rho}}^{\prime}\\\
Q_{\hat{\rho}\hat{\rho}}&=l_{\hat{\rho}\hat{\rho}}+F_{\hat{\rho}}^{T}V_{\hat{\rho}\hat{\rho}}^{\prime}F_{\hat{\rho}}+V_{\hat{\rho}}^{\prime}\cdot
F_{\hat{\rho}\hat{\rho}}\\\
Q_{uu}&=l_{uu}+F_{u}^{T}V_{\hat{\rho}\hat{\rho}}^{\prime}F_{u}+V_{\hat{\rho}}^{\prime}\cdot
F_{uu}\\\
Q_{\hat{\rho}u}&=l_{\hat{\rho}u}+F_{\hat{\rho}}^{T}V_{\hat{\rho}\hat{\rho}}^{\prime}F_{u}+V_{\hat{\rho}}^{\prime}\cdot
F_{\hat{\rho}u}\\\ \end{split}$

where the notation $(\cdot)^{\prime}$ 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 $\delta u$

$\delta u^{*}=\arg\min_{\delta u}Q(\delta\hat{\rho},\delta
u)=-Q_{uu}^{-1}\left(Q_{u}+Q_{u\hat{\rho}}\delta\hat{\rho}\right)$ (32)

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 $\alpha$ to update the policy, that is

$u_{\text{new}}=u-\alpha
Q_{uu}^{-1}Q_{u}-Q_{uu}^{-1}Q_{u\hat{\rho}}\delta\hat{\rho}$

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. Tassa, Erez, and Todorov (2012); Yakowitz and Rutherford
(1984).

Figure 2: Rotor controls and particle distribution moments for transport by
rotors in free space as shown in Fig. 1. Left, Top: Rotor strengths for the
left rotor ($\gamma_{L}$) and right rotor ($\gamma_{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.

##  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 Pozrikidis (1992), whose stream function is given by

$\psi(\mathbf{x})=-\gamma\log|\mathbf{x}-\mathbf{x}_{r}|$ (33)

where $\mathbf{x}=(x,y)$ is the position a point in the fluid,
$\mathbf{x}_{r}=(x_{r},y_{r})$ is the position of the rotlet, and $\gamma$ is
the strength of the rotlet. Physically, $\gamma$ 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:

$\mathbf{u}(\mathbf{x})=-\sum_{i=1}^{n_{r}}\left(\gamma_{i}\hat{k}\times\frac{~{}\mathbf{x}-\mathbf{x}_{i}~{}}{r_{i}^{2}}\right)$
(34)

where $\mathbf{x}_{i}$ is the location of the $i$-th rotlet and
$r_{i}=|x-x_{i}|$. Clearly, this results in a flow with a singularity at
$\mathbf{x}_{r}$, circular streamlines around the singularity with
counterclockwise flow for positive $\gamma$, 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 $\gamma_{L}$
and $\gamma_{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
$\Sigma(0)=0.025I_{2}$, where $I_{2}$ is the $2\times 2$ identity matrix. That
is, $\rho(x(0))=\mathcal{N}\big{(}(1,1),0.025\,I_{2}\big{)}$. 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

$m_{2}^{ij}=\sigma^{ij}+m_{1}^{i}m_{1}^{j}$ (35)

where the $\sigma^{ij}$ are the elements of the covariance matrix $\Sigma$, 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 $\Delta t=0.005$. The Perron-Frobenius
operators are computed using Eq. 14 with this trajectory data and a $25\times
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.3.

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 $\gamma_{L}$ and $\gamma_{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’).

###  IV.1 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.4 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\times 50$ grid
over $[-2,2]^{2}$ and the constant function, $\psi=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=\frac{1}{m}\Psi_{X}\Psi_{Y}^{T}$. 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, $\rho(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}$,

$\bar{g}(f_{X},f_{Y})=\frac{1}{m}\sum_{i=1}^{m}\mathrm{sign}(f_{X}(x_{i}))\mathrm{sign}(f_{Y}(y_{i}))$
(36)

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
$\bar{g}=0.9904$.

$t=0.00$ |  $t=1.25$ |  $t=2.50$ |  $t=3.75$ |  $t=5.00$
---|---|---|---|---

Figure 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 $f_{X}$, as
approximated by the left singular vector of
$A=\frac{1}{m}\Psi_{X}\Psi_{Y}^{T}$. The black line in the first image depicts
the $f_{X}=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.

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 Inanc, Shadden, and Marsden (2005); Senatore and Ross (2008);
Ramos  _et al._ (2018); Krishna, Song, and Brunton (2022) 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 vehicle, 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 time-varying 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 Inanc, Shadden, and Marsden (2005);
Krishna, Song, and Brunton (2022), 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 $\rho(x(0))=\mathcal{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
$\bar{g}(f_{X}^{2},f_{Y}^{2})=0.9978$ and
$\bar{g}(f_{X}^{3},f_{Y}^{3})=0.9928$ where $f_{X}^{2},f_{Y}^{2}$ and
$f_{X}^{3},f_{Y}^{3}$ refer to the partition functions given by the second and
third singular vectors, respectively.

Figure 4: Finite-time coherent sets produced by the optimal rotor-driven flow
field with initial density $\rho(x(0))=\mathcal{N}((-0.5,1),0.025I_{2})$ .
Figs. (a) and (b) show the $10^{4}$ data points at the initial and final time,
colored according to the 3nd left singular vector of $A$, along with the
$f_{X}=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.

##  V Transport by rotors near an infinite plane wall

Figure 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
$\rho(x(0))=\mathcal{N}((1,1),0.025I_{2})$ 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.

(a)

(b)

(c)

Figure 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, $J_{w}$ normalized by the cost for the
free space case, $J_{0}$.

Figure 7: Finite-time coherent sets for the case shown in Fig. 5 with an
infinite plane wall at $y=-1.25$. Left: the initially uniform data set colored
according to the partition function $f_{X}$, as approximated by the left
singular vector of $A=\frac{1}{m}\Psi_{X}\Psi_{Y}^{T}$. Right: The same points
at the final time. The black line shows the $f_{X}=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.

For the case of a rotlet located at a point $\mathbf{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 Pozrikidis (1992);
Ranger (1980)

$\begin{split}\psi_{w}(\mathbf{x})&=\gamma\Bigg{(}-\log|\mathbf{x}-\mathbf{x}_{r}|+\log|\mathbf{x}-\mathbf{x}_{\mathrm{im}}|\\\
&\qquad\qquad-\frac{2(y-w)(y-y_{\mathrm{im}})}{r_{\mathrm{im}}^{2}}\Bigg{)}\end{split}$
(37)

where
$\mathbf{x}_{\mathrm{im}}=(x_{\mathrm{im}},y_{\mathrm{im}})=\left(x_{r},2w-y_{r}\right)$
is the location of an image singularity which has the effect of making the
fluid velocity vanish at the plane wall. Similarly,
$r_{\mathrm{im}}=|\mathbf{x}-\mathbf{x}_{\mathrm{im}}|$. Therefore, the flow
produced for this case is $\mathbf{u}_{w}=(u_{w},v_{w})$ where

$\displaystyle u_{w}$
$\displaystyle=\gamma\Bigg{(}-\frac{y-y_{r}}{r^{2}}-\frac{y-y_{\mathrm{im}}}{r_{\mathrm{im}}^{2}}-\frac{2(y-w)}{r_{\mathrm{im}}^{2}}$
(38a)
$\displaystyle\qquad\qquad+\frac{4(y-w)(y-y_{\mathrm{im}})^{2}}{r_{\mathrm{im}}^{4}}\Bigg{)}$
$\displaystyle v_{w}$
$\displaystyle=\gamma\Bigg{(}\frac{x-x_{r}}{r^{2}}-\frac{x-x_{r}}{r_{\mathrm{im}}^{2}}$
(38b)
$\displaystyle\qquad\qquad-\frac{4(y-w)(x-x_{\mathrm{im}})(y-y_{\mathrm{im}})}{r_{\mathrm{im}}^{4}}\Bigg{)}$

and the velocity field for multiple rotlets above a plane wall can be found by
summing the individual velocity fields as in Eq. 34.

Figure 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
$\rho(x(0))=\mathcal{N}((1,1),0.025I_{2})$ 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.

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
$\rho(x(0))=\mathcal{N}\big{(}(1,1),0.025\,I_{2}\big{)}$, 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.

(a)

(b)

(c)

Figure 9: Comparison of boundary effects due to an 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, $J_{c}$ normalized by the cost
for the free space case, $J_{0}$.

##  VI Transport by rotors within a circular boundary

For the case of a rotlet positioned at a point $\mathbf{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 Pozrikidis (1992); Ranger
(1980); Meleshko and Aref (1996)

$\begin{split}\psi_{c}=\gamma\Bigg{(}-\log|&\mathbf{x}-\mathbf{x}_{r}|+\log|\mathbf{x}-\mathbf{x}_{\mathrm{im}}|+\log\frac{R_{r}}{a}\\\
&-\frac{1}{2}\left(\frac{R^{2}-a^{2}}{x_{\mathrm{im}}^{2}}\right)\left(\frac{a^{2}}{R_{r}}-\frac{R^{2}}{a^{2}}\right)\Bigg{)}\end{split}$
(39)

where $R=|\mathbf{x}|$ and $R_{r}=|\mathbf{x}_{r}|$ are the radial distances
from the center of the circle to the evalutation point and to the rotlet,
respectively, and
$\mathbf{x}_{\mathrm{im}}=\frac{a^{2}}{R_{r}^{2}}\mathbf{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
$\mathbf{u}_{c}=(u_{c},v_{c})$ where

$\displaystyle u_{c}$
$\displaystyle=\gamma\Bigg{(}-\frac{y-y_{r}}{r^{2}}+\frac{y-y_{\mathrm{im}}}{r_{\mathrm{im}}^{2}}$
(40a)
$\displaystyle\qquad\qquad+\left(\frac{a^{2}}{R_{r}^{2}}-2\frac{R^{2}}{a^{2}}+1\right)\frac{y}{r_{\mathrm{im}}^{2}}$
$\displaystyle\qquad\qquad-\frac{(y-y_{\mathrm{im}})(R^{2}-a^{2})\left(\frac{a^{2}}{R_{r}}-\frac{R^{2}}{a^{2}}\right)}{r_{\mathrm{im}}^{4}}\Bigg{)}$
$\displaystyle v_{c}$
$\displaystyle=\gamma\Bigg{(}\frac{x-x_{r}}{r^{2}}-\frac{x-x_{\mathrm{im}}}{r_{\mathrm{im}}^{2}}$
(40b)
$\displaystyle\qquad\qquad-\left(\frac{a^{2}}{R_{r}^{2}}-2\frac{R^{2}}{a^{2}}+1\right)\frac{x}{r_{\mathrm{im}}^{2}}$
$\displaystyle\qquad\qquad+\frac{(x-x_{\mathrm{im}})(R^{2}-a^{2})\left(\frac{a^{2}}{R_{r}}-\frac{R^{2}}{a^{2}}\right)}{r_{\mathrm{im}}^{4}}\Bigg{)}.$

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 circular
boundary on the transport problem. As before, the rotor positions of $(-1,0)$
and $(1,0)$, initial density of
$\rho(x(0))=\mathcal{N}\big{(}(1,1),0.025\,I_{2}\big{)}$, 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.

Figure 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 $f_{X}$, as approximated by the
left singular vector of $A=\frac{1}{m}\Psi_{X}\Psi_{Y}^{T}$. Right: The same
points at the final time. The black line shows the $f_{X}=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. 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.

$t=0.00$ |  $t=1.25$ |  $t=2.50$ |  $t=3.75$ |  $t=5.00$
---|---|---|---|---

Figure 11: Controlled transport of two distributions of fluid particles by two
micro-rotors fixed at $(-1,0)$ and $(1,0)$ from initial densities
$\rho^{A}(x(0))=\mathcal{N}((1,1),0.025I_{2})$ and
$\rho^{B}(x(0))=\mathcal{N}((b,1),0.025I_{2})$ to a target mean at $(-1,-1)$
(green circle). 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$.
White filled circles indicate the sample means of the two distributions.
Streamlines depict the flow field produced by the rotor control at the instant
shown.

##  VII Transport of two densities

Figure 12: Rotor strengths for controlled transport for varying initial
distributions, $\rho^{B}(x(0))=\mathcal{N}((b,1),0.025I_{2})$ with
$\rho^{A}(x(0))=\mathcal{N}((1,1),0.025I_{2})$, corresponding to the sequences
shown in Fig. 11.

Figure 13: Finite-time coherent sets for the the transport of two
distributions, as shown in Fig. 11. Left column: the initially uniform data
set colored according to the partition function $f_{X}$, as approximated by
the left singular vector of $A=\frac{1}{m}\Psi_{X}\Psi_{Y}^{T}$. Right: The
same points at the final time. The black line shows the $f_{X}=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$.

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
$\hat{\rho}$. Here we consider an augmented state containing the projection
coefficients of the two density functions. Denoting these two density
functions as $\rho^{A}$ and $\rho^{B}$, and their corresponding projection
coefficients by $\hat{\rho}^{A}$ and $\hat{\rho}^{B}$, the augmented state for
this case is $[(\hat{\rho}^{A})^{T},(\hat{\rho}^{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 $\rho^{A}$ and
$\rho^{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,
$\rho^{A}(x(0))=\mathcal{N}\big{(}(1,1),0.025\,I_{2}\big{)}$, and consider a
second distribution starting from an initial density of
$\rho^{B}(x(0))=\mathcal{N}\big{(}(b,1),0.025\,I_{2}\big{)}$, 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_{1}^{A}(5)=m_{1}^{B}(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.

Fig. 11 shows snapshots from the evolution of the particle distributions for
the flow induced by the rotors controlled using the strengths determined from
the DDP algorithm for four different initial distributions $\rho^{B}$ with the
initial $x$-coordinate of the mean being $b=0.5$, $b=0.0$, $b=-0.5$, and
$b=-1.0$, respectively, on the rows. It can be seen that the flow produced in
the case where $b=0.5$ is qualitatively similar to the case of controlling the
density $\rho^{A}$ alone, as was considered in Sec. IV. For the next two cases
of $b=0$ and $b=-0.5$, we see that as the initial distribution of $\rho^{B}$
starts farther from the initial distribution of $\rho^{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 $\rho^{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 $\rho^{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. 12.

Fig. 13 shows the coherent sets at the initial and final time for the cases
shown in Fig. 11 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 robotsBuzhardt and Tallapragada (2019), 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 Chen, Georgiou, and Pavon (2021). 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.

## References

  * Iverson and Garimella (2008) B. D. Iverson and S. V. Garimella, “Recent advances in microscale pumping technologies: a review and evaluation,” Microfluidics and nanofluidics 5, 145–174 (2008).
  * Zhou  _et al._ (2016) C. Zhou, H. Zhang, Z. Li, and W. Wang, “Chemistry pumps: a review of chemically powered micropumps,” Lab on a Chip 16, 1797–1811 (2016).
  * Xu  _et al._ (2017) T. Xu, W. Gao, L.-P. Xu, X. Zhang, and S. Wang, “Fuel-free synthetic micro-/nanomachines,” Advanced Materials 29, 1603250 (2017).
  * Wang and Fu (2018) Y.-N. Wang and L.-M. Fu, “Micropumps and biomedical applications–a review,” Microelectronic Engineering 195, 121–138 (2018).
  * Aref (1984) H. Aref, “Stirring by chaotic advection,” Journal of fluid mechanics 143, 1–21 (1984).
  * Aref  _et al._ (2017) H. Aref, J. R. Blake, M. Budišić, S. S. Cardoso, J. H. Cartwright, H. J. Clercx, K. El Omari, U. Feudel, R. Golestanian, E. Gouillart, _et al._ , “Frontiers of chaotic advection,” Reviews of Modern Physics 89, 025007 (2017).
  * Ottino (1989) J. M. Ottino,  _The kinematics of mixing: stretching, chaos, and transport_ , Vol. 3 (Cambridge university press, 1989).
  * Chwang and Wu (1975) A. T. Chwang and T. Y.-T. Wu, “Hydromechanics of low-reynolds-number flow. part 2. singularity method for stokes flows,” Journal of Fluid mechanics 67, 787–815 (1975).
  * Pozrikidis (1992) C. Pozrikidis,  _Boundary integral and singularity methods for linearized viscous flow_ (Cambridge university press, 1992).
  * Buzhardt and Tallapragada (2023) J. Buzhardt and P. Tallapragada, “Controlled density transport using perron frobenius generators,” arXiv preprint arXiv:2304.13829 (2023).
  * Meleshko and Aref (1996) V. Meleshko and H. Aref, “A blinking rotlet model for chaotic advection,” Physics of Fluids 8, 3215–3217 (1996).
  * Van Der Woude  _et al._ (2007) D. Van Der Woude, H. Clercx, G. Van Heijst, and V. Meleshko, “Stokes flow in a rectangular cavity by rotlet forcing,” Physics of Fluids 19, 083602 (2007).
  * Lushi and Vlahovska (2015) E. Lushi and P. M. Vlahovska, “Periodic and chaotic orbits of plane-confined micro-rotors in creeping flows,” Journal of Nonlinear Science 25, 1111–1123 (2015).
  * Tallapragada and Sudarsanam (2019) P. Tallapragada and S. Sudarsanam, “Chaotic advection and mixing by a pair of microrotors in a circular domain,” Physical Review E 100, 062207 (2019).
  * Delmotte (2019) B. Delmotte, “Hydrodynamically bound states of a pair of microrollers: A dynamical system insight,” Physical Review Fluids 4, 044302 (2019).
  * Mathew  _et al._ (2007) G. Mathew, I. Mezić, S. Grivopoulos, U. Vaidya, and L. Petzold, “Optimal control of mixing in stokes fluid flows,” Journal of Fluid Mechanics 580, 261–281 (2007).
  * D’Alessandro, Dahleh, and Mezic (1999) D. D’Alessandro, M. Dahleh, and I. Mezic, “Control of mixing in fluid flow: A maximum entropy approach,” IEEE Transactions on Automatic Control 44, 1852–1863 (1999).
  * Cortelezzi, Adrover, and Giona (2008) L. Cortelezzi, A. Adrover, and M. Giona, “Feasibility, efficiency and transportability of short-horizon optimal mixing protocols,” Journal of Fluid Mechanics 597, 199–231 (2008).
  * Lin, Thiffeault, and Doering (2011) Z. Lin, J.-L. Thiffeault, and C. R. Doering, “Optimal stirring strategies for passive scalar mixing,” Journal of Fluid Mechanics 675, 465–476 (2011).
  * Hassanzadeh, Chini, and Doering (2014) P. Hassanzadeh, G. P. Chini, and C. R. Doering, “Wall to wall optimal transport,” Journal of fluid mechanics 751, 627–662 (2014).
  * Zhang and Balasuriya (2020) L. Zhang and S. Balasuriya, “Controlling trajectories globally via spatiotemporal finite-time optimal control,” SIAM Journal on Applied Dynamical Systems 19, 1609–1632 (2020).
  * Lasota and Mackey (1994) A. Lasota and M. C. Mackey,  _Chaos, Fractals, and Noise : Stochastic Aspects of Dynamics_ (Springer, 1994).
  * Otto and Rowley (2021) S. E. Otto and C. W. Rowley, “Koopman operators for estimation and control of dynamical systems,” Annual Review of Control, Robotics, and Autonomous Systems 4, 59–87 (2021).
  * Kaiser, Kutz, and Brunton (2020) E. Kaiser, J. N. Kutz, and S. L. Brunton, “Data-driven approximations of dynamical systems operators for control,” The Koopman Operator in Systems and Control: Concepts, Methodologies, and Applications , 197–234 (2020).
  * Rowley  _et al._ (2009) C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S. Henningson, “Spectral analysis of nonlinear flows,” Journal of fluid mechanics 641, 115–127 (2009).
  * Mezić (2013) I. Mezić, “Analysis of fluid flows via spectral properties of the koopman operator,” Annual Review of Fluid Mechanics 45, 357–378 (2013).
  * Froyland, González-Tokman, and Watson (2016) G. Froyland, C. González-Tokman, and T. M. Watson, “Optimal mixing enhancement by local perturbation,” SIAM Review 58, 494–513 (2016).
  * Froyland and Santitissadeekorn (2017) G. Froyland and N. Santitissadeekorn, “Optimal mixing enhancement,” SIAM Journal on Applied Mathematics 77, 1444–1470 (2017).
  * Klünker, Padberg-Gehle, and Thiffeault (2022) A. Klünker, K. Padberg-Gehle, and J.-L. Thiffeault, “Open-flow mixing and transfer operators,” Philosophical Transactions of the Royal Society A 380, 20210028 (2022).
  * Sinha, Vaidya, and Rajaram (2016) S. Sinha, U. Vaidya, and R. Rajaram, “Operator theoretic framework for optimal placement of sensors and actuators for control of nonequilibrium dynamics,” Journal of Mathematical Analysis and Applications 440, 750–772 (2016).
  * Brockett (2007) R. W. Brockett, “Optimal control of the Liouville equation,” AMS IP Studies in Advanced Mathematics 39, 23 (2007).
  * Brockett (2012) R. W. Brockett, “Notes on the control of the Liouville equation. control of partial differential equations, lecture notes in mathematics,” (2012).
  * Grover and Elamvazhuthi (2018) P. Grover and K. Elamvazhuthi, “Optimal perturbations for nonlinear systems using graph-based optimal transport,” Communications in Nonlinear Science and Numerical Simulation 59, 197–215 (2018).
  * Elamvazhuthi and Grover (2018) K. Elamvazhuthi and P. Grover, “Optimal transport over nonlinear systems via infinitesimal generators on graphs.” Journal of Computational Dynamics 5 (2018).
  * Cvitanovic  _et al._ (2005) P. Cvitanovic, R. Artuso, R. Mainieri, G. Tanner, G. Vattay, N. Whelan, and A. Wirzba, “Chaos: classical and quantum,” ChaosBook. org (Niels Bohr Institute, Copenhagen 2005) 69, 25 (2005).
  * Ulam (1960) S. M. Ulam, “A collection of mathematical problems,” New York 29 (1960).
  * Dellnitz, Froyland, and Junge (2001) M. Dellnitz, G. Froyland, and O. Junge, “The algorithms behind GAIO—set oriented numerical methods for dynamical systems,” in  _Ergodic theory, analysis, and efficient simulation of dynamical systems_ (Springer, 2001) pp. 145–174.
  * Klus, Koltai, and Schütte (2016) S. Klus, P. Koltai, and C. Schütte, “On the numerical approximation of the Perron-Frobenius and Koopman operator,” Journal of Computational Dynamics 3, 51–79 (2016).
  * Williams, Kevrekidis, and Rowley (2015) M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, “A data–driven approximation of the koopman operator: Extending dynamic mode decomposition,” Journal of Nonlinear Science 25, 1307–1346 (2015).
  * Huang and Vaidya (2018) B. Huang and U. Vaidya, “Data-driven approximation of transfer operators: Naturally structured dynamic mode decomposition,” in  _2018 Annual American Control Conference (ACC)_ (IEEE, 2018) pp. 5659–5664.
  * Goswami, Thackray, and Paley (2018) D. Goswami, E. Thackray, and D. A. Paley, “Constrained Ulam dynamic mode decomposition: Approximation of the Perron-Frobenius operator for deterministic and stochastic systems,” IEEE control systems letters 2, 809–814 (2018).
  * Korda and Mezić (2018) M. Korda and I. Mezić, “Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control,” Automatica 93, 149–160 (2018).
  * Goswami and Paley (2017) D. Goswami and D. A. Paley, “Global bilinearization and controllability of control-affine nonlinear systems: A koopman spectral approach,” in  _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_ (IEEE, 2017) pp. 6107–6112.
  * Klus  _et al._ (2020) S. Klus, F. Nüske, S. Peitz, J.-H. Niemann, C. Clementi, and C. Schütte, “Data-driven approximation of the Koopman generator: Model reduction, system identification, and control,” Physica D: Nonlinear Phenomena 406, 132416 (2020).
  * Peitz, Otto, and Rowley (2020) S. Peitz, S. E. Otto, and C. W. Rowley, “Data-driven model predictive control using interpolated koopman generators,” SIAM Journal on Applied Dynamical Systems 19, 2162–2193 (2020).
  * Bruder, Fu, and Vasudevan (2021) D. Bruder, X. Fu, and R. Vasudevan, “Advantages of bilinear koopman realizations for the modeling and control of systems with unknown dynamics,” IEEE Robotics and Automation Letters 6, 4369–4376 (2021).
  * Folkestad and Burdick (2021) C. Folkestad and J. W. Burdick, “Koopman NMPC: Koopman-based learning and nonlinear model predictive control of control-affine systems,” in  _2021 IEEE International Conference on Robotics and Automation (ICRA)_ (IEEE, 2021) pp. 7350–7356.
  * Happel and Brenner (1983) J. Happel and H. Brenner,  _Low Reynolds number hydrodynamics: with special applications to particulate media_ , Vol. 1 (Springer Science & Business Media, 1983).
  * Dellnitz and Junge (1999) M. Dellnitz and O. Junge, “On the approximation of complicated dynamical behavior,” SIAM Journal on Numerical Analysis 36, 491–515 (1999).
  * Froyland, Santitissadeekorn, and Monahan (2010) G. Froyland, N. Santitissadeekorn, and A. Monahan, “Transport in time-dependent dynamical systems: Finite-time coherent sets,” Chaos: An Interdisciplinary Journal of Nonlinear Science 20, 043116 (2010).
  * Williams, Rypina, and Rowley (2015) M. O. Williams, I. I. Rypina, and C. W. Rowley, “Identifying finite-time coherent sets from limited quantities of lagrangian data,” Chaos: An Interdisciplinary Journal of Nonlinear Science 25 (2015).
  * Allshouse and Peacock (2015) M. R. Allshouse and T. Peacock, “Lagrangian based methods for coherent structure detection,” Chaos: An Interdisciplinary Journal of Nonlinear Science 25 (2015).
  * Hadjighasem  _et al._ (2017) A. Hadjighasem, M. Farazmand, D. Blazevski, G. Froyland, and G. Haller, “A critical comparison of lagrangian methods for coherent structure detection,” Chaos: An Interdisciplinary Journal of Nonlinear Science 27 (2017).
  * Hofer and Tibken (1988) E. Hofer and B. Tibken, “An iterative method for the finite-time bilinear-quadratic control problem,” Journal of optimization Theory and applications 57, 411–427 (1988).
  * Tassa, Erez, and Todorov (2012) Y. Tassa, T. Erez, and E. Todorov, “Synthesis and stabilization of complex behaviors through online trajectory optimization,” in  _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_ (IEEE, 2012) pp. 4906–4913.
  * Yakowitz and Rutherford (1984) S. Yakowitz and B. Rutherford, “Computational aspects of discrete-time optimal control,” Applied Mathematics and Computation 15, 29–45 (1984).
  * Inanc, Shadden, and Marsden (2005) T. Inanc, S. C. Shadden, and J. E. Marsden, “Optimal trajectory generation in ocean flows,” in  _Proceedings of the 2005, American Control Conference, 2005._ (IEEE, 2005) pp. 674–679.
  * Senatore and Ross (2008) C. Senatore and S. D. Ross, “Fuel-efficient navigation in complex flows,” in  _2008 American Control Conference_ (IEEE, 2008) pp. 1244–1248.
  * Ramos  _et al._ (2018) A. Ramos, V. García-Garrido, A. Mancho, S. Wiggins, J. Coca, S. Glenn, O. Schofield, J. Kohut, D. Aragon, J. Kerfoot, _et al._ , “Lagrangian coherent structure assisted path planning for transoceanic autonomous underwater vehicle missions,” Scientific reports 8, 4575 (2018).
  * Krishna, Song, and Brunton (2022) K. Krishna, Z. Song, and S. L. Brunton, “Finite-horizon, energy-efficient trajectories in unsteady flows,” Proceedings of the Royal Society A 478, 20210255 (2022).
  * Ranger (1980) K. Ranger, “Eddies in two dimensional stokes flow,” International Journal of Engineering Science 18, 181–190 (1980).
  * Buzhardt and Tallapragada (2019) J. Buzhardt and P. Tallapragada, “Dynamics of groups of magnetically driven artificial microswimmers,” Physical Review E 100, 033106 (2019).
  * Chen, Georgiou, and Pavon (2021) Y. Chen, T. T. Georgiou, and M. Pavon, “Controlling uncertainty,” IEEE Control Systems Magazine 41, 82–94 (2021).