% -*- Mode: LaTeX; Default-character-style: (:SWISS :ROMAN :LARGE) -*-

\documentstyle[12pt,ocite]{article}

\renewcommand{\baselinestretch}{1.5}
\setlength{\oddsidemargin}{.25in}
\setlength{\evensidemargin}{.25in}
\setlength{\topmargin}{0pt}
\setlength{\textheight}{8in}
\setlength{\textwidth}{6.5in}

\newcommand{\plr}{{\sc plr}}
\newcommand{\Plr}{{\sc Plr}}
\newcommand{\qr}{{\sc qr}}
\newcommand{\Qr}{{\sc Qr}}
\newcommand{\bounder}{{\sc bounder}}
\newcommand{\Bounder}{{\sc Bounder}}

\begin{document}

\begin{titlepage}
\begin{center}
{\large\bf Qualitative Analysis by Piecewise Linear Approximation}\\

Elisha Sacks%
\footnote{This research was supported (in part) by National Institutes
of Health Grant No.\ R01 LM04493 from the National Library of
Medicine.}\\ 
MIT Laboratory for Computer Science\\
545 Technology Square, Room 370\\
Cambridge, MA 02139, USA
\end{center}

\begin{abstract}
\noindent
This paper describes a computer program called \plr\ that derives the
qualitative behavior of ordinary differential equations.  Current
qualitative reasoning programs derive the abstract behavior of a system
by simulating hand-crafted ``qualitative'' versions of the differential
equations that characterize it and summarizing the results.  \Plr\ infers
more detailed information by constructing and analyzing piecewise
linear approximations of the original equations.  The analysis employs
the phase space representation of dynamic systems theory.  \Plr\
constructs a phase diagram for a system of piecewise linear equations
by partitioning phase space into regions on which the system is linear,
analyzing the linear systems, and combining the results.  It pastes
together the local analyses into a global phase diagram by determining
which sequences of regions trajectories can traverse.  The current
implementation of \plr\ only handles second-order systems, but the
method extends to higher-order systems.  As evidence that \plr\ can
provide useful information to modelers, I present its analyses of the
Lienard and van der Pol equations.
\end{abstract}

\end{titlepage}

\section{Introduction}

Scientists and engineers model many physical systems with sets of
ordinary differential equations.  They predict the behavior of the
systems by analyzing the associated equations.  The differential
equations for typical systems contain unknown or variable parameters.
Rather than treat every instance separately, people analyze the
parameterized equations abstractly and instantiate the results with
appropriate numbers.  The analysis ends at the abstract stage if the
parameter values are unknown or irrelevant.  This approach avoids
redundancy, provides global insight, and facilitates design.  For
example, the parameterized equation $y'(t)=ay(t)$ describes the class
of one-tank devices with instantaneous mixing.  The solution,
$y(t)=y_0e^{at}$, and the physical constraint $y_0>0$ imply that $y$
increases toward infinity if $a>0$, is constant if $a=0$, and decreases
asymptotically to zero if $a<0$.  An engineer can design a specific
one-tank device by choosing an appropriate value of $a$.

This paper describes a computer program called \plr\ that derives the
qualitative behavior of second-order ordinary differential equations:
local properties in interesting regions as well as global properties such
as stability, periodicity, limit cycles, and asymptotic behavior.  (The
method extends to higher-order systems.)  For linear systems, simple
mathematical tools exist for obtaining this information.  Nonlinear
systems, however, generally require extremely sophisticated techniques
and often resist analysis altogether.  Artificial intelligence researchers
have responded to this problem by replacing differential equations with
weaker models and inference methods.  I find it more promising to
construct and analyze piecewise linear approximations of intractable
equations.  The approximations are accurate enough to reproduce the
essential properties of their nonlinear prototypes, yet simple enough to
be analyzed completely and efficiently.

The next section outlines the phase space representation.  The
following two sections describe the \plr\ methodology and demonstrate
its capabilities with two examples that frequently arise in nonlinear
oscillators: the Lienard equation and the van der Pol equation.  Both
examples are simple enough for mathematicians to analyze directly.
The solutions afford a standard against which to measure \plr.  Every
stage of piecewise analysis requires sophisticated inequality reasoning.
Section~\ref{s-bounder} presents an inequality prover that does the job
quickly and effectively.  The final two sections contain a review of
previous work and conclusions. 

\section{Phase Space}

\Plr\ uses the standard phase space representation for systems of
first-order differential equations, which stresses the long-term
qualitative form of solutions for all initial values, while downplaying
exact solutions.  It converts higher-order equations to this format by
introducing new variables as synonyms for higher derivatives.   The
phase space for a system
\begin{equation}\label{e-system}
x_i'=f_i(x_1,\ldots,x_n);\; i=1,\ldots,n
\end{equation}
is the Cartesian product of the $x_i$'s domains.  Points in phase space
represent states of the system.  Curves on which
equation~(\ref{e-system}) is satisfied, called {\em trajectories,}
represent solutions.  A phase diagram for a system depicts its phase
space and trajectories graphically.  The topological and geometric
properties of trajectories characterize the qualitative behavior of
solutions.  For example, a point trajectory, called a {\em fixed point,}
indicates a constant solution, whereas a closed curve indicates a
periodic solution.  A fixed point is {\em stable\/} if every nearby
trajectory asymptotically approaches it and {\em unstable\/} otherwise.

\section{The \plr\ Methodology}

\Plr\ analyzes intractable nonlinear equations by approximating them
with piecewise linear equations and constructing phase diagrams of
the results.  It can approximate many nonlinear systems automatically,
including both examples in this paper, by replacing nonlinear terms with
piecewise linear terms.  It approximates each nonlinear term $g(x)$ by
connecting its consecutive extrema and zeros with straight lines. For
example, it approximates the Lienard equation
\begin{equation}\label{e-lienard}
x''+x'+x^2+2x=0
\end{equation}
with the piecewise linear equations 
\begin{equation}\label{e-llienard}
\left\{\begin{array}{l}
        x'=y\\
        y'=f(x)-y
       \end{array}\right. 
\;\mbox{with}\;
f(x)=
\left\{\begin{array}{ccc}
	x+2 & \mbox{for} & x\leq-1\\
	-x & \mbox{for} & x>-1
       \end{array}\right.
\end{equation}
by introducing $y$ as a synonym for $x'$ and replacing the nonlinear
term $x^2+2x$ with two lines (Figure~\ref{f-llienard}).  The user must
approximate systems that \plr\ cannot handle.  Other analysis systems
always require their users to specify approximations.  Users of standard
numerical packages must choose appropriate algorithms, error margins,
initial guesses, and step sizes.  Similarly, de Kleer and
Brown\cite{bobrow} note that qualitative reasoning requires users to
derive the confluences for systems by themselves.  Like those systems,
\plr\ focuses on analysis, but unlike them, it also has significant model
construction capabilities. 

\begin{figure}[htbp]
\caption{Piecewise linearization of $x^2+2x$}\label{f-llienard}
\end{figure}

\Plr\ constructs a phase diagram for a parameterized piecewise linear
system in two analysis stages: local and global.  In local analysis, it
partitions phase space into regions on which the system is linear and
constructs phase diagrams for the linear systems.  It employs the
standard theory of linear equations\cite{hirsch74} to derive the phase
diagram for a linear system from the disposition in the complex plane
of its fixed point and eigenvalues.  It splits the diagram into one or
more regions on which all solutions have identical qualitative
properties.  For example, the piecewise linear Lienard
system~(\ref{e-llienard}) has two linear regions: $x\leq-1$ and $x>-1$.
The phase diagrams of the corresponding linear systems appear in
Figure~\ref{f-lphase1}.  The left-hand linear system has a saddle at
$(-2,0)$; the right-hand system has a spiral sink at $(0,0)$.  \Plr\
splits the left-hand region into subregions A through D along the
asymptotes of the saddle, denoted by dashed lines.  It leaves the
right-hand region, E, intact.

\begin{figure}[htbp]
\caption{Local phase diagrams for (a) left-hand and (b) right-hand 
           linear systems of equation~(\protect\ref{e-llienard})}
           \label{f-lphase1}
\end{figure}

Global analysis infers the joint phase diagram for a system from the
local phase diagrams through a combination of algebraic and geometric
reasoning.  It concatenates the relevant portions of the individual
diagrams and tests whether trajectories can cross the boundaries
between pairs of adjoining regions.  It summarizes the results in a
transition graph whose nodes and links represent regions and possible
transitions.  Whenever the out-degree of a region exceeds one, \plr\
splits it into subregions of degree one by case analysis.  Each
subregion consists of all points through which trajectories cross into a
specific successor.  Subregions inherit all the predecessors of their
parent region.  \Plr\ then rules out impossible transitions, as illustrated
below.  Every walk through the transition graph denotes a trajectory in
the joint phase diagram.  Loops denote trajectories that remain in one
region forever, whereas longer cycles denote trajectories that
continually shift between a sequence of regions.

Figure~\ref{f-lphase2} contains the phase diagram and transition graph
for the piecewise linear Lienard equation.   In the phase diagram, thick
lines, dashed lines, and dots represent boundaries between linear
regions, asymptotes, and fixed points.  In the transition graph, named
nodes represent regions and unnamed ones represent fixed points.  The
contents of a fixed point indicate its type: a pair of arrows pointing
inward or outward for a sink or source, one pair pointing in each
direction for a saddle, and a double circle for a center.  The node at
the top right of the figure stands for the spiral sink at the origin and
the node at the middle left for the saddle at $(-2,0)$.

Case analysis splits region~E into E1, E2, and E3 corresponding to
trajectories that remain in E forever, cross into B, and cross into C.
Trajectories in regions~C and~D cannot cross into any other region.
Trajectories in region~A cross into region~E because they increase
toward infinity in the $x$ direction.  If the height, $h$, at which a
trajectory crosses into~E is less than $h_1$, it enters E1 and remains
there forever, spiraling around the origin.  For $h$ between $h_1$ and
$h_2$, the trajectory crosses into region E2 then into region~B then
back into E1 for good.  It cannot cross from region~B to E2 or E3
because the upper boundary of~B intersects the boundary of E below
$h_1$.  Hence, \plr\ rules out these transitions.  For $h$ greater than
$h_2$, the trajectory crosses from region~A to E3 then enters
region~C and remains there.  We can verify these results by comparing
Figure~\ref{f-lphase2} with Figure~\ref{f-lphase3}, the phase diagram
for the original Lienard equation~(\ref{e-lienard}), as given by Brauer
and Nohel\cite{brauer1}. 

\begin{figure}[htbp]
\caption{(a) Transition graph and (b) phase diagram for the piecewise
           linear Lienard equation}\label{f-lphase2}
\end{figure}

\begin{figure}[htbp]
\caption{Phase diagram for the actual Lienard equation.}\label{f-lphase3} 
\end{figure}

\section{The Van der Pol Equation}
This section demonstrates \plr's analysis of the van der Pol equation,
which often arises in oscillatory dynamic systems.  Figure~\ref{f-circuit}
depicts a simple example from network theory: a capacitor, an
inductor, and a nonlinear resistor connected in series.  
\begin{figure}[b]
\caption{A circuit governed by van der Pol's equation}\label{f-circuit}
\end{figure}
By Kirchoff's laws, the current through the circuit, $I$, obeys the
equation 
\begin{equation}\label{e-vanderpol}
i''+\frac{k}{L}(i^2-1)i'+\frac{1}{LC}i=0
\end{equation}
with $C$ the capacitance, $L$ the inductance, and $k$ a positive scaling
factor.  Intuitively, the system oscillates because the nonlinear
resistor adds energy to the circuit at low currents and drains energy at
high currents.  \Plr\ obtains a piecewise linear approximation of
equation~(\ref{e-vanderpol}) by replacing the nonlinear resistor model 
$i^3/3-i$ with the piecewise linear one
\begin{equation}
f(i)=
\left\{\begin{array}{cll}
	  2(i+\sqrt{3})/3(\sqrt{3}-1) & \mbox{for} & i<-1\\
	  -2i/3 &\mbox{for} & -1\leq i\leq 1\\
	  2(i-\sqrt{3})/3(\sqrt{3}-1) &\mbox{for} & i>1
       \end{array}\right.
\end{equation}
shown in Figure~\ref{f-vanderpol}.  The analysis of the resulting
equations 
\begin{equation}\label{e-lvanderpol}
\left\{\begin{array}{l} 
          i'=y\\ 
          \displaystyle y'=-\frac{1}{LC}i-\frac{kf'(i)}{L}y
\end{array}\right.
\end{equation}
follows the general pattern described in the previous section, although
the symbolic parameters, $k$, $L$, and $C$, complicate the process
somewhat.  The fact that $f'(i)$ is undefined on the boundaries
between linear regions does not pose a problem; it suffices for the
analysis that the one-sided derivatives exist.  \Plr\ must consider two
cases, depending on whether the eigenvalues are real or complex.  I
will only discuss real eigenvalues; the complex case is similar, but
simpler.

\begin{figure}[htbp]
\caption{Piecewise linear approximation of $i^3/3-i$}\label{f-vanderpol}
\end{figure}

In the case of real eigenvalues, equation~(\ref{e-lvanderpol}) has three
linear regions: $i<-1$, $|i|\leq 1$, and $i>1$.  The outer two regions
have the same linear system.  It has negative eigenvalues, while the
middle system has positive ones.  Both have a fixed point at $(0,0)$.
Figure~\ref{f-vphase1} depicts the phase diagrams for the two linear
systems along with their regions of applicability.  As with the Lienard
equation, \plr\ infers the joint phase diagram from the individual ones
by constructing a transition graph (Figure~\ref{f-vphase2}).  The
regions are labeled A--J with regions C and E split into C1 and C2 and
E1 and E2 by the $i$-axis.  Every walk in the transition graph
terminates in the cycle: C2-D-E2-H.  This shows that all trajectories
spiral around the origin, but tells nothing of whether they move inward,
move outward, or wobble around.  In a future publication, I explain
how \plr\ proves that they converge to a unique limit cycle. 

\begin{figure}[htbp]
\caption{Local phase diagrams for the piecewise linear van der Pol
           equation (inapplicable regions shaded): (a) $|i|\leq1$ (b)
           $|i|>1$.} \label{f-vphase1} 
\end{figure}

\begin{figure}[htbp]
\caption{Phase diagram and transition graph for the piecewise linear
           van der Pol equation} \label{f-vphase2}
\end{figure}

\section{Inequality Reasoning}\label{s-bounder}

Every stage of \plr\ requires inequality reasoning: determining the signs
of eigenvalues, detecting transitions between linear regions, and ruling
out transitions in case analysis.  \Plr\ uses my \bounder\ inequality
reasoner\cite{sacks87a}.  \Bounder\ proves inequalities between
extended elementary functions over finite sets of constraints.  It can
also derive upper and lower bounds for symbolic expressions over
constraint sets.  Previous inequality algorithms perform well on some
subset of the extended elementary functions, but poorly elsewhere.
\Bounder\ maintains a hierarchy of increasingly complex algorithms to
avoid this limitation.  When one fails to resolve an inequality, it tries
the next.  This strategy resolves more inequalities than any single
algorithm.  It also performs well on hard problems without wasting time
on easy ones.  The current hierarchy consists of four algorithms:
bounds propagation, substitution, derivative inspection, and iterative
approximation.  Propagation is an extension of interval arithmetic that
takes linear time, but ignores constraints between variables and
multiple occurrences of variables.  The remaining algorithms consider
these factors, but require exponential time.  Substitution is a new,
provably correct algorithm for utilizing constraints between variables.
The final two algorithms analyze constraints created by multiple
occurrences of variables.  Inspection examines the signs of partial
derivatives.  Iteration is based on several earlier algorithms from
interval arithmetic.

\section{Previous Work}

Commonly used techniques for analyzing nonlinear differential equations
fall into the following categories: explicit solutions, dynamic systems
theory, numeric simulation, and qualitative reasoning.  Few systems
have explicit solutions.  Even when a modeler obtains an explicit
solution, he must often expend considerable effort to infer its
qualitative properties because of its complexity.  Dynamic systems
theory provides powerful tools for deriving qualitative behavior.  At
present, though, these tools apply to a limited number of systems and
require a high level of mathematical sophistication of their users.  \Plr\
invokes dynamic systems theory when possible, but relies mainly on
analysis of piecewise linear models.

Simulation yields low-level, numeric data about systems for specific
parameter values and initial conditions.  Modelers must interpret the
data and generalize the results.  This process becomes difficult for
systems containing many parameters.  Important properties of the model
can be missed in interpretation for lack of raw data.  For example,
discontinuities and extrema might fall between the observed or
simulated points, while asymptotes may exist beyond their range.  One
can also overlook important properties due to the sheer volume of raw
data.  Generalization can fail too, since a model need not behave in a
certain manner for {\em all\/} parameter values just because it does so
for certain ones.

Qualitative reasoning\cite{bobrow} (\qr) derives the abstract behavior
of dynamic systems by simulating hand-crafted ``qualitative'' versions
of their differential equations and summarizing the results.  In its
current form, \qr\ falls far short of telling modelers what they need to
know about nonlinear systems.  It can only provide extremely abstract
descriptions, such as ``the quantity $f$ increases for a while, reaches
a maximum, and decreases thereafter.''  More information is required to
analyze actual devices, physical laws, and other models:  local
properties in interesting regions such as estimates of maxima, minima,
and rates of change as well as global properties such as stability,
periodicity, limit cycles, and asymptotic behavior.  \Qr\ abstracts away
the details required to derive this information by representing dynamic
systems with confluences instead of differential equations.  It cannot
even express many important functional properties, such as linearity,
exponential decay, asymptotic approach, oscillation, damped
oscillation, stability, and limit cycles.

\Qr\ also generates spurious behaviors.  One cause, described by
Kuipers\cite{Kuipers85a}, is the local character of its analysis.  In
addition, the abstract nature of confluences introduces ambiguities that 
differential equations preclude.  For example, the equation $y'=y-y^2\/$
implies that $y'$ is negative whenever $y$ exceeds $1$, whereas the
corresponding confluence leaves the sign completely ambiguous.
Consequently, \qr\ concludes that $y$ can increase toward infinity,
even though it is bounded from above.  Kuipers\cite{Kuipers85b} notes
that this type of ambiguity crops up in almost every clinical system of
second-order or higher.  The same result holds for other domains. 

Raiman\cite{Raiman86} addresses a special case of this problem by
incorporating assertions of the form ``quantity $a$ is negligible in
relation to quantity $b$'' into \qr.  His extension does not solve our
example because neither $y^2$ nor $y$ is negligible with respect to the
other.  Nor does similar work by Mavrovouniotis and
Stephanopoulos\cite{MandM87}.

We have seen that commonly used analysis techniques for nonlinear
dynamic systems are inadequate.  Simulation yields low-level, numeric
data about individual systems, whereas qualitative reasoning focuses
on the abstract behavior of extremely general systems.  Some
researchers are interested in \qr\ as a tool for cognitive simulation, in
which case its analysis may be appropriate.  Nevertheless, it is too
weak a tool for expert-level reasoning about complex models.  \Plr\
provides an intermediate level of analysis that does the job.  It might
be possible to attain this level with an extended version of \qr\ that
included a richer set of confluences and stronger analysis algorithms.  I
find the piecewise linear approach more promising, although \plr\ takes
ideas from \qr\ as well.

\section{Conclusions}

This paper has described piecewise linear reasoning, a technique for
deriving the properties of nonlinear dynamic systems from their
piecewise linear approximations.  Although the current implementation
of \plr\ only handles second-order systems, the techniques for
constructing and analyzing piecewise linear equations carry over
directly to higher-order systems.  The purpose of \plr\ is to enable
computers to analyze dynamic systems qualitatively.  The fact that
engineers normally find piecewise linear equations adequate for
modeling devices provides empirical support for the \plr\ paradigm.  As
more concrete evidence, I have demonstrated its analyses of the
Lienard and van der Pol equations.  Although mathematicians can
handle both equations analytically, the prospect of encoding the
prerequisite creativity, sophistication, and knowledge is daunting.  \plr\
provides an algorithmic alternative. It can also help understand
nonlinear systems that defy known analytic techniques. 

\bibliography{z:>medg>bibcite>bib}
\bibliographystyle{unsrt}
\end{document}