\documentstyle[IEEEtran]{article}
\begin{document}
\tolerance 10000

\pagestyle{myheadings}

\markright{APIC'95, El Paso,
Extended Abstracts,
A Supplement to the international journal of {\rm Reliable
Computing}\ \ \ \ \ \ \ \ \ \ \  \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ }

\title{Adding Interval Constraints to the Moore-Skelboe Global
Optimization Algorithm}
\author{H.M. Chen and M.H. van Emden}

\date{}

\maketitle

\begin{abstract}
Interval constraints are a promising addition to global optimization
with interval analysis \cite{rtrkn88,hnsn92}. We explain the
principle of interval constraints and the application of this method
to global
optimization. We then show how the Moore-Skelboe global optimization
algorithm can be made more efficient by adding interval constraints.
This yields a simple algorithm in BNR Prolog \cite{BNR88,bnldr94}
and achieves a factor of
two to five
in function evaluations on the unconstrained optimization
examples of \cite{rtrkn88}.
\end{abstract} 
 

\auffil{The authors are with the University of Victoria, Victoria,
B.C.,
Canada. The authors gratefully acknowledge support from the Natural
Science and Engineering Research Council of Canada.}

\section{Introduction}

Although interval constraints, as a method for solving numerical problems,
has much in common with interval arithmetic, there are some interesting
differences between the two approaches.
In this paper we show that one of these differences
can be used to improve the Moore-Skelboe algorithm
 \cite{rtrkn88} for
global optimization. This algorithm is unique in nonlinear optimization
by giving an interval for a global minimum while not requiring
differentiability, or even continuity, in the objective function.

We briefly introduce interval constraints and explain by means
of an example how this
method applies to global optimization. We then present our modification
of the Moore-Skelboe algorithm followed by
a discussion of its performance on the unconstrained global
optimization examples in \cite{rtrkn88}.

\section{Constraint Processing}

Constraint processing \cite{brnng77,ssmnstl80}
is the general approach in computer
problem solving where one expects the problem to be solved by merely
entering the constraints between the variables of the problem without
the need for any algorithm in addition to a general-purpose
constraint solver that comes with the system.
Of course, in interesting problems the
number of constraints is large enough to require an algorithm to
generate them. This is also the case in the global optimization problems
on which we report in this paper.

A common form of constraints are {\em range} constraints: where the
constraint solver updates the sets of possible values of the
variables in view of the relations constraining these values.
In problems where the variables range over the real
numbers, an efficient way to represent such sets of possible values
is as intervals with floating point numbers as bounds. The primitive
constraining relations include the usual numerical relations
$=$, $\leq$, $\geq$, and so on. They also include
the sum relation ($\hbox{sum}(x,y,z)$ iff $x+y = z$), the product relation
($\hbox{prod}(x,y,z)$ iff $x \times y = z$) and the exponentiation
relation ($\hbox{exp}(x,n,y)$ iff $x^n = y$ for positive integer
exponents).

Typical constraint-solving software allows new constraints to be defined
in terms of primitive or previously defined constraints.

\section{The Consistency Method for Interval Constraints}

The consistency method consists in removing
from the intervals of some of the variables certain values
that are not consistent with the constraints.
It was first
used in discrete combinatorial problems \cite{mntn74,wltz75,mckw77}.
Here is an example of how it works with interval constraints.

Suppose
that initial conditions or other constraints determine that the values
for the variables $x$, $y$, and $z$ have to lie in the intervals
$[3,5]$, $[1,3]$, and $[2,5]$. It is then clear that all three intervals
contain values for their variables that are not consistent with,
for instance, the
constraint $x+y=z$. For example, $z$ has to be at least 4.
Interestingly, the
constraint not only causes the interval for $z$ to change, but also
those for $x$ and $y$.
As a result, the constraint causes the above intervals to narrow to
$[3,4]$, $[1,2]$, and $[4,5]$.
Thus the {\em sum} constraint includes the addition operation and its
inverse. Similarly for the product constraint.

Such considerations may also yield the conclusion that an interval is
empty. In that case, when the usual precaution of outward rounding is
observed, a proof of nonexistence of a solution is obtained. In this way,
interval constraints shares with interval analysis the property of being
able to find complete sets of solutions to nonlinear
equations as well as global optima.

The purpose of computation in the consistency method is to use all
constraints to narrow as much as possible the intervals for all of the
variables. The narrower the intervals, the more information they convey
about the solution, if one exists. It is quite common to obtain
intervals that are not much wider than is allowed by the precision of
the floating point system. It is also quite common to obtain intervals
that are so wide as not to contain any useful information about possibly
existing solutions. In the application reported here we obtained
intervals not exceeding the tolerances specified.

We have used the BNR Prolog system \cite{BNR88,bnldr94},
which embeds interval constraints
and a consistency algorithm in the Prolog language \cite{ckp72,ssh86}.

\section{Examples of the Role of Interval Constraints in Global
Minimization}

The Moore-Skelboe algorithm for global
minimization subdivides the
area of interest in successively smaller boxes (i.e. Cartesian products
of intervals)
and computes an interval
for the objective function over each of these.
During the execution of the algorithm a list of disjoint
boxes exists of which the union is guaranteed to contain
the global minima.
At any time during execution we know that global minima have to lie
between the lowest upper bound in the list and the lowest lower bound.
The algorithm terminates as soon as these bounds differ less than
a predetermined tolerance.

We have programmed in BNR Prolog a
close counterpart of the Moore-Skelboe algorithm.
With the algorithm as stated in \cite{rtrkn88} we were unable to avoid
memory overflow in the examples cited. But with some obvious and minor
modifications, but without using constraints,
we obtained results much like the ones reported in
\cite{rtrkn88}; see the table below.

The Moore-Skelboe algorithm only uses information about the objective
function obtained by evaluation of an inclusion function for it over
each of
the boxes in the list. With interval constraints we can obtain more
information about the values of the objective function.
By including such constraints
in the Moore-Skelboe algorithm, we increase its effectiveness.
Before discussing the algorithm,
we show by an example
how interval constraints can give such additional information.

To find the global optimum with interval constraints, the objective
function is translated to a set of primitive constraints. For example,
suppose we wish to find the global minimum of
$$f(x) =
  x_1^2(4 + x_1^2(-2.1 + (1/3)x_1^2))
  + 4x_2^2(x_2^2-1) + x_1x_2.
$$
This is example 3 on page 93 of Ratschek and Rokne \cite{rtrkn88}.
The BNR Prolog system for interval constraints allows one to
enter the definition of the relation \verb$f$ defined as
\begin{verbatim}
f(X1,X2,Y)
:- Y is X1**2*(4+X1**2*(-2.1+X1**2/3))
   +4*X2**2*(X2**2-1)+X1*X2.
\end{verbatim}
Read this as: ``The relation \verb$f$ holds between \verb$X1$, \verb$X2$,
and \verb$Y$ if \verb$Y$ equals the polynomial shown.''
A constraint containing an arithmetic expression is translated to an
equivalent set of primitive constraints, typically introducing auxiliary
variables.

We continue with example 3 from Ratschek and Rokne. They start with the
independent variables in the intervals $[-2.5,2.5]$ and $[-1.5,1.5]$.
After suitable type declarations, we can now enter the
BNR Prolog command
\begin{verbatim}
-2.5 =< X1,X1 =< 2.5, -1.5 =< X2,X2 =< 1.5,
 f(X1,X2,Y)
\end{verbatim}
using the definition for \verb$f$ and setting up constraints for
\verb$X1$ and \verb$X2$ as shown.
By leaving \verb$Y$ initially unconstrained, it gets an interval of
possible values equal to the result of an inclusion function
\cite{rtrkn88} or interval extension \cite{hnsn92}.
In this example it is $[-70,40]$.

Now, BNR Prolog allows one to enter additional constraints. So why
not probe the large interval for \verb$Y$
obtained from the inclusion function
by means of a constraint on \verb$Y$?
In this way we find, for example, that \verb$Y$ cannot be less than $-13$,
because the constraint system
\begin{verbatim}
-2.5 =< X1,X1 =< 2.5, -1.5 =< X2,X2 =< 1.5,
 f(X1,X2,Y), Y < -13
\end{verbatim}
is found not to have any solutions. This
establishes $-13$ as a lower bound for the global minimum of $f$
within the given intervals for $x_1$ and $x_2$. When
we enter instead $-12$, the constraint system leaves open the
possibility of solutions. This points up an important limitation of the
consistency method: all it does is remove inconsistent values. Hence the
{\em absence} of solutions can be reported with certainty; failure to
determine absence of solutions does not imply the {\em presence} of a
solution. But of course the fact that $f$ is everywhere
defined on the given initial intervals for $x_1$ and $x_2$
allows us to obtain an upper bound to the global minimum by arbitrarily
selecting point intervals for $x_1$ and $x_2$
and again invoking the interval extension of $f$.
This is done in the following BNR Prolog command, where we choose the
midpoints of the previous intervals as the arbitrary point intervals for
\verb$X1$ and \verb$X2$:
\begin{verbatim}
X1 == 1.23456789, X2 == 0.987654321,
f(X1,X2,Y).
\end{verbatim}
which results in the interval
$$[3.52202582505467, 3.52202582505469]$$ for
\verb$Y$. The fact that $f$ is defined for these $x$-values
implies that the upper bound of the interval for \verb$Y$ is an
upper bound for the global minimum. Thus the interval constraint system
has improved the interval $[-70,40]$,
obtained from the inclusion function used, to, say, $[-13,4]$.
The Moore-Skelboe algorithm relies on the intervals given by an
inclusion function. Our modification is to use instead an improved
interval by the same method as in this example.

Trying the above command with \verb$Y < -13$ was of course a lucky
guess. The point is, that with any number equal to about $-13$ or less,
a single constraint evaluation can establish it as a lower bound. This
is what we do in our modification of the Moore-Skelboe algorithm.

To further improve the results we have to find a systematic way for
subdividing the initial search space. For that we adopt the method of
the Moore-Skelboe algorithm.

\section{Modification of the Moore-Skelboe Algorithm}

In our version of the Moore-Skelboe algorithm,
we exploit the fact that in BNR Prolog all computation is performed
via interval constraints, so that the extra constraint on the value of
the objective function, as discussed above, can be inserted easily. 
\\
\begin{center}{   
          \bf  The Modified Algorithm\\}
\end {center}


 input: \\ \hspace*{0.2in}
        (1) a box $X = (X_{1},  \cdots , X_{m})$, \\ \hspace*{0.4in}
                  where $X_{i} = [a_{i}, b_{i}] (i = 1, \cdots , m)$. \\ \hspace*{0.2in}
        (2) an inclusion function $F$ of $f$. \\ 
 output:\\ \hspace*{0.2in}
         $L = ( (V_{1}, y_{1}),   \cdots ,  (V_{n}, y_{n}) )$, where\\ \hspace*{0.4in}
             $V_{i}$ is a box, $y_{i}=F(V_{i}), i=1, \cdots , n$, \\ \hspace*{0.4in}
             $lb(y_{1}) \leq lb(y_{2}) \leq \cdots \leq lb(y_{n})$. \\ 

 {\bf 1:} \\ \hspace*{0.2in}
          $y = F(X)$ ; $L = ((X, y))$ ; \\ \hspace*{0.2in}
          $f_{min} = ub(F(c))$ , where $c = mid(X)$.

 {\bf 2:} \\ \hspace*{0.2in}
          remove the 1st element $(V, y)$ of $L = ((V, y), \cdots )$, \\ \hspace*{0.4in}
              where $V = (Z_{1}, \cdots , Z_{m})$.

 {\bf 3:} \\ \hspace*{0.2in}
          find $Z_{k}$ such that \\ \hspace*{0.4in}
              $width(Z_{k}) = max(width(Z_{i})) (i=1, \cdots , m)$; \\ \hspace*{0.2in}
          bisect $Z_{k}$ into $Z_{k1}$, $Z_{k2}$, and obtain subboxes\\ \hspace*{0.4in}
             $V_{1} = (Z_{1}, \cdots , Z_{k-1}, Z_{k1}, Z_{k+1}, \cdots, Z_{m})$ \\ \hspace*{0.4in}
             $V_{2} = (Z_{1}, \cdots , Z_{k-1}, Z_{k2}, Z_{k+1}, \cdots, Z_{m})$ 

 {\bf 4:} \\ \hspace*{0.2in}
         $v_{1} = F(V_{1})$; $v_{2} = F(V_{2})$; \\ \hspace*{0.2in}
         enter $(V_{1}, v_{1}), (V_{2}, v_{2})$ into $L$ such that \\ \hspace*{0.2in}
         the lower bounds of the 2nd members of the \\ \hspace*{0.2in}
         pairs of $L$ do not decrease.

 {\bf 5:} \\ \hspace*{0.2in}
         $f_{min} = min(f_{min}, ub(F(c_{1})), ub(F(c_{2})))$, \\ \hspace*{0.4in}
              where $c_{1}=mid(V_{1}), c_{2}=mid(V_{2})$. \\ \hspace*{0.2in}
         for each $(V, v)$ of $L$,there are two different cases: \\ \hspace*{0.35in}
            (1) discard it if interval constraint $v \leq f_{min}$ \\ \hspace*{0.55in}
                can not be consistently added to constraint system; \\ \hspace*{0.35in}
            (2) keep it in $L$ and add $v \leq f_{min}$ to\\ \hspace*{0.55in}
                constraint system if $v \leq f_{min}$ \\ \hspace*{0.55in}
                can be consistently added. 
                  
 {\bf 6:} \\ \hspace*{0.2in}
          if termination criteria hold then halt \\ \hspace*{0.3in}
                     else goto step $2$. \\

                  

\section{Computational Results}
The above algorithm was
tested on the three examples used by Ratschek and Rokne for testing
the Moore-Skelboe
algorithm. To compare our modification with the original version,
we implemented both in BNR-Prolog. The following
tables show the results.
Here follows an explanation of the symbols used in the table.

{\halign{\indent#\hfil&\quad#\hfil\cr
$D_{1} \times D_{2}$ &Initial box \cr
                       &(Step 1 of the algorithm).\cr
$\epsilon$             &Intended absolute accuracy \cr 
                       &(Termination criterion(3.6)\cite{rtrkn88}).\cr
$N$       &Number of inclusion function \cr 
          &evaluations till termination when\cr
          &Moore-Skelboe algorithm is used. \cr
$N'$      &Number of inclusion function \cr 
          &evaluations till termination when\cr
          &the modified algorithm is used. \cr
$y$*      &Lower bound of the global \cr
          &minimum.\cr
($X$*, $y$*)&Leading pair at termination, \cr
          &$X$*$=X_{1}$* $\times X_{2}$*.\cr}}

\pagebreak
\begin{center}{   
          \bf  Example 1\\}

\begin{tabular}{|l|l|l|} \hline\hline
           & $D_{1}$   & [-100,100]       \\ \cline{2-3}
 Input     & $D_{2}$   & [-100,100]       \\ \cline{2-3}
           & $\epsilon$& 1.0e-6           \\ \hline
           & $N$       & 143              \\ \cline{2-3}
           & $y$*      &                  \\ \cline{2-3}
 Moore-    &           & -9999.0000167    \\ \cline{2-3}
 Skelboe   & $X_{1}$*  & [-100,-99.9999]  \\ \cline{2-3}
           & $X_{2}$*  & [-100,-99.9999]  \\ \hline
           & $N'$      & 57               \\ \cline{2-3}
 Modified  & $y$*      & -9999.000015693  \\ \cline{2-3}
 Moore-    &           & -9999.000016667    \\ \cline{2-3}
 Skelboe   & $X_{1}$*  & [-100,-99.9999]  \\ \cline{2-3}
           & $X_{2}$*  & [-100,-99.9999]  \\ \hline
 Comparison& $N/N'$    & 2.5              \\ \hline\hline
\end{tabular}
\end {center}

\begin{center}{   
          \bf  Example 2\\}

\begin{tabular}{|l|l|l|} \hline\hline
           & $D_{1}$   & [0,1]            \\ \cline{2-3}
 Input     & $D_{2}$   & [0,1]            \\ \cline{2-3}
           & $\epsilon$& 1.0e-6           \\ \hline
           & $N$       & 81               \\ \cline{2-3}
           & $y$*      &                  \\ \cline{2-3}
 Moore-    &           & 1.0              \\ \cline{2-3}
 Skelboe   & $X_{1}$*  & [0,9.5367e-7]    \\ \cline{2-3}
           & $X_{2}$*  & [0,9.5367e-7]    \\ \hline
           & $N'$      & 37               \\ \cline{2-3}
 Modified  & $y$*      & 1.00000055190    \\ \cline{2-3}
 Moore-    &           & 1.0              \\ \cline{2-3}
 Skelboe   & $X_{1}$*  & [0,5.5189e-7]    \\ \cline{2-3}
           & $X_{2}$*  & [0,4.6928e-7]    \\ \hline
 Comparison& $N/N'$    & 2.1              \\ \hline\hline
\end{tabular}
\end {center}

\begin{center}{   
          \bf  Example 3\\}

\begin{tabular}{|l|l|l|} \hline\hline
           & $D_{1}$   & [-2.5,2.5]         \\ \cline{2-3}
 Input     & $D_{2}$   & [-1.5,1.5]         \\ \cline{2-3}
           & $\epsilon$& 1.0e-1             \\ \hline
           & $N$       & 1171               \\ \cline{2-3}
           & $y$*      &                    \\ \cline{2-3}
 Moore-    &           & -1.1313            \\ \cline{2-3}
 Skelboe   & $X_{1}$*  & [0.0391, 0.0781]   \\ \cline{2-3}
           & $X_{2}$*  & [-0.7239,-0.6979]  \\ \hline
           & $N'$      & 223                \\ \cline{2-3}
 Modified  & $y$*      & -1.0316            \\ \cline{2-3}
 Moore-    &           & -1.1258            \\ \cline{2-3}
 Skelboe   & $X_{1}$*  & [-0.1563,-0.1172]  \\ \cline{2-3}
           & $X_{2}$*  & [0.6094, 0.6563]   \\ \hline
 Comparison& $N/N'$    & 5.2                \\ \hline\hline
\end{tabular}
\end {center}

\section{Conclusions}
We combined interval constraints, as implemented in BNR Prolog, with
the older method of interval analysis, as embodied in the Moore-Skelboe
algorithm for global optimization. We ran the resulting program on
several examples. We reported a comparison on the three examples we
found in the literature where the performance of the Moore-Skelboe
algorithm is given. We found a factor of two to five improvement in
the number of function evaluations.

We also attempted to
implement the original
Moore-Skelboe algorithm in BNR Prolog, making only such modifications as
to allow us to get results at all. The performance we found was very
nearly the same as reported by Ratschek and Rokne.

We feel this result suggests investigating other numerical applications
of interval constraints. As interval constraints translate {\em any} form of
computation to constraint solving, {\em constrained} global optimization
is an attractive area for such new applications.

\begin{thebibliography}{99}

\bibitem{bnldr94}
F. Benhamou and W.~J. Older,
``Applying interval arithmetic to real, integer, and {B}oolean
  constraints,
{\em Journal of Logic Programming} 
(to appear).

\bibitem{BNR88}
BNR, {\bf BNR Prolog user guide and reference manual},
1988.

\bibitem{brnng77}
A. Borning,
``Thinglab --- an object-oriented system for building simulations using
  constraints'',
{\em Proc. 5th Int. Joint Conf. on Artificial Intelligence}, 1977.

\bibitem{ckp72}
A.~Colmerauer, H.~Kanoui, R.~Pas\'{e}ro, and P.~Roussel,
{\bf Un Esyst\`{e}me de communication homme-machine en fran\c{c}ais},
Universit\'{e} d'Aix-Marseille II, 
Groupe d'Intelligence Artificielle, Technical Report, 1972.

\bibitem{hnsn92}
E. Hansen,
{\bf Global Optimization Using Interval Analysis},
Marcel Dekker, 1992.

\bibitem{mckw77}
A. K. Mackworth, ``Consistency in networks of relations'',
{\em Artificial Intelligence}, 1977, Vol. 8, pp. 99--118.

\bibitem{mntn74}
Ugo Montanari,
``Networks of constraints: Fundamental properties and applications to
  picture processing'',
{\em Information Science}, 1974, Vol. 7, No. 2, pp. 95--132.

\bibitem{rtrkn88}
H.~Ratschek and J.~Rokne,
{\bf New Computer Methods for Global Optimization}.
Ellis Horwood/John Wiley, 1988.

\bibitem{ssh86}
L. Sterling and E. Shapiro,
{\bf The Art of {P}rolog}, MIT Press, 1994.

\bibitem{ssmnstl80}
G. J. Sussman and G. L. Steele Jr.,
``Constraints --- a language for expressing almost-hierarchical
  descriptions'',
{\em AI Journal}, 1980, Vol. 14, pp. 1--39.

\bibitem{wltz75}
D.~Waltz,
``Understanding line drawings in scenes with shadows'',
In: P. H. Winston (editor), {\bf The Psychology of Computer
  Vision}, McGraw-Hill, 1975, pp. 19--91.

\end{thebibliography}

\end{document}

