3 USAGE

DONLP2 is written as a selfcontained system of subprograms. The user issues a

call donlp2

supplying necessary information in a series of userwritten subprograms and in some common areas described below. A rudimentary main-program donlp2main.f is contained in the distribution of the code.

3.1 PROBLEM DESCRIPTION

The user has to define the problem through function evaluation routines. This may be done in two ways: Method one gives any function and gradient evaluation code individually. In this case the parameter BLOC has to be set to FALSE. This is the default.

  1. EF(X,FX) returns FX=f(x) given x as input. Arguments DOUBLE PRECISION X(*), FX

  2. EGRADF(X,GRADF) returns GRADF= Ñ f(x) given x as input. Arguments DOUBLE PRECISION X(*),GRADF(*).

  3. EH(I,X,HXI) returns the value of the i-th equality constraint as hxi given i and x as input. Arguments INTEGER I; DOUBLE PRECISION X(*),HXI.

  4. EGRADH(I,X,GRADHI) returns Ñhi(x) as GRADHI given i and x as input. Arguments INTEGER I ; DOUBLE PRECISION X(*),GRADHI(*).

  5. EG(I,X,GXI) returns gi(x) as GXI given i and x as input. Arguments INTEGER I; DOUBLE PRECISION X(*),GXI. Bound constraints have to be programmed here too. The evaluation of their gradients however may be much simplified. See below concerning the usage of GUNIT.

  6. EGRADG(I,X,GRADGI) returns Ñgi(x) as GRADGI given i and x as input. Arguments INTEGER I ;DOUBLE PRECISION X(*),GRADGI(*). For bound constraints, no evaluations must be coded here if GUNIT is initialized properly.

  7. SETUP0 has to set some problem specific data in DONLP2's common areas as described below. SETUP0 is the first subprogram called by DONLP2 prior to any other computation.

  8. SETUP may set user specific data (not transparent to DONLP2) as well as parameters of DONLP2 e.g.in order to supersede DONLP2's standard parameter setting. 
    !!! SETUP is called after SETUP0 and after DONLP2's standard initialization, which is done in O8ST.

  9. SOLCHK may add additional final computations, graphical output using e.g. information from ACCINF (see below) and so on. This distribution contains a solchk.f with an empty body.

donlp2usrfc.f in this directory contains a rudimentary version of these routines which easily may be completed by the user. The bodies of SETUP, EG, EGRADG, EH, EGRADG and SOLCHK may be empty of course. In that case unconstrained minimization using the BFGS-method is done. The correctness of EGRADF, EGRADH and EGRADG can to some extent be checked by running testgrad.f (in this directory) with the set of user routines .

Alternatively the user might prefer or may be forced to evaluate the problem functions by a black-box external routine. In this case he has to set

BLOC=.TRUE.

within SETUP0 and to use a subroutine

EVAL EXTERN(MODE)

This routine has to take XTR, which is a copy of X, from the O8FUEXTR common area as input and to compute , according to the setting of MODE and ANALYT either FU(0:NH+NG) or both FU(0:NH+NG) and FUGRAD(I,0:NH+NG), I=1,N and then return to the user. The meaning of MODE is as follows:

  1. MODE=0 compute only those constraints which depend linearly on exactly one variable, i.e. fixed variables conditions and bound constraints. These constraints are identified by GUNIT(1,I)=1.

  2. MODE=1 compute only FU, i.e. function values.

  3. MODE=3 compute function and gradient values FU and FUGRAD.

The routines SETUP0, SETUP and SOLCHK must be present in this case too, SETUP and SOLCHK may have an empty body of course.

DONLP2 has a built-in feature for doing gradients numerically. This can be used for both methods of problem description. This is controlled by the variables ANALYT , EPSFCN, TAUBND and DIFFTYPE. If ANALYT=.TRUE. then DONLP2 uses the values from the EGRAD.. routines or from FUGRAD, according to the setting of BLOC. If ANALYT=.FALSE., then numerical differentiation is done internally using a method depending on DIFFTYPE.

  1. DIFFTYPE=1 Use the ordinary forward difference quotient with discretization stepsize 0.1EPSFCN1/2
    componentwise relative.

  2. DIFFTYPE=2 Use the symmetric difference quotient with discretization stepsize 0.1EPSFCN1/3 componentwise relative

  3. DIFFTYPE=3 Use a sixth order approximation computing a Richardson extrapolation of three symmetric difference quotient values. This uses a discretization stepsize 0.01EPSFCN1/7.

Here

EPSFCN

is the expected relative precision of the function evaluation. The precision obtained in the gradients is then of the order of EPSFCN1/2, EPSFCN2/3, EPSFCN6/7 respectively. Since numerical differentiation can also occur if variables are on their bounds, the user must supply a further parameter

TAUBND

which gives a positive value by which any bound present may be violated. The user is warned that numerical differentiation uses n, 2n, 6n additional function evaluations for a single gradient for DIFFTYPE=1,2,3 respectively.

The sequence of routines described above has to be present in any case. (If EVAL EXTERN is used, then of course the bodies of EF, EGRADF, EH, EGRADH, EG, EGRADG are to be empty. This distribution contains such a file donlp2dummyfu.f .) Bound constraints are treated in a special way by DONLP2. The initial point is corrected automatically to fit the bounds. Thereafter any succeeding point is projected with respect to the bounds. The presence of bound constraints is indicated using the GUNIT-array. This is a descriptive array for f, h1 ,..., he , g1 ,..., gm (with e=NH and m=NG) (exactly in that order) with the following meaning:

  1. GUNIT(1,j)=-1: j-th function is "general".
  2. GUNIT(1,j)=1: j-th function depends linearly on X(GUNIT(2,j)) alone with GUNIT(3,j) as partial derivative, that is this function has the form gunit(3,j)*x(gunit(2,j)) + const.

E. g. let

                                  3x7 - 5.6 ³ 0

be the 20-th inequality constraint ( the call EG(20,X,GXI) gives GXI=3x7 - 5.6 ). This will be expressed as

GUNIT(1,NH+20)=1
GUNIT(2,NH+20)=7
GUNIT(3,NH+20)=3

The evaluation of these bound constraints however must be done in EG and EH in the usual way, e.g.

.................................
IF ( I .EQ. 20 ) THEN
  GXI=3.D0*X(7)-5.6D0
  RETURN
ENDIF
.................................

GUNIT must be given completely from (I,0) (corresponding to f ), over (I,1,...NH), (corresponding h) until (I,NH+NG), I=1,2,3. If a user doesn't want to use this feature , she(he) simple writes

DO I=0,NRESM
  GUNIT(1,I)=-1
  GUNIT(2,I)=0
  GUNIT(3,I)=0
ENDDO

within SETUP0. This however may be inefficient, since now bound-constraints are treated like general nonlinear ones and may also become violated(!).

Within SETUP0, the following data must be given:

  1. N number of variables (£ 300 at present. check O8PARA.INC)
     
  2. NH number of equality constraints (may be zero).
     
  3. NG number of inequality constraints including bounds. (may be zero). NH+NG £ NRESM £ 1800 by default setting in O8PARA.INC at present.
  4. NAME Problem identifier, 40 characters maximum. The leading 8 characters (if any) have to be alphanumeric, the first one alphabetic. Nonalphabetic characters are interpreted as 'X'.
  5. X Initial guess (also holds the current solution).
  6. TAU0 !!! This parameter should be carefully selected by the user: It gives a (universal) bound describing how much the unscaled penalty-term (the l1-norm of the constraint violation) may deviate from zero. DONLP2 assumes that within the region described by

                                    NH                 NG
                                    å |hi(x)|   -  å min{0, gi(x)} £  TAU0
                                    i=1                  i=1

    all functions may be evaluated safely. The initial guess however may violate these requirements. In that case an initial feasibility improvement phase is run by DONLP2 until a point is found, that fits them. This is done by scaling f by SCF=0, such that the pure unscaled penalty term is minimized. TAU0 may be chosen as large as a user may want. A small TAU0 diminishes the efficiency of DONLP2, because the iterates then will follow the boundary of the feasible set closely. These remarks do not apply for bound constraints, since DONLP2 treats these by a projection technique. Bounds always remain satisfied. In the example HS114 TAU0=1 works by good luck only. Contrary, a large TAU0 may degrade the reliability of the code (outdoors the wolfes are howling.)
  7. DEL0. In the initial phase of minimization a constraint is considered binding if

                                   gi(x) / max{1,|| Ñgi(x)||} £ DEL0 

    Good values are between 1.d0 and 1.d-2 . If equal to zero , DEL0 is set equal to TAU0. If DEL0 is chosen too small, then identification of the correct set of binding constraints may be delayed. Contrary, if DEL0 is too large, then the method will often escape to the full regularized SQP method, using individual slack variables for any active constraint, which is quite costly. For well scaled problems DEL0=1 is reasonable. However, sometimes it should be much smaller (compare the list of testresults in the accompanying paper).
  8. NRESET. If there are more then NRESET steps using "small" stepsizes and therefore small corrections, a "restart" of the accumulated quasi-Newton-update is tried. NRESET is internally bounded by 4 from below. Good values are between N and 3*N. In the standard setting NRESET is bounded by N. Override this by using SETUP if desired.
  9. ANALYT. ANALYT=.TRUE., if the gradients are given by analytical expressions. With ANALYT=.FALSE. donlp2 relaxes its termination criteria using the variable
  10. EPSDIF. DONLP2 assumes relative precision of EPSDIF in the gradients if ANALYT is set to .FALSE. 
  11. EPSFCN. This is used only if ANALYT=.FALSE.. In this case it must give the relative precision of the function evaluation routine. Since precision of the numerical gradients is lower than this, using the method with very crude function values makes no sense.
  12. DIFFTYPE Type of numerical differentiation to be used, if any. See above.
  13. PROU. FORTRAN unit number for output protocol 
  14. MEU. FORTRAN unit number for messages concerning "special events" (PROU and MEU must be different!!!!)
  15. SILENT. logical. If .TRUE., neither PROU nor MEU are written. DONLP2 returns silently with its results in the appropriate common-areas giving no form of output. See the description of this areas below. The user might extract his desired information from these.
  16. TAUBND. Amount by which bounds may be violated if numerical differentiation is used. 
  17. COLD. If true, then any information is computed afresh . Otherwise, the penalty-weights and quasinewton-udates as stored in the common-fields are used for the current run. This is useful if a parametric problem is to be solved.
  18. GUNIT see above

Default settings are

 TAUBND=1.D0
 EPSFCN=1.D-16
 EPSDIF=1.D-14
 SILENT=.FALSE.
 COLD=.TRUE.

As an example for coding the problem case 114 from Hock and Schittkowski (Lecture notes on economics and mathematical systems 187, alkylation problem from Bracken and McCormick) is given here.

The problem is to minimize

f(x) = 5.04x1 + 0.035x2 + 10x3 + 3.36x5 - 0.063x4x7

subject to

h1(x) = 1.22x4 - x1 -x5 = 0
h2(x) = 98000x3/(x4x9 + 1000x3) - x6 = 0
h3(x) = (x2 + x5)/x1 - x8 = 0
g1(x) = 35.82 - 0.222x10 - bx9 ³ 0, b = 0.9
g2(x) = -133 + 3x7 - ax10 ³ 0, a = 0.99
g3(x) = -g1(x) + x9(1/b - b) ³  0,
g4(x) = -g2(x) + (1/a - a)x10 ³  0,
g5(x) = 1.12x1 + 0.13167x1x8 - 0.00667x1x82 - ax4 ³  0,
g6(x) = 57.425 + 1.098x8 - 0.038x82 + 0.325x6 - ax7 ³  0,
g7(x) = -g5(x) + (1/a - a)x4 ³  0,
g8(x) = -g6(x) + (1/a - a)x7 ³  0

and the bounds

0.00001   £    x1    £ 2000
0.00001   £    x2    £ 16000
0.00001   £    x3     £ 120
0.00001   £    x4    £ 5000
0.00001   £    x5    £ 2000
         85   £    x6    £ 93
         90   £    x7    £ 95
           3   £    x8    £ 12
         1.2  £    x9     £ 4
          45  £    x10  £ 162.

This may be coded as follows (given as alkylati.f in the directory examples): (For easier handling, INCLUDE-files are used here, which are to be found in this distribution too.)

           SUBROUTINE SETUP
C  NO SPECIAL USER DATA
           RETURN
           END
           SUBROUTINE SETUP0
           INCLUDE 'O8COMM.INC'
           INTEGER I,J
           DOUBLE PRECISION XST0(10)
           DATA (XST0(I),I=1,10)
        F            /1745.D0,12.D3,11.D1,3048.D0,1974.D0,89.2D0,92.8D0,8.D0,
        1           3.6D0,145.D0/
C  NAME IS IDENT OF THE EXAMPLE/USER AND CAN BE SET AT USERS WILL
C  THE FIRST CHARACTER MUST BE ALPHABETIC. 40 CHARACTERS MAXIMUM
           NAME='ALKYLATION'
C X IS INITIAL GUESS AND ALSO HOLDS THE CURRENT SOLUTION
C PROBLEM DIMENSION N=DIM(X), NH=DIM(H), NG=DIM(G)
           N=10
           NH=3
           NG=28
           ANALYT=.TRUE.
           EPSDIF=1.D-16
           EPSFCN=1.D-16
           PROU=10
           MEU=20
C DEL0 AND TAU0: SEE BELOW
           DEL0=1.D0
           TAU0=1.D0
           TAU=0.1D0
           DO I=1,N
           X(I)=XST0(I)
           ENDDO
C GUNIT-ARRAY, SEE BELOW
           DO J=0,11
           GUNIT(1,J)=-1
           GUNIT(2,J)=0
           GUNIT(3,J)=0
           ENDDO
           DO J=12,31
           GUNIT(1,J)=1
           IF ( J .LE. 21 ) THEN
              GUNIT(2,J)=J-11
              GUNIT(3,J)=1
           ELSE
              GUNIT(2,J)=J-21
              GUNIT(3,J)=-1
           ENDIF
           ENDDO
           RETURN
           END

C   OBJECTIVE FUNCTION
           SUBROUTINE EF(X,FX)
           INCLUDE 'O8FUCO.INC'
           DOUBLE PRECISION FX,X(*)
           ICF=ICF+1
           FX=5.04D0*X(1)+.035D0*X(2)+10.D0*X(3)+3.36D0*X(5)-.063D0*X(4)*X(7)
           RETURN
           END
C    GRADIENT OF OBJECTIVE FUNCTION
           SUBROUTINE EGRADF(X,GRADF)
           INCLUDE 'O8FUCO.INC'
           INTEGER J
           DOUBLE PRECISION X(*),GRADF(*),A(10)
           SAVE A
           DATA A/5.04D0,0.035D0,10.D0,0.D0,3.36D0,5*0.D0/
           ICGF=ICGF+1
           DO     100       J=1,10
           GRADF(J)=A(J)
100     CONTINUE
           GRADF(4)=-0.063D0*X(7)
           GRADF(7)=-0.063D0*X(4)
           RETURN
           END

C    COMPUTE THE I-TH EQUALITY CONSTAINT, VALUE IS HXI
           SUBROUTINE EH(I,X,HXI)
           INCLUDE 'O8FUCO.INC'
           DOUBLE PRECISION X(*),HXI
           INTEGER I
           CRES(I)=CRES(I)+1
           GOTO (100,200,300),I
100     CONTINUE
           HXI=1.22D0*X(4)-X(1)-X(5)
           RETURN
200     CONTINUE
           HXI=9.8D4*X(3)/(X(4)*X(9)+1.D3*X(3))-X(6)
           RETURN
300     CONTINUE
           HXI=(X(2)+X(5))/X(1)-X(8)
           RETURN
           END
C    COMPUTE THE GRADIENT OF THE I-TH EQUALITY CONSTRAINT
           SUBROUTINE EGRADH(I,X,GRADHI)
           INCLUDE 'O8FUCO.INC'
           INTEGER I,J
           DOUBLE PRECISION X(*),GRADHI(*),T,T1
           CGRES(I)=CGRES(I)+1
           DO     100    J=1,10
           GRADHI(J)=0.D0
100     CONTINUE
           GOTO(200,300,400),I
200     CONTINUE
           GRADHI(1)=-1.D0
           GRADHI(4)=1.22D0
           GRADHI(5)=-1.D0
           RETURN
300     CONTINUE
           T=9.8D4/(X(4)*X(9)+1.D3*X(3))
           T1=T/(X(4)*X(9)+1.D3*X(3))*X(3)
           GRADHI(3)=T-1.D3*T1
           GRADHI(4)=-X(9)*T1
           GRADHI(9)=-X(4)*T1
           GRADHI(6)=-1.D0
           RETURN
400     CONTINUE
           GRADHI(1)=-(X(2)+X(5))/X(1)**2
           GRADHI(2)=1.D0/X(1)
           GRADHI(5)=GRADHI(2)
           GRADHI(8)=-1.D0
           RETURN
           END

C    COMPUTE THE I-TH INEQUALITY CONSTAINT, BOUNDS INCLUDED

           SUBROUTINE EG(I,X,GXI)
           INCLUDE 'O8FUCO.INC'
           INTEGER I,K
           DOUBLE PRECISION GXI
           DOUBLE PRECISION X(NX),OG(10),UG(10),A,B,C,D,T
           SAVE A,B,C,D,UG,OG
           DATA A/.99D0/, B/.9D0/, C/2.01010101010101D-2/, D/
         1              2.11111111111111D-1/
           DATA OG/2.D3,16.D3,1.2D2,5.D3,2.D3,93.D0,95.D0,12.D0,4.D0,162.D0/
           DATA UG/5*1.D-5,85.D0,90.D0,3.D0,1.2D0,145.D0/
           IF ( GUNIT(1,I+NH) .EQ. -1 ) CRES(I+NH)=CRES(I+NH)+1
           IF(I .GT. 8)       GOTO 500 
           K=(I+1)/2
           GOTO(100,200,300,400) K
100     CONTINUE
           T=35.82D0-.222D0*X(10)-B*X(9)
           IF(K+K .EQ. I) T=-T+X(9)*D
           GXI=T
           RETURN
200     CONTINUE
           T=-133.D0+3.D0*X(7)-A*X(10)
           IF(K+K .EQ. I)          T=-T+C*X(10)
           GXI=T
           RETURN
300     CONTINUE
           T=1.12D0*X(1)+.13167D0*X(1)*X(8)-.00667D0*X(1)*X(8)**2-A*X(4)
           IF(K+K .EQ. I)         T=-T+C*X( 4)
           GXI=T
           RETURN
400     CONTINUE
           T=57.425D0+1.098D0*X(8)-.038D0*X(8)**2+.325D0*X(6)-A*X(7)
           IF(K+K .EQ. I)       T=-T+C*X( 7)
           GXI=T
           RETURN
500     CONTINUE
           IF(I .GT. 18)    GXI=OG(I-18)-X(I-18)
           IF(I .LE. 18)     GXI=X(I-8)-UG(I-8)
           END
C    COMPUTE THE GRADIENT OF THE I-TH INEQUALITY CONSTAINT
           SUBROUTINE EGRADG(I,X,GRADGI)
           INCLUDE 'O8FUCO.INC'
           INTEGER I,J,K
           DOUBLE PRECISION X(NX) ,GRADGI(NX),A,B,C,D
           SAVE A,B,C,D
           DATA A/.99D0/, B/.9D0/, C/1.01010101010101D0 /, D/
         1               1.11111111111111D0/
           DO 50 J=1,NX
           GRADGI(J)=0.D0
50       CONTINUE
           IF(I .GT. 8) GOTO 500
           K=(I+1)/2
           GOTO(100,200,300,400) K
100     CONTINUE
           IF(K+K .EQ. I)      GOTO 150
           GRADGI(9)=-B
           GRADGI(10)=-.222D0
           RETURN
150     CONTINUE
           GRADGI( 9)= D
           GRADGI(10)= .222D0
           RETURN
200     CONTINUE
           IF(K+K .EQ. I)     GOTO 250
           GRADGI( 7)=3.D0
           GRADGI(10)=-A
           RETURN
250     CONTINUE
           GRADGI( 7)=-3.D0
           GRADGI(10)= C
           RETURN
300     CONTINUE
           GRADGI( 1)=1.12D0+.13167D0*X(8)-.00667D0*X(8)**2
           GRADGI( 4)= -A GRADGI( 8)=.13167D0*X(1)-.01334D0*X(1)*X(8)
           IF(K+K .NE. I)      RETURN
           GRADGI( 1)=-GRADGI(1)
           GRADGI( 8)=-GRADGI(8)
           GRADGI( 4)= C
           RETURN 
400     CONTINUE
           GRADGI( 6)=.325D0
           GRADGI( 7)= -A
           GRADGI( 8)=1.098D0-.076D0*X(8)
           IF(K+K .NE. I)     RETURN
           GRADGI( 6)=-.325D0
           GRADGI( 7)=C
           GRADGI( 8)=-GRADGI(8)
           RETURN
500    CONTINUE
           IF(I .GT. 18)     GRADGI(I-18)=-1.D0
           IF(I .LE. 18)     GRADGI(I-8)=1.D0
           RETURN
           END

If the user wants valid statistics concerning the number of function evaluations etc. she(he) has to increase the counters in the named common block O8CNT as done in the example above.

  1. ICF counts the number of function calls of the objective,

  2. CRES(1:NH) those of the equality and

  3. CRES(NH+1:NG+NH of the inequality constraints.

  4.  ICGF, and CGRES(1:NH+NG) are the counters of the gradient calls.
DONLP2 does no evaluation statistics for its own.

3.2 ADDITIONAL FEATURES

There are additional possibilities to enhance the performance of the code:

  1. Affine linear functions.
    If some function is affine-linear in x (i.e. has constant gradient) then the user might set GCONST(I)=.TRUE. in SETUP0, with I in the range of 0 to NH+NG. In that case the gradient is evaluated once only and the corresponding function is excluded from the computation of the Maratos-correction. Also the evaluation of this gradient is counted once only. Feasibility for affine linear functions other than bounds and fixed variables is not maintained.
  2. Parametric problems.
    If the user has to solve a problem in a parametric fashion, she/he may do this using COLD=.FALSE. in a subsequent call of SETUP0 , thereby avoiding the standard initialization of DONLP2 (done in O8ST). In that case the old Quasi-Newton-update and the old penalty-weights are used as initial guesses. The files USER1.INC, parmain.f and parametric.f show how this may be done.
  3. Parameter settings.
    Any of the parameters of DONLP2 (e.g. the parameters controlling amount of output TE0, TE1, TE2, TE3 ) may be changed at the users will within SETUP by redefining the values located in the corresponding common blocks described below. This must be done within SETUP whose call follows that of the standard initialization routine O8ST. In SETUP use INCLUDE 'O8COMM.INC' in order to access these parameters.
  4. Error return from function evaluation.
    The user has the possibility to prevent DONLP2 from leaving the domain of definition of one of her(his) functions using the error-indicators of the COMMON-block O8ERR. The code proposes changes of the current solution using some formula

                                                      xnew = xcurrent + s d + s 2dd

    with correction vectors d and dd. These directions might well be "downhill" but too large, hence xnew leaving the domain of definition of some function involved. The user may check the actual parameter of the function evaluation and set FFUERR=.TRUE. (in the objective function evaluation) or CFUERR(I)=.TRUE. in one of the constraint function evaluations and return to the optimizer then. In this case DONLP2 tries to avoid the problem by reducing s. If such an error occurs with the current x, i.e. s = 0, the run is terminated with an appropriate error message.
  5. Scaling.
    The user has the possibility to prescribe an internal scaling of x, setting the array XSC appropriately in SETUP0 (must be done there, if any). The initial value given and the function and gradient evaluations however are always in the original unscaled variables. The external variables are then XSC(I)*X(I), with X(I) as internal variables. The first internal variable is obtained by dividing the given initial values X(I) from SETUP0 or BLOCKDATA by XSC(I).
  6. Changing termination criteria.
    Termination criteria can be changed in various ways. Seven variables enter these criteria (see below), namely

                            EPSX, DELMIN, SMALLW, EPSDIF, NRESET, NUMSM , EPSPHI

    The first five have been mentioned already. The remaining two may be changed from their default values MAX(10,N) and 1000*EPSMAC in order to avoid many steps with only very marginally changing objective function values in the terminal phase of optimization. (Such a situation occurs e.g. for illconditioned cases). The user is warned that being too sloppy with EPSPHI might result in premature termination and solutions far from optimum, if the problem is nonconvex.