Adaptive Detection of Instabilities: An Experimental Feasibility Study.
K. Krischer$^2$
G. Fl\"atgen$^2$
J.S. Anderson$^3$ and I.G. Kevrekidis$^{2,3}$
$^1$ Instituto Tecnol\'ogico de Celaya
Depto. de Ingenier\'{\i}a Qu\'{\i}mica
Celaya
Gto. 38010 M\'exico
$^2$ Fritz-Haber-Institut der Max-Planck-Gesellschaft
Faradayweg 4-6
14195 Berlin
Germany
$^3$ Department of Chemical Engineering
Princeton
NJ 08544 U.S.A.
We present an example of the practical implementation of a protocol for experimental bifurcation detection based on {on-line} identification and feedback control ideas. % The current experimental practice for the detection of bifurcations involves a cumbersome approach typically requiring long observation times in the vicinity of marginally stable solutions, as well as frequent re-settings of the experiment for the detection of turning-point or subcritical bifurcations. % The approach exemplified here addresses these issues drawing from ideas of numerical bifurcation detection. % The idea is to couple the experiment with an on-line computer-assisted identification/feedback protocol so that the closed-loop system will converge to the open-loop bifurcation points. % We demonstrate the applicability of this instability detection method by real-time, computer-assisted detection of period doubling bifurcations of an electronic circuit; the circuit implements an analog realization of the R\"ossler system. % The method succeeds in locating the bifurcation points even in the presence of modest experimental uncertainties, noise and limited resolution. % The results presented here include bifurcation detection experiments that rely on measurements of a single state variable and delay-based phase space reconstruction, as well as an example of tracing entire segments of a codimension-1 bifurcation boundary in two parameter space. \leftline{{PACS numbers: 05.45.,82.40}} \leftline{{Keywords:} Bifurcation Detection, Nonlinear Systems, Adaptive Control.}
adbifphysd8.tex

\documentstyle[12pt]{article} 
\include{./psfig}
\begin{document}

\title{Adaptive Detection of Instabilities: An Experimental
Feasibility Study.}
\author{R. Rico-Mart\'{\i}nez$^{1,2}$\thanks{
Corresponding author. E-mail: ramiro@arnold.princeton.edu, Telephone
+52-46117802, Fax +52-46117744},
K. Krischer$^2$, G. Fl\"atgen$^2$, \and J.S. Anderson$^3$ and
I.G. Kevrekidis$^{2,3}$\\
$^1$ Instituto Tecnol\'ogico de Celaya, Depto. de Ingenier\'{\i}a 
Qu\'{\i}mica,\\ Celaya, Gto. 38010 M\'exico \\
$^2$ Fritz-Haber-Institut der Max-Planck-Gesellschaft,\\
Faradayweg 4-6, 14195 Berlin, Germany\\
$^3$ Department of Chemical Engineering, Princeton University, 
\\Princeton, NJ 08544 U.S.A.\\}


\maketitle

\begin{abstract}

We present an example of the practical implementation of a
protocol for experimental
bifurcation detection based on {\it on-line}
identification and feedback control ideas.
%
The current experimental practice for the detection
of bifurcations involves a cumbersome
approach typically requiring long observation times in the
vicinity of marginally stable solutions, as well as
frequent re-settings of the experiment for the
detection of turning-point or subcritical bifurcations.
%
The approach exemplified here addresses these issues drawing from
ideas of numerical bifurcation detection.
%
The idea is to couple the experiment with an on-line computer-assisted
identification/feedback protocol so that the closed-loop
system will converge to the open-loop bifurcation points.
%
We demonstrate the applicability of
this instability detection method by real-time, computer-assisted
detection of period doubling bifurcations of an electronic circuit;
the circuit implements an analog realization of
the R\"ossler system.
%
The method succeeds in locating the bifurcation points even
in the presence of modest experimental uncertainties, noise and limited
resolution.
%
The results presented here include bifurcation detection
experiments that rely on measurements
of a single state variable and delay-based phase space reconstruction, as
well as an example of tracing entire segments of a codimension-1
bifurcation boundary in two parameter space.

\leftline{{\bf PACS numbers: 05.45.,82.40}}

\leftline{{\bf Keywords:} Bifurcation Detection, Nonlinear Systems, Adaptive
Control.}
\end{abstract}
\newpage


\section{{\bf Introduction}}
\thispagestyle{empty}
\noindent
Over the last 30 years increasing efforts have been directed towards
the understanding of experimentally observed complex dynamical behavior.
%
The development of Nonlinear Dynamics as a branch of applied 
mathematics provides us the set of tools 
necessary for understanding the response
of nonlinear dynamical systems 
as operating parameters are varied.
%
Characterization of the nonlinear
dynamical behavior of experimental systems, however, remains a difficult task.
%
This is mainly because real systems cannot be manipulated as 
easily or as cleanly as a model would on
the computer when its bifurcation behavior is traced.

Typically, an  experimental bifurcation diagram is
obtained via repeated application of a {\it set-and-observe}
procedure. The operating parameter of interest
is set and the system dynamics allowed to evolve asymptotically to stationary 
(e.g. steady or oscillatory) behavior.  
%
The operating parameter setting is then changed (e.g. incremented slightly)
and the experiment is repeated; if a qualitative difference is detected
between the asymptotic dynamics reached for two consecutive parameter
settings, a bifurcation is inferred to have occurred somewhere in-between.
%
A similar approach  can be applied to the {\it computational} detection
of bifurcations if one 
is only allowed to use forward-in-time integration of the model equations 
(i.e. set parameter values and initial conditions and run a simulation to
stationarity).
%
In both cases there are serious drawbacks.
%
For example, in the investigation of a period doubling bifurcation
of a periodic orbit, as the 
critical value of the search parameter is approached the dynamics will
slow down considerably, and the time required to determine whether ultimately 
single loop or double loop oscillations are
observed increases dramatically.
%
Furthermore, if one wishes to detect a subcritical bifurcation, or
a saddle-node bifurcation, once one steps over the critical parameter
value, the
system will move far away from the phase space region of interest.
%
The experiment will then need to be reset to reach the vicinity of the
stable solution located just before the critical parameter value
that was overstepped. 
%

In previous work \cite{ciss,prl} we proposed a protocol for the
experimental detection of bifurcations motivated by the procedure
of computational construction of bifurcation diagrams from mathematical
models.
%
Numerical bifurcation theory in a computational
environment allows the bifurcation location via a simultaneous
search for criticality in (phase)$\times$(parameter) space.
%
An augmented set of equations (e.g. steady state equations 
plus a criticality condition) is typically formed, and its fixed point 
(steady state values at criticality {\it and} critical
parameter value) is converged upon quadratically using a Newton-Raphson
contraction mapping for the augmented system of equations.
%
Obviously in this case the states do not evolve on the computer
based on the system
dynamics, but based on the computational contraction mapping dynamics.

The proposed protocol for the detection of bifurcations seeks an
``experimental-computational" compromise:
an experiment (i.e. states evolving based on system dynamics) is coupled
to a parameter evolution law which is based on criticality conditions and
feedback control. The augmented, closed loop system,
constructed in this
manner has a {\it stable} steady state at what used to be
the open loop system's {\it marginally stable} steady state.
%
The computational building blocks that
constitute the parameter evolution law are based on  the
{\it on-line} identification of
a ``local nonlinear" dynamic model, obtained via least-squares regression
\cite{PTVF},
used to estimate criticality on-line, and simple feedback control
algorithms that use the bifurcation parameter as the actuator in order to
drive the system towards the bifurcation point.
%
As it will be illustrated below,
this algorithm enhances our ability
to efficiently detect instabilities, trace segments of the
bifurcation diagram (``operating diagram"),
and, in general, characterize the system dynamics in multiparameter space.
%
The ``experimental" system used in our illustration is an analog realization
of the R\"ossler system. For this system, we converge experimentally
on the period doubling
bifurcation of simple periodic orbits in one and two-parameters. The
local adaptive model used to infer the presence of the bifurcation is a 2D nonlinear
map. Our illustrations are based on (a) full state vector measurement; and 
(b) reconstruction of the state-space based on time series of a single state 
measurement and the use of delay coordinates.

An alternative motivation for studying this type of problem is the
typical Ziegler-Nichols (ZN) tuning procedure for PID (proportional-
integral-derivative) controllers
\cite{ziegler},
which motivated the work of O'Neil, Lyberatos and Svoronos
\cite{Lyb} on adaptive determination of bifurcation points.
%
In the ``textbook" ZN procedure the ``experiment" is the 
closed loop system with a proportional controller; the gain of
this controller is the ``bifurcation parameter".
%
The closed-loop system loses stability when the proportional gain $K$
reaches the value $K_u$. 
%
This critical parameter value is precisely the ZN ``ultimate gain".
%
In practice, one obtains $K_u$ by slowly increasing the gain of 
the proportional controller
until one observes incipient oscillations.
%
Knowledge of the location of the bifurcation and the frequency of the oscillations
dictates then the ZN-tuning settings for a PID controller;
it is interesting that this practical controller tuning approach is based
on experimental bifurcation detection.
%
%the tuning of the controller by selecting
%the gain $K$ as a function of
%the distance from the nominal operation point to the ``bifurcation" point.

The paper is organized as follows: first, we describe
our system (the electronic R\"ossler circuit) and the bifurcation 
detection protocol before presenting
our experimental results. We then provide some concluding remarks. 

\begin{figure}[t]
\centerline{\psfig{file=./schemeanalog.eps,width=5.0in}}
\caption{Schematic of the experimental setup.}
\label{schem}
\end{figure}

\section{{\bf The R\"ossler Circuit}}

The R\"ossler model was proposed in the early 1970's \cite{rossler}
as one of ``the simplest" nonlinear vectorfields 
capable of generating chaotic behavior.
%
It consists of the following set of three ordinary differential equations:

\begin{eqnarray}
{{dx}\over{dt}} & = & -y -x \\
{{dy}\over{dt}} & = & x + ay \\
{{dz}\over{dt}} & = & b + z(x-c) 
\end{eqnarray}

This system has been extensively studied. For $b=2$ and $c=4$ and as $a$
is increased, it exhibits a steady state
undergoing a Hopf bifurcation followed by a cascade of period doubling
bifurcations of the periodic orbits culminating in an apparently chaotic 
attractor around $a=0.396$.

\begin{figure}[t]
\centerline{\psfig{file=./rosscomp.eps,width=5.0in}}
\caption{Comparison of attractors for the circuit (left) and the
numerically computed trajectory obtained integrating the set of ODEs
of the R\"ossler system. The comparison is made for a double-loop
trajectory observed at $a=0.35$, $b=2$ and $c=4$.}
\label{rosscomp}
\end{figure}

In the implementation of the protocol for the adaptive detection of
experimental bifurcations, we use as model system an analog realization
of the R\"ossler system. The R\"ossler equations
are implemented via a sequence of feedback loops of electronic signals.
The circuit was manufactured at the FHI labs\footnote{
A detailed diagram of the circuit is available in the web in
the following address
http://w3.rz-berlin.mpg.de/pc/spatdyn/spatdyn.cgi?Publications}
and includes an encasing
that allows both manual and digital manipulation of all three system parameters. 
%
In our experiments the parameter dynamics were computer-generated  in closed
loop, and were communicated to the device through its digital I/O channels. 
%
The circuit is coupled (Fig.~\ref{schem})
with a VMEbus System (32bit backplane system) consisting of a
VMEbus CPU-board Baja40 from Heurikon Division controlled
via a Motorola 68040 32-bit microprocessor chip running internally
with 66 MHz, externally with 33 MHz, 32 MByte RAM and
several I/O interfaces and a VMEbus IP carrier board with
four IndustryPack slots. One slot was used
with a TIP501 optically isolated 16 channel
16 bit ADC module from TEWS Datentechnik, input voltage +/- 10 volt,
14.5 usec data acquisition and conversion time,
accuracy +/- 4 LSB, linearity +/- 4 LSB, and the other slot with a
TIP550 optically isolated 4 channel 12 bit DAC from TEWS
Datentechnik, output voltage +/- 10 volt, 13 usec settling time,
accuracy +/- 1 LSB, linearity +/- 1 LSB. The array uses an
operating system VxWorks 5.3.1 (Tornado) from Wind River System.
The data acquisition
components allow us to specify the accuracy of the measurements from
11 to 16 bits.


The dynamic response of the motivating set of differential equations
is well reproduced. Figure~\ref{rosscomp} presents a comparison of a couple
of projections of the double-loop limit cycle predicted and observed at $a=0.35$,
$b=2$ and $c=4$. The deviation of the signal given by the circuit with respect
to the one obtained by integrating the differential equations is minimal,
both in the actual values of the variables as well as in the period of the
oscillations. Similar results are observed within the range of operation
of the circuit: $a \in [0,1]$, $b\in [0,10]$, $c\in [0,10]$ and values of
the state variables inside the range [-10,10].


\section{{\bf The Bifurcation Detection Protocol}}

The building blocks that constitute the bifurcation detection
protocol are schematically depicted in Figure \ref{schemecontrol}.
%
They involve four different modules: one that
identifies {\it on-line} an approximate model of the system, a second one
that uses the identified model to estimate the location of the
bifurcation point in phase$\times$parameter space; this second module
relies in the information provided by a test function (third
module) that contains the
conditions defining the type of bifurcation sought. 
%
The fourth and final module uses the estimated critical parameter value
(our estimate of the open-loop critical conditions) to devise a policy that
will bring and keep the system closer to open-loop criticality.
%
The repeated application of the algorithm, composed by these four
modules, drives the system to the true open-loop bifurcation point,
and holds it there.
%
In devising the algorithm, we assume
that at least one of the states of the experimental system can be
observed (measured), and that at least one operating parameter is available to
be manipulated (the bifurcation will be located with respect to this parameter).
%
Furthermore, we assume that the open-loop system is in the neighborhood of a
local bifurcation; thus its dynamics can be captured via a low-order
nonlinear model. 
%
Using such a low-order approximation to model the dynamics of the
system is motivated (and can be justified) by the existence 
of a normal form for the bifurcation.
%
The discrete-time model used in the system identification module of the
algorithm
((1) in Fig.~\ref{schemecontrol})
takes the form of a low order polynomial on the observed variables and the
parameter, motivated from Taylor series expansions in the neighborhood
of the bifurcation point, and the theory of Center Manifolds/Normal Forms.
%
The polynomial model is of the type:

\begin{eqnarray}
x_{i}(k+1) &\equiv& F_{i}[{{\bf x}}(k);p(k)] \nonumber \\
& = & a_{i} + \sum_{j=1}^n b_{i,j}x_j(k) +  c_{i}p(k) +
\sum_{j=1}^n \sum_{l=1}^n d_{i,j,l}x_j(k)x_l(k)\nonumber \\
&  & + \sum_{j=1}^n  e_{i,j}x_j(k)p(k)+
\sum_{j=1}^n \sum_{l=1}^n \sum_{m=1}^n f_{i,j,l,m}x_j(k)x_l(k)x_m(k)
\quad   (i = 1, \cdots,n)  \nonumber  
\end{eqnarray}

\begin{figure}[t]
\centerline{\psfig{file=./schemecontrol.eps,width=5.0in}}
\caption{Schematic diagram of the algorithm for the experimental
detection of bifurcations.}
\label{schemecontrol}
\end{figure}

\noindent
where ${\bf x}$ represents the $n$-dimensional
vector of state variables, $p$ is the parameter, $a,b,c,d,e$ and $f$
are the polynomial coefficients, and $k$ is the iteration counter.
%
Note that the polynomial is linear in the parameter and contains cubic terms
on the state variables. The normal forms of
the most common bifurcations are linear
in the parameter, and the normal form of the
Hopf bifurcation has cubic terms in the state variables \cite{Wig}.
%
The polynomial can be simpler for certain specific types of bifurcation;
the form presented here is appropriate for the detection of
local codimension-1 bifurcations.
%
For the experimental system, we rely in Poincar\'e sections 
and/or reconstructions based on delayed measurements to form the vector 
of states.
%
The identification module uses a least-squares algorithm \cite{PTVF} to fit the
coefficients of the polynomial from experimental observations;
other algorithms, such as the projection method, can of course also be used. 
%
The model is updated at every step as new observations become available.
%
In order to prevent ill-conditioning of the least-squares algorithm during
the experimental run, a well-known problem in the
{\it on-line} identification of experimental systems,
the parameter signal must be persistently excited
(here, by addition of random noise) \cite{astrom,G&S}.

The second module of our algorithm ((2) in Fig.~\ref{schemecontrol})
uses the polynomial model identified
by the first module to estimate the location of the open-loop bifurcation.
%
This is a computational step that
involves the simultaneous solution of the steady state equations
along with a criticality condition appropriate for the bifurcation sought
((3) in Fig.~\ref{schemecontrol}).
%
This is accomplished using an algebraic equation solver (usually a contraction
mapping like a Newton-Raphson). 
%
For example, for a period-doubling bifurcation, 
which we detected in the R\"ossler circuit, 
the bifurcation is located at the point where one Floquet multiplier
(an eigenvalue of the linearization around the fixed point of the
open-loop Poincar\'e map model) is equal to -1. 
%
Thus, a suitable test function, for
a two-variable discrete-time system, will be $det[{\bf{I}}+{\bf{J}}]= 0$,
where $\bf{I}$ is the identity matrix and
${\bf{J}}=[\nabla_x{\bf F}[{\bf x};p]]$ is the Jacobian of the
identified model at the steady state $({\bf x},p)$.

In order to find the location of the critical 
state ($\{{\bf \hat{x}}_{cr},\hat{p}_{cr}\}$),
one must solve for $\{\bf{x},p\}$ such that (period doubling case):

\begin{eqnarray}
x_i-F_{i}[{\bf{x}};p]& = & 0 \quad   (i = 1, \cdots,n) \nonumber \\
det[\bf{I}+\bf{J}] & = & 0  \nonumber  
\end{eqnarray}

The final module ((4) in Fig.~\ref{schemecontrol})
uses the identified model and the estimated location of
the bifurcation ($\{{\bf \hat{x}}_{cr},\hat{p}_{cr}\}$) to find a parameter
variation protocol that will drive the system to the bifurcation point.
This critical step can be accomplished using several alternative ``feedback
laws". Previously \cite{ciss,prl},
we had illustrated the use of one of the simplest:
a policy that will bring the system to criticality in the
minimal possible number of steps (dead-beat policy).
In the runs presented below, however, we rely
on ``softer" policies in order
to overcome problems associated with the local character of the model
identified, the limited resolution of measurements and
the need for continuous excitation of the control signal associated with
the ill-conditioning of the least squares identification procedure.
%
In addition, the control policy should be possible to be
evaluated and implemented very fast (almost instantaneously) in real time.
%
This is because we have ignored in our experiment the delays
associated with I/O tasks as well as those associated with
the computation of the control action - we have effectively
assumed that computation and implementation of control action
is instantaneous. 
%
This avoids a number of complications and possible instabilities 
having to do with physical delays in the control loop.
%
These are
important issues, and will be taken up in future work; they are,
however, more ``implementation difficulties", and do not directly
pertain to the mathematical underpinnings of our procedure.


With these restrictions in mind, we explored a couple of 
plausible two parameter variation
protocols: the first policy attempts to minimize the distance from the
critical point in the next iteration, penalizing large
parameter variations; 
the second protocol is based on decreasing the distance to the critical
state (in phase space) in the next iteration by a fixed fraction.
%
In the notation used above, the first protocol results from:

\begin{equation}
\min_{p(k)} \sum_{i=1}^n (x_i(k+1)-\hat{x}_{i,cr})^2 + W (p(k)-\hat{p}_{cr})^2
\end{equation}

Note that, since the identified model (used to predict
${\bf x}(k+1)$) is linear in the parameter $p$, the solution of
this minimization problem renders a linear equation for $p$.  
%
The parameter $W$ is a penalty weight, which we usually set in the range [0.1,10].
Since the model form used is only valid locally,
close to the open-loop bifurcation, its predictions away from the
bifurcation are not accurate. 
%
In the presence
of modeling errors, a dead-beat policy may drive
the system farther away from the bifurcation instead of bringing it
closer. 
%
The introduction of the weight $W$ gives us a mechanism to 
moderate excessive control action, and thus, to possibly 
overcome the effect of 
such modeling errors away from the bifurcation point.
%
A more appropriate implementation of this policy
would seek only to penalize actions that bring the system
outside the range of validity of the model. 
%
That is, the penalty should only be applied if the prescribed control action
is larger that the estimated range of validity of the identified model.
We are currently exploring this issue.

The second policy is motivated by the form of transients approaching a
steady-state near a Hopf bifurcation. In such a scenario the dynamical
system will slowly converge to a marginally stable steady state by
(possibly fast, but slowly decaying) spiraling towards it. 
%
Thus, this second policy attempts to control the
speed of the approach by requiring that the distance from the critical
state location in the next iteration becomes a fraction of the {\it current}
distance from  criticality. In mathematical terms we solve for $p(k)$ such
that:

\begin{equation}
\sum_{i=1}^n (x_i(k+1)-\hat{x}_{i,cr})^2 = C^2 \sum_{i=1}^n
(x_i(k)-\hat{x}_{i,cr})^2
\end{equation}


\noindent
%
Here $C$ is the desired fractional reduction, in the next iteration,
of the distance from estimated open-loop criticality. 
%
We typically set $C$ in the range [0.9,0.95].
Since the model to predict $x_i(k+1)$ is linear in $p$, the solutions
of this equation are the roots of a quadratic equation that can be solved
explicitly. The root selected is the one closer to the current value of
the parameter.

Both control policies must be applied repeatedly in order to drive the
closed loop 
system to the point that, for the open-loop system,
was a marginally stable solution. Although the number
of iterations
to reach the bifurcation may be large (e.g. more than 100
iterations), we have found that, even with these
simple policies, the closed-loop
system performs well, finding the bifurcation for every experiment
independently of the initial condition used.

By using the identification and control methods in the manner described, we 
circumvent the long settling periods associated with open-loop
near-critical dynamics. As previously illustrated in \cite{prl}
through computational modeling, this
methodology is also capable of driving the system back to steady 
state after (mildly) overshooting the critical value of the
bifurcation parameter even for 
``hard" bifurcations such as turning points or subcritical Hopf bifurcations.
%
In the following section, we 
illustrate our approach by locating
period-doubling bifurcations of 
periodic orbits for the R\"ossler circuit device.

\begin{figure}[t]
\centerline{\psfig{file=./fig1-16bitrun10.eps,width=5in}}
\caption{Experimental results for the location of the period doubling
bifurcation of the period-1 limit cycle in the R\"ossler circuit. After
a short transient (50 steps) the initial model is properly identified and the
circuit is driven to the bifurcation point. The parameters $b$ and
$c$ were set at 2.0 and 4.0 respectively.}
\label{fig1}
\end{figure}


\section{{\bf Experimental Runs.}}

A period-doubling cascade of periodic orbits (or fixed points of
maps or Poincar\'e maps) is perhaps one of the
best understood routes to deterministic chaos. 
%
The R\"ossler system exhibits this
type of transitions at $b=2$ and $c=4$ for increasing values of the
remaining parameter $a$. The first of such period doubling bifurcations
is observed for $a\approx 0.33485$. It is precisely this bifurcation that
we set to locate using the protocol described above.

Poincar\'e sections provide a natural framework for observing/reporting
period doubling of periodic orbits.
The construction of our approximate models is performed on
such a (transversally intersecting the flow) Poincar\'e section.
Initially, we measure all three state variables from
the circuit $(x,y,z)$ and define the Poincar\'e section as
the plane $x=0.2$, crossed towards increasing values of $x$. 
%
For this Poincar\'e section
the critical point is located at $y=-4.1409$ and $z=0.44601$.
%
In the proposed protocol, the model to be identified involves two state
variables (the coordinates $(y,z)$ on the Poincar\'e section) and one parameter
($a$). Thus, our model will try to predict where the trajectory
intersects the plane at the next iteration ($y_{k+1},z_{k+1}$)
based on the current
location ($y_k,z_k$) and the value of the parameter ($a_k=p_k$). For
the experiment, the intersection point on the surface defining the
Poincar\'e surface
was obtained via fast (between 200 to 1000 Hz) sampling of all variables and
linear interpolation between the sampled points once the crossing with
the surface was detected.

\begin{figure}[t]
\centerline{\psfig{file=./fig2-16bitrun14.eps,width=5in}}
\caption{Experimental run exhibiting the effect of the lower bound
for the control parameter. The parameter is initially ``trapped"
near the lower bound until a better model is identified thus restarting
the movement towards the bifurcation point. The
parameters $b$ and $c$ were set as indicated in Figure 4.}
\label{fig2}
\end{figure}

For the circuit, the period of the oscillations before the bifurcation
is approximately 0.62 s. In comparison with the R\"ossler ODE model, the
arbitrary units of time translate as seconds for the circuit device and the
arbitrary units of the state variables and
parameters as voltages measured in Volts.
%
Using a 16-bit data acquisition card, the measurement resolution for
the state variables is approximately 0.0006 volts and
0.00003 volts for the parameter $a$
which has a smaller range of variation (voltage). 
%
As we mentioned earlier,
in order to successfully identify a varying {\it on-line} model, one should
apply persistent excitation to the actuator (the parameter signal).
This is done in the
form of random noise with maximum amplitude of 0.003 volts. The excitation
signal prevents the system from reaching stationary states for which the
measurement will be effectively constant.
%
Such measurements will render the matrix used in the identification 
computations severely ill-conditioned, and the approximate model
will deteriorate. 
%
With persistent excitation signals the system
will be always ``on the move" allowing the collection of information-rich data
resulting in more accurate models. 
%
Since, however the final goal is to
drive the system to the open-loop bifurcation, very large excitation signals are 
undesirable because they work against the control actions implemented. 
%
The excitation level selected
represents a compromise between these competing effects. Depending on the
system under consideration, the magnitude of the excitation signal
should be adjusted.
%
In our experience with several different simulated systems, we have found that
the amplitude of the excitation signal must be at least one order of magnitude
larger than the resolution available in the measurement.

The use of this relatively large excitation signal combined with the soft
control policies described previously constrains the selection 
of the appropriate convergence criterion. 
%
The approach to the
critical state may be slow, possibly also ``erratic" due to the
effect of the excitation; on the other hand, the system may become trapped
and the limited
resolution may translate in biased models and incorrect criticality
predictions. 

\begin{figure}[t]
\centerline{\psfig{file=./fig4-16bitrun13.eps,width=5in}}
\caption{Experimental run exhibiting bursting.
The adaptive bifurcation driver settles for a while in a point near the
critical coordinates, without achieving convergence,
until near-singularity of the matrix used
for the identification produces a sudden burst that drives the
circuit far away. The
parameters $b$ and $c$ were set as indicated in Figure 4.}
\label{fig3}
\end{figure}

In order to prevent the latter occurrence, we opted for a
two-level convergence criterion: (1) we require that the norm of
the deviation from criticality in the extended
phase$\times$parameter space is sufficiently small (i.e. less than
0.008 volts) for a relatively small
number of consecutive iterations (typically 5); and (2) after complying
with the first criterion, we restart the system identification
procedure by gathering enough data
to construct a ``fresh" model, using the excitation signal
without control, and then restart the feedback control until the first criterion
is again satisfied. 
%
Overall convergence is declared only
if two consecutive restarts reach the same point within the tolerance
indicated above. The first criterion does not excessively burden the
controller performance, while the second seeks to prevent declaring
convergence at a possibly unconverged,
incorrect prediction of the location of the bifurcation.
%
The tolerance to declare convergence (0.008 volts in our case)
should clearly be larger than the set
maximum amplitude of the random noise (excitation
signal, 0.003 volts); a reasonable rule-of-thumb is to
set this tolerance half-an-order of magnitude above the maximum
amplitude of the random noise.
%
We note, however, that the achievable accuracy in locating the bifurcation
is limited, not because of the convergence criterion, but most
importantly because of the excitation signal, model sensitivity and,
ultimately, limited resolution in the measurements.

With the feedback protocols implemented, we observed, during our runs,
that large parameter variations (``control actions") 
arose frequently.  Since such
variations may drive the system outside the range physically accessible to
the circuit device (i.e. negative parameter values), 
or outside the range for
which sustained oscillations are observed (thus, preventing the
trajectories from intersecting the plane used in the definition of
the Poincar\'e section) we opted for establishing hard bounds
(``saturation") for the values that the parameter can take. 
%
If the control
action prescribed goes beyond these bounds, the parameter is
set to the value of the bound that has been surpassed (the actuator saturates).
For the run at $b=2$ and $c=4$, the minimum value allowed for $a$ was 0.12
and the maximum $0.5$. Below
$a\approx 0.125$ the only attractor of the system is a steady state,
while above $0.5$ the chaotic attractor observed in the range [0.396,0.5]
disappears via global bifurcations and the state variables of the
circuit tend to saturate. 

With the hardware available to us, the evaluation of the protocol,
including identification of an updated model, location of criticality
and evaluation of the control policy, involved between 30 and 60 milliseconds,
depending on the number of evaluations required for some of the subroutines
such as the algebraic solver to find criticality. The effect of such ``large"
delay (about one tenth of the period of oscillation in the worst case)
was often detrimental to the performance of the protocol, preventing the accurate
location of the bifurcation. In order to overcome this situation, we 
modified the evaluation order of the subroutines. That is, we first
calculate and implement a control policy using one-iteration-old models,
then we use the collected data to update the model and prediction of
criticality for the next iteration.
%
This strategy was successful, since the evaluation of the
control policies proposed require only the solution of linear or quadratic
equations and thus take little time (less than 3 ms), and the slow approach
to the critical state makes consecutive models very similar most of the
time (depending on the conditioning of the matrix involved in the
identification of the model).

Figure ~\ref{fig1} illustrates the performance of the protocol.
The parameter was initially set at $a=0.42$. At this parameter
setting, the open-loop circuit exhibits a chaotic attractor. 
%
The experiment is left to evolve in the vicinity
of this parameter value for about 30 iterations, while data are gathered
to construct an initial model. 
%
After approximately 70 iterations, the system
has been driven to the vicinity of the bifurcation. The last 50 iterations,
for which the response of the circuit is practically constant, represents
the period needed to declare convergence with one restart as described above.
For this run the critical state was predicted at $y=-4.1127$, $z=0.4537$
and $a=0.3389$.
%
This compares favorably with the location of the bifurcation
for the R\"ossler ODE system. 
%
As we indicated before, we anticipated that the
estimation of the bifurcation location will not be very accurate
because of the persistent excitation signal and soft control policies. 
%
A more representative result should be reported
after gathering statistics over a number of runs. On average, over 15 runs,
each run with different initial condition and initial parameter
value, the
bifurcation is predicted to be located at $a=0.3309$ with a standard deviation
of $0.0082$ volts, with the critical state located at $y=-4.0964$
and $z=0.4552$ and standard deviations of $0.0134$ and $0.0014$ volts
respectively. 
%
The deviation of this average is well within the range
allowed for declaring convergence in our runs and the standard deviation
commensurate with the excitation signal and convergence criteria.
We note that all our runs did converge to the critical state
albeit occasionally at a larger number of iterations than for Figure~\ref{fig1}.

\begin{figure}[t]
\centerline{\psfig{file=./fig6-11bitrun2.eps,width=5in}}
\caption{Experimental run with lower resolution. For a resolution based
on an 11 bit A/D converter, the bifurcation is also properly
located. The parameters $b$ and $c$ were set as indicated in Figure 4.}
\label{fig4}
\end{figure}

Figure~\ref{fig2} presents another run that illustrates the role of
hard, saturation bounds limiting the values that the parameter can take
during the search for criticality. 
%
After approximately 40 steps, an incorrect criticality prediction translates 
into a drastic control action.
Since the value prescribed for the parameter would be below the preset lower
bound, the parameter is set to the bound. The performance of the controller
deteriorates and the system is trapped near the lower bound for about 40
iterations before restarting its movement towards the critical
state.

Figure~\ref{fig3} shows an example of the phenomenon
commonly referred to as ``bursting".
The system starts its gradual approach toward criticality, but since
the coordinates of the points in the trajectory do not change much, the
matrix used for the identification becomes progressively more ill-conditioned,
eventually rendering a ``bad" model that, in turn, dictates a large parameter
excursion (``burst"). The trajectory to criticality exhibits, in this
example, a couple of instances of such bursting on its route to the
eventual correct detection of the critical fixed point.

For those runs initialized at parameter values lower than the bifurcation value,
as the one illustrated in
Figure~\ref{fig3}, the convergence to the critical point was
consistently declared for slightly lower values of the
parameter than the true critical parameter value. 
%
In contrast, when the initialization took place at conditions 
for which the open loop system exhibited apparently chaotic oscillations (i.e.
Figure~\ref{fig1}), the collection of information-rich data allowed
the accurate estimation of the critical fixed point location early on, 
resulting in faster convergence towards the bifurcation.

One potential limitation of the application of this strategy may be the
available measurement resolution. The examples presented above rely on a 
16 bit A/D card. However, as illustrated in Figure~\ref{fig4}, the
bifurcation was successfully located also for lower resolution, equivalent
to that provided by a standard 11 bit A/D converter (i.e. one order
of magnitude less in resolution for the parameter and state variables).

Although the examples discussed until this
point involve measurements of all the states of
the system, it is also possible to use only a single measurement time series.
%
Figure~\ref{fig5} illustrates the performance of the protocol for
the location of the same bifurcation using only measurements of the
variable $x$ and relying on delayed measurements to reconstruct the phase space.
The delay ($\tau$) for the reconstruction was set to 40 ms. In practice one
could seek to find an optimal value of the delay using, for example,
mutual information measures; this was not done here.
The Poincar\'e section was set as before using the plane $x(t)=0.2$.
The critical state is located for the R\"ossler system at
$x(t-\tau)=-1.1701$ and $x(t-2\tau)=-2.2020$. On average over 5 runs,
the bifurcation is predicted at $a=0.3235$ for
$x(t-\tau)=-1.1143$ and $x(t-2\tau)=-2.1320$, with standard deviations
of $0.0080$, $0.0042$ and $0.0010$ volts respectively. The trajectory to reach
the critical state, however, exhibits a larger number of
instances of bursting and saturation of the parameter and thus requires
a larger number of iterations as illustrated
in Figure~\ref{fig5}. This deterioration in performance of the protocol
can be attributed to the increased sensitivity of the identification
resulting from the use of a sub-optimal delay and the smaller effective
resolution of the sampled variable.


\begin{figure}[t]
\centerline{\psfig{file=./fig8b1meas0.25.eps,width=5in}}
\caption{Experimental run using  measurements of only one
variable and a delay-based reconstruction of
the phase-space as the basis to take the Poincar\'e section.
The delay was set to 40 ms, two delays were used and
the Poincar\'e section was set for the current value of the
measured variable (x) equal to 0.2.
The parameters $b$ and $c$ were set as indicated in Figure 4.}
\label{fig5}
\end{figure}

A final test of the capabilities of the proposed protocol for the
detection of bifurcations comes by attempting to draw, experimentally,
a two-parameter bifurcation diagram. 
%
That is, we attempt to trace automatically
a segment of our codimension-1 bifurcation boundary in two parameters. 
%
We vary the (second) parameter $b$ in the range [1,3]. 
%
Starting at $b=3.0$, we locate the value of $a$ for which the 
bifurcation is observed and then
decrease $b$ by a fixed amount (0.05 volts). 
%
Figure~\ref{fig6}
presents the comparison of these runs for our circuit device,
and the computed boundary for the R\"ossler ODE system. 
%
The results presented are the
average over five sweeps, three
calculated decreasing $b$ and two increasing it. The previously found
bifurcation point was used as the starting point for the new search
with a different value of $b$.
Figure~\ref{fig7} presents the comparison of the location
of the (open loop critical) fixed points for the Poincar\'e section.
%
The boundaries for the R\"ossler ODE system were
calculated using AUTO \cite{Doe}.
%
As can be seen in Figures~\ref{fig6} and ~\ref{fig7}
the segment is accurately located, as are
the critical states. For both plots, standard deviations have been
omitted; they are comparable to the previously indicated
values for the single parameter runs.

\begin{figure}[t]
\centerline{\psfig{file=./comptwoparbif.eps,width=5in,angle=270}}
\caption{Comparison of the experimentally
obtained two parameter boundary and the numerically calculated
from the R\"ossler System. The parameter $b$ was changed gradually
and the bifurcation point obtained from the previous value used as
starting point to locate the new bifurcation. The experimental
results are averaged over 5 runs, three runs using a variation
for $b$ decreasing and two for $b$ increasing.
The parameter $c$ was set to 4.0.}
\label{fig6}
\end{figure}

\section{{\bf Conclusions}}

We have presented a successful implementation of
an algorithm for efficiently driving an experimental system to 
(what at open loop would have been) criticality.
%
The control protocol was applied to an experiment in
the form of an electronic circuit device that
implemented the dynamics of the R\"ossler system.
%
The proposed algorithm is based on adaptive system identification, 
separation of time scales close to criticality, and on the
exploitation of the identified local model to search for the
criticality conditions; the estimation
of open-loop criticality is used to devise a control policy that drives the
closed-loop system towards the bifurcation.
%
The protocol seeks, in a sense, to mimic the procedure used in the numerical
construction of bifurcations from mathematical models:
the critical state of the open-loop system becomes a {\it stable} steady state of 
the {\it augmented} system in (phase) $\times$ (operating parameter) space.
%
We illustrated the capabilities of the algorithm in the location of
period doubling bifurcations of periodic orbits.
%
The bifurcation was successfully located even in the presence of limited
resolution.
%
We also illustrated the applicability of this approach when only
one state variable can be measured by reconstructing the phase space
using delayed measurements.
%
Our final illustration involved the tracing of a segment of the
two parameter bifurcation boundary.

Our selection of an electronic circuit device with R\"ossler dynamics,
to demonstrate the experimental feasibility of the protocol for
the automatic detection of bifurcations, was intended to provide
a testing ground for assessing the performance of the protocol.
%
The circuit device provides us with a highly reproducible and easily 
manipulated experiment with a well charted dynamical behavior.
%
We chose the period doubling bifurcations of the oscillations
of the circuit as our representative codimension-1 bifurcation.
%
The low-dimensionality of the model used crucially depends on the 
time-scale gap between the long time scale associated with
the critical eigenvalue(s) and the time scale of the ``next least
stable" eigenvalue. 
%
It is important that all operations performed in our system 
(I/O as well as software computations) are fast compared to 
this ``next least slow" time scale.
%
Since the return time for our Poincar\'e map was roughly 600 ms,
our relatively short (approximately 40 ms) 
action delay could be considered ``approximately instantaneous".
%
An alternative that would speed up the protocol evaluation is the use of
recursive least squares techniques \cite{astrom,G&S}
that may reduce significantly
the time involved in the identification of a local model; however, this
alternative resulted in more sensitive and prone-to-burst trajectories,
even though we attempted several different strategies in the assignment of
forgetting factors for the identification.
%
Other alternatives, such as the use of projection algorithms \cite{G&S},
are also feasible, but were not investigated.

Perhaps the most limiting aspect for the widespread application of the
proposed protocol involves the high sensitivity of the strategy
as a result of the combination of the control policies, excitation
signal and on-line identification of a model with only local validity.
%
Such sensitivity translates in erratic and slow motion towards the
bifurcation, frequent bursts, trapping of trajectories near hard bounds
and ill-conditioning of the identification
process.

\begin{figure}[t]
\centerline{\psfig{file=./comptwoparloc.eps,width=5in,angle=270}}
\caption{Comparison of the experimentally
converged location of the marginally stable fixed points
and the numerically calculated ones 
from the R\"ossler system. The parameter $b$ was changed as indicated
in the previous figure and $c$ was set to 4.0.}
\label{fig7}
\end{figure}


%
We are currently working on posing the entire problem (system, identification
and control) as an augmented dynamical system, for which the existence 
of an overall converged state can be demonstrated, and its stability quantified.
%
Such an approach could also allow us to explore issues related to the
performance in the presence of multiple nearby bifurcations of the same type.

Extensions of the algorithm are straightforward: higher codimension
bifurcations can be detected by solving the corresponding (more than one)
criticality conditions for the (open-loop) marginally stable steady state and the
(more than one) critical parameters.
%
In principle one might even be able to accommodate unstable, nonbifurcating
parts of the spectrum on the right-hand plane.
%
Several important elements can be improved upon:
(a) Formulation of the criticality conditions in a less
naive and computationally more efficient manner, along lines of modern
numerical bifurcation theory (e.g. \cite{Doe}
and references therein); (b) Usage of more sophisticated
(e.g. receding horizon) feedback control laws \cite{rawlings};
and (c) Improvement of the type of local model used.
%
While in this problem we used polynomial models based on Taylor-series
expansions and normal forms, it is obvious that other
types of nonlinear dynamic models (e.g. models based on artificial neural
networks) can also be used in this context; we, among many others, have 
demonstrated that such models are capable of  accurately predicting
the bifurcation structure of a nonlinear system when trained with time series
data \cite{greece,bulsari}.

While we are continuing theoretical and computational work along the directions 
described above, we are also implementing this approach in the
study and characterization of the bifurcations of electrochemical
systems \cite{Wol}.
%

It may be interesting to note that this approach can also be coupled with
large scale scientific
simulators, for which the relevant long term dynamics are low-dimensional,
in order to chart their bifurcation scenario \cite{copeland} (as
opposed to writing a bifurcation algorithm for them from scratch).
%
The strategy of seeking and stabilizing a fixed point with one additional
condition (here, open-loop criticality), falls in the same category as
extremum seeking controllers \cite{krstic}; controller developments along
these lines would find direct application in our protocols, and vice versa.
%
We believe that an approach like the one described here can eventually become a
standard real-time tool which will assist experimentalists in
efficiently detecting bifurcations and
tracing operating diagrams; such a tool could drastically
enhance the experimental characterization of nonlinear dynamics
in multiparameter space.
%

\section{{\bf Acknowledgements}}

We are grateful to Professor Ertl for his interest in 
and generous support of this work.
We thank  Heinz Junkes and Georg
Heyne from the FHI technical staff for technical support.
We gratefully acknowledge the support of CONACYT and COSNET (RRM),
NATO (IGK, KK), and AFOSR (IGK) for partial funding of this research. 
IGK and RRM gratefully acknowledge
the support from the Alexander von Humboldt Foundation.

\vspace{.1in}
\begin{thebibliography}{99}
%
\bibitem{ciss}
J. S. Anderson, I.G. Kevrekidis, R. Rico-Mart\'{\i}nez and
K. Krischer, Use of adaptation and feedback in the experimental detection
of bifurcations, {\em Proc. 30th Annual Conf. on Information Sciences and
Systems}, {\bf 1}:155--160, 1996.
%
\bibitem{prl}
J. S. Anderson, S. Y. Shvartsman, G. Fl\"atgen, I. G. Kevrekidis,
R. Rico-Mart\'{\i}nez and K. Krischer, An adaptive method for the
experimental detection of instabilities, {\it Physical Review Letters},
{\bf 82}:532-535, 1999.
%
\bibitem{astrom}
K.J. Astr\"om and B. Wittenmark, {\em Adaptive Control}, Addison-Wesley,
New York, NY, 1989.
%
\bibitem{copeland}
G. S. Copeland, I.G. Kevrekidis and R. Rico-Mart\'{\i}nez,
Adaptive detection of instabilities and nonlinear analysis of a
reduced-order model for flutter and rotating stall in turbomachinery,
{\it Proceedings of the IEEE International Conference of Control
Applications} (CCA, Hawaii), Paper TuP3-1, 1146-1150, 1999.
%
\bibitem{Doe}
E. Doedel, H.B. Keller, and J.P. Kern\'evez,
Numerical analysis and control of bifurcation problems:
(I) Bifurcation in finite dimensions, {\it Int. J. of Bifurcation
and Chaos}, {\bf 1}:493--520, 1991.
%

\bibitem{rawlings}
J. W. Eaton and J. B. Rawlings, Feedback control of nonlinear
processes using on-line optimization techniques, {\it Computers chem. Engng},
{\bf 14}:469-479, 1990.
%
\bibitem{G&S}
G.C. Goodwin and K.S. Sin, {\em Adaptive Filtering, Prediction,
and Control}, Prentice-Hall, Englewood Cliffs, NJ, 1984.
%
\bibitem{Wol}
K. Krischer, 
Principles of temporal and spatial pattern formation in electrochemical 
systems, {\it Mod. Aspects Electrochem.},
{\bf 32}:1--142, 1999.
%
\bibitem{Lyb}
D.G. O'Neil, G. Lyberatos, and S.A. Svoronos,
Adaptive Determination of Bifurcation Points, Presentation at
{\it AIChE Annual Meeting}, Paper 172j, Washington, D.C., 1988.
%
\bibitem{PTVF}
W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery,
{\em Numerical Recipes in C}, 2nd ed., Cambridge University Press,
New York, NY, 1992.
%
\bibitem{greece}
R. Rico-Mart\'{\i}nez, J.S. Anderson and I.G. Kevrekidis,
Continuous-time nonlinear signal processing: A neural network
based approach for gray box identification,
{\it Proc. 1994 IEEE Workshop on Neural Networks for Signal Processing},
596--605, 1994.
%
\bibitem{bulsari}
R. Rico-Mart\'{\i}nez, I.G. Kevrekidis and K. Krischer,
Nonlinear system identification using neural networks: dynamics and
instabilities, Chapter 16 in the book {\em Neural Networks for Chemical
Engineers}, A.B. Bulsari (ed.), Elsevier Science B.V., 1995.
%
\bibitem{rossler}
O. E. R\"ossler, 
An equation for continuous chaos, {\it Phys. Lett.}, {\bf57A}:397--398, 1976.
%
\bibitem{krstic}
H-H Wang, S. Yeung and M. Krstic,
Experimental application of extremum seeking on an axial flow compressor,
{\it IEEE Trans. on Control System Technology},
{\bf 8}:300-309, 2000.
%
\bibitem{Wig}
S. Wiggins, {\em Introduction to Applied Nonlinear Dynamical 
Systems and Chaos}, Springer-Verlag, New York, NY, 1990.
%
\bibitem{ziegler}
J. G. Ziegler and N. B. Nichols, Optimum Settings for Automatic Controllers.
{\it Trans. ASME}, {\bf 64}:759, 1942.
%

\end{thebibliography}



\end{document}

psfig.tex
% Psfig/TeX Release 1.8
% dvips version
%
% All psfig/tex software, documentation, and related files
% in this distribution of psfig/tex are 
% Copyright 1987, 1988, 1991 Trevor J. Darrell
%
% Permission is granted for use and non-profit distribution of psfig/tex 
% providing that this notice is clearly maintained. The right to
% distribute any portion of psfig/tex for profit or as part of any commercial
% product is specifically reserved for the author(s) of that portion.
%
% *** Feel free to make local modifications of psfig as you wish,
% *** but DO NOT post any changed or modified versions of ``psfig''
% *** directly to the net. Send them to me and I'll try to incorporate
% *** them into future versions. If you want to take the psfig code 
% *** and make a new program (subject to the copyright above), distribute it, 
% *** (and maintain it) that's fine, just don't call it psfig.
%
% Bugs and improvements to trevor@media.mit.edu.
%
% Thanks to Greg Hager (GDH) and Ned Batchelder for their contributions
% to the original version of this project.
%
% Modified by J. Daniel Smith on 9 October 1990 to accept the
% %%BoundingBox: comment with or without a space after the colon.  Stole
% file reading code from Tom Rokicki's EPSF.TEX file (see below).
%
% More modifications by J. Daniel Smith on 29 March 1991 to allow the
% the included PostScript figure to be rotated.  The amount of
% rotation is specified by the "angle=" parameter of the \psfig command.
%
% Modified by Robert Russell on June 25, 1991 to allow users to specify
% .ps filenames which don't yet exist, provided they explicitly provide
% boundingbox information via the \psfig command. Note: This will only work
% if the "file=" parameter follows all four "bb???=" parameters in the
% command. This is due to the order in which psfig interprets these params.
%
%  3 Jul 1991	JDS	check if file already read in once
%  4 Sep 1991	JDS	fixed incorrect computation of rotated
%			bounding box
% 25 Sep 1991	GVR	expanded synopsis of \psfig
% 14 Oct 1991	JDS	\fbox code from LaTeX so \psdraft works with TeX
%			changed \typeout to \ps@typeout
% 17 Oct 1991	JDS	added \psscalefirst and \psrotatefirst
%

% From: gvr@cs.brown.edu (George V. Reilly)
%
% \psdraft	draws an outline box, but doesn't include the figure
%		in the DVI file.  Useful for previewing.
%
% \psfull	includes the figure in the DVI file (default).
%
% \psscalefirst width= or height= specifies the size of the figure
% 		before rotation.
% \psrotatefirst (default) width= or height= specifies the size of the
% 		 figure after rotation.  Asymetric figures will
% 		 appear to shrink.
%
% \psfigurepath#1	sets the path to search for the figure
%
% \psfig
% usage: \psfig{file=, figure=, height=, width=,
%			bbllx=, bblly=, bburx=, bbury=,
%			rheight=, rwidth=, clip=, angle=, silent=}
%
%	"file" is the filename.  If no path name is specified and the
%		file is not found in the current directory,
%		it will be looked for in directory \psfigurepath.
%	"figure" is a synonym for "file".
%	By default, the width and height of the figure are taken from
%		the BoundingBox of the figure.
%	If "width" is specified, the figure is scaled so that it has
%		the specified width.  Its height changes proportionately.
%	If "height" is specified, the figure is scaled so that it has
%		the specified height.  Its width changes proportionately.
%	If both "width" and "height" are specified, the figure is scaled
%		anamorphically.
%	"bbllx", "bblly", "bburx", and "bbury" control the PostScript
%		BoundingBox.  If these four values are specified
%               *before* the "file" option, the PSFIG will not try to
%               open the PostScript file.
%	"rheight" and "rwidth" are the reserved height and width
%		of the figure, i.e., how big TeX actually thinks
%		the figure is.  They default to "width" and "height".
%	The "clip" option ensures that no portion of the figure will
%		appear outside its BoundingBox.  "clip=" is a switch and
%		takes no value, but the `=' must be present.
%	The "angle" option specifies the angle of rotation (degrees, ccw).
%	The "silent" option makes \psfig work silently.
%

% check to see if macros already loaded in (maybe some other file says
% "\input psfig") ...
\ifx\undefined\psfig\else\endinput\fi

%
% from a suggestion by eijkhout@csrd.uiuc.edu to allow
% loading as a style file:
\edef\psfigRestoreAt{\catcode`@=\number\catcode`@\relax}
\catcode`\@=11\relax
\newwrite\@unused
\def\ps@typeout#1{{\let\protect\string\immediate\write\@unused{#1}}}
\ps@typeout{psfig/tex 1.8}

%% Here's how you define your figure path.  Should be set up with null
%% default and a user useable definition.

\def\figurepath{./}
\def\psfigurepath#1{\edef\figurepath{#1}}

%
% @psdo control structure -- similar to Latex @for.
% I redefined these with different names so that psfig can
% be used with TeX as well as LaTeX, and so that it will not 
% be vunerable to future changes in LaTeX's internal
% control structure,
%
\def\@nnil{\@nil}
\def\@empty{}
\def\@psdonoop#1\@@#2#3{}
\def\@psdo#1:=#2\do#3{\edef\@psdotmp{#2}\ifx\@psdotmp\@empty \else
    \expandafter\@psdoloop#2,\@nil,\@nil\@@#1{#3}\fi}
\def\@psdoloop#1,#2,#3\@@#4#5{\def#4{#1}\ifx #4\@nnil \else
       #5\def#4{#2}\ifx #4\@nnil \else#5\@ipsdoloop #3\@@#4{#5}\fi\fi}
\def\@ipsdoloop#1,#2\@@#3#4{\def#3{#1}\ifx #3\@nnil 
       \let\@nextwhile=\@psdonoop \else
      #4\relax\let\@nextwhile=\@ipsdoloop\fi\@nextwhile#2\@@#3{#4}}
\def\@tpsdo#1:=#2\do#3{\xdef\@psdotmp{#2}\ifx\@psdotmp\@empty \else
    \@tpsdoloop#2\@nil\@nil\@@#1{#3}\fi}
\def\@tpsdoloop#1#2\@@#3#4{\def#3{#1}\ifx #3\@nnil 
       \let\@nextwhile=\@psdonoop \else
      #4\relax\let\@nextwhile=\@tpsdoloop\fi\@nextwhile#2\@@#3{#4}}
% 
% \fbox is defined in latex.tex; so if \fbox is undefined, assume that
% we are not in LaTeX.
% Perhaps this could be done better???
\ifx\undefined\fbox
% \fbox code from modified slightly from LaTeX
\newdimen\fboxrule
\newdimen\fboxsep
\newdimen\ps@tempdima
\newbox\ps@tempboxa
\fboxsep = 3pt
\fboxrule = .4pt
\long\def\fbox#1{\leavevmode\setbox\ps@tempboxa\hbox{#1}\ps@tempdima\fboxrule
    \advance\ps@tempdima \fboxsep \advance\ps@tempdima \dp\ps@tempboxa
   \hbox{\lower \ps@tempdima\hbox
  {\vbox{\hrule height \fboxrule
          \hbox{\vrule width \fboxrule \hskip\fboxsep
          \vbox{\vskip\fboxsep \box\ps@tempboxa\vskip\fboxsep}\hskip 
                 \fboxsep\vrule width \fboxrule}
                 \hrule height \fboxrule}}}}
\fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% file reading stuff from epsf.tex
%   EPSF.TEX macro file:
%   Written by Tomas Rokicki of Radical Eye Software, 29 Mar 1989.
%   Revised by Don Knuth, 3 Jan 1990.
%   Revised by Tomas Rokicki to accept bounding boxes with no
%      space after the colon, 18 Jul 1990.
%   Portions modified/removed for use in PSFIG package by
%      J. Daniel Smith, 9 October 1990.
%
\newread\ps@stream
\newif\ifnot@eof       % continue looking for the bounding box?
\newif\if@noisy        % report what you're making?
\newif\if@atend        % %%BoundingBox: has (at end) specification
\newif\if@psfile       % does this look like a PostScript file?
%
% PostScript files should start with `%!'
%
{\catcode`\%=12\global\gdef\epsf@start{%!}}
\def\epsf@PS{PS}
%
\def\epsf@getbb#1{%
%
%   The first thing we need to do is to open the
%   PostScript file, if possible.
%
\openin\ps@stream=#1
\ifeof\ps@stream\ps@typeout{Error, File #1 not found}\else
%
%   Okay, we got it. Now we'll scan lines until we find one that doesn't
%   start with %. We're looking for the bounding box comment.
%
   {\not@eoftrue \chardef\other=12
    \def\do##1{\catcode`##1=\other}\dospecials \catcode`\ =10
    \loop
       \if@psfile
	  \read\ps@stream to \epsf@fileline
       \else{
	  \obeyspaces
          \read\ps@stream to \epsf@tmp\global\let\epsf@fileline\epsf@tmp}
       \fi
       \ifeof\ps@stream\not@eoffalse\else
%
%   Check the first line for `%!'.  Issue a warning message if its not
%   there, since the file might not be a PostScript file.
%
       \if@psfile\else
       \expandafter\epsf@test\epsf@fileline:. \\%
       \fi
%
%   We check to see if the first character is a % sign;
%   if so, we look further and stop only if the line begins with
%   `%%BoundingBox:' and the `(atend)' specification was not found.
%   That is, the only way to stop is when the end of file is reached,
%   or a `%%BoundingBox: llx lly urx ury' line is found.
%
          \expandafter\epsf@aux\epsf@fileline:. \\%
       \fi
   \ifnot@eof\repeat
   }\closein\ps@stream\fi}%
%
% This tests if the file we are reading looks like a PostScript file.
%
\long\def\epsf@test#1#2#3:#4\\{\def\epsf@testit{#1#2}
			\ifx\epsf@testit\epsf@start\else
\ps@typeout{Warning! File does not start with `\epsf@start'.  It may not be a PostScript file.}
			\fi
			\@psfiletrue} % don't test after 1st line
%
%   We still need to define the tricky \epsf@aux macro. This requires
%   a couple of magic constants for comparison purposes.
%
{\catcode`\%=12\global\let\epsf@percent=%\global\def\epsf@bblit{%BoundingBox}}
%
%
%   So we're ready to check for `%BoundingBox:' and to grab the
%   values if they are found.  We continue searching if `(at end)'
%   was found after the `%BoundingBox:'.
%
\long\def\epsf@aux#1#2:#3\\{\ifx#1\epsf@percent
   \def\epsf@testit{#2}\ifx\epsf@testit\epsf@bblit
	\@atendfalse
        \epsf@atend #3 . \\%
	\if@atend	
	   \if@verbose{
		\ps@typeout{psfig: found `(atend)'; continuing search}
	   }\fi
        \else
        \epsf@grab #3 . . . \\%
        \not@eoffalse
        \global\no@bbfalse
        \fi
   \fi\fi}%
%
%   Here we grab the values and stuff them in the appropriate definitions.
%
\def\epsf@grab #1 #2 #3 #4 #5\\{%
   \global\def\epsf@llx{#1}\ifx\epsf@llx\empty
      \epsf@grab #2 #3 #4 #5 .\\\else
   \global\def\epsf@lly{#2}%
   \global\def\epsf@urx{#3}\global\def\epsf@ury{#4}\fi}%
%
% Determine if the stuff following the %%BoundingBox is `(atend)'
% J. Daniel Smith.  Copied from \epsf@grab above.
%
\def\epsf@atendlit{(atend)} 
\def\epsf@atend #1 #2 #3\\{%
   \def\epsf@tmp{#1}\ifx\epsf@tmp\empty
      \epsf@atend #2 #3 .\\\else
   \ifx\epsf@tmp\epsf@atendlit\@atendtrue\fi\fi}


% End of file reading stuff from epsf.tex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% trigonometry stuff from "trig.tex"
\chardef\letter = 11
\chardef\other = 12

\newif \ifdebug %%% turn me on to see TeX hard at work ...
\newif\ifc@mpute %%% don't need to compute some values
\c@mputetrue % but assume that we do

\let\then = \relax
\def\r@dian{pt }
\let\r@dians = \r@dian
\let\dimensionless@nit = \r@dian
\let\dimensionless@nits = \dimensionless@nit
\def\internal@nit{sp }
\let\internal@nits = \internal@nit
\newif\ifstillc@nverging
\def \Mess@ge #1{\ifdebug \then \message {#1} \fi}

{ %%% Things that need abnormal catcodes %%%
	\catcode `\@ = \letter
	\gdef \nodimen {\expandafter \n@dimen \the \dimen}
	\gdef \term #1 #2 #3%
	       {\edef \t@ {\the #1}%%% freeze parameter 1 (count, by value)
		\edef \t@@ {\expandafter \n@dimen \the #2\r@dian}%
				   %%% freeze parameter 2 (dimen, by value)
		\t@rm {\t@} {\t@@} {#3}%
	       }
	\gdef \t@rm #1 #2 #3%
	       {{%
		\count 0 = 0
		\dimen 0 = 1 \dimensionless@nit
		\dimen 2 = #2\relax
		\Mess@ge {Calculating term #1 of \nodimen 2}%
		\loop
		\ifnum	\count 0 < #1
		\then	\advance \count 0 by 1
			\Mess@ge {Iteration \the \count 0 \space}%
			\Multiply \dimen 0 by {\dimen 2}%
			\Mess@ge {After multiplication, term = \nodimen 0}%
			\Divide \dimen 0 by {\count 0}%
			\Mess@ge {After division, term = \nodimen 0}%
		\repeat
		\Mess@ge {Final value for term #1 of 
				\nodimen 2 \space is \nodimen 0}%
		\xdef \Term {#3 = \nodimen 0 \r@dians}%
		\aftergroup \Term
	       }}
	\catcode `\p = \other
	\catcode `\t = \other
	\gdef \n@dimen #1pt{#1} %%% throw away the ``pt''
}

\def \Divide #1by #2{\divide #1 by #2} %%% just a synonym

\def \Multiply #1by #2%%% allows division of a dimen by a dimen
       {{%%% should really freeze parameter 2 (dimen, passed by value)
	\count 0 = #1\relax
	\count 2 = #2\relax
	\count 4 = 65536
	\Mess@ge {Before scaling, count 0 = \the \count 0 \space and
			count 2 = \the \count 2}%
	\ifnum	\count 0 > 32767 %%% do our best to avoid overflow
	\then	\divide \count 0 by 4
		\divide \count 4 by 4
	\else	\ifnum	\count 0 < -32767
		\then	\divide \count 0 by 4
			\divide \count 4 by 4
		\else
		\fi
	\fi
	\ifnum	\count 2 > 32767 %%% while retaining reasonable accuracy
	\then	\divide \count 2 by 4
		\divide \count 4 by 4
	\else	\ifnum	\count 2 < -32767
		\then	\divide \count 2 by 4
			\divide \count 4 by 4
		\else
		\fi
	\fi
	\multiply \count 0 by \count 2
	\divide \count 0 by \count 4
	\xdef \product {#1 = \the \count 0 \internal@nits}%
	\aftergroup \product
       }}

\def\r@duce{\ifdim\dimen0 > 90\r@dian \then   % sin(x+90) = sin(180-x)
		\multiply\dimen0 by -1
		\advance\dimen0 by 180\r@dian
		\r@duce
	    \else \ifdim\dimen0 < -90\r@dian \then  % sin(-x) = sin(360+x)
		\advance\dimen0 by 360\r@dian
		\r@duce
		\fi
	    \fi}

\def\Sine#1%
       {{%
	\dimen 0 = #1 \r@dian
	\r@duce
	\ifdim\dimen0 = -90\r@dian \then
	   \dimen4 = -1\r@dian
	   \c@mputefalse
	\fi
	\ifdim\dimen0 = 90\r@dian \then
	   \dimen4 = 1\r@dian
	   \c@mputefalse
	\fi
	\ifdim\dimen0 = 0\r@dian \then
	   \dimen4 = 0\r@dian
	   \c@mputefalse
	\fi
%
	\ifc@mpute \then
        	% convert degrees to radians
		\divide\dimen0 by 180
		\dimen0=3.141592654\dimen0
%
		\dimen 2 = 3.1415926535897963\r@dian %%% a well-known constant
		\divide\dimen 2 by 2 %%% we only deal with -pi/2 : pi/2
		\Mess@ge {Sin: calculating Sin of \nodimen 0}%
		\count 0 = 1 %%% see power-series expansion for sine
		\dimen 2 = 1 \r@dian %%% ditto
		\dimen 4 = 0 \r@dian %%% ditto
		\loop
			\ifnum	\dimen 2 = 0 %%% then we've done
			\then	\stillc@nvergingfalse 
			\else	\stillc@nvergingtrue
			\fi
			\ifstillc@nverging %%% then calculate next term
			\then	\term {\count 0} {\dimen 0} {\dimen 2}%
				\advance \count 0 by 2
				\count 2 = \count 0
				\divide \count 2 by 2
				\ifodd	\count 2 %%% signs alternate
				\then	\advance \dimen 4 by \dimen 2
				\else	\advance \dimen 4 by -\dimen 2
				\fi
		\repeat
	\fi		
			\xdef \sine {\nodimen 4}%
       }}

% Now the Cosine can be calculated easily by calling \Sine
\def\Cosine#1{\ifx\sine\UnDefined\edef\Savesine{\relax}\else
		             \edef\Savesine{\sine}\fi
	{\dimen0=#1\r@dian\advance\dimen0 by 90\r@dian
	 \Sine{\nodimen 0}
	 \xdef\cosine{\sine}
	 \xdef\sine{\Savesine}}}	      
% end of trig stuff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\def\psdraft{
	\def\@psdraft{0}
	%\ps@typeout{draft level now is \@psdraft \space . }
}
\def\psfull{
	\def\@psdraft{100}
	%\ps@typeout{draft level now is \@psdraft \space . }
}

\psfull

\newif\if@scalefirst
\def\psscalefirst{\@scalefirsttrue}
\def\psrotatefirst{\@scalefirstfalse}
\psrotatefirst

\newif\if@draftbox
\def\psnodraftbox{
	\@draftboxfalse
}
\def\psdraftbox{
	\@draftboxtrue
}
\@draftboxtrue

\newif\if@prologfile
\newif\if@postlogfile
\def\pssilent{
	\@noisyfalse
}
\def\psnoisy{
	\@noisytrue
}
\psnoisy
%%% These are for the option list.
%%% A specification of the form a = b maps to calling \@p@@sa{b}
\newif\if@bbllx
\newif\if@bblly
\newif\if@bburx
\newif\if@bbury
\newif\if@height
\newif\if@width
\newif\if@rheight
\newif\if@rwidth
\newif\if@angle
\newif\if@clip
\newif\if@verbose
\def\@p@@sclip#1{\@cliptrue}


\newif\if@decmpr

%%% GDH 7/26/87 -- changed so that it first looks in the local directory,
%%% then in a specified global directory for the ps file.
%%% RPR 6/25/91 -- changed so that it defaults to user-supplied name if
%%% boundingbox info is specified, assuming graphic will be created by
%%% print time.
%%% TJD 10/19/91 -- added bbfile vs. file distinction, and @decmpr flag

\def\@p@@sfigure#1{\def\@p@sfile{null}\def\@p@sbbfile{null}
	        \openin1=#1.bb
		\ifeof1\closein1
	        	\openin1=\figurepath#1.bb
			\ifeof1\closein1
			        \openin1=#1
				\ifeof1\closein1%
				       \openin1=\figurepath#1
					\ifeof1
					   \typeout{Error, File #1 not found}
						\if@bbllx\if@bblly
				   		\if@bburx\if@bbury
			      				\def\@p@sfile{#1}%
			      				\def\@p@sbbfile{#1}%
							\@decmprfalse
				  	   	\fi\fi\fi\fi
					\else\closein1
				    		\def\@p@sfile{\figurepath#1}%
				    		\def\@p@sbbfile{\figurepath#1}%
						\@decmprfalse
	                       		\fi%
			 	\else\closein1%
					\def\@p@sfile{#1}
					\def\@p@sbbfile{#1}
					\@decmprfalse
			 	\fi
			\else
				\def\@p@sfile{\figurepath#1}
				\def\@p@sbbfile{\figurepath#1.bb}
				\@decmprtrue
			\fi
		\else
			\def\@p@sfile{#1}
			\def\@p@sbbfile{#1.bb}
			\@decmprtrue
		\fi}

\def\@p@@sfile#1{\@p@@sfigure{#1}}

\def\@p@@sbbllx#1{
		%\ps@typeout{bbllx is #1}
		\@bbllxtrue
		\dimen100=#1
		\edef\@p@sbbllx{\number\dimen100}
}
\def\@p@@sbblly#1{
		%\ps@typeout{bblly is #1}
		\@bbllytrue
		\dimen100=#1
		\edef\@p@sbblly{\number\dimen100}
}
\def\@p@@sbburx#1{
		%\ps@typeout{bburx is #1}
		\@bburxtrue
		\dimen100=#1
		\edef\@p@sbburx{\number\dimen100}
}
\def\@p@@sbbury#1{
		%\ps@typeout{bbury is #1}
		\@bburytrue
		\dimen100=#1
		\edef\@p@sbbury{\number\dimen100}
}
\def\@p@@sheight#1{
		\@heighttrue
		\dimen100=#1
   		\edef\@p@sheight{\number\dimen100}
		%\ps@typeout{Height is \@p@sheight}
}
\def\@p@@swidth#1{
		%\ps@typeout{Width is #1}
		\@widthtrue
		\dimen100=#1
		\edef\@p@swidth{\number\dimen100}
}
\def\@p@@srheight#1{
		%\ps@typeout{Reserved height is #1}
		\@rheighttrue
		\dimen100=#1
		\edef\@p@srheight{\number\dimen100}
}
\def\@p@@srwidth#1{
		%\ps@typeout{Reserved width is #1}
		\@rwidthtrue
		\dimen100=#1
		\edef\@p@srwidth{\number\dimen100}
}
\def\@p@@sangle#1{
		%\ps@typeout{Rotation is #1}
		\@angletrue
%		\dimen100=#1
		\edef\@p@sangle{#1} %\number\dimen100}
}
\def\@p@@ssilent#1{ 
		\@verbosefalse
}
\def\@p@@sprolog#1{\@prologfiletrue\def\@prologfileval{#1}}
\def\@p@@spostlog#1{\@postlogfiletrue\def\@postlogfileval{#1}}
\def\@cs@name#1{\csname #1\endcsname}
\def\@setparms#1=#2,{\@cs@name{@p@@s#1}{#2}}
%
% initialize the defaults (size the size of the figure)
%
\def\ps@init@parms{
		\@bbllxfalse \@bbllyfalse
		\@bburxfalse \@bburyfalse
		\@heightfalse \@widthfalse
		\@rheightfalse \@rwidthfalse
		\def\@p@sbbllx{}\def\@p@sbblly{}
		\def\@p@sbburx{}\def\@p@sbbury{}
		\def\@p@sheight{}\def\@p@swidth{}
		\def\@p@srheight{}\def\@p@srwidth{}
		\def\@p@sangle{0}
		\def\@p@sfile{} \def\@p@sbbfile{}
		\def\@p@scost{10}
		\def\@sc{}
		\@prologfilefalse
		\@postlogfilefalse
		\@clipfalse
		\if@noisy
			\@verbosetrue
		\else
			\@verbosefalse
		\fi
}
%
% Go through the options setting things up.
%
\def\parse@ps@parms#1{
	 	\@psdo\@psfiga:=#1\do
		   {\expandafter\@setparms\@psfiga,}}
%
% Compute bb height and width
%
\newif\ifno@bb
\def\bb@missing{
	\if@verbose{
		\typeout{psfig: searching \@p@sbbfile \space  for bounding box}
	}\fi
	\no@bbtrue
	\epsf@getbb{\@p@sbbfile}
        \ifno@bb \else \bb@cull\epsf@llx\epsf@lly\epsf@urx\epsf@ury\fi
}	
\def\bb@cull#1#2#3#4{
	\dimen100=#1 bp\edef\@p@sbbllx{\number\dimen100}
	\dimen100=#2 bp\edef\@p@sbblly{\number\dimen100}
	\dimen100=#3 bp\edef\@p@sbburx{\number\dimen100}
	\dimen100=#4 bp\edef\@p@sbbury{\number\dimen100}
	\no@bbfalse
}
% rotate point (#1,#2) about (0,0).
% The sine and cosine of the angle are already stored in \sine and
% \cosine.  The result is placed in (\p@intvaluex, \p@intvaluey).
\newdimen\p@intvaluex
\newdimen\p@intvaluey
\def\rotate@#1#2{{\dimen0=#1 sp\dimen1=#2 sp
%            	calculate x' = x \cos\theta - y \sin\theta
		  \global\p@intvaluex=\cosine\dimen0
		  \dimen3=\sine\dimen1
		  \global\advance\p@intvaluex by -\dimen3
% 		calculate y' = x \sin\theta + y \cos\theta
		  \global\p@intvaluey=\sine\dimen0
		  \dimen3=\cosine\dimen1
		  \global\advance\p@intvaluey by \dimen3
		  }}
\def\compute@bb{
		\no@bbfalse
		\if@bbllx \else \no@bbtrue \fi
		\if@bblly \else \no@bbtrue \fi
		\if@bburx \else \no@bbtrue \fi
		\if@bbury \else \no@bbtrue \fi
		\ifno@bb \bb@missing \fi
		\ifno@bb \ps@typeout{FATAL ERROR: no bb supplied or found}
			\no-bb-error
		\fi
		%
%\ps@typeout{BB: \@p@sbbllx, \@p@sbblly, \@p@sbburx, \@p@sbbury} 
%
% store height/width of original (unrotated) bounding box
		\count203=\@p@sbburx
		\count204=\@p@sbbury
		\advance\count203 by -\@p@sbbllx
		\advance\count204 by -\@p@sbblly
		\edef\ps@bbw{\number\count203}
		\edef\ps@bbh{\number\count204}
		%\ps@typeout{ psbbh = \ps@bbh, psbbw = \ps@bbw }
		\if@angle 
			\Sine{\@p@sangle}\Cosine{\@p@sangle}
	        	{\dimen100=\maxdimen\xdef\r@p@sbbllx{\number\dimen100}
					    \xdef\r@p@sbblly{\number\dimen100}
			                    \xdef\r@p@sbburx{-\number\dimen100}
					    \xdef\r@p@sbbury{-\number\dimen100}}
%
% Need to rotate all four points and take the X-Y extremes of the new
% points as the new bounding box.
                        \def\minmaxtest{
			   \ifnum\number\p@intvaluex<\r@p@sbbllx
			      \xdef\r@p@sbbllx{\number\p@intvaluex}\fi
			   \ifnum\number\p@intvaluex>\r@p@sbburx
			      \xdef\r@p@sbburx{\number\p@intvaluex}\fi
			   \ifnum\number\p@intvaluey<\r@p@sbblly
			      \xdef\r@p@sbblly{\number\p@intvaluey}\fi
			   \ifnum\number\p@intvaluey>\r@p@sbbury
			      \xdef\r@p@sbbury{\number\p@intvaluey}\fi
			   }
%			lower left
			\rotate@{\@p@sbbllx}{\@p@sbblly}
			\minmaxtest
%			upper left
			\rotate@{\@p@sbbllx}{\@p@sbbury}
			\minmaxtest
%			lower right
			\rotate@{\@p@sbburx}{\@p@sbblly}
			\minmaxtest
%			upper right
			\rotate@{\@p@sbburx}{\@p@sbbury}
			\minmaxtest
			\edef\@p@sbbllx{\r@p@sbbllx}\edef\@p@sbblly{\r@p@sbblly}
			\edef\@p@sbburx{\r@p@sbburx}\edef\@p@sbbury{\r@p@sbbury}
%\ps@typeout{rotated BB: \r@p@sbbllx, \r@p@sbblly, \r@p@sbburx, \r@p@sbbury}
		\fi
		\count203=\@p@sbburx
		\count204=\@p@sbbury
		\advance\count203 by -\@p@sbbllx
		\advance\count204 by -\@p@sbblly
		\edef\@bbw{\number\count203}
		\edef\@bbh{\number\count204}
		%\ps@typeout{ bbh = \@bbh, bbw = \@bbw }
}
%
% \in@hundreds performs #1 * (#2 / #3) correct to the hundreds,
%	then leaves the result in @result
%
\def\in@hundreds#1#2#3{\count240=#2 \count241=#3
		     \count100=\count240	% 100 is first digit #2/#3
		     \divide\count100 by \count241
		     \count101=\count100
		     \multiply\count101 by \count241
		     \advance\count240 by -\count101
		     \multiply\count240 by 10
		     \count101=\count240	%101 is second digit of #2/#3
		     \divide\count101 by \count241
		     \count102=\count101
		     \multiply\count102 by \count241
		     \advance\count240 by -\count102
		     \multiply\count240 by 10
		     \count102=\count240	% 102 is the third digit
		     \divide\count102 by \count241
		     \count200=#1\count205=0
		     \count201=\count200
			\multiply\count201 by \count100
		 	\advance\count205 by \count201
		     \count201=\count200
			\divide\count201 by 10
			\multiply\count201 by \count101
			\advance\count205 by \count201
			%
		     \count201=\count200
			\divide\count201 by 100
			\multiply\count201 by \count102
			\advance\count205 by \count201
			%
		     \edef\@result{\number\count205}
}
\def\compute@wfromh{
		% computing : width = height * (bbw / bbh)
		\in@hundreds{\@p@sheight}{\@bbw}{\@bbh}
		%\ps@typeout{ \@p@sheight * \@bbw / \@bbh, = \@result }
		\edef\@p@swidth{\@result}
		%\ps@typeout{w from h: width is \@p@swidth}
}
\def\compute@hfromw{
		% computing : height = width * (bbh / bbw)
	        \in@hundreds{\@p@swidth}{\@bbh}{\@bbw}
		%\ps@typeout{ \@p@swidth * \@bbh / \@bbw = \@result }
		\edef\@p@sheight{\@result}
		%\ps@typeout{h from w : height is \@p@sheight}
}
\def\compute@handw{
		\if@height 
			\if@width
			\else
				\compute@wfromh
			\fi
		\else 
			\if@width
				\compute@hfromw
			\else
				\edef\@p@sheight{\@bbh}
				\edef\@p@swidth{\@bbw}
			\fi
		\fi
}
\def\compute@resv{
		\if@rheight \else \edef\@p@srheight{\@p@sheight} \fi
		\if@rwidth \else \edef\@p@srwidth{\@p@swidth} \fi
		%\ps@typeout{rheight = \@p@srheight, rwidth = \@p@srwidth}
}
%		
% Compute any missing values
\def\compute@sizes{
	\compute@bb
	\if@scalefirst\if@angle
% at this point the bounding box has been adjsuted correctly for
% rotation.  PSFIG does all of its scaling using \@bbh and \@bbw.  If
% a width= or height= was specified along with \psscalefirst, then the
% width=/height= value needs to be adjusted to match the new (rotated)
% bounding box size (specifed in \@bbw and \@bbh).
%    \ps@bbw       width=
%    -------  =  ---------- 
%    \@bbw       new width=
% so `new width=' = (width= * \@bbw) / \ps@bbw; where \ps@bbw is the
% width of the original (unrotated) bounding box.
	\if@width
	   \in@hundreds{\@p@swidth}{\@bbw}{\ps@bbw}
	   \edef\@p@swidth{\@result}
	\fi
	\if@height
	   \in@hundreds{\@p@sheight}{\@bbh}{\ps@bbh}
	   \edef\@p@sheight{\@result}
	\fi
	\fi\fi
	\compute@handw
	\compute@resv}

%
% \psfig
% usage : \psfig{file=, height=, width=, bbllx=, bblly=, bburx=, bbury=,
%			rheight=, rwidth=, clip=}
%
% "clip=" is a switch and takes no value, but the `=' must be present.
\def\psfig#1{\vbox {
	% do a zero width hard space so that a single
	% \psfig in a centering enviornment will behave nicely
	%{\setbox0=\hbox{\ }\ \hskip-\wd0}
	%
	\ps@init@parms
	\parse@ps@parms{#1}
	\compute@sizes
	%
	\ifnum\@p@scost<\@psdraft{
		%
		\special{ps::[begin] 	\@p@swidth \space \@p@sheight \space
				\@p@sbbllx \space \@p@sbblly \space
				\@p@sbburx \space \@p@sbbury \space
				startTexFig \space }
		\if@angle
			\special {ps:: \@p@sangle \space rotate \space} 
		\fi
		\if@clip{
			\if@verbose{
				\ps@typeout{(clip)}
			}\fi
			\special{ps:: doclip \space }
		}\fi
		\if@prologfile
		    \special{ps: plotfile \@prologfileval \space } \fi
		\if@decmpr{
			\if@verbose{
				\typeout{psfig: including \@p@sfile.Z \space }
			}\fi
			\special{ps: plotfile "`zcat \@p@sfile.Z" \space }
		}\else{
			\if@verbose{
				\typeout{psfig: including \@p@sfile \space }
			}\fi
			\special{ps: plotfile \@p@sfile \space }
		}\fi
		\if@postlogfile
		    \special{ps: plotfile \@postlogfileval \space } \fi
		\special{ps::[end] endTexFig \space }
		% Create the vbox to reserve the space for the figure
		\vbox to \@p@srheight true sp{
			\hbox to \@p@srwidth true sp{
				\hss
			}
		\vss
		}
	}\else{
		% draft figure, just reserve the space and print the
		% path name.
		\if@draftbox{		
			% Verbose draft: print file name in box
			\hbox{\frame{\vbox to \@p@srheight true sp{
			\vss
			\hbox to \@p@srwidth true sp{ \hss \@p@sfile \hss }
			\vss
			}}}
		}\else{
			% Non-verbose draft
			\vbox to \@p@srheight true sp{
			\vss
			\hbox to \@p@srwidth true sp{\hss}
			\vss
			}
		}\fi	



	}\fi
}}
\psfigRestoreAt




nlin0202057.extracted_bib
\begin{thebibliography}{99}
%
\bibitem{ciss}
J. S. Anderson, I.G. Kevrekidis, R. Rico-Mart\'{\i}nez and
K. Krischer, Use of adaptation and feedback in the experimental detection
of bifurcations, {\em Proc. 30th Annual Conf. on Information Sciences and
Systems}, {\bf 1}:155--160, 1996.
%
\bibitem{prl}
J. S. Anderson, S. Y. Shvartsman, G. Fl\"atgen, I. G. Kevrekidis,
R. Rico-Mart\'{\i}nez and K. Krischer, An adaptive method for the
experimental detection of instabilities, {\it Physical Review Letters},
{\bf 82}:532-535, 1999.
%
\bibitem{astrom}
K.J. Astr\"om and B. Wittenmark, {\em Adaptive Control}, Addison-Wesley,
New York, NY, 1989.
%
\bibitem{copeland}
G. S. Copeland, I.G. Kevrekidis and R. Rico-Mart\'{\i}nez,
Adaptive detection of instabilities and nonlinear analysis of a
reduced-order model for flutter and rotating stall in turbomachinery,
{\it Proceedings of the IEEE International Conference of Control
Applications} (CCA, Hawaii), Paper TuP3-1, 1146-1150, 1999.
%
\bibitem{Doe}
E. Doedel, H.B. Keller, and J.P. Kern\'evez,
Numerical analysis and control of bifurcation problems:
(I) Bifurcation in finite dimensions, {\it Int. J. of Bifurcation
and Chaos}, {\bf 1}:493--520, 1991.
%

\bibitem{rawlings}
J. W. Eaton and J. B. Rawlings, Feedback control of nonlinear
processes using on-line optimization techniques, {\it Computers chem. Engng},
{\bf 14}:469-479, 1990.
%
\bibitem{G&S}
G.C. Goodwin and K.S. Sin, {\em Adaptive Filtering, Prediction,
and Control}, Prentice-Hall, Englewood Cliffs, NJ, 1984.
%
\bibitem{Wol}
K. Krischer, 
Principles of temporal and spatial pattern formation in electrochemical 
systems, {\it Mod. Aspects Electrochem.},
{\bf 32}:1--142, 1999.
%
\bibitem{Lyb}
D.G. O'Neil, G. Lyberatos, and S.A. Svoronos,
Adaptive Determination of Bifurcation Points, Presentation at
{\it AIChE Annual Meeting}, Paper 172j, Washington, D.C., 1988.
%
\bibitem{PTVF}
W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery,
{\em Numerical Recipes in C}, 2nd ed., Cambridge University Press,
New York, NY, 1992.
%
\bibitem{greece}
R. Rico-Mart\'{\i}nez, J.S. Anderson and I.G. Kevrekidis,
Continuous-time nonlinear signal processing: A neural network
based approach for gray box identification,
{\it Proc. 1994 IEEE Workshop on Neural Networks for Signal Processing},
596--605, 1994.
%
\bibitem{bulsari}
R. Rico-Mart\'{\i}nez, I.G. Kevrekidis and K. Krischer,
Nonlinear system identification using neural networks: dynamics and
instabilities, Chapter 16 in the book {\em Neural Networks for Chemical
Engineers}, A.B. Bulsari (ed.), Elsevier Science B.V., 1995.
%
\bibitem{rossler}
O. E. R\"ossler, 
An equation for continuous chaos, {\it Phys. Lett.}, {\bf57A}:397--398, 1976.
%
\bibitem{krstic}
H-H Wang, S. Yeung and M. Krstic,
Experimental application of extremum seeking on an axial flow compressor,
{\it IEEE Trans. on Control System Technology},
{\bf 8}:300-309, 2000.
%
\bibitem{Wig}
S. Wiggins, {\em Introduction to Applied Nonlinear Dynamical 
Systems and Chaos}, Springer-Verlag, New York, NY, 1990.
%
\bibitem{ziegler}
J. G. Ziegler and N. B. Nichols, Optimum Settings for Automatic Controllers.
{\it Trans. ASME}, {\bf 64}:759, 1942.
%

\end{thebibliography}