# Stochastic Thermodynamics of Thermophoresis

Mingnan Ding1<EMAIL_ADDRESS>Jun Wu1 Xiangjun Xing1,2,3<EMAIL_ADDRESS>1Wilczek Quantum Center, School of Physics and Astronomy, Shanghai Jiao Tong
University, Shanghai, 200240 China
2T.D. Lee Institute, Shanghai Jiao Tong University, Shanghai, 200240 China
3Shanghai Research Center for Quantum Sciences, Shanghai 201315, China

###### Abstract

We study stochastic thermodynamics of Brownian motion in a temperature
gradient. Using the covariant formulation of over-damped Ito-Langevin theory
and the associated Fokker-Planck theory, we construct excess entropy
production both at trajectory level and at ensemble level, and relate it to
the log-ratio between path probabilities of the forward process and the
backward process. We prove both a second-law like inequality and transient
fluctuation theorem for the excess entropy production. We also simulate a
Brownian particle inside a Lennard-Jones fluid and find that the fluctuation
theorem is valid both in the over-damped regime and in the under-damped
regime.

Introduction After several decades of hard work, two types of Fluctuation
theorems (FTs) [1, 4, 3, 5, 2, 6] have been firmly established. Firstly the
steady state FT (SSFT) characterizes the large deviation behaviors [4, 8, 7,
1] of entropy production (EP) rate in non-equilibrium steady states. Secondly,
various transient FTs [3, 2, 5] characterize the fluctuations of dissipated
work in non-equilibrium processes starting from and ending at equilibrium
states. These FTs constitute the backbone of the emerging field of stochastic
thermodynamics [5, 9, 10], which provides a unified view of non-equilibrium
fluctuations.

Many physical systems and most biological systems are however embedded in
dissipative backgrounds and undergo non-stationary processes. They therefore
do not fall into either of above mentioned two categories. Remarkably, a
series of seminal theoretical works [11, 12, 13, 14, 15, 16, 17] indicate that
in Markov processes with even variables, the total EP may be broken into two
parts, each satisfying a separate FT. For non-conservative Langevin systems,
it may be understood that these two parts of EP are due to respectively the
dissipative work of the conservative and non-conservative components of forces
[20]. This decomposition of EP is also intimately related to the Glansdorff-
Prigogine stability theory of non-equilibrium steady states [21, 18, 19, 20],
as well as the steady-state thermodynamics of Sasa and Tasaki [22, 20].
Following these earlier literatures, it is therefore natural to call them
respectively the excess EP and the house-keeping EP 111They are called
respectively the non-adiabatic EP and adiabatic EP by Esposito and Van den
Broeck [14, 15, 16].. We note that except for a few partial results [23, 24],
there has been no systematic verification of these theories, either
experimentally or numerically.

Brownian motion in temperature gradient constitutes one of the simplest
examples of systems embedded in dissipative backgrounds. Without any confining
potential, the particle exhibits directed motion, a phenomenon called the
Ludwig-Soret effect [25, 26, 27], or thermophoresis [29, 30, 28], which has
many important applications [31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42,
43, 44, 45]. One interesting recent variation is the micro-thermophoresis [31,
32], where micro-scale temperature gradient is produced by laser sources.
Whilst the thermodynamical theory of thermophoresis was understood long ago
[27, 28], its physical mechanism at microscopic level remains illusive [29].
Some recent studies on Brownian dynamics in temperature gradient [46, 47, 48,
49] yields disagreement between EP in under-damped theory and that in over-
damped theory, i.e., the so-called entropy anomaly [46, 47]. The Langevin
equations in these works however are constructed using the traditional ad hoc
Langevin approach, and do not capture the complex microscopic physics of
thermophoresis [29]. Furthermore, there has been no derivation or verification
of FT in the presence of temperature gradient. A systematic theory of
stochastic thermodynamics for thermophoresis is therefore still lacking.

In the present Letter, we shall develop such a theory, and make contact with
thermodynamics of thermophoresis [28, 29]. Firstly we assume that the dynamics
of fluid particles is much faster than that of the Brownian particles, so that
modeling of collisions with fluid particles as white noises can be justified.
Furthermore, to avoid the confusing issue of entropy anomaly, we take the
over-damped Langevin approach. Because the environment does not have a fixed
temperature, the conventional formulation of stochastic thermodynamics in
terms of heat and work is not directly applicable. We define and study the
excess EP both at the trajectory level and at the ensemble level, and prove
rigorously the Clausius inequality and the FT for excess EP. The FT is further
verified using molecular dynamic simulation of a colloid confined in liquid
with a temperature gradient. Remarkably, the FT remains valid even deep in the
under-damped regime. We explain this surprising result using coarse-graining
of stochastic thermodynamics.

Figure 1: Brownian particle in temperature gradient and external potential.
The dots show the fluid temperature determined numerically. The straight-line
is the linear fit of temperature, whereas the green shade gives the error bar.

Langevin dynamics Consider the over-damped Langevin equation:

$\displaystyle
dx^{i}+\left[\frac{T({\boldsymbol{x}})}{\gamma({\boldsymbol{x}})}\partial_{i}U({\boldsymbol{x}},\lambda)-\partial_{i}\frac{T({\boldsymbol{x}})}{\gamma({\boldsymbol{x}})}\right]dt=\sqrt{\frac{2T({\boldsymbol{x}})}{\gamma({\boldsymbol{x}})}}dW^{i},\quad$
(1)

which, as illustrated in Fig. 1, describes a Brownian particle in a
temperature gradient. Here $T({\boldsymbol{x}}),\gamma({\boldsymbol{x}})$ are
respectively temperature and friction coefficient, and
$D({\boldsymbol{x}})=T({\boldsymbol{x}})/\gamma({\boldsymbol{x}})$ can be
understood as the diffusion constant, all of which are position dependent.
$dW^{i}(t)$ is the Wiener noise satisfying $dW^{i}(t)dW^{k}(t)=\delta_{ik}dt$,
whereas the product in the r.h.s. of Eq. (1) is interpreted in Ito’s sense.
Equation (1) is mathematically equivalent to the Fokker-Planck equation:

$\displaystyle\partial_{t}\,p=\partial_{i}\frac{T({\boldsymbol{x}})}{\gamma({\boldsymbol{x}})}\left(\partial_{i}+\left(\partial_{i}U\right)\right)p\equiv\mathcal{L}_{\rm
FP}\,p,$ (2)

where $p({\boldsymbol{x}},t)$ is the pdf of $x$, and $\mathcal{L}_{\rm FP}$ is
the Fokker-Planck operator. Equation (2) admits a steady state solution
$p_{\rm ss}({\boldsymbol{x}})=e^{-U({\boldsymbol{x}},\lambda)}$, where
$U({\boldsymbol{x}},\lambda)$ is the generalized potential. The relation
between $U({\boldsymbol{x}},\lambda)$ and the external potential
$V({\boldsymbol{x}},\lambda)$ encodes the rich microscopic physics of
thermophoresis, which may be measured in experiments or computed in
simulations. $\lambda$ is the control parameter, which may represent the
location or strength of potential.

Note that ${\boldsymbol{x}}$ and $\lambda$ are both even under time-reversal.
Equations (1) and (2) therefore are in the covariant form studied in Ref. [50,
51], with a diagonal but position dependent kinetic matrix
$L^{ij}({\boldsymbol{x}})=(T({\boldsymbol{x}})/\gamma({\boldsymbol{x}}))\delta_{ij}$.
They also satisfy conditions of detailed balance (Eqs. (2.31) of Ref. [51]),
which guarantee that two-time pdfs of position satisfy:
$p_{\lambda}({\boldsymbol{x}}_{2},t_{2};{\boldsymbol{x}}_{1},t_{1})=p_{\lambda}({\boldsymbol{x}}_{1},t_{1};{\boldsymbol{x}}_{2},t_{2})$.
The probability current associated with Eq. (2) is

$\displaystyle
j^{i}({\boldsymbol{x}},t,\lambda)=-\frac{T({\boldsymbol{x}})}{\gamma({\boldsymbol{x}})}(\partial_{i}+(\partial_{i}U))p({\boldsymbol{x}},t),$
(3)

which vanishes identically in the steady state.

Excess Clausius inequality We define the excess change of environmental
entropy at the trajectory level:

$\displaystyle d\mathcal{S}_{\rm env}^{\rm
ex}=-d_{{\boldsymbol{x}}}U\equiv-\partial_{i}U\circ dx^{i},$ (4)

where $d_{{\boldsymbol{x}}}U\equiv
U({\boldsymbol{x}}+d{\boldsymbol{x}})-U({\boldsymbol{x}})$ denotes
differential of $U$ due to change of $x$, and $\circ$ is the product in
Stratonovich’s sense. Note that repeated indices are summed over. Averaging
this over both the noises and the pdf 222This is the meaning of
$\langle\\!\langle\,\cdot\rangle\\!\rangle$ below., we obtained the
corresponding quantity at the ensemble level:

$\displaystyle dS_{\rm env}^{\rm ex}=\langle\\!\langle d\mathcal{S}_{\rm
env}^{\rm
ex}\rangle\\!\rangle=dt\int_{{\boldsymbol{x}}}(\partial_{i}U)\frac{T({\boldsymbol{x}})}{\gamma({\boldsymbol{x}})}(\partial_{i}p+p\,\partial_{i}U).$
(5)

For details see Supplementary Information (SI) Sec. I. The differential of
system entropy is

$\displaystyle dS_{\rm sys}=-d\int_{{\boldsymbol{x}}}p\,\log
p=-dt\int_{{\boldsymbol{x}}}\log p\,\mathcal{L}_{\rm FP}p.$ (6)

The excess entropy production 333This is called the non-adiabatic entropy
production by Esposito and van der Broeck. Here we follow the terminology of
Glansdorff and Prigogine. at ensemble level is defined as the sum of Eqs. (5)
and (6), which turns out to be non-negative and vanishes only in the steady
state:

$\displaystyle dS^{\rm ex}$ $\displaystyle=$ $\displaystyle dS_{\rm env}^{\rm
ex}+dS_{\rm sys}$ (7) $\displaystyle=$ $\displaystyle
dt\int_{\boldsymbol{x}}\frac{T({\boldsymbol{x}})}{p\gamma({\boldsymbol{x}})}\left(e^{-U}\partial_{i}p\,e^{U}\right)^{2}\geq
0.$

This is called the second law of steady states by Hatano and Sasa [22, 11].

We may also define a dimensionless functional

$\Phi=\int_{{\boldsymbol{x}}}p({\boldsymbol{x}})[U({\boldsymbol{x}},\lambda)+\log
p({\boldsymbol{x}})],$ (8)

which, for fixed $\lambda$, is minimized at the steady state. This is the
analogue of minimum free energy principle in classical thermodynamics.

It is important to emphasize that at the level of our description, the
environment consists of the ambient fluid as well as the thermostats who are
used to generate the temperature gradient. The environment therefore does not
have a fixed temperature. As a consequence, the even though Eqs. (1) and (2)
have detailed balance, the steady state $p_{\rm
ss}({\boldsymbol{x}})=e^{-U({\boldsymbol{x}},\lambda)}$ is of truly non-
equilibrium nature. Additionally, $d\mathcal{S}^{\rm ex}_{\rm env}$ as defined
by Eq. (7) is not the entropy change of the environment, because it does not
take into account the dissipation due to heat transport in the ambient fluid,
which is non-vanishing even if the Brownian particle stand still. For the
special case of vanishing temperature gradient, however, $T$ is constant, and
we may rewrite the generalized potential as $\beta(H-F)$, where $H$ is the
potential of mean force and $F$ is the equilibrium free energy. Equation (4)
then reduces to $-\beta d_{{\boldsymbol{x}}}H=-\beta\,d\bar{}\hskip
1.00006pt\mathcal{Q}$, where $d\mathcal{Q}$ is the heat absorbed by the
particle [52]. Equation (7) then reduces to $dS_{\rm sys}-\beta d\bar{}\hskip
1.00006ptQ\geq 0$, which is the usual Clausius inequality for systems in
contact with a single heat bath [51]. Also in this case, Eq. (8) reduces to
$\beta F[p]$, where $F[p]$ is the non-equilibrium free energy.

Excess fluctuation theorem Consider a forward process with dynamic protocol
$\lambda(t)$ defined in the time interval $0\leq t\leq\tau$, and initial
parameter $\lambda(0)=\lambda_{i}$. The backward process is defined by the
protocol $\tilde{\lambda}(t)=\lambda(\tau-t)$, and initial parameter
$\tilde{\lambda}(0)=\lambda(\tau)=\lambda_{f}$. In both the forward and
backward processes, the temperature gradient remains fixed, whereas the system
starts from steady states:

$\displaystyle
p_{F}({\boldsymbol{x}},t=0)=e^{-U({\boldsymbol{x}},\lambda_{i})},$ (9a)
$\displaystyle
p_{B}({\boldsymbol{x}},t=0)=e^{-U({\boldsymbol{x}},\lambda_{f})}.$ (9b)

Since Eq. (1) is in the covariant form, we may invoke Eq. (A27) of Ref. [51]
and obtain the following relation:

$\displaystyle\log\frac{P_{F}({\boldsymbol{x}}_{1}|{\boldsymbol{x}};dt)}{P_{B}({\boldsymbol{x}}|{\boldsymbol{x}}_{1};dt)}=-d_{{\boldsymbol{x}}}U=d\mathcal{S}_{\rm
env}^{\rm ex},$ (10)

where
$p_{F}({\boldsymbol{y}}|{\boldsymbol{x}};dt),p_{B}({\boldsymbol{y}}|{\boldsymbol{x}};dt)$
are respectively the short-time transition probabilities for the forward and
backward processes. Using Eqs. (9) and (10), we can prove the following for a
dynamic path ${\boldsymbol{\gamma}}={\boldsymbol{x}}(t)$ in the forward
process and the reversed path
$\tilde{\boldsymbol{\gamma}}={\boldsymbol{x}}(\tau-t)$ in the backward process
444It is important to distinguish $\gamma$ the friction coefficient and
${\boldsymbol{\gamma}}$ (boldface) the dynamic trajectory.:

$\displaystyle\log\frac{P_{F}[{\boldsymbol{\gamma}}]}{P_{B}[\tilde{\boldsymbol{\gamma}}]}$
$\displaystyle=$ $\displaystyle
U({\boldsymbol{x}}_{f},\lambda_{f})-U({\boldsymbol{x}}_{i},\lambda_{i})-\medint\int_{\boldsymbol{\gamma}}d_{{\boldsymbol{x}}}U$
$\displaystyle=$ $\displaystyle\Sigma^{\rm ex}_{\rm
F}[{\boldsymbol{\gamma}}]=-\Sigma^{\rm ex}_{\rm
B}[\tilde{\boldsymbol{\gamma}}]=\medint\int_{\boldsymbol{\gamma}}d_{\lambda}U,$

where
${\boldsymbol{x}}_{i}={\boldsymbol{x}}(0),{\boldsymbol{x}}_{f}={\boldsymbol{x}}(\tau)$
are respectively the initial and final state of the forward trajectory.
Detailed proof of Eq. (Stochastic Thermodynamics of Thermophoresis) is
supplied in Sec. II of SI. The r.h.s. of Eq. (Stochastic Thermodynamics of
Thermophoresis) is the sum of change of stochastic entropy of the system and
the excess entropy change of the environment along the trajectory, and may be
understood as the excess EP along the path ${\boldsymbol{\gamma}}$ in the
forward process, or negative the excess EP along $\tilde{\boldsymbol{\gamma}}$
in the backward process.

We can now follow the standard procedure of stochastic thermodynamics and
define the pdfs of excess EP of the forward and backward processes:

$\displaystyle p_{F}(\Sigma^{\rm ex})\equiv\int
D\Gamma\,p_{F}(\Gamma)\,\delta\left(\Sigma^{\rm
ex}-\log\frac{p_{F}[\Gamma]}{p_{B}[\tilde{\Gamma}]}\right),\quad\quad$ (12a)
$\displaystyle p_{B}(\Sigma^{\rm ex})\equiv\int
D\tilde{\Gamma}\,p_{B}(\tilde{\Gamma})\,\delta\left(\Sigma^{\rm
ex}+\log\frac{p_{F}[\Gamma]}{p_{B}[\tilde{\Gamma}]}\right),$ (12b)

and prove the following transient fluctuation theorem:

$\displaystyle\log\frac{p_{F}(\Sigma^{\rm ex})}{p_{B}(-\Sigma^{\rm
ex})}={\Sigma^{\rm ex}},$ (13)

which shall be called the excess FT. For vanishing temperature gradient, Eq.
(Stochastic Thermodynamics of Thermophoresis) reduces to
$\beta(\mathcal{W}-\Delta F)$, and Eq. (13) reduces to the usual Crooks FT.

Numerical simulation To test our theory, we simulate a colloid immersed in a
2d fluid, which consists of 2809 particles with unit mass confined in a
$200\times 200$ square and interacting via truncated Lennard-Jones (LJ)
potential. As schematized in Fig. 1, a cold bath with dimensionless
temperature $T_{h}=5.8$ and a hot heat bath with $T_{c}=4.2$ are introduced
near the left and right boundaries respectively. Further, periodic boundary
condition is imposed on the vertical direction, whereas reflection boundary
condition is imposed on horizontal direction, so that there is no macroscopic
transport of fluid particles in the steady state. Details of simulation
methods are explained in Sec. III of SI.

We analyze the distribution of velocity of the fluid particles for given
horizontal coordinate $x$. As demonstrated in detail in Sec. IV of SI, the pdf
of $v_{y}$ is Gaussian, with a variance linear in $x$ in the bulk, whereas the
pdf of $v_{x}$ exhibits a small yet systematic deviation from Gaussian. We use
several distinct methods to compute the local temperature as a function of
$x$, all of which give a temperature profile:
$T({\boldsymbol{x}})=-0.0083\,x+5.005$, as shown in Fig. 3 of SI. See also
Fig. 1 for an illustration.

A Brownian particle with triple the linear size and ten times the mass of
fluid particles is introduced into the fluid. It interacts with the fluid via
the truncated LJ potential, as well as an anisotropic harmonic potential:

$\displaystyle
V({\boldsymbol{x}})=\frac{1}{2}\kappa_{x}(x-x_{0})^{2}+\frac{1}{2}\kappa_{y}y^{2}.$
(14)

We analyze the velocity pdf of Brownian particle in the steady state. As shown
in detail in Sec. VA of SI, the pdf of $v_{x}$ is independent of $\kappa_{x}$,
and exhibits small deviation from Gaussian. The variances of $v_{x},v_{y}$ are
however consistent with those of fluid particles, see Fig. 4 (b) and (c) of
SI.

The steady state pdf of the position of Brownian particle is used to compute
the generalized potential via $U({\boldsymbol{x}},\lambda)=-\log p_{\rm
SS}({\boldsymbol{x}},\lambda)$. In Sec. VB of SI, we show that for $V\neq 0$,
the following approximation is good:

$\displaystyle
U({\boldsymbol{x}},\lambda)=\beta({\boldsymbol{x}})V({\boldsymbol{x}},\lambda)+S_{T}\,{\boldsymbol{x}}\cdot\nabla
T+b(\lambda).$ (15)

Whilst the constant $b(\lambda)$ only serves to normalize the steady state
pdf, the parameter $S_{T}$, which is positively valued in our simulation,
turns out to be the Soret coefficient [29]. To see this, we recall in the
thermodynamic theory of thermophoresis, where there is no external potential,
the probability current of the Brownian particle can be written as [29]

$j=-D\,\nabla p-D\,S_{T}\,T^{\prime}\,p.$ (16)

This must agree with Eq. (3), the current given by FP theory. Indeed,
substituting Eq. (15) with $V=0$ into Eq. (3) and noticing $T/\gamma=D$, we
obtain Eq. (16). Note that in principle the Soret coefficient $S_{T}$, same as
the diffusion coefficient, depends on local temperature. In the region of
simulation, the overall variation of temperature is less than 8%, hence
$S_{T}$ can be treated approximately as a constant. For the special case
$V=0$, Eq. (15) is unbounded, and hence the steady state is not normalizable.
As shown in Sec. VC of SI, the Brownian particle drifts towards lower
temperature with approximately constant speed. Finally as detailed in Sec. VD
of SI, we also compute the local friction constant $\gamma({\boldsymbol{x}})$
using two different methods, and obtain consistent result. Note however,
$\gamma({\boldsymbol{x}})$ does not enter Eq. (4), hence numerical value of
$\gamma({\boldsymbol{x}})$ is not needed for verification of excess FT.

Figure 2: (a), (b), (c): Distributions of excess EP. $F,B$ means forward and
backward respectively. In (a),
$x_{0}=-10+20t/\tau,\kappa_{x}=0.01,\kappa_{y}=0.01$; In (b),
$x_{0}=0,\kappa_{x}=0.01+0.1t/\tau,\kappa_{y}=0.1$; in (c),
$x_{0}=-10+20t/\tau,\kappa_{x}=0.01+0.1t/\tau,\kappa_{y}=0.1$. (d):
verification of FT.

Verification of FT We simulate three types of forward process as well as
corresponding backward process, and calculate the excess EP defined by Eq.
(Stochastic Thermodynamics of Thermophoresis), where the generalized potential
is calculated using the fitting formula (15). Throughout the simulation we fix
$\kappa_{y}=0.1$, so that the Brownian particle is localized in $y$ direction.
In the first protocol, we change the equilibrium position $x_{0}$ with
$\kappa_{x}$ fixed. In the second protocol, we fix $x_{0}$ but tune
$\kappa_{x}$. In the third protocol we tune $\kappa_{x},x_{0}$ simultaneously.
The pdfs of $\Sigma^{\rm ex}$ are shown in Fig. 2(a), (b), (c), respectively,
whereas in Fig. 2(d), the log ratio $\log p_{F}(\Sigma^{\rm
ex})/p_{B}(-\Sigma^{\rm ex})$ is plotted against $\Sigma^{\rm ex}$. As shown
there, all data collapse to the straight-line as predicted by the excess FT
Eq. (13).

Figure 3:  Verification of excess FT in under-damped regime. In the first two
protocols, $x_{0}=-10+20t/\tau,\kappa_{x}=0.1,\kappa_{y}=0.1$; in the
remaining two protocols,
$x_{0}=-10+20t/\tau,\kappa_{x}=0.1+0.1t/\tau,\kappa_{y}=0.1$.

The velocity auto-correlation function (VACF) of the Brownian particle are
shown in Fig. 5 of SI. It is seen that the dynamics is over-damped for
$\kappa\leq 0.05$ and under-damped otherwise. Hence in the second and third
protocols, as we tune up the spring constant $\kappa_{x}$ from $0.01$ to
$0.11$, the over-damped Langevin equation (1) is not valid in the later part
of the forward process. Nonetheless, the excess FT seems hold to high
precision, as can be seen from Fig. 2(d). As shown in Fig. 3, we also simulate
two other protocols where the system is deep in the over-damped regime in the
entire process, and find that the excess FT (13) remains valid as well.

The validity of the excess FT in the under-damped regime may be understood
using the stochastic thermodynamic theory developed in Ref. [20], and the
coarse-graining method of stochastic thermodynamics developed in Ref. [53] and
[54]. In principle, we may construct a fine-grained stochastic thermodynamics
which includes both the fluid particles and the Brownian particle, with
temperature gradient taken into account as some kind of non-conservative
forces. For this system, local detailed balance holds, and the theory
developed in Ref. [20] is applicable. Furthermore, all variables except the
position of Brownian particle are decoupled from the control parameter. We can
the apply the coarse-grained theory developed in Ref. [54], integrate our all
variables of fluid particles, as well as the velocity of Brownian particle.
This yields the over-damped theory which was studied in the present work. As
proved in Ref. [54], the excess FT remains valid in the coarse-grained theory.

The authors acknowledge support from Chenxing Post Doctoral Incentive Project,
NSFC grant 11674217, and Shanghai Municipal Science and Technology Major
Project via grant 2019SHZDZX01.

## References

  * [1] Gallavotti, Giovanni, and Ezechiel Godert David Cohen. “Dynamical ensembles in nonequilibrium statistical mechanics.” Physical review letters 74.14 (1995): 2694.
  * [2] Jarzynski, C. Equalities and inequalities: Irreversibility and the second law of thermodynamics at the nanoscale. Annu. Rev. Condens. Matter Phys., 2(1), 329-351. (2011)
  * [3] Crooks, G. E. Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences. Physical Review E, 60(3), 2721. (1999)
  * [4] Searles, D. J. and Evans, D. J. Fluctuation theorem for stochastic systems. Physical Review E, 60(1), 159 (1999)
  * [5] Seifert, U. Stochastic thermodynamics, fluctuation theorems and molecular machines. Reports on progress in physics, 75(12), 126001 (2012).
  * [6] Ciliberto, Sergio. Experiments in stochastic thermodynamics: Short history and perspectives. Physical Review X 7.2 (2017): 021051.
  * [7] Lebowitz, Joel L., and Herbert Spohn. A Gallavotti–Cohen-type symmetry in the large deviation functional for stochastic dynamics. Journal of Statistical Physics 95 (1999): 333-365.
  * [8] Van Zon, R., and Cohen, E. G. D. Stationary and transient work-fluctuation theorems for a dragged Brownian particle. Physical Review E, 67(4), 046102. (2003)
  * [9] Peliti, Luca, and Simone Pigolotti. Stochastic Thermodynamics: An Introduction. Princeton University Press, 2021.
  * [10] Van den Broeck, Christian. “Stochastic thermodynamics: A brief introduction.” Phys. Complex Colloids 184 (2013): 155-193.
  * [11] Hatano, T. and Sasa, S. I. Steady-state thermodynamics of Langevin systems. Physical review letters, 86(16), 3463. (2001)
  * [12] Speck, Thomas, and Udo Seifert.“Integral fluctuation theorem for the housekeeping heat.” Journal of Physics A: Mathematical and General 38.34 (2005): L581.
  * [13] Chernyak, Vladimir Y., Michael Chertkov, and Christopher Jarzynski. “Path-integral analysis of fluctuation theorems for general Langevin processes.” Journal of Statistical Mechanics: Theory and Experiment 2006.08 P08001 (2006).
  * [14] Esposito, M. and Van den Broeck, C. Three detailed fluctuation theorems. Physical review letters, 104(9), 090601. (2010)
  * [15] Esposito, M. and Van den Broeck, C. Three faces of the second law. I. Master equation formulation. Physical Review E, 82(1), 011143. (2010)
  * [16] Van den Broeck, C. and Esposito, M. Three faces of the second law. II. Fokker-Planck formulation. Physical Review E, 82(1), 011144. (2010)
  * [17] Ge, Hao, and Hong Qian. “Physical origins of entropy production, free energy dissipation, and their mathematical representations.” Physical Review E 81.5 (2010): 051133.
  * [18] Qian, Hong. “Mesoscopic nonequilibrium thermodynamics of single macromolecules and dynamic entropy-energy compensation.” Physical Review E 65.1 (2001): 016102.
  * [19] Qian, Hong. “Entropy production and excess entropy in a nonequilibrium steady-state of single macromolecules.” Physical Review E 65.2 (2002): 021111.
  * [20] Mingnan Ding, Fei Liu and Xiangjun Xing. “A Unified Theory for Thermodynamics and Stochastic Thermodynamics of Nonlinear Langevin Systems Driven by Non-conservative Forces” Physical Review Research 4.4 (2022): 043125.
  * [21] Glansdorff, Paul, and Ilya Prigogine. Thermodynamic theory of structure, stability and fluctuations. J. Willey & Sons, 1971.
  * [22] Sasa, Shin-ichi, and Hal Tasaki. “Steady state thermodynamics.” Journal of statistical physics 125.1 (2006): 125-224.
  * [23] Trepagnier, E. H., Jarzynski, C., Ritort, F., Crooks, G. E., Bustamante, C. J. and Liphardt, J. Experimental test of Hatano and Sasa’s nonequilibrium steady-state equality. Proceedings of the National Academy of Sciences, 101(42), 15038-15041.(2004)
  * [24] Thalheim, Tobias, et al. “Jarzynski Equality for Soret Equilibria–Virtues of Virtual Potentials.” arXiv preprint arXiv:2009.00950 (2020).
  * [25] Ludwig C 1859 Sitz. Ber. Akad. Wiss. Wien 20 539.
  * [26] Soret C 1879 Arch. Sci. Phys. Nat. 2 48.
  * [27] de Groot S R 1945 L’Effet Soret, Diffusion Thermique Dans Les Phasescondensées (Amsterdam: North-Holland).
  * [28] De Groot, Sybren Ruurds, and Peter Mazur. Non-equilibrium thermodynamics. Courier Corporation, 2013.
  * [29] Piazza, Roberto, and Alberto Parola. “Thermophoresis in colloidal suspensions.” Journal of Physics: Condensed Matter 20.15 (2008): 153102.
  * [30] Piazza, Roberto. “Thermophoresis: moving particles with thermal gradients.” Soft Matter 4.9 (2008): 1740-1744.
  * [31] Jerabek-Willemsen, Moran, et al. “MicroScale Thermophoresis: Interaction analysis and beyond.” Journal of Molecular Structure 1077 (2014): 101-113.
  * [32] Jerabek-Willemsen, Moran, et al. “Molecular interaction studies using microscale thermophoresis.” Assay and drug development technologies 9.4 (2011): 342-353.
  * [33] Thamdrup, Lasse H., Niels B. Larsen, and Anders Kristensen. “Light-induced local heating for thermophoretic manipulation of DNA in polymer micro-and nanochannels.” Nano letters 10.3 (2010): 826-832.
  * [34] Baaske, Philipp, et al. “Optical thermophoresis for quantifying the buffer dependence of aptamer binding.” Angewandte Chemie International Edition 49.12 (2010): 2238-2241.
  * [35] Wen, Y., Liu, Q. and Liu, Y. “Temperature gradient-driven motion and assembly of two-dimensional (2D) materials on the liquid surface: a theoretical framework and molecular dynamics simulation. ” Physical Chemistry Chemical Physics, 22(41), 24097-24108.(2020)
  * [36] Becton, M. and Wang, X. Thermal gradients on graphene to drive nanoflake motion. Journal of chemical theory and computation, 10(2), 722-730. (2014)
  * [37] Dehbi, A. A stochastic Langevin model of turbulent particle dispersion in the presence of thermophoresis. International Journal of Multiphase Flow, 35(3), 219-226.(2009)
  * [38] Pakravan, H. A. and Yaghoubi, M. Combined thermophoresis, Brownian motion and Dufour effects on natural convection of nanofluids. International Journal of Thermal Sciences, 50(3), 394-402.(2011)
  * [39] Helden, L., Eichhorn, R. and Bechinger, C. Direct measurement of thermophoretic forces. Soft matter, 11(12), 2379-2386. (2015)
  * [40] Braibanti, M., Vigolo, D. and Piazza, R. Does thermophoretic mobility depend on particle size? Physical review letters, 100(10), 108303.(2008)
  * [41] Rurali, R. and Hernandez, E. R. Thermally induced directed motion of fullerene clusters encapsulated in carbon nanotubes. Chemical Physics Letters, 497(1-3), 62-65. (2010)
  * [42] Schoen, P. A., Walther, J. H., Poulikakos, D. and Koumoutsakos, P. Phonon assisted thermophoretic motion of gold nanoparticles inside carbon nanotubes. Applied Physics Letters, 90(25), 253116. (2007)
  * [43] Duhr, S. and Braun, D. Thermophoretic depletion follows Boltzmann distribution. Physical review letters, 96(16), 168301. (2006)
  * [44] Duhr, S. and Braun, D. Why molecules move along a temperature gradient. Proceedings of the National Academy of Sciences, 103(52), 19678-19682. (2006)
  * [45] Kroy, Klaus, Dipanjan Chakraborty, and Frank Cichos. “Hot microswimmers.” The European Physical Journal Special Topics 225 (2016): 2207-2225.
  * [46] Celani, A., Bo, S., Eichhorn, R. and Aurell, E. Anomalous thermodynamics at the microscale. Physical review letters, 109(26), 260603. (2012)
  * [47] Marino, Raffaele, Ralf Eichhorn, and Erik Aurell. “Entropy production of a Brownian ellipsoid in the over-damped limit.” Physical Review E 93.1 (2016): 012132.
  * [48] Sancho, J. M. Brownian colloids in underdamped and overdamped regimes with nonhomogeneous temperature. Physical Review E, 92(6), 062110. (2015)
  * [49] Polettini, M. Diffusion in nonuniform temperature and its geometric analog. Physical Review E, 87(3), 032126. (2013)
  * [50] Ding, M., Tu, Z. and Xing, X. Covariant formulation of nonlinear Langevin theory with multiplicative Gaussian white noises. Physical Review Research, 2(3), 033381. (2020).
  * [51] Ding, M. and Xing, X. (2022). Covariant nonequilibrium thermodynamics from Ito-Langevin dynamics. Physical Review Research, 4(3), 033247.
  * [52] Ding, M., Tu, Z. and Xing, X. Strong coupling thermodynamics and stochastic thermodynamics from the unifying perspective of time-scale separation. Physical Review Research, 4(1), 013015. (2022).
  * [53] Mingnan Ding, Xiangjun Xing, “Multi-scale Projection Operator Method and Coarse-graining of Fokker-Planck Theory”, to appear in Physical Review Research.
  * [54] Mingnan Ding, Xiangjun Xing, “Thermodynamics and Stochastic Thermodynamics of Nonlinear Langevin Systems Driven by Non-conservative Forces”, submitted to Physical Review E.


# Determining Follower Aircraft’s Optimal Trajectory in Relation to a Dynamic
Formation Ring

Carl A. Gotwald & Michael D. Zollars
Department of Aeronautics and Astronautics
Air Force Institute of Technology
Wright-Patterson AFB  OH 45433
<EMAIL_ADDRESS><EMAIL_ADDRESS>Isaac E. Weintraub
Controls Science System Center
Wright-Patterson AFB  OH 45433
<EMAIL_ADDRESS>U.S. Government work not protected by U.S.
copyright

###### Abstract

The specific objective of this paper is to develop a tool that calculates the
optimal trajectory of the follower aircraft as it completes a formation
rejoin, and then maintains the formation position, defined as a ring of
points, until a fixed final time. The tool is designed to produce optimal
trajectories for a variety of initial conditions and leader trajectories.
Triple integrator dynamics are used to model the follower aircraft in three
dimensions. Control is applied directly to the rate of acceleration. Both the
follower’s and leader’s velocities and accelerations are bounded, as dictated
by the aircraft’s performance envelope. Lastly, a path constraint is used to
ensure the follower avoids the leader’s jet wash region. This optimal control
problem is solved through numerical analysis using the direct orthogonal
collocation solver GPOPS-II. Two leader trajectories are investigated,
including a descending spiral and continuous vertical loops. Additionally, a
study of the effect of various initial guesses is performed. All trajectories
displayed a direct capture of the formation position, however changes in
solver initial conditions demonstrate various behaviors in how the follower
maintains the formation position. The developed tool has proven adequate to
support future research in crafting real-time controllers capable of
determining near-optimal trajectories.

###### Contents

  1. 1 Introduction
  2. 2 Problem Formulation
    1. 2.1 Assumptions
    2. 2.2 Reference Frames
    3. 2.3 Mathematical Formulation
    4. 2.4 Discrete Hamiltonian Derivation
    5. 2.5 Test Cases
  3. 3 Solution Methodology
    1. 3.1 Direct Collocation Solver
    2. 3.2 NLP Solver
    3. 3.3 Initial Guess Generation
  4. 4 Results
    1. 4.1 Spiral Trajectory
    2. 4.2 Loop Trajectory
    3. 4.3 Initial Guess Comparison
    4. 4.4 Analysis of Optimally
  5. 5 Conclusions

##  1 Introduction

Autonomous formation control is a growing area of interest for future Air
Force operations. Autonomous aircraft are increasingly used across the battle
space to avoid putting personnel at risk. Thus, the need for autonomous
aircraft to handle formation tasks normally requiring a pilot has grown. One
model for this interaction consists of a leader aircraft who is accompanied by
an autonomous follower aircraft. The ability for the follower aircraft to
autonomously rejoin and maintain a designated formation position would free
the leader aircraft to focus on other duties, leading to increased mission
effectiveness. The leader aircraft often has primary responsibility for
interactions outside of the formation, so having an autonomous wingman which
can execute without additional oversight from the lead aircraft in an
efficient manner would decrease the leader workload. The scenario investigated
in this paper focused on the response of a wingman rejoining from an arbitrary
starting location to a defined formation position and maintaining that
position until a fixed final time. This same behavior could readily be applied
to the scenario of an autonomous aircraft attempting to rapidly attain an
attack position on an enemy aircraft, along with many other combat and
training applications.

Figure 1: Fighting Wing Position [1]

Extensive research in maintaining a desired formation position, defined as a
single location relative to a leader aircraft, has been conducted in articles
such as [2], [3], and [4]. The author of [5] presented a model predictive
controller which added the ability to arrive at a target waypoint with a
desired velocity vector orientation. To extend this concept beyond single-
point formation positions, further research was conducted which crafted a
control architecture allowing a follower aircraft to rejoin to a ring of
points defined relative to a leader aircraft [6]. This ring of formation
points was inspired by the traditional formation position ”fighting wing” as
seen in Figure 1 referenced from [1]. Future research aims to dictate a
trajectory in real-time for which an autonomous wingman would attain and
maintain a formation position defined in reference to a maneuvering leader
while minimizing a desired cost function. In order to attain an approximate
minimum in real-time may require the use of different control laws depending
on the initial states of the leader and follower aircraft. Thus, a state space
of boundaries could be defined to dictate when certain control laws would be
optimal. Further, future research aims to create a strategy which will
optimally select from these predefined control laws and implement them in
real-time. Since it is unlikely a fully optimal solution can be attained in
real-time for this problem, the control strategy will aim to approximate an
optimal trajectory. Therefore, the ability to calculate the degree of the
error of this approximation is needed. This paper specifically aimed to create
a validation tool that will determine the optimal path for future comparison
to the approximated trajectory generated by the real-time control strategy.

One specific scenario was considered. The follower aircraft being tasked to
rejoin to a defined formation ring in minimum time, and then maintaining that
formation position with minimum deviation over a fixed time duration. During
this task, the leader aircraft maneuvered in all three axes, similar to the
behavior dictated by typical formation roles. The leader’s maneuvers were
determined a priori for these simulations, and were treated as a given set of
parameters for each test case. To simulate the follower aircraft’s flight
envelope, constraints were applied to the magnitude of both the follower’s
velocity and acceleration. The leader aircraft maintained its own constraints
along the predetermined path. Two cases were analyzed and evaluated to assess
the tools effectiveness at determining the optimal path. The analysis of this
problem is organized in the following sections: Section 2 details the problem
formulation, Section 3 outlines the solution methodology, Section 4 includes
discussion of the results, and Section 5 includes the conclusions and future
recommendations.

##  2 Problem Formulation

###  2.1 Assumptions

To aid in the initial creation of the validation tool, a number of simplifying
assumptions were used. The follower aircraft was modeled with triple
integrator dynamics, with control applied as a three-dimensional jerk vector.
This model assumed perfect knowledge of the aerodynamic forces and thrust
interactions. Full knowledge of the leader’s path was provided to the solver
for both cases, including three-dimensional position, velocity, acceleration,
and jerk vectors defined in the inertial reference frame. No exogenous inputs
were considered within the tool. Therefore, deterministic knowledge of all the
follower states was assumed, which included three-dimensional position,
velocity, and acceleration.

Both the leader and follower aircraft were modeled as point masses about their
center of gravity. For any images created throughout this report, the velocity
vector was used as a reference to determine the pitch and the yaw of the
aircraft. However, without modeling the full Euler angles of the aircraft the
ability to determine the roll angle was limited. To aid in graphics
development, roll was visually estimated.

In formulating the optimal control problem, the existence of an optimal
solution was assumed. A scenario without a solution could be easily crafted,
such as the simple case where a trailing follower cannot accelerate fast
enough to catch a higher performance leader aircraft in the two-minute
simulation time. This assumption was satisfied within this study by choosing
appropriate initial conditions and using matching aircraft performance
envelopes. This assumption could easily be relaxed while still providing
optimal solutions in future work. Additionally, the test cases presented
leader trajectories such that the opportunity existed for an optimal rejoin
and formation position maintenance.

Lastly, the aircraft performance envelope was modeled with a minimum velocity
and acceleration to prevent aerodynamic stall, a maximum velocity to limit
maximum dynamic pressure forces, and a maximum acceleration to limit maximum
structural forces.

###  2.2 Reference Frames

This problem required two primary reference frames: the inertial reference
frame and the leader reference frame. Both are shown in Fig 2, reproduced from
[6].

Figure 2: The formation ring is defined in the leader frame as shown by the
red circle, in trail of the leader. The follower strives to rejoin by reaching
the ring and then sustaining that position thereafter. [6]

The inertial reference frame was defined as a traditional north, east, down
right-handed coordinated system with an origin fixed to the starting north and
east location of the follower aircraft at zero feet mean sea level (MSL). The
body fixed reference frame of interest was fixed to the center of gravity of
the leader aircraft. This leader reference frame was defined with a
longitudinal x axis for which positive extended out the nose of the aircraft,
and a lateral y axis with positive defined as pointing out the right wing. The
third axis of the leader frame was defined by the cross product of the x and y
axis, using the right-hand rule to determine the sign. The rotation matrix
between the two reference frames was defined by the leader’s flight path
angle, $\gamma$, and the course angle $\chi$. The applicable direction cosine
matrix to rotate a vector from the leader frame to the inertial frame was
given by the following equation referenced in [7].

$R_{L}^{i}=\begin{bmatrix}\cos{\chi}&-\sin{\chi}&0\\\
\sin{\chi}&\cos{\chi}&0\\\
0&0&1\end{bmatrix}\begin{bmatrix}\cos{\gamma}&0&\sin{\gamma}\\\ 0&1&0\\\
-\sin{\gamma}&0&\cos{\gamma}\end{bmatrix}$ (1)

Additionally, the rotation matrix in (1) satisfies the relations
$(R_{i}^{L})^{-1}=[R_{L}^{i}]^{\top}=R_{L}^{i}$ and $\text{det}(R_{i}^{L})=1$.
$R_{i}^{L}$ defines the rotation from the inertial frame to the leader frame.
A superscript L on any vector denotes a vector defined in the leader reference
frame, while all other vectors are defined in the inertial frame. The
subscript $l$ and $f$ represent the leader and follower vectors respectively.
Lastly, the vector $p_{d}$ represent the vector LF in the inertial frame as
shown in Figure 2. This formulation was chosen to match the research in [6].
The follower and leader states and control jerk vector are defined in relation
to the inertial coordinate system as follows:

$p_{i}=[p_{i_{x}},p_{i_{y}},p_{i_{z}}]^{\top}$; $LF=p_{d}=p_{f}-p_{l}$ | Position
---|---
$v_{i}=[v_{i_{x}},v_{i_{y}},v_{i_{z}}]^{\top}$ | Velocity
$a_{i}=[a_{i_{x}},a_{i_{y}},a_{i_{z}}]^{\top}$ | Acceleration
$u_{i}=[u_{x},u_{y},u_{z}]^{\top}$ | Control

###  2.3 Mathematical Formulation

As stated previously, the goal of this validation tool is to determine the
optimal path for a follower aircraft to rejoin to a formation position,
defined by the set of points contained on the formation ring in minimum time,
and subsequently maintain the formation position with least deviation until a
fixed final time. This problem was crafted into a two-phase dynamic optimal
control problem based on the GPOPS-II standard problem formulation [8].

The cost function for phase 1 dictates a minimum time to execute the initial
rejoin from a predefined starting location to a point contained on the
formation ring. The Mayer form of a minimum time cost function was used as the
final time of phase 1. The follower’s dynamics were modeled by three-
dimensional triple integrator dynamics, with the three-element control applied
directly to the follower’s jerk vector. The scenario begins at time zero with
a known initial state of both the leader and follower. A path constraint was
enforced to prevent the follower from entering the leader’s jet wash, a safety
hazard to the follower aircraft. The jet wash was modeled as an infinite
cylinder centered about the $x^{L}$-axis centered at $(y^{L},z^{L})=(0,0)$
with a radius of $R_{jw}$. Note that $t_{f}^{(i)}$ indicates the final time of
phase $i$. Thus, phase 1 included the terminal condition requiring the
follower to attain a formation position at $t_{f}^{(1)}$ which was enforced by
the functions:

$f_{1}(p_{d})=(p_{d_{x}}-[R_{L}^{i}c^{L}]_{x})^{2}$
---
$f_{2}(p_{d})=(p_{d_{y}}-[R_{L}^{i}c^{L}]_{y})^{2}+(p_{d_{z}}-[R_{L}^{i}c^{L}]_{z})^{2}-R^{2}_{ring}$
(2)

The equations in (2) provided a way to characterize the distance from the
follower aircraft to the formation ring. $f_{1}(p_{d}(t_{f}^{(1)}))$ ensured
the follower aircraft attained the same $x^{L}$-coordinate in the leader
reference frame as the tip of the vector $c^{L}$, which was used to define the
center point of the formation ring. $f_{2}(p_{d}(t_{f}^{(1)}))$ was used to
ensure that the $y^{L}$ and $z^{L}$ coordinates of the follower in relation to
the leader were at a distance from the ring center equal to the radius of the
ring. Simply put, when $f_{1}=f_{2}=0$ the follower was established in a valid
formation position defined by the ring of points of radius $R_{ring}$ with
center defined by the vector $c^{L}$. In order to permit small numerical
errors, a tolerance of $\pm 10$ meters was allowed as the maximum value of
each function at the terminal condition. Lastly, constraints were enforced on
the magnitude of the follower’s velocity and acceleration vectors for the
duration of the phase. The magnitude of the follower’s velocity vector was
constrained by a minimum and maximum velocity, as dictated by the aircraft’s
stall speed and maximum speed. The magnitude of the acceleration vector was
constrained by a minimum acceleration to avoid zero g flight, and a maximum
acceleration to prevent overloading the aircraft structure. Both of these
constraints were enforced as path constraints throughout the duration of the
phase.

Minimizing the cost function of phase 2 required the follower aircraft to
minimize the deviation from the formation ring for the duration of the phase.
This deviation was characterized by squaring and summing $f_{1}$ and $f_{2}$
as defined during phase 1, and integrating the resulting values over the
duration of the phase. Phase 2 used the same follower dynamics along with the
same path, state, and control constraints. The second phase also contained the
requirement for the position and velocity to match across phases with a
tolerance of $\pm 0.1$ m, ensuring a continuous trajectory throughout the
entire scenario. Additionally, the condition for the free final time of phase
1 to match the initial time of phase 2 was enforced. The final time of phase 2
was fixed, and was chosen to ensure that the follower could attain the
formation ring during phase 1 before the fixed final time of phase 2.

Since the optimization software required one cost function for the overall
problem, the cost functions from both phase 1 and phase 2 were summed to
produce a cumulative cost function. In order to account for a difference in
the magnitudes between the cost functions of each phase, a weighting parameter
$\beta$ was included in the summation.

Overall: | $\min\limits_{u(t)}$ | $J=\beta J_{1}+(1-\beta)J_{2}$ | Cost Function
---|---|---|---
Phase 1: | $\min\limits_{u(t)}$ | $J_{1}=t_{f}^{(1)}$ | Cost Function
| S.T. | $\dot{p_{f}}=v_{f}$ | Dynamics
|  | $\dot{v_{f}}=a_{f}$ | Dynamics
|  | $\dot{a_{f}}=u$ | Dynamics
|  | $t_{0}=0\text{, }p_{f}(t_{0})=p_{f_{0}}\text{, }v_{f}(t_{0})=v_{f_{0}}\text{, }a_{f}(t_{0})=a_{f_{0}}$ | Boundary Condition
|  | $-\epsilon_{form}\leq f_{1}(p_{d}(t_{f}^{(1)}))=(p_{d_{x}}(t_{f}^{(1)})-[R_{L}^{i}(\gamma,\chi)c^{L}]_{x})^{2}\leq\epsilon_{form}$ | Terminal Condition
|  | $-\epsilon_{form}\leq f_{2}(p_{d}(t_{f}^{(1)}))=(p_{d_{y}}(t_{f}^{(1)})-[R_{L}^{i}(\gamma,\chi)c^{L}]_{y})^{2}$ |
|  | $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +(p_{d_{z}}(t_{f}^{(1)})-[R_{L}^{i}(\gamma,\chi)c^{L}]_{z})^{2}-R_{ring}^{2}\leq\epsilon_{form}$ | Terminal Condition
|  | $-([R_{i}^{L}(\gamma,\chi)p_{d}]_{y})^{2}-([R_{i}^{L}(\gamma,\chi)p_{d}]_{z})^{2}+R_{jw}^{2}\leq 0$ | Jet Wash
|  | $v_{f_{min}}\leq|v_{f}|\leq v_{f_{max}}\text{, }a_{f_{min}}\leq|a_{f}|\leq a_{f_{max}}$ | Constraints
Phase 2: | $\min\limits_{u(t)}$ | $J_{2}=\int\limits_{t_{0}^{(2)}}^{t_{f}^{(2)}}\\{f_{1}(p_{d}(t))^{2}+f_{2}(p_{d}(t))^{2}\\}dt$ | Cost Function
| S.T. | $\dot{p_{f}}=v_{f}$ | Dynamics
|  | $\dot{v_{f}}=a_{f}$ | Dynamics
|  | $\dot{a_{f}}=u$ | Dynamics
|  | $t_{f}^{(2)}=120$ | Boundary Condition
|  | $t_{o}^{(2)}=t_{f}^{(1)}\text{, }p_{f}(t_{0}^{(2)})=p_{f}(t_{f}^{(1)})$ | Phase Continuity
|  | $v_{f}(t_{0}^{(2)})=v_{f}(t_{f}^{(1)})\text{, }a_{f}(t_{0}^{(2)})=a_{f}(t_{f}^{(1)})$ | Phase Continuity
|  | $-([R_{i}^{L}(\gamma,\chi)p_{d}]_{y})^{2}-([R_{i}^{L}(\gamma,\chi)p_{d}]_{z})^{2}+R_{jw}^{2}\leq 0$ | Jet Wash
|  | $v_{f_{min}}\leq|v_{f}|\leq v_{f_{max}}\text{, }a_{f_{min}}\leq|a_{f}|\leq a_{f_{max}}$ | Constraints
Table 1: Initial Conditions $p_{f_{0}}=\begin{bmatrix}0\text{ ft}\\\ 0\text{ ft}\\\ 20,000\text{ ft}\end{bmatrix}$ | $v_{f_{0}}=\begin{bmatrix}450\text{ kts}\\\ 0\text{ kts}\\\ 0\text{ kts}\end{bmatrix}$ | $a_{f_{0}}=\begin{bmatrix}0.07\text{ g}\\\ 0\text{ g}\\\ 1\text{ g}\end{bmatrix}$ | $c^{L}=\begin{bmatrix}-700ft\\\ 0\\\ 0\end{bmatrix}$
---|---|---|---
$R_{ring}$ = 700 ft |  $R_{jw}$=5 ft | $\beta=0.9$ |  $\epsilon_{form}=10$ m
$V_{min}$ = 200 kts |  $V_{max}$ = 700 kts |  $A_{min}$ = 0.25 g |  $A_{max}$ = 7 g

A summary of the problem’s mathematical formulation is presented on the
following page.

###  2.4 Discrete Hamiltonian Derivation

A derivation of the discrete Hamiltonian is provided and was used to ensure
the validity of the optimal solution. The continuous Hamiltonian was defined
as $\mathcal{H}=\mathcal{L}+\bar{\lambda}^{\text{T}}\mathbf{F}$. The
Lagrangian of the cost function was defined as
$\mathcal{L}$$=(1-\beta)\\{f_{1}^{2}+f_{2}^{2}\\}$. $\bar{\lambda}$ contained
a column vector of the 9 costates. Finally, $\mathbf{F}$ was a column vector
of the right-hand side of the 9 dynamics equations. Since the tool produces a
solution at discrete points in time, the discrete form of the Hamiltonian
equation was used, defined as
$\mathcal{H}$${}_{d}=\mathcal{L}$${}_{k}+\bar{\lambda}_{k}^{\text{T}}\mathbf{F}_{k}$.
Here the subscript $k$ denotes the value of each function at the point
$t_{k}$. For a problem in which the Hamiltonian is not an explicit function of
time, the value of the Hamiltonian will be constant on an extremal path [9].
This criterion is applicable to both phase 1 and phase 2.

###  2.5 Test Cases

A summary of the initial conditions and parameters is presented in Table 1.
Both aircraft were modeled as a typical 7g capable fighter aircraft. The
formation ring center was set at 700 ft directly behind the leader aircraft.
The radius of the formation ring was set to 700 ft. This placed the formation
ring points approximately 1,000 ft from the leader aircraft on a 45-degree
aspect, directly correlating to the fighting wing position discussed
previously. The jet wash cylinder was set 5 feet wide. The follower aircraft
was initialized in 1 NM trail of the leader aircraft, with a 1,000 ft lateral
offset to the leader’s right wing. The simulation was run for 2 minutes. The
velocity was limited to a minimum of 200 knots to prevent stall, and a maximum
of 700 knots. The maximum acceleration was limited to 7 g’s, with a minimum of
0.25 g’s.

Figure 3: Initial Guess Provided by Algorithm

In order to assess the capabilities of the validation tool, two test cases
were evaluated. Both aircraft began on a north heading in level flight at 450
knots and 20,000 ft MSL. The first leader trajectory consisted of a constant
left bank, descending spiral. The leader descended with constant airspeed at
approximately 3.25 degrees of pitch losing 6,500 ft over the 2-minute
simulation. The second case used the same initial conditions followed by the
leader performing a number of continuous loops. The loops were approximately
12,000 ft tall starting and ending at the initial altitude. The loops were
executed on a constant heading. Again, a 2-minute simulation time was used.
Due to the minimum time formulation of the problem, the optimal solution was
expected to take a direct path to the formation ring, with allowances for
maintaining the ring position throughout phase 2.

##  3 Solution Methodology

###  3.1 Direct Collocation Solver

This problem was solved with the direct orthogonal collocation solver GPOPS-II
discussed in [8]. GPOPS-II uses a polynomial approximation of the state and
control vectors. The solver transcribed the continuous optimal control problem
by interpolating these approximations at LGR based collocation points as
discussed in [10]. Gaussian quadrature was used to approximate the cost
function, and a differential matrix based on the Lagrange polynomial basis
allowed the dynamics to be enforced exactly at each collocation point.
Additionally, the boundary and path constraints were enforced at each
collocation point. GPOPS-II then formatted the problem for hand-off to a
nonlinear program (NLP) solver. The NLP solver used both the gradient and
Hessian of the transcribed problem. These functions were computed using the
MATLAB program ADiGator discussed in [11].

Scaling was used to aid in solver convergence, which simply consisted of
normalizing to 100s of meters of distance. The NLP solver requires bounds to
constrain the feasible region for all variables. The position state was
bounded to the max distance the follower aircraft could travel straight along
the coordinate direction at maximum velocity over the two-minute simulation
time. The velocity and acceleration components were bounded by $\pm V_{max}$
and $\pm A_{max}$ respectively. An acceleration rate limit of $\pm 5$ g/s was
used to bound the control components. Unless specified elsewhere, default
settings were used as provided with GPOPS-II.

Once the first solution was returned from the NLP solver, the error in the
enforcement of the dynamics at the midpoints between collocation points was
estimated by GPOPS-II. The maximum magnitude of this error was used as the
termination criteria for the algorithm, with a tolerance of $10^{-3}$. If the
error exceeded the tolerance, GPOPS-II modified the mesh of each phase, and
the number of collocation points included on each mesh interval. The hp-LiuRao
adaptive mesh refinement method provided with GPOPS-II was used. This process
was then repeated until the mesh tolerance was satisfied. GPOPS-II reported
the final solution with values of the states, controls, and costate estimates
at each collocation point and one interpolated point at the end of each of the
two phases.

###  3.2 NLP Solver

The NLP Solver IPOPT was used as provided with GPOPS-II. IPOPT used an
interior point boundary value method which was described in [12]. A tolerance
of $10^{-5}$ was used for the determination of convergence within IPOPT. The
linear solver MA57 was used by IPOPT to solve the system of equations on each
NLP iteration.

###  3.3 Initial Guess Generation

GPOPS-II required an initial guess for the states and controls over the
duration of the problem. In order to align with future work, the initial guess
of the follower’s states and controls were generated from a control strategy
algorithm. This algorithm was currently under development with the intent of
creating the ability to approximate the optimal trajectory in real-time.
Details of this algorithm will be documented in a future publication. The
algorithm used a control strategy which selected from simple maneuver
primitives based on various decision criteria informed by the relation between
the leader and follower aircraft at a given instance in time. While this
algorithm is still under development, it has a basic capability such that the
follower aircraft can track to a single point defined in reference to the
leader aircraft.

Using this basic capability, the algorithm was run such that the follower
tracked to and maintained a point that would coincide with a single point
located on the formation ring. This algorithm did not allow for the tracking
of a set of points at the time of this writing, as the validation tool does.
For both cases the initial guess was generated by having the follower aircraft
attempt to track to the point that is longitudinally in line with the
projection of the right-wing tip of the leader. While this algorithm was still
under development, it provided a sufficient initial guess such that the
optimal control solver was able to converge to an optimal trajectory in all
cases.

In order to assess the impact various initial guesses may have had on the
convergence of the optimal control solver, the spiral case was evaluated with
4 different algorithm results. The algorithm was set to track to 4 various
points along the formation ring. The four points were located along the
projection of the left wing, along the projection of the right wing, and at
the top and bottom of the ring when referencing the $z^{L}$ axis. The
algorithm results were then evaluated on their time to rejoin to the formation
ring and the magnitude of the final cost function for comparison. An example
of an initial guess for the spiral trajectory can be seen in figure 3.

Figure 4: Position Response for Case 1

##  4 Results

###  4.1 Spiral Trajectory

For the descending spiral trajectory, GPOPS-II reported an optimal solution
for which the states and controls are depicted in Figures 4, 5, 6, and 7. The
position response showed the follower aircraft in black, the leader aircraft
in blue, and a maroon diamond represented the point in the trajectory when the
phase transition occurs. As expected, the follower took a direct path to the
formation ring during phase 1, while maintaining the formation ring nearly
perfectly throughout the duration of phase 2. It should be noted that some
motion along the formation ring was observed throughout phase 2, but as this
incurs no penalty within the cost function this does not preclude this path
from being optimal. The velocity and acceleration response showed the aircraft
attempting to achieve a maximum velocity in order to minimize the rejoin time.
The limits on the jerk control vector could also be seen by the linear changes
in the acceleration components. The optimal control is very noisy during phase
1. However, this was the main reason to use triple integrator dynamics as this
noise simply represents a g onset rate which wouldn’t be prohibitively
difficult to implement on a real aircraft.

The path constraints imposed on the magnitudes of velocity and acceleration
can be seen in Figures 8 and 9. These plots showed the follower aircraft
achieving its maximum velocity at the maximum acceleration rate throughout
phase 1, reducing the rejoin time to the greatest extent possible. The rejoin
time reduction was further enabled by the ability of the follower to
maneuvering along the formation ring and dissipate excess energy during phase
2 with no additional cost penalty. This reduction in rejoin time would not
have been possible if the formation position was represented by a single point
as the follower would have to arrive at that point with its velocity already
matched to the leader’s velocity in order to minimize deviation during phase
2.

Figure 5: Velocity Response for Case 1 Figure 6: Acceleration Response for
Case 1 Figure 7: Optimal Control for Case 1 Figure 8: Magnitude of Velocity
for Case 1 Figure 9: Magnitude of Acceleration for Case 1 Figure 10: Deviation
from Formation Ring in Case 1

The final plot seen in figure 10 for case 1 shows the progression of the ring
deviation functions throughout the simulation. As expected, the follower
rapidly reduced the distance to zero during phase 1, and maintained a zero
deviation throughout phase 2 even as the leader continued to maneuver.

Figure 11: Position Response for Case 2

###  4.2 Loop Trajectory

For the leader trajectory consisting of continuous vertical loops, the optimal
path can be seen in Fig 11. This leader trajectory was chosen for study due to
its significant differences from the descending spiral. The leader trajectory
contained significant vertical deviations along with large changes in the
leader’s velocity. However, similar behavior of the optimal trajectory was
seen across the two cases. During phase 1, the follower again accelerated at
the maximum rate to maximum velocity in a direct path to the formation ring.
The follower was then able to maintain the ring with little deviation for the
duration of phase 2.

Figures 12-16 show the other state responses for the follower’s trajectory.
The maximum velocity attained during phase 1 can be clearly seen, and the
large variations in the vertical direction match with the loop maneuvers of
the leader. Acceleration rate limiting was also observed in the same fashion
as the spiral trajectory. The follower still maneuvered along the formation
ring during phase 2 while dissipating excess energy from the rejoin, but with
little deviation, as can be seen in Figure 17.

Table 2: Initial Guess Cost Comparison | Left Wing | Right Wing | Top | Bottom
---|---|---|---|---
$t_{f}^{(1)}$ | 15.2319 s | 15.2320 s | 15.2422 s | 15.2362 s
$J$ | 16.1472 | 16.1471 | 16.1221 | 16.1221

###  4.3 Initial Guess Comparison

The second phase of this problem allows for motion along the formation ring
without penalty. Thus, more than one path could be considered equally optimal
during the second phase. A study of this behavior was conducted by varying the
initial guess given to GPOPS-II. Figure 18 shows the resulting trajectories
that were reported as optimal by the software for the 4 tested initial
conditions. As predicted, the trajectories were nearly identical during phase
one as they took a direct path to minimize the time of the rejoin. During the
later portion of phase 2 the trajectories all settled to stable locations
along the ring with only small deviations observed between the paths. However,
when the follower was dissipating excess energy carried from the first phase,
the paths diverged.

Table 2 shows the final times of phase 1 and overall cost function values for
each initial condition. It is notable that the cost functions are nearly
identical even with the path deviations that occurred at the beginning of
phase 2. While the various initial conditions caused the solver to converge to
somewhat different phase 2 paths, the solver still converged to paths that
were equally optimal in regards to the cost function that was studied. This
result brings interest to future studies which could include further penalties
in the cost function such as limiting control usage or energy expenditure.
Figure 19 shows the similarity in the velocity profiles between each initial
condition. This figure uses the same color scheme as Figure 18. Similar
results were also seen for the other states and the control vector.

Figure 12: Velocity Response for Case 2 Figure 13: Acceleration Response for
Case 2 Figure 14: Optimal Control for Case 2 Figure 15: Magnitude of Velocity
for Case 2 Figure 16: Magnitude of Acceleration for Case 2 Figure 17:
Deviation from Formation Ring in Case 2

Figure 18: Initial Condition Trajectory Comparisons

In all cases the tool produced feasible trajectories with no violation of the
path constraints. This was observed for the jet wash constraint by noting that
the value of the constraint never took on a positive value. However, the
maximum value of the jet wash constraint, $-8.54\times 10^{-5}$, occurred
during the spiral trajectory with the goal position at the top of the ring So,
the optimal trajectory did closely approach the jet wash cylinder, but never
violated the path constraint.

Figure 19: Initial Condition Velocity Comparisons Figure 20: Discrete
Hamiltonian

###  4.4 Analysis of Optimally

The discrete Hamiltonian is presented for both leader trajectories in Figure
20. As previously discussed, the Hamiltonian should be constant during phase 1
and 2. Both figures show that using the costate estimates provided by GPOPS-II
did yield Hamiltonians which were consistent with these conditions. There was
some variation near the phase switch likely caused by numerical error.
However, the magnitude of this error was not observed to have significantly
affected the validity of the solutions.

These results demonstrate the initial ability of the tool to produce optimal
trajectories for a variety of complex trajectories.

##  5 Conclusions

This paper demonstrated the initial capability of a validation tool which can
find the optimal path for a follower aircraft to rejoin to a set of dynamic
formation points in minimum time and maintain that position for a given leader
trajectory. Of note was the ability of the follower to use the increased
flexibility of the formation ring to decrease the rejoin time during the first
phase with no penalty during the second phase of the trajectory. It allows for
a follower to maintain maximum flexibility while maintaining formation
position. This in turn allows the leader aircraft freedom to perform
aggressive maneuvering if required. Additionally, the solver showed consistent
behavior for two drastically different leader trajectories. The optimal
trajectory worked to quickly rejoin to the formation ring, and maintained the
ring with nearly zero deviation during phase 2, regardless of the leader
maneuvering laterally or vertically. Additionally, the validity of the
solutions was shown through the consistent behavior of the discrete
Hamiltonian. This tool should prove useful in future research for finding the
optimal path to evaluate the accuracy of a control strategy.

Creating this validation tool uncovered a number of lessons learned for future
use. The tool produced a noisy control response, which necessitated the use of
triple integrator dynamics to yield realistic acceleration and velocity
vectors. While choosing different initial conditions did produce somewhat
different optimal trajectories, they did not differ significantly in their
rejoin time or total amount of deviation from the formation ring. This result
demonstrated that there was a family of optimal solutions which were all
approximately equal in value for the given cost function. Future research
could consider adding additional cost parameters to penalize things such as
excess control effort or energy expenditure.

A number of areas are planned for future research. The first area was planned
to test additional leader trajectories to continue building trust in the
tool’s ability to generate optimal trajectories. A more complicated dynamics
model will be incorporated for the follower’s dynamics, incorporating aircraft
aerodynamics. Different formation shapes can be defined by modifying the cost
functions accordingly. However, the next major area of research will consist
of determining control laws and strategies which can approximate these optimal
trajectories, but can be computed at a speed which will allow for real-time
implementation.

###### Acknowledgements.

The authors thank the Air Force Research Laboratories for funding this
project.

## References

  * [1] T.-C. F. Fundamentals, “Air education and training command manual (aetcman) 11-251,” Tech. Rep., 2017.
  * [2] L. E. Buzogany, M. Pachter, and J. J. D’Azzo, “Automatic control of aircraft in formation flight,” in _Guidance Navigation and Control Conference_. Monterey, CA, USA: AIAA, 1993\.
  * [3] M. Pachter, J. J. D’Azzo, and J. L. Dargan, “Automatic formation flight control,” _Journal of Guidance, Control, and Dynamics_ , vol. 17, no. 6, pp. 1380–1383, 1994.
  * [4] V. P. Reyna, “Automation of formation flight control,” Wright Patterson AFB, OH, 1994.
  * [5] I. E. Weintraub, Z. J. L. Demers, J. Shroyer, N. P. Ritsema, and D. P. Sluss, “Circular rejoin in 3d using bézier paths,” in _IEEE Conference on Control Technology and Applications_ , Montreal, Canada, 2020.
  * [6] D. Tran, D. Casbeer, E. Garcia, I. E. Weintraub, and D. Milutinović, “Ring formation maneuver: Double-integrator kinematics with input saturation,” _Journal of Guidance, Control, and Dynamics_ , vol. 44, no. 11, pp. 2091–2098, 2021.
  * [7] R. W. Beard and T. W. McLain, _Small Unmanned Aircraft: Theory and Practice_. Princeton, NJ: Princeton University Press, 2012.
  * [8] M. A. Patterson and A. V. Rao, “GPOPS - II: A MATLAB software for solving multiple-phase optimal control problems using hp-adaptive gaussian quadrature collocation methods and sparse nonlinear programming,” _ACM Transactions on Mathematical Software_ , vol. 41, no. 1, pp. 1–37, 2014.
  * [9] D. Kirk, _Optimal Control Theory: An Introduction_. Englewood Cliffs, NJ: Prentice-Hall, 1970, Section 5.1, pp.188.
  * [10] D. Garg, M. Patterson, W. W. Hager, A. V. Rao, D. A. Benson, and G. T. Huntington, “A unified framework for the numerical solution of optimal control problems using pseudospectral methods,” _Automatica_ , vol. 46, no. 11, pp. 1843–1851, 2010.
  * [11] M. J. Weinstein and R. A. V., “Algorithm 984: ADiGator, a toolbox for the algorithmic differentiation of mathematical functions in matlab using source transformation via operator overloading,” _ACM Transactions on Mathematical Software_ , vol. 44, no. 2, pp. 1–25, 2018.
  * [12] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” _Mathematical Programming_ , vol. 106, pp. 25–57, 2006.

Carl GotwaldOfficialPhotoSeptember22.png received his B.S. from the United
States Air Force Academy. He holds M.S. degrees from both Mississippi State
University and Air University, and is currently pursuing a Ph.D. at the Air
Force Institute of Technology. He works as a test pilot with extensive
experience testing remotely piloted aircraft and instructing at the United
States Air Force Test Pilot School. His current research interests lie in
optimal control and real-time simulation.

Michael ZollarsZollars_Picture.png received a B.S. in Mechanical Engineering
from the Pennsylvania State University in 2003 and an M.S. and Ph.D. in
Aeronautical Engineering from the Air Force Institute of Technology in 2007
and 2018 respectively. He is a professor with the Department of Aeronautics
and Astronautics at the Air Force Institute of Technology. His research
interests include optimal control for vehicle trajectory analysis and dynamics
and control of aircraft systems.

Isaac E. WeintraubIsaac_Photo.png (S’09-M’15-SM’21) holds a Ph.D. from The Air
Force Institute of Technology (2021), an M.S. in Electrical Engineering from
University of Texas at Arlington (2011), and a B.S. in Mechanical Engineering
from Rose-Hulman Institute of Technology (2009). He is currently an
Electronics Engineer with the Control Science Center, Air Force Research
Laboratory, Wright-Patterson Air Force Base, Dayton, OH, USA. His research
interests are in automation and control of aerospace systems for defense
applications.