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.
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.
- EF(X,FX) returns FX=f(x) given x as input. Arguments DOUBLE PRECISION X(*), FX
- EGRADF(X,GRADF) returns GRADF= Ñ f(x) given x as input. Arguments DOUBLE PRECISION X(*),GRADF(*).
- 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.
- EGRADH(I,X,GRADHI) returns Ñhi(x) as GRADHI given i and x as input. Arguments INTEGER I ; DOUBLE PRECISION X(*),GRADHI(*).
- 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.
- 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.
- 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.
- 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.
- 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:
- 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.
- MODE=1 compute only FU, i.e. function values.
- 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.
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:
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
.................................
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:
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.
There are additional possibilities to enhance the performance of the code: