text:
\section{Introduction}

The concept of input-to-state stability (ISS), which unified Lyapunov and input-output approaches, plays a foundational role in nonlinear control theory \cite{Son08}. It is central for robust stabilization of nonlinear systems \cite{FrK08,KKK95}, design of nonlinear observers \cite{ArK01}, analysis of large-scale networks \cite{JTP94,DRW07} etc.

The success of ISS methodology for ordinary differential equations and time-delay systems, see e.g. \cite{Tee98, PeJ06, PKJ08, Krs08d}, as well as importance of robust control for distributed parameter systems, motivated a rapid development of the ISS theory for abstract infinite-dimensional systems \cite{JLR08, KaJ11, DaM13b, MiI15b, MiW17b, Mir16} and more specifically of partial differential equations (PDEs) \cite{MaP11, PrM12, MiI15b, AVP16, CDP17, PiO17, TPT17} during the last decade.
Important results achieved include characterizations of ISS and local ISS for broad classes of nonlinear infinite-dimensional systems \cite{MiW17b, Mir16}, nonlinear small-gain theorems for interconnections of $n\in\N$ infinite-dimensional systems \cite{KaJ11, MiI15b}, applications of ISS Lyapunov theory to analysis and control of various classes of PDE systems \cite{PrM12, AVP16, TPT17, CDP17, PiO17} etc.
However, most of these papers are devoted to PDEs with distributed inputs.
\textit{In this work we study input-to-state stability (ISS) of nonlinear parabolic partial differential equations (PDEs) with boundary disturbances on multidimensional spatial domains.} This question naturally arises in such fundamental problems of PDE control as robust boundary stabilization of PDE systems, design of robust boundary observers, stability analysis of cascades of parabolic and ordinary differential equations (ODEs) etc.

It is well known, that PDEs with boundary disturbances can be viewed as evolution equations in Banach spaces with unbounded input (disturbance) operators. This makes the analysis of such systems much more involved than stability analysis of PDEs with distributed disturbances (which are described by bounded input operators), even in the linear case.

At the same time, ISS of linear parabolic systems w.r.t. boundary disturbances has been studied in several recent papers using different methodologies \cite{AWP12, KaK16b, KaK17a, JNP18, ZhZ17}.
In \cite{AWP12} the authors have shown that parabolic systems with boundary disturbances are ISS w.r.t. $C^1$ norm of disturbances (known as $D^1$-ISS, see \cite[p. 190]{Son08}) provided the system is uniformly asymptotically stable if disturbances are set to zero.
In \cite{ZhZ17} linear and nonlinear parabolic systems with Neumann or Robin boundary inputs over one-dimensional spatial domain have been studied by means of Lyapunov methods combined with novel Poincare-like inequalities. As a result several criteria for ISS of such systems w.r.t. $C$-norm of boundary inputs have been obtained.

In \cite{KaK16b, KaK17a, KaK17b} linear parabolic PDEs with Sturm-Liouville operators over 1-dimensional spatial domain have been treated by using two different methods: (i) the spectral decomposition of the solution, and (ii) the approximation of the solution by means of a  finite-difference scheme.
This made possible to avoid differentiation of boundary disturbances, and to obtain ISS of classical solutions w.r.t. $L^\infty$ norm of disturbances, as well as in weighted $L^2$ and $L^1$ norms.
An advantage of these methods is that this strategy can be applied also to other types of linear evolution PDEs. At the same time, for multidimensional spatial domains, the computations can become quite complicated.

In \cite{JNP18, JSZ17} ISS of linear boundary control systems has been approached by means of methods of semigroup and admissibility theory. In particular, using the theory of Orlicz spaces and $H_\infty$ calculus, deep relations between ISS and integral input-to-state stability (iISS) have been achieved. This method can be applied to a large class of linear control systems (in particular, the results on equivalence between ISS and integral ISS have been shown in \cite{JSZ17} not only for parabolic systems, but also for the systems governed by analytic semigroups). At the same time, the lack of proper generalizations of admissibility to general nonlinear systems makes this method hard (if possible) to apply to nonlinear distributed parameter systems.

\textit{In this paper, we develop a novel method for investigation of parabolic PDEs with Dirichlet boundary disturbances.} In contrast to previous results, we do not restrict ourselves to linear equations over 1-dimensional spatial domains. Our results are valid for a class of nonlinear equations over multidimensional bounded domains with a sufficiently smooth boundary. Our method is based on the concept of monotone control systems introduced in \cite{AnS03} and inspired by the theory of monotone dynamical systems pioneered by M. Hirsch in 1980-s in a series of papers, beginning with \cite{Hir82}. For the introduction to this theory, a reader may consult \cite{Smi95}. Monotone control systems appear frequently in applications, see e.g. in chemistry and biology \cite{LAS07, Smi95}.
An early effort to use monotonicity methods to study of ISS of infinite-dimensional systems has been made in \cite{DaM10} to show ISS of a particular class of parabolic systems with distributed inputs and Neumann boundary conditions.

Our strategy is as follows.

First, we show in Section~\ref{sec:ISS_monotone_control_systems} that \textit{a monotone control system (and in particular a nonlinear parabolic PDE system with boundary disturbances) is ISS if and only if it is ISS over a much smaller class of inputs (e.g. constant inputs)}. This is achieved by proving that, for any given disturbance, there is a larger constant disturbance which leads to the larger deviation from the origin, and hence, in a certain sense, the constant disturbances are the "worst case" ones.

Next, in Section~\ref{ISS_nonlinear_parabolic_PDEs}, using maximum and comparison principles for nonlinear parabolic operators, we show that \textit{under certain regularity assumptions nonlinear parabolic equations with Dirichlet boundary disturbances are monotone control systems}.
In turn, this nice property helps us to show in Section~\ref{ISS_nonlinear_parabolic_PDEs}
that \textit{ISS of the original nonlinear parabolic PDE with constant \textit{boundary disturbances} is equivalent to ISS of a closely related nonlinear parabolic PDE with constant \textit{distributed disturbances}}.
The latter problem has been studied extensively in the last years, and a number of powerful results are available for this class of systems.
In particular, in \cite{MaP11}, constructions of strict Lyapunov functions for certain nonlinear parabolic systems have been provided. In \cite{MiI15b, DaM13}, ISS and integral ISS small-gain theorems for nonlinear parabolic systems interconnected via spatial domain have been proved which give powerful tools to study stability of large-scale parabolic systems on the basis of knowledge of stability of its components.
Results in this paper make it possible to use this machinery for analysis of ISS of nonlinear parabolic PDEs with boundary inputs.

Finally, in Section~\ref{sec:Lin_Parabolic_Systems}, we apply above results to linear parabolic problems, which are of specific interest. The obtained ISS criteria can be applied to the problem of ISS stabilization of linear parabolic systems by means of a PDE backstepping method (see for instance \cite{KrS08, SmK10}) in the presence of actuator disturbances, see \cite{MKK17} for more details.

Next, we introduce some notation used throughout these notes.
By $\R_+$ we denote the set of nonnegative real numbers. For $z \in \R^n$ the Euclidean norm of $z$ is denoted by $|z|$.
For any open set $G\subset \R^n$ we denote by $\partial G$ the boundary of $G$; by $\overline{G}$ the closure of $G$ and by $\mu(G)$ the Lebesgue measure of $G$.
Also for such $G$ and any $p\in[1,+\infty)$ we denote
by $L^p(G)$ the space of Lebesgue measurable functions $y$ with $\|y\|_p=\Big(\int_G |y(z)|^p dz\Big)^{1/p}$,
by $L^\infty(G)$ the space of Lebesgue measurable functions $y$ with $\|y\|_\infty=\esssup_{z\in G}|y(z)|$,
and by $H^k(G)$ the set of $g\in L^2(G)$, which possess weak derivatives up to the $k$-th order, all of which belong to $L^2(G)$.
$C^k(G)$ consists of $k$ times continuously differentiable functions defined on $G$ and $C_0^1(G)$ consists of a functions from $C^1(G)$ which have a compact support.
$H^{1}_0(G)$ is a closure of $C_0^1(G)$ in the norm of $H^1(G)$.
If $I \subset \R_+$, then $C^{1,2}(I\times G)$ means the set of functions mapping $I\times G$ to $\R$, which are continuously differentiable w.r.t. the first argument and possess continuous second derivatives w.r.t. the second arguments.

Also we will use the following classes of comparison functions.

\begin{equation*}
\begin{array}{ll}
{\K} &:= \left\{\gamma:\R_+\rightarrow\R_+\left|\ \gamma\mbox{ is continuous, strictly} \right. \right. \\
&\phantom{aaaaaaaaaaaaaaaaaaa}\left. \mbox{ increasing and } \gamma(0)=0 \right\}, \\
{\K_{\infty}}&:=\left\{\gamma\in\K\left|\ \gamma\mbox{ is unbounded}\right.\right\},\\
{\LL}&:=\left\{\gamma:\R_+\rightarrow\R_+\left|\ \gamma\mbox{ is continuous and strictly}\right.\right.\\
&\phantom{aaaaaaaaaaaaaaaa} \text{decreasing with } \lim\limits_{t\rightarrow\infty}\gamma(t)=0\},\\
{\KL} &:= \left\{\beta:\R_+\times\R_+\rightarrow\R_+\left|\ \beta \mbox{ is continuous,}\right.\right.\\
&\phantom{aaaaaa}\left.\beta(\cdot,t)\in{\K},\ \beta(r,\cdot)\in {\LL},\ \forall t\geq 0,\ \forall r >0\right\}. \\
\end{array}
\end{equation*}

\textit{Due to the page limits we omit most of the proofs of the results obtained in the article. Please refer to the full version of the paper for the proofs and further applications of the main results \cite{MKK17}.}

\section{Monotonicity of control systems}
\label{sec:Framework_Monotonicity}

We start with a definition of a control system.
\begin{definition}
\label{Steurungssystem}
Consider a triple $\Sigma=(X,\Uc,\phi)$, consisting of
\begin{enumerate}[(i)]
 \item A normed linear space $(X,\|\cdot\|_X)$, called the {state space}, endowed with the norm $\|\cdot\|_X$.
 \item A set of  input values $U$, which is a nonempty subset of a certain normed linear space.

 \item A normed linear space of inputs $\Uc \subset \{f:\R_+ \to U\}$ endowed with the norm $\|\cdot\|_{\Uc}$.
   We assume that $\Uc$ satisfies \textit{the axiom of shift invariance}, which states that for all $u \in \Uc$ and all $\tau\geq0$ the time
shift $u(\cdot + \tau)$ is in $\Uc$.

\item A family of nonempty sets $\{\Uc(x) \subset \Uc: x \in X\}$, where $\Uc(x)$ is a set of admissible inputs for the state $x$.

\item A transition map $\phi:\R_+ \times X \times \Uc \to X$, defined for any $x\in X$ and any $u\in \Uc(x)$ on a certain subset of $\R_+$.
\end{enumerate}
The triple $\Sigma$ is called a (forward-complete) control system, if the following properties hold:
\begin{sysnum}
 \item \textit{Forward-completeness}: for every $x\in X$, $u\in\Uc(x)$ and for all $t \geq 0$ the value
$\phi(t,x,u) \in X$ is well-defined.
 \item\label{axiom:Identity} \textit{The identity property}: for every $(x,u) \in X \times \Uc(x)$
          it holds that $\phi(0,x,u)=x$.
 \item \label{axiom:Causality}\textit{Causality}: for every $(t,x,u) \in \R_+ \times X \times
          \Uc(x)$, for every $\tilde{u} \in \Uc(x)$, such that $u(s) =
          \tilde{u}(s)$, $s \in [0,t]$ it holds that $\phi(t,x,u) = \phi(t,x,\tilde{u})$.

  \item \label{axiom:Cocycle} \textit{The cocycle property}: for all $t,h \geq 0$, for all
                  $x \in X$, $u \in \Uc(x)$ we have $u(t+\cdot) \in \Uc(\phi(t,x,u))$ and
\[
\phi(h,\phi(t,x,u),u(t+\cdot))=\phi(t+h,x,u).
\]
\end{sysnum}
\end{definition}
In the above definition $\phi(t,x,u)$ denotes the state of a system at the moment $t \in
\R_+$ corresponding to the initial condition $x \in X$ and the input $u \in \Uc(x)$.
A pair $(x,u)\in X\times\Uc(x)$ is referred to as an \textit{admissible pair}.

\begin{remark}\em
\label{rem:Compatibility_Conditions}
For wide classes of systems, in particular for ordinary differential equations, one can assume that $\Uc(x)=\Uc$ for all $x\in X$, that is every input is admissible for any state. On the other hand, the classical solutions of PDEs with Dirichlet boundary inputs have to satisfy compatibility conditions (see Section~\ref{ISS_nonlinear_parabolic_PDEs}), and hence for such systems $\Uc(x)\neq \Uc$ for any $x\in X$.
Another class of systems for which one cannot expect that $\Uc(x)=\Uc$ for all $x\in X$ are differential-algebraic equations
(DAEs), see e.g. \cite{KuM06, KrT15}.
\end{remark}

\begin{definition}
A subset $K \subset X$ of a normed linear space $X$ is called a \textit{positive cone} if $K \cap (-K)=\{0\}$
and for all $a \in \R_+$ and all $x,y \in K$ it follows that $ax \in K$; $x+y \in K$.
\end{definition}

\begin{definition}
A normed linear space $X$ together with a cone $K \subset X$ is called an \textit{ordered normed linear space} (see \cite{Kra64}), which we denote $(X,K)$ with an order $\leq$ given by $x \leq y\ \Iff \  y-x \in K$. Analogously $x \geq y\ \Iff \  x-y \in K$.
\end{definition}

\begin{definition}
\index{control system!ordered}
We call a control system $\Sigma=(X,\Uc,\phi)$ ordered, if $X$ and $\Uc$ are ordered normed linear spaces.
\end{definition}

An important for applications subclass of control systems are monotone control systems:
\begin{definition}
An ordered control system $\Sigma=(X,\Uc,\phi)$ is called monotone, provided for all $t \geq 0$, all $x_1, x_2 \in X$ with $x_1 \leq x_2$ and all $u_1 \in\Uc(x_1), u_2 \in \Uc(x_2)$ with $u_1 \leq u_2$ it holds that
$\phi(t,x_1,u_1)  \leq \phi(t,x_2,u_2)$.
\end{definition}

To treat situations when the monotonicity w.r.t. initial states is not available, the following definition is useful:
\begin{definition}
An ordered control system $\Sigma=(X,\Uc,\phi)$ is called monotone w.r.t. inputs, provided for all $t \geq 0$, all $x \in X$ and all $u_1,u_2 \in \Uc(x)$ with $u_1 \leq u_2$ it holds that $ \phi(t,x,u_1)  \leq \phi(t,x,u_2)$.
\end{definition}

\section{Input-to-state stability of monotone control systems}
\label{sec:ISS_monotone_control_systems}

Next we introduce the notion of input-to-state stability, which will be central in this paper.
\begin{definition}
\label{def:ISS_Uc}
Let $\Sigma=\left(X,\mathcal{U},\phi \right)$ be a control system. Let $\Uc_{c}$ be a subset of $\mathcal{U}$.
System $\Sigma$ is called \textit{input-to-state stable (ISS) with respect to inputs in $\Uc_{c}$} if there exist functions $\beta \in \KL$, $\gamma \in \K$ such that for every $x \in X$ for which $\Uc(x ) \bigcap \Uc_{c}$ is non-empty, the following estimate holds for all $u \in \Uc(x ) \bigcap \Uc_{c}$ and $t\ge 0$:
\begin{eqnarray}
\| \phi(t,x,u) \|_X \leq \beta(\| x \|_X,t) + \gamma( \|u\|_{\Uc}).
\label{eq:ISS_estimate}
\end{eqnarray}
If $\Sigma$ is ISS w.r.t. inputs from $\Uc$, then $\Sigma$ is called \textit{input-to-state stable (ISS)}.
\end{definition}

For applications the following notion, which is stronger than ISS is of importance:
\begin{definition}
\label{def:exp-ISS_Uc}
Let $\Sigma=\left(X,\mathcal{U},\phi \right)$ be a control system. Let $\Uc_{c}$ be a subset of $\mathcal{U}$.
System $\Sigma$ is called \textit{exponentially input-to-state stable (exp-ISS) with respect to inputs in $\Uc_{c}$} if there exist constants $M,a$ and $\gamma\in\Kinf$ such that for every $x \in X$ for which $\Uc(x ) \bigcap \Uc_{c}$ is non-empty, the following estimate holds for all $u \in \Uc(x ) \bigcap \Uc_{c}$ and $t\ge 0$:
\begin{eqnarray}
\| \phi(t,x,u) \|_X \leq Me^{-at}\| x \|_X + \gamma(\|u\|_{\Uc}).
\label{eq:exp-ISS_estimate}
\end{eqnarray}
If $\Sigma$ is exp-ISS w.r.t. inputs from $\Uc$, then $\Sigma$ is called \textit{exponentially input-to-state stable}.
If in addition $\gamma$ can be chosen to be linear, then $\Sigma$ is exp-ISS with a linear gain function.
\end{definition}

We are also interested in the stability properties of control systems in absence of inputs.
\begin{definition}
\label{def:0UGAS_Uc}
A control system $\Sigma=\left(X,\mathcal{U},\phi \right)$ is {\it globally asymptotically
stable at zero uniformly with respect to the state} (0-UGAS), if there
exists a $ \beta \in \KL$, such that for all $x \in X$: $0 \in \Uc(x)$ and for all $
t\geq 0$ it holds that
\begin{equation}
\label{UniStabAbschaetzung}
\left\| \phi(t,x,0) \right\|_{X} \leq  \beta(\left\| x \right\|_{X},t) .
\end{equation}
\end{definition}

It may be hard to verify the ISS estimate \eqref{eq:ISS_estimate} for all admissible pairs of states and inputs.
Therefore a natural question appears: to find a smaller set of admissible states and inputs, so that validity of the ISS estimate \eqref{eq:ISS_estimate} implies ISS of the system for all admissible pairs (possibly with larger $\beta$ and $\gamma$).
For example, for wide classes of systems it is enough to check ISS estimates on the properly chosen dense subsets of the space of admissible pairs (a so-called "density argument", see e.g. \cite[Lemma 2.2.3]{Mir12}). As Propositions~\ref{prop:ISS_Monotone_Systems}, \ref{prop:ISS_Monotone_states_inputs_Systems} show, for monotone control systems such sets can be much more sparse.

We start with the case, when all $(x,u) \in X \times \Uc$ are admissible pairs.
\begin{proposition}
\label{prop:ISS_Monotone_Systems}
Let $(X,\Uc,\phi)$ be a control system that is monotone with respect to inputs  with $\Uc(x)=\Uc$ for all $x\in X$. Furthermore, let $\Uc_{c}$ be a subset of $\mathcal{U}$ and let the following two conditions hold:
\begin{enumerate}
 \item[(i)] There exists $\rho \in\Kinf$ so that for any $x_-,x,x_+ \in X$ satisfying $x_-\leq x \leq x_+$ it holds that
\begin{eqnarray*}
\|x\|_X \leq \rho(\|x_-\|_X + \|x_+\|_X).
\end{eqnarray*}
 \item[(ii)] There exists $\eta \in \Kinf$ so that for any $u \in \Uc$ there are $u_-,u_+ \in \Uc_{c}$, satisfying $u_-\leq u \leq u_+$ and $\|u_-\|_{\Uc} \leq \eta(\|u\|_{\Uc})$ and $\|u_+\|_{\Uc} \leq \eta(\|u\|_{\Uc})$.
\end{enumerate}
Then $(X,\Uc,\phi)$ is ISS if and only if $(X,\Uc,\phi)$ is ISS w.r.t. inputs in $\Uc_{c}$.

Moreover, if $\rho$ is linear, then $\Sigma$ is exp-ISS if and only if $\Sigma$ is exp-ISS w.r.t. inputs in $\Uc_{c}$.
If additionally $\Sigma$ is exp-ISS w.r.t. inputs in $\Uc_{c}$ with a linear gain function and $\eta$ is linear, then $\Sigma$ is exp-ISS
with a linear gain function.
\end{proposition}

Here, as for most results in this paper we omit the proofs due to the page limit constraints. \textit{We refer to \cite{MKK17} for the proofs.}
\begin{example}
\label{rem:Results_in_ODE_context}
Consider ODEs of the form
\begin{eqnarray}
\dot{x}=f(x,u),
\label{eq:ODEsys}
\end{eqnarray}
with $X=\R^n$ with an order induced by the cone $\R^n_+$, $\Uc:=L^\infty(\R_+,\R^m)$ with the order induced by the cone $L^\infty(\R_+,\R^m_+)$ and $\Uc(x)=\Uc$ for all $x\in\R^n$.
Under assumptions that $f$ is Lipschitz continuous w.r.t. the first argument uniformly w.r.t. the second one and that \eqref{eq:ODEsys} is forward complete,
\eqref{eq:ODEsys} defines a control system $(X,\Uc,\phi)$, where $\phi(t,x_0,u)$ is a state of \eqref{eq:ODEsys} at the time $t$ corresponding to $x(0)=x_0$ and input $u$.

Assume that \eqref{eq:ODEsys} is monotone w.r.t. inputs with such $X$ and $\Uc$ and consider $\Uc_{c}:=\{u\in\Uc:\exists k\in\R^m:\ u(t)=k \text{ for a.e. } t\in\R_+\}$. It is easy to verify that assumptions of Proposition~\ref{prop:ISS_Monotone_Systems} are fulfilled, and hence \eqref{eq:ODEsys}
is ISS iff it is ISS w.r.t. the inputs with constant in time controls.

This may simplify analysis of ISS of monotone ODE systems since input-to-state stable ODEs with constant inputs have some specific properties,
see e.g. \cite[pp. 205--206]{Son08}.
\end{example}

As we argued in Remark~\ref{rem:Compatibility_Conditions}, for many systems $\Uc(x) \neq \Uc$ for some $x\in X$. For such systems Proposition~\ref{prop:ISS_Monotone_Systems} is inapplicable, due to the fact that existence of inputs $u_-,u_+ \in \Uc_{c}$ for a given initial condition $x$ and for an input $u$, satisfying assumption (ii) of Proposition~\ref{prop:ISS_Monotone_Systems} does not guarantee that the pairs $(x,u_-)$ and $(x,u_+)$ are admissible, which is needed for the proof of Proposition~\ref{prop:ISS_Monotone_Systems}.
However, if $(X,\Uc,\phi)$ is in addition monotone w.r.t. states, the following holds:
\begin{proposition}
\label{prop:ISS_Monotone_states_inputs_Systems}
Let $\Sigma=\left(X,\mathcal{U},\phi \right)$ be a control system that is monotone with respect to states and inputs.
Assume that $\Uc_{c}\subset\mathcal{U}$ and let the following two conditions hold:
\begin{itemize}
   \item[(i)]There exists $\rho \in \Kinf $ so that for any $x_{-} ,x,x_{+} \in X$ satisfying $x_{-} \le x\le x_{+} $ it holds that
\begin{equation} \label{GrindEQ__19_}
\left\| x\right\| _{X} \le \rho \left(\left\| x_{-} \right\| _{X} +\left\| x_{+} \right\| _{X} \right).
\end{equation}
\item[(ii)] There exist $\eta ,\xi \in \Kinf $ so that for every $x\in X$, $u\in \Uc(x)$ and for every $\varepsilon >0$ there exist $x_{-} ,x_{+} \in X$ and $u_{-} \in \Uc(x_{-} ) \bigcap \Uc_{c}$, $u_{+} \in \Uc(x_{+}) \bigcap \Uc_{c}$, satisfying $x_{-} \le x\le x_{+} $, $u_{-} \le u\le u_{+} $, so that the estimates
\begin{equation} \label{GrindEQ__20_}
\max \left(\left\| u_{-} \right\| _{\mathcal{U}} ,\left\| u_{+} \right\| _{\mathcal{U}} \right)\le \eta \left(\left\| u\right\| _{\mathcal{U}} +\varepsilon \right),
\end{equation}
\begin{equation} \label{GrindEQ__21_}
\max \left(\left\| x_{-} \right\| _{X} ,\left\| x_{+} \right\| _{X} \right)\le \xi \left(\left\| x\right\| _{X} +\left\| u\right\| _{\mathcal{U}} +\varepsilon \right).
\end{equation}
\end{itemize}
hold. Then $\Sigma$ is ISS $\Iff$ $\Sigma$ is ISS w.r.t. inputs in $\Uc_{c}$.

Moreover, if $\rho$ and $\xi$ are linear, then $\Sigma$ is exp-ISS if and only if $\Sigma$ is exp-ISS w.r.t. inputs in $\Uc_{c}$.
If additionally $\Sigma$ is exp-ISS w.r.t. inputs in $\Uc_{c}$ with a linear gain function and $\eta$ is linear, then $\Sigma$ is exp-ISS with a linear gain function.
\end{proposition}

\section{Input-to-state stability of nonlinear parabolic equations}
\label{ISS_nonlinear_parabolic_PDEs}

In this section we apply results from Section~\ref{sec:ISS_monotone_control_systems}
 to nonlinear parabolic equations with boundary inputs.
Our strategy is as follows: first we formulate in Proposition~\ref{prop:New_Comparison_Principle} the comparison principle for nonlinear parabolic operators. Afterwards we use this principle in Theorem~\ref{thm:Monotonicity_Parabolic_Systems} to show that under certain technical conditions the nonlinear parabolic systems with boundary disturbances give rise to monotone control systems. Finally, in Theorems~\ref{thm:ISS_Parabolic_Systems_Characterization}, \ref{thm:exp-ISS_Parabolic_Systems_Characterization} we show our main result that \textit{a nonlinear parabolic system with boundary disturbances is ISS if and only if a related system with distributed disturbances is ISS for constant inputs}. The precise statements of these results follow.

Let $G\subset \mathbb{R} ^{n} $ be an open bounded region, let $T>0$ be a constant and denote $D:=(0,T)\times G$.
Let $CL(D)$ denote the class of functions $x\in C^{0} \left(\overline{D}\right)\bigcap C^{1,2} (D)$.
Denote for each $t\in[0,T]$ and each $x\in CL(D)$ a function $x[t]:\R_+ \to C^0(\overline{G})$ by $x[t]:=x(t,\cdot)$.

Consider the operator $L$ defined for 	 $x\in CL(D)$ by
\begin{eqnarray}
\label{GrindEQ__1_}
(Lx)(t,z)&:=&\frac{\partial \, x}{\partial \, t} (t,z)-\sum _{i,j=1}^{n}a_{i,j} (z)\frac{\partial ^{2} \, x}{\partial \, z_{i} \, \partial \, z_{j} } (t,z) \nonumber\\
&& \qquad -f\big(z,x(t,z),\nabla x(t,z)\big),
\end{eqnarray}
where $(t,z)\in D$, $a_{i,j} \in C^{0} (\overline{G})$ for $i,j=1,...,n$ and $f:G\times \mathbb{R} \times \mathbb{R} ^{n} \to \mathbb{R} $ is a continuous function. The operator $L$ is called uniformly parabolic, if there exists a constant $K>0$ so that for all $\xi =(\xi _{1} ,...,\xi _{n} )\in \mathbb{R} ^{n} $ it holds that
\begin{equation} \label{GrindEQ__2_}
\sum _{i,j=1}^{n}a_{i,j} (z)\xi _{i} \xi _{j}  \ge K\left|\xi \right|^{2}  \quad \mbox{for all} \ z\in \overline{G}.
\end{equation}

We need the following proposition, based upon a classical comparison principle from \cite[Theorem 16, p. 52]{Fri83}.
\begin{proposition}
\label{prop:New_Comparison_Principle}
Let $L$ be uniformly parabolic. Assume that for every bounded set $W\subset \mathbb{R} $ there is a constant $k>0$ such that for every $w_1,w_2\in W$, $z\in G$, $\xi \in \mathbb{R} ^{n}$ with $w_1>w_2$ it holds that
\begin{equation}
\label{GrindEQ__3_}
f(z,w_1,\xi )-f(z,w_2,\xi ) < k(w_1-w_2).
\end{equation}
Let $x,y\in CL(D)$ be so that
\begin{eqnarray}
(Ly)(t,z)&\ge& 0\ge (Lx)(t,z), \quad \mbox{for all} \ (t,z)\in D,\label{GrindEQ__4_}\\
y(0,z)&\ge& x(0,z), \quad \mbox{for all} \ z\in G, \label{GrindEQ__5_}\\
y(t,z)&\ge& x(t,z), \quad \mbox{for all} \ (t,z)\in [0,T]\times \partial G. \label{GrindEQ__6_}
\end{eqnarray}
Then $y(t,z)\ge x(t,z)$ for all $(t,z)\in \overline{D}$.
\end{proposition}

Since our intention is to analyze forward complete systems, we introduce some more notation.
Let $CL$ denote the class of functions $x\in C^{0} \left(\mathbb{R} _{+} \times \overline{G}\right)\bigcap C^{1,2} ((0,+\infty )\times G)$.

Now we apply the established results to analyze the initial boundary value problem:
\begin{eqnarray}
(Lx)(t,z)&=&0, \quad \mbox{for all }(t,z)\in (0,+\infty )\times G, \label{GrindEQ__22_}\\
x(0,z)&=&x_{0} (z), \quad \mbox{for all }z\in G, \label{GrindEQ__23_}\\
x(t,z)&=&u(t,z), \quad \mbox{for all }(t,z)\in \mathbb{R} _{+} \times \partial G,  \label{GrindEQ__24_}
\end{eqnarray}
where $x_{0} \in C^{0} (\overline{G})$, $L$ is the uniformly parabolic operator defined by \eqref{GrindEQ__1_} with $f\in C^{0} (\overline{G}\times \mathbb{R} \times \mathbb{R} ^{n} )$.

In this section we assume that the space of input values is $U=C^{0} (\partial G)$, endowed with the standard sup-norm
and that
$\mathcal{U}:=\left\{\, u\in C^{0} (\mathbb{R} _{+} ;U):\, \, u \mbox{ is bounded}\right\}$, endowed with
\begin{itemize}
  \item the partial order $\le $ for which $u\le v$ iff $u(t,z)\le v(t,z)$ for all $(t,z)\in \mathbb{R} _{+} \times \partial G$,
 \item the norm
$\left\| u\right\| _{\mathcal{U}} =\mathop{\sup }\limits_{t\ge 0} \left\| u[t]\right\| _{U} =\mathop{\sup }\limits_{z\in \partial G,\ t\ge 0} \left|u(t,z)\right|$.
\end{itemize}
In the sequel we will need also a subspace of $\Uc$ consisting of constant in time and space inputs:
\begin{eqnarray}
\Uc_{c}:=\{u \in \Uc:\exists k \in\R:\ u(t,z)=k \ \forall (t,z)\in \mathbb{R} _{+} \times \partial G\}.
\label{eq:Uc_const_def}
\end{eqnarray}

Define for $x\in C^0(G)$ the standard $L^p(G)$-norm as $\left\|x \right\| _{p}:=\Big(\int_G |x(z)|^p dz\Big)^{1/p}$.
The Euclidean distance between $w\in G$ and $W \subset \overline{G}$ is denoted by $\rho(w,W):=\inf_{s\in W}|w-s|$.

The following assumptions will be instrumental:
\begin{itemize}
   \item[(H1)] There exists a linear space $X\subseteq C^{0} (\overline{G})$, containing the functions $\{x\in C^{0} (\overline{G}):\exists k\in\R \mbox{ s.t. } x(\cdot)=k\}$, such that for each $x_{0} \in X$ which is constant on $\partial G$, there exists a set of inputs $\Uc(x_{0} )\subseteq \left\{\, v\in \Uc\, :\, v(0,z)=x_{0} (z)\, \mbox{for}\, \, z\in \partial G\, \right\}$, which contains constant in time and space inputs
\[
\left\{\, v\in \Uc\, :\, v(t,z)=x_{0} (z)\, \mbox{for}\, \, z\in \partial G\, \mbox{and}\, \, t\ge 0\, \right\}
\]
with the following property: for every $x_{0} \in X$, $u\in \Uc(x_{0} )$ there exists a solution $x\in CL$ of the initial boundary
value problem \eqref{GrindEQ__22_}, \eqref{GrindEQ__23_}, \eqref{GrindEQ__24_} for which $x[t]\in X$ for all $t\ge 0$.

 \item[(H2)] Assume that for every bounded set $W\subset \mathbb{R} $ there exists $k>0$ such that for every $w_1,w_2\in W$, $(t,z)\in D$, $\xi \in \mathbb{R} ^{n} $ with $w_1>w_2$ inequality \eqref{GrindEQ__3_} holds. Moreover, the function $\overline{G}\times \mathbb{R} \times \mathbb{R} ^{n} \ni (z,w,\xi )\mapsto f(z,w,\xi )\in \mathbb{R} $ is continuously differentiable w.r.t. $w\in \mathbb{R} $ and $\xi \in \mathbb{R} ^{n} $.

 \item[(H3)] For every $\delta >0$, $a\in \mathbb{R} $, $x \in X$ there exists a continuous function $k:\overline{G}\to [0,1]$ with $k(z)=1$ for $z\in \partial G$, $k(z)=0$ for all $z\in G$ with $\rho(z,\partial G)\ge \delta $ and such that $f\in X$ where $f(z)=(1-k(z))x(z)+a\, k(z)$.
\end{itemize}

\begin{remark}\em
Note that (H3) is a condition on the geometry of the boundary of $G$ which is automatically satisfied when $G$ is an open interval in $\R$.
\end{remark}

We equip $X$ in (H1) with the partial order $\le $ for which $x\le y$ iff $x(z)\le y(z)$ for all $z\in \overline{G}$.
For existence theorems, which can be used to verify (H1), we refer to \cite[Chapters 3, 7]{Fri83}.

The next result assures monotonicity of the initial boundary value problem \eqref{GrindEQ__22_}, \eqref{GrindEQ__23_}, \eqref{GrindEQ__24_}.
\begin{theorem}
\label{thm:Monotonicity_Parabolic_Systems}
Suppose that assumptions (H1), (H2) hold and let $p\in[1,+\infty]$. Let us endow the linear space $X$ in (H1) with the standard $L^{p} (G)$-norm, which we denote by $\left\| x\right\| _{p}$.
Then:
\begin{itemize}
   \item[(i)] Initial boundary value problem \eqref{GrindEQ__22_}, \eqref{GrindEQ__23_}, \eqref{GrindEQ__24_} gives rise to
 the monotone control system $\Sigma=\left(X,\mathcal{U},\phi \right)$, where $\phi$ is the solution map of \eqref{GrindEQ__22_}, \eqref{GrindEQ__23_}, \eqref{GrindEQ__24_}.
 \item[(ii)] If additionally (H3) holds, then conditions (i) and (ii) of Proposition~\ref{prop:ISS_Monotone_states_inputs_Systems} hold with $\rho$, $\eta$, $\zeta$ being linear functions and $\Uc_{c}$ given by \eqref{eq:Uc_const_def}.
\end{itemize}
\end{theorem}

We continue to assume that the axioms (H1) and (H2) hold and that $X$ is as in (H1).

Consider now the following equations:
\begin{align}
\frac{\partial y}{\partial t}(t,z) - \sum_{i,j=1}^n  & a_{ij}(z)\frac{\partial^2 y}{\partial z_i \partial z_j}(t,z) \nonumber\\
&-  f\Big(z,y(t,z)+v(t,z), \nabla y(t,z)\Big) =0,
\label{eq:Nonlinear_Parabolic_Equation_transformed}
\end{align}
where $t >0$, $z \in G$, together with homogeneous Dirichlet boundary conditions
\begin{eqnarray}
y(t,z) = 0, \quad z \in \partial G,\ t \geq 0.
\label{eq:Dirichlet_transformed}
\end{eqnarray}
The state space of \eqref{eq:Nonlinear_Parabolic_Equation_transformed} is a linear space
\[
Y:=\{y\in X: y(z)=0,\ z\in\partial{G}\}
\]
and the input $v$ belongs to the space of constant in time and space inputs
\[
\Vc:=\{v\in C^0(\R_+\times \overline{G},\R): \exists k\in\R: v(t,z)=k,\ (t,z)\in\R_+\times \overline{G}\}.
\]
Denote the solutions of \eqref{eq:Nonlinear_Parabolic_Equation_transformed}, \eqref{eq:Dirichlet_transformed}, corresponding to the initial condition $y \in Y$ and input $v$ by $\phi_y(t,y,v)$.
It is easy to see that for any $x\in X$ and any $v \in \Vc$ for which $v|_{\partial G} \in \Uc(x)$  it holds that
\begin{eqnarray}
\phi_y(t,y,v)=\phi(t,x,v|_{\partial G}) -v,\quad \mbox{where}\ y=x-v,
\label{eq:phi_phyy_relation}
\end{eqnarray}
and $y \in Y$ since $v\in X$, $X$ is a linear space and $x(z)=v(z)$ for $z\in\partial G$.

Now we are able to state our main result, showing that \textit{ISS of nonlinear parabolic systems w.r.t. a boundary input can be reduced to the problem of ISS of a parabolic system with a distributed and constant input, which is conceptually much simpler.}
\begin{theorem}
\label{thm:ISS_Parabolic_Systems_Characterization}
Suppose that assumptions (H1), (H2) and (H3) hold and let $p\in[1,+\infty]$. Let us endow the linear spaces $X,Y$ with the norm $\left\| \cdot \right\| _{p}$. The following statements are equivalent:

\begin{itemize}
 \item[(i)] The system \eqref{GrindEQ__22_}, \eqref{GrindEQ__23_}, \eqref{GrindEQ__24_}  with the state space $X$ is ISS w.r.t. inputs of class $\Uc$.
 \item[(ii)] The system \eqref{GrindEQ__22_}, \eqref{GrindEQ__23_}, \eqref{GrindEQ__24_} with the state space $X$ is ISS w.r.t. constant in time and space inputs of class $\Uc$.
 \item[(iii)] The system \eqref{eq:Nonlinear_Parabolic_Equation_transformed}, \eqref{eq:Dirichlet_transformed}
 with the state space $Y$ is ISS w.r.t. inputs in $\Vc$.
\end{itemize}
\end{theorem}

Similarly, the following result on exp-ISS property holds:
\begin{theorem}
\label{thm:exp-ISS_Parabolic_Systems_Characterization}
Suppose that assumptions (H1), (H2) and (H3) hold and let $p\in[1,+\infty]$. Let us endow the linear spaces $X,Y$ with the norm $\left\| \cdot \right\| _{p}$.
The following statements are equivalent:

\begin{itemize}
 \item[(i)] The system \eqref{GrindEQ__22_}, \eqref{GrindEQ__23_}, \eqref{GrindEQ__24_}  with the state space $X$ is exp-ISS w.r.t. inputs of class $\Uc$.
 \item[(ii)] The system \eqref{GrindEQ__22_}, \eqref{GrindEQ__23_}, \eqref{GrindEQ__24_} with the state space $X$ is exp-ISS w.r.t. constant in time inputs of class $\Uc$.
 \item[(iii)] The system \eqref{eq:Nonlinear_Parabolic_Equation_transformed}, \eqref{eq:Dirichlet_transformed}
 with the state space $Y$ is exp-ISS w.r.t. inputs in $\Vc$.
\end{itemize}
\end{theorem}

\section{ISS of linear parabolic systems with boundary inputs}
\label{sec:Lin_Parabolic_Systems}

In this section, we apply the above results to the linear problems with boundary disturbances over multidimensional spatial domains, which are of specific interest.
We continue to assume that $G\subset \R^n$ is an open connected and bounded set with the smooth boundary
and $D:=(0,+\infty)\times G$.

Consider the linear heat equation with a Dirichlet boundary input:
\begin{equation}
\label{eq:IBVP_Operator_L_Linear}
\begin{array}{l}
\frac{\partial x}{\partial t}(t,z) = \Delta x(t,z) + ax(t,z),\quad (t,z) \in D,  \\
x(t,z) = u(t,z), \quad (t,z)\in (0,+\infty) \times \partial G,\\
x(0,z) = x_0(z), \quad z \in \overline{G},
\end{array}
\end{equation}
where $\Delta$ is a Laplacian and $a \in \R$.

\textit{We are going to show that under reasonable assumptions ISS of the system \eqref{eq:IBVP_Operator_L_Linear}
is equivalent to its uniform global asymptotic stability in absence of inputs.}

The input $u:\R_+ \times\partial G\to \R$ is the trace of a function $\nu\in C^{0}(\R_+\times\overline{G}) \bigcap C^{1,2}((0,+\infty)\times G)$.
Defining the function
\vspace{-1.5mm}
\begin{eqnarray}
f(t,z):=\Delta \nu(t,z) - \frac{\partial \nu}{\partial t}(t,z),\quad t\geq0,\ z\in G
\label{eq:f_def}
\end{eqnarray}
\vspace{-0.5mm}
and using the transformation
\vspace{-1.5mm}
\begin{eqnarray}
x(t,z)=e^{at}\big(y(t,z)+\nu(t,z)\big),\quad t\geq 0 ,\ z\in \overline{G}
\label{eq:Aux_Transformation}
\end{eqnarray}
\vspace{-0.5mm}
we are in a position to study an equivalent to \eqref{eq:IBVP_Operator_L_Linear} initial boundary value problem
\vspace{-1.5mm}
\begin{equation}
\label{eq:IBVP_Operator_L_Linear_modified}
\begin{array}{l}
\frac{\partial y}{\partial t}(t,z) = \Delta y(t,z) + f(t,z),\quad (t,z) \in D, \\
y(t,z) = 0, \quad (t,z)\in (0,+\infty) \times \partial G,\\
y(0,z) = x(0,z)-\nu(0,z), \quad z \in \overline{G}.
\end{array}
\end{equation}
Let $m>0$ be the smallest integer for which $m\geq \frac{1+[n/2]}{2}$. \cite[Theorem 6, p. 365]{Eva98} in conjunction with \cite[Theorem 4, p. 288]{Eva98}
and \cite[Theorem 6, p. 270]{Eva98} guarantees that if
\begin{itemize}
 \item[(p-i)] $y[0]\in H^{2m+1}(G)$, $\frac{d^k}{dt^k}(f[t]) \in L^2(0,T;H^{2m-2k}(G))$ for every $k=0,\ldots,m$ and $T>0$,
 \item[(p-ii)] $g_i \in H^1_0(G)$ for $i=0,\ldots,m$, where $g_0:=y[0]$, $g_m:=\frac{d^{m-1}}{dt^{m-1}}f[0]+\Delta g_{m-1}$,
\end{itemize}
then the initial boundary value problem \eqref{eq:IBVP_Operator_L_Linear_modified} has a unique solution $y\in CL$.

Therefore, using the transformation \eqref{eq:Aux_Transformation}, we conclude that for every $x_0\in H^{2m+1}(G)$ and for every input $u\in \R_+\times \partial G \to \R$ being the trace of a function $\nu\in C^{0}(\R_+\times\overline{G}) \bigcap C^{1,2}((0,+\infty)\times G)$ and satisfying (p-i), (p-ii),
also the initial boundary value problem \eqref{eq:IBVP_Operator_L_Linear} has a unique solution $x\in CL$.
Hence, \eqref{eq:IBVP_Operator_L_Linear} defines a control system with
\begin{itemize}
 \item $X:=H^{2m+1}(G)$ with $\|\cdot\|_p$-norm, for any fixed $p\in[1,+\infty]$.

 \item  $\Uc(x)$ being the set of all inputs $u:\R_{+} \times \partial G\to \R$ which are traces of functions $\nu\in C^{0}(\R_+\times\overline{G}) \bigcap C^{1,2}((0,+\infty)\times G)$ satisfying $\frac{d^{k} }{d\, t^{k} } \left(f[t]\right)\in L^{2} \big(0,T;H^{2m-2k} (G)\big)$ for every $k=0,...,m$ and $T>0$ and $g_{i} \in H_{0}^{1} (G)$ for $i=0,...,m$, where $f$ is defined by
 \eqref{eq:f_def}, $g_{0} :=x-v[0]$, $g_{1} :=f[0]+\Delta g_{0} $,\dots , $g_{m} :=\frac{d^{m-1} }{d\, t^{m-1} } \left(f[0]\right)+\Delta g_{m-1} $,

 \item $\phi(t,x_0,u)$ being the unique solution $x[t]$ of the initial boundary value problem \eqref{eq:IBVP_Operator_L_Linear}.
\end{itemize}
Notice that $X$ contains the constant functions. We conclude that $\Sigma$ satisfies (H1). Clearly, (H2) holds for $\Sigma$ as well.

Hence we obtain from Theorem~\ref{thm:exp-ISS_Parabolic_Systems_Characterization}:
\begin{corollary}
\label{cor:ISS_Linear_Heat_Equation}
Assume that $G \subset \R^n$ is an open bounded set with a smooth boundary for which Assumption (H3) holds.
Then $\Sigma=(X,\Uc,\phi)$ is exp-ISS with a linear gain function iff $\Sigma$ is 0-UGAS.
\end{corollary}

\begin{proof}
Clearly, if $\Sigma=(X,\Uc,\phi)$ is exp-ISS, then $\Sigma=(X,\Uc,\phi)$ is 0-UGAS.

Now assume that $\Sigma=(X,\Uc,\phi)$ is 0-UGAS and let us prove the converse implication.
\eqref{eq:IBVP_Operator_L_Linear} is a problem \eqref{GrindEQ__1_}, corresponding to the operator
$Lx:= \frac{\partial x}{\partial t} - \Delta x - ax$.

According to Theorem~\ref{thm:exp-ISS_Parabolic_Systems_Characterization}, exp-ISS of \eqref{eq:IBVP_Operator_L_Linear} is equivalent to exp-ISS of the system
\begin{eqnarray}
\frac{\partial y}{\partial t} = \Delta y +ay+av,\quad t >0, \ z \in G,
\label{eq:Nonlinear_Parabolic_Equation_transformed-2}
\end{eqnarray}
with homogeneous Dirichlet boundary condition \eqref{eq:Dirichlet_transformed} and constant inputs $v \in \R$.

In order to prove the claim, we are going to use the semigroup approach. Consider three cases: $p\in (1,+\infty)$, $p=1$ and $p=+\infty$.
The operator $A_p$, $p\in (1,+\infty)$ with a domain of definition $D(A_p):=W^{2,p}(G) \bigcap W^{1,p}_0(G)$ and $A_p x:=\Delta x + ax$, for $x\in D(A_p)$ generates an analytic semigroup over $L^p(G)$, see \cite[Theorem 3.6, p. 215]{Paz83}.

For $p=1$ the operator $A_1$ with a domain of definition $D(A_1):= \{x\in W^{1,1}_0(G):\Delta x \in L^1(G)\}$ and $A_1 x:=\Delta x + ax$, for $x\in D(A_1)$ generates an analytic semigroup over $L^1(G)$, see \cite[Theorem 3.10, p. 218]{Paz83}.
Analogously, one can define an operator $A_\infty:=\Delta + aI$ with a certain domain of definition $D(A_\infty)$, which generates an analytic semigroup over $C_0(\overline{G}):=\{x\in C(\overline{G}): x(z)=0 \ \forall z\in\partial G\}$, see \cite[Theorem 3.7, p. 217]{Paz83}.

Now, since we assume that $\Sigma$  is 0-UGAS in the norm $\|\cdot \|_p$, then also the operator $A_p$ generates 0-UGAS (exponentially stable) $C_0$-semigroup $T_p$, which follows since $X$ is dense in the $L^p(G)$ for $p\in[1,\infty)$ and in $C_0(\overline{G})$ for $p=\infty$ and since $T_p$ is a semigroup of bounded operators (hence the norm of the operator $T_p(t)$ with a domain restricted to $Y$ is equal to the norm of $T_p(t)$ as an operator on $X$).

Since $T_p$ is an exponentially stable semigroup, \cite[Proposition 3]{DaM13} ensures that \eqref{eq:Nonlinear_Parabolic_Equation_transformed-2} is exp-ISS with a linear gain function.
An inspection of Theorem~\ref{thm:exp-ISS_Parabolic_Systems_Characterization} shows that \eqref{eq:IBVP_Operator_L_Linear} is exp-ISS with a linear gain function as well.
\end{proof}

\begin{remark}\em
\label{rem:Linfinity}
Note that the operator $\Delta + aI$ in the above proof does not generate a strongly continuous semigroup over the space $L^\infty(G)$, see \cite[p.217]{Paz83} or \cite[Lemma 2.6.5, Remark 2.6.6]{CaH98}, therefore it is of importance to define $A_\infty$ as a generator of strongly continuous semigroup over $C_0(\overline{G})$.
\end{remark}

\begin{remark}\em
\label{rem:Practical_applications}
There are different ways to ensure that $\Sigma$ is 0-UGAS. For $p\in(1,+\infty)$ a usual way would be to construct a Lyapunov functional.
On the other hand, since $T_p$ is an analytic semigroup for any $p\in[1,+\infty]$,
$T_p$ is exponentially stable (i.e. 0-UGAS) if and only if the spectrum of $A_p$ lies in $\{z\in \C:\text{Re}z<0\}$, see e.g. \cite[p.387]{Tri75}.
\end{remark}

\section{Conclusions}

We presented a new technique for analyzing ISS of linear and nonlinear parabolic equations with boundary inputs.
We prove that parabolic equations with Dirichlet boundary inputs are monotone control systems and use this fact to transform the parabolic system with boundary disturbances into a related system with distributed constant disturbances. We show that ISS of the original equation is equivalent to ISS of the transformed system.
Analysis of ISS of the transformed system is much easier to perform, for example, by means of Lyapunov methods.

Although in this paper we concentrate on parabolic scalar equations and study properties of classical solutions of such equations, the scheme which we have developed here can be useful for other classes of monotone control systems: monotone parabolic systems, ordinary differential equations, ODE-heat cascades, some classes of time-delay systems. We expect that a big part of our analysis can be transferred to study properties of mild solutions of parabolic systems. Finally, some of our results can be used to study ISS of monotone systems of a general nature.

\vspace{-2mm}


text:
\section{Introduction}
Since the beginning of last century, radar technology has been continuously developed and benefiting us in many different ways: From large scale surveillance radar systems that are crucial for aircraft safety and space exploration \cite{Skolnik2001}, to small scale systems for monitoring insects \cite{Oneal2004}. Considering the limit of radar tracking technique, it is intuitive to ask: How small can an object be accurately tracked by a radar system? Can it be as small as a tracer particle with a size comparable to a grain of sand? If this challenge can be solved to a satisfactory level, radar technology can help us greatly in a better monitoring, understanding, and predicting granular flow problems that widely exist in nature, industry and our daily lives, ranging from geo-science (e.g. landslide, debris flow), through chemical and civil engineering (e.g. pile drilling), to space exploration (e.g. landing on an asteroid) \cite{Duran2000,pg13}. The reason behind is that the mobility of single particles in a granular material can influence dramatically its collective behavior, owing to its discrete nature as well as heterogeneous distributions of force chains inside \cite{Bi2011}.

In the past decades, there have been substantial progresses in imaging granular particles \cite{Amon2017}. Due to the fact that most of the granular materials are opaque, optical means for imaging particles in three dimensions (3D), such as refractive index matching \cite{Dijksman2017}, are very limited. Instead, X-ray tomography \cite{Athanassiadis2014} and Magnetic Resonance Imaging \cite{Stannarius2017} have both been used to identify the internal structures of granular materials. However, the limited time resolution of scanning technique as well as the huge amount of data to be processed hinder the investigation of granular dynamics. For the investigation of granular dynamics, alternative approaches including radar systems for tracer particles have also been discussed \cite{Hill2010}.

Recently, we introduced a small scale continuous-wave (CW) radar system working at X-band to track a spherical object with a size down to $5$\,mm \cite{Ott2017}. In comparison to other techniques, the continuous trajectory of a tracer particle obtained by the radar system helps in deciphering granular dynamics greatly. Here, we characterize the uncertainty of this system, discuss possible sources of error and ways to improve the accuracy.

\section{Experimental set-up and tracking algorithm}
\begin{figure}[!ht]

\centering
\includegraphics[width=0.85\columnwidth]{fig1.eps}
\caption{Block diagram of the radar system (a) and a sketch showing the configuration of the transmission (Tx.) and receiving (Rx.) antennae as well as the tracer holding device.}
\label{fig:setup}
\end{figure}

Figure~\ref{fig:setup} shows the block diagram of the radar system (a) together with a sketch of the set-up used to test the accuracy of reconstructed trajectories (b). The bi-static radar system operates at 10\,GHz (X-band) with one transmission (Tx.) antenna pointing in the direction of gravity (defined as $-z$ direction) and three receiving (Rx.) antennae mounted symmetrically around the $z$ axis. Polarized electromagnetic (EM) waves, after being scattered by the tracer particle, are captured by the Rx. antennae. With the help of an IQ mixer, which compares the phase shift between the Tx. and Rx. antennae, the change of the absolute traveling distance for the $i$th antenna $L_i = l_0 + l_i$ can be obtained, where $l_0$ and $l_i$ are the distance between Tx. Ant. and the target and the distance from the target to the Rx. Ant, respectively. Subsequently, a transformation matrix $\vec T$ is applied to the distance vector $\vec L = (L_1,L_2,L_3)$ to reconstruct the tracer trajectory in a 3D Cartesian system. The smallest spherical object being traceable by the system was found to be $\sim 5$\,mm, in agreement with the prediction of radar equation. More detailed descriptions of the radar system can be found in \cite{Ott2017}.

IQ-Mixers play an essential role in the accurate ranging of a target. The LO and RF inputs correspond to the signal sent to the Tx. antenna [$a\cos(2\pi f_{\rm 0}t)$] and that received by a Rx. antenna [$b\cos(2\pi ft + \theta)$], where $a$ and $b$ are the magnitudes of the corresponding signals, $f_{\rm 0}$ and $f$ are the transmitted and received signal frequencies. The output signals of the IQ-Mixer are

\begin{eqnarray}
I=\frac{ab}{2}\cos[2\pi(f_{\rm 0}-f)t - \theta], \nonumber \\
Q=\frac{ab}{2}\sin[2\pi(f_{\rm 0}-f)t - \theta]. \label{eq:iq2}
\end{eqnarray}

\noindent Subsequently, the relative movement of the tracer is obtained from the phase shift of $I+Qi$ in a complex plane. If $L_i$ varies with a distance of one wavelength, the vector $I+Qi$ rotates $2\pi$. As IQ mixers provide analogue signals that representing the mobility of the tracer, the time resolution of the radar system is only limited by the analogue-digital (AD) converter.

Although distance measures rely only on the phase information, the sensitivity and accuracy of the system depend on IQ signal strength. In order to have a sufficient signal to noise ratio, the directions of the horn antennae (Dorado GH-90-20) are adjusted with the help of a laser alignment and range meter (Umarex GmbH, Laserliner) to face the target area. According to the specification of the antenna, the main lobe of its radiation pattern has an opening angle of $\sim15$ degrees. Thus, we estimate the field of `view' of the radar system has a volume of about $30 {\rm cm} \times 30 {\rm cm} \times 30 {\rm cm}$, taking into account the average working distance of the antennae. The distance between each antenna and the center of the coordinate system is also measured by the laser meter during the adjustment process. The polarization of the antennae are adjusted to maximize the raw I and Q signals. The whole system is covered with microwave absorbers (Eccosorb AN-73) to reduce clutter and unwanted noises from the surrounding. In addition, the container for granular materials and the holder of the tracer are made of expanded polystyrene (EPS), which is transparent to EM waves.

A metallic sphere with a diameter $d=$10\,mm is used as the tracer. It is initially held by a thin thread wrapped around and released by gently pulling the thread such that the initial velocity of the falling sphere is close to $0$. This design enables a defined and reproducible initial condition for a comparison among various experimental runs. The raw IQ signals from the AD converter (NI DAQPAD-6015) are recorded and further processed with a Matlab program to obtain the reconstructed trajectories.

\section{Data analysis and error correction}

\begin{figure}
\centering
\vskip -2em
\includegraphics[width=0.85\columnwidth]{fig2.eps}
\caption{Raw (continuous lines) and corrected (open symbols) signals representing a free-falling sphere from a height of $27$\,cm. Red (dark red), green (dark green) and blue (dark blue) curves (points) correspond to the results from channel 1, 2, and 3, respectively. For a better visibility, the offsets of the raw signals are removed. }
\label{fig:raw}
\end{figure}

As Fig.\,\ref{fig:raw} shows, the raw IQ signals are typically not ideal in the sense that the in-phase (I) and quadrature (Q) signals are not always orthogonal with each other. This mismatch may arise from the DC offsets of either I or Q signal, gain and phase imbalance. How to correct such kind of errors has been extensively discussed in, for instance \cite{Churchill1981} or \cite{Huang2002}, particularly along with the development of telecommunication and non-invasive motion detecting techniques \cite{Park2007}. The distortions are typically attributed to device imperfections as well as clutter. However, for the system being used here, there are additional errors arising from the mobility of the tracer itself, which can not be readily corrected with a calibration of the hardware. Moreover, distortion may also arise from the interaction of the scattered signal from the tracer with that from static objects that are not completely transparent to EM waves. In that case, the existence of `mirrored' particles may lead to additional uncertainty.

Here, we use the following approach to correct IQ mismatch arising from multiple sources of errors. It works best when the object moves in a distance covering multiple wavelengths. As illustrated in Fig.\,\ref{fig:corr}, the correction process is composed of the following steps: First, we identify the time segment of the raw data $V_{\rm raw}$ that contains the movement of the tracer particle via finding the start and end of the fluctuations. Second, the peaks (red circles) and valleys (blue circles) of individual fluctuations are determined by finding the local extreme values in the selected data. Third, the bias error $V_{\rm bias}$ [green line in (a)] is estimated as the mean value of the spline fits of peaks and valleys (dashed lines). In order to avoid unrealistic extrapolations, the bias error starts to vary only from the first peak. Fourth, the bias error is removed and the corrected signal $V_{\rm raw - bias}$ is segmented by zero crossings. Finally, the data in individual segments are rescaled by local maxima and minima to correct gain mismatch.

As shown in Fig.\,\ref{fig:raw}, this approach can effectively find time dependent correction factors due to tracer movement. For the corrected data of channel 3 (dark blue circles), there exists a slight deviation from a perfect circle, indicating the existence of a small phase error. This arises presumably from the fact that perfect polarization cannot be achieved for all three Rx. antennae. Further investigations are needed to check whether this error can be avoided by using circular polarized EM waves or by correcting the phase error between I and Q signals in the Matlab program.

\begin{figure
\centering
\includegraphics[width=0.85\columnwidth]{fig3.eps}
\caption{Process for IQ mismatch correction. (a) A representative raw signal with peaks and valleys marked with red and blue open circles. From an average of both spline fits for the peaks and valleys, the bias error (green line) for the raw signal as a function of time is estimated. (b) Bias corrected signal time dependent rescaling factors (green line) for the correction of gain error. (c) Corrected signal for further analysis.}
\label{fig:corr}
\end{figure}

After the correction of IQ mismatch, the corresponding phase angles are obtained by $\phi = \arctan(Q/I)$. Because $\phi$ is a modulo of $2\pi$, a further correction on the phase jump is needed to obtain the continuous phase shift $\Phi$. In this step, a threshold is introduced to determine whether a phase jump occurs or not and in which direction the jump takes place. As the phase shift of the $i$th channel $\Phi_i \propto L_i$, the variation of $\Phi$ with time (see the blue curve in Fig.\,\ref{fig:dist}) indicates that the target object moves initially slow and accelerates while moving away from the antennae. As demonstrated by a comparison between corrected $\Phi$ and uncorrected $\Phi_{\rm uncorr}$ phase shift, the aforementioned correction method can effectively reduce unrealistic fluctuations in the reconstructed curves. As shown in the inset of Fig.\,\ref{fig:dist}, the distance $L_1$ obtained from Rx. antenna represents exactly what is expected: The object falls freely with a growing velocity and bounces back when reaching the container bottom, suggesting that the coefficient of restitution, which measures the relative rebound over impact velocities, can be determined with the radar system. In comparison to the standard high speed imaging technique \cite{Mueller2016}, the radar tracking technique requires less data collection and processing efforts.

\section{Validation of reconstructed trajectory}

\begin{figure
\vskip -3em
\centering
\includegraphics[width=0.85\columnwidth]{fig4.eps}
\caption{Arc-tangential demodulation process to obtain the traveling distance from the Tx. to a Rx. antenna. The red open symbols correspond to the outcome from the demodulation and the blue curve represents the continuous phase shift that scales with the traveling distance of an EM wave. The gray curve corresponds to the $\Phi$ without correcting IQ mismatch. Inset shows an example of the variation of $L_{\rm 1}$ over a longer time.}
\label{fig:dist}
\end{figure}

From the measured distances $\vec{L} \equiv (L_1, L_2, L_3)$, the reconstructed trajectory can be obtained with a coordinate transformation

\begin{equation}
\label{eq:reconst}
\left(
\begin{array}{c}
x \\
y \\
z
\end{array}
\right)=\frac{\vec{r}-\vec{L}}{\vec{T}},
\end{equation}

\noindent where the vector $\vec{r}$ is chosen to be $0$ as it contributes only to a constant offset to the reconstructed trajectory, the transformation matrix reads

\begin{equation}
\label{eq:trmat}
\vec{T} \equiv \left(
\begin{array}{ccc}
\sin{\theta_{\rm 1}}\cos{\phi_{\rm 1}} & \sin{\theta_{\rm 1}}\sin{\phi_{\rm 1}} & 1+\cos{\theta_{\rm 1}} \\
\sin{\theta_{\rm 2}}\cos{\phi_{\rm 2}} & \sin{\theta_{\rm 2}}\sin{\phi_{\rm 2}} & 1+\cos{\theta_{\rm 2}} \\
\sin{\theta_{\rm 3}}\cos{\phi_{\rm 3}} & \sin{\theta_{\rm 3}}\sin{\phi_{\rm 3}} & 1+\cos{\theta_{\rm 3}}
\end{array}\right)
\end{equation}

\noindent with $\theta_i$ and $\phi_i$ the tilt and azimuth angles of the $i$th antenna, respectively. The transformation matrix is determined from a calibration process using the same tracer particle moving in a known trajectory. A detailed description of the calibration and coordinate transformation processes can be found in \cite{Ott2017}.

\begin{figure
\vskip -4em
\centering
\includegraphics[width=0.85\columnwidth]{fig5.eps}
\caption{A comparison of reconstructed free-falling curves at various initial falling heights. The solid line corresponds to the expected free-falling curve for the largest falling height. Note that the curves for various $H$ are shifted to have the initial falling position $z=0$\,cm, and only the trajectories before the first bouncing with the container bottom are shown except for $H=22.37$\,m. For each curve, one over 15 data points are shown here for a better visibility.}
\label{fig:varHeight}
\end{figure}

Finally, we compare the reconstructed trajectories of a free-falling object from different initial falling heights with the expected parabolic curve. As shown in Fig.\,\ref{fig:varHeight}, the falling curves agree with the expected curve well, demonstrating that, after a proper correction of IQ mismatch, the radar system can be used for particle tracking. Further tests with the Styrofoam container filled with expanded polypropylene (EPP) particles also show that this system can be readily used for the investigation of granular dynamics. Further investigations will focus on particle tracking with various types of granular materials, particularly how to deal with distorted signals arising from the multiple scattering of the surrounding particles.

\section{CONCLUSION}

To summarize, this investigation suggests that advances in radar tracking technology can be helpful in the investigation of granular dynamics. Using an X-band continuous wave radar system, we are able to track a centimeter sized metallic object in 3D, which enables, for instance, a measurement of the coefficient of restitution of the particle. In comparison to other particle imaging techniques already being used for granular particles \cite{Amon2017}, continuous-wave radar tracking has the advantage of high time resolution and low data collection and processing requirements. With the rapid development of radar technology, this approach is also expected to be more cost effective and accurate.

Moreover, we show that the accuracy of the radar tracking technique depends strongly on a proper correction of  IQ mismatch, which arises predominately from the mobility of the tracer itself. A practical approach has been proposed to correct the instantaneously changing bias as well as gain errors in the raw IQ signals. Finally, we validate this approach through an analysis on the reconstructed trajectories of a free-falling sphere.

\section*{Acknowledgment}

We acknowledge Felix Ott for his preliminary work on the experimental set-up and Klaus Oetter for technical support. Helpful discussions with Valentin Dichtl, Simeon V\"olkel and Ingo Rehberg are gratefully acknowledged. This work is partly supported by German Research Foundation through Grant No.~HU1939/4-1.

redpajama_set_name: RedPajamaArXiv


text:
\section{Evidential Value in Forensic Statistics}\label{EVFS}

At some crime scene a trace has been found that links a suspect to
the crime. In the court case the prosecutor puts forward the
hypothesis $H_p$ that the suspect is the donor of the trace. The
defendant claims the hypothesis $H_d$ holds, which states that an
unknown person, not the suspect, is the donor of the trace. The
juror (judge, jury) has to decide in favor of $H_p$ or $H_d.$ An
important current scientific approach to such criminal court cases
is via the so-called {\em Bayesian Paradigm of Forensic
Statistics.}

Within this paradigm the juror has to construct a prior opinion
about $H_p$ and $H_d.$ This means that the juror has to decide
beforehand, before seeing the evidence, how likely the hypothesis
of the prosecutor is in comparison to the hypothesis of the
defendant. This prior opinion might be based on e.g. the number of
possible offenders, and it may be formulated in terms of the prior
odds in favor of the hypothesis of the prosecutor, namely
$$ P(H_p)\,/\,P(H_d).$$
The evidence in such a court case consists of the trace found at
the crime scene and characteristics of the suspect. Let us denote
it by $E.$ The forensic expert has to determine now the
probability that a randomly chosen person would leave a trace like
the one found, at the crime scene. This probability is denoted by
$P(E\,|\,H_d).$ Likewise he has to determine $P(E\,|\,H_p),$ the
probability that the suspect would leave a trace like the one
found, at the crime scene. The ratio
$$P(E\,|\,H_p)\,/\,P(E\,|\,H_d)$$ is called the likelihood ratio.
Multiplying the prior odds and the likelihood ratio the juror
obtains the so-called posterior odds in favor of the hypothesis of
the prosecutor $$P(H_p\,|\,E)\,/\,P(H_d\,|\,E),$$ i.e., the odds
in favor of $H_p$ after having seen the evidence. The juror has to
base his decision on these posterior odds. In summary, the
Bayesian Paradigm of Forensic Statistics reads as follows
\begin{equation}\label{BayesParadigm}
 \underbrace{\frac{P(H_p)}{P(H_d)}}_{prior\, odds}\ \
\overbrace{\frac{P(E\,|\,H_p)}{P(E\,|\,H_d)}}^{likelihood\,
ratio}\, =\
\underbrace{\frac{P(H_p\,|\,E)}{P(H_d\,|\,E)}}_{posterior\,
odds}\,.
\end{equation}\\
The validity of equation (\ref{BayesParadigm}) may be checked
straightforwardly by applying the definition of conditional
probability, which is $$P(A\,|\,B) = P(A \cap B)\,/\,P(B),$$ where
$A \cap B$ is the intersection of $A$ and $B.$ Since the
likelihood ratio in (\ref{BayesParadigm}) may be interpreted as
the weight that the evidence should have in the decision of the
juror, it is often called the evidential value in forensic
science.

The evidence $E$ is viewed here as a realization of a random
mechanism, both under $H_d$ and $H_p.$ In case this random
mechanism produces outcomes via probability density functions
$f(E\,|\,H_p)$ and $f(E\,|\,H_d),$ the probabilities in the
likelihood ratio or evidential value are replaced by the
corresponding probability density functions, resulting in
\begin{equation}\label{BayesParadigmdensity}
 \underbrace{\frac{P(H_p)}{P(H_d)}}_{prior\, odds}\ \
\overbrace{\frac{f(E\,|\,H_p)}{f(E\,|\,H_d)}}^{likelihood\,
ratio}\, =\
\underbrace{\frac{P(H_p\,|\,E)}{P(H_d\,|\,E)}}_{posterior\,
odds}\,.
\end{equation}\\

\section{Modelling Fabrication of Data Underlying an ANOVA Study}\label{MFD}

In Analysis of Variance the basic assumption is that all
observations may be viewed as realizations of independent normally
distributed random variables with the same variance $\sigma^2$ and
with means that depend on the values of some categorical
covariates. Let $I$ be the total number of cells that are defined
via these categorical covariates, and let the number of
observations per cell be the same, namely $n.$ The random
variables denoting the observations are then
\begin{equation}\label{model}
X_{ij}= \mu_i +\varepsilon_{ij}\,,\quad i=1,\dots,I,\, j=1,\dots,
n.
\end{equation}
The cell means $\mu_i$ are unknown real numbers, and the
measurement errors $\varepsilon_{ij}$ are independent, normally
distributed random variables with mean 0 and variance $\sigma^2.$

If authors are fiddling around with data and are fabricating and
falsifying data, they tend to underestimate the variation that the
data should show due to the randomness within the model. Within
the framework of the above ANOVA case, we model this by
introducing dependence between the normal random variables
$\varepsilon_{ij},$ which represent the measurement errors.
Actually, we assume that the measurement errors in any cell have
correlation coefficient $\rho$ with respect to the corresponding
measurement errors in the other cells. More precisely formulated,
we assume that the correlations between the random variables
$\varepsilon_{ij}$ no longer all vanish, but satisfy
\begin{equation}\label{dependence}
\rho(\varepsilon_{ij}, \varepsilon_{hj})= \rho,\quad
j=1,\dots,n,\, 0\leq i \neq h \leq I,
\end{equation}
with all other correlations still being equal to 0. In the sequel
we restrict attention to nonnegative values of $\rho$ and we
exclude $\rho=1$ for technical reasons, so $0 \leq \rho < 1.$ We
note that under the standard assumptions of ANOVA $\rho=0$ holds.
Furthermore, we note that within cells observations may be
renumbered in order to get the structure (\ref{dependence}).
Nevertheless, we still assume (\ref{model}) to hold and the
measurement errors to be normally distributed with mean 0 and
variance $\sigma^2.$

A way in which fabrication of measurement errors may take place is
by copying some of them. This might be modelled as follows. Let
$U_j\,,\, j=1, \dots, n,$ and $V_{ij}\,,\, i=1, \dots, I,\ j=1,
\dots, n,$ be independent and identically distributed normal
random variables with mean 0 and variance $\sigma^2.$ Independent
of these, let the random indicators $\Delta_{ij}\,,\, i=1, \dots,
I,\ j=1, \dots, n,$ be independent and identically distributed
Bernoulli random variables with $P(\Delta_{ij}=1)=\sqrt \rho$ and
$P(\Delta_{ij}=0)=1- \sqrt \rho.$ Then
\begin{equation}\label{structure}
\varepsilon_{ij}=\Delta_{ij} U_j + (1-\Delta_{ij}) V_{ij}\,,\quad
i=1, \dots, I,\ j=1, \dots, n,
\end{equation}
satisfy (\ref{dependence}) and (\ref{model}). Note that for $0\leq
i \neq h \leq I$ we have $\varepsilon_{ij}= \varepsilon_{hj}= U_j$
with probability ${\sqrt \rho}^2 = \rho$ then, and the measurement
errors satisfy (\ref{dependence}).

Finally, we note that (\ref{dependence}) is just one possible way
to model dependence, and that the actual way in which fabrication
has been implemented, might lead to quite different dependence
structures. However, this model will come close to some types of
fabrication and falsification.

\section{Evidential Value for Fabrication of Data Underlying an ANOVA Study}\label{EVFD}

Consider a study in a scientific research paper. The data in this
study are analyzed by ANOVA and presented via the sample averages
of the cells and the values of some F-statistics. The underlying
data themselves are not published and are not available. The
conclusion of this study is that the $I$ cells can be grouped into
$K$ groups of $I_k$ cells ($\sum_{k=1}^K I_k=I$), such that
(possibly after renumbering of the cells) group $k$ consists of
cells $i=L_{k-1}+1, \dots, L_k,\, k=1,\dots,K,$ with $0=L_0 < L_1<
\dots < L_K=I,\ L_k -L_{k-1}=I_k,$ and such that for each group
the population cell means are the same, i.e.,
\begin{equation}\label{conclusion}
\mu_i=\nu_k,\quad i=L_{k-1}+1,\dots,L_k,\quad k=1,\dots,K,
\end{equation}
for some values $\nu_k,\, k=1,\dots,K.$

There are two hypotheses to be formulated about the data
underlying the ANOVA study. The hypothesis $H_p$ of fabrication of
the data underlying the results presented in the paper, is
$0<\rho<1.$ The other hypothesis $H_d$ represents the situation
that data have been collected according to (\ref{model}) with
independent $X_{ij},$ i.e., $\rho=0.$ We want to determine the
evidential value of the ANOVA study, i.e., of the sample means of
the cells and the published F-statistics, with respect to these
hypotheses $H_p$ and $H_d.$

To this end we first note that the sample averages in the cells,
\begin{equation}\label{sampleaverage}
X_{i\cdot}=\frac 1n \sum_{j=1}^n X_{ij},\quad i=1,\dots,I,
\end{equation}
have a joint $I$-dimensional multivariate normal distribution.
Actually, the dependence structure (\ref{dependence}) implies
\begin{equation}\label{multinormal}
\begin{pmatrix}
X_{1\,\cdot} \cr \cdot \cr \cdot \cr X_{I\,\cdot}\cr
\end{pmatrix}
\sim \cal N \left(
\begin{pmatrix}
\mu_1 \cr \cdot \cr \cdot \cr \mu_I\cr
\end{pmatrix}
,\, \mbox{$\sigma^2 n^{-1}$}
\begin{pmatrix} 1& \rho & \cdot & \rho \cr \rho & 1 & \cdot & \rho \cr
\cdot & \cdot & \cdot & \cdot \cr \rho & \rho & \cdot & 1 \cr
\end{pmatrix}
\right).
\end{equation}
In stead of assuming (\ref{dependence}), we could have started
right away from (\ref{multinormal}).

Since the inverse of the covariance matrix in (\ref{multinormal})
equals
\begin{equation}\label{inverse}
\mbox{\Large{$\frac n{\sigma^2 (1-\rho)(1+(I-1)\rho)}$}}
\begin{pmatrix} 1+(I-2)\rho& -\rho & \cdot & -\rho \cr -\rho & 1+(I-2)\rho & \cdot & -\rho \cr
\cdot & \cdot & \cdot & \cdot \cr -\rho & -\rho & \cdot &
1+(I-2)\rho \cr
\end{pmatrix}
\end{equation}
and the determinant of $n\sigma^{-2}$ times this covariance matrix
equals
\begin{eqnarray}\label{determinant}
\lefteqn{
 \left|
\begin{array}{cccc}
    1 & \rho & \cdot & \rho \\
    \rho    & 1 & \cdot & \rho \\
    \cdot & \cdot & \cdot & \cdot \\
    \rho & \rho & \cdot & 1
\end{array} \right|
= \left|
\begin{array}{cccc}
    1 & \rho & \cdot & \rho \\
    \rho -1    & 1 -\rho & \cdot & 0 \\
    \cdot & \cdot & \cdot & \cdot \\
    \rho-1 & 0 & \cdot & 1-\rho
\end{array} \right| }\\
&&= \left|
\begin{array}{cccc}
    1+(I-1)\rho & \rho & \cdot & \rho \\
    0    & 1 -\rho & \cdot & 0 \\
    \cdot & \cdot & \cdot & \cdot \\
    0 & 0 & \cdot & 1-\rho
\end{array} \right| =  \left(1+(I-1)\rho \right)
(1-\rho)^{I-1},\nonumber
\end{eqnarray}
(\ref{multinormal}) and (\ref{conclusion}) entail that the joint
density of $X_{1\cdot},\dots,X_{I\cdot}$ equals
\begin{eqnarray}\label{jointdensity1}
\lefteqn{\left(\frac n{2\pi \sigma^2}\right)^{I/2}
\left[(1+(I-1)\rho)(1-\rho)^{I-1}\right]^{-1/2}} \nonumber
\\
 && \exp\left(-\frac n{2\sigma^2} \left[
 \frac 1{1-\rho}\sum_{k=1}^K\sum_{i=L_{k-1}+1}^{L_k}
 \left(X_{i\cdot}-\nu_k\right)^2
 \right. \right.  \\
&&  \qquad \qquad \qquad \left. \left.
  -\, \frac \rho {(1+(I-1)\rho)(1-\rho)} \left(\sum_{k=1}^K
 \sum_{i=L_{k-1}+1}^{L_k}\left(X_{i\cdot}-\nu_k\right)\right)^2
\right]\right). \nonumber
\end{eqnarray}
This density depends on the parameters $\rho, \sigma^2,
\nu_1,\dots, \nu_K.$ If the underlying data would be available
their mean square error
\begin{equation}\label{MSE}
\frac 1{I(n-1)} \sum_{i=1}^I
\sum_{j=1}^n\left(X_{ij}-X_{i\cdot}\right)^2
\end{equation}
would be the proper unbiased estimator of $\sigma^2.$ The
distribution of this estimator depends on $\rho,$ but its mean
does not. Furthermore, standard ANOVA theory shows that this
estimator is independent of the exponent in (\ref{jointdensity1}).
Since the underlying data are not available, the value of the
parameter $\sigma^2$ should be retrieved from the values of the
F-statistics given. For a method to do this that does not depend
on $\rho,$ see the next section. Let us call the resulting
estimate ${\hat \sigma}_n^2,$ and let us denote the density from
(\ref{jointdensity1}) with $\sigma$ replaced by ${\hat \sigma}_n$
by $f_n(X_{1\cdot},\dots, X_{I\cdot}; \nu_1,\dots,\nu_K,\rho).$

The hypothesis $H_p$ of fabrication of the data corresponds to the
parameter values $0<\rho<1$ and $\nu_1,\dots,\nu_K$ arbitrary, and
the hypothesis $H_d$ of proper data corresponds to the parameter
values $\rho = 0$ and $\nu_1,\dots,\nu_K$ arbitrary. The
evidential value
$$\frac{f(E\,|\,H_p)}{f(E\,|\,H_d)}$$ from
(\ref{BayesParadigmdensity}) in favor of $H_p$ versus $H_d$
becomes in this case (cf. Zhang (2009), Bickel (2012))
\begin{equation}\label{evidentialvalue1}
\mathbb V = \frac{\sup_{0<\rho<1,\, \nu_1,\dots,\nu_K \in \mathbb
R}f_n(X_{1\cdot},\dots, X_{I\cdot}; \nu_1,\dots,\nu_K,\rho)}
{\sup_{\nu_1,\dots,\nu_K \in \mathbb R}f_n(X_{1\cdot},\dots,
X_{I\cdot}; \nu_1,\dots,\nu_K, 0)}.
\end{equation}
Straightforward computation shows that for any $\rho$
\begin{equation}\label{nus}
\sup_{\nu_1,\dots,\nu_K \in \mathbb R}f_n(X_{1\cdot},\dots,
X_{I\cdot}; \nu_1,\dots,\nu_K,\rho)
\end{equation}
is attained at
$$\nu_k= {\bar X}_k =\frac 1{L_k -L_{k-1}}
\sum_{i=L_{k-1}+1}^{L_k} X_{i\cdot}\,,\quad k=1,\dots,K.$$ This
implies that the evidential value from (\ref{evidentialvalue1})
reduces to
\begin{equation}\label{evidentialvalue2}
\mathbb V = \sup_{0<\rho<1} \chi_n(\rho)
\end{equation}
with
\begin{eqnarray}\label{jointdensity2}
\lefteqn{\chi_n(\rho) =
\left[(1+(I-1)\rho)(1-\rho)^{I-1}\right]^{-1/2}}
\\
 &&  \qquad \quad \exp\left(-\frac {n \rho}{2{\hat \sigma}_n^2(1-\rho)}\
 \sum_{k=1}^K \sum_{i=L_{k-1}+1}^{L_k}\left(X_{i\cdot}-{\bar X}_k\right)^2
 \right). \nonumber
\end{eqnarray}
We need the additional notation
\begin{eqnarray}\label{notation}
S_n & = & \frac n{I{\hat \sigma}_n^2}\ \sum_{k=1}^K
\sum_{i=L_{k-1}+1}^{L_k
}\left(X_{i\cdot}-{\bar X}_k\right)^2,
\\
{\hat \rho}_n & = & \frac 12 (1-S_n)\left[1+\sqrt{1-\frac
{4S_n}{(I-1)(1-S_n)^2}}\right].
\end{eqnarray}
In Proposition \ref{computationV} of the Appendix the following is
shown.
\begin{itemize}
\item
If $$S_n \geq \frac {\sqrt I -1}{\sqrt I +1}$$ holds, then the
evidential value from (\ref{evidentialvalue1}) and
(\ref{evidentialvalue2}) reduces to $\mathbb V = 1.$
\item
If $$S_n < \frac {\sqrt I -1}{\sqrt I +1}$$ holds, then ${\hat
\rho}_n$ is well-defined and the evidential value from
(\ref{evidentialvalue1}) and (\ref{evidentialvalue2}) reduces to
\begin{equation}\label{evidentialvalue5}
\mathbb V = \max\left\{\chi_n({\hat \rho}_n), 1\right\}.
\end{equation}
\end{itemize}

\section{Estimating $\sigma^2$ from F-Statistics}\label{EFS}

Table 1 in Stapel, Koomen and Van der Pligt (1996) presents the
sample means in a three-way layout ANOVA study with a $3 \times 2
\times 2$ design.

\medskip
\begin{tabular}{|l|l|l|l|}\hline
 Prime type              & Positive & Negative & Irrelevant \\\hline
 Impersonal / Memory     & 2.3      & 3.5      & 2.9        \\
 Impersonal / Impression & 3.4      & 2.5      & 2.9        \\
 Personal / Memory       & 3.3      & 2.3      & 2.9        \\
 Personal / Impression   & 3.5      & 2.5      & 3.0        \\\hline
\end{tabular}
\medskip

The estimate of the error variance $\sigma^2$ is not given.
It should be possible to retrieve this estimate from the value of
any F-statistic. On page 441 of ibid. the value of the F-statistic
for testing three-way interactions is given, namely,
$F(2,326)=3.21.$ We assume that the 338 observations are
approximately uniformly distributed over the 12 cells. This yields
an average of 28.17 observations per cell. Applying e.g. Table
4.5.2 (Analysis of Variance of the Three-Way Layout with $M$
Observations per Cell) of Scheff\'e (1959) we obtain by some
computation that the mean square error for interaction equals
7.769. Dividing this by the value 3.21 of the F-statistic we get
2.420 as the mean square for error, i.e., the estimate for
$\sigma^2.$ However, this is {\em not} the value that we would
have gotten, would we have used the underlying observations, since
the cell means, which are used in the above computation, are given
in very low precision.

In an ANOVA of the upper half of Table 1 in ibid. the two way
interaction terms are tested by an F-statistic with value
$F(2,164)=14.28.$ By Table 4.3.1 of Scheff\'e (1959) a similar
computation as above yields 1.095 as the value of the mean square
for error, based on 170 observations.

Note that a value like 2.3 for a cell mean implies that the actual
value of the cell mean lies in the interval $[2.25, 2.35).$ Using
this rounding off property we may conclude that the first three
F-values given on page 442 of ibid., which have 1 and 164 degrees
of freedom, imply that the value of the mean square for error,
based on 170 observations, lies in the interval $[0.918, 1.218].$
Note that 1.095 belongs to this interval. Averaging the values of
the mean square for error that we get from the last four F-values,
we obtain 1.047 as our estimate.

The F-values presented on page 442 of ibid. that are based on the
second half of Table 1 of ibid., namely $F(2,162)=11.49$ and $F(1,
162)=23.00,$ yield 1.223 and 1.217 as value of the mean square for
error, based on the remaining 168 observations. Averaging yields
1.220. Pooling 1.047 and 1.220 we obtain $${\hat \sigma}_n^2 =
1.134 $$ as our final estimate for $\sigma^2.$ Note that this
deviates considerably from the value 2.42, which has been obtained
from the F-value 3.21 for three-way interaction. Let us presume
here that this is a misprint and that this F-value should have
been something like 6.9.

In order to take care of the rounding off of the values of the
cell means, we have adapted Table 1 in a direction that increases
the double sum in (\ref{jointdensity2}) as much as possible and
that should decrease the evidential value. The resulting table is
given below.

\medskip
\begin{tabular}{|l|l|l|l|}\hline
 Prime type              & Positive & Negative & Irrelevant \\\hline
 Impersonal / Memory     & 2.25     & 3.55     & 2.85       \\
 Impersonal / Impression & 3.35     & 2.55     & 2.85       \\
 Personal / Memory       & 3.25     & 2.25     & 2.85       \\
 Personal / Impression   & 3.55     & 2.55     & 3.05       \\\hline
\end{tabular}
\medskip

Analyzing the same F-statistics as above and performing the same
computations we see that the F-statistics for interactions yield
exactly the same values for the mean square for error. Only the
three F-statistics of the type $F(1,164)$ yield different values.
Averaging the four values for the mean square for error that we
get out of the four F-values related to the upper half of the
table, we arrive at 1.117. The F-statistics for the second half of
the table yield the same estimate 1.220. Pooling 1.117 and 1.220
we obtain $${\hat \sigma}_n^2 = 1.168 $$ as our final estimate for
$\sigma^2$ based on our version of Table 1 of ibid.

\section{Computing Evidential Value}\label{CEV}

Let us group the cells of the tables in the preceding section into
three groups, namely the groups corresponding to the covariate
Prime type with the first two cells in the row Impersonal / Memory
interchanged; $I=12, K=3, I_1=I_2=I_3=4.$ According to the social
psychology theory as put forward in Stapel, Koomen and Van der
Pligt (1996), the participants within these groups should have
similar mean scores. By (\ref{evidentialvalue2}) through
(\ref{evidentialvalue5}) we may compute the evidential value
$\mathbb V$ in favor of the hypothesis $H_p$ that these data have
been fabricated in some way resulting in (\ref{multinormal}) with
$0<\rho<1.$ For the first table from the preceding section, i.e.,
Table 1 from ibid., this yields
$$\mathbb V = 56.88$$
and for the second, adapted table from the preceding section this
yields $$\mathbb V = 1.92.$$

\section{Interpreting Evidential Value}\label{IEV}

With the evidential value $\mathbb V$ defined as in
(\ref{evidentialvalue2}) through (\ref{evidentialvalue5}) the
Bayesian paradigm for criminal court cases
(\ref{BayesParadigmdensity}) becomes
\begin{equation}\label{BayesParadigmV}
 \underbrace{\frac{P(H_p)}{P(H_d)}}_{prior\, odds}\ \
\overbrace{\mathbb V}^{evidential\, value}\, =\
\underbrace{\frac{P(H_p\,|\,E)}{P(H_d\,|\,E)}}_{posterior\,
odds}\,.
\end{equation}
An important principle in criminal court cases is `in dubio pro
reo', which means that in case of doubt the accused is favored. In
science one might argue that the leading principle should be `in
dubio pro scientia', which should mean that in case of doubt a
publication should be withdrawn. Within the framework of this
paper this would imply that if the posterior odds in favor of
hypothesis $H_p$ of fabrication equal at least 1, then the
conclusion should be that $H_p$ is true. So an ANOVA study for
which
\begin{equation}\label{BayesParadigmV2}
 \underbrace{\frac{P(H_p)}{P(H_d)}}_{prior\, odds}\ \
\overbrace{\mathbb V}^{evidential\, value}\, =\
\underbrace{\frac{P(H_p\,|\,E)}{P(H_d\,|\,E)}}_{posterior\, odds}
> 1
\end{equation}
holds, should be rejected and disqualified scientifically.

We conclude with some notes.
\begin{itemize}
\item
ANOVA studies are based on the assumption of normality. Often this
assumption is not satisfied, but the technique is still applied.
This is the case in Stapel et al. (1996), since in Table 1 of
ibid. the measurements are averages of two 7 point Likert scales,
which hardly behave like normal random variables. However, in view
of the central limit theorem cell means like in our basic model
(\ref{multinormal}) behave approximately like (jointly
multivariate) normal random variables.
\item
Note that (\ref{evidentialvalue5}) implies
$$\mathbb V \geq 1.$$
Consequently, within this framework there does not exist
exculpatory evidence. This is reasonable since bad science cannot
be compensated by very good science. It should be very good
anyway.
\item
When a paper contains more than one study based on independent
data, then the evidential values of both studies can and may be
combined into an overall evidential value by multiplication in
order to determine the validity of the whole paper; see the
preceding item.
\item
One may wonder if the way in which the mean square error
(\ref{MSE}) is retrieved from the values of F-statistics,
interferes with the randomness in (\ref{evidentialvalue1}). As
mentioned in Section \ref{EVFD} standard ANOVA theory shows that
this estimator is independent of the exponent in
(\ref{jointdensity1}) and hence (\ref{evidentialvalue1}), provided
the underlying data have a normal distribution; see also item 1.
\end{itemize}

\section{Evidential Value for Fabrication of Data Underlying an ANOVA Study Based on an Alternative Dependence Structure}\label{EVFDA}

In this Section we present an analysis as in Sections \ref{MFD}
and \ref{EVFD}, but under a different dependence structure. Given
the group structure of the cells as presented in the first
paragraph of Section \ref{EVFD} we assume the existence of
$\rho_1,\dots,\rho_K \in [0,1]$ such that
\begin{equation}\label{dependencea} \rho(\varepsilon_{ij},
\varepsilon_{hj})= \rho_k,\quad j=1,\dots,n,\, L_{k-1}+1\leq i
\neq h \leq L_k,\quad k=1,\dots,K,
\end{equation}
hold with all other correlations being equal to 0. This implies
independence between different groups of cells. We note that
(\ref{dependencea}) is just another possible way to model
dependence, and we note again that the actual way in which
fabrication has been implemented, might lead to quite different
dependence structures.

We reconsider the ANOVA study presented via the sample averages of
the cells and the values of some F-statistics. Again the
underlying data themselves are not published and are not
available, and the conclusion of this study is given by
(\ref{conclusion}). There are two hypotheses to be formulated
about the data underlying the ANOVA study. The hypothesis $H_p$ of
fabrication of the data underlying the results, is that at least
one of the $\rho_k$s is positive. The other hypothesis $H_d$
represents the situation that data have been collected according
to (\ref{model}) with independent $X_{ij},$ i.e.,
$\rho_1=\dots=\rho_K=0.$ We want to determine the evidential value
of the ANOVA study, i.e., of the sample means of the cells and the
published F-statistics, with respect to these hypotheses $H_p$ and
$H_d.$ Here the evidential value is defined analogously to
(\ref{evidentialvalue1}) with the supremum taken over $0<\rho_k
<1,\, k=1,\dots,K.$

The sample averages in the cells, $X_{i\cdot}$ from
(\ref{sampleaverage}), have a joint $I$-dimensional multivariate
normal distribution with
\begin{equation}\label{multinormala}
\begin{pmatrix}
X_{L_{k-1}+1\,\cdot} \cr \cdot \cr \cdot \cr X_{L_k\,\cdot}\cr
\end{pmatrix}
\sim \cal N \left(
\begin{pmatrix}
\nu_k \cr \cdot \cr \cdot \cr \nu_k\cr
\end{pmatrix}
,\, \mbox{$\sigma^2 n^{-1}$}
\begin{pmatrix} 1& \rho_k & \cdot & \rho_k \cr \rho_k & 1 & \cdot & \rho_k \cr
\cdot & \cdot & \cdot & \cdot \cr \rho_k & \rho_k & \cdot & 1 \cr
\end{pmatrix}
\right)
\end{equation}
for each $k=1, \dots, K$ and with independence between groups with
different indices $k.$ This entails that the joint density of
$X_{1\cdot},\dots,X_{I\cdot}$ equals
\begin{eqnarray}\label{jointdensity1a}
\lefteqn{\left(\frac n{2\pi \sigma^2}\right)^{I/2} \prod_{k=1}^K
\left[(1+(I_k-1)\rho_k)(1-\rho_k)^{I_k-1}\right]^{-1/2}} \nonumber
\\
 && \exp\left(-\frac n{2\sigma^2} \sum_{k=1}^K\left[
 \frac 1{1-\rho_k}\sum_{i=L_{k-1}+1}^{L_k}
 \left(X_{i\cdot}-\nu_k\right)^2
 \right. \right.  \\
&&  \qquad \qquad \qquad \left. \left.
  -\, \frac {\rho_k} {(1+(I_k-1)\rho_k)(1-\rho_k)} \left(
 \sum_{i=L_{k-1}+1}^{L_k}\left(X_{i\cdot}-\nu_k\right)\right)^2
\right]\right). \nonumber
\end{eqnarray}
This density depends on the parameters $\rho_1,\dots, \rho_K,
\sigma^2, \nu_1,\dots, \nu_K.$ Again, we write ${\hat
\sigma}_n^2,$ for the estimate of $\sigma^2.$

With the notation
\begin{eqnarray}\label{notationa}
\lefteqn{\chi_{n,k}(\rho) =
\left[(1+(I_k-1)\rho)(1-\rho)^{I_k-1}\right]^{-1/2}} \nonumber
\\
 &&  \qquad \quad \exp\left(-\frac {n \rho}{2{\hat \sigma}_n^2(1-\rho)}\
 \sum_{i=L_{k-1}+1}^{L_k}\left(X_{i\cdot}-{\bar X}_k\right)^2
 \right), \nonumber \\
S_{n,k} & = & \frac n{I_k{\hat \sigma}_n^2}\
\sum_{i=L_{k-1}+1}^{L_k }\left(X_{i\cdot}-{\bar X}_k\right)^2,
\\
{\hat \rho}_{n,k} & = & \frac 12 (1-S_{n,k})\left[1+\sqrt{1-\frac
{4S_{n,k}}{(I_k-1)(1-S_{n,k})^2}}\right]\,{\bf 1}_{[S_{n,k} <
(\sqrt{I_k}-1)/(\sqrt{I_k}+1)]}  \nonumber
\end{eqnarray}
Proposition \ref{computationV} of the Appendix shows
\begin{equation}
\mathbb V = \prod_{k=1}^K \chi_{n,k}\left({\hat
\rho}_{n,k}\right).
\end{equation}

Computing this evidential value for Table 1 in Stapel, Koomen and
Van der Pligt (1996), i.e., for the first table of Section
\ref{EFS}, we obtain $$\mathbb V = 14.49.$$ The adapted table,
namely the second table of Section \ref{EFS}, yields
$$\mathbb V = 1.28.$$

Here and in (\ref{evidentialvalue1}) we have defined the
evidential value in the presence of the nuisance parameters
$\nu_1,\dots, \nu_K$ by replacing these parameters by their
maximum likelihood estimators. An alternative approach is to
compute the evidential value keeping these parameters fixed, and
to subsequently minimize the resulting evidential value over these
nuisance parameters; in formula
\begin{equation}\label{evidentialvalue6}
\tilde {\mathbb V} = \inf_{\nu_1,\dots,\nu_K \in \mathbb R}
\frac{\sup_{0<\rho_k<1, k=1,\dots,K}f_n(X_{1\cdot},\dots,
X_{I\cdot}; \nu_1,\dots,\nu_K,\rho_1,\dots,\rho_K)}
{f_n(X_{1\cdot},\dots, X_{I\cdot}; \nu_1,\dots,\nu_K, 0,\dots,0)},
\end{equation}
where $f_n(X_{1\cdot},\dots, X_{I\cdot};
\nu_1,\dots,\nu_K,\rho_1,\dots,\rho_K)$ is the density as given in
(\ref{jointdensity1a}) with $\sigma$ replaced by ${\hat
\sigma}_n.$ In fact, both definitions of evidential value yield
the same value in the situation of this Section \ref{EVFDA}, as is
shown in Theorem \ref{V=V}.