\documentclass[doublecol]{epl2} 
% or \documentclass[page-classic]{epl2} for one column style
%\documentstyle[a4paper,12pt]{article}
%\documentclass[12pt]{article}

%%\usepackage{jas}
%\usepackage{a4}
%
\usepackage{psfig,epsfig} 


%%\usepackage{amsfonts}
\usepackage{amsmath}
%\usepackage{graphicx}
%\usepackage{latexsym}
\usepackage{bm}
%% bold math
%\usepackage{ulem}
%\usepackage{a4wide,amsmath,epsfig,fancyheadings,times,eepic,psfig,pstricks}
\usepackage{amssymb}
%
%\usepackage{footnpag}
%\usepackage{rotating}
%\usepackage{euscript}
%\usepackage{tabularx}
%%\usepackage{minitoc}

%\setcounter{secnumdepth}{2}
%\setcounter{tocdepth}{1}

%\tolerance=8000
%\pretolerance=4000

%\newenvironment{listdot}{\begin{list}{$\bullet$}{\setlength{\itemsep}
%{0ex} \setlength{\parsep}{0ex} \leftmargin 1em}}{ \end{list}}
%\newenvironment{listcar}{\begin{list}{$\square$}{\setlength{\itemsep}
%{0ex} \setlength{\parsep}{0ex} \leftmargin 1em}}{ \end{list}}

%\oddsidemargin 0cm

%\textwidth 17 cm
%\textheight 23 cm
%\topmargin 0cm
%\usepackage{endf}

\newcommand{\scal}[2]{{< #1 \mid #2>}}
\newcommand{\ext}[2]{{\mid\, #1 >< #2\mid}}
\newcommand{\bra}[1]{<\!{#1}\!\mid}
\newcommand{\ket}[1]{\mid\!{#1}\!>}
\newcommand{\Deta}{{\mathbf{\delta} \vect{\eta }} }
\newcommand{\Dfi}{{\delta \varphi } }
%%\newcommand{\vv}[1]{\boldsymbol{#1}}

%\renewcommand{\thesubsection}{\alph{subsection}}

% met lews ref entre () in lieu of []
%\usepackage[round]{natbib}

% bib refs as [1] etc
\bibliographystyle{unsrt}
%\bibliographystyle{jas}


%\draft
\title{Feedback characteristics of nonlinear dynamical systems}

\author{A. Lahellec\inst{1}\thanks{E-mail: \email{alain@lmd.jussieu.fr}}
 \and S.Hallegatte\inst{2} \and
J-Y. Grandpeix\inst{1} \and P.~Dumas\inst{3} \and S.~Blanco\inst{4}}
\shortauthor{A. Lahellec \etal}

\institute{
 \inst{1}Laboratoire de M\'et\'erologie Dynamique, CNRS, Universit\'e~Paris-VI\\
 \inst{2}Centre International de Recherche sur l'Environnement et le
D\'eveloppement, ENPC, Paris\\
 \inst{3}Environment Research and Teaching Institute, Ecole Normale Sup\'erieure, Paris\\
 \inst{4}Laboratoire LAPLACE, Universit\'e~Paul-Sabatier, Toulouse
}

\pacs{05.45.-a}{Nonlinear dynamics and chaos }
\pacs{92.70.Gt}{Climate dynamics}

%\renewcommand{\baselinestretch}{3}
%\normalsize
%\setlength{\topmargin}{.3 in}

\abstract{
We propose a method to extend the concept of feedback
gain to nonlinear models. The method is designed to dynamically characterise a feedback
mechanism along system natural trajectory.
\revision{ The numerical efficiency of the method is proved using the~Lorenz~(1963)
classical model. Finally, a simple climate model of water-vapour feedback
shows how nonlinearity impacts feedback intensity along the seasonal cycle}.}
%
%\bigskip
%\small
%\begin{tabular}{l}
%\hline
%{\it Corresponding author's address}:\\
%Alain Lahellec, \\
%Laboratoire de M\'et\'eorologie Dynamique,\\
%UNIVERSITE Pierre et Marie Curie \\
%Tour 55 - 45 / 3\`eme etage - Boite 99\\
%4 place Jussieu\\
% F-75252 PARIS CEDEX 05\\
% France\\
%E-mail: alain@lmd.jussieu.fr  fax:(33) 1 44 27 62 72
%\end{tabular}
%\large
%\normalsize
\begin{document}
%
\maketitle
%
\section{Introduction}
%---------------------------------------------------------------
Since its introduction by Bode~\cite{Bode45}, the concept of feedback gain has
proved a very efficient tool in the engineering fields.
Originally elaborated in the frequency domain for
linear systems, it further allowed for the boom in techniques of optimal
control that have indeed been extended to nonlinear systems in the 60s --
robust filtering, optimal or adaptive control, etc (e.g.~\cite{gs96}, and the references
therein). It also inspired a variety of research applications for controlling or
synchronising chaotic systems using feedback laws~(see~\cite{ott02}).

Our concern here is different: the concept of feedback is thought
as a way to gain insight into the behaviour of a system.
In this vein, feedback is frequently mentioned
in biology, in chemistry, in physics, and in economics~(e.g.~the Keynes feedback, in~\cite{Chiarella}),
but usually on a qualitative basis only: one will essentially speak
of positive (destabilising) or negative (stabilising) effect of a feedback mechanism.
Manabe and Wetherald~\cite{manabe67} originated a breakthrough in climate sciences
when they proposed an analysis of the Earth climate
allowing to assess in terms of quantitative feedback gain the relative importance of
various contributing phenomena.
%
From then on, the study of global warming, caused by the
increasing concentration of greenhouse gases like \chem{CO_2}, has been systematically analysed
in terms of feedback processes (water vapour, clouds, albedo ...), 
as in~\cite{Bony06,Cess89,Hansen84,Schneider72} among many others. Nevertheless, there
is still ambiguity about a general definition of feedback gain, as analysed
in~\cite{Bates07}, where four definitions are found.

Considering the role of water vapour feedback (WVF) in global warming,
the direct climate response to an increase in \chem{CO_2} concentration is an increase
in the mean atmospheric temperature ($T_{\mbox{\tiny{atm}}}$), due to change in long wave
radiation budget. Following this direct effect, the feedback process
runs as follows: since climatic mechanisms maintain relative humidity unchanged, the
atmospheric  warming translates into a larger water vapour content in the atmosphere.
Water vapour being a strong greenhouse gas, the initial $T_{\mbox{\tiny{atm}}}$
perturbation is amplified.
From the comparison between the equilibrium response obtained in climate models when
including that feedback and the similar response when artificially suppressing that
feedback, it is commonly said that ``WVF increases by 60 percent the
warming due to carbon dioxide perturbation''.

A major limitation of this approach is that it ignores the dynamical aspects of Bode's original
concept since only the difference between two systems at equilibrium is considered.
Indeed, for systems that take very long to reach a new
equilibrium -- as the climate system -- the dynamics are of primary importance.
% 
In~\cite{HLG05}, the authors explained how to retrieve the dynamics of feedback for systems between
two equilibrium. Regarding the \chem{CO_2} increase,
they showed with an idealised model that the atmospheric WVF response is negative on the short term
(damping effect), and only becomes positive (amplifying effect) after a build-up period of four to seven
years, where the 60 percent increase is finally reached at the perturbed equilibrium. 

This equilibrium approach is nevertheless a severe limitation. Indeed, climate is driven by natural
forcing (e.g. solar forcing) and can hardly be considered as
in equilibrium. The purpose of this article is to extend the
feedback concept to the dynamical nonlinear system.

We shall first introduce a generic method to define a feedback loop in a state-space model.
Then, system response to external perturbation is characterised by a set of dynamical functions
related to feedback gain. We observe that when the system simplifies into autonomous
conditions, these functions are the same as those introduced by Bode.
% new
\revision{ A first application on the Lorenz~(1963) model is carried out to demonstrate the numerical
efficiency of the method on a strongly nonlinear model. Finally, we come back to
the water-vapour feedback and, using a very simple nonlinear model,
we compare the equilibrium response with more realistic trajectories}.
%
\subsection{Formal background}
Consider a state-space model driven by some external forcing~$\vect{h}(t)$:
%
\begin{equation}
\partial_t \vect{\eta} (t) =  \vect{g}(\vect{\eta}(t),\vect{h}(t));\qquad \vect{\eta}(t=0)=\vect{\eta}_{o},  \\
\label{state_ini}\end{equation}
%
where $\vect{\eta}$ is the state vector. Once a reference trajectory~($\mathcal{T}^{\mbox{\tiny{ref}}}$) is followed, perturbations
to system~(\ref{state_ini}) are solutions of the so-called Tangent Linear System (TLS)~:
%
\begin{equation}
\partial_t \Deta(t) =  A(t)\, \Deta(t) + B(t)\, \delta\vect{h}(t),\quad t\in[0,T];
\label{state}
\end{equation}
%
In this equation, $\Deta$ is the deviation of the perturbed system from the reference
trajectory, caused by a perturbation $\delta \vect{h}(t)$ of the reference forcing function,
and $A(t)=\partial_\eta\vect{g}(t)$ and~$B(t)=\partial_h\vect{g}(t)$ are the Jacobian matrices.
We call Circulating Tangent Linear System~(CTLS) the non-autonomous TLS~(\ref{state}),
meaning that it is defined along~$\mathcal{T}^{\mbox{\tiny{ref}}}$.
The classical solution of~(\ref{state})~reads
%
\begin{equation}
\left \lbrace
\begin{array}{l}
\Deta(t) = \Phi (t,0)\Deta(0) + \int_0^{t} {\Phi (t,\tau)B(\tau)\delta\vect{h}(\tau)} d\tau \\
\partial_t \Phi (t,\tau) = A(t) \Phi (t,\tau);\qquad  \Phi (\tau,\tau) = I \, ,
\end{array}
\right .
\label{propag}\end{equation}
%
where $\Phi(t,0)$ is the state-transition matrix of the perturbed system and
$\Deta(0)$ a perturbation of initial conditions~\cite{Luenberger}.
%
\revision{Such a system does not exhibit any explicit structure in which a feedback loop can be
specified. We now introduce a way to provide a state model with a feedback structure.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
\section{Feedback concepts}
%---------------------------------------------------------------
\subsection{Feedback loop}
A feedback-loop is a circular structure that links
some chosen physical quantity to the rest of the system. A general method to provide
a model with such a structure can be stated as follows:
%
\begin{enumerate}
\item choose an ``output'' variable from the model that will be used as
a test-variable -- {\it \revision{in the WVF case in this letter,
we choose $T_s$, the atmospheric surface temperature, as the test-variable}};
\item add a perturbed version of this variable to the system -- {\it say} $\varphi=T_s+ \textsl{perturbation}$;
\item close the feedback-loop: the new variable~$\varphi$ is introduced back in the model
where the mechanisms are depending on its perturbed version
-- \revision{{\it in equation~(\ref{Qequ}), the water content is chosen as depending on $\varphi$ instead of $T_s$};}
\end{enumerate}
%
Accordingly, the introduction of a feedback-loop in the original system corresponds to the following
formal procedure:
%
\begin{itemize} 
\item a scalar variable $\varphi$ is added to the system with equation
$\varphi=f(\eta)$ -- see eq.~(\ref{Lorenz}) for instance;
\item the primitive model~(\ref{state_ini}) is modified to be sensitive to $\varphi$
in such a way that original
and modified models are mathematically equivalent:\\
$\partial_t \vect{\eta} (t) =
  \vect{g'}(\vect{\eta}(t),\varphi(t),\vect{h}(t))=\vect{g}(\vect{\eta}(t),\vect{h}(t))$;
\end{itemize}
%
Now, if a scalar perturbation $u(t)$ is applied to the $\varphi$ equation only
($\varphi'(t)= f(\eta) + u(t) $), the system will respond to that change
\revision{within a feedback loop: perturbation $u$ on $\varphi$ 
$\rightarrow$ selected mechanisms $\rightarrow$
rest~of~the~system $\rightarrow$ response $\delta\varphi$~(see fig.\ref{fbloop}).}
%
\begin{figure}[ht]
\centerline{\psfig{figure=other/deuxg.eps,%
bbllx=0pt,bblly=0pt,bburx=391pt,bbury=160pt,width=8cm,height=3.3cm,clip=}}
\caption{Feedback structure using an extra-variable $\varphi$ (left panel) and
the classical block-diagram (right panel).}
\label{fbloop}
\end{figure}
%

\revision{The opening of the loop means preventing
the system to respond on $\varphi$ (cross in fig.~\ref{fbloop}): $\delta\varphi^o=u$.
The classical response would be $\delta\varphi=\frac{1}{1-g}u$ in a symbolic form.}

The difference between perturbation $u(t)$ and response in $\varphi(t)$
is commonly understood as the {\it feedback effect} of the system on $\varphi$.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
\subsection{A general definition of dynamical feedback}
% -------------------------------------------------------
Following the preceding procedure,
the original tangent system~(\ref{state}) is modified into:
%
\begin{equation}
\left \lbrace
\begin{array}{l}
\partial_t \Deta(t) = A^{\flat}(t)\, \Deta(t)\, + \ket{b(t)} \Dfi (t) \\
\Dfi (t) = \bra{c(t)} \Deta(t) + u(t)\, .  
\end{array}
\right .
\label{basemod}\end{equation}
%
 The column and row matrices
$|\,b\!\!>=\partial_{\varphi}\vect{g'}$ and $<\!\!c\,|=\partial_{\eta}f$\\ formalise the structure of
the feedback loop.
For system~(\ref{basemod}) to be mathematically equivalent to~(\ref{state}), the following
equality $A(t) = A^{\flat}(t) + \ext{b}{c}(t)$ must hold. Opening the loop corresponds to the
cancelling of $\bra{c}$, that is disallowing the system perturbation to influence $\varphi$.

It is easily checked that system~(\ref{basemod}) takes the form of~(\ref{state})
when $\Dfi$ has been eliminated. Correspondingly, we define the state sensitivity
$\vect{s}(t,0)$ as the state response in~(\ref{basemod})
using the solution~(\ref{propag}) -- setting $\Deta(0)=0$:
%
\begin{equation}
\vect{s}(t,0) =\Deta(t) =  \int_{0}^{t} {\Phi (t,\tau)\ket{b(\tau)} u(\tau)}\, d\tau\, .
\label{repeta}
\end{equation}
%
We call~$\sigma^u(t,0)=\Dfi(t)$ the response in $\varphi$,
and $\varrho^u(t,0) = \sigma^u(t,0)-u(t)$ the {\it feedback effect} of system on $\varphi$,
which is the deviation of the response from the perturbing function starting at $\tau=0$:
%
\begin{equation}
\varrho^u(t,0) = \,\bra{c(t)}\int_{0}^{t} {\Phi (t,\tau)\ket{b(\tau)} u(\tau)}\, d\tau\, .
\label{repg}
\end{equation}
%
This function can be computed along $\mathcal{T}^{\mbox{\tiny{ref}}}$ from the system obtained
by time-differentiating eq.~(\ref{repeta}):
%
\begin{equation}
\left \lbrace
\begin{array}{l}
\partial_t \vect{s}(t,0) = A^{\flat} \vect{s}(t,0)  + \ket{b} \varrho^u(t,0)+\ket{b}u(t)\\
\varrho^u(t,0)  =  \,\bra{c} \vect{s}(t,0)   ;\quad \vect{s}(0,0)=\vect{0}
\end{array}
\right.
\label{ptsig}
\end{equation}
%

In \cite{Cacuci81a}, Cacuci has shown how to obtain such systems
using G\^ateaux differentiation, and called ``general sensitivity''
the response $\Dfi$ to $u$.

It is worth noticing that the CTLS integration of~(\ref{ptsig}) is carried out simultaneously 
with the reference trajectory, that is, {\it without} applying any actual
numerical perturbation that would change the trajectory. This is admissible as far as perturbation
amplitude stays within the limits of linearity of the perturbed system. Also, the CTLS provides
a \revision{definite} advantage concerning the generic features of the results, compared to finite difference
methods using finite perturbation amplitudes.
%
%...............................
\subsection{The autonomous limit}
%--------------------------------
On short segments of trajectory, the constant Jacobian
matrices define a local TLS, with a simpler state transition
matrix: $\Phi^s(t,\tau) = \exp A (t-\tau)$. For such autonomous system, a Laplace transformation
can be applied to the convolution product in~(\ref{repg}) to obtain:
%
\begin{equation}
\tilde{\varrho}^u _s(\mu) = \,\bra{c}(\mu I - A)^{-1}\ket{b} \tilde{u}(\mu)\, ,
\label{repeflap}
\end{equation}
%
where $\tilde{f}(\mu)$ is the Laplace transform of $f(t)$, with $\mu$
as the Laplace variable. We call {\it standing} this function corresponding to the dynamical
response of the local autonomous system to perturbation. When the system reaches an equilibrium,
the CTLS becomes truly autonomous, and \revision{the} standing response identifies with
the circulating one, a property that allowed the dynamical analysis in~\cite{HLG05}
simply using the inversion of eq.~(\ref{repeflap}).

To retrieve the feedback gain, one can directly apply a Laplace transformation to
the autonomous version of eq.~(\ref{basemod}) to get
%
\begin{equation}
\left \lbrace
\begin{array}{l}
(\mu I-A^{\flat}) \tilde{\vect{s}}(\mu) =\ket{b}\tilde{\sigma}(\mu) \\
\tilde{\sigma}(\mu)=\,\bra{c}\tilde{\vect{s}}(\mu)  + \tilde{u}(\mu) \, ,
\end{array}
\right .
\end{equation}
%
which, after elimination of $ \tilde{\vect{s}}$, leads to:
%
\begin{equation}
 ( 1-\bra{c}(\mu I-A^{\flat})^{-1}\ket{b})\tilde{\sigma}(\mu)=\tilde{u}(\mu)\, ,\\
\label{sltbor}\end{equation}
%
The form $(1-\tilde{g}(\mu))\tilde{\sigma}(\mu)=\tilde{u}(\mu)$ of eq.~(\ref{sltbor}) is the one
used by Bode as the definition of feedback gain in the Laplace domain:
\begin{equation}
\tilde{g}(\mu)=\,\bra{c}(\mu I-A^{\flat})^{-1}\ket{b}
\label{sltborg}\end{equation}
%
In this form, it is easily seen that setting $\bra{c}=0$ corresponds to the open-loop mode,
that is, to a null feedback gain.

\bigskip
We may conclude that the two feedback functions $\sigma$ and $\varrho$ are
generalising the \revision{autonomous} ones with the  following correspondence:
%
\begin{equation}
\left \lbrace
\begin{array}{l}
\sigma^u(t,0)\longleftrightarrow\tilde{\sigma}(\mu)=\frac{1}{1-\tilde{g}(\mu)}\tilde{u}(\mu)\\
\varrho^u(t,0)\longleftrightarrow\tilde{\varrho}(\mu)=\tilde{u}(\mu)-\tilde{\sigma}(\mu)=\frac{\tilde{g}(\mu)}{1-\tilde{g}(\mu)}\tilde{u}(\mu)\end{array}
\right .
\label{coresp}\end{equation}
%
%
%...............................................................
\subsection{Dynamical feedback characteristics}
%...............................................................
The standard definition of feedback characteristic functions $\sigma^{\Upsilon}$ and
$\varrho^{\Upsilon}$ considers the particular
perturbing step-function $u(\tau)=\Upsilon(\tau)$, whose Laplace transform
is $\tilde{\Upsilon}(\mu)=\frac{1}{\mu}$ in~(\ref{coresp}). Correspondingly, the
feedback-effect characteristic function $\varrho^{\Upsilon}$ is given by
eq.~(\ref{repg}) setting $u(\tau)=1$. This form is referred to as the
``F4 sensitivity-altering'' feedback in~\cite{Bates07}.

Unlike the linear feedback mechanism which is completely characterised by autonomous functions,
the nonlinear case is depending both on a trajectory and a starting point.
It may happen that, in practical applications, those primitives are known:
in an economic-model, a shock (i.e. on oil prices) arises at a known date; in climate models,
a volcano eruption is also chronologised.
In general the feedback response may strongly depend on the initial
time of perturbation: this pertains to the richness of nonlinear systems behaviour.
A comparison between Standing-(\ref{repeflap}) and Circulating-(\ref{repg}) TLS-functions~ is
presented in the application section.

As a final formal remark, the Dirac impulse $u(\tau)= \delta(\tau)$ may also be useful
as a perturbing function because
the corresponding $\varrho^{\delta}(t,\tau)$ kernel can be convoluted with any
perturbing function $v(t)$ to produce the model response:
$\delta\varphi(t) = \bra{c(t)}\int_0 ^t\varrho^{\delta}(t,\tau) v(\tau) d\tau+v(t)$. 
%
Because in that case $\varrho^{\delta}(t,\tau)$ is needed as
a function of $\tau$, the feedback system~(\ref{ptsig}) shall be integrated in its adjoint-system
form, see~\cite{Cacuci84}. That feedback form corresponds to the ``F3 stability-altering''
class in~\cite{Bates07}.
%
\subsection{Circulating feedback-gain function}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Contrary to the linear case, one cannot retrieve the feedback gain of nonlinear systems
from neither $\varrho$ nor $\sigma$, as in~(\ref{coresp}). A particular system has to be built.
Let us consider the following system:
%
\begin{equation}
\left \lbrace
\begin{array}{ll}
\partial_t \vect{s}(t,0) = A^{\flat} \vect{s}(t,0)  + \ket{b} \Upsilon(t)  \\
g(t,0)  =  \bra{c} \vect{s}(t,0) \, ,
\end{array}
\right .
\label{psig}\end{equation}
%
In the autonomous case, that system simplifies into:
$\tilde{g}(\mu) = \mu^{-1}\bra{c}(\mu I -A^{\flat})^{-1}\ket{b}$,
that is recognised as Bode's feedback gain.
System~(\ref{psig}) is the so-called ``open-loop'' mode which provides for the general
feedback-gain function in its circulating form.
%% %end ancien
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Application to the Lorenz model}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We illustrate the method using the Lorenz~\cite{Lorenz63} model
to prove the numerical efficiency of the method on a well known nonlinear system.
%\revision
%In particular, standing and circulating
%responses are compared to exhibit their fundamental difference.
%
The system -- modified into feedback form -- reads
%
\begin{equation}
\quad
\left \lbrace
\begin{array}{l}
\partial_t x = s  (y - x)\\
\partial_t y = (r  x - y - xz)\\
\partial_t z\ = -b z + \varphi \\
\varphi = xy   \, ,
\end{array}
\right .
\label{Lorenz}
\end{equation}
%
with the set of parameters $s=10, r=28,b=\frac{8}{3}$ that leads to
a chaotic behaviour~(Figure~\ref{lor4} displays the $x,y$ trajectory
with the two wings of the attractor).

%
\begin{figure}[ht]
\centerline{\epsfig{figure=PS/xyrepcolokn.eps,%
bbllx=49pt,bblly=49pt,bburx=409pt,bbury=301pt,width=8cm,height=5cm,clip=}}
%bbllx=49pt,bblly=49pt,bburx=409pt,bbury=301pt,width=13.5cm,height=6.25cm,clip=}}
\caption{Test variable in the Lorenz model and feedback effect CTLS-function.}
\label{lor1}
\end{figure}
%
To build the feedback form~(\ref{Lorenz}), we have chosen
$\varphi=xy$ as a nonlinear test variable, but any other term could have been chosen,
since we are not interested in physical analysis here\footnote{\revision{nevertheless, a physical interpretation
can be found in {\it http://www.lmd.jussieu.fr/ZOOM/doc/examples.html}}}. Also, the comments on
dynamics will only consider numerical aspects of the results.
Figure~\ref{lor1} shows the evolution of $\varphi$ and CTLS~$\varrho^{\Upsilon}(t,0)$ function
along the trajectory, starting at the beginning of an
orbital point~A and going after four revolutions through a transition point~B (in the following,
$\Upsilon$ is dropped from symbol $\varrho^{\Upsilon}$).

%
\begin{figure}[hb!]
\centerline{\epsfig{figure=PS/ctlscolok.eps,%
bbllx=50pt,bblly=50pt,bburx=408pt,bbury=300pt,width=8cm,height=4.7cm,clip=}}
%bbllx=50pt,bblly=50pt,bburx=408pt,bbury=300pt,width=13cm,height=6.5cm,clip=}}
\caption{Long term feedback effect~(log scale). \revision{The $x$ trajectory is
rotated along the Lyapunov exponent:~$x\rightarrow x+\lambda_L t$, with~$\lambda_L \simeq 0.95$,
to ease curve collocation}.}
\label{loa2}
\end{figure}
%
Because of the symmetry of the system across $z$, the test variable $\varphi=xy$ does not
clearly reveal the transition, whereas the CTLS feedback-effect does.
Between points A and B on the figure, the interaction
mode responds with a zigzag at each orbit, increasing regularly in amplitude.
After B the response at transition amplifies.

This systematic amplification before transitions is seen on 
fig.~\ref{loa2} which displays the long-term tendency of $\varrho(t,0)$.
%
The chaotic nature of the model is correctly characterised by the average linear increase
of the log-CTLS-effect which is found to follow a local Lyapunov exponent.
Moreover, the oscillation is similarly increasing in amplitude, an aspect that is not revealed
by the Lyapunov analysis.

%
\begin{figure}[ht]
\centerline{\epsfig{figure=PS/bigctlscolb.eps,%
bbllx=50pt,bblly=50pt,bburx=408pt,bbury=300pt,width=8cm,height=5cm,clip=}}
%bbllx=50pt,bblly=50pt,bburx=408pt,bbury=300pt,width=13.5cm,height=6.75cm,clip=}}
\caption{Same as fig.~\ref{loa2} with initial conditions out of the attractor.}
\label{loabis}
\end{figure}
%
Figure~\ref{loabis} illustrates what happens when starting far from the attractor.
Clearly, the feedback-effect is very different depending on whether the current point
is lying within the attractor or not:
the average increase of $\varrho(t,0)$ reduces each time
the circulating point leaves the attractor (see times $0,~35,~48,~73$).
The attractor can therefore be characterised by an enhanced action of the
$xy$ mode interaction.
%
\subsection{Short-term dynamical response}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The asymptotic behaviour of feedback function shows coherence with Lyapunov stability analysis,
but feedback analyses are further concerned with short and medium-term dynamics.
To account for the Lyapunov divergence of the Lorenz system, we will remove this divergence
by applying a factor of $\exp(-\lambda_{\mbox{\tiny{L}}} t)$ on dynamical functions.

Let now zoom on the four orbits sequence going to transition. In fig.~\ref{lor4}, the
$x,y$ trajectory is plot with error-bars. Error-bars represent the
CTLS-effect amplitude~(log-scaled), plot as blue-vertical when {\it above} the mean Lyapunov
increase, green-horizontal when {\it under}.
%
\begin{figure}[ht]
\centerline{\epsfig{figure=PS/refig5.eps,%
bbllx=60pt,bblly=65pt,bburx=552pt,bbury=750pt,width=5.5cm,height=8cm,angle=-90,clip=}}
%bbllx=49pt,bblly=49pt,bburx=410pt,bbury=300pt,width=13cm,height=8cm,clip=}}
\caption{Lorenz attractor in the $x,y$ plane; error bars:~CTLS-effect (scaled, see text).}
\label{lor4}
\end{figure}
%
Both orbits are followed clockwise. The first (inner) orbit shows the {\it under} effect only,
when the point is orbiting around the left-wing point of steady convection~(at the centre of orbits).

Between points 1 and 2 on the fourth orbit, the {\it above} reaction is active when the movement decelerates
toward perigee (close to central point [0,0]), and when it accelerates heading to the
apogee~(3,4). This suggests that the {\it above} effect
is determined by a repulsion from central point. The transition initiates with a strong {\it under}
reaction after the apogee~(5), and the circulating point decelerates again from the central point
repulsive effect, but with an offset that knocks the state on to the right wing.

\revision{Such dynamical characteristics show how different regions of the attractor
influence the behaviour of the system. The feedback formalism just brings another approach,
simple and general, to the numerous and sophisticated tools developed for nonlinear analysis.}
%
\section{\revision{Water vapour feedback model}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\revision{
Coming back to the water vapour feedback effect, we have built the most simple model
to analyse the slow response of climate to a \chem{CO_2} perturbation that matches the results
obtained in~\cite{HLG05}. The purpose of the exercise is to give an idea of what importance
nonlinear mechanisms might have on WV feedback.

We consider the evolution of a 30~m thick mixed ocean layer, irradiated by a mean solar flux
$F_{SW}=343$~W~m$^{-2}$. The energy loss is from the LW radiative flux only:
%
\begin{equation}
C_o \frac{d T_s}{dt}= F_{SW} - F_{LW}\;,
\label{wv}\end{equation}
%
with $C_o=1.257~10^8~JK^{-1}m^{-2}$.

To simulate the effect of greenhouse gases (WV or \chem{CO_2}) on radiation, we use an equivalent radiation altitude $z_e$.
Seen from space, this altitude gives an equivalent outcoming blackbody radiation temperature due
to both surface and low-level layers of the atmosphere (that contain a substantial amount of water
vapour):
%
\begin{eqnarray}
F_{LW} &=&\sigma(T_s-z_e\gamma_z)^4\;,
\label{FLW}\end{eqnarray}
%
where $\sigma=5.67~10^{-8}W m^{-2} K^{-4}$ is the Stephan-Boltzmann's constant and $\gamma_s=6.5~Kkm^{-1}$
the fixed atmospheric temperature lapse-rate~({\it i.e.}~$\gamma_s=\partial_z T_{\mbox{\tiny{atm}}}$). 
The specific humidity at surface $Q$ is given by the Clausius-Clapeyron relation:
%
\begin{eqnarray}
Q &=& f(T_s) = \mu\;\exp{\displaystyle[\frac{L_v}{R_m}(\frac{1}{T_0}-\frac{1}{T_s})]}\,,
\label{Qequ}
\end{eqnarray}
%
with $\mu=11.05$, $\frac{L_v}{R_m}=5423~K$ and $T_0=273.15~$K.

We now determine the dependence of $z_e$ in $Q$. When the ocean temperature T$_s$ increases due to
a \chem{CO_2} concentration increase, more water is evaporated in the atmosphere, and $z_e$ increases: $z_e= z_0+\alpha Q$.
Parameters $z_0=1$~km and $\alpha=0.02$ have been calibrated from the more complete model in~\cite{HLG05}, as
they are accounting for the elevation of $z_e$ due to the amount of water vapour in low atmosphere, which is
non trivially related to $Q$.

To build the WV feedback loop, we define the new variable as $\varphi=T_s$, and let the atmospheric
water vapour depend on it in eq.~(\ref{Qequ}): $Q = f(\varphi)$. In this new structure of the system,
$\varphi$ is the air-surface temperature, and can now be perturbed independently from $T_s$.

The WV feedback loop is~:
$\varphi\uparrow\; \rightarrow Q\uparrow\;\rightarrow z_e\uparrow\; \rightarrow LW\downarrow\, \rightarrow T_s\uparrow\;\rightarrow
\varphi$ increases more (positive feedback). The open-loop mode is obtained when $Q$ does not respond
to change in $\varphi$.

We first consider a $T_s$ trajectory starting out of equilibrium from 16 $^\circ$C, and slowly reaching
equilibrium at 17~$^\circ$C (red points-curve in fig.~\ref{wvfb}). At the final equilibrium,
the standing response to a perturbation in $\varphi$ is given by the blue curve, an exponential
response that reaches 0.43 with an e-folding time of 423 days.
%
\begin{figure}[ht]
\centerline{\epsfig{figure=PS/wv2repeff.eps,%
bbllx=60pt,bblly=65pt,bburx=560pt,bbury=750pt,width=5.5cm,height=8cm,angle=-90,clip=}}
%bbllx=49pt,bblly=49pt,bburx=410pt,bbury=300pt,width=13cm,height=8cm,clip=}}
\caption{Water vapour feedback model responses.}
\label{wvfb}
\end{figure}
%
The corresponding feedback gain is 30~\% as in \cite{HLG05}~(the model is calibrated for
that purpose). This means that if a step in \chem{CO_2} concentration is responsible for
a final increase of 1~K in the model temperature, then an immediate fraction of~0.7~K
is due to the direct \chem{CO_2} effect, while the WV feedback is gradually adding an extra
0.3~K~({\it i.e. }$0.7\times0.43$).

The response time is nevertheless much shorter than the several years in~\cite{HLG05},
this being due to the absence of delay in the coupling of the ocean to the atmosphere in this crude model,
where the response time is proportional to $C_o$  -- or to the ocean layer
thickness\footnote{the autonomous response can be analytically obtained from the application
of~eq.(\ref{repeflap}) and inversion. One gets the asymptotic feedback gain:
 $g_\infty=\alpha\gamma_z\frac{L_v}{R_m}\frac{Q(T_s)}{T_s^2}$
and an e\_folding time of $\tau=\frac{1}{(1-g_\infty) T^3(z_e)}\frac{C_o}{4\sigma\gamma_z}$
for the response.}.

The CTLS response is the green curve
just above with the same asymptotic response of 0.43.
Following the trajectory, the circulating response is very close to the standing one, meaning
that the nonlinearity of the model is marginally active.

We now consider a cyclic seasonal SW forcing with an added amplitude of $\pm$100~W~m$^{-2}$; the response
is the blue oscillating curve and the thin red line is the surface temperature trajectory, finally
oscillating between 9~$^\circ$C and 25~$^\circ$C. The nonlinear effect of the model is illustrated by a
stronger WV feedback in the mean.  Moreover, it is even stronger during autumn ($\varrho\simeq 0.52$)
than in spring~(0.44).

If this minimalist model of climate where realistic, it would tell that in response
to a constant \chem{CO_2} step, fall temperatures would be amplified more that spring temperatures.
Such asymmetry proceeds from the nonlinearity of the model.}
%
\section{Final remarks}
%%%%%%%%%%%%%%%%%%%%%%%%%%%
We have shown how to extend to nonlinear systems the original Bodes'approach of feedback analysis.
Using the \revision{Tangent Linear System} along trajectories gives rise to a natural extension of linear
concepts and allows to exhibit the dynamical aspects of feedback characteristic functions.

In practical applications, the addition of an extra variable \revision{to the basic state model} as we did
is a simple way to unambiguously build a feedback-loop structure in a model,
including explicit reference to the open-loop mode.
Moreover, this new variable offers flexibility in implementing meaningful feedback loops.
In gridded models for instance, a meaningful response relies in the
average of an ensemble of grid variables, because the space discretisation is guided
by numerical rather than physical concern. A combination of
state-space-averaging methods with feedback characterisation would bring important gain
in gridded-systems analysis.

As conclusion, the proposed CTLS approach offers a technique to apply dynamical feedback concepts
to nonlinear models while reflecting some generic features of the results.
Moreover, the dynamical aspects of feedback mechanisms
address the basic issues in the behaviour of complex systems, especially the effect of interdependence of parts
in the functional complexity of the whole system.

\acknowledgments
The first author wishes to thank J.L.~Dufresne, M.~Ghil,  J.I.~Yano and
two anonymous referees for useful comments, and especially his daughter
Julie Lahellec for the English version.

%
%\bibliography{/home/al1/Cahiers/Gains/Articles/biblio}

\begin{thebibliography}{10}

\bibitem{Bode45}
\Name{Bode H.W.}
\Book{Network Analysis and Feedback Amplifier Design}.
%\newblock {\em Network Analysis and Feedback Amplifier Design}.
\Publ{Van Nostrand, New York}
\Year{1945}

\bibitem{gs96}
\Editor{Gupta M. \and Sinha N.}
\newblock {\em {Intelligent Control Systems:~Theory and Applications}}.
\Publ{IEEE Press}
\Year{1996}


\bibitem{ott02}
\Name{Ott Ed.}
\Book {Chaos in Dynamical Systems}.
\Publ{Cambridge University Press}
\Year{2002}

\bibitem{Chiarella}
\Name{Chiarella C. \and Flaschel P.}
\Book {Keynesian Monetary Growth: Macro Foundations}.
\Publ{Cambridge University Press}
\Year{2000}

\bibitem{manabe67}
\Name{Manabe S. \and Wetherald R.}
%\newblock \protect{Thermal Equilibrium of the Atmosphere with a given
%  Distribution of Relative Humidity}.
\REVIEW{J. Atmos. Sci.}{24}{1967}{241}
%, 24(3):241--259, 1967.

\bibitem{Bony06}
\Name{Bony S. \etal}
%, R.~Colman, V.M. Kattsov, R.P. Allan, C.S. Bretherton, J-L. Dufresne,
%  A.~Hall, S.~Hallegatte, M.M. Holland, W.~Ingram, D.A. Randall, B.J. Soden,
%  G.~Tselioudis, and M.J. Webb.
%\newblock {How Well do we Understand and Evaluate Climate Change Feedback
%  Processes~?}
\REVIEW{Journal of Climate}{19}{2006}{3445}
%, 19:3445--3482, 2006.

\bibitem{Cess89}
\Name{Cess R. D. \etal}
% G.L. Potter, J.P. Blanchet, G.J. Boer, S.J. Ghan, J.T. Kiehl, T.H.
%  Le, X.Z. Li, X.Z. Liang, J.F.B. Mitchell, J.J. Morcrette, D.A. Randall, M.R.
%  Riches, E.~Roeckner, U.~Schlese, A.~Slingo, K.E. Taylor, W.M. Washington,
%  R.T. Wetherald, and I.~Yagai.
%\newblock Interpretation of cloud-climate feedback as produced by 14
%  atmospheric general circulation models.
\REVIEW{Science}{245}{1989}{513}
%--516, 1989.

\bibitem{Hansen84}
\Name{Hansen J. E. \etal}
%, D.~Lacis, D.~Rind, G.~Rusell, P.~Stone, I.~Fung, R.~Ruedy, and
%  J.~Lerner.
%\newblock \protect{Climate Sensitivity: Analysis of Feedback Mechanisms}.
\Editor{Hansen J. E. \and Takahashi T.}
\Book{Climate Processes and Climate Sensitivity}
\Vol{5}
\Page{130}
%--162. Geophysical monograph {\bf
%  29},
\Editor{AGU}
\Year{1984}

\bibitem{Schneider72}
\Name{Schneider E. K.}
%\newblock Cloudiness as a global climatic feedback mechanism - effects on
%  radiation balance and surface-temperature of variations in cloudiness.
\REVIEW{J. Atmos. Sci.}{29}{1972}{1413}
%, 29(8):1413--1422, 1972.

\bibitem{Bates07}
\Name{Bates J.R.}
%\newblock Some considerations of the concept of climate feedback.
\REVIEW{Quart. J. R. Met. Soc.}{133}{2007}{545}
%, 133(624):545--560, 2007.

\bibitem{HLG05}
\Name{Hallegatte S., Lahellec A., \and Grandpeix J.-Y.}
%\newblock \protect{An elicitation of the dynamic nature of water feedback in
%  Climate Change using a~1D model}.
\REVIEW{J. Atmos. Sci.}{63}{2006}{1878}
%, 63:1878--1894, 2006.

\bibitem{Lorenz63}
\Name{Lorenz E. N.}
%\Book{Deterministic nonperiodic flow}.
\REVIEW{J. Atmos. Sci.}{20}{1963}{130}
%, 20:130--141, 1963.

\bibitem{Luenberger}
\Name{Luenberger D. G.}
\Book{Introduction to Dynamic Systems}
\Publ{John Wiley \& Sons}
\Year{1979}

\bibitem{Cacuci81a}
\Name{Cacuci D. G.}
%\newblock Sensitivity theory for non-linear system: {I}. {N}onlinear
%  functionnal analysis approach.
\REVIEW{J. Math. Phys.}{22}{1981}{2794}
%, 22:2794--2802, 1981.

\bibitem{Cacuci84}
\Name{Cacuci D. G. \and Hall M. C. G.}
%\newblock Efficient estimation of feedback effects with application to climate
%  models.
\REVIEW{J. Atmos. Sci.}{41}{1984}{2063}
%, 41:2063--2068, 1984.

\end{thebibliography}

\end{document}
