\documentclass[a4paper]{article}
\usepackage{amsmath}
\usepackage{url}
\title{Rdonlp2 - an R interface to DONLP2}
\date{\today (Version 0.3-1)}
\author{Ryuichi Tamura(\texttt{ry.tamura @ gmail.com})}
%%\VignetteIndexEntry{Rdonlp2 tutorial}
\begin{document}
\maketitle

\tableofcontents
\section{Abount This Package}

Rdonlp2 is an R package to use Peter Spellucci's DONLP2 from R. DONLP2
is the copyrighted software written by Peter Spellucci for solving
nonlinear programming problems.

DONLP2 is available from:
\begin{itemize}
 \item \url{http://plato.la.asu.edu/donlp2.html}
\end{itemize}
Rdonlp2 is a wrapper for ANSI-C version of DONLP2(also called DONLP3):
\begin{itemize}
 \item \url{http://plato.asu.edu/ftp/donlp2/donlp2_intv_dyn.tar.gz}
\end{itemize}
Current Rdonlp2 package is available from:
\begin{itemize}
 \item \url{http://arumat.net/Rdonlp2/Rdonlp2_0.3-1.zip}(Windows Binary)
 \item \url{http://arumat.net/Rdonlp2/Rdonlp2_0.3-1_R_i386-apple-darwin8.8.1.tar.gz}(OSX Universal Binary)
 \item \url{http://arumat.net/Rdonlp2/Rdonlp2_0.3-1.tar.gz}(Source File)
\end{itemize}
Since Rdonlp2 is simply wrapper program, user is required to refer PDF
manual included in \texttt{donlp2\_intv\_dyn.tar.gz} for the detail
of the algorithm.

As for condition of use, Please refer Section\ref{copyright}.

\section{Problem Definition}

R Package Rdonlp2 solves following nonlinear minimization problem with
linear, nonlinear constraints as well as parameter bounds:
\begin{align*}
 \min_x f(x)& \quad\text{subject to }\quad x\in \mathcal{S}\\
           \mathcal{S} \in& \quad\{ \quad x\in R^n, \\
                          & \qquad x_l\le x\le x_u, \\
                          & \qquad b_l \le Ax \le b_u, \\
                          & \qquad c_l \le c(x) \le c_u \quad\},
\end{align*}
where, $f(x)$ is a continuous function, $x_l,x_u$ are bounds for
parameters($x$), $b_l, b_u$ are bounds for linear combinations
$Ax$(linear constraints), and $c_l,c_u$ are bounds for nonlinear
function $c(x)$(nonlinear constraints). To describe equality constraint
or parameter constancy, let lower and upper bounds for constraint be
equal.

\subsection{R function donlp2()}

Rdonlp provides single function for this problem:
\begin{verbatim}
 donlp2 <- function(par, fn,
                   par.upper=rep(+Inf, length(par)),
                   par.lower=rep(-Inf, length(par)),

                   A = NULL,
                   lin.upper=rep(+Inf, length(par)),
                   lin.lower=rep(-Inf, length(par)),

                   nlin = list(),
                   nlin.upper=rep(+Inf, length(nlin)),
                   nlin.lower=rep(-Inf, length(nlin)),

                   control=donlp2.control(),
                   control.fun=function(lst){return(TRUE)},
                   env=.GlobalEnv, name="Rdonlp2")
\end{verbatim}
where,
\begin{description}
 \item{\texttt{fn}} - the objective function to be minimized. Currently,
	    \texttt{fn} must take only one argument, and the parameter
	    vector(\texttt{par}) will be passed to \texttt{fn} during
	    the optimization. The first element of return value must be
	    the evaluated value.
	    
 \item{\texttt{par}} - parameter vector(vector object)

 \item{\texttt{par.upper}, \texttt{par.lower}} - upper and lower bounds for
	    parameter vector, respectively. Their length must equal to
	    \texttt{length(par)}. If some elements are unbounded,
	    specify \texttt{+Inf} or \texttt{-Inf} explicitly.

 \item{\texttt{A}} - the matrix object that represents linear
	    constraints. Its columns must be equal to
	    \texttt{length(par)}, and its rows must be equal to the
	    number of linear constraints.
 \item{\texttt{lin.upper}, \texttt{lin.lower}} - upper and lower bounds
	    for linear constraints, respectively. Their length must equal
	    to the number of linear constraints. If some elements are unbounded,
	    specify \texttt{+Inf} or \texttt{-Inf} explicitly.

 \item{\texttt{nlin}} - list object whose elements are functions that
	    represents nonlinear constraints.  rule for argument and
	    return value is the same as \texttt{fn}, i.e., these
	    functions take only one arugument(\texttt{par}), and return
	    a vector object whose first element is the evaluated value.

 \item{\texttt{lin.upper}, \texttt{lin.lower}} - upper and lower bounds
	    for nonlinear constraints, respectively. Their length must
	    equal to \texttt{length(nlin)}. If some elements are
	    unbounded, specify \texttt{+Inf} or \texttt{-Inf}
	    explicitly.
 \item{\texttt{control}} - control parameters that defines the behavior
	   of DONLP2. See below for details.
 \item{\texttt{control.fun}} - \texttt{donlp2()} reports a group of
	   optimization parameters in every iteration(see below for
	   details). This (read-only) information can be available
	   within \texttt{control.fun()}, in which user can decide
	   whether the optimization should be iterrupted. Set its return
	   value to \texttt{FALSE} to interrupt \texttt{donlp2()}.
\item{\texttt{env}} - the environment in which objective, constraint,
	   control functions are evaluated.
\item{\texttt{name}} - an character object that specify file
	   name(without extension) of
	   output file. DONLP2 can output following 2 files(\texttt{name}.pro and
	   \texttt{name}.mes) in working directory which
           contain detailed information during the optimization.
	      \begin{itemize}
	       \item \texttt{name}.pro: the results of optimization and other
		     information during the optimazation will be written. the
		     latter depends on the values of
		     \texttt{te0},\texttt{te1},\texttt{te2},\texttt{te3}.
	       \item \texttt{name}.mes: other message from DONLP2 (warnings for
		     ill-conditions, etc) will be written. This will be helpful
		     if your R code does not work correctly.
	      \end{itemize}
\end{description}

\section{Details}

\subsection{Initial Values}

Although intial values(given to \texttt{par}) should be carefully
determined by user, DONLP2 has the feature to correct initial values
\textbf{automatically} that violate given constraints.

\subsection{Objective Function and its Gradients}\label{fngr}

Objective function should be implemented so that parameter vector is
the only argument and the first element of return value is
numeric. Typically, it looks like:
\begin{verbatim}
objective.fun <- function(par){
  # calculation with par, and
  # results are stored to ans
  :
  ans # return value
}

:

ret <- donlp2(par=par, fn=objective.fun, ....)
\end{verbatim}

User can implement gradient function to improve accuracy and
efficiency. Gradient function should be implement such that parameter
vector is the only argument and return value is vector of length equal
to \texttt{length(par)}. Then, assign it to the attibute \texttt{gr} of
objective function:
\begin{verbatim}
# par is vector of length n
grad.fun <- function(par){
  
  c(v1, v2, ..., vn)
}
# assign grad.fun to objective.fun's gr attribute
attr(objective.fn, "gr") <- grad.fun
\end{verbatim}


\subsection{Bounds and Equality Constraints}

Bounds for parameter, linear and nonlinear constraints are given as
vector of appropriate
length(\texttt{par.upper},\texttt{par.lower},\texttt{lin.upper},\texttt{lin.lower},\texttt{nlin.upper},\texttt{nlin.lower}). If
some parameter or constraints are bounded from below (above), then
specify \texttt{+Inf}(\texttt{-Inf}), respectively.

Set upper and lower bounds to be equal if equality constraint is
needed.
\begin{verbatim}
# par[1]<0, 0<par[2]<1, par[3]>1
par.lower <- c(-Inf, 0,    1)
par.upper <- c(   0, 1, +Inf)

# two linear constraints on two parameters
# (1) par[1]+par[2]=0
# (2) par[1]-2*par[2]+10>0
lin.lower <- c(0,  -10)
lin.upper <- c(0, +Inf)
\end{verbatim}


\subsection{Linear Constraints}

Bounds arguments and the coefficient matrix \texttt{A} represents the
linear constraints. Every row of \texttt{A} stands for linear
combination of parameters:
\begin{verbatim}
# two linear constraints on two parameters
# (1) par[1]+par[2]=0
# (2) par[1]-2*par[2]+10>0
lin.lower <- c(0,  -10)
lin.upper <- c(0, +Inf)

A = rbind( c(1, 1),   # 1st linear compination
           c(1,-2) )  # 2nd linear combination
\end{verbatim}

\subsection{Nonlinear Constraints and its Gradients}

Nonlinear constraints are represented as their bounds given to
\texttt{nlin.upper} and \texttt{nlin.lower}, and user-defined
function(and gradients). The way to implement function and gradients are
the same as objective function and its gradients(see Section
\ref{fngr}):
\begin{verbatim}
# Nonlinear constraints: 1 constraint on 2 parameters
# par[1]*par[2] = 1
nlcon1 <- function(par){
  par[1]*par[2]
}
nlcon1.gr <- function(par){
  c(par[2], par[1])
}
attr(nlcon1, "gr")<-nlcon1.gr

nlin.upper = c(1)
nlin.lower = c(1)
:
ret <- donlp2(par, fn, 
              nlin=list(nlcon1), 
              nlin.upper=nlin.upper,nlin.lower=nlin.lower,.....)
\end{verbatim}
All the nonlinear constraint function are collected in a
list(\texttt{nlin}).

\subsection{Numerical Gradients}\label{numdif}

User can omit the implementation of gradients. In this case one of 3
algorithm fro numerical differentiation provided by DONLP2 will be
performed. If there are \texttt{n} parameters,
\begin{enumerate}
 \item the forward difference: requires \texttt{n} additional
       evaluations of function(\texttt{difftype=1}).
 \item the central difference: requires \texttt{2n} additional
       evaluations of function(\texttt{difftype=2}).
 \item a sixth order approximation computing a Richardson extrapolation
       of the three symmetric difference: requires \texttt{6n}
       additional evaluations of function(\texttt{difftype=3}).
\end{enumerate}
By default, Rdonlp2 uses 3rd algorithm(most acculate, but quite
costly). User can change this by the control variable(below)
\texttt{difftype}.

Currently, if user want to use analytical gradients, he \textbf{must}
implement \textbf{all} of the gradients for \textbf{both} objective
function and nonlinear constraint functions. If one of them are not
implemented, Rdonlp2 gives up using any gradient function and uses
numerical method instead.

\subsection{Control Variables}

User can control the behavior of \texttt{donlp2()} by
\texttt{donlp2.control()}.  \texttt{donlp2.control()}returns a group of
default control parameters as list object, so user change some of them
by giving \texttt{tag=value} pairs as arguments. Currently following
tags(control variables) are available (values in () are the defaults).
\begin{itemize}
 \item Settings
       \begin{itemize}
	\item \texttt{iterma} (\texttt{4000}) - maximum number of iterations
	\item \texttt{nstep} (\texttt{20}) - maximum number of tries in the
	      backtracking allowed.
	\item \texttt{fnscale}(\texttt{1}) - set \texttt{-1} for
	      maximization instead of minimization. values other than
	      \texttt{1,-1} are not recommended.
       \end{itemize}
 \item Tunings and erfomance of the optmization
       \begin{itemize}
	\item \texttt{tau0} (\texttt{1.0}) - the positive amount how much any
	      constaint other than abound can be violated. A small
	      \texttt{tau0} diminishes the efficiency of DONLP2, while a large
	      \texttt{tau0} may degarde the reliability of the code.
	\item \texttt{tau} (\texttt{0.1}) - gives a weight between
	      descent for \texttt{fn} and infeasibility and is also used
	      as a safety parameter for chosing the penalty weigths. It
	      can be chosen larger zero at will, but useful values are
	      between 0.1 and 1. The smaller tau, the more may
	      \texttt{fn} be scaled down. Tau is also used as an
	      additive increase for the penalty-weights. Therefore it
	      should not be chosen too large, since that degrades the
	      performance.
	\item \texttt{del0} (\texttt{1.0}) - The positive amount by which
	      constraints are considered binding. If too small, the
	      indentification of correct sets of binding constraints may be
	      delayed. If too large, DONLP2 will escape to the full
	      regularlized SQP method(quite costly).  Good values are in [0.01,
	      1.0]
       \end{itemize} 
 \item Termination criteria
       \begin{itemize}
	\item \texttt{epsx}(\texttt{1e-5}) - successful temination is
	      indicated if the Kuhn-Tucker criteria are satisfied within
	      the value.
	\item \texttt{delmin}(\texttt{0.1*del0}) - constraints are
	      considered as sufficiently satisfied if absolute values of
	      their violation are less than the value.
	\item \texttt{epsdif}(\texttt{1e-8}) - relative precision in the
	      gradients.
	\item \texttt{nreset.multiplier} (\texttt{1}) - maximum number of steps
	      (\texttt{nreset.multiplier} times \texttt{n}) until a ``restart''
	      of the accumulated quasi-newton-update is tried. Value
	      should be integer between 1 and 4.
       \end{itemize}
 \item Numerical differentiation
       \begin{itemize}
	\item \texttt{difftype} (\texttt{3}) - See Section\ref{numdif}.
	\item \texttt{taubnd} (\texttt{0.1}) - The positive amount by which
	      bounds may be violated if numerical differention is used.
	\item \texttt{epsfcn} (\texttt{1e-16}) - relative precision of the
	      function evaluation routine.
	\item \texttt{hessian} (\texttt{FALSE}) - if \texttt{TRUE}, caliculate numeric Hessian matrix evaluated at the optimum by numerical differentiation specified in \texttt{difftype}.
       \end{itemize}
 \item Ouputs on console or files
       \begin{itemize}
	\item \texttt{report}(\texttt{FALSE}) - if \texttt{TRUE}, a list
	      object which contains detailed information will be passed
	      to \texttt{control.fun()}. See Section\ref{taglist}.
	\item \texttt{rep.freq}(\texttt{1}) - the frequency of
	      report. the report will be passed to \texttt{control.fun}
	      every \texttt{rep.freq} iterations.
	\item \texttt{te0} (\texttt{TRUE}) - if \texttt{TRUE} one-line-output
	      for every step is printed on R console.
	\item \texttt{te1} (\texttt{FALSE}) - if \texttt{TRUE}
	      post-mortem-dump of accumlated information is printed on R
	      console. Note in Rdonlp2 the same information will be
	      passed to \texttt{control.fun()}(See Section\ref{accinf}).
	\item \texttt{te2} (\texttt{FALSE}) - if \texttt{TRUE}, more detailed
	      information is ``pretty-printed'' on R console.
	\item \texttt{te3} (\texttt{FALSE}) - if \texttt{TRUE}, and output file is specified in \texttt{donlp2}, also print the gradients and Newton-Raphson update in upper trianglar matrix in the \texttt{.pro} file.
	\item \texttt{silent} (\texttt{FALSE}) - if \texttt{TRUE}, \texttt{donlp2()} runs quietly, i.e. nothing is output to R console and \texttt{.pro, .mes} files are never created even if you specify file name in the \texttt{name} argument of \texttt{donlp2().}
	\item \texttt{intakt} (\texttt{TRUE}) - if \texttt{TRUE},
          various information(depends on \texttt{te0,te1,te2}) from current iteration step is output to R console.
       \end{itemize}
\end{itemize}

Some of parameteters listed above are not well documented in this
tutorial. Please refer the original PDF manual included in DONLP
distribution for details.

\section{Information during the Optimization}\label{accinf}

User may want to know what happens during the optimation, and if some
parameters are 'undesirable' he may stop the execution. Rdonlp2
provides the way to access the information via
\texttt{control.fun(lst)}. On each iterantion, 35 parameters that
tell us how the optimization is working are passed to
\texttt{control.fun(lst)}. The last four parameters are not reported
until the optimization has finished. Parameters are collected into single
list object, so user can easily access some of them by specifying the
\texttt{tag}s. Complete tag list is given Section \ref{taglist}.

\texttt{control.fun()} should return \texttt{TRUE} if user want to
continue and \texttt{FALSE} if user want to interrupt.

\begin{verbatim}
## Example 1
## keep track of current lagrange multiplier values
mycontrol <- function(lst){
 print(lst$u)
 TRUE # return TRUE to continue execution
}

# tell donlp2 to use mycontrol()
donlp2(.....,control.fun=mycontrol,....)

## Example 2(useless example)
## force to terminate optimization after 10 iterations
mycontrol2 <- function(lst){
  lst$step.nr <= 10 # return FALSE when step.nr>11
}
# tell donlp2 to use mycontrol2()
donlp2(.....,control.fun=mycontrol2,....)
\end{verbatim}

\subsection{Tag list}\label{taglist}

Some of tags listed here are not well documented in this
tutorial. Please refer the original PDF manual included in DONLP
distribution for details.

\begin{itemize}
 \item \texttt{par} - current value of parameter vector.
 \item \texttt{u} - current value of langrange multipliers.
 \item \texttt{w} - current value of penalty terms.
 \item \texttt{gradf} - current gradient vector.
 \item \texttt{step.nr} - step number(total number of iterations when finished).
 \item \texttt{fx} - current value of \texttt{fn}.
 \item \texttt{scf} - scaling of \texttt{fn}.
 \item \texttt{psi} - the weighted penalty term.
 \item \texttt{upsi} - the unweighted penalty term(L1 norm of constraint
       vector).
 \item \texttt{del.k.1} - bound for currently active constraints.
 \item \texttt{b2n0} - weighted L2 norm of projected gradients.
 \item \texttt{b2n} - L2 norm of gradients based on \texttt{del.k.1}
 \item \texttt{nr} - number of binding constraints.
 \item \texttt{sing} - value other than \texttt{-1} indicates working
       set is singular
 \item \texttt{umin} - infinity norm of negative part of inequalities
       multipliers.
 \item \texttt{not.used} - always \texttt{0}(curretnly not used). 
 \item \texttt{cond.r} - condition number of diagonal part of qr
       decomposition of normalized gradients of binding constraints.
 \item \texttt{cond.h} - condition number of diagonal of cholesky factor
       of updated full Hessian.
 \item \texttt{scf0} - the relative damping of tangential component if
       \texttt{upsi > tau0/2}.
 \item \texttt{xnorm} - L2 norm of \texttt{par}.
 \item \texttt{dnorm} - unscaled L2 norm of \texttt{d}, correction from
       eqp/qp subproblem.
 \item \texttt{phase} - 
       \begin{itemize}
	\item \texttt{-1}: infeasibility improvement phase.
	\item \texttt{0}: initial optimization.
	\item \texttt{1}: binding constraints unchanged
	\item \texttt{2}: \texttt{d} small, maratos correction in use
       \end{itemize}
 \item \texttt{c.k} - number of decreases of penalty weights
 \item \texttt{wmax} - infinity norm of weights
 \item \texttt{sig.k} - stepsize from unidimensional
       minimization(backtracking).
 \item \texttt{cfincr} - number of objective function evaluations for
       stepsize algorithm.
 \item \texttt{dirder} - scaled directional derivative of penalty
       function along \texttt{d}.
 \item \texttt{dscal} - scaling factor for \texttt{d}.
 \item \texttt{cosphi} - cosine of arc between \texttt{d} and previous
       \texttt{d}. 
 \item \texttt{violis} - number of constraints not binding at
       \texttt{par} but hit during line search.
 \item \texttt{hesstype} - type of update for Hessian:
       \begin{itemize}
	\item \texttt{1}: normal P\&M-BFGS update
	\item \texttt{0}: update suppressed
	\item \texttt{-1}: restart with scaled unit matrix
	\item \texttt{2}: standard BFGS
	\item \texttt{3}: BFGS modified by Powell's Device
       \end{itemize}
 \item \texttt{modbfgs} - modification factor for damping the projector
       int the BFGS or pantoja-mayne update.
 \item \texttt{modnr} -modification factor for damping the
       quasi-newton-relation in BFGS.
 \item \texttt{qpterm} -
       \begin{itemize}
	\item \texttt{0}: if \texttt{sing=-1}: temination indicator of
	      the QP solver
	\item \texttt{1}: successful
	\item \texttt{-1}: \texttt{tau} becomes larger than
	      \texttt{tauqp} without slack variables becoming
	      sufficiently small.
	\item \texttt{-2}: infeasible QP problem(theoretically impossible).
       \end{itemize}
 \item \texttt{tauqp} - weight of slack variables in QP solver
 \item \texttt{infeas} - L1 norm of slack variables in QP solver
\end{itemize}

\section{Value from donlp2()}

The return value from \texttt{donlp2()} is the list object with
38(35+3; if \texttt{hessian=FALSE}(default)) or 39(35+4; if
\texttt{hessian=TRUE}) elements with specified tags. The 35
pareameters are identical to those listed in Section \ref{taglist},
and the rest parameters are:
\begin{itemize}
 \item \texttt{nr.update}: the approximated newton-raphson updates in upper triangular matrix.
 \item \texttt{hessian}(if \texttt{hessian=TRUE} in \texttt{donlp2.control()}): numeric Hessian matrix evaluated at the final value \texttt{par}.
 \item \texttt{runtime}: the elapsed time for the optimization.
 \item \texttt{message}: the termination message. See \ref{termreason}.
\end{itemize}

\subsection{The termination Reason}\label{termreason}

When the optimization finishes, DONLP2 returns one of 19 messages listed
below. They are classified to following 3 groups, last of which user
need to decide the result is 'reasonable' and accestable.

\begin{itemize}
 \item 1.-10. : irregular case 
 \item 11.-13.: successful
 \item 14.-19.: successful, but the precision may be very poor.
\end{itemize}

\begin{enumerate}
 \item \texttt{"constraint evaluation returns error with current point"}
 \item \texttt{"objective evaluation returns error with current point"}
 \item \texttt{"qpsolver: working set singular in dual extended qp "}
 \item \texttt{"qpsolver: extended qp-problem seemingly infeasible "}
 \item \texttt{"qpsolver: no descent for infeas from qp for tau=tau\_max"}
 \item \texttt{"qpsolver: on exit correction small, infeasible point"}
 \item \texttt{"stepsizeselection: computed d from qp not a dir. of descent"}
 \item \texttt{"more than maxit iteration steps"}
 \item \texttt{"stepsizeselection: no acceptable stepsize in [sigsm,sigla]"}
 \item \texttt{"small correction from qp, infeasible point"}
 \item \texttt{"kt-conditions satisfied, no further correction computed"}
 \item \texttt{"computed correction small, regular case "}
 \item \texttt{"stepsizeselection: x almost feasible, dir. deriv. very small"}
 \item \texttt{"kt-conditions (relaxed) satisfied, singular point"}
 \item \texttt{"very slow primal progress, singular or illconditoned problem"}
 \item \texttt{"more than nreset small corrections in x "}
 \item \texttt{"correction from qp very small, almost feasible, singular "}
 \item \texttt{"numsm small differences in penalty function,terminate"}
 \item \texttt{"user required termination"}
\end{enumerate}

Some of messages listed above are not well documented in this
tutorial. Please refer the original PDF manual included in DONLP
distribution for details.


\section{Examples}

Example C source files are included in the original DONLP2
distribution. In this section, we rewrite some of them in R and show you
how to code constrained optimization problem with Rdonlp2.

\subsection{examples/simple.c: linear constraint}

\begin{equation}
 \min_{x,y}x^2+y^2\quad\text{subject to}\quad 0< x,y< 100,\quad x+y=1
\end{equation}
with intial value: $(x,y)=(-10,10)$. Note that initial value is
infeasible because $x=-10\notin (0,100)$.

This problem has 2 parameters, 2 parameter bounds, 1 linear equality
constraint. R script looks like:
\begin{verbatim}
p <- c(-10,10)
par.l <- c(0,0); par.u <- c(100,100)

lin.u <- 1; lin.l <- 1
A <- t(c(1,1))

fn <- function(x){
  x[1]^2+x[2]^2
}
ret <- donlp2(p, fn, par.lower=par.l, par.upper=par.u,
              A=A, lin.u=lin.u, lin.l=lin.l, name="simple")
	
\end{verbatim}
Note that \texttt{A} must be represented as row vector(\texttt{1x2})
with single linear constraint. Also we use numerical gradients.

Since control variables are all default values(specifically
\texttt{te0=TRUE} and \texttt{silent=FALSE}), running the script outputs
detailed information on console:

\begin{verbatim}
    1 fx=   0.000000e+00 upsi=  5.9e+01 b2n= -1.0e+00 umi=  0.0e+00 nr   1 si-1
    2 fx=   0.000000e+00 upsi=  2.0e+01 b2n= -1.0e+00 umi=  0.0e+00 nr   2 si-1
    3 fx=   1.000000e+00 upsi=  0.0e+00 b2n=  4.4e-16 umi=  0.0e+00 nr   2 si-1
simple

     n=         2    nlin=         1    nonlin=         0

  epsx= 1.000e-05 sigsm= 3.293e-10

startvalue
  5.0000000e+01   1.0000000e+01 

  eps=  1.08e-19  tol=  0.00e+00 del0=  1.00e+00 delm=  1.00e-06 tau0=  1.00e+00
  tau=  1.00e-01   sd=  1.00e-01   sw=  2.27e-13  rho=  1.00e-06 rho1=  1.00e-10
 scfm=  1.00e+04  c1d=  1.00e-02 epdi=  1.00e-16
  nre=         4 anal=         0
taubnd=  1.00e+00 epsfcn=  1.00e-16 difftype=3

 termination reason:
 KT-conditions satisfied, no further correction computed
 evaluations of f                            3
 evaluations of grad f                       2
 evaluations of constraints                  6
 evaluations of grads of constraints         0
 final scaling of objective           1.000000e+00
 norm of grad(f)                      1.414214e+00
 lagrangian violation                 4.718134e-14
 feasibility violation                0.000000e+00
 dual feasibility violation           0.000000e+00
 optimizer runtime sec's              0.000000e+00


 optimal value of f =   5.00000000000000e-01

 optimal solution  x =
  4.99999999999991e-01  5.00000000000009e-01

  multipliers are relativ to scf=1
  nr.    constraint     multiplier norm(grad) or 1 
    1   5.0000000e-01    0.0000000e+00 
    2   9.9500000e+01    0.0000000e+00 
    3   5.0000000e-01    0.0000000e+00 
    4   9.9500000e+01    0.0000000e+00 
    5   0.0000000e+00    1.0000000e+00   1.3906925e-309
    6   0.0000000e+00    0.0000000e+00 

 evaluations of restrictions and their gradients
 (     6,     0)
 last estimate of cond.nr. of active gradients   1.414e+00

 last estimate of cond.nr. of approx.  hessian   1.000e+00
iterative steps total               3
# of restarts                       0
# of full regular updates           1
# of updates                        1
# of regularized full SQP-steps     0
\end{verbatim}

If user want to access parameter values, simply
\begin{verbatim}
> ret$par
[1] 0.5 0.5
\end{verbatim}

\subsection{examples/simple2.c: nonlinear constraint}

Second example shows optimization problem with parameter bounds and
nonlinear constraint.
\begin{equation}
 \min_{x,y}x^2+y^2\quad\text{subject to}\quad -100< x,y< 100,\quad xy=2
\end{equation}
with intial value: $(x,y)=(10,10)$. We give gradient functions for
objective and nonlinear constraint(\texttt{dfn()} and
\texttt{dnlcon()}, respectively).
\begin{verbatim}
p <- c(10,10)
par.l <- c(-100,-100); par.u <- c(100,100)
nlin.l <- nlin.u <- 2

fn <- function(x){
  x[1]^2+x[2]^2
}
dfn <- function(x){
  c(2*x[1], 2*x[2])
}
attr(fn, "gr") <- dfn

nlcon <- function(x){
  x[1]*x[2]
}
dnlcon <- function(x){
  c(x[2], x[1])
}
attr(nlcon, "gr") <- dnlcon

ret<-donlp2(p, fn, par.u=par.u, par.l=par.l,
            nlin=list(nlcon), nlin.u=nlin.u, nlin.l=nlin.l)
\end{verbatim}

\subsection{example/hs211.c}

This problem uses all types of constraints:
\begin{align*}
 \min_{x_i, i=1\ldots 10} &5.04x_1+0.035x_2+10x_3+3.36x_5-0.063x_4x_7\\
\text{subject to:}&\\
h_1(x) &= 1.22x_4 - x_1 - x_5 = 0\\
h_2(x) &= 98000x_3/(x_4x_9+1000x_3)-x_6=0\\
h_3(x) &= (x_2+x_5)/x_1-x_8 = 0\\
g_1(x) &= 35.82-0.222x_{10}-bx_9\ge 0, b=0.9\\
g_2(x) &= -133+3x_7 -ax_10 \ge 0, a=0.99\\
g_3(x) &= -g_1(x) + x_9(1/b-b)\ge 0\\
g_4(x) &= -g_2(x) +(1/a-a)x_{10} \ge 0 \\
g_5(x) &= 1.12x_1+0.13167x_1x_8 - 0.00667x_1x_8^2-ax_4 \ge 0\\
g_6(x) &= 57.425 + 1.098x_8 - 0.038x_8^2 + 0.325x_6 - ax_7 \ge 0\\
g_7(x) &= -g_5(x) + (1/a-a)x_4 \ge 0\\
g_8(x) &= -g_6(x) + (1/a-a)x_7 \ge 0\\
0.00001 &\le x_1 \le 2000\\
0.00001 &\le x_2 \le 16000\\
0.00001 &\le x_3 \le 120\\
0.00001 &\le x_4 \le 5000\\ 
0.00001 &\le x_5 \le 2000\\
85&\le x_6 \le 93\\
90&\le x_7 \le 95\\
3&\le x_8 \le 12\\
1.2&\le x_9 \le 4\\
145&\le x_{10} \le 162
\end{align*}
We have 10 parameter bounds, 5 linear constraints(arranging terms):
\begin{align*}
h_1 &\rightarrow 1.22x_4 - x_1 - x_5 = 0\\
g_1 &\rightarrow -0.222x_10-bx_9\ge -35.82, b=0.9\\
g_2 &\rightarrow 3x_7 -ax_10 \ge 133, a=0.99\\
g_3 &\rightarrow x_9(1/b-b+b)+0.222x_{10}\ge 35.82\\
g_4 &\rightarrow -3x_7+(1/a-a+a)x_{10} \ge -133,
\end{align*}
1 of which($h_1$) is equality constraint, and 6 nonlinear
constraints:
\begin{align*}
h_2 &\rightarrow 98000x_3/(x_4x_9+1000x_3)-x_6=0\\
h_3 &\rightarrow (x_2+x_5)/x_1-x_8 = 0\\
g_5 &\rightarrow 1.12x_1+0.13167x_1x_8 - 0.00667x_1x_8^2-ax_4 \ge 0\\
g_6 &\rightarrow 57.425 + 1.098x_8 - 0.038x_8^2 + 0.325x_6 - ax_7 \ge 0\\
g_7 &\rightarrow -g_5(x) + (1/a-a)x_4 \ge 0\\
g_8 &\rightarrow -g_6(x) + (1/a-a)x_7 \ge 0,
\end{align*}
2 of which($h_2$, $h_3$) are equality constraint.

R code will be follows:

\begin{verbatim}
a <- 0.99; b <- 0.9
##
## Objective Function
##
fn <- function(par){
  x1 <- par[1]; x2 <- par[2]; x3 <- par[3]; x4 <- par[4]; x5 <- par[5]
  x6 <- par[6]; x7 <- par[7]; x8 <- par[8]; x9 <- par[9]; x10 <- par[10]

  5.04*x1 + 0.035*x2 + 10*x3 +3.36*x5 - 0.063*x4*x7
}
##
## Parameter Bounds
##
par.l <- c(rep(1e-5, 5), 85, 90, 3, 1.2, 145)
par.u <- c(2000, 16000, 120, 5000, 2000, 93, 95, 12, 4, 162)

##
## Constraints
##
linbd  <- matrix(0, nr=5, nc=2)
nlinbd <- matrix(0, nr=6, nc=2)

## linear equality 
linbd[1,] <- c(0,0)          # h1
linbd[2,] <- c(-35.82, Inf)  # g1
linbd[3,] <- c(133, Inf)     # g2
linbd[4,] <- c(35.82,Inf)    # g3
linbd[5,] <- c(-133, Inf)    # g4

## nonlinear equality
h2 <- function(par){
  x1 <- par[1]; x2 <- par[2]; x3 <- par[3]; x4 <- par[4]; x5 <- par[5]
  x6 <- par[6]; x7 <- par[7]; x8 <- par[8]; x9 <- par[9]; x10 <- par[10]

  98000*x3/(x4*x9+1000*x3)-x6
}
nlinbd[1,] <- c(0,0)

## nonlinear equality
h3 <- function(par){
  x1 <- par[1]; x2 <- par[2]; x3 <- par[3]; x4 <- par[4]; x5 <- par[5]
  x6 <- par[6]; x7 <- par[7]; x8 <- par[8]; x9 <- par[9]; x10 <- par[10]

  (x2+x5)/x1 - x8
}
nlinbd[2,] <- c(0,0)

## nonlinear inequality
g5 <- function(par){
  x1 <- par[1]; x2 <- par[2]; x3 <- par[3]; x4 <- par[4]; x5 <- par[5]
  x6 <- par[6]; x7 <- par[7]; x8 <- par[8]; x9 <- par[9]; x10 <- par[10]

  1.12*x1 + 0.13167*x1*x8 - 0.00667*x1*x8^2 - a*x4
}
nlinbd[3,] <- c(0,Inf)

## nonlinear inequality
g6 <- function(par){
  x1 <- par[1]; x2 <- par[2]; x3 <- par[3]; x4 <- par[4]; x5 <- par[5]
  x6 <- par[6]; x7 <- par[7]; x8 <- par[8]; x9 <- par[9]; x10 <- par[10]

  57.425 + 1.098*x8 - 0.038*x8^2 + 0.325*x6 - a*x7
}
nlinbd[4,] <- c(0,Inf)

## nonlinear inequality
g7 <- function(par){
  x1 <- par[1]; x2 <- par[2]; x3 <- par[3]; x4 <- par[4]; x5 <- par[5]
  x6 <- par[6]; x7 <- par[7]; x8 <- par[8]; x9 <- par[9]; x10 <- par[10]

  -g5(par) + (1/a-a)*x4
}
nlinbd[5,] <- c(0,Inf)

## nonlinear inequality
g8 <- function(par){
  x1 <- par[1]; x2 <- par[2]; x3 <- par[3]; x4 <- par[4]; x5 <- par[5]
  x6 <- par[6]; x7 <- par[7]; x8 <- par[8]; x9 <- par[9]; x10 <- par[10]

  -g6(par) + (1/a-a)*x7
}
nlinbd[6,] <- c(0,Inf)

## A is 5(linear constraints) x 10(params) matrix
A <- rbind(c(-1, 0, 0, 1.22,-1, 0,  0, 0,        0,         0), #h1
           c( 0, 0, 0,    0, 0, 0,  0, 0,       -b,    -0.222), #g1
           c( 0, 0, 0,    0, 0, 0,  3, 0,        0,        -a), #g2
           c( 0, 0, 0,    0, 0, 0,  0, 0,(1/b-b)+b,     0.222), #g3
           c( 0, 0, 0,    0, 0, 0, -3, 0,        0,(1/a-a)+a))  #g4

## initial values
p0 <- c(1745, 12e3, 11e1, 3048, 1974, 89.2, 92.8, 8, 3.6, 145)

## control variables
cntl <- donlp2.control(del0=0.2, tau0=1.0, tau=0.1, taubnd=5e-6)

## start constrained optimization
ret <- donlp2(par=p0, fn=fn,
              par.u=par.u, par.l=par.l,
              A=A,
              lin.u=linbd[,2], lin.l=linbd[,1],
              nlin=list(h2,h3,g5,g6,g7,g8),
              nlin.upper=nlinbd[,2], nlin.lower=nlinbd[,1], name="alkylation",
              control=cntl)
\end{verbatim}
Since gradient functions are not implemented in the program and the
numerical differentiation algorithm is \texttt{difftype=3}(default), it
takes about 0.90 sec to finish the computation on my machine(Intel
CoreDuo 2G, Memory 2G), while original C version does within 0.1 sec.

\section{Bugs}

DONLP2 provides much more 'fine-tuning' parameters than those exported
to Rdonlp2. So if you want other parameters that should be exported,
please e-mail me. Also, any comments of suggestions are highly welcome.

\section{Copyright}\label{copyright}

\paragraph*{Original DONLP2:}

\begin{verbatim}
/* Conditions of use:                                                        */
/* 1. donlp2 is under the exclusive copyright of P. Spellucci                */
/*    (e-mail:spellucci@mathematik.tu-darmstadt.de)                          */
/*    "donlp2" is a reserved name                                            */
/* 2. donlp2 and its constituent parts come with no warranty, whether ex-    */
/*    pressed or implied, that it is free of errors or suitable for any      */
/*    specific purpose.                                                      */
/*    It must not be used to solve any problem, whose incorrect solution     */
/*    could result in injury to a person , institution or property.          */
/*    It is at the users own risk to use donlp2 or parts of it and the       */
/*    author disclaims all liability for such use.                           */
/* 3. donlp2 is distributed "as is". In particular, no maintenance, support  */
/*    or trouble-shooting or subsequent upgrade is implied.                  */
/* 4. The use of donlp2 must be acknowledged, in any publication which       */ 
/*    contains                                                               */
/*    results obtained with it or parts of it. Citation of the authors name  */
/*    and netlib-source is suitable.                                         */
/* 5. The free use of donlp2 and parts of it is restricted for research      */
/*    purposes                                                               */
/*    commercial uses require permission and licensing from P. Spellucci.    */
\end{verbatim}

\paragraph*{Rdonlp2}

Copyright (C) 2007 Ryuichi Tamura(ry.tamura at gmail.com). You may
re-distribute or modify this library under GNU LGPL ver.2.

\end{document}