<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">\documentclass[12pt]{article}

%=====================================================
% Include common packages, environments, counters etc
\input{common/common_latex_packages}
\input{common/common_latex_memo_format}
\input{common/common_latex_index_setup}
\input{common/common_latex_definitions}

\begin{document}

\newcommand{\ivp}{\mbox{{I}\hspace{-.008in}{v}\hspace{-.011in}{P} }}
\newcommand{\ivpp}{{I}\hspace{-.008in}{v}\hspace{-.011in}{P}}

\begin{center}

% Title:

{\LARGE The Interval Programming Model for Multi-objective \\
Decision Making}\\

\medskip

% Authors and addresses:

\footnotesize

\mbox{\large Michael R. Benjamin}\\
MIT Computer Science and Artificial Intelligence Laboratory, 32-220 \\
Cambridge MA 02139 \\
mikerb@csail.mit.edu

\normalsize

\end{center}

\noindent
The interval programming model (\ivpp) is a mathematical programming
model for representing and solving multi-objective optimization
problems.  The central characteristic of the model is the use of
piecewise linearly defined objective functions and a solution method
that searches through the combination space of pieces rather than
through the actual decision space. The piecewise functions typically
represent an approximation of some underlying function, but this
concession is balanced on the positive side by relative freedom from
function form assumptions as well as the assurance of global optimality.
In this paper the model and solution algorithms are described, and the
applicability of IvP to certain applications are
discussed.

\medskip
\noindent
{\em Key words:} branch and bound, artificial intelligence,
multi-objective optimization, behavior-based control, multiple
criteria, piecewise linear, decision making, nonlinear, convexity


%===============================================================
\section{Introduction}
%===============================================================
\label{intro}

Interval programming (\ivpp) is a model and set of algorithms for
representing and solving multi-objective optimization problems. The
basic idea, as with other mathematical programming models, is to
provide a specific set of mathematical constructs for representing a
class of problems, and a set of algorithms that generate useful
results by exploiting the nature of those constructs. The aim of this
paper is to describe the constructs used in interval programming,
along with the algorithms, and describe why this particular
arrangement is uniquely useful for a certain class of applications.  A
common thread in the applications motivating the development of IvP is
the need to fully automate a decision (as in a robot or autonomous
vehicle), or provide a decision to a human in a critical window of
time (as in a submarine control center) with little or no room for
human interaction. Furthermore, the environment is typically composed
of distinct components, or ``sub-situations'' that, in isolation, are
fairly familiar and easy to handle. Their combination however
typically creates a completely unique and unfamiliar overall situation
where the decision-making challenge is in finding the right balance
between the concerns of each sub-situation.

The central characteristic of the model is the use of piecewise
linearly defined objective functions that may represent only an
approximation of an underlying objective function. The solution method
searches through the combination space of pieces, one from each
function, rather than through the actual decision space. By performing
search using the piecewise linear functions rather than the actual
underlying function, the search algorithm is freed from function form
assumptions, and global optimality (modulo the error between the two
function representations) can be guaranteed. The problem of producing
a piecewise linear function that adequately represents the underlying
function is indeed a significant part of the overall solution process,
and the notion of adequacy may be subjective. However, the piecewise
approximation may be produced either off-line and outside the critical
decision window, or it may be done in real-time by a module that has
deep insight concerning the particular kind of function it is tasked
with approximating and can exploit this insight to create sufficiently
accurate piecewise approximations sufficiently fast. The distinct
separation of these two aspects of the solution process, the
construction of the piecewise functions, and searching through the
combination of pieces, is a central design decision that was made not
only in the interest of balancing speed and accuracy tradeoffs, but
also in the interest of providing a flexible model without function
form assumptions to accommodate an eclectic collection of objective
function producing modules.

%===============================================================
\subsection{Decision Automation and the Problem of State Space Growth}

A common method for decision making automation is to develop a policy
that associates a prescribed decision with each possible distinct
situation or state (the terms are used interchangeably
here). Automation then is primarily a matter of determining which
situation applies at the moment. However, the overhead of this method
grows exponentially as each new possible sub-situation (e.g. aspect
of the environment) is accounted for in the policy mapping situations
to actions. This problem is compounded in applications where the
objectives associated with sub-situations also vary in their relative
priority since the possible combinations of priority values
distinguish additional unique overall situations.

The problems associated with an exponential growth of state space can
be countered by dealing with separate sub-situations
independently. Rather than form a policy based on all possible
configurations in which sub-situations A and B may simultaneously come
to be, a policy is separately formed for each. This approach is more
modular and grows linearly in complexity with each new
sub-situation. The separate policies, which typically disagree in
their prescribed actions, ultimately need to be dealt with in
coordination during the final decision making process. Effective
coordination is difficult, if not impossible, if the policies do not
indicate alternatives in addition to the single best decision
associated with a situation. If a policy rates {\em all} alternative
decisions, it is, in effect, providing an objective function, and the
coordination of policies to find the best overall decision can be
characterized as solving a multi-objective optimization problem. A
primary motivation for developing the \ivp model was to use it in this
way to counter the problem of exponential state space growth.

%===============================================================
\subsection{Decision Automation and the Problem of State Space Uncertainty}

In addition to using objective functions to coordinate policies
between two or more situations that arise simultaneously, they can
also be used to coordinate policies when there is uncertainty as to
which situation aptly describes the world facing the decision maker.
Loosely speaking, in coordinating two or more {\em simultaneous}
situations the aim is to ``kill two birds with one stone'', and when
coordinating two or more {\em uncertain situations}, the aim is to
``hedge one's bet''.

Situation uncertainty is not exclusively a matter of the decision
maker's ability to sense the present state of the world without
error. It also includes uncertainty as to how the world is changing
during and shortly after the decision will be made. Any configuration of
sub-situations, each with a measure of uncertainty, is in itself a
situation, but the number of such potential combinations prohibits
developing a policy a priori for each unique overall situation. Choosing
an action that will effectively hedge bets between two potential
situations also depends on rating alternative actions in addition to the
single best action for each situation, i.e., providing objective
functions. In some cases the alternatives from each show that it is
simply not possible to hedge bets, allowing the decision maker to reap
the benefits of being decisive. Applying multi-objective optimization in
this manner was another primary motivation for the development of \ivpp.

%===============================================================
\subsection{The Scope of This Paper}

The primary focus of this paper is to clearly define what constitutes
an interval programming problem, and provide an algorithm for solving
them. Five algorithms are provided in progression, with relative
empirical results reported on ``randomly'' generated problems. The
methods to generate these random problems are presented in the Online
Supplement. Among the five algorithm versions, some are completely
dominated by others, and some may be the best for only certain types
of problems.  By ``dominated'' we mean that there is a faster
algorithm for virtually all types of conceivable problem
instances. Some of the dominated algorithms are simpler to state
however, and are provided for illustrative purposes. This paper does
not include a discussion of all currently available algorithm
versions, primarily since further versions are tied closely to
assumptions of problem type, which are tied to the specifics of an
application. This paper does not address the use of
multi-objective decision making in autonomous vehicle control systems
or in multi-objective action selection in behavior-based systems,
although it has been used in this way effectively and is a primary
motivation for the development of IvP. See \cite{bcl06}, \cite{bng06},
\cite{pir98}, \cite{rie99}, and \cite{ros97} for more on this
topic. For more on behavior based control see \cite{ark87} and
\cite{bro86}, and for examples of behavior-based control applications
see \cite{ba98}, \cite{bng06}, \cite{bl00}, \cite{srk99} and
\cite{vwlbb00}.

%===============================================================
\section{The Interval Programming Model}
%===============================================================
\label{ipdef}	

The interval programming model consists of a set of constructs for
representing multi-objective optimization problems. The primary aim
in deciding on these constructs is to strike the right balance between
constructs that allow for exploitation by a set of optimization
algorithms, and at the same time concede as little as possible in the
ability to express problems within the applications of interest. The
constructs and assumptions are set forth here, and the algorithms are
addressed in the next section.

%===============================================================
\subsection{The Decision Space}

The decision space is the set of all possible decisions that could
result from a decision making problem. Typically there are a number of
decision variables, each with their own set of possible values, and
the decision space is the set formed by the Cartesian product of each
of the variables. In the \ivp model, the assumption is made that for
each decision variable, the set of possible values is finite,
uniformly discrete, with a known upper and lower bound. Thus the
overall decision space is also finite and uniformly discrete.

The discrete assumption was considered acceptable since decision
variables typically correspond to decisions in the real world that
have limited precision in their execution. For example, if the
decision variable corresponds to an actuator controlling the heading
of robot that has a precision of $0.5$ degrees, it makes little sense
to produce a decision of $90.003$ degrees. Furthermore, with a
discrete domain, a globally optimal point can always be found in a
finite number of steps. Brute-force search may be prohibitively
expensive when searching through the overall decision space, but is
useful at times in conjunction with other methods. Proving global
optimality for a given search method can be reduced to proving how
each discrete point is implicitly ruled out from contention. The
uniform spacing assumption is purely a convenience that the algorithms
may exploit in readily knowing the value of the next higher or lower
point in the domain.

%===============================================================
\subsection{Interval Programming Functions}
\label{ivp_func}

Objective functions used in interval programming problems are
piecewise defined functions, but with a few strategic restrictions on
representation: an \ivp function is piecewise defined such that each
point in the decision space is contained in at most one piece, and
each piece is given by a set of {\em {boundary intervals}}, one for
each decision variable, and a linear {\em {interior function}}
evaluating each point in the piece. If a point is in the decision
space is not contained in any piece, it is considered infeasible, as
discussed below in Section \ref{feasible_space}.  Note that the
uniformly spaced discrete decision space discussed above does not
imply that the piecewise function has uniformly shaped or distributed
pieces.

Piecewise functions have the ability to serve as a common
representation scheme that can level the differences between modules
that produce objective functions. Each may in fact serve as an
approximation to some other underlying function, perhaps an analytical
expression or a block of software capable of evaluating a point in
the decision space. The time and process used to create sufficiently
accurate piecewise functions is outside the scope here, but discussed
in \cite{ben02a}.

%===============================================================
\subsection{Interval Programming Problems}
\label{ivpprob}

An interval programming problem consists of a collection of \ivp
functions, each with an associated priority weighting. Each function
typically corresponds to an aspiration of the decision maker, or
autonomous agent, and maps each point in the decision space to a value
that reflects the degree to which that decision supports the
corresponding aspiration. The priority weightings reflect the degree
to which the decision maker is willing to trade off achievement in one
aspiration for another based on the overall context at the moment.

For a problem defined over a decision space with $n$ decision
variables $(x_1, \ldots, x_n)$, and having $k$ objective functions
$f_1(x_1,\ldots, x_n), \ldots, f_k(x_1,\ldots, x_n)$, with $k$
priority weights $(w_1, \ldots, w_k)$, the general form is given by:

\medskip

\begin{minipage}[b]{0.85\textwidth}
\centering  
\begin{tabular}{lcl}
maximize  &amp; &amp; $w_1f_1(x_1,\ldots, x_n)+\ldots+w_kf_k(x_1,\ldots, x_n)$ \\ 
such that &amp; &amp; $f_i$ is an \ivp piecewise defined function,  \\
          &amp; &amp; $w_i$ is a positive number. \\  
\end{tabular}
\end{minipage}

\vspace{0.25in}

\noindent
The $k$ objective functions are effectively combined into a single
objective function, which begs the question as to whether or not this
constitutes ``multi-objective'' optimization. This term is applied
here to discern a subclass of single-objective optimization problems
where the single objective function to be optimized is composed of
components that are themselves meaningful objective functions.

In non-interactive applications, where the user is not engaged in a
cycle of what-ifs that alter the preference structure between
objectives, there indeed must be only one collective objective
function. For this reason, Pareto optimality (see \cite{mie99},
\cite{kr93}, or \cite{par64} for further on Pareto optimality), and
the Pareto optimal set of solutions, does not play a prominent role in
the search algorithms, only in that the solution algorithms provide a
globally optimal solution that is indeed Pareto optimal.

%===============================================================
\subsection{The Feasible Space}
\label{feasible_space}

Generally speaking, the {\em {feasible space}} is a subset of the
decision space carved out by a set of explicit constraint
constructs. No such constructs are utilized in interval programming
problems. However, constraints can be effectively embedded within an
objective function by leaving points in the decision space unmapped,
i.e., not contained in any piece, for at least one objective function.
This effect is a by-product of the search process described in the
next section.


%The use of explicit linear
%constraints in a modified form of the interval programming model is
%outside the scope here, but discussed in \cite{ben02a}.

%Might want to add references to LP or IP regarding linear constraints
%also to ground the reader (and to produce some references!!)


%Although a single point in the decision space uniquely identifies a 
%combination of $k$ pieces, the opposite is not true. 

%To qualify a candidate combination of $k$ pieces as the ``best'' solution, 
%they must share a common point in the decision space 

%Since a combination of $k$ pieces may, however, be associated with
%more than one point in the decision space, the notion of an {\em
%intersection} of two or more pieces is fundemental. The intersection
%of two pieces is a new piece that contains only those points common to
%both pieces and has a new interior function formed by the sum of the
%interior functions from each piece.

%================================================================
\section{Solving Interval Programming Problems}
%================================================================
\input{sec_alg1}

%================================================================
\subsection{Overview}

The {\em {solution space}} for a particular problem is none other than
the decision space itself. However, since each point in the feasible
decision space is contained in exactly one piece from each of the $k$
objective functions, the solution space can also be thought of as the
set of possible combinations of $k$ pieces. The search algorithms in
this section are based on this observation, using branch and bound to
search over a tree formed from the set of all such possible
combinations. Although this set can be larger than the actual decision
space, effective pruning on the corresponding search tree can greatly
reduce the number of nodes expanded, i.e., combinations
considered. Furthermore, with correct pruning of the search tree,
global optimality is preserved since all points in the feasible space
are implicitly considered.

%A set of $k$ pieces, however, may not correlate to a unique point in
%the decision space, i.e., they may share zero, one or many points. To
%compare two candidate sets of $k$ pieces, a set must map to a single
%numerical value. The definitions of piece intersection and maximal 
%value described below serve this purpose. Furthermore, 

%To make this connection between the two notions of solution space, two
%operations on pieces must be well defined: the intersection of two
%pieces, and the single maximum value of the piece.



%================================================================
\subsection{Piece Intersection and Piece Maximal Value}
\label{int_maximal_val}

Two important operations need to be well defined for comparing two
sets of $k$ pieces in the search algorithms. One is the definition of
the {\em intersection} of two pieces, and the other is the {\em
maximal value} of a piece. The intersection of two pieces $p_i$ and
$p_j$, denoted $p_i \cap p_j$, is the piece $p_k$ that has new
boundary intervals given by the pairwise intersection of each of its
intervals, and a new interior function equal to the sum of its two
interior functions. The example in Figure \ref{piece_int} illustrates
this for two pieces defined over the decision variables $x$ and $y$.

\newcommand{\pcone}{
\small
      \renewcommand{\arraystretch}{1.2}
      \begin{tabular}{|ccc|} \hline
      $4$ &amp; $\leq x \leq$ &amp; $10$ \\
      $5$ &amp; $\leq y \leq$ &amp; $40$ \\ \hline \hline 
      \multicolumn{3}{|c|}{$3x - 8y + 18$} \\ \hline 
      \end{tabular}
\normalsize}

\newcommand{\pctwo}{
\small
      \renewcommand{\arraystretch}{1.2}
      \begin{tabular}{|ccc|} \hline
      $8$ &amp; $\leq x \leq$ &amp; $20$ \\
      $1$ &amp; $\leq y \leq$ &amp; $25$ \\ \hline \hline 
      \multicolumn{3}{|c|}{$5x + 6y - 11$} \\ \hline 
      \end{tabular}
\normalsize}

\newcommand{\pcint}{
\small
      \renewcommand{\arraystretch}{1.2}
      \begin{tabular}{|ccc|} \hline
      $8$ &amp; $\leq x \leq$ &amp; $10$ \\
      $5$ &amp; $\leq y \leq$ &amp; $25$ \\ \hline \hline 
      \multicolumn{3}{|c|}{$8x - 2y + 7$} \\ \hline 
      \end{tabular}
\normalsize}

\begin{figure}[ht]
  \subfigure[one piece, $p_i$]{
    \begin{minipage}[b]{0.33\textwidth}
    \centering \pcone
    \end{minipage}}%
  \subfigure[another piece, $p_j$]{
    \begin{minipage}[b]{0.33\textwidth}
    \centering \pctwo
    \end{minipage}}%
  \subfigure[the intersection, $p_k$]{
    \begin{minipage}[b]{0.33\textwidth}
    \centering \pcint
    \end{minipage}}%
\vspace{-.2in}
\caption{The intersection of two two-dimensional pieces with 
linear interior functions. \label{piece_int}}
\end{figure}

When one or more of the intervals from $p_i$ and $p_j$ do not overlap,
$p_k$ is said to be {\em empty}, and the intersection of an empty
piece with any other piece is also an empty piece.

The maximal value of a piece, denoted ${p_i}\to${\tt{maxval}} is the
greatest value of all points contained in the piece when evaluated by
the interior function. For example, with a linear interior function,
this is obtained readily, i.e., independent of the number of points in
the piece, by examining the slope of each variable. For $p_k$ above,
$x$ has a positive slope and $y$ has a negative slope, so the function
is evaluated at $x=10$ and $y=5$ to give a maximal value of $77$. A
set of $k$ pieces can therefore be compared against another set by the
maximal value of the piece representing the intersection of all pieces
in the set. The maximal value of an empty piece is undefined.


%To qualify a candidate combination of $k$ pieces as the ``best'' solution, 
%they must share a common point in the decision space 

%Since a combination of $k$ pieces may, however, be associated with
%more than one point in the decision space, the notion of an {\em
%intersection} of two or more pieces is fundamental. The intersection
%of two pieces is a new piece that contains only those points common to
%both pieces and has a new interior function formed by the sum of the
%interior functions from each piece.

%So, to compare one combination of $k$ pieces to another, the two
%maximum values of the piece representing the intersection of the $k$
%pieces

%================================================================
\subsection{The Search Tree}
\label{search_tree}
\newcommand{\nodebox}{\texttt{nodeBox}}
\newcommand{\nodepc}{\texttt{nodePiece}}
\newcommand{\nodepiece}{\texttt{nodePiece}}
\newcommand{\ripal}{\texttt{RIPAL}}
\newcommand{\ipal}{\texttt{IPAL}}

Each level of the search tree corresponds to a single objective
function, and each node within the level corresponds to a single piece
within the objective function. Figure \ref{combo-tree} depicts the
situation for a problem with $k=3$ objective functions, each with $m$
pieces. 

\begin{figure}[ht]
  \centering 
  \includegraphics[width=5.5in]{localfigs/search_tree.png}
\caption{The search tree for $k=3$ objective functions with $m$ pieces
         each.  \label{combo-tree}} 
\end{figure}

Each leaf node corresponds to a unique combination of $k$ pieces. In
what follows, the piece associated with a particular node is referred
to as the \nodepc{}, the expression {\tt{p(i, j)}} denotes the $i$th
piece in the $j$th objective function, and {\tt {universe}} refers to
a special piece with boundary intervals spanning the entire decision
space and an interior function of $f(x_1, \ldots, x_n) = 0$. Each
node, in addition to its \nodepc, has a \nodebox, which is the
intersection of its \nodepc{} and its parent's \nodebox, with the root
node having a \nodebox{} set to the {\tt {universe}}. (For nodes on
the first level, therefore, the \nodepc{} and \nodebox{} are
equivalent in both boundary intervals and interior function.)

%----------------------------------------------------------------
\subsubsection{Full Tree Traversal}
\label{full_tree}

The algorithm in Figure \ref{ipal-v1} expands and traverses the entire
search tree. The main function, \ipal{} (Interval Programming
ALgorithm), sets up the root node and launches the tree traversal with
the call to the recursive function \ripal{}. The variable
{\tt{bestBox}}, initially null, is the best working solution, i.e.,
the intersection of the best combination of $k$ pieces. The variable
{\tt{bestVal}} is the maximal value of {\tt{bestBox}}. For simplicity
in the pseudo-code, both are assumed to globally accessible rather
than passing them as arguments to \ripal{}.  The value of
{\tt{bestVal}} cannot be meaningfully initialized, so when
{\tt{bestBox}} is null, the first non-empty leaf node will become the
working solution, regardless of its maximal value.
\begin{figure}[ht]
  \subfigure{
    \label{fig_ipalA}  %% label for first subfigure
    \begin{minipage}[b]{.35\textwidth}
    \centering \ipalA	
    \end{minipage}}%
  \subfigure{
    \label{fig_ripalA}  %% label for first subfigure
    \begin{minipage}[b]{.65\textwidth}
    \centering \ripalA	
    \end{minipage}}
\caption{\ipal{} and \ripal{} (version-1): Full search tree expansion.
\label{ipal-v1}}
\end{figure}

The first part of \ripal{} (lines 0-5) checks for the boundary, i.e.,
leaf node, condition (line 0) and potentially (lines 1-2) updates the
best solution (lines 3-4). The second part (lines 6-8) iterates
through the pieces associated with that level and builds a new node to
branch on. This algorithm performs no pruning whatsoever and will
expand all $m^k$ leaf nodes, thus ensuring global optimality.

%----------------------------------------------------------------
\subsubsection{Simple Pruning of the Search Tree}
\label{simple-prune}

The above algorithm will evolve over the coming sections, while
preserving the guarantee of global optimality. In pruning the search
tree, the time saved by pruning is balanced against the time needed to
recognize a valid pruning opportunity. The first pruning strategy
comes from slightly altering the first version of \ripal{} in Figure
\ref{ipal-v1} so that the check for an empty \nodebox{} is done at
internal nodes rather than exclusively at leaf nodes. The revised
algorithm is given in Figure \ref{ipal-v2}.

\begin{figure}[H]
    \begin{minipage}[b]{.99\textwidth}
    \centering \ripalB	
    \end{minipage}%
\caption{\ripal{} (version-2): Pruning nodes with an empty \nodebox.
\label{ipal-v2}}
\end{figure}

The basic observation is that if a node has an empty \nodebox, the
\nodebox{} for all its leaf children will also be empty. The check for
non-emptiness on line 7 replaces the same check done on line 1 in the
first version of \ripal{} in Figure \ref{ipal-v1}. In practice, this
results in an enormous amount of pruning at very little cost. Not all
pruning strategies are such clear winners. 

%----------------------------------------------------------------
\subsubsection{Empirical Results}
\label{emp_trial1}
A comparison of \ripal{}, between version-1 and version-2 is shown in
Figure \ref{trial1ab} below. Each data point represents the average
solution time of 15 randomly generated problems in 2 dimensions.
Note the difference in time scale between plots.

\begin{figure}[ht]
  \subfigure[\ivp problems with 3 objective functions.]{
    \label{trial1a}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
    \centering \includegraphics[width=3.0in]{experiments/trial1A.png}
    \end{minipage}}%
  \subfigure[\ivp problems with 4 objective functions.]{
    \label{trial1b}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
    \centering \includegraphics[width=3.0in]{experiments/trial1B.png}
    \end{minipage}}%
\caption{Comparison of \ripal: version-1 vs. version-2.  
\label{trial1ab}}
\end{figure}

Each problem tested was generated on a set of ramdomly parameterized
{\em ring functions}, described in the supplement, where the motivation
for this choice of experimentation is also discussed.  The particulars
of the two trials in Figure \ref{trial1ab}, are detailed further in
the supplement
%Section \ref{appx_trial1}, 
along with the tabular data used to
generate these plots, with high, low, and standard deviation results
for each data set representing a point in a plot. This is also true
for the empirical results reported in later sections.

Overall, two types of pruning are possible: a) pruning empty nodes as
done in the above revised algorithm, and b) pruning non-empty nodes
based on a calculated upper bound on solutions below a given
node. Pruning strategies discussed after this point utilize a grid
structure introduced below in Section \ref{ivp_grid}. Once this has
been defined, further versions of {\tt{IPAL}} and \ripal{} are
provided that are capable of further pruning.

%================================================================
\section{Grid Structures for Facilitating Search}
%================================================================
\label{ivp_grid}

%================================================================
\subsection{Overview}

An instance of a grid structure described in this section can be
associated with each objective function in an \ivp problem to
facilitate branch and bound search. It is utilized in the search
algorithms to retrieve two types of information at each internal,
i.e., non-leaf, node in the search tree. In considering whether to
prune any of the $m$ children of an internal node, a child node may be
pruned either because (a) its \nodepiece{} does not intersect the
\nodebox{} in question, or (b) an upper bound on the its leaf nodes is
inferior to the current working solution. The use of a grid structure,
introduced in this section, is central to efficient implementation of
both of the above types of pruning.

Although the grids are constructed and populated during the process of
creating a piecewise defined approximation of an underlying function,
which is outside the scope of this paper, grid configuration choices
have a significant impact on solution time. For this reason two 
grid aspects of objective function creation are discussed below, grid
initialization, and semi-uniform objective functions.


%================================================================
\subsection{The Basic Grid Structure and Contained Information}
\label{grid_build}

Grids are composed of a multi-dimensional array of ``grid elements'',
with two pieces of information associated with each element:
{\tt{gridVal}} and {\tt{gridList}}.  The latter is an initially empty
linked list. When a piece from an objective function is added to a
grid, the set of grid elements related to the piece is determined, and
then the piece is added to the linked list of each such grid
element. The value of {\tt{gridVal}} indicates the maximum value
across all pieces contained in the gridList.  When a piece from an
objective function is added to a grid element, the maximum value
within the piece (always at a corner) is first determined, and this
value is compared to the current {\tt{gridVal}}, and updated if
greater.

In the example in Figure \ref{gridex-one}, there are 23 pieces in an
example objective function, and 9 grid elements in the grid shown in
Figure \ref{gridex-three} with a constant function interiors shown for
each piece.
\begin{figure}[H]
  \subfigure[Example 2D objective function with 23 pieces]{
    \label{gridex-one}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
    \centering \includegraphics[width=2.05in]{localfigs/ipal_grid_example1.png}
    \end{minipage}}%
  \subfigure[Grid containing the function pieces]{
    \label{gridex-three}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
    \centering \includegraphics[width=2.05in]{localfigs/ipal_grid_example2.png}
    \end{minipage}}%
\caption{Example grid with 9 grid elements and 23 rectilinear pieces. \label{gridex}}
\end{figure}
The top shaded grid element in Figure \ref{gridex-three}, has the
{\tt gridList} containing pieces $\{b, c, d, g, h\}$, and the rightmost
shaded grid element has the {\tt gridList} $\{h, i, e, m, n\}$. Note that
piece $h$ belongs to both gridLists.  The {\tt gridVal} associated with the
latter grid element would be $12$, the maximum of $\{4, 12, 2, 2, 12\}$.

Each objective function may use a grid structure most suitable
for its needs. The number of pieces and distribution of pieces can be
quite different between objective functions, and therefore the number
and shape of grid elements may be quite different as well to allow
each grid structure to be efficient.  The expression {\tt grid[i]}
refers to the grid associated with objective function $i$.


%================================================================
\subsection{Initializing the Grid Element Layout}
\label{grid_init_gel}
The initialization of the grid associated with objective function $j$
takes two arguments as in:

\medskip

$${\tt{grid[j]{\to}initialize(universe, \; grid\_element\_box)}}$$

\medskip

The {\tt universe} is given by the Cartesian product of each
variable's domain. It is equivalent to the domain of the corresponding
objective function, and is also identical between objective functions
within an \ivp problem.  The second argument, {\tt
grid\_element\_box}, indicates the size of of each grid element. All
grid elements within a grid will be the same size, provided that they
divide evenly into the universe. For example, consider a 2D universe
with both variable domains given by the interval $[0, 200]$, and a 2D
{\tt grid\_element\_box} with a length of 67 in both dimensions.  The
resulting grid is shown in Figure \ref{gelex_onea} with 9 identical
grid elements.
\begin{figure}[ht]
  \subfigure[Uniform grid elements]{
    \label{gelex_onea}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
      \centering \includegraphics[width=1.85in]{localfigs/ipal_gel_example1.png}
  \end{minipage}}%
  \subfigure[Non-uniform grid elements]{
    \label{gelex_oneb}  %% label for second subfigure
    \begin{minipage}[b]{.50\textwidth}
      \centering \includegraphics[width=1.85in]{localfigs/ipal_gel_example2.png}
  \end{minipage}}
  \caption{Two different universes result in different 
    layouts with the same requested grid element size.}
  \label{gelex_one}
\end{figure}

If the universe domains were instead $[0, 180]$, the higher indexed
grid elements would be truncated, resulting in a non-uniform grid, as
shown in Figure \ref{gelex_oneb}.  The uniformity of all but possibly
the highest indexed grid elements is a key property that allows for an
efficient determination of which grid elements intersect a query box.
In the example in Figure \ref{gelex_oneb}, the 9 grid elements could
have been shaped in a more equitable way, which generally speaking,
results in a more efficient grid.  However, the initialization
function respects the requested {\tt grid\_element\_box} since it
typically is coordinated with the shape of the objective function
pieces it contains, as discussed next.

%================================================================
\subsection{Coordinating a Grid with a Semi-uniform Objective Function}
\label{grid_coord}

Semi-uniform \ivp functions, are created by allocating an initial
portion of the total number of pieces to represent the underlying
function in a uniform piecewise manner. This idea is depicted in
Figures \ref{gridex_fivea} and \ref{gridex_fiveb}. The piece
distribution in Figure \ref{gridex_fivea} is perfectly uniform, using
only 9 pieces initially, with the remaining 53 pieces created by
splitting the original 9 pieces and their descendents, as shown in
Figure \ref{gridex_fiveb}. 
\begin{figure}[ht]
  \subfigure[Uniform 2D \ivp function]{
    \label{gridex_fivea}  %% label for first subfigure
    \begin{minipage}[b]{.33\textwidth}
    \centering \includegraphics[width=1.7in]{localfigs/ipal_grid_example5A.png}
    \end{minipage}}%
  \subfigure[Semi-uniform \ivp function]{
    \label{gridex_fiveb}  %% label for second subfigure
    \begin{minipage}[b]{.33\textwidth}
    \centering \includegraphics[width=1.7in]{localfigs/ipal_grid_example5B.png}
    \end{minipage}}%
  \subfigure[\ivp function aligned with grid]{
    \label{gridex_fivec}  %% label for third subfigure
    \begin{minipage}[b]{.33\textwidth}
    \centering \includegraphics[width=1.7in]{localfigs/ipal_grid_example5C.png}
    \end{minipage}}%
\caption{Aligning the grid with the initial uniformity of an \ivp function.
\label{gridex_five}}
\end{figure}
The latter non-uniform pieces are allocated strategically to areas of
the function domain that may need more pieces to adequately represent
the underlying function. Grid efficiency is improved if the grid
element layout is then aligned with the initial uniform pieces in the
\ivp function (Figure \ref{gridex_fivec}).

The alignment ensures that no piece belongs to more than one grid
element, and thus ensures the smallest possible {\tt gridList}
associated with each grid element. The alignment also improves the
accuracy of the {\tt gridVal} associated with each element, especially
for pieces with linear interior functions. When a piece is added to
the grid, the upper bound for the grid element is based on the maximal
value of the linear interior function of the piece. If a portion of
the piece lies outside the grid element, then the bound may not be
accurate. This situation is eliminated when no piece lies in more than
one grid element.

For all empirical results reported in this paper, objective functions
are semi-uniform, and aligned with the grid (see the supplement
for more on the degree of semi-uniformity for specific trials). The
use of a semi-uniform function allows the contributor of each
objective function to retain significant latitude in deciding what
configuration and amount of pieces is sufficient to represent its
underlying function. In certain application domains, the source of
each objective function can vary significantly, and this latitude may
be important. If, however, it makes sense to further coordinate the
objective functions within a problem to have identical grid element
layouts and initial uniformity, the search time may be dramatically
reduced since each piece is prone to intersect with less pieces from
other objective functions, thus reducing the branching factor of the
search tree.

%================================================================
\section{Intersection Pruning Using \ivp Grid Structures}
%================================================================

%================================================================
\subsection{Overview}
For a given internal node in the search tree, any of its $m$ children
may be pruned if the child's \nodepiece{} does not intersect the
\nodebox, as discussed in Section \ref{simple-prune}. Although it is
possible to simply check each of the $m$ children to see if they
intersect the \nodebox{} (as in Figure \ref{ipal-v2}), in practice,
the great majority do not intersect.  For a given internal node, a
grid structure associated with the next level is used to return a list
of pieces that are ``close'' to a given \nodebox{} and therefore
``probably'' intersect the \nodebox.  A full intersection test is
reserved only for pieces on this list, and any consideration
whatsoever for the other ``non-close'' pieces is avoided.

%================================================================
\subsection{Retrieving Intersection Information From a Grid Structure}
\label{grid_int}

A grid structure may be queried to determine which of the stored boxes
intersect with a provided {\em{query box}}. The syntax for this query
is:

$${\tt{boxset = grid[i]{\to}getBoxes(qbox)}}$$

\medskip

\noindent
where {\tt grid[i]} is the grid associated with objective function $i$,
the argument {\tt{qbox}}{} is a query box, and {\tt{boxset}} is the resulting set (linked
list) of boxes intersecting the query box.  In the example depicted in
Figure \ref{gridex-five}, the query box intersects one grid element
containing boxes $\{h,i,e,m,n\}$. Individual pairwise intersection tests
are performed on these five, with the returned
list being $\{m,n\}$.
\begin{figure}[H]
  \subfigure[Query box intersects one grid element]{
    \label{gridex-five}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
    \centering \includegraphics[width=1.8in]{localfigs/ipal_grid_example3.png}
    \end{minipage}}%
  \subfigure[Query box intersects two grid elements]{
    \label{gridex-six}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
    \centering \includegraphics[width=1.8in]{localfigs/ipal_grid_example4.png}
    \end{minipage}}%
\caption{Retrieving intersection information from a grid given a query box.
\label{gridex-query}}
\end{figure}
In the second case, in Figure \ref{gridex-six}, the query box
intersects two grid elements. The concatenation of lists from these
two grid elements is $\{g,h,k,l,h,i,e,m,n\}$, which contains piece $h$
twice. Thus duplicates must be removed from this list, and, in
practice, this is done before the pairwise intersection tests are
performed.  In this simple case, the cost in removing duplicates is
negligible, but in problems with more pieces, and higher dimensions
(with usually coarser grids), the removal of duplicates may consume a
significant amount of time. Thus grids that ensure that no piece lies
in more than one grid element (discussed in Section \ref{grid_coord})
may have a noticeably more efficient implementation of this function by
obviating the check for duplicates.

%================================================================
\subsection{Search Tree Pruning Using Grid Intersection Information}
\label{ipal_with_grid}

The versions of {\tt IPAL} and {\tt RIPAL} given in Figure
\ref{ipal-v3} use the grid associated with each of the $k$ objective
functions to avoid as many of the $m$ intersection tests as possible
at each node.  The grid associated with each objective function is
assumed to have been a priori constructed by the modules that produced
the objective function.  As before there are $k$ objective functions,
and $m$ pieces in each function.
\begin{figure}[ht]
  \subfigure[IPAL (unchanged)]{
    \label{fig_ipalB}  %% label for first subfigure
    \begin{minipage}[b]{.35\textwidth}
      \centering \ipalA
  \end{minipage}}%
  \subfigure[RIPAL (version-3)]{
    \label{fig_ripalD}  %% label for first subfigure
    \begin{minipage}[b]{.65\textwidth}
      \centering \ripalD
  \end{minipage}}
  \caption{\ipal{} and \ripal{} (version-3): Utilizing
    intersection information from grid structures.}
  \label{ipal-v3}
\end{figure}

In this third version of \ripal, lines 0-4 remain unchanged from
Figure \ref{ipal-v2}. Lines 7-10 now reflect that iteration through a
linked list is conducted rather than the array of pieces of fixed
length $m$ as before. This linked list is the subset of the $m$ pieces
in the objective function at {\tt{level+1}} determined, by a call to
the grid in Line 5, to intersect the {\tt nodeBox}. An empty returned
{\tt boxset}, or the end of the linked list is detected when {\tt
ibox$=$NULL}. 

%================================================================
\subsection{Empirical Results}
\label{emp_trial2}

A comparison of \ripal{}, between version-2 and version-3 is shown in
Figure \ref{trial2ab} below.  Each data point represents the average
solution time of 15 randomly generated problems in 2 dimensions. This
trial is similar to the previous trial shown in Figure \ref{trial1ab}
in that the x-axis shows a progression in the number of pieces in each
objective function. Note that, as well as version-2 performed,
relatively speaking, in the previous trial, it is completely dominated
by version-3, which uses the grid structure to obviate many
intersection tests. All solution times using version-3 and later
algorithm versions using grids do not include the time taken to
populate the grid data structures. Note the shape of the version-2
curve, and that, without the use of the grid to obviate intersection
tests, $m$ tests are needed at each node, and the algorithm has
$\Omega{(m^2)}$ performance.

\begin{figure}[H]
  \subfigure[\ivp problems with 2 objective functions.]{
    \label{trial2a}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
    \centering \includegraphics[width=3.0in]{experiments/trial2A.png}
    \end{minipage}}%
  \subfigure[\ivp problems with 3 objective functions.]{
    \label{trial2b}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
    \centering \includegraphics[width=3.0in]{experiments/trial2B.png}
    \end{minipage}}%
  \caption{Comparison of \ripal: version-2 vs. version-3. 
  \label{trial2ab}}
\end{figure}

The particulars of the two trials in Figure \ref{trial2ab}, are
detailed further in the supplement
%Section \ref{appx_trial2}, 
along with the tabular data used to generate these plots, with high,
low, and standard deviation results for each data set representing a
point in a plot.

%================================================================
\section{Upper Bound Pruning Using \ivp Grid Structures}
%================================================================
\label{grid_bound}

%================================================================
\subsection{overview}

Grid structures are also used to quickly obtain an upper bound on
solution for a particular node before any of its children are tested
for intersection.

%================================================================
\subsection{Retrieving Upper Bound Information From a Grid Structure}
\label{bound_alg}

The second primary use of the grid structures is to find an upper
bound associated with a particular search node. Recall the situation
shown in Figure \ref{combo-tree}, where a particular node has $m$
children.  If an upper bound on all leaf children is shown to be
poorer than a previously found solution, backtracking can be invoked
immediately, precluding any further intersection tests and branching
below the current node.  

The tighter the bound, the more likely that pruning opportunities will
be identified and taken advantage of. However, typically, the better
the quality of the bound, the more expensive the cost of each node
that is actually expanded. The bounding method presented here will be
referred to with the following syntax:

$${\tt{bound = grid[i]{\to}getBound(qbox, covered)}},$$

\medskip

\noindent
where, as before, {\tt{qbox}} is a query box, and {\tt{i}} is an index
indicating a particular objective function.  The argument {\tt
covered} is a boolean return value set to be true if the query box
intersects at least one grid element with a non-empty {\tt gridList},
and false otherwise.  The determination of the upper bound is given by
the greatest {\tt{gridVal}} of those grid elements intersecting the
query box. In the example in Figure \ref{upbound-fig} below, the query
box intersects the two shaded grid elements, which have a
{\tt{gridVal}} of $22$ and $91$ respectively, thus giving a upper
bound of $91$.

\begin{figure}[H]
  \begin{minipage}[c]{.99\textwidth}
    \centering
    \includegraphics[height=2.25in]{localfigs/ipal_grid_example6.png}
  \end{minipage}%
  \caption{An example set of boxes from an objective function, with
the {\tt{maxval}} shown for each box.}
  \label{upbound-fig}
\end{figure}

This bound is extremely quick to obtain because the value associated
with each grid element is set once, when the grid is built and
populated. For a problem with $k$ objective functions, an upper bound
for a node at level $i$ is the sum of the {\tt{maxVal}} of the
{\tt{nodeBox}} and the upper bound on the portion of the tree below
level $i$ given by:

$$\sum_{j=i+1}^k {\tt{grid[j]{\to}getBound(qbox, covered)}},$$ 
\medskip

\noindent
where the {\tt{qbox}} is the {\tt{nodeBox}} of the node in question.
Note that the bound is progressively more likely to be tighter as the
size of the query box shrinks. The {\tt{nodeBox}} does indeed tend to
shrink deeper in the tree, but this is tempered by the fact that
pruning is also less rewarding deeper in the tree.

%================================================================
\subsection{Search Tree Pruning Using Grid Upper Bound Information}

To utilize the bounding information from the grid structures,
version-4 of \ripal{} is listed below in Figure \ref{ipal-v4}. The
primary difference is the addition of lines 9-15 between lines 8 and 9
of version-3 listed in Figure \ref{ipal-v3}.  
\begin{figure}[ht]
  \begin{minipage}[b]{.99\textwidth}
    \centering \ripalEEE
  \end{minipage}%
  \caption{A version of {\tt RIPAL} utilizing the bounding information
    from the grid of each objective function.}
  \label{ipal-v4}
\end{figure}
In determining the upper bound, rather then simply iterating through
all lower levels in lines 12-14, a while loop is used to abort the 
bounding process if the {\tt newNodeBox} is not covered by one of the
objective functions.

%================================================================
\subsection{Empirical Results}
\label{emp_trial3}

A comparison of \ripal{}, between version-3 and version-4 is shown in
the three pairs of trials in Figures \ref{trial3ab} thru
\ref{trial3ef} below. In the first trial pair, in Figure
\ref{trial3ab}, each data point represents the average solution time
of 60 randomly generated problems with 10 objective functions, first
in (a) 2 dimensions, and then in (b) 4 dimensions.  The x-axis
shows a progression in the number of pieces in each of the objective
functions, with a limit of 30,000 pieces versus the 20,000 limit in
Figure \ref{trial2ab}. In this case, however, the search tree is
significantly larger primarily due to the ten objective functions in
each test problem. With 10,000 pieces, there are $10^{40}$ possible
leaf nodes.

In this trial set, with 10 objective functions in each problem, the
number of local optima is likely to be greater, and the solution time
variance between similar problems was observed to be much higher.  See
the tables in the supplement for variance data.  For this
reason, the test set size was chosen to be large enough to see a
largely monotonic growth in time as the number of pieces
increase. This trial pair shows the effectiveness of using the bounding
technique described in the previous section. For problems with 30,000
pieces, version-4 is roughly 17 times faster in the 2D trials, and 12
times faster in the 4D trials. This ratio of relative effectiveness
appears to grow as the number of pieces grow.


\begin{figure}[H]
  \subfigure[\ivp problems with 2 dimensions, 10 objective functions.]{
    \label{trial3a}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
      \centering \includegraphics[width=3.0in]{experiments/trial3A.png}
  \end{minipage}}%
  \subfigure[\ivp problems with 4 dimensions, 10 objective functions]{
    \label{trial3b}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
      \centering \includegraphics[width=3.0in]{experiments/trial3B.png}
  \end{minipage}}%
  \caption{Comparison of \ripal: Version-3 versus Version-4.}
  \label{trial3ab}
\end{figure}

\noindent
In the second trial set, in Figure \ref{trial3cd}, the test files are
nearly identical to the previous set, but have 20 objective functions
each instead. As with the problems in the first trial set with 10
objective functions, the effectiveness of the bounding technique used
in version-4 is apparent. However, the relative effectiveness is not
as stark as with the problems with only 10 objective functions.  For
problems with 30,000 pieces, version-4 is roughly 13 times faster in
the 2D trials, and 3 times faster in the 4D trials. In fact, for 4D
problems with 6,000 pieces or less, version-3 is actually faster! This
crossover is easier to discern in the 
supplement, in the table 
%table \ref{table3c}, in the appendix,
corresponding to Figure \ref{trial3d}.

\begin{figure}[H]
  \subfigure[\ivp problems with 2 dimensions, 20 objective
  functions.]{
    \label{trial3c}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
      \centering \includegraphics[width=3.0in]{experiments/trial3C.png}
  \end{minipage}}%
  \subfigure[\ivp problems with 4 dimensions, 20 objective functions.]{
    \label{trial3d}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
      \centering \includegraphics[width=3.0in]{experiments/trial3D.png}
  \end{minipage}}%
  \caption{Comparison of \ripal: Version-3 versus Version-4.}
  \label{trial3cd}
\end{figure}

\noindent
The diminished relative effectiveness between version-3 and version-4
between the two trial sets can be accounted for in the bounding
algorithm. The accuracy of the bound calculated at each search node
is not only diminished, but also more expensive to calculate as the
number of objective functions grow. At some point, the overhead of
attempting to bound the solution does not justify the benefits of
the pruning resulting from the bounding. 

To further demonstrate this relationship, the problems in the third
pair of trials, in Figure \ref{trial3ef}, vary instead in the number
of objective functions while holding the number of pieces constant at
4,000. In each trial pair there is a cross-over point where version-4,
with bounding, no longer out-performs version-3. In the 2D trial, this
roughly occurs in problems with 50 or more objective functions. In the
4D trial, this happens much sooner, at roughly 20 objective functions.
(Note the data points in each plot at 20 objective functions. These
points can also be seen in the trial pair in Figure \ref{trial3cd},
where the number of pieces equals 4,000.)

\begin{figure}[H]
  \subfigure[\ivp problems with 2 dimensions, 4000 pieces.]{
    \label{trial3e}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
      \centering \includegraphics[width=3.0in]{experiments/trial3E.png}
  \end{minipage}}%
  \subfigure[\ivp problems with 4 dimensions, 4000 pieces.]{
    \label{trial3f}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
      \centering \includegraphics[width=3.0in]{experiments/trial3F.png}
  \end{minipage}}%
  \caption{Comparison of {\tt RIPAL}, version-3 vs. version-4 as the 
      number of objective functions grow.}
  \label{trial3ef}
\end{figure}

These results indicate the likely advantage of mixing the strategy of
bounding, e.g., not bounding at all nodes as in version-4, and not
forgoing all bounding as in version-3. Mixing strategies is outside
the scope of this paper. The particulars of the three pairs of trials
in Figures \ref{trial3ab} thru \ref{trial3ef}, are detailed further in
the supplement
%Section \ref{appx_trial3}, 
along with the tabular data used to generate these plots, with high, low, and
standard deviation results for each data set representing a point in a
plot.

%================================================================
\section{Solving IvP Problems Using an Initial Solution}
%================================================================

%================================================================
\subsection{Overview}

An application of the branch and bound search method typically stands
to benefit from the use of some heuristic that provides an initial
solution that is likely to be pretty good for some reason or another.
By measuring this solution prior to the traversal of the search tree,
the hope is that the initial best-so-far solution is high enough that
bounding and pruning will be effective earlier in the search tree.
Otherwise, with no initial solution, it would suffer through an
initial progression of weak and moderate solutions before the better
solutions would allow the bounding methods to identify pruning
opportunities higher up in the tree where pruning is most rewarding.

To accommodate an initial solution, the algorithm {\tt IPAL} in Figure
\ref{fig_ipalA} and \ref{fig_ipalB}, is altered as follows:

\medskip
\begin{figure}[ht]
  \begin{minipage}[b]{.99\textwidth}
    \centering \ipalB	
  \end{minipage}
  \caption{A revised version of \ipal{} with an ``initial solution'' 
    parameter.}
  \label{ipalisol}
\end{figure}

\noindent
The parameter {\tt initialBestBox} could be provided by any one of
possible heuristic algorithms, like the one provided below. If this
parameter is {\tt NULL}, the branch and bound search method simply
proceeds as normal, with the first best-so-far solution being the
first leaf node representing a combination of intersecting pieces from
the $k$ objective functions.


%================================================================
\subsection{An Initial-Solution Heuristic}

For a problem with $k$ objective functions, the initial-solution
heuristic described here will evaluate $k$ points in the decision
space, where each point is the maximal point associated with one of
the objective functions. To obtain the latter information, the grid
population method described in Section \ref{grid_build} is augmented
slightly by keeping a running indicator of the maximal point of any
piece in the grid associated with an objective function as it is
populated. The maximal value of each piece is compared to this as each
piece is added. The expression for obtaining this maximal point from
the grid associated with objective function $i$ is:

\medskip
$${\tt{box = grid[i]{\to}getMaxPoint()}}.$$
\medskip

\noindent
Given a maximal point for one of the $k$ objective functions, the
heuristic will need to evaluate this point from the perspective of the
other objective functions. The expression for evaluating a point for
an objective function $i$ is:

\medskip
$${\tt{val = grid[i]{\to}evaluatePoint(box, covered)}}.$$
\medskip

\noindent
The argument, {\tt covered}, is a boolean return value set to be true
if the query box intersects at least one grid element with a non-empty
{\tt gridList}, and false otherwise.  As mentioned in Section
\ref{feasible_space}, if a point in the decision space is not
contained within a piece of a particular objective function, this
indicates that the point is infeasible. When a box is not covered, the
numerical return value, {\tt val}, is undefined.

The initial-solution heuristic, utilizing the above two expressions,
is provided below in Figure \ref{isolalg}.
\begin{figure}[ht]
  \begin{minipage}[b]{.99\textwidth}
    \centering \isolone
  \end{minipage}
  \caption{A heuristic for (usually) finding a good initial solution.} 
  \label{isolalg}
\end{figure}
It is possible that this heuristic returns a null solution if all
points individually optimal by one of the $k$ objective functions are
infeasible, i.e., uncovered, from the perspective of one of the other
$k-1$ functions. For the purposes of testing the effectiveness of this
heuristic in the following empirical tests, its use prior to the
invocation of {\tt IPAL(initialSolution)}, will be referred to as
``version-5'' of the search algorithm.

%================================================================
\subsection{Empirical Results}
\label{emp_trial4}

A comparison of \ripal{}, between version-4 and version-5 is shown in
Figures \ref{trial4ab} thru \ref{trial4ef} below. In the first pair
of trials shown in Figure \ref{trial4ab}, the test problems are
the same used for comparing version-3 with version-4 reported
previously in Figure \ref{trial3ab}. Each data point represents the
average solution time of 60 randomly generated problems with 10
objective functions, first in (a) 2 dimensions, and then in (b) 4
dimensions, with the x-axis showing a progression in the number of
pieces in each of the objective functions.

\begin{figure}[H]
  \subfigure[\ivp problems with 2 dimensions, 10 objective functions.]{
    \label{trial4a}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
      \centering \includegraphics[width=3.0in]{experiments/trial4A.png}
  \end{minipage}}%
  \subfigure[\ivp problems with 4 dimensions, 10 objective functions.]{
    \label{trial4b}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
      \centering \includegraphics[width=3.0in]{experiments/trial4B.png}
  \end{minipage}}%
  \caption{Comparison of \ripal: Version 4 versus Version 5 as the number
of pieces grow.
  Trials were run on the same problems reported in Figure \ref{trial3ab}.}
  \label{trial4ab}
\end{figure}

\noindent
In the second pair of trials shown in Figure \ref{trial4cd}, the test
problems are the same used for comparing version-3 with version-4
reported previously in Figure \ref{trial3ef}. Each data point
represents the average solution time of 80 randomly generated problems
with 4000 pieces in 2 dimensions, and 18 randomly generated problems
in 4 dimensions, with the x-axis showing a progression in the number
of objective functions in each problem.


\begin{figure}[H]
  \subfigure[\ivp problems with 2 dimensions, 4000 pieces.]{
    \label{trial4c}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
      \centering \includegraphics[width=3.0in]{experiments/trial4C.png}
  \end{minipage}}%
  \subfigure[\ivp problems with 4 dimensions, 4000 pieces.]{
    \label{trial4d}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
      \centering \includegraphics[width=3.0in]{experiments/trial4D.png}
  \end{minipage}}
  \caption{Comparison of \ripal: version-4 versus version-5 as the 
    number of objective functions grow. Trials were run on the same 
    problems reported in Figure \ref{trial3ef}.}
    \label{trial4cd}
\end{figure}

\noindent
These results indicate that the use of this initial-solution heuristic
is advantageous across a wide variety of problem types ranging in not
only piece count, but also in the number of objective functions. In
the last pair of trials shown in Figure \ref{trial4ef}, version-5 is
compared with version-3 instead, on the same problems reported in
Figure \ref{trial3ef}. With the use of the initial-solution heuristic,
the cross-over point where version-3 outperforms version-5 happens at
roughly 60 objective functions, versus 50 previously. The effect in
the 4D case is not as apparent.

\begin{figure}[H]
  \subfigure[\ivp problems with 2 dimensions, 4000 pieces.]{
    \label{trial4e}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
      \centering \includegraphics[width=3.0in]{experiments/trial4E.png}
  \end{minipage}}%
  \subfigure[\ivp problems with 4 dimensions, 4000 pieces.]{
    \label{trial4f}  %% label for first subfigure
    \begin{minipage}[b]{.50\textwidth}
      \centering \includegraphics[width=3.0in]{experiments/trial4F.png}
  \end{minipage}}%
  \caption{Comparison of \ripal: version-3 versus version-5 as the 
    number of objective functions grow. Trials were run on the same 
    problems reported in Figure \ref{trial3ef}.}
  \label{trial4ef}
\end{figure}

\noindent
The particulars of the three pairs of trials in Figures \ref{trial4ab}
thru \ref{trial4ef}, are detailed further in 
the supplement 
%Section \ref{appx_trial4}, 
along with the tabular data used to generate these
plots, with high, low, and standard deviation results for each data
set representing a point in a plot.

%================================================================
\section{The Effect of Dimension Size on the Solution Algorithms}
%================================================================

The versions of \ipal{} and \ripal{} presented, up through version-5,
make no reference to {\em dimension}, i.e., the number of decision
variables over which the objective functions are defined. (True not
only in pseudo-code but the also in the C++ code used to generate
empirical results.) The issue of problem dimension is embedded in the
implementation of the two key operations, intersection and maximal
value, described in Section \ref{int_maximal_val}. Grids have also
been generalized to handle objective functions with arbitrary
dimension.

Although the algorithms make no explicit mention of {\em dimension},
trials 3F and 4F in Figures \ref{trial3ef} and \ref{trial4ef}, show
that the 4D problems take longer to solve than their 2D
counterparts. This is due to more than one issue. The simplest issue
is that the core operations of intersection and finding the maximal
value of a piece contain a number of steps that grows linearly in the
number of variables.  Broadly speaking, there are at least four
influences on overall run-time of the search algorithm:

\begin{itemize}
\item
{\em The number of nodes}: The number of nodes with a non-empty {\tt
nodeBox}. See Section \ref{search_tree}. This is primarily determined
by the number of objective functions, and the number of pieces in each
function, but also the average number of pieces from one objective
function that intersect a given piece from another.  In functions
defined over a larger domain, determined in part by dimension, the
number of pieces needed to accurately represent the underlying
function may increase.

\item
{\em The cost of the intersection and maxval operations}: The basic
tree traversal operations. See Section \ref{int_maximal_val}. The number
of steps in each grows linearly in the number of dimensions. 

\item
{\em The cost of executing the grid functions}: The primary grid
functions are {\tt getBoxes} and {\tt getBound} described in Sections
\ref{grid_int} and \ref{bound_alg} respectively. The number
of steps in each grows linearly in the number of dimensions, for grids
with identical numbers of grid elements.

\item
{\em The accuracy of the grid function results}: In functions defined
over a larger domain, determined in part by dimension, each grid will
will be defined over a greater expanse of the domain. For the
function, {\tt getBound}, this means that the upper bound will likely
deviate more from, i.e., be less representative of, pieces not near
the piece contributing to the upper bound. For the function {\tt
getBoxes}, this means that the query box may be tested for
intersection with a greater number of pieces that don't actually 
intersect.

\end{itemize}

\noindent
For these reasons, the empirical results, as dimensions vary, are not
as conslusive. A full study of this issue is outside the scope of this
paper.

%================================================================
\section{Conclusions}
%================================================================
\label{summary}

The interval programming model is a mathematical programming model for
representing and solving multi-objective optimization problems. The
central characteristic of the model is the use of piecewise linearly
defined objective functions that may represent only an approximation
of an underlying objective function. The solution method searches
through the combination space of pieces, one from each function,
rather than through the actual decision space. By performing search
using the piecewise linear functions rather than the actual underlying
function, the search algorithm is freed from function form
assumptions, and global optimality (modulo the error between the two
function representations) can be guaranteed. In this paper, the model
was defined, and a progression of solution algorithms was
provided. Two of the algorithms outperform all the others on all test
problems, and outperform each other on certain test problems.

The problem of producing a piecewise linear function that adequately
represents the underlying function is indeed a significant part of the
overall solution process, and the notion of adequacy may be
subjective. However, the piecewise approximation may be produced
either off-line and outside the critical decision window, or it may be
done in real-time by a module that knows about the particular kind of
function it is tasked with approximating and uses this insight to
create sufficiently accurate piecewise approximations sufficiently
fast. The issue of how such functions are created in particular
applications was outside the scope of this paper. 

The distinct separation of these two aspects of the solution process,
the construction of the piecewise functions, and searching through the
combination of pieces, is a central design decision that was made not
only in the interest of balancing speed and accuracy tradeoffs, but
also in the interest of providing a flexible model without function
form assumptions to accommodate an eclectic collection of objective
function producing modules. Work is ongoing to further the
capabilities of the algorithms to solve larger and more widely varied
classes of problems, as well as work motivated by particular
applications that focus on a specific subclass of problems with
features that can be exploited, such as problems with very large sets
of objective functions in low dimensions. 

\section*{Acknowledgements}

The author gratefully acknowledges the Office of Naval Research (ONR)
for support of this work, in particular Dr. Don Wagner ONR 311, and
the In-House Laboratory Independent Research (ILIR) program through
ONR and the Naval Undersea Warfare Center, Division Newport RI. The
author also gratefully acknowledges the influence of Pascal van
Hentenryck, John Savage and Philip Klein of the Brown University
Computer Science Department in Providence Rhode Island.


\bibliographystyle{plain} 
\bibliography{biblio}

%\input{sec_appendix}

\end{document}
</pre></body></html>