#  Skew Dyck paths without up–down–left

Helmut Prodinger  Helmut Prodinger, Department of Mathematical Sciences,
Stellenbosch University, 7602 Stellenbosch, South Africa, and NITheCS
(National Institute for Theoretical and Computational Sciences), South Africa
<EMAIL_ADDRESS>
###### Abstract.

Skew Dyck paths without up–down–left are enumerated. In a second step, the
number of contiguous subwords ‘up–down–left’ are counted. This explains and
extends results that were posted in the Encyclopedia of Integer Sequences.

###### 1991 Mathematics Subject Classification:

05A15

##  1\. Introduction

The entry A128729 in [5] is a bit mysterious. It presents some results about
skew Dyck paths without a (contiguous) sequence ‘up–down–left,’ but without
justification, and without reference to a research paper. The present note is
designed to fill these gaps. First, skew Dyck paths are Dyck paths with
additional left steps $(-1,-1)$ but without overlapping itself. We found it
more convenient [4] to replace a left step by a red down-step $(1,-1)$,
labelled red. So we might say that up–down–red is forbidden. In a second step,
it won’t be forbidden, but counted how often this happens. If it does not
happen, it is the same as ‘forbidden’. Probably the first paper that deals
with skew Dyck paths is [1].

The Figure 1 describes skew Dyck paths. One can compute more: Paths that end
at level $k$, by any of the 3 types of steps. Note that up–red and red–up are
forbidden, as they would overlap itself.

Figure 1. Three layers of states according to the type of steps leading to
them (up, down-black, down-red).

The next step is to take care of up–down–red. Such a red step is depicted in
Figure 2 in a special color.

Figure 2. Four layers of states according to the type of steps leading to
them.

##  2\. The magenta colored steps are forbidden

We introduce 4 sequences of generating functions, according to the 4 layers in
the defining graph. The variable $z$ counts the steps. The recursions can be
seen immediately from the graph:

$\displaystyle f_{n+1}$ $\displaystyle=zf_{n}+zg_{n}+zh_{n},\ f_{0}=1,$
$\displaystyle g_{n}$ $\displaystyle=zf_{n+1},$ $\displaystyle h_{n}$
$\displaystyle=zg_{n+1}+zh_{n+1}+zk_{n+1},$ $\displaystyle k_{n}$
$\displaystyle=zh_{n+1}+zk_{n+1}.$

And now we use a second variable $u$ to be able to work with bivariate
generating functions;

$\displaystyle\sum_{n\geq 0}u^{n+1}f_{n+1}$ $\displaystyle=\sum_{n\geq
0}zu^{n+1}f_{n}+\sum_{n\geq 0}zu^{n+1}g_{n}+\sum_{n\geq 0}zu^{n+1}h_{n},$
$\displaystyle\sum_{n\geq 0}u^{n+1}g_{n}$ $\displaystyle=\sum_{n\geq
0}zu^{n+1}f_{n+1},$ $\displaystyle\sum_{n\geq 0}u^{n+1}h_{n}$
$\displaystyle=\sum_{n\geq 0}zu^{n+1}g_{n+1}+\sum_{n\geq
0}zu^{n+1}h_{n+1}+\sum_{n\geq 0}zu^{n+1}k_{n+1},$ $\displaystyle\sum_{n\geq
0}u^{n+1}k_{n}$ $\displaystyle=\sum_{n\geq 0}zu^{n+1}h_{n+1}+\sum_{n\geq
0}zu^{n+1}k_{n+1}.$

Translating this into bivariate generating functions $F(u)=F(u,z)$ etc., leads
to

$\displaystyle F(u)$ $\displaystyle=1+zuF(u)+zuG(u)+zuH(u),$ $\displaystyle
uG(u)$ $\displaystyle=zF(u)-z,$ $\displaystyle uH(u)$
$\displaystyle=zG(u)-zg_{0}+zH(u)-zh_{0}+zK(u)-zk_{0},$ $\displaystyle uK(u)$
$\displaystyle=zH(u)-zh_{0}+zK(u)-zk_{0}.$

Solving the system leads to

$\displaystyle F(u)$
$\displaystyle=\frac{{u}^{2}{z}^{2}g_{0}+{u}^{2}{z}^{2}h_{0}+{u}^{2}{z}^{2}k_{0}-u{z}^{3}+2zu-{u}^{2}+{u}^{2}{z}^{2}-u{z}^{3}g_{0}-{z}^{4}}{-{u}^{2}+z{u}^{3}+2zu-{u}^{2}{z}^{2}-u{z}^{3}-{z}^{4}},$
$\displaystyle G(u)$
$\displaystyle=\frac{{z}^{2}\left(uzg_{0}+uzh_{0}+uzk_{0}+2zu-{z}^{2}g_{0}-{u}^{2}\right)}{-{u}^{2}+z{u}^{3}+2zu-{u}^{2}{z}^{2}-u{z}^{3}-{z}^{4}},$
$\displaystyle H(u)$
$\displaystyle=\frac{-z\left(u{z}^{2}-{z}^{3}-g_{0}u-h_{0}u-k_{0}u+zg_{0}+z{u}^{2}g_{0}+{u}^{2}zh_{0}+{u}^{2}zk_{0}+u{z}^{2}h_{0}+u{z}^{2}k_{0}-{z}^{3}g_{0}\right)}{-{u}^{2}+z{u}^{3}+2zu-{u}^{2}{z}^{2}-u{z}^{3}-{z}^{4}},$
$\displaystyle K(u)$
$\displaystyle=\frac{-\left({u}^{2}zk_{0}+{u}^{2}zh_{0}-h_{0}u+u{z}^{2}g_{0}-k_{0}u+u{z}^{2}k_{0}+u{z}^{2}h_{0}-zg_{0}+{z}^{3}+{z}^{3}g_{0}+k_{0}{z}^{3}+h_{0}{z}^{3}\right)z}{-{u}^{2}+z{u}^{3}+2zu-{u}^{2}{z}^{2}-u{z}^{3}-{z}^{4}},$

One cannot just plug in $u=0$ to identify the constants. However,

$-{u}^{2}+z{u}^{3}+2zu-{u}^{2}{z}^{2}-u{z}^{3}-{z}^{4}=z(u-u_{1})(u-u_{2})(u-u_{3}),$

and the bad factors $(u-u_{2})(u-u_{3})$ can be cancelled out, both, in the
denominator and the numerator. This is essential for the kernel method [3].
This leads to

$\displaystyle F(u)$
$\displaystyle=\frac{z^{2}g_{0}+z^{2}h_{0}+z^{2}k_{0}-1+z^{2}}{z(u-u_{1})},$
$\displaystyle G(u)$ $\displaystyle=\frac{-z}{(u-u_{1})},$ $\displaystyle
H(u)$ $\displaystyle=\frac{-z(g_{0}+h_{0}+k_{0})}{(u-u_{1})},$ $\displaystyle
K(u)$ $\displaystyle=\frac{-z(h_{0}+k_{0})}{(u-u_{1})}.$

Now we can plug in $u=0$ and solve to get

$\displaystyle g_{0}$ $\displaystyle=\frac{z}{u_{1}},$ $\displaystyle h_{0}$
$\displaystyle=\frac{1-z^{2}-zu_{1}}{zu_{1}},$ $\displaystyle k_{0}$
$\displaystyle=\frac{1-z^{2}-zu_{1}}{u_{1}(u_{1}-z)}.$

In total

$\displaystyle 1+g_{0}+h_{0}+k_{0}$
$\displaystyle=\frac{u_{1}-z^{2}u_{1}-z^{3}}{zu_{1}(u_{1}-z)}=\frac{1-zu_{1}}{z^{2}}$
$\displaystyle=1+z^{2}+2z^{4}+6z^{6}+20z^{8}+71z^{10}+262z^{12}+994z^{14}+3852z^{16}+\cdots.$

The last simplification was done by noticing that
$u_{1}=\textsf{RootOf}(-u^{2}+zu^{3}+2zu-u^{2}z^{2}-uz^{3}-z^{4},u)$ and using
Maple’s evala command. (One can also use the algeqtoseries command to go from
the algebraic equation to the series expansion.)

We further find

$F+G+H+K=-\frac{u_{1}(1-zu_{1})}{z^{2}(u-u_{1})}$

and

$[u^{k}](F+G+H+K)=\frac{1-zu_{1}}{z^{2}u_{1}^{k}}.$

We summarize the results of this section:

######  Theorem 1.

The generating function of skew Dyck paths without up–down–left ending on
level $k$ is given by

$\frac{1-zu_{1}}{z^{2}u_{1}^{k}}$

where

$u_{1}=\frac{1}{z}-z-z^{3}-2z^{5}-6z^{7}-20z^{9}-71z^{11}-262z^{13}-\cdots$

is the solution of the algebraic equation
$-u^{2}+zu^{3}+2zu-u^{2}z^{2}-uz^{3}-z^{4}=0$ that does not contain imaginary
numbers in its series expansion around $z=0$. Equivalently, the other
solutions $u_{2}$, $u_{3}$ are such that $1/(u-u_{2})$ and $1/(u-u_{3})$ have
no power series expansion around $(0,0)$.

Noticing that one can only return to level 0 in an even number of steps, one
can write an alternative form for $f_{0}+g_{0}+h_{0}+k_{0}$:

RootOf $\displaystyle(z^{2}u^{3}-z(2-z)u^{2}+(1-z^{2})u-1+z+z^{2},u)$
$\displaystyle=1+z+2z^{2}+6z^{3}+20z^{4}+71z^{5}+262z^{6}+994z^{7}+3852z^{8}+\cdots.$

We show how the algebraic equation for $u_{1}$ is transformed into one for
$\frac{1-zu_{1}}{z^{2}}$; in the last step, $z^{2}=Z$ was substituted:

$\displaystyle
u_{1}\longrightarrow-u^{2}+zu^{3}+2zu-u^{2}z^{2}-uz^{3}-z^{4}=0,$
$\displaystyle
zu_{1}\longrightarrow-(zu)^{2}+(zu)^{3}+2z^{2}(zu)-z^{2}(zu)^{2}-z^{4}(zu)-z^{6}=0,$
$\displaystyle-
zu_{1}\longrightarrow-(-zu)^{2}-(-zu)^{3}-2z^{2}(-zu)-z^{2}(-zu)^{2}+z^{4}(-zu)-z^{6}=0,$
$\displaystyle
1-zu_{1}\longrightarrow-(1-zu-1)^{2}-(1-zu-1)^{3}-2z^{2}(1-zu-1)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-z^{2}(1-zu-1)^{2}+z^{4}(1-zu-1)-z^{6}=0,$
$\displaystyle\qquad=2(1-zu)^{2}-(1-zu)-(1-zu)^{3}+z^{2}-z^{2}(1-zu)^{2}+z^{4}(1-zu)-z^{4}-z^{6}=0,$
$\displaystyle\frac{1-zu_{1}}{z^{2}}=U_{1}\longrightarrow
2z^{4}U^{2}-z^{2}U-z^{6}U^{3}+z^{2}-z^{6}U^{2}+z^{6}U-z^{4}-z^{6}=0,$
$\displaystyle U_{1}\longrightarrow
2z^{2}U^{2}-U-z^{4}U^{3}+1-z^{4}U^{2}+z^{4}U-z^{2}-z^{4}=0,$ $\displaystyle
U_{1}\longrightarrow 2ZU^{2}-U-Z^{2}U^{3}+1-Z^{2}U^{2}+Z^{2}U-Z-Z^{2}=0.$

So, writing $S(z)=f_{0}+g_{0}+h_{0}+k_{0}$, the algebraic equation is

$z^{2}S^{3}-z(2-z)S^{2}+(1-z^{2})S-1+z+z^{2}=0,$

this is the form given in A128729 [5], and $z$ counts only the half-length of
the skew Dyck path, and Maple’s gfun translates that into a differential
equation

$31z-8-15zS-(2z-1)(44z^{3}+15z^{2}-48z+8)S^{\prime}-z(11z^{2}+16z-4)(2z-1)^{2}S{{}^{\prime\prime}}=0,$

from which gfun derives a linear recursion with polynomial coefficients:

$-44n(n+1)s_{n}-2(n+1)(10n-7)s_{n+1}\\\
+3(115+106n+23n^{2})s_{n+2}-32(n+4)(n+3)s_{n+3}+4(n+5)(n+4)s_{n+4}=0,$

where

$S(z)=\sum_{n\geq 0}s_{n}z^{n}.$

### Asymptotics

The technique is _singularity analysis of generating functions_ , as described
in [2].

We consider $eq:=z^{2}S^{3}-z(2-z)S^{2}+(1-z^{2})S-1+z+z^{2}$ and find the
dominant singularity $(z_{0},S_{0})$. For that, we consider the equation

$\frac{eq}{dS}=3z^{2}S^{2}-2z(2-z)S+1-z^{2}=0$

and thus $S_{0}=\frac{z+1}{3z}$. Plugging this into $eq$, we find
$z_{0}=\frac{2}{11}(3\sqrt{3}-4)$ as closest solution to the origin. From
this, $S_{0}=1+\frac{\sqrt{3}}{2}$. Expanding $eq$ locally,

$z-z_{0}\sim\frac{216-129\sqrt{3}}{121}\Big{(}S-1-\frac{\sqrt{3}}{2}\Big{)}^{2},$

or

$(z_{0}-z)\Big{(}8+\frac{43}{9}\sqrt{3}\Big{)}=\Big{(}2+\frac{8}{9}\sqrt{3}\Big{)}\Big{(}1-\frac{z}{z_{0}}\Big{)}\sim\Big{(}S-1-\frac{\sqrt{3}}{2}\Big{)}^{2}.$

This gives us the local expansion around the square-root singularity:

$S\sim
1+\frac{\sqrt{3}}{2}-\sqrt{2+\frac{8}{9}\sqrt{3}}\;\sqrt{1-\frac{z}{z_{0}}}.$

Therefore

$[z^{n}]S\sim\sqrt{2+\frac{8}{9}\sqrt{3}}\;\frac{1}{2\sqrt{\pi}}\Big{(}2+\frac{3}{2}\sqrt{3}\Big{)}^{n}n^{-3/2}.$

##  3\. Counting the up–down–red configurations

This is not too different from what we did before, but now we use another
variable, $t$, which is attached to the magenta colored edges. In this way,
the exponent of $t$ refers to the number of up–down–red configurations;

$\displaystyle f_{n+1}$ $\displaystyle=zf_{n}+zg_{n}+zh_{n},\ f_{0}=1,$
$\displaystyle g_{n}$ $\displaystyle=zf_{n+1},$ $\displaystyle h_{n}$
$\displaystyle=zg_{n+1}+zh_{n+1}+zk_{n+1},$ $\displaystyle k_{n}$
$\displaystyle=ztg_{n+1}+zh_{n+1}+zk_{n+1}.$

We use again a second variable $u$ to be able to work with bivariate
generating functions;

$\displaystyle\sum_{n\geq 0}u^{n+1}f_{n+1}$ $\displaystyle=\sum_{n\geq
0}zu^{n+1}f_{n}+\sum_{n\geq 0}zu^{n+1}g_{n}+\sum_{n\geq 0}zu^{n+1}h_{n},$
$\displaystyle\sum_{n\geq 0}u^{n+1}g_{n}$ $\displaystyle=\sum_{n\geq
0}zu^{n+1}f_{n+1},$ $\displaystyle\sum_{n\geq 0}u^{n+1}h_{n}$
$\displaystyle=\sum_{n\geq 0}zu^{n+1}g_{n+1}+\sum_{n\geq
0}zu^{n+1}h_{n+1}+\sum_{n\geq 0}zu^{n+1}k_{n+1},$ $\displaystyle\sum_{n\geq
0}u^{n+1}k_{n}$ $\displaystyle=\sum_{n\geq 0}ztu^{n+1}g_{n+1}+\sum_{n\geq
0}zu^{n+1}h_{n+1}+\sum_{n\geq 0}zu^{n+1}k_{n+1}.$

Translating this into bivariate generating functions $F(u)=F(u,z)$ etc., leads
to

$\displaystyle F(u)$ $\displaystyle=1+zuF(u)+zuG(u)+zuH(u),$ $\displaystyle
uG(u)$ $\displaystyle=zF(u)-z,$ $\displaystyle uH(u)$
$\displaystyle=zG(u)-zg_{0}+zH(u)-zh_{0}+zK(u)-zk_{0},$ $\displaystyle uK(u)$
$\displaystyle=ztG(u)-ztg_{0}+zH(u)-zh_{0}+zK(u)-zk_{0}.$

Solving the system leads to

$\displaystyle F(u)$
$\displaystyle=\frac{\mathscr{F}}{2zu-z^{4}-u^{2}+z^{4}t-uz^{3}+zu^{3}-u^{2}z^{2}},$
$\displaystyle G(u)$
$\displaystyle=\frac{\mathscr{G}}{2zu-z^{4}-u^{2}+z^{4}t-uz^{3}+zu^{3}-u^{2}z^{2}},$
$\displaystyle H(u)$
$\displaystyle=\frac{\mathscr{H}}{2zu-z^{4}-u^{2}+z^{4}t-uz^{3}+zu^{3}-u^{2}z^{2}},$
$\displaystyle K(u)$
$\displaystyle=\frac{\mathscr{K}}{2zu-z^{4}-u^{2}+z^{4}t-uz^{3}+zu^{3}-u^{2}z^{2}},$

with
$\mathscr{F}={u}^{2}{z}^{2}h_{0}+{z}^{3}utg_{0}+2zu-{u}^{2}+{z}^{4}t-{z}^{4}+{u}^{2}{z}^{2}k_{0}+{u}^{2}{z}^{2}-u{z}^{3}+{u}^{2}{z}^{2}g_{0}-u{z}^{3}g_{0}$,
$\mathscr{G}={z}^{3}uh_{0}+{z}^{4}g_{0}t+{z}^{3}uk_{0}+2u{z}^{3}+u{z}^{3}g_{0}-{z}^{4}g_{0}-{u}^{2}{z}^{2}$,
$\mathscr{H}={z}^{2}tg_{0}-{z}^{3}utg_{0}-{z}^{4}g_{0}t+{z}^{4}-u{z}^{3}-{u}^{2}{z}^{2}h_{0}-{z}^{4}t+uzg_{0}+uzh_{0}+uzk_{0}-{z}^{2}g_{0}-{u}^{2}{z}^{2}k_{0}-{u}^{2}{z}^{2}g_{0}-{z}^{3}uh_{0}-{z}^{3}uk_{0}+{z}^{4}g_{0}$,
$\mathscr{K}=uztg_{0}-{z}^{2}{u}^{2}tg_{0}-{z}^{2}tg_{0}+{z}^{4}g_{0}t+{z}^{4}h_{0}t+{z}^{4}k_{0}t-{z}^{4}-{z}^{3}ut+uzh_{0}+uzk_{0}+{z}^{2}g_{0}-{z}^{4}g_{0}-u{z}^{3}g_{0}-{z}^{4}h_{0}-{z}^{3}uh_{0}-{u}^{2}{z}^{2}h_{0}-{z}^{4}k_{0}-{z}^{3}uk_{0}-{u}^{2}{z}^{2}k_{0}+{z}^{4}t$.

The denominator $2zu-z^{4}-u^{2}+z^{4}t-uz^{3}+zu^{3}-u^{2}z^{2}$ still
factors as $z(u-u_{1})(u-u_{2})(u-u_{3})$ where $u_{1}$, $u_{2}$, $u_{3}$ now
depend on $t$. Again, the factor $(u-u_{2})(u-u_{3})$ can be cancelled out:

$\displaystyle F(u)$
$\displaystyle=\frac{z^{2}g_{0}+z^{2}h_{0}+z^{2}k_{0}-1+z^{2}}{z(u-u_{1})},$
$\displaystyle G(u)$ $\displaystyle=\frac{-z}{(u-u_{1})},$ $\displaystyle
H(u)$ $\displaystyle=\frac{-z(g_{0}+h_{0}+k_{0})}{(u-u_{1})},$ $\displaystyle
K(u)$ $\displaystyle=\frac{-z(tg_{0}+h_{0}+k_{0})}{(u-u_{1})}.$

Now we can plug in $u=0$ and compute the constants:

$\displaystyle g_{0}$ $\displaystyle=\frac{z}{u_{1}},$ $\displaystyle h_{0}$
$\displaystyle=\frac{1-z^{2}-zu_{1}}{zu_{1}},$ $\displaystyle k_{0}$
$\displaystyle=\frac{(1-z^{2}-zu_{1})(tu_{1}-tz+z)}{zu_{1}(u_{1}+tz-z)}.$

In total

$1+g_{0}+h_{0}+k_{0}=\frac{-tzu_{1}^{2}+tu_{1}+z^{3}t-z^{2}u_{1}-z^{3}+u_{1}}{zu_{1}(u_{1}+tz-z)}=\frac{1-zu_{1}}{z^{2}}.$

The series starts as

$1+z^{2}+(2+t)z^{4}+(4t+6)z^{6}+(16t+20)z^{8}+(71+64t+2t^{2})z^{10}+(262+261t+20t^{2})z^{12}+\cdots.$

One can compute the triple generating function, where $z$ refers to the length
of the walk, $u$ to the level at the end, and $t$ to the number of up–down–red
configurations:

$F(u)+G(u)+H(u)+K(u)=\frac{z^{4}t-u_{1}^{2}+tz^{2}u_{1}^{2}+z^{3}u_{1}+z^{2}u_{1}^{2}-2u_{1}tz-t^{2}z^{4}}{u_{1}(u_{1}+tz-z)z(u-u_{1})}=\frac{u_{1}(zu_{1}-1)}{z^{2}(u-u_{1})}.$

Consequently

$[u^{k}](F(u)+G(u)+H(u)+K(u))=\frac{1-zu_{1}}{z^{2}u_{1}^{k}}.$

Since only an even number of steps brings us to level 0, we have

$z^{2}R^{3}-z(2-z)R^{2}+(1-z^{2})R-1+z+z^{2}-tz^{2}=0$

for $1+g_{0}+h_{0}+k_{0}$ where $z$ counts the half-length of the skew Dyck
path. This is the form given in sequence A128728 [5].

## References

  * [1] Emeric Deutsch, Emanuele Munarini, Simone Rinaldi.  Skew Dyck paths.  J. Stat. Plann. Infer., 140 (8) (2010) 2191–2203.
  * [2] P. Flajolet and A. M. Odlyzko.  Singularity analysis of generating functions.  SIAM J. Discret. Math. 3:216–240, 1990.
  * [3] H. Prodinger,  The kernel method: A collection of examples,  Sém. Lothar. Combin., B50f (2004), 19 pages.
  * [4] H. Prodinger, Partial skew Dyck paths—a kernel method approach, preprint, 2021.
  * [5] N. J. A. Sloane et al., The On-line Encyclopedia of Integer Sequences, 2022.


# Parameter Space Design of Repetitive Controllers for Satisfying a Robust
Performance Requirement

Burak Demirel1, and Levent Güvenç12  2 Corresponding author. 1The authors are
with the MEKAR (Mechatronics Research) Labs and Automotive Control and
Mechatronics Research Center, Department of Mechanical Engineering, Istanbul
Technical University, Istanbul TR-34437, Turkey (e-mail<EMAIL_ADDRESS>guvencl@itu.edu.tr).

###### Abstract

A parameter space procedure for designing chosen parameters of a repetitive
controller to satisfy a robust performance criterion is presented. Using this
method, low order robust repetitive controllers can be designed and
implemented for plants that possibly include time delay, poles on the
imaginary axis and discontinuous weights. A design and simulation study based
on a high speed atomic force microscope position control example is utilized
to illustrate the method presented in this paper.

###### Index Terms:

Repetitive control; Parameter space method; Robust performance; Atomic force
microscope (AFM) control

##  I Introduction

Repetitive controllers are used to accurately track a periodic reference
signal or to reject a periodic disturbance with a known period by introducing
a highly frequency selective gain through a positive feedback loop which
contains a time delay element as this is a generator of periodic signals. The
delay time is equal to the known period of the repetitive reference (or
disturbance) signal. Repetitive control system is a special type of servo-
system but its basic structure is based on the Internal Model Principle of
Francis and Wonham [1]. Significant improvements in the tracking accuracy or
disturbance rejection characteristics of systems subject to periodic exogenous
signals can be achieved using repetitive control. The idea of repetitive
control was first created by Inoue et al [2] to replace conventional motion
control techniques in the control of a proton synchrotron magnet power supply.
Until recently, it has been widely utilized in many application areas
including control of hard disc drives [3], control of optical disc drives [4],
control of noncircular tuning [5], trajectory control of industrial robot arms
[6, 7], motor speed control [8], high precision rotational control [9],
control of material testing machine [10], control of cold rolling process
[11], suppression of torque vibration in motors [12], reduction of waveform
distortion in PWM inverter or UPS [13, 14, 15] and accurate position control
of piezoelectric actuators [16, 17].

The earlier papers in the literature have generally focused on the stability
analysis in both continuous time [18, 19] and discrete time systems [20]. Tsao
and Tomizuka [5] have analyzed the robust stability of repetitive control
systems applied to plants with unstructured modeling error. In order to
achieve a specified level of nominal performance, Srinivasan et al [21] have
utilized the Nevanlinna-Pick interpolation method to modify repetitive
controller design by means of optimizing a measure of stability robustness.
Peery and Özbay have modified $\mathcal{H}_{\infty}$optimal design approach
presented in [22] and then applied the extension of this methodology based on
Youla parameterization to repetitive control systems in [23]. Moon et al [4]
have developed a robust design methodology for parametric uncertainty in
interval plants under repetitive control. Similarly, Roh and Chung [24] have
created a new synthesis method based on Kharitonov’s theorem for repetitive
control systems with uncertain parameters. Weiss et al [25, 26, 27] have made
a stability and robustness analysis for MIMO repetitive control systems based
on $\mathcal{H}_{\infty}$ control theory. $\mu$ analysis has been used for
assessing stability and performance robustness of SISO continuous time
repetitive control systems by Güvenç [28]. $\mu$ synthesis has been applied to
sampled data repetitive controller design by Li and Tsao [29].

The repetitive controller design approach presented in this paper is a
continuation of the work presented in Aksun Güvenç and Güvenç [30] on
repetitive controller design based on mapping the nominal performance and
robust stability frequency domain constraints to controller parameter space
where a servo-hydraulic material testing machine application was used. This
work, in contrast, treats the robust performance constraint. Additionally, the
efficiency of methodology is illustrated by using a high speed AFM scanner
application. Moreover, the repetitive controller design approach presented in
this paper is significantly different from those of the abovementioned
references including the application of $\mathcal{H}_{\infty}$ methods. The
significant advantages of the approach here in comparison with
$\mathcal{H}_{\infty}$ methods are: (i) the ease of visualization due to the
graphical representation of the solution in the parameter space approach and
the capability and ease of doing multi-objective optimization by simply
intersecting solution regions for different objectives, (ii) the determination
of a solution region rather than one specific solution for the control system
satisfying a frequency domain constraints (this makes it easier to design non-
fragile controllers as changes in controller parameters will not violate the
chosen objectives so long as the parameters are within the solution region),
(iii) the determination of controller parameters that guarantee robust
performance, (iv) being able to treat plants with time delay and poles on the
imaginary axis, (v) not having to use rational, continuous weights in the
robust performance specifications, and (vi) obtaining fixed structure low
order repetitive controller filters that are easily implementable. There are
also some shortcomings of the proposed design method in comparison to the
methods that exist in the literature including $\mathcal{H}_{\infty}$ methods
such as: i) the method can simultaneously accommodate the design of only two
controller parameters due to its graphical display of the solution region, ii)
the method does not result in a single analytical solution and the methods
used do not look mathematically elegant as a constructive frequency-by-
frequency design approach is used.

It is difficult to apply standard robust control methods like
$\mathcal{H}_{\infty}$ control to repetitive controller design for robust
performance as the repetitive control system is infinite dimensional due to
the presence of the inherent time delay in the controller. Robust control
methods such as $\mathcal{H}_{\infty}$ optimal control have been extended to
infinite dimensional systems and applied to repetitive control (see [21] and
[23], for example). However, very high order weighting functions need to be
used in the robust controller synthesis. Consequently, the resulting
repetitive controller filters also have high order. Model order reduction
techniques are used to reduce the order of the repetitive controller filters
in an actual implementation. Some of the most powerful characteristics of the
proposed method are that the weights used in the design do not need to be
continuous functions of frequency and that plants can have time delay and/or
poles on the imaginary axis because the computations are naturally carried out
only at the frequencies of interest. Secondly, the choice of the frequency
grid used is not a problematic issue for the repetitive control design
procedure presented here as the main design frequencies are exactly known and
are the fundamental frequency of the periodic exogenous signal (reference or
disturbance) with the period $\tau_{d}$ and its harmonics. The largest
harmonic frequency considered is chosen to be close to the bandwidth of the
repetitive control system which is limited by the bandwidth of the actuator
used in the implementation. The method presented here is for SISO systems;
however, it can be used to design controllers for MIMO systems where one loop
at a time design is possible.

### Paper Organization:

The organization of the rest of the paper is as follows: Section 2 gives some
basic information on robust repetitive controller design. In Section 3, the
technique of mapping robust performance frequency domain specifications into
repetitive controller parameter space is explained in detail. Then, a
numerical example of a high speed AFM scanner position control application is
utilized in order to demonstrate the effectiveness of the proposed method in
Section 4. The paper ends with conclusions in Section 5.

##  II Repetitive Control Basics

Consider the repetitive control structure shown in Fig. 1 where $G_{n}$ is the
nominal model of the plant, $\Delta_{m}$ is the normalized unstructured
multiplicative model uncertainty, $W_{T}$ is the multiplicative uncertainty
weighting function and $\tau_{d}$ is the period of the periodic exogenous
signal. $q(s)$ and $b(s)$ are filters used for tuning the repetitive
controller. Repetitive control systems can track periodic signals very
accurately and can reject periodic disturbances very satisfactorily. This is
due to the fact that the positive feedback loop in Fig. 1 is a generator of
periodic signals with period $\tau_{d}$ for $q(s)=1$. A low pass filter with
unity d.c. gain is used for $q(s)$ for robustness of stability [18, 25].

Figure 1: Repetitive control structure.

The repetitive controller design involves the design of the two filters $q(s)$
and $b(s)$ seen in Fig. 1. In the frequency domain, the ideal low-pass filter
$q(j\omega)$ would be 1 in the frequency range of interest and zero at higher
frequencies. This is not possible and $q(j\omega)$ will have negative phase
angle which will make $q(j\omega)$ differ from 1, resulting in reduced
accuracy. So as to improve the accuracy of the repetitive controller, a small
time advance is customarily incorporated into $q(s)$ to cancel out the
negative phase of its low-pass filter part within its bandwidth. This small
time advance can easily be absorbed by the much larger time delay $\tau_{d}$
corresponding to the period of the exogenous input signal and does not
constitute an implementation problem (see Fig. 2).

The main objective of the usage of the dynamic compensator $b(s)$ is improving
the relative stability, the transition response and the steady state accuracy
in combination with the low-pass filter $q(s)$. Consider the function of
frequency given by

$\displaystyle\mathcal{R}(\omega)\triangleq\bigg{|}q(j\omega)\Big{[}1-b(j\omega)\frac{G(j\omega)}{1+G(j\omega)}\Big{]}\bigg{|}\;,$
(1)

which is called the regeneration spectrum in [19].
$\mathcal{R}(\omega)<1-\epsilon$ for all $\omega\in[0,\infty)$ and some
positive $\epsilon$ is a sufficient condition for stability [19]. Moreover,
$\mathcal{R}(\omega)$ can be utilized to obtain a good approximation of the
locus of the dominant characteristic roots of the repetitive control system
for large time delay, thus resulting in a measure of relative stability, as
well. Accordingly, the compensator $b(s)$ is designed to approximately invert
$G/(1+G)$ within the bandwidth of in an effort to minimize
$\mathcal{R}(\omega)$. The dynamic compensator $b(s)$ can be selected as only
a small time advance or time advance multiplied by a low-pass filter in order
to minimize $\mathcal{R}(\omega)$. In order to make $\mathcal{R}(\omega)<1$,
the time advance in the filter $b(s)$ is chosen to cancel out the negative
phase of $G/(1+G)$. This small time advance can easily be absorbed by the much
larger time delay $\tau_{d}$ corresponding to the period of the exogenous
input signal and does not constitute an implementation problem (see Fig. 2).

The $q(s)$ and $b(s)$ filters are thus expressed as

$\displaystyle q(s)=q_{p}(s)e^{\tau_{q}s}\;,$ (2)

and

$\displaystyle b(s)=b_{p}(s)e^{\tau_{b}s}\;.$ (3)

Figure 2: Modified repetitive control structure.

The time advances $\tau_{q}$ and $\tau_{b}$ are firstly chosen to decrease the
magnitude of $\mathcal{R}(\omega)$ given in Eq. (1). Then, the design focuses
on pairs of chosen parameters in $q_{p}(s)$ or $b_{p}(s)$ to satisfy a
frequency domain bound on the robust performance criterion. If $L$ denotes the
loop gain of a control system, its sensitivity and complementary sensitivity
transfer functions are

$\displaystyle S\triangleq\frac{1}{1+L}\;,$ (4)

and,

$\displaystyle T\triangleq\frac{L}{1+L}\;.$ (5)

The parameter space design, presented in the following section, aims at
satisfying the condition

$|W_{S}S|+|W_{T}T|<1\;,\qquad\forall\omega\in[0,\infty)\;,$ (6)

which is similar to satisfying the robust performance requirement
$\parallel|W_{S}S|+|W_{T}T|\parallel_{\infty}<1$ where $W_{S}$ and $W_{T}$ are
sensitivity and complementary sensitivity function weights.

The loop gain of the repetitive control system seen in Figs. 1 and 2 is given
by

$\displaystyle
L=G_{n}\Big{(}1+\frac{q_{p}}{1-q_{p}e^{(-\tau_{d}+\tau_{q})s}}b_{p}e^{(-\tau_{d}+\tau_{q}+\tau_{b})s}\Big{)}\;.$
(7)

The robust performance design requires

$\displaystyle|W_{S}(\omega)S(j\omega)|+|W_{T}(\omega)T(j\omega)|=\bigg{|}\frac{W_{S}(\omega)}{1+L(j\omega)}\bigg{|}+\bigg{|}\frac{W_{T}(\omega)L(j\omega)}{1+L(j\omega)}\bigg{|}<1\;,\quad\forall\omega\in[0,\infty)\;,$
(8)

or equivalently,

$\displaystyle|W_{S}(\omega)|+|W_{T}(\omega)L(j\omega)|<|1+L(j\omega)|\;,\quad\forall\omega\in[0,\infty)\;,$
(9)

to be satisfied.

##  III Mapping Robust performance Frequency Domain Specifications into
Repetitive Controller Parameter Space

In the present section, a repetitive controller design procedure based on
mapping the robust performance frequency domain performance specification
given in Eq. (9) with an equality sign into the chosen repetitive controller
parameter plane at a chosen frequency is described.

Consider the robust performance problem given in Fig. 3 illustrating Eq. (9)
with an equality sign (called the robust performance point condition). Apply
the cosine rule to the triangle with vertices at the origin, $-1$ and $L$ in
Fig. 3 to obtain

$\displaystyle(|W_{S}(\omega)|+|W_{T}(\omega)L(j\omega)|)^{2}=|L(j\omega)|^{2}+1^{2}+2|L(j\omega)|\cos{\theta_{L}}\;.$
(10)

Equation (10) is a quadratic equation in $|L(j\omega)|$ and its solutions are

$\displaystyle|L(j\omega)|=\frac{(-\cos{\theta_{L}}+|W_{S}(\omega)||W_{T}(\omega)|)\pm\sqrt{\Delta_{M}}}{1-|W_{T}(\omega)|^{2}}\;,$
(11)

where

$\displaystyle\Delta_{M}=\cos^{2}{\theta_{L}}+|W_{S}(\omega)|^{2}+|W_{T}(\omega)|^{2}-2|W_{S}(\omega)|^{2}|W_{T}(\omega)|^{2}\cos{\theta_{L}}-1\;.$
(12)

Only, positive and real solutions for $|L|$ are allowed, i.e., $\Delta_{M}\geq
0$ in Eq. (11) must be satisfied. The point condition solution procedure is
outlined below.

M1. Define the set of frequencies to be used as

$\displaystyle\Omega=\big{\\{}\omega_{1},\omega_{2},\ldots,\omega_{n};\omega_{n+1},\omega_{n+2},\ldots,\omega_{m};\omega_{m+1},\omega_{m+2},\ldots,\omega_{l}\big{\\}},\omega_{k}=\frac{2\pi
k}{\tau_{d}},k=1,2,\ldots,l$

where $\omega_{1}=2\pi/\tau_{d}$ is the frequency of the periodic exogenous
input and $\omega_{k}=2\pi k/\tau_{d}$ is the chosen bandwidth of repetitive
control (limited by the bandwidth of the actuator used). Frequencies
$\omega_{m+1}$ to $\omega_{l}$ are high frequencies where significant model
uncertainty exists ($\omega_{m+1}>10\omega_{n}$) and the intermediate
frequencies $\omega_{n+1}$ to $\omega_{m}$.

_Remark:_ It should be noted that the inherent time delay in a repetitive
control system will improve performance only at the fundamental frequency
$\omega_{1}=2\pi/\tau_{d}$ and its harmonics. Repetitive control will worsen
performance at frequencies between the fundamental frequencies and harmonics.
For this reason, repetitive control is only used in the presence of a periodic
external input (reference or disturbance) as it will result in degraded
performance for non-periodic external inputs. For that reason, the weights
($W_{S}$ and $W_{T}$) can be assumed to be zero outside the finite set
$\Omega$ given in step M1. Once the design is complete, the designer checks
the $|S|$ plot at low frequencies and the $|T|$ plot at high frequencies to
make sure that the magnitude envelopes corresponding to intermediate
frequencies (between the harmonics) are at an acceptable level. Another
approach will be to specify weights for intermediate frequencies in between
the fundamental frequency and harmonics. A major deficiency of parameter space
methods is that one needs to use a large number of frequencies in the set
$\Omega$ to make sure that the robust performance condition will not be
violated at frequencies outside of $\Omega$. In the case of repetitive
control, this problem is less severe at low frequencies where the designer is
interested mainly in reducing the sensitivity function at the fundamental
frequency and its harmonics. A large number of frequency points can be used at
higher frequencies above the bandwidth of the repetitive controlled system.

Figure 3: Illustration of the point condition for robust performance.

M2. Choose a specific frequency value $\omega=\omega_{i}\in\Omega$,
$i=1,2,\ldots,l$ from set $\Omega$ in step M1.
$|W_{S}(\omega)|$,$|W_{T}(\omega)|$ and $|G(\omega)|$ at a frequency $\omega$
are known at this point.

M3. Let $\theta_{L}\in[0,2\pi]$. Evaluate $\Delta_{M}$ by using Eq. (12) and
select the active range of $\theta_{L}$ where $\Delta_{M}\geq 0$ is satisfied.
For all values of $\theta_{L}$ in the active range:

M3a. Evaluate $|L|$ by using Eq. (11). Keep only the positive solutions.

M3b. Evaluate $L=|L|e^{j\theta_{L}}$.

M3c. Solve for the corresponding repetitive controller filters
$q_{p}(j\omega)$ and $b_{p}(j\omega)$ at the chosen frequency $\omega$ by
utilizing

$\displaystyle
q_{p}(j\omega)=\frac{L(j\omega)-G(j\omega)}{L(j\omega)-G(j\omega)[1-b(j\omega)]}e^{(\tau_{d}-\tau_{q})j\omega}\;,$
(13)

and

$\displaystyle
b_{p}(j\omega)=[L(j\omega)-G(j\omega)]\Big{[}\frac{1-q(j\omega)e^{-\tau_{d}j\omega}}{q(j\omega)e^{-\tau_{d}j\omega}}\Big{]}e^{-\tau_{b}j\omega}\;.$
(14)

M3d. Using the specific structure of $q_{p}(j\omega)$ or $b_{p}(j\omega)$,
back solve for the two chosen controller parameters within them. For example,
$q_{p}(s)$ and $b_{p}(s)$ can be chosen as a multiplication of the second
order controllers given by

$\displaystyle
q_{p}(s)=\prod_{i=1}^{n}\frac{q_{5i}s^{2}+q_{4i}s+q_{3i}}{q_{2i}s^{2}+q_{1i}s+q_{0i}},\quad\text{and}\quad~{}b_{p}(s)=\prod_{i=1}^{n}\frac{b_{5i}s^{2}+b_{4i}s+b_{3i}}{b_{2i}s^{2}+b_{1i}s+b_{0i}}\;.$
(15)

There are six tunable parameters for $n=1$ in Eq. (15) which can be used to
represent different types of controllers. These six tunable parameters are
$q_{5i}$, $q_{4i}$, $q_{3i}$, $q_{2i}$, $q_{1i}$ and $q_{0i}$ for $q_{p}(s)$.
For $n=1$, the controller structure in (15) consists of some well-known
controller types such as proportional-integral-derivative (PID), lead-lag
controller, first or second order filters as illustrated in Table 1. If the
performance of the filters which are utilized in the repetitive controller for
$n=1$ is unsatisfactory, $n$ can be increased and new higher order filters can
be synthesized. For the filter structure choice in Eq. (15), the back solution
procedure uses

$\displaystyle\begin{array}[]{rcl}\mathrm{Re}[q_{p}(j\omega)]=\frac{(q_{3i}-q_{5i}\omega^{2})(q_{0i}-q_{2i}\omega^{2})+q_{1i}q_{4i}\omega^{2}}{(q_{0i}-q_{2i}\omega^{2})^{2}+(q_{1i}\omega)^{2}}\;,\\\
\mathrm{Im}[q_{p}(j\omega)]=\frac{q_{4i}\omega(q_{0i}-q_{2i}\omega^{2})-q_{1i}\omega(q_{3i}-q_{5i}\omega^{2})}{(q_{0i}-q_{2i}\omega^{2})^{2}+(q_{1i}\omega)^{2}}\;,\end{array}$
(18)

and

$\displaystyle\begin{array}[]{rcl}\mathrm{Re}[b_{p}(j\omega)]=\frac{(b_{3i}-b_{5i}\omega^{2})(b_{0i}-b_{2i}\omega^{2})+b_{1i}b_{4i}\omega^{2}}{(b_{0i}-b_{2i}\omega^{2})^{2}+(b_{1i}\omega)^{2}}\;,\\\
\mathrm{Im}[b_{p}(j\omega)]=\frac{b_{4i}\omega(b_{0i}-b_{2i}\omega^{2})-b_{1i}\omega(b_{3i}-b_{5i}\omega^{2})}{(b_{0i}-b_{2i}\omega^{2})^{2}+(b_{1i}\omega)^{2}}\;.\end{array}$
(21)

M4. The solution in step M3 above results in a closed curve which is plotted
for solving two of the twelve parameters $q_{0i}$ to $q_{5i}$ and $b_{0i}$ to
$b_{5i}$. Plot the closed curve obtained in the chosen controller parameter
space. Either the inside (drawn with a solid boundary) or outside (drawn with
a dashed boundary) of this curve is a solution of Eq. (8) at the chosen
frequency (see the ellipses in Fig. 5, for example). The region obtained is
the point condition solution in the chosen repetitive controller parameter
plane at the frequency chosen in step M2.

M5. Go back to step M2 and repeat the procedure at a different frequency until
all frequencies in set $\Omega$ are used.

M6. Plot the intersection of all point condition solutions for all frequencies
in set $\Omega$. This is the overall solution region for the robust
performance requirement.

As the solution procedure only uses frequency response data and is numerical
in nature, plants with time delay or poles on the imaginary axis and
discontinuous weights do not pose any problems. Note that solution regions for
nominal performance $|W_{S}S|<1$ for all $\omega\in[0,\infty)$ and for robust
stability $|W_{T}T|<1$ for all $\omega\in[0,\infty)$ can easily be obtained
using the algorithm above by setting $W_{S}=0$ and $W_{T}=0$, respectively. It
is then possible to concentrate on nominal performance at low frequencies,
robust performance at intermediate frequencies and robust stability at high
frequencies, obtaining three solution regions. The overall solution region in
the controller parameter space is then determined by the intersection of all
three regions for nominal performance, robust performance and robust
stability. This procedure is easily programmable and quickly results in a
controller parameter space representation of the solution. The controller
parameter space presentation obtained offers the ease of visualization of
parameter space methods (see Fig. 5) when one accepts the shortcoming of
treating only two controller parameters at a time. Multi-objective design can
easily be formed in parameter space as it amounts simply to intersection of
individual solution regions. It is also possible to determine the final design
(or tuning point) by optimizing some other criteria, such as nominal time
domain performance within the solution region obtained. In contrast to
$\mathcal{H}_{\infty}$ optimal control synthesis, there is no relationship
between the order of repetitive control filters and the complexity of weights
in this proposed method. The main strength of this method is that low-order,
easily implementable repetitive control filters are specified from the
beginning.

TABLE I: Controller Coefficients Table Control Action | $n$ | $q_{5i}$ | $q_{4i}$ | $q_{3i}$ | $q_{2i}$ | $q_{1i}$ | $q_{0i}$
---|---|---|---|---|---|---|---
P | 1 | 0 | 0 | $K$ | 0 | 0 | 1
PD | 1 | 0 | $KT_{d}$ | $K$ | 0 | 0 | 1
PI | 1 | 0 | $K$ | $KT_{i}$ | 0 | 1 | 0
PID | 1 | $KT_{d}$ | $K$ | $KT_{i}$ | 0 | 1 | 0
Lag ($\beta>1$) | 1 | 0 | $KT$ | $K$ | 0 | $\beta T$ | 1
Lead ($0<\alpha<1$) | 1 | 0 | $KT$ | $K$ | 0 | $\alpha T$ | 1
$1^{st}$ Order Filter | 1 | 0 | 0 | $K$ | 0 | $\tau$ | 1
$2^{nd}$ Order Filter | 1 | 0 | 0 | $K\omega^{2}$ | 1 | $2\zeta\omega$ | $\omega^{2}$

It is possible that for certain data sets $|W_{S}|$, $|W_{T}|$, $G$, $\omega$;
no solutions to the solution procedure outlined above exist. Nonexistence of a
solution for a specific frequency $\omega$ could be because of nonexistence of
a positive $\Delta_{M}$ in Eq. (12) or nonexistence of a positive solution
$|L|$ in Eq. (11). Nonexistence of a solution usually results from a weight
$|W_{S}|$ or $|W_{T}|$ that is too restrictive. The solution procedure, which
is programmed in an interactive fashion, results in no solution points in this
case. Then, the user will know that his robust performance requirement at that
frequency was too restrictive and has the choice of relaxing this requirement.
Note that solutions might exist at all individual frequencies, however; their
intersection in Step M6 resulting in the overall solution region, might still
be empty. In that case, the user must change the sensitivity and complementary
sensitivity weights at the problematic frequencies.

##  IV Numerical Example

In this part of the paper, the high speed atomic force microscope (AFM)
scanner which is designed and modeled in [31] is utilized as a numerical
example to explain the methodology of the multi-objective parameter space
approach for SISO repetitive controller design. The second order and fourth
order mathematical models of this high speed AFM scanner are given in [31]. In
this example, the fourth order model is used because it includes the first
mode of the piezoelectric stack in the vertical direction. The transfer
function of the AFM scanner is given by

$\displaystyle
G(s)=\frac{K(s^{2}+2\zeta_{2}\omega_{2}s+\omega_{2}^{2})}{(s^{2}+2\zeta_{1}\omega_{1}s+\omega_{1}^{2})(s^{2}+2\zeta_{3}\omega_{3}s+\omega_{3}^{2})}\;,$
(22)

where $K=1\times 10^{12}$nm/V includes the power amplifier and sensor gain.
The system seen in (22) has two resonant frequencies and one anti-resonant
frequency. The numerical values of these frequencies are given as $f_{1}=40.9$
kHz, $f_{2}=41.6$ kHz and $f_{3}=120$ kHz and can be seen in Fig. 4. The
numerical values of the relative damping coefficients are given as
$\zeta_{1}=0.016$, $\zeta_{2}=0.016$ and $\zeta_{3}=0.17$. The dynamic
compensator $b(s)$ is chosen as a pure time advance as

$\displaystyle b(s)=b_{p}(s)e^{\tau_{b}}=e^{3\times 10^{-6}s}\;.$ (23)

The low-pass filter $q(s)$ is chosen as

$\displaystyle
q(s)=q_{p}(s)e^{\tau_{q}s}=\frac{a_{0}}{s^{2}+a_{1}s+a_{0}}e^{7.5\times
10^{-6}s}\;.$ (24)

The parameters of $q_{p}(s)$ given in Eq. (24) are chosen as
$q_{51}=q_{41}=0$, $q_{21}=1$, $q_{31}=q_{01}=a_{0}$ and $q_{11}=a_{1}$ in the
general form (15) in order to obtain unity d.c. gain. Phase advance is also
added to this low-pass filter phase cancellation. Thus, a decrease in the
steady state error is aimed. The region in the $q_{01}-q_{11}$ controller
parameter space are computed for three cases which are respectively the
nominal performance at low frequencies ($W_{T}=0$), robust performance at
intermediate frequencies and robust stability at high frequencies ($W_{S}=0$).

Figure 4: The Bode magnitude plot of high speed AFM-scanner with the mapping frequencies for the nominal performance, robust performance and robust stability. TABLE II: Desired Sensitivity Magnitude Upper Bounds at $\tau_{d}=0.0005sec$ Frequency Range | $k$ | $f=k/\tau_{d}(kHz)$ | $W_{S}$ | $W_{T}$
---|---|---|---|---
Low (NP) | 1 | 2 | 500 | 0
Low (NP) | 2 | 4 | 225 | 0
Low (NP) | 3 | 6 | 115 | 0
Low (NP) | 4 | 8 | 75 | 0
Intermediate (RP) | 40 | 80 | 3.3 | 0.001
Intermediate (RP) | 50 | 100 | 4.5 | 0.045
Intermediate (RP) | 55 | 110 | 4.5 | 0.001
Intermediate (RP) | 60 | 120 | 1.5 | 0.005
Intermediate (RP) | 70 | 140 | 1.5 | 0.01
High (RS) | 80 | 160 | 0 | 0.05
High (RS) | 90 | 180 | 0 | 0.05
High (RS) | 100 | 200 | 0 | 0.05

(a) $|W_{S}S|<1$ and $|W_{T}T|<1$, $\forall\omega\in[0,\infty)$

(b) $|W_{S}S|+|W_{T}T|<1$, $\forall\omega\in[0,\infty)$

Figure 5: Parameter space region of the low-pass filter $q_{p}(s)$.

(a) Output signal of the controlled system

(b) Error signal of the controlled system

Figure 6: Simulation results for triangular wave input at $2$ kHz.

The sensitivity constraints are specified at a set of discrete frequencies.
The periodic input command of the high speed AFM scanner has a period of
$\tau_{d}$ sec. The specific numerical values of the chosen weights used in
the computation of the controller parameters are seen in Table 2. The
frequencies corresponding to the weights in Table 2 are shown with dots in the
Bode magnitude plot of Fig. 4. The overall region calculated for nominal
performance, robust performance and stability robustness can be seen in Fig.
5b. The method in Aksun Güvenç and Güvenç [30] where the solution region is
obtained by intersecting the nominal performance and stability robustness
plots is presented in Fig. 5a for allowing easy comparison with an existing
method. Note that the sufficient stability condition $\mathcal{R}(\omega)<1$
for all $\omega\in[0,\infty)$ was also mapped and nominal stability is thus
satisfied for the two solution regions shown in Fig. 5. Also note that the
results in reference [30] are for the more conservative problem where nominal
performance and robust stability solution regions are found separately. This
conservatism is reduced by the method in the present paper where the robust
performance condition in Eq. (6) is handled directly. The method in this paper
is more general in nature and contains the method of [30] as a special case.

The mathematical model of the high-speed AFM scanner cannot be fitted very
well for frequencies above $160$ kHz. A uniform weight for the robust
stability requirement for frequencies above this value is chosen here as
$|T|<0.05$ for $f\geq 160$ kHz. This corresponding discontinuous weight
$W_{T}$ has been shown graphically with red-colored cross sign in Fig. 4. The
relative multiplicative error $|(G-G_{n})/G_{n}|$ has to be below the weight
specified in the stability robustness considerations given in Fig 4. The
intersection of the regions, which are calculated in order for the nominal
performance, the robust performance and the robust stability requirements, in
the $q_{01}-q_{11}$ controller parameter space is filled with green color. The
designation procedure is concluded by choosing a point in the controller
parameter plane given in Fig. 5. The solution within this region is chosen
arbitrarily in this example and is point is marked with a cross in Fig. 5. The
simulation result for a triangular wave input with the period $2$ kHz and
amplitude can be seen in Fig. 6. This result shows the effectiveness of the
repetitive controller in decreasing the steady state error while tracking a
periodic input signal.

##  V Conclusion

A multi-objective parameter space repetitive controller design procedure for
satisfying a robust performance objective was presented here. The main idea
was to use a simple easily implementable structure for the repetitive
controller filters and compute solution regions in the chosen controller
parameter space where frequency domain specifications on the nominal
performance at low frequencies ($W_{T}=0$), robust performance at intermediate
frequencies and robust stability at high frequencies ($W_{S}=0$) are
satisfied. The abovementioned method is well suited to the structure of a
repetitive control system with discrete frequencies of interest and the
computations were also quite fast. The proposed method is successfully applied
to the infinite dimensional nature of the repetitive control system with its
inherent time delay. The effectiveness of the proposed method was demonstrated
by carrying out a design and simulation study for high speed AFM scanner
position control problem.

## References

  * [1] B. A. Francis and W. M. Wonham, “The internal model principle for linear multivariable regulators,” _Applied Mathematics and Optimization_ , vol. 2, pp. 170–194, 1975.
  * [2] T. Inoue, S. Iwai, and M. Nakano, “High accuracy control of a proton synchrotron magnet power supply,” in _Proceedings of the $8^{th}$ IFAC World Congress_, ser. Part 3, 1981, pp. 3137–3142.
  * [3] H. Fujimoto, F. Kawakami, and S. Kondo, “Switching based repetitive control for hard disk drives: experiments on RIDR and RPTC methods,” in _Proceedings of the 2004 IEEE International Conference on Control Applications_ , vol. 1, 2004, pp. 99–104.
  * [4] J. H. Moon, M. N. Lee, and M. Chung, “Repetitive control for the track-following servo system of an optical disk drive,” _IEEE Transactions on Control Systems Technology_ , vol. 6, no. 5, pp. 663–667, 1998\.
  * [5] T. C. Tsao and M. Tomizuka, “Robust adaptive and repetitive digital tracking control and application to hydraulic servo for noncircular machining,” _ASME Journal of Dynamical Systems, Measurement and Control_ , vol. 116, no. 1, pp. 24–32, 1994.
  * [6] T. Omata, S. Hara, and M. Nakano, “Nonlinear repetitive control with application to trajectory control of manipulators,” _Journal of Robotic Systems_ , vol. 4, no. 5, pp. 631–652, 1987.
  * [7] N. Sadegh, R. Horowitz, W. W. Kao, and M. Tomizuka, “A unified approach to the design of adaptive and repetitive controllers for robotic manipulators,” _ASME Journal of Dynamical Systems, Measurement and Control_ , vol. 112, no. 4, pp. 618–629, 1990.
  * [8] F. Kobayashi, S. Hara, and H. Tanaka, “Reduction of motor speed fluctuation using repetitive control,” in _Proceedings of the $29^{th}$ IEEE Conference on Decision and Control_, 1990, pp. 1697–1702.
  * [9] R. F. Fung, J. S. Huang, and C. C. Chien, “Design and application of a continuous repetitive controller for rotating mechanisms,” _International Journal of Mechanical Sciences_ , vol. 42, no. 9, pp. 1805–1819, 2000.
  * [10] F. R. Shaw and K. Srinivasan, “Discrete-time repetitive control systems design using the regeneration spectrum,” _ASME Journal of Dynamical Systems, Measurement and Control_ , vol. 115, no. 2A, pp. 228–237, 1993.
  * [11] S. S. Garimella and K. Srinivasan, “Application of repetitive control to eccentricity compensation in rolling,” _ASME Journal of Dynamical Systems, Measurement and Control_ , vol. 118, no. 4, pp. 657–664, 1996.
  * [12] S. Tinghsu, S. Hattori, M. Ishida, and T. Hori, “Suppression control method for torque vibration of AC motor utilizing repetitive controller with fourier transform,” _IEEE Transactions on Industry Application_ , vol. 38, no. 5, pp. 1316–1325, 2002.
  * [13] C. Rech, H. Pinheriro, and H. Grundling, “Comparison of digital control techniques with repetitive integral action for low cost PWM inverters,” _IEEE Transactions on Power Electronics_ , vol. 18, no. 1, pp. 401–410, 2003\.
  * [14] K. Zhang, Y. Kang, and J. Xiong, “Direct repetitive control of SPWM inverter for UPS purpose,” _IEEE Transactions on Power Electronics_ , vol. 18, no. 3, pp. 784–792, 2003.
  * [15] K. Zhou and D. Wang, “Digital repetitive controlled three-phase PWM rectifier,” _IEEE Transactions on Power Electronics_ , vol. 18, no. 1, pp. 309–316, 2003.
  * [16] G. S. Choi, Y. A. Lim, and G. H. Choi, “Tracking position control of piezoelectric actuators for periodic reference inputs,” _Mechatronics_ , vol. 12, no. 5, pp. 669–684, 2002.
  * [17] G. H. Choi, H. O. Jong, and G. S. Choi, “Repetitive tracking control of a course-fine actuator,” in _Proceedings of IEEE/ASME International Conference Advanced Intelligent Mechatronics_ , 1999, pp. 335–340.
  * [18] S. Hara, Y. Yamamoto, T. Omata, and M. Nakano, “Repetitive control systems: A new type servo system for periodic exogenous signals,” _IEEE Transactions on Automatic Control_ , vol. 33, no. 7, pp. 657–667, 1988.
  * [19] K. Srinivasan and F. R. Shaw, “Analysis and design of repetitive control systems using the regeneration spectrum,” _ASME Journal of Dynamical Systems, Measurement and Control_ , vol. 113, no. 2, pp. 216–222, 1991.
  * [20] M. Tomizuka, T. C. Tsao, and K. K. Chew, “Analysis and synthesis of discrete-time repetitive controllers,” _ASME Journal of Dynamical Systems, Measurement and Control_ , vol. 111, no. 3, pp. 353–358, 1989.
  * [21] K. Srinivasan, H. Özbay, and I. S. Jung, “A design procedure for repetitive control systems,” in _Proceedings of ASME International Mechanical Engineer’s Congress and Exposition_ , 1995, pp. 581–593.
  * [22] H. Özbay, “${H}_{\infty}$ optimal controller design for a class of distributed parameter space,” _International Journal of Control_ , vol. 58, no. 4, pp. 739–782, 1993.
  * [23] T. Peery and H. Özbay, “${H}_{\infty}$ optimal repetitive controller design for stable plants,” _ASME Journal of Dynamical Systems, Measurement and Control_ , vol. 119, no. 3, pp. 541–547, 1997.
  * [24] C. L. Roh and M. J. Chung, “Design of repetitive control system for an uncertain plant,” _Electronics Letter_ , vol. 31, no. 22, pp. 1959–1960, 1995.
  * [25] G. Weiss and M. Häfele, “Repetitive control of MIMO systems using ${H}_{\infty}$ design,” _Automatica_ , vol. 35, no. 7, pp. 1185–1199, 1999\.
  * [26] G. Weiss, Q.-C. Zhong, T. C. Green, and J. Liang, “${H}_{\infty}$ repetitive control of DC-AC converters in microgrids,” _IEEE Transactions on Power Electronics_ , vol. 19, no. 1, pp. 219–230, 2004.
  * [27] Q.-C. Zhong, J. Liang, G. Weiss, C. Feng, and T. C. Green, “${H}_{\infty}$ control of the neutral point in four-wire three-phase DC-AC converters,” _IEEE Transactions on Industrial Electronics_ , vol. 53, no. 5, pp. 1594–1602, 2006.
  * [28] L. Güvenç, “Stability and performance robustness analysis of repetitive control systems,” _ASME Journal of Dynamical Systems, Measurement and Control_ , vol. 118, no. 3, pp. 593–597, 1996.
  * [29] J. Li and T. C. Tsao, “Robust performance repetitive control design using structured singular values,” in _Proceedings of ASME International Mechanical Engineer’s Congress and Exposition_ , 1998.
  * [30] B. A. Güvenç and L. Güvenç, “Robust repetitive controller design in parameter space,” _ASME Journal of Dynamical Systems, Measurement and Control_ , vol. 128, no. 2, pp. 406–413, 2006.
  * [31] G. Schitter, K. J. A. ström, B. E. DeMartini, P. J. Thurner, K. L. Turner, and P. K. Hansma, “Design and modeling of a high-speed AFM-scanner,” _IEEE Transactions on Control Systems Technology_ , vol. 15, no. 5, pp. 906–915, 2007.


§ DYNAMICS OF INFALLING CLOUDS

§.§ Cloud Evolution and Terminal Velocities

A falling cloud growing via accretion can be described by the following set of differential equations:

\begin{align}
\frac{\dd z}{\dd t}    &= v \label{eq:diff_eq1} \\
\frac{\dd (mv)}{\dd t} &= mg - \frac{1}{2}\rho_{\rm hot} v^2 C_0 A_{\rm cross} \label{eq:diff_eq2} \\
\frac{\dd m}{\dd t}    &= \frac{m}{t_{\rm grow}}
\label{eq:diff_eq3}
\end{align}

where $z$, $v$ and $m$ represent the distance fallen, velocity, and mass of the cloud respectively, $t_{\rm grow} \equiv m/\Dot{m}$ is the growth timescale (which we discuss in Section <ref>), $g$ is the gravitational acceleration, $C_0$ is the drag coefficient (geometry dependent; of order unity here), $\rho_{\rm hot}$ is the density of the background medium, and $A_{\rm cross}$ is the cross-sectional area which the cloud presents to the background flow. We shall see that it is important to distinguish $A_{\rm cross}$ from $A_{\rm cloud}$, the overall surface area of the cloud.  We shall also see that $t_{\rm grow}$ is roughly independent of mass growth, so that from equation (<ref>), mass growth is nearly exponential. Note that equation (<ref>) assumes steady growth and omits terms which contribute to cloud destruction. Thus, it does not apply to clouds which are losing rather than gaining mass. In this paper, we focus on scenarios where clouds survive and grow, which is the novel feature in our new model (previous works, e.g. , have looked at scenarios with significant mass loss). In Section <ref> we will quantify the criterion for cloud survival. In this work, we only consider the hydrodynamic case and leave investigation of other factors such as magnetic fields, externally driven turbulence, and cosmic rays to future work.

The terms on the right hand side in the momentum equation (equation (<ref>)) represent the gravitational and hydrodynamic drag forces. In standard models, these two terms are assumed to balance one another in steady-state, giving the hydrodynamic drag terminal velocity

\begin{align}
v_{\rm T,drag} = \sqrt{\frac{2mg}{\rho_{\rm hot} C_0 A_{\rm cross}}} \simeq \sqrt{\frac{2 \chi L g}{C_0}}
\label{eq:vtdrag}
\end{align}

for a falling cloud with volume $\sim A_{\rm cross} L$ and $\chi = \rho_{\rm cloud}/\rho_{\rm hot}$. The hydrodynamic drag time (momentum divided by the drag force) is given by $t_{\rm drag} \sim \chi L/v$. In fact, this gives the terminal velocity only if the left hand side of equation (<ref>) vanishes, $\dot{p} = m \dot{v} + \dot{m} v = 0 \Rightarrow \dot{v} =0$, which is correct only if cloud mass does not evolve so $\dot{m}=0$.  If $\Dot{m} > 0$, i.e. the cloud grows by accreting mass from the background, then from momentum conservation, since the background gas is at rest and has zero initial momentum, this will slow down the cloud. In the limit that the hydrodynamic drag term is small compared to $\Dot{m}v$,

\begin{align}
v = \frac{m}{\Dot{m}}(g - \Dot{v}).
\label{eq:nonconstantv}
\end{align}

Thus, if $\dot{v} \ll g$, there is a second terminal velocity

\begin{align}
v_{\rm T,grow} = \frac{mg}{\Dot{m}} = g t_{\rm grow}.
\label{eq:vtgrow}
\end{align}

The two terminal velocities in equations (<ref>) and (<ref>) represent regimes where the cloud acceleration under gravity is predominantly balanced by either hydrodynamic drag or the momentum transfer from background accretion respectively. We can separate them by considering the ratio $t_{\rm grow}/t_{\rm drag}$. When this ratio is large, gravity is balanced by drag. Conversely, when this ratio is small, gravity is balanced by accretion. The transition between the two is marked by where $t_{\rm grow} \sim t_{\rm drag}$.

The terminal velocity (normalized by the drag terminal velocity) as a function of the ratio of the growth time $t_{\rm grow}$ and the drag time $t_{\rm drag}$. The dashed lines show the corresponding values for the terminal velocities when assumed to be set by either drag or growth. As the ratio increases, there is a smooth transition from $v_{\rm T,grow}$ to $v_{\rm T,drag}$.

We can illustrate this by solving equations (<ref>) – (<ref>) numerically for a constant $t_{\rm grow}$. The result is shown in Fig. <ref>, where we plot the terminal velocity as a function of $t_{\rm grow}/t_{\rm drag}$. We can see that when $t_{\rm grow} \ll t_{\rm drag}$, the terminal velocity follows $v_{\rm T,grow}$, and when $t_{\rm grow} \gg t_{\rm drag}$, the terminal velocity follows $v_{\rm T,drag}$ as expected.

Which regime is more realistic? Let us first explore this for the idealized case of spherical clouds. For a spherical cloud, the growth time is given by <cit.>:

\begin{equation}
t_{\rm grow} \equiv \frac{m}{\dot m}\sim \frac{\rho_{\rm cold}r^3}{\rho_{\rm hot}A_{\rm cloud} v_{\rm mix}} \sim \chi \frac{r}{v_{\rm mix}}.
\label{eq:tgrow-approx}
\end{equation}

This seems long: if $v_{\rm mix} \sim c_{\rm s,cold}$ (a reasonable estimate; see Section <ref> of  and Sections 4.6 & 5.3.3 of ), then $t_{\rm grow} \sim \chi t_{\rm sc}$, where $t_{\rm sc}$ is the sound crossing time across the cloud. By contrast, the hydrodynamic drag time for a spherical cloud (as mentioned previously) is:

\begin{equation}
t_{\rm drag} \sim \chi \frac{r}{v},
\end{equation}

which is much shorter, since $t_{\rm drag}/t_{\rm grow} \sim v_{\rm mix}/v \sim c_{\rm s,cold}/c_{\rm s,hot} \sim \chi^{-1/2} \ll 1$, if we assume the the virial velocity to be a characteristic infall speed, $v \sim v_{\rm vir} \sim c_{\rm s,hot}$. The fact that $t_{\rm drag} \ll t_{\rm grow}$ makes physical sense. The hydrodynamic drag time is also the timescale for a cloud to sweep up its own mass in hot gas ($\rho_{\rm hot}A_{\rm cross}vt_{\rm drag} \sim \rho_{\rm hot}r^3\chi \sim m$). Even if all this mass is incorporated into the cloud, then at best $t_{\rm grow} \sim t_{\rm drag}$. In fact, only a small fraction of this gas is actually incorporated into the cloud, so that $t_{\rm grow} \gg t_{\rm drag}$. This suggests that hydrodynamic drag is the main drag force, which results in a terminal velocity given by equation (<ref>).

However, as previously mentioned, clouds in a shearing wind do not remain spherical; they develop extended cometary tails (as seen both in simulations and observations). This change in geometry – and in particular the large increase in surface area – is crucial for enabling momentum transfer via mass growth. In hydrodynamic drag, $F_{\rm drag} \sim \rho_{\rm hot} v^2 A_{\rm cross}$, the area $A_{\rm cross}\approx \pi r^2$ is the cross-sectional area the cloud presents to the wind. Thus, $F_{\rm drag}$ remains roughly constant during cloud evolution. By contrast, in $\dot{m} \sim \rho_{\rm hot} A_{\rm cloud} v$, the area $A_{\rm cloud}$ is the surface area of the cloud available for mixing. In a cometary structure, this is dominated by the sides of the cylinder, so that $A_{\rm cloud} \sim 2 \pi r L$, where $L$ is the length of the tail. Thus, $\dot{m} \propto A_{\rm cloud} \propto L$ increases as a cloud develops a cometary tail. It is this increase in $\dot{m}$, and thus the effective momentum transfer rate $F_{\rm grow} \sim \dot{p}_{\rm grow} \sim \dot{m} v$, compared to a constant $F_{\rm drag}$, which causes mass growth to dominate momentum transfer: $t_{\rm grow} \sim \chi r/v_{\rm mix}$ is roughly constant, while $t_{\rm drag} \sim m v/F_{\rm drag} \sim \rho_{\rm cloud} A_{\rm cross} L v/(\rho_{\rm hot} A_{\rm cross} v^{2}) \sim \chi L/v$ increases as the mass of the cloud increases. In particular,

\begin{equation}
\frac{t_{\rm grow}}{t_{\rm drag}} \sim \frac{r}{L} \frac{v}{v_{\rm mix}} \sim \frac{r}{L} {\chi^{1/2}}.
\label{eq:tgrow_tdrag}
\end{equation}

In cloud crushing simulations, the tail grows during the process of entrainment to a length $L \sim vt_{\rm drag} \sim \chi r$ during the `tail formation' phase <cit.>, so that $t_{\rm grow}/t_{\rm drag} \sim \chi^{-1/2} \ll 1$. The continuous shear for infalling clouds can lead to even more extended tails since the cloud does not entrain, so $t_{\rm grow}/t_{\rm drag} \ll 1$ is easily satisfied[Shorter entrainment times than $t_{\rm drag}$ have been observed in cloud crushing simulations <cit.>.].

Finally, it is important to realize that there is a third timescale in the problem, the free-fall time $t_{\rm ff}\sim v_{\rm vir}/g$. This sets the evolutionary lifetime available to clouds, before they fall to the halo center. Clouds will not grow significantly (and reach the terminal velocity $v_{\rm T,grow}$ given by equation (<ref>)), unless $t_{\rm grow} < t_{\rm ff}$.
Indeed, $t_{\rm grow} < t_{\rm ff}$ is required for a subvirial terminal velocity. We can show this by recalling that $F_{\rm grav} \sim m g \sim m v_{\rm vir}/t_{\rm ff}$, while the drag force from mass growth is $F_{\rm grow} \sim \dot{m}v \sim m v/t_{\rm grow}$. At the terminal velocity $v_T$, we have $F_{\rm grav} \sim F_{\rm grow}$, so that:

\begin{equation}
f_{\rm sub-vir} \equiv \frac{v_{\rm T,grow}}{v_{\rm vir}} \sim \frac{v_{\rm T}}{c_{\rm s,hot}} \sim \frac{t_{\rm grow}}{t_{\rm ff}}.
\label{eqn:alpha}
\end{equation}

This is useful because $f_{\rm sub-vir}$ – infall velocities, normalized to the virial velocity – can be measured observationally. Indeed, $f_{\rm sub-vir} < 1$, sub-virial infall velocities, is commonly observed in LRGs <cit.> and galaxy clusters <cit.>, much lower than predicted terminal velocities from hydrodynamic drag models <cit.>. Our models can explain these puzzling observations, as we describe in Section <ref>. It also allows for testable predictions. Since $f_{\rm sub-vir}$ is measured and $t_{\rm ff}$ is known from the density profile, we can predict $t_{\rm grow} \approx f_{\rm sub-vir} t_{\rm ff}$ from kinematic observations, assuming that clouds have reached terminal velocity. This can be compared with predictions for $t_{\rm grow}$ from equations (<ref>) and (<ref>), given measured or inferred cloud and background hot gas properties. Lastly, the mass growth that a cloud experiences is $m/m_0 \sim \exp(t_{\rm ff}/t_{\rm grow}) \sim {\rm exp} (f_{\rm sub-vir}^{-1})$. Thus, a measurement of sub-virial velocities directly constrains the degree to which mixing and cooling enhances cool gas infall to the central galaxy. Significantly sub-virial infall implies that cold clouds grow considerably before reaching the halo center. These analytical estimates can be compared to measurements of the mass infall rate <cit.>. In Section <ref>, we will also show the rather remarkable result that in an isothermal atmosphere with constant gravity, $f_{\rm sub-vir}$ is fixed by geometry, specifically the scaling between cloud mass and area ($\alpha$ in equation (<ref>)), independent of all other properties of the system. For our infalling clouds, we find $f_{\rm sub-vir} \approx 0.6$.

§.§ Cloud Growth

Previous models of infalling clouds have considered the interplay between gravity and hydrodynamic drag forces, assuming a fixed cloud mass <cit.>. However, a fixed cloud mass is unrealistic due to various processes that trigger mixing with the hot background gas or shred the cloud. Mass evolution therefore cannot remain static; clouds should either be destroyed ($\dot{m} < 0$), or grow ($\dot{m} > 0$) over time.

In the absence of cooling, clouds moving relative to a background medium are destroyed by hydrodynamic instabilities on the cloud crushing timescale <cit.>

\begin{align}
t_{\rm cc} \sim \sqrt{\chi} \frac{r}{v},
\end{align}

where $\chi$ is the ratio of the cloud density to the background density, $r$ is the cloud radius, and $v$ is the magnitude of the relative velocity between the cloud and the background. This cloud crushing timescale reflects the destruction of the cloud via internal shocks induced inside the cloud due to its velocity with respect to the medium it is moving through (assuming that this velocity is supersonic with respect to the sound speed within the cloud), and is roughly the same timescale on which surface instabilities such as the Kelvin-Helmholtz and Rayleigh-Taylor instabilities grow to the cloud scale <cit.>. This destructive fate can however be counteracted by mass growth due to cooling. In wind tunnel simulations of `cloud crushing', <cit.> found that in order for cold gas to survive, cooling needs to be strong enough to satisfy the criterion

\begin{align}
t_{\rm cool,mix} < t_{\rm cc},
\label{eq:t_cool_mix}
\end{align}

where $t_{\rm cool,mix}$ is the cooling time of the mixed gas, defined as $T_{\rm mix} \sim \sqrt{T_{\rm cloud} T_{\rm hot}}$ (in the spirit of , see also  for an alternative derivation). That is, if the cooling time of the mixed gas is shorter than the initial cloud crushing time, then cold gas survives and is eventually entrained in the hot background wind.

However, infalling clouds have an important aspect that differentiates them from clouds in a wind – gravity. Clouds encountering a hot wind gradually entrain in the wind, so that shear eventually drops to zero if the cloud manages to survive until entrainment. The cloud thus encounters destructive forces for a limited period of time. By contrast, clouds in a gravitational field will always keep falling and shearing against the background gas. Thus, the survival criterion is different, and more stringent; we discuss this in Section <ref>.

Assuming cloud survival, let us quantify the timescale on which clouds grow. We first derive some scaling relations, before deriving numerical expressions. For now, we ignore fudge factors (due to geometry, etc.) which can be up to an order of magnitude. As in equation (<ref>), the mass growth rate of a cloud can be written as

\begin{align}
\Dot{m} \sim \rho_{\rm hot} A_{\rm cloud} v_{\rm mix},
\label{eq:mdot}
\end{align}

where $\rho_{\rm hot}$ is the density of the hot background medium, $A_{\rm cloud}$ is the effective surface area of the cloud[The effective surface area corresponds to the (smoothed) enveloping area of the cloud and not the (non-convergent) surface area of the cold gas. See <cit.> for further discussion of this distinction.] and $v_{\rm mix}$ is the the velocity corresponding to the mass flux from the hot background onto this surface. As above, if we write $m \sim \rho_{\rm cold} A_{\rm cloud} r$, this gives:

\begin{equation}
t_{\rm grow} \sim \chi \frac{r}{v_{\rm mix}}.
\label{eq:tgrow-simple}
\end{equation}

Plane parallel simulations of mixing layers <cit.> show:

\begin{equation}
v_{\rm mix} \sim u^{\prime 3/4} \left( \frac{r}{t_{\rm cool}} \right)^{1/4}
\sim v_{\rm shear}^{3/5} v_0^{3/20} \left( \frac{r}{t_{\rm cool}} \right)^{1/4}
\label{eq:vmix-simple}
\end{equation}

where $t_{\rm cool}$ is the cooling time in cold gas (the minimum cooling time in the mixing layer, a convention we adopt henceforth) and $u^{\prime}$ is the peak turbulent velocity in the mixing layer (usually in intermediate temperature gas). Note that while the first step in equation (<ref>), i.e., $v_{\rm mix}(u')$, is generally valid, we have used the scaling $u^{\prime} \propto v_{\rm shear}^{4/5}$ for relating $u'$ to the parameters of the setup. This scaling was found numerically in <cit.> for plane-parallel mixing layers and we have written it here as $u^{\prime} \sim v_{\rm shear}^{4/5} v_0^{1/5}$ to preserve dimensionality ($v_0$ simply encodes normalization). If we set $v_{\rm shear} \sim v_{\rm T,grow} \sim g t_{\rm grow}$, this yields:

\begin{equation}
t_{\rm grow} \sim \chi^{5/8} \frac{r^{15/32} t_{\rm cool}^{5/32}}{g^{3/8} v_0^{3/32}}; \ \ {v_{\rm T,grow} < c_{\rm s,hot}}.
\label{eq:tgrow_nonumbers}
\end{equation}

While the above scalings focus on the subsonic and transonic cases, large enough clouds can reach velocities exceeding the sound speed of the hot gas. In such a case, the turbulent mixing velocity saturates and stops scaling with the cloud velocity <cit.>, changing the above scalings. In this case, from equations (<ref>) and (<ref>), we obtain:

\begin{equation}
t_{\rm grow} \propto \frac{\chi  r^{3/4} }{{\rm c_{\rm s,hot}^{3/5}} t_{\rm cool}^{1/4}}.
\label{eq:tgrow-sat-simple}
\end{equation}

We now give numerical expressions, which are calibrated to simulations.
For cooling dominated regimes (defined below), <cit.> found that $v_{\rm mix}$ in turbulent mixing layers follows

\begin{equation}
v_{\rm mix} \approx
9.5 {\mathrm{~km}\,\mathrm{s}^{-1}}
\left(\frac{u^{\prime}}{50 \, {\rm km \, s^{-1}}}\right)^{3/4}
\left(\frac{L_{\rm turb}}{100\,{\rm pc}}\right)^{1/4}
\left(\frac{t_{\rm cool}}{0.03\,{\rm Myr}}\right)^{-1/4},
\label{eq:vmix}
\end{equation}

where $L_{\rm turb}$ is the outer scale of the turbulence.
Note that equation (<ref>) only applies in the `fast cooling' (${\rm Da}_{\rm mix} \equiv L_{\rm turb}/(u^{\prime} t_{\rm cool,mix}) > 1$, where ${\rm Da}_{\rm mix}$ is the Damköhler number; ) regime, where the cooling time
is much smaller than the turbulent mixing time $L_{\rm turb}/u'$. As we will discuss below, however, this is always true for surviving clouds.

<cit.> note that $u'$ is geometry dependent, but find for shearing layers that

\begin{equation}
u' \approx 50{\mathrm{~km}\,\mathrm{s}^{-1}}
\mathcal{M}^{4/5}
\left(\frac{c_{\rm s,hot}}{150 {\mathrm{~km}\,\mathrm{s}^{-1}}} \right)^{4/5}
\left(\frac{t_{\rm cool}}{0.03\,{\rm Myr}}\right)^{-0.1} ,
\label{eq:uprime}
\end{equation}

for $\chi\gtrsim 100$ and $\mathcal{M}\equiv v_{\rm shear} / c_{\rm s,hot}$. From equation (<ref>), we can approximate $v_{\rm mix} \sim c_{\rm s,cold}$ for quick estimates. While <cit.> only considered mixing layers with subsonic to transonic velocity shears, <cit.> found that beyond $\mathcal{M} = 1$, $u'$ in the mixing region stops scaling with $\mathcal{M}$ and saturates. We include this in our model by setting $\mathcal{M} \rightarrow \min(1,\mathcal{M})$. We find good evidence for this in our simulations.

Equations (<ref>) and (<ref>) assume fully developed turbulence. When a cloud falls from rest however, there is a transient period when turbulence is developing. We hence set a time dependent weight factor $w_{\rm kh}(t)$ to account for the initial onset of turbulence. Turbulence develops over the timescale for the development of the Kelvin-Helmholtz instability; on the scale of the cloud $t_{\rm kh} = f_{\rm kh} t_{\rm cc}$ where $f_{\rm kh}$ is some constant of proportionality <cit.>. We use the simplest ansatz that

\begin{equation}
{v}_{\rm mix} \rightarrow w_{\rm kh}(t) v_{\rm mix};  \ \
w_{\rm kh}(t) = \min\left(1, \frac{t}{f_{\rm kh} t_{\rm cc}}\right),
\label{eq:w_kh}
\end{equation}

which amounts to $v_{\rm in}$ growing linearly with time over the instability growth time, until fully developed and capped at unity. We will justify this ansatz in our simulations. Since $t_{\rm cc}$ is changing over time, we note that $t/t_{\rm cc} \propto vt \sim ~z$, where $z$ is distance the cloud has fallen. We find in our simulations that $f_{\rm kh} \sim 5$ for a constant background and $\sim 1$ for a stratified background. In a more realistic setting with less idealized initial conditions, this time-dependent weight factor might not be necessary as the initial mixing can be already seeded from the outflowing section (assuming $v<v_{\rm esc}$), extrinsic turbulence, or cooling induced pulsations <cit.>.

What is an appropriate scaling relation for the effective cloud surface area $A_{\rm cloud}$? In cloud crushing simulations, areal growth follows two phases <cit.>. In the `tail-formation' phase, surface area growth is dominated by the formation of a cometary tail, with $A_{\rm cloud} \propto L \propto m$, where $L$ is the length of the tail. The stretching of the cloud means that the area to mass ratio $A_{\rm cloud}/m \approx $ constant, rather than $A_{\rm cloud}/m \propto m^{-1/3}$, as for fixed geometry. Once the tail grows to a length $L \sim \chi r$ (the hydrodynamic drag length), the cloud becomes entrained in the wind from efficient momentum transfer, and due to lack of shear the tail no longer grows. The cloud surface area thereafter scales roughly as $A_{\rm cloud} \propto (m/\rho_{\rm cloud})^{2/3}$, as one would expect for a monolithic cloud.

However, our falling clouds do not get entrained - rather the opposite in fact, as they start at rest and accelerate until reaching some terminal velocity. This means they start `entrained' and then begin to shear against background gas. They never leave the `tail-formation' phase, since there is a constant velocity difference between the cloud and background medium. The cloud sees a continuous headwind which drives turbulence, mixing, and lengthening. Instead of $A_{\rm cloud} \propto m/\rho_{\rm cloud}$ or $A_{\rm cloud} \propto (m/\rho_{\rm cloud})^{2/3}$, we assume that $A_{\rm cloud} \propto (m/\rho_{\rm cloud})^{\alpha}$, where $\alpha$ is a growth scaling exponent between 2/3 and 1. Physically, this is because both mass growth onto the surface of the cloud and a lengthening tail are concurrent processes. We will demonstrate that this is a good assumption for the mass growth of the falling clouds in our simulations, where we find $\alpha \approx$ 5/6. The cloud surface area is thus

\begin{align}
A_{\rm cloud} \approx A_{\rm cloud,0} \left( \frac{m}{m_0} \frac{\rho_{\rm cloud,0}}{\rho_{\rm cloud}} \right)^{\alpha},
\label{eq:acloud}
\end{align}

where $A_{\rm cloud,0}$, $\rho_{\rm cloud,0}$ and $m_0$ are the initial cloud surface area, density and mass respectively. Note that since $\dot{m} \propto m^{\alpha}$ where $\alpha = 5/6$ is close to 1, the growth is close to exponential[Similar scalings $\dot{m} \propto A \propto m^{\alpha}$, where $\alpha \approx 0.8$, are seen in simulations of cloud growth when clouds are embedded in a turbulent medium <cit.>. This super-Euclidean scaling can be understood as the outcome of the fractal nature of the mixing surface, where area $A \propto m^{D/3}$, where $D$ is the fractal dimension <cit.>. In their mixing layer simulations, <cit.> measure $D \approx 2.5$, which gives $\alpha \approx D/3 = 5/6$, consistent with the above.]. The cloud density $\rho_{\rm cloud}$ changes because the ambient pressure increases as the cloud falls in a stratified medium, compressing the cloud.

Using equations (<ref>), (<ref>), (<ref>) and (<ref>), we can write the growth time $t_{\rm grow} \sim m/\Dot{m}$ as

\begin{align}
t_{\rm grow} = \frac{t_{\rm grow,0}}{w_{\rm kh}(t)}
\left( \frac{c_{\rm s,150}}{v} \right)^{3/5}
\left( \frac{t_{\rm cool}}{t_{\rm cool,0}} \right)^{1/4}
\left( \frac{m}{m_0} \frac{\rho_{\rm hot,0}}{\rho_{\rm hot}} \right)^{1-\alpha},
\label{eq:tgrow_scaling}
\end{align}

where $c_{\rm s,150} = 150\,{\rm km\,s}^{-1}$ is the sound speed of gas at $10^6$ K and the initial growth time $t_{\rm grow,0}$ is given by

\begin{align}
t_{\rm grow,0} \approx 35\,{\rm Myr }\,
\left( \frac{f_{\rm A}}{0.23} \right)
\left( \frac{\chi}{100} \right)
\left( \frac{r}{r_{100}} \right)
\left( \frac{L_{\rm turb}}{L_{100}} \right)^{-1/4}
\left( \frac{t_{\rm cool,0}}{0.03\,{\rm Myr}} \right)^{1/4}.
\label{eq:tgrow_scaling2}
\end{align}

where $r_{100} = L_{100} = 100$ pc and $r$ is the initial cloud size. We will assume generally that $L_{\rm turb} \sim r$ (since the hydrodynamic instabilities which drive turbulence and mixing have an outer scale set by cloud size). We have included an unknown normalization factor $f_{\rm A}$ to account for uncertainties arising from geometrical differences between the single mixing layers in <cit.> and our cloud setup here, the use of the initial size of the sphere as a characteristic scale (see discussion at the end of Section <ref>), and any other simplifying assumptions we might have made. We find in our simulations that $f_{\rm A} \sim 0.23$.
We can simplify equations (<ref>) and (<ref>) by ignoring the weak mass and hot gas density dependence, and setting $L_{\rm turb} \sim r$, to obtain:

\begin{align}
t_{\rm grow}     = \frac{35\,{\rm Myr }}{w_{\rm kh}(t)} \,
\left(\frac{f_{\rm A}}{0.23}\right) \left( \frac{c_{\rm s,150}}{v} \right)^{3/5}
\left( \frac{\chi}{100} \right)
\left( \frac{r}{r_{100}} \right)^{3/4}
\left( \frac{t_{\rm cool}}{0.03\,{\rm Myr}}\right)^{1/4}.
\label{eq:tgrow_terminal}
\end{align}

Equations (<ref>) and (<ref>) should be used when evaluating $t_{\rm grow}$ if the velocity $v(t)$ varies with time (i.e., when solving equations (<ref>) – (<ref>). However, a key quantity is the growth time at the terminal velocity $v=g t_{\rm grow}$, which we shall see determines whether the cloud can survive (Section <ref>). Inserting $v = g t_{\rm grow}$ into equation (<ref>), setting $w_{\rm kh}(t)=1$, and using $f_A=0.23$, we obtain the numerical version of equation (<ref>):

\begin{align}
t_{\rm grow}     = 40\,{\rm Myr }\,
\left( \frac{g}{g_{\rm fid}}\right)^{-3/8}
\left( \frac{\chi}{100} \right)^{5/8}
\left( \frac{r}{r_{100}} \right)^{15/32}
\left( \frac{t_{\rm cool}}{0.03\,{\rm Myr}} \right)^{5/32},
\label{eq:tgrow_vgt}
\end{align}

where $g_{\rm fid} = 10^{-8}\,{\rm cm\,s}^{-2}$. On the other hand, for supersonic speeds, as we have discussed, the turbulent mixing velocity saturates and stops scaling with the cloud velocity <cit.>. Setting $v \sim gt_{\rm grow}$ to $ v \sim c_{\rm s,hot}$ instead in equation (<ref>), we find the numerical version of equation (<ref>):

\begin{align}
t_{\rm grow}     = 35\,{\rm Myr }\,
\left( \frac{c_{\rm s,hot}}{c_{s,150}} \right)^{-3/5}
\left( \frac{\chi}{100} \right)
\left( \frac{r}{r_{100}} \right)^{3/4}
\left( \frac{t_{\rm cool}}{0.03\,{\rm Myr}} \right)^{1/4}.
\label{eq:tgrow_sat}
\end{align}

§.§ Cloud Survival

The model we have presented only accounts for mass growth of the cloud and does not include processes that result in mass loss. In addition, the initial onset and development of turbulence is only very crudely incorporated. The absence of these refinements mean that we should expect differences between model predictions and simulations, certainly for clouds that are losing mass, and at early times even for clouds that do survive and grow. We leave the inclusion and refinement of these components for future work, as we find that the model as presented works well for surviving clouds. Since the key assumption of our model is that the cloud is growing, we now discuss when this is a valid assumption.

As we previously discussed, clouds placed in a wind tunnel encountering a hot wind can survive if $t_{\rm cool, mix} < t_{\rm cc}$ (equation (<ref>)).
Physically, $t_{\rm cool,mix}$ can be understood as the time it takes gas to cool in the downstream tail region of the cloud. Even if the initial pristine cloud material does not survive, if mixed gas can cool and survive, then the cold gas mass will increase. Since this mixed gas in the tail is cooling from the background, it is much more entrained in the wind than the initial cloud and hence able to survive – once the cold gas is entrained, it is no longer subject to destruction by shear.

The `usual' survival criterion $t_{\rm cool,mix} < t_{\rm cc}$ above is certainly a necessary condition for survival. If no gas can cool before the cloud is completely disrupted, the cloud cannot survive. However, this criterion is not a sufficient one. This is because the physical process associated with $t_{\rm cc}$ is not simply surface evaporation. If this were so, then the above criterion would indeed be sufficient as any mixing would lead to a net increase in cloud mass. Instead, the entire cloud is disrupted <cit.>. Hence, as we shall see, it is not enough that mixed gas can cool faster than the cloud crushing time.

Compared to a wind tunnel setup, the considerations for an infalling cloud are different. Since the cloud's velocity increases instead, and there is no entrainment, $t_{\rm cc}$ decreases over time.
The only way for cold gas to survive is if it is produced at a rate faster than it is destroyed:

\begin{align}
t_{\rm grow} < f_{\rm S} t_{\rm cc},
\label{eq:survive}
\end{align}

where $f_{\rm S}$ is some constant[Although we find that a constant factor is sufficient for our purposes, this coefficient has been found to vary in supersonic flows. For example, <cit.> found that in the cloud crushing setup with a supersonic wind, $f_{\rm S}$ scales as $\sqrt{1+\mathcal{M}_{\rm hot}}$ where $\mathcal{M}_{\rm hot}$ is the Mach number of the hot medium <cit.>. However, we mostly probe the subsonic to transonic regime.] factor of order unity, which we shall calibrate in simulations. It encodes the fact cloud destruction takes place over several cloud crushing times <cit.>. In evaluating $t_{\rm cc} \sim \chi^{1/2} r/v$, the cloud radius is evaluated at its initial value. As in wind tunnel experiments, this turns out to be a very good approximation, since the cloud grows mostly in the streamwise direction. If the velocity is evaluated at the terminal velocity $v_{\rm T} \sim g t_{\rm grow}$, then equation (<ref>) is equivalent to:

\begin{equation}
\frac{g t_{\rm grow}^{2}}{\chi^{1/2} r} < f_{\rm S}.
\label{eq:survive-terminal}
\end{equation}

As we have seen, there are two regimes for $t_{\rm grow}$, subsonic and supersonic infall. The criterion for subsonic infall is $t_{\rm grow} < t_{\rm ff}$ (equation (<ref>)). Using equation (<ref>), and assuming $t_{\rm ff} \sim c_{\rm s,hot}/g$, this can be rewritten as $r< r_{\rm sonic}$, where

\begin{equation}
r_{\rm sonic} \sim 150\,{\rm pc }\,
\left(\frac{t_{\rm cool}}{0.03\,{\rm Myr}}\right)^{-1/3}
\left( \frac{g}{g_{\rm fid}} \right)^{-4/3}
\left( \frac{\chi}{100} \right)^{-4/3}
\left( \frac{c_{\rm s,hot}}{c_{\rm s,150}} \right)^{32/15}.
\label{eq:rsonic}
\end{equation}

Thus, clouds must be smaller than some critical radius to fall at sub-virial velocities. In this regime, ($v_{\rm T} \lsim c_{\rm s,hot}$), $t_{\rm grow}$ is given by equation (<ref>), and the survival criterion, equation (<ref>), becomes:

\begin{eqnarray}
t_{\rm cool} < 5 \times 10^{-3} \, {\rm Myr}
\left(\frac{f_{\rm S}} {2}\right)^{16/5}
\left(\frac{r}{r_{100}}\right)^{1/5}
\left( \frac{g}{g_{\rm fid}} \right)^{-4/5}
\left( \frac{\chi}{100} \right)^{-12/5} .
\label{eq:survival_criterion}
\end{eqnarray}

Note that equation (<ref>) is almost independent of cloud size. Indeed, $t_{\rm grow}/t_{\rm cc} \propto g t_{\rm grow}^2/r \propto r^{-1/16}$, i.e., a very weak scaling. We shall verify this in Section <ref>.

Is it possible for clouds to survive in the supersonic regime ($r> r_{\rm sonic}$)? This requires $t_{\rm ff} < t_{\rm grow} < f_S t_{\rm cc}$. This in turn requires that clouds be smaller than some critical size $r_{\rm SS}$, since $t_{\rm grow}/t_{\rm cc} \propto gt_{\rm grow}^{2}/r \propto r^{1/2}$ in the supersonic regime (using $t_{\rm grow} \propto r^{3/4}$ from equation (<ref>)). Thus, supersonic infall and survival requires:

\begin{equation}
r_{\rm sonic} < r < r_{\rm SS} ,
\label{eq:survive_supersonic}
\end{equation}

where $r_{\rm SS}$ is given by:

\begin{equation}
r_{\rm SS} = 100\,{\rm pc }\,
\left(\frac{t_{\rm cool}}{0.03\,{\rm Myr}}\right)^{-1}
\left( \frac{g}{g_{\rm fid}} \right)^{-2}
\left( \frac{\chi} {100}\right)^{-2}
\left( \frac{c_{\rm s,hot}}{c_{\rm s,150}} \right)^{12/5}.
\label{eq:rSS}
\end{equation}

Note that equation (<ref>) can only be fulfilled if $r_{\rm sonic}/r_{\rm SS} < 1$, where:

\begin{equation}
\frac{r_{\rm sonic}}{r_{\rm SS}} \sim 1.5
\left(\frac{t_{\rm cool}}{0.03\,{\rm Myr}}\right)^{2/3}
\left( \frac{g}{g_{\rm fid}} \right)^{2/3}
\left( \frac{\chi} {100}\right)^{2/3}
\left( \frac{c_{\rm s,hot}}{c_{\rm s,150}} \right)^{-4/15}.
\end{equation}

Cloud survival for subsonic and supersonic infall for different cloud          sizes and cooling times. Survival is mostly sensitive to the latter.

Figure <ref> shows the survival criteria above (equations (<ref>) and (<ref>)) for $g = g_{\rm fid}$, $\chi=100$, $c_{\rm s,hot}=c_{\rm s,150}$, and $f_{\rm S}=2$. It is clear that survival is mostly independent of cloud size and depends instead on the cooling time.

In practice, the subsonic case is of most interest. There, clouds must satisfy $t_{\rm grow} < {\rm min}(t_{\rm ff}, t_{\rm cc})$, which translates into a maximum allowed cloud size (equation (<ref>)) and a maximum allowed cooling time in cold gas (equation (<ref>)). The latter criterion is quite stringent. Since the dependence on size in equation (<ref>) is weak, under isobaric conditions $t_{\rm cool} \propto 1/P$, we can translate equation (<ref>) into a critical pressure. For $f_{\rm S} = 2$, and ignoring the size dependence, equation (<ref>) is equivalent to:

\begin{align}
P > 3000 \,{\rm k}_B\,{\rm K}\,{\rm cm}^{-3}
\left(\frac{g}{g_{\rm fid}}\right)^{4/5} \left(\frac{\chi}{100}\right)^{12/5} ,
\label{eq:safe_zone}
\end{align}

where the RHS is the critical pressure $P_{\rm crit}$ above which a falling cloud can survive. We can also write equation (<ref>) in terms of the cooling time of the hot gas $t_{\rm cool,hot} \sim \chi^2 t_{\rm cool} [\Lambda(T_{\rm cold})/\Lambda(T_{\rm hot})]$ and the free fall time $t_{\rm ff} \sim c_{\rm s,hot}/g$ to obtain:

\begin{align}
\frac{t_{\rm cool,hot}}{t_{\rm ff}} \lesssim 1 \left( \frac{\chi}{100} \right)^{-2/5} \left( \frac{\Lambda(T_{\rm cold})}{\Lambda(T_{\rm hot})} \right),
\label{eq:tcool_over_tff}
\end{align}

where we have ignored the weak dependence on $g$, $t_{\rm cool}/t_{\rm ff} \propto g^{1/5}$. This is similar to the criterion ($t_{\rm cool,hot}/t_{\rm ff} < 1$, where $t_{\rm cool}$ is evaluated at one scale height) for precipitation out of a thermally unstable background medium in a plane parallel atmosphere[It is somewhat more stringent than the requirement for thermal instability in spherical systems ($t_{\rm cool,hot}/t_{\rm ff} < 10$; ), where the gravitational acceleration $g$ and hence $t_{\rm ff}$ varies as a function of radius. However, it has been shown that there is no geometrical difference in cold gas condensation in plane parallel and spherical geometries; the apparent difference arises from definitional differences in where $t_{\rm cool}/t_{\rm ff} $ is evaluated and cold gas is located <cit.>.] <cit.>. Since all our analytics and simulations are in the framework of plane parallel systems, the numerical factor in equation (<ref>) will likely change in spherical systems. Equation (<ref>) has the very interesting implication that clouds which condense via thermal instability are able to survive subsequent infall, as long as they are below the critical size given by equation (<ref>). Note that the physics of stratified thermal instability which leads to the $t_{\rm cool,hot}/t_{\rm ff} < 1$ criterion – overstable gravity waves driven by cooling – is quite different from what we have discussed here, so it is non-trivial (perhaps coincidental) that both thermal instability and falling cloud survival have similar criteria.