# Controlled density transport using Perron Frobenius generators

Jake Buzhardt and Phanindra Tallapragada The authors are with the Department
of Mechanical Engineering, Clemson University, Clemson, SC, USA.
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We consider the problem of the transport of a density of states from an
initial state distribution to a desired final state distribution through a
dynamical system with actuation. In particular, we consider the case where the
control signal is a function of time, but not space; that is, the same
actuation is applied at every point in the state space. This is motivated by
several problems in fluid mechanics, such as mixing and manipulation of a
collection of particles by a global control input such as a uniform magnetic
field, as well as by more general control problems where a density function
describes an uncertainty distribution or a distribution of agents in a multi-
agent system. We formulate this problem using the generators of the Perron-
Frobenius operator associated with the drift and control vector fields of the
system. By considering finite-dimensional approximations of these operators,
the density transport problem can be expressed as a control problem for a
bilinear system in a high-dimensional, lifted state. With this system, we
frame the density control problem as a problem of driving moments of the
density function to the moments of a desired density function, where the
moments of the density can be expressed as an output which is linear in the
lifted state. This output tracking problem for the lifted bilinear system is
then solved using differential dynamic programming, an iterative trajectory
optimization scheme.

##  I Introduction

In this paper, we consider the problem of controlled density transport, where
given an initial distribution of states specified by a density function, we
seek to determine a control sequence to drive this initial distribution to a
desired final distribution. We consider the case where a common control signal
is applied to the entire distribution of states. This differs from the usual
formulation of swarm control and optimal transport problems, where typically
each agent can select a control input independently, making the control signal
a function of the states and time. This problem of density transport is
motivated by problems of manipulation of a large collection of agents using a
uniform control signal [1, 2]. The transport of density also has relevance to
the propagation of an uncertainty distribution arising due to uncertainty in
the initial state or of a model parameter through an otherwise deterministic
control system (see, e.g. [3, 4, 5]). We formulate and solve this problem
using an operator theoretic approach, specifically using the generator of the
Perron-Frobenius operator.

In recent years the operator theoretic approach to dynamical systems and
control has gained significant research attention [6, 7, 8]. A dynamical
system can be framed in terms of such an operator either by considering the
evolution of observable functions of the state using the Koopman operator or
by considering the evolution of densities of states using the Perron-Frobenius
operator [9]. The interest in these approaches is primarily due to the fact
that these operators allow for a linear, although typically infinite
dimensional representation of a nonlinear system. The linearity of these
operators is useful from an analytical perspective, as it allows for the use
of linear systems techniques such as the analysis of eigenvalues and
eigenfunctions, but also from a computational perspective, as in many cases a
useful approximation for these operators can be found by considering a finite
dimensional approximation in which the operator is represented as a matrix
acting on coordinates corresponding to a finite set of a set of dictionary
functions [10, 11, 12].

In applications in control systems, much of the recent work has been on
developing methods involving the Koopman operator [13, 8], as the
transformation to a space of observable functions can be viewed as a nonlinear
change of coordinates which maps the system to a higher dimensional space
where the dynamics are (approximately) linear [14]. This makes the numerical
approximation of the operator particularly amenable to linear control methods,
such as the linear quadratic regulator (LQR) and model predictive control
(MPC) [14, 15, 16]. On the other hand, the Perron-Frobenius operator
propagates densities of states forward in time along trajectories of the
system, which can have multiple interpretations in the controlled setting. For
example, the Perron-Frobenius operator and the Liouville equation, the related
PDE formulation, have been used to determine controls for agents in an
ensemble or swarm formulation [17, 18]. Such formulations are closely related
to optimal transport problems which also involve driving an initial
distribution to a desired final distribution (see, e.g., [18, 19, 20]).
Formulations involving the Perron-Frobenius operator have also been used in
the context of fluid flows to study the transport of distributions of fluid
particles and to detect invariant or almost invariant sets [21, 22].

Our approach involves first obtaining a finite dimensional approximation of
the Perron-Frobenius generators associated with the drift and control vector
fields of the system, which allow us to represent the density transport
dynamics as a bilinear system in a lifted state. With this system, we frame
the density control problem as a problem of driving moments of the density
function to the moments of a desired density function, where the moments of
the density can be expressed as an output which is linear in the lifted state.
This output tracking problem for the lifted bilinear system is then solved
using differential dynamic programming (DDP), an iterative trajectory
optimization scheme.

##  II Preliminaries

Consider first the autonomous dynamical system on a measure space
$(\mathbb{X}\subset\mathbb{R}^{n},\mathcal{A},\mu)$ with a $\sigma$-algebra
$\mathcal{A}$ on $\mathbb{X}$ and $\mu$ a measure on
$(\mathbb{X},\mathcal{A})$,

$\dot{x}=f(x)$ (1)

and denote the associated time-$t$ flow from an initial state $x_{0}$ as
$\Phi^{t}(x_{0})$, where $x\in\mathbb{X}$ is the state. The Perron-Frobenius
operator $\mathcal{P}^{t}:L^{1}(\mathbb{X})\mapsto L^{1}(\mathbb{X})$
associated with the flow map $\Phi^{t}$ is defined as

$\int_{\mathbb{A}}\left[\mathcal{P}^{t}\,\rho\right](x)\,dx=\int_{(\Phi^{t})^{-1}(\mathbb{A})}\rho(x)\,dx$
(2)

for any $\mathbb{A}\in\mathcal{A}$, assuming that the relevant measure $\mu$
is absolutely continuous with respect to the Lebesgue measure and can thus be
expressed in terms of a density $\rho$ (i.e., $d\mu(x)=\mu(dx)=\rho(x)dx$). It
can be shown that the family of these operators $\\{\mathcal{P}^{t}\\}_{t\geq
0}$ form a semigroup, (see [9]). The generator of this semigroup is known as
the Liouville operator, denoted $\mathcal{L}$, or Perron-Frobenius generator
and expresses the deformation of the density $\rho$ under infinitesimal action
of the operator $\mathcal{P}^{t}$ [23, 9]. That is,

$\frac{d\rho}{dt}=\mathcal{L}\rho=-\nabla_{x}\cdot(\rho f)$ (3)

Alternatively, the action of the generator can be written in terms of the
Perron-Frobenius operator 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$ (4)

where $\mathcal{I}$ is the identity operator.

###### Lemma 1

Suppose the Liouville operator associated with a vector field
$f_{1}:\mathbb{X}\mapsto\mathbb{R}^{n}$ is denoted by $\mathcal{L}_{1}$ and
the Liouville operator associated with the vector field
$f_{2}:\mathbb{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}$.

###### Proof:

The proof is a direct consequence of Eq. 3. Suppose $f(x)=f_{1}(x)+f_{2}(x)$.
Then
$\mathcal{L}\rho=-\nabla_{x}\cdot(\rho(f_{1}+f_{2}))=-\nabla_{x}\cdot(\rho
f_{1})-\nabla_{x}\cdot(\rho f_{2})=(\mathcal{L}_{1}+\mathcal{L}_{2})\rho.$ ∎

The Koopman operator $\mathcal{K}^{t}:L^{\infty}(\mathbb{X})\mapsto
L^{\infty}(\mathbb{X})$ propagates observable functions forward in time along
trajectories of the system and is defined as

$[\mathcal{K}^{t}h](x)=[h\circ\Phi^{t}](x)$ (5)

where $h(x)$ is an observable. The Koopman and Perron-Frobenius operators are
adjoint to one another,

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

Figure 1: Moment propagation of proposed method for a Duffing oscillator with
sinusoidal forcing. Red points show trajectories from initial conditions
sampled from the initial density, $\rho(x(0))\sim\mathcal{N}([-0.5;1],0.05I)$.
Red circle and red ellipse show the sample mean and 2$\sigma$ sample
covariance ellipse, respectively. The black circle and black ellipse are the
predicted mean and $2\sigma$ covariance ellipse.

##  III Numerical approximation of the Perron-Frobenius operator and generator

Here, we implement extended dynamic mode decomposition (EDMD) [11] for the
computation of the Perron-Frobenius operator, which we outline below, largely
following [12].

We begin by selecting a dictionary $\mathbb{D}$ of $k$ scalar-valued basis
functions, $\mathbb{D}=\\{\psi_{1},\psi_{2},\dots,\psi_{k}\\}$, where
$\psi_{i}:\mathbb{X}\mapsto\mathbb{R}$ for $i=1,\dots,k$, and denote by
$\mathbb{V}$ the function space spanned by the elements of $\mathbb{D}$. We
then collect trajectory data 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})$.

We then approximate the observable function $h$ and density $\rho$ in Eq. 6 by
their projections onto $\mathbb{V}$. That is,

$\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:\mathbb{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, we have

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

Then replacing $[\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}$ gives a least-squares problem for the matrix $P$

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

where $\Psi_{X}$,$\Psi_{Y}\in\mathbb{R}^{k\times m}$ are matrices with columns
formed by evaluating $\Psi$ 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}$ (13)

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 Perron Frobenius generator can be approximated
based on the limit definition of the generator in Eq. 4. as

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

where $I_{k}$ is the $k\times k$ identity matrix. Just as the matrix operator
$P$ approximates the propagation of a density function $\rho$ by advancing the
projection coordinates $\hat{\rho}$ forward for a finite time, the
approximation of the generator allows us to approximate 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}\,.$ (15)

###  III-A Extension to controlled systems

In the context of applying the Koopman operator to control systems, several
recent works have noted the usefulness of formulating the problem in terms of
the Koopman generator, rather than the Koopman operator itself (see [24, 25,
26] and others), which typically results in a lifted system that is bilinear
in the control and lifted state. This approach allows for a better
approximation of the effects of control, especially for systems in control-
affine form

$\dot{x}=f(x)+\sum_{i=1}^{n_{c}}g_{i}(x)u_{i}$ (16)

as it expresses the effect of the control vector fields $g_{i}$ in a way that
is also dependent on the lifted state. Here we apply a similar approach to the
density transport problem, expressed in terms of the Perron-Frobenius
generator. As shown in Ref. [26] 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. 3. This leads to density transport dynamics of the
following form

$\frac{d}{dt}\rho(x)=(\mathcal{L}_{0}\rho)(x)+\sum_{i=1}^{n_{c}}u_{i}(\mathcal{B}_{i}\rho)(x)$
(17)

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}$
(18)

where the matrices $L_{0}$ and $B_{i}$ are the matrix approximations of the
operators in Eq. 17.

###  III-B Propagation of moments

In order to make the problem of controlling a density function using a finite
dimensional control input well-posed, we formulate the problem as a control
problem of a finite number of outputs. In particular, we will describe the
density function, in terms of a finite number of its moments. For a scalar
$x$, recall that the $i^{\text{th}}$ raw moment, $m_{i}$ is defined as
$m_{i}=\int_{\mathbb{X}}x^{i}\rho(x)dx$ and the $i^{\text{th}}$ central
moment, $\mu_{i}$ about the mean $m_{1}$ is
$\mu_{i}=\int_{\mathbb{X}}(x-m_{1})^{i}\rho(x)dx$.

Given a projection of the density as in Eq. 10, the mean is approximated as

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

which simply indicates that the mean of the density can be written as a
summation of the means of the dictionary functions, weighted by the projection
coefficients. For higher moments, if the central moment is considered, it will
be polynomial in the projection coefficients due to its dependence on the
mean, whereas the raw moments remain linear in the projection coefficients.
For this reason, we choose to work with the raw moments, as the central
moments can also be expressed in terms of the raw moments. 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$
(20)

where, again, superscripts $i$ and $j$ are coordinate indices.

Figure 2: Controlled transport through the Duffing system from an initial
density $\rho(x(0))\sim\mathcal{N}([0;1],0.025I)$

toward the equilibrium point at $(0,1)$ (green circle). Red circle and red
ellipse show the sample mean and 2$\sigma$ sample covariance ellipse,
respectively. The black circle and black ellipse are the predicted mean and
$2\sigma$ covariance ellipse.

Figure 3: Moment propagation for a forced Duffing oscillator. Left, top: First
raw moment (mean). Left, bottom: sinusoidal control signal. Right: 2nd raw
moment.

###  III-C Numerical example

To illustrate the ability of the proposed framework to propagate density
functions forward in time, we consider the propagation of an initial density
for a forced Duffing oscillator system, given by

$\frac{d}{dt}\begin{pmatrix}x_{1}\\\
x_{2}\end{pmatrix}=\begin{pmatrix}x_{2}\\\ x_{1}-x_{1}^{3}+u\end{pmatrix}$
(21)

where $u$ is the control input. For the purpose of this simulation, we set
$u(t)=\sin(4\pi t)$ and the prediction results are shown in Figs. 1 and 3. For
the generator calculation, a dictionary of Gaussian radial basis functions is
used where the centers lie on an evenly spaced $30\times 30$ grid ranging from
$-2.5$ to $2.5$ in $x_{1}$ and $x_{2}$. The operators are approximated using
data collected from short time trajectories with $\Delta t=0.005$ for a
$50\times 50$ grid of initial conditions on the same region. The predicted
moment is compared to the sample moment obtained from $1000$ trajectories from
initial conditions sampled according to the initial density. We see that the
moment propagation of the proposed method is good for at least 3 seconds,
which motivates the use of this method in a control formulation, as detailed
in the following sections. Also shown for comparison in Fig. 3 is a linear
prediction, which is computed by propagating the initial Gaussian through a
linearization of Eq. 21, where the linearization is re-computed at each
timestep about the predicted mean, as is commonly done in the a priori
prediction step of an extended Kalman filter.

##  IV Control formulation

We have shown in Sec. III-A and III-B that the problem of steering a density
$\rho$ to a desired density can be expressed as an output tracking problem on
a lifted, bilinear system given by Eq. 18, where the projection coefficients
$\hat{\rho}$ can be interpreted as the lifted state. Then, if the raw moments
are taken to be the relevant output, the output, $y$ is linear in the lifted
state, $y=C\hat{\rho}$, where the elements of the output matrix $C$ are given
by rewriting Eqs. 19, 20 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})$
(22a) $\displaystyle\mathrm{s.t.}\qquad$
$\displaystyle\hat{\rho}_{t+1}=F(\hat{\rho}_{t},u_{t})$ (22b) $\displaystyle
y_{t}=C\hat{\rho}_{t}$ (22c)

where $H$ is the number of timesteps in the time horizon and Eq. 22b
represents the discrete time version of Eq. 18.

In particular, for output tracking, we consider a quadratic cost of the form

$\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}$
(23) $\displaystyle l_{H}(\hat{\rho}_{H})$
$\displaystyle=(y_{H}-y^{\mathrm{ref}}_{H})^{T}S_{H}(y_{H}-y^{\mathrm{ref}}_{H})$
(24)

where $S$ and $R$ are weighting matrices which define the relative penalty on
tracking error and control effort, respectively. This cost is quadratic in
$\hat{\rho}_{t}$ (with a 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
[27]. For this reason, we solve the optimal control problem using differential
dynamic programming (DDP) [28, 29], which is closely related to the method of
iterative LQR.

Figure 4: Controlled transport of fluid particles, driven by two rotlets, or
micro-rotors, in a Stokes flow from an initial density,
$\rho(x(0))\sim\mathcal{N}([1;1],0.025I)$ toward a target mean (green circle)
at $(-1,-1)$. The rotors are located at $(-1,0)$ and $(0,1)$, as indicated by
the black circle-cross. Red circle and red ellipse show the sample mean and
2$\sigma$ sample covariance ellipse, respectively. The black circle and black
ellipse are the predicted mean and $2\sigma$ covariance ellipse. Gray
streamlines indicate the flow field produced by the rotlets at the instant
shown.

##  V Validation and Results

Here we consider two examples of the density control formulation using Perron
Frobenius generators.

###  V-A Example 1: Forced Duffing oscillator

Figure 5: Control of raw moments for the forced Duffing system. The black line
indicates the target (reference). Left, bottom: control from differential
dynamic programming for the generator system. Right: 2nd raw moment. Predicted
values are from the Perron-Frobenius generator computation. True values are
given by the sample moment.  Figure 6: Comparison of the proposed Perron
Frobenius generator DDP with standard DDP for the forced Duffing system.
Shaded region indicates one standard deviation

As a first example, we consider the forced Duffing oscillator of Eq. 21, and
we use DDP to determine a control sequence to steer a Gaussian initial density
$\rho(x(0))\sim\mathcal{N}([1;1],0.025I)$ toward the equilibrium point at
$(1,0)$ over a time horizon of 2s. The generators are approximated using the
same data as described in Sec. III-C and the same timestep of $0.005$s is used
for DDP. In differential dynamic programming, the in-horizon cost weights on
the reference error of the first two moments and control effort are all set to
unity. The terminal cost on the reference error in the moments is set to
$1000$, and only the first two moments are considered. The target raw second
moment is computed from the desired mean, with the desired variance taken to
be zero. The results of this computation are shown in Fig. 2-5.

In Fig. 6, we show a comparison of the performance of the proposed controller,
labelled ‘PF-DDP’ with a standard DDP controller, which computes a control on
the Duffing system directly (rather than a lifted state), with the initial
condition being the mean. That is, the standard DDP controller acts on the
mean as if it were a deterministic initial condition of the system. The
comparison shown is the mean and one standard deviation region from
trajectories from a set of 500 initial conditions sampled from the initial
distribution, with each of the respective controllers applied. We see that the
proposed controller moves the mean of the distribution closer to the
reference, while maintaining nearly the same variance as the standard DDP
controller.

###  V-B Example 2: Rotor-driven Stokes flow

As a second example, we consider the problem of steering a distribution of
fluid particles in a Stokes flow, where the flow is produced by two micro-
rotors. The micro-rotors are modelled as rotlets, the Stokes flow singularity
associated with a point torque in the fluid. For a collection of $N_{r}$
rotlets, this flow is given by

$\frac{d\mathbf{x}}{dt}=\sum_{i=1}^{N_{r}}\frac{T_{i}\times(\mathbf{x}-\mathbf{\bar{x}_{i}})}{\|\mathbf{x}-\mathbf{\bar{x}}_{i}\|^{3}}$
(25)

where $\mathbf{\bar{x}}_{i}$ and $T_{i}$ denote the location and torque of the
$i$-th rotor, respectively. We consider the case where there are two such
rotors lying in the $x_{1}$-$x_{2}$ plane, located at
$\mathbf{\bar{x}_{1}}=(-1,0)$ and $\mathbf{\bar{x}_{2}}=(1,0)$ and the
controls for the problem are taken to be torques $u_{1}=T_{1}$, $u_{2}=T_{2}$,
where the direction of these torques is taken to be normal to the
$x_{1}$-$x_{2}$ plane (in the positive $x_{3}$ direction.

For this example, the task is to drive a Gaussian initial density
$\rho(x(0))\sim\mathcal{N}([1;1],0.025I)$ toward a target mean at $(-1,-1)$
over a time horizon of 2s. Again, we consider a timestep of $\Delta t=0.005$s
for both the computation of the generators and for DDP. The in horizon costs
weights for the mean error are set to 2, while the weights for the second
moment error and control effort are unity. The terminal cost weight on the
mean error is 1000 and the terminal cost weight on the second moment error is
500. Results from this computation are shown in Figs. 4 and 7.

We see in Fig. 7 that DDP yields a control which drives the mean to the target
by first giving a significant counterclockwise torque on the right rotor to
drive the density to a point between the rotors, at which point the left rotor
is initiated to drive the flow with a clockwise torque, pulling the density
mean near to the target.

This example demonstrates an alternative physical interpretation of the
density transport problem, where the density represents a distribution of
fluid particles. This also demonstrates the effectiveness of the proposed
method on a system which is linear in the controls, but in which the control
vector fields are nonlinear.

Figure 7: Control of raw moments for the rotlet-driven Stokes flow. The black
line indicates the target (reference). Left, top: 1st raw moment. Left,
bottom: control inputs $u_{L}$ and $u_{R}$ for the left and right rotor.
Right: 2nd raw moment. Predicted values are from the Perron-Frobenius
generator computation. True values are given by the sample moment.

##  VI Conclusion

In this work, we have studied the problem of transporting density functions of
states through a controlled dynamical system. This problem formulation has
applications both in fluid mechanics and in the control of uncertain systems.
Our approach is based on approximations of the Perron-Frobenius operator,
whereby we show that approximations of this operator and its generator can be
used to model the density transport dynamics as a high-dimensional system
which is bilinear in the lifted state and the control. We demonstrated this
approach on two examples, a forced Duffing system, in which the density can
have the interpretation as an uncertainty in the initial state and on a rotor
driven Stokes flow, in which the density formulation takes on the fluid
mechanics interpretation of describing a distribution of fluid particles.
Future work in these areas could include extending the proposed control
formulation for use in a constrained model predictive control framework for
uncertain systems or by studying the fluid transport by more realistic
biological microswimmers or artificial microrobots [30].

## References

  * [1] A. Becker, G. Habibi, J. Werfel, M. Rubenstein, and J. McLurkin, “Massive uniform manipulation: Controlling large populations of simple robots with a common input signal,” in 2013 IEEE/RSJ international conference on intelligent robots and systems, pp. 520–527, IEEE, 2013.
  * [2] J. Buzhardt and P. Tallapragada, “Dynamics of groups of magnetically driven artificial microswimmers,” Physical Review E, vol. 100, no. 3, p. 033106, 2019.
  * [3] A. Mesbah, S. Streif, R. Findeisen, and R. D. Braatz, “Stochastic nonlinear model predictive control with probabilistic constraints,” in 2014 American control conference, pp. 2413–2419, IEEE, 2014.
  * [4] G. I. Boutselis, Y. Pan, and E. A. Theodorou, “Numerical trajectory optimization for stochastic mechanical systems,” SIAM Journal on scientific computing, vol. 41, no. 4, pp. A2065–A2087, 2019.
  * [5] Y. Aoyama, A. D. Saravanos, and E. A. Theodorou, “Receding horizon differential dynamic programming under parametric uncertainty,” in 2021 60th IEEE Conference on Decision and Control, pp. 3761–3767, IEEE, 2021.
  * [6] M. Budišić, R. Mohr, and I. Mezić, “Applied Koopmanism,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 22, no. 4, p. 047510, 2012.
  * [7] S. L. Brunton, M. Budišić, E. Kaiser, and J. N. Kutz, “Modern Koopman theory for dynamical systems,” SIAM Review, vol. 64, no. 2, pp. 229–340, 2022.
  * [8] S. E. Otto and C. W. Rowley, “Koopman operators for estimation and control of dynamical systems,” Annual Review of Control, Robotics, and Autonomous Systems, vol. 4, pp. 59–87, 2021.
  * [9] A. Lasota and M. C. Mackey, Chaos, Fractals, and Noise : Stochastic Aspects of Dynamics.  Springer, 1994.
  * [10] M. Dellnitz and O. Junge, “On the approximation of complicated dynamical behavior,” SIAM Journal on Numerical Analysis, vol. 36, no. 2, pp. 491–515, 1999.
  * [11] 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, vol. 25, pp. 1307–1346, 2015\.
  * [12] S. Klus, P. Koltai, and C. Schütte, “On the numerical approximation of the Perron-Frobenius and Koopman operator,” Journal of Computational Dynamics, vol. 3, no. 1, pp. 51–79, 2016.
  * [13] 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, pp. 197–234, 2020.
  * [14] M. Korda and I. Mezić, “Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control,” Automatica, vol. 93, pp. 149–160, 2018.
  * [15] J. L. Proctor, S. L. Brunton, and J. N. Kutz, “Generalizing Koopman theory to allow for inputs and control,” SIAM Journal on Applied Dynamical Systems, vol. 17, no. 1, pp. 909–930, 2018.
  * [16] M. O. Williams, M. S. Hemati, S. T. Dawson, I. G. Kevrekidis, and C. W. Rowley, “Extending data-driven Koopman analysis to actuated systems,” IFAC-PapersOnLine, vol. 49, no. 18, pp. 704–709, 2016.
  * [17] R. W. Brockett, “Notes on the control of the Liouville equation. control of partial differential equations, lecture notes in mathematics,” 2012.
  * [18] K. Elamvazhuthi and P. Grover, “Optimal transport over nonlinear systems via infinitesimal generators on graphs.,” Journal of Computational Dynamics, vol. 5, 2018.
  * [19] Y. Chen, T. T. Georgiou, and M. Pavon, “Optimal transport in systems and control,” Annual Review of Control, Robotics, and Autonomous Systems, vol. 4, pp. 89–113, 2021.
  * [20] A. Halder and R. Bhattacharya, “Geodesic density tracking with applications to data driven modeling,” in 2014 American Control Conference, pp. 616–621, IEEE, 2014.
  * [21] P. Tallapragada and S. D. Ross, “A set oriented definition of finite-time Lyapunov exponents and coherent sets,” Comm. in Nonlinear Science and Numerical Simulation, vol. 18, no. 5, pp. 1106–1126, 2013.
  * [22] G. Froyland, N. Santitissadeekorn, and A. Monahan, “Transport in time-dependent dynamical systems: Finite-time coherent sets,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 20, no. 4, p. 043116, 2010\.
  * [23] 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), vol. 69, p. 25, 2005.
  * [24] 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), pp. 6107–6112, IEEE, 2017.
  * [25] 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, vol. 406, p. 132416, 2020.
  * [26] S. Peitz, S. E. Otto, and C. W. Rowley, “Data-driven model predictive control using interpolated koopman generators,” SIAM Journal on Applied Dynamical Systems, vol. 19, no. 3, pp. 2162–2193, 2020.
  * [27] E. Hofer and B. Tibken, “An iterative method for the finite-time bilinear-quadratic control problem,” Journal of optimization Theory and applications, vol. 57, pp. 411–427, 1988.
  * [28] 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, pp. 4906–4913, IEEE, 2012.
  * [29] S. Yakowitz and B. Rutherford, “Computational aspects of discrete-time optimal control,” Applied Mathematics and Computation, vol. 15, no. 1, pp. 29–45, 1984.
  * [30] J. Buzhardt and P. Tallapragada, “Optimal trajectory tracking for a magnetically driven microswimmer,” in 2020 American Control Conference (ACC), pp. 3211–3216, IEEE, 2020.