12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412 |
- COBYLA
- ~~~~~~
- Here is a single-precision Fortran implementation of the algorithm for
- constrained optimization that is the subject of the report I have written on
- "A direct search optimization method that models the objective and constraint
- functions by linear interpolation". This report has the number DAMTP 1992/NA5,
- University of Cambridge, and it has been published in the proceedings of the
- conference on Numerical Analysis and Optimization that was held in Oaxaca,
- Mexico in January, 1992, which is the book "Advances in Optimization and
- Numerical Analysis" (eds. Susana Gomez and Jean-Pierre Hennart), Kluwer
- Academic Publishers (1994).
- The instructions for using the Fortran code are given in the comments of
- SUBROUTINE COBYLA, which is the interface between the user and the main
- calculation that is done by SUBROUTINE COBYLB. There is a need for a linear
- programming problem to be solved subject to a Euclidean norm trust region
- constraint. Therefore SUBROUTINE TRSTLP is provided too, but you may have some
- software that you prefer to use instead. These 3 subroutines are separated by
- lines of hyphens below. Further, there follows the main program, the CALCFC
- subroutine and the output that are appropriate to the numerical examples that
- are discussed in the last section of DAMTP 1992/NA5. Please note, however,
- that some cosmetic restructuring of the software has caused the given output
- to differ slightly from Table 1 of the report.
- There are no restrictions on the use of the software, nor do I offer any
- guarantees of success. Indeed, at the time of writing this note I had applied
- it only to test problems that have up to 10 variables.
- Mike Powell (May 7th, 1992).
- -------------------------------------------------------------------------------
- SUBROUTINE COBYLA (N,M,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W,IACT)
- DIMENSION X(*),W(*),IACT(*)
- C
- C This subroutine minimizes an objective function F(X) subject to M
- C inequality constraints on X, where X is a vector of variables that has
- C N components. The algorithm employs linear approximations to the
- C objective and constraint functions, the approximations being formed by
- C linear interpolation at N+1 points in the space of the variables.
- C We regard these interpolation points as vertices of a simplex. The
- C parameter RHO controls the size of the simplex and it is reduced
- C automatically from RHOBEG to RHOEND. For each RHO the subroutine tries
- C to achieve a good vector of variables for the current size, and then
- C RHO is reduced until the value RHOEND is reached. Therefore RHOBEG and
- C RHOEND should be set to reasonable initial changes to and the required
- C accuracy in the variables respectively, but this accuracy should be
- C viewed as a subject for experimentation because it is not guaranteed.
- C The subroutine has an advantage over many of its competitors, however,
- C which is that it treats each constraint individually when calculating
- C a change to the variables, instead of lumping the constraints together
- C into a single penalty function. The name of the subroutine is derived
- C from the phrase Constrained Optimization BY Linear Approximations.
- C
- C The user must set the values of N, M, RHOBEG and RHOEND, and must
- C provide an initial vector of variables in X. Further, the value of
- C IPRINT should be set to 0, 1, 2 or 3, which controls the amount of
- C printing during the calculation. Specifically, there is no output if
- C IPRINT=0 and there is output only at the end of the calculation if
- C IPRINT=1. Otherwise each new value of RHO and SIGMA is printed.
- C Further, the vector of variables and some function information are
- C given either when RHO is reduced or when each new value of F(X) is
- C computed in the cases IPRINT=2 or IPRINT=3 respectively. Here SIGMA
- C is a penalty parameter, it being assumed that a change to X is an
- C improvement if it reduces the merit function
- C F(X)+SIGMA*MAX(0.0,-C1(X),-C2(X),...,-CM(X)),
- C where C1,C2,...,CM denote the constraint functions that should become
- C nonnegative eventually, at least to the precision of RHOEND. In the
- C printed output the displayed term that is multiplied by SIGMA is
- C called MAXCV, which stands for 'MAXimum Constraint Violation'. The
- C argument MAXFUN is an integer variable that must be set by the user to a
- C limit on the number of calls of CALCFC, the purpose of this routine being
- C given below. The value of MAXFUN will be altered to the number of calls
- C of CALCFC that are made. The arguments W and IACT provide real and
- C integer arrays that are used as working space. Their lengths must be at
- C least N*(3*N+2*M+11)+4*M+6 and M+1 respectively.
- C
- C In order to define the objective and constraint functions, we require
- C a subroutine that has the name and arguments
- C SUBROUTINE CALCFC (N,M,X,F,CON)
- C DIMENSION X(*),CON(*) .
- C The values of N and M are fixed and have been defined already, while
- C X is now the current vector of variables. The subroutine should return
- C the objective and constraint functions at X in F and CON(1),CON(2),
- C ...,CON(M). Note that we are trying to adjust X so that F(X) is as
- C small as possible subject to the constraint functions being nonnegative.
- C
- C Partition the working space array W to provide the storage that is needed
- C for the main calculation.
- C
- MPP=M+2
- ICON=1
- ISIM=ICON+MPP
- ISIMI=ISIM+N*N+N
- IDATM=ISIMI+N*N
- IA=IDATM+N*MPP+MPP
- IVSIG=IA+M*N+N
- IVETA=IVSIG+N
- ISIGB=IVETA+N
- IDX=ISIGB+N
- IWORK=IDX+N
- CALL COBYLB (N,M,MPP,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W(ICON),
- 1 W(ISIM),W(ISIMI),W(IDATM),W(IA),W(IVSIG),W(IVETA),W(ISIGB),
- 2 W(IDX),W(IWORK),IACT)
- RETURN
- END
- C------------------------------------------------------------------------------
- SUBROUTINE COBYLB (N,M,MPP,X,RHOBEG,RHOEND,IPRINT,MAXFUN,
- 1 CON,SIM,SIMI,DATMAT,A,VSIG,VETA,SIGBAR,DX,W,IACT)
- DIMENSION X(*),CON(*),SIM(N,*),SIMI(N,*),DATMAT(MPP,*),
- 1 A(N,*),VSIG(*),VETA(*),SIGBAR(*),DX(*),W(*),IACT(*)
- C
- C Set the initial values of some parameters. The last column of SIM holds
- C the optimal vertex of the current simplex, and the preceding N columns
- C hold the displacements from the optimal vertex to the other vertices.
- C Further, SIMI holds the inverse of the matrix that is contained in the
- C first N columns of SIM.
- C
- IPTEM=MIN0(N,5)
- IPTEMP=IPTEM+1
- NP=N+1
- MP=M+1
- ALPHA=0.25
- BETA=2.1
- GAMMA=0.5
- DELTA=1.1
- RHO=RHOBEG
- PARMU=0.0
- IF (IPRINT .GE. 2) PRINT 10, RHO
- 10 FORMAT (/3X,'The initial value of RHO is',1PE13.6,2X,
- 1 'and PARMU is set to zero.')
- NFVALS=0
- TEMP=1.0/RHO
- DO 30 I=1,N
- SIM(I,NP)=X(I)
- DO 20 J=1,N
- SIM(I,J)=0.0
- 20 SIMI(I,J)=0.0
- SIM(I,I)=RHO
- 30 SIMI(I,I)=TEMP
- JDROP=NP
- IBRNCH=0
- C
- C Make the next call of the user-supplied subroutine CALCFC. These
- C instructions are also used for calling CALCFC during the iterations of
- C the algorithm.
- C
- 40 IF (NFVALS .GE. MAXFUN .AND. NFVALS .GT. 0) THEN
- IF (IPRINT .GE. 1) PRINT 50
- 50 FORMAT (/3X,'Return from subroutine COBYLA because the ',
- 1 'MAXFUN limit has been reached.')
- GOTO 600
- END IF
- NFVALS=NFVALS+1
- CALL CALCFC (N,M,X,F,CON)
- RESMAX=0.0
- IF (M .GT. 0) THEN
- DO 60 K=1,M
- 60 RESMAX=AMAX1(RESMAX,-CON(K))
- END IF
- IF (NFVALS .EQ. IPRINT-1 .OR. IPRINT .EQ. 3) THEN
- PRINT 70, NFVALS,F,RESMAX,(X(I),I=1,IPTEM)
- 70 FORMAT (/3X,'NFVALS =',I5,3X,'F =',1PE13.6,4X,'MAXCV =',
- 1 1PE13.6/3X,'X =',1PE13.6,1P4E15.6)
- IF (IPTEM .LT. N) PRINT 80, (X(I),I=IPTEMP,N)
- 80 FORMAT (1PE19.6,1P4E15.6)
- END IF
- CON(MP)=F
- CON(MPP)=RESMAX
- IF (IBRNCH .EQ. 1) GOTO 440
- C
- C Set the recently calculated function values in a column of DATMAT. This
- C array has a column for each vertex of the current simplex, the entries of
- C each column being the values of the constraint functions (if any)
- C followed by the objective function and the greatest constraint violation
- C at the vertex.
- C
- DO 90 K=1,MPP
- 90 DATMAT(K,JDROP)=CON(K)
- IF (NFVALS .GT. NP) GOTO 130
- C
- C Exchange the new vertex of the initial simplex with the optimal vertex if
- C necessary. Then, if the initial simplex is not complete, pick its next
- C vertex and calculate the function values there.
- C
- IF (JDROP .LE. N) THEN
- IF (DATMAT(MP,NP) .LE. F) THEN
- X(JDROP)=SIM(JDROP,NP)
- ELSE
- SIM(JDROP,NP)=X(JDROP)
- DO 100 K=1,MPP
- DATMAT(K,JDROP)=DATMAT(K,NP)
- 100 DATMAT(K,NP)=CON(K)
- DO 120 K=1,JDROP
- SIM(JDROP,K)=-RHO
- TEMP=0.0
- DO 110 I=K,JDROP
- 110 TEMP=TEMP-SIMI(I,K)
- 120 SIMI(JDROP,K)=TEMP
- END IF
- END IF
- IF (NFVALS .LE. N) THEN
- JDROP=NFVALS
- X(JDROP)=X(JDROP)+RHO
- GOTO 40
- END IF
- 130 IBRNCH=1
- C
- C Identify the optimal vertex of the current simplex.
- C
- 140 PHIMIN=DATMAT(MP,NP)+PARMU*DATMAT(MPP,NP)
- NBEST=NP
- DO 150 J=1,N
- TEMP=DATMAT(MP,J)+PARMU*DATMAT(MPP,J)
- IF (TEMP .LT. PHIMIN) THEN
- NBEST=J
- PHIMIN=TEMP
- ELSE IF (TEMP .EQ. PHIMIN .AND. PARMU .EQ. 0.0) THEN
- IF (DATMAT(MPP,J) .LT. DATMAT(MPP,NBEST)) NBEST=J
- END IF
- 150 CONTINUE
- C
- C Switch the best vertex into pole position if it is not there already,
- C and also update SIM, SIMI and DATMAT.
- C
- IF (NBEST .LE. N) THEN
- DO 160 I=1,MPP
- TEMP=DATMAT(I,NP)
- DATMAT(I,NP)=DATMAT(I,NBEST)
- 160 DATMAT(I,NBEST)=TEMP
- DO 180 I=1,N
- TEMP=SIM(I,NBEST)
- SIM(I,NBEST)=0.0
- SIM(I,NP)=SIM(I,NP)+TEMP
- TEMPA=0.0
- DO 170 K=1,N
- SIM(I,K)=SIM(I,K)-TEMP
- 170 TEMPA=TEMPA-SIMI(K,I)
- 180 SIMI(NBEST,I)=TEMPA
- END IF
- C
- C Make an error return if SIGI is a poor approximation to the inverse of
- C the leading N by N submatrix of SIG.
- C
- ERROR=0.0
- DO 200 I=1,N
- DO 200 J=1,N
- TEMP=0.0
- IF (I .EQ. J) TEMP=TEMP-1.0
- DO 190 K=1,N
- 190 TEMP=TEMP+SIMI(I,K)*SIM(K,J)
- 200 ERROR=AMAX1(ERROR,ABS(TEMP))
- IF (ERROR .GT. 0.1) THEN
- IF (IPRINT .GE. 1) PRINT 210
- 210 FORMAT (/3X,'Return from subroutine COBYLA because ',
- 1 'rounding errors are becoming damaging.')
- GOTO 600
- END IF
- C
- C Calculate the coefficients of the linear approximations to the objective
- C and constraint functions, placing minus the objective function gradient
- C after the constraint gradients in the array A. The vector W is used for
- C working space.
- C
- DO 240 K=1,MP
- CON(K)=-DATMAT(K,NP)
- DO 220 J=1,N
- 220 W(J)=DATMAT(K,J)+CON(K)
- DO 240 I=1,N
- TEMP=0.0
- DO 230 J=1,N
- 230 TEMP=TEMP+W(J)*SIMI(J,I)
- IF (K .EQ. MP) TEMP=-TEMP
- 240 A(I,K)=TEMP
- C
- C Calculate the values of sigma and eta, and set IFLAG=0 if the current
- C simplex is not acceptable.
- C
- IFLAG=1
- PARSIG=ALPHA*RHO
- PARETA=BETA*RHO
- DO 260 J=1,N
- WSIG=0.0
- WETA=0.0
- DO 250 I=1,N
- WSIG=WSIG+SIMI(J,I)**2
- 250 WETA=WETA+SIM(I,J)**2
- VSIG(J)=1.0/SQRT(WSIG)
- VETA(J)=SQRT(WETA)
- IF (VSIG(J) .LT. PARSIG .OR. VETA(J) .GT. PARETA) IFLAG=0
- 260 CONTINUE
- C
- C If a new vertex is needed to improve acceptability, then decide which
- C vertex to drop from the simplex.
- C
- IF (IBRNCH .EQ. 1 .OR. IFLAG .EQ. 1) GOTO 370
- JDROP=0
- TEMP=PARETA
- DO 270 J=1,N
- IF (VETA(J) .GT. TEMP) THEN
- JDROP=J
- TEMP=VETA(J)
- END IF
- 270 CONTINUE
- IF (JDROP .EQ. 0) THEN
- DO 280 J=1,N
- IF (VSIG(J) .LT. TEMP) THEN
- JDROP=J
- TEMP=VSIG(J)
- END IF
- 280 CONTINUE
- END IF
- C
- C Calculate the step to the new vertex and its sign.
- C
- TEMP=GAMMA*RHO*VSIG(JDROP)
- DO 290 I=1,N
- 290 DX(I)=TEMP*SIMI(JDROP,I)
- CVMAXP=0.0
- CVMAXM=0.0
- DO 310 K=1,MP
- SUM=0.0
- DO 300 I=1,N
- 300 SUM=SUM+A(I,K)*DX(I)
- IF (K .LT. MP) THEN
- TEMP=DATMAT(K,NP)
- CVMAXP=AMAX1(CVMAXP,-SUM-TEMP)
- CVMAXM=AMAX1(CVMAXM,SUM-TEMP)
- END IF
- 310 CONTINUE
- DXSIGN=1.0
- IF (PARMU*(CVMAXP-CVMAXM) .GT. SUM+SUM) DXSIGN=-1.0
- C
- C Update the elements of SIM and SIMI, and set the next X.
- C
- TEMP=0.0
- DO 320 I=1,N
- DX(I)=DXSIGN*DX(I)
- SIM(I,JDROP)=DX(I)
- 320 TEMP=TEMP+SIMI(JDROP,I)*DX(I)
- DO 330 I=1,N
- 330 SIMI(JDROP,I)=SIMI(JDROP,I)/TEMP
- DO 360 J=1,N
- IF (J .NE. JDROP) THEN
- TEMP=0.0
- DO 340 I=1,N
- 340 TEMP=TEMP+SIMI(J,I)*DX(I)
- DO 350 I=1,N
- 350 SIMI(J,I)=SIMI(J,I)-TEMP*SIMI(JDROP,I)
- END IF
- 360 X(J)=SIM(J,NP)+DX(J)
- GOTO 40
- C
- C Calculate DX=x(*)-x(0). Branch if the length of DX is less than 0.5*RHO.
- C
- 370 IZ=1
- IZDOTA=IZ+N*N
- IVMC=IZDOTA+N
- ISDIRN=IVMC+MP
- IDXNEW=ISDIRN+N
- IVMD=IDXNEW+N
- CALL TRSTLP (N,M,A,CON,RHO,DX,IFULL,IACT,W(IZ),W(IZDOTA),
- 1 W(IVMC),W(ISDIRN),W(IDXNEW),W(IVMD))
- IF (IFULL .EQ. 0) THEN
- TEMP=0.0
- DO 380 I=1,N
- 380 TEMP=TEMP+DX(I)**2
- IF (TEMP .LT. 0.25*RHO*RHO) THEN
- IBRNCH=1
- GOTO 550
- END IF
- END IF
- C
- C Predict the change to F and the new maximum constraint violation if the
- C variables are altered from x(0) to x(0)+DX.
- C
- RESNEW=0.0
- CON(MP)=0.0
- DO 400 K=1,MP
- SUM=CON(K)
- DO 390 I=1,N
- 390 SUM=SUM-A(I,K)*DX(I)
- IF (K .LT. MP) RESNEW=AMAX1(RESNEW,SUM)
- 400 CONTINUE
- C
- C Increase PARMU if necessary and branch back if this change alters the
- C optimal vertex. Otherwise PREREM and PREREC will be set to the predicted
- C reductions in the merit function and the maximum constraint violation
- C respectively.
- C
- BARMU=0.0
- PREREC=DATMAT(MPP,NP)-RESNEW
- IF (PREREC .GT. 0.0) BARMU=SUM/PREREC
- IF (PARMU .LT. 1.5*BARMU) THEN
- PARMU=2.0*BARMU
- IF (IPRINT .GE. 2) PRINT 410, PARMU
- 410 FORMAT (/3X,'Increase in PARMU to',1PE13.6)
- PHI=DATMAT(MP,NP)+PARMU*DATMAT(MPP,NP)
- DO 420 J=1,N
- TEMP=DATMAT(MP,J)+PARMU*DATMAT(MPP,J)
- IF (TEMP .LT. PHI) GOTO 140
- IF (TEMP .EQ. PHI .AND. PARMU .EQ. 0.0) THEN
- IF (DATMAT(MPP,J) .LT. DATMAT(MPP,NP)) GOTO 140
- END IF
- 420 CONTINUE
- END IF
- PREREM=PARMU*PREREC-SUM
- C
- C Calculate the constraint and objective functions at x(*). Then find the
- C actual reduction in the merit function.
- C
- DO 430 I=1,N
- 430 X(I)=SIM(I,NP)+DX(I)
- IBRNCH=1
- GOTO 40
- 440 VMOLD=DATMAT(MP,NP)+PARMU*DATMAT(MPP,NP)
- VMNEW=F+PARMU*RESMAX
- TRURED=VMOLD-VMNEW
- IF (PARMU .EQ. 0.0 .AND. F .EQ. DATMAT(MP,NP)) THEN
- PREREM=PREREC
- TRURED=DATMAT(MPP,NP)-RESMAX
- END IF
- C
- C Begin the operations that decide whether x(*) should replace one of the
- C vertices of the current simplex, the change being mandatory if TRURED is
- C positive. Firstly, JDROP is set to the index of the vertex that is to be
- C replaced.
- C
- RATIO=0.0
- IF (TRURED .LE. 0.0) RATIO=1.0
- JDROP=0
- DO 460 J=1,N
- TEMP=0.0
- DO 450 I=1,N
- 450 TEMP=TEMP+SIMI(J,I)*DX(I)
- TEMP=ABS(TEMP)
- IF (TEMP .GT. RATIO) THEN
- JDROP=J
- RATIO=TEMP
- END IF
- 460 SIGBAR(J)=TEMP*VSIG(J)
- C
- C Calculate the value of ell.
- C
- EDGMAX=DELTA*RHO
- L=0
- DO 480 J=1,N
- IF (SIGBAR(J) .GE. PARSIG .OR. SIGBAR(J) .GE. VSIG(J)) THEN
- TEMP=VETA(J)
- IF (TRURED .GT. 0.0) THEN
- TEMP=0.0
- DO 470 I=1,N
- 470 TEMP=TEMP+(DX(I)-SIM(I,J))**2
- TEMP=SQRT(TEMP)
- END IF
- IF (TEMP .GT. EDGMAX) THEN
- L=J
- EDGMAX=TEMP
- END IF
- END IF
- 480 CONTINUE
- IF (L .GT. 0) JDROP=L
- IF (JDROP .EQ. 0) GOTO 550
- C
- C Revise the simplex by updating the elements of SIM, SIMI and DATMAT.
- C
- TEMP=0.0
- DO 490 I=1,N
- SIM(I,JDROP)=DX(I)
- 490 TEMP=TEMP+SIMI(JDROP,I)*DX(I)
- DO 500 I=1,N
- 500 SIMI(JDROP,I)=SIMI(JDROP,I)/TEMP
- DO 530 J=1,N
- IF (J .NE. JDROP) THEN
- TEMP=0.0
- DO 510 I=1,N
- 510 TEMP=TEMP+SIMI(J,I)*DX(I)
- DO 520 I=1,N
- 520 SIMI(J,I)=SIMI(J,I)-TEMP*SIMI(JDROP,I)
- END IF
- 530 CONTINUE
- DO 540 K=1,MPP
- 540 DATMAT(K,JDROP)=CON(K)
- C
- C Branch back for further iterations with the current RHO.
- C
- IF (TRURED .GT. 0.0 .AND. TRURED .GE. 0.1*PREREM) GOTO 140
- 550 IF (IFLAG .EQ. 0) THEN
- IBRNCH=0
- GOTO 140
- END IF
- C
- C Otherwise reduce RHO if it is not at its least value and reset PARMU.
- C
- IF (RHO .GT. RHOEND) THEN
- RHO=0.5*RHO
- IF (RHO .LE. 1.5*RHOEND) RHO=RHOEND
- IF (PARMU .GT. 0.0) THEN
- DENOM=0.0
- DO 570 K=1,MP
- CMIN=DATMAT(K,NP)
- CMAX=CMIN
- DO 560 I=1,N
- CMIN=AMIN1(CMIN,DATMAT(K,I))
- 560 CMAX=AMAX1(CMAX,DATMAT(K,I))
- IF (K .LE. M .AND. CMIN .LT. 0.5*CMAX) THEN
- TEMP=AMAX1(CMAX,0.0)-CMIN
- IF (DENOM .LE. 0.0) THEN
- DENOM=TEMP
- ELSE
- DENOM=AMIN1(DENOM,TEMP)
- END IF
- END IF
- 570 CONTINUE
- IF (DENOM .EQ. 0.0) THEN
- PARMU=0.0
- ELSE IF (CMAX-CMIN .LT. PARMU*DENOM) THEN
- PARMU=(CMAX-CMIN)/DENOM
- END IF
- END IF
- IF (IPRINT .GE. 2) PRINT 580, RHO,PARMU
- 580 FORMAT (/3X,'Reduction in RHO to',1PE13.6,' and PARMU =',
- 1 1PE13.6)
- IF (IPRINT .EQ. 2) THEN
- PRINT 70, NFVALS,DATMAT(MP,NP),DATMAT(MPP,NP),
- 1 (SIM(I,NP),I=1,IPTEM)
- IF (IPTEM .LT. N) PRINT 80, (X(I),I=IPTEMP,N)
- END IF
- GOTO 140
- END IF
- C
- C Return the best calculated values of the variables.
- C
- IF (IPRINT .GE. 1) PRINT 590
- 590 FORMAT (/3X,'Normal return from subroutine COBYLA')
- IF (IFULL .EQ. 1) GOTO 620
- 600 DO 610 I=1,N
- 610 X(I)=SIM(I,NP)
- F=DATMAT(MP,NP)
- RESMAX=DATMAT(MPP,NP)
- 620 IF (IPRINT .GE. 1) THEN
- PRINT 70, NFVALS,F,RESMAX,(X(I),I=1,IPTEM)
- IF (IPTEM .LT. N) PRINT 80, (X(I),I=IPTEMP,N)
- END IF
- MAXFUN=NFVALS
- RETURN
- END
- C------------------------------------------------------------------------------
- SUBROUTINE TRSTLP (N,M,A,B,RHO,DX,IFULL,IACT,Z,ZDOTA,VMULTC,
- 1 SDIRN,DXNEW,VMULTD)
- DIMENSION A(N,*),B(*),DX(*),IACT(*),Z(N,*),ZDOTA(*),
- 1 VMULTC(*),SDIRN(*),DXNEW(*),VMULTD(*)
- C
- C This subroutine calculates an N-component vector DX by applying the
- C following two stages. In the first stage, DX is set to the shortest
- C vector that minimizes the greatest violation of the constraints
- C A(1,K)*DX(1)+A(2,K)*DX(2)+...+A(N,K)*DX(N) .GE. B(K), K=2,3,...,M,
- C subject to the Euclidean length of DX being at most RHO. If its length is
- C strictly less than RHO, then we use the resultant freedom in DX to
- C minimize the objective function
- C -A(1,M+1)*DX(1)-A(2,M+1)*DX(2)-...-A(N,M+1)*DX(N)
- C subject to no increase in any greatest constraint violation. This
- C notation allows the gradient of the objective function to be regarded as
- C the gradient of a constraint. Therefore the two stages are distinguished
- C by MCON .EQ. M and MCON .GT. M respectively. It is possible that a
- C degeneracy may prevent DX from attaining the target length RHO. Then the
- C value IFULL=0 would be set, but usually IFULL=1 on return.
- C
- C In general NACT is the number of constraints in the active set and
- C IACT(1),...,IACT(NACT) are their indices, while the remainder of IACT
- C contains a permutation of the remaining constraint indices. Further, Z is
- C an orthogonal matrix whose first NACT columns can be regarded as the
- C result of Gram-Schmidt applied to the active constraint gradients. For
- C J=1,2,...,NACT, the number ZDOTA(J) is the scalar product of the J-th
- C column of Z with the gradient of the J-th active constraint. DX is the
- C current vector of variables and here the residuals of the active
- C constraints should be zero. Further, the active constraints have
- C nonnegative Lagrange multipliers that are held at the beginning of
- C VMULTC. The remainder of this vector holds the residuals of the inactive
- C constraints at DX, the ordering of the components of VMULTC being in
- C agreement with the permutation of the indices of the constraints that is
- C in IACT. All these residuals are nonnegative, which is achieved by the
- C shift RESMAX that makes the least residual zero.
- C
- C Initialize Z and some other variables. The value of RESMAX will be
- C appropriate to DX=0, while ICON will be the index of a most violated
- C constraint if RESMAX is positive. Usually during the first stage the
- C vector SDIRN gives a search direction that reduces all the active
- C constraint violations by one simultaneously.
- C
- IFULL=1
- MCON=M
- NACT=0
- RESMAX=0.0
- DO 20 I=1,N
- DO 10 J=1,N
- 10 Z(I,J)=0.0
- Z(I,I)=1.0
- 20 DX(I)=0.0
- IF (M .GE. 1) THEN
- DO 30 K=1,M
- IF (B(K) .GT. RESMAX) THEN
- RESMAX=B(K)
- ICON=K
- END IF
- 30 CONTINUE
- DO 40 K=1,M
- IACT(K)=K
- 40 VMULTC(K)=RESMAX-B(K)
- END IF
- IF (RESMAX .EQ. 0.0) GOTO 480
- DO 50 I=1,N
- 50 SDIRN(I)=0.0
- C
- C End the current stage of the calculation if 3 consecutive iterations
- C have either failed to reduce the best calculated value of the objective
- C function or to increase the number of active constraints since the best
- C value was calculated. This strategy prevents cycling, but there is a
- C remote possibility that it will cause premature termination.
- C
- 60 OPTOLD=0.0
- ICOUNT=0
- 70 IF (MCON .EQ. M) THEN
- OPTNEW=RESMAX
- ELSE
- OPTNEW=0.0
- DO 80 I=1,N
- 80 OPTNEW=OPTNEW-DX(I)*A(I,MCON)
- END IF
- IF (ICOUNT .EQ. 0 .OR. OPTNEW .LT. OPTOLD) THEN
- OPTOLD=OPTNEW
- NACTX=NACT
- ICOUNT=3
- ELSE IF (NACT .GT. NACTX) THEN
- NACTX=NACT
- ICOUNT=3
- ELSE
- ICOUNT=ICOUNT-1
- IF (ICOUNT .EQ. 0) GOTO 490
- END IF
- C
- C If ICON exceeds NACT, then we add the constraint with index IACT(ICON) to
- C the active set. Apply Givens rotations so that the last N-NACT-1 columns
- C of Z are orthogonal to the gradient of the new constraint, a scalar
- C product being set to zero if its nonzero value could be due to computer
- C rounding errors. The array DXNEW is used for working space.
- C
- IF (ICON .LE. NACT) GOTO 260
- KK=IACT(ICON)
- DO 90 I=1,N
- 90 DXNEW(I)=A(I,KK)
- TOT=0.0
- K=N
- 100 IF (K .GT. NACT) THEN
- SP=0.0
- SPABS=0.0
- DO 110 I=1,N
- TEMP=Z(I,K)*DXNEW(I)
- SP=SP+TEMP
- 110 SPABS=SPABS+ABS(TEMP)
- ACCA=SPABS+0.1*ABS(SP)
- ACCB=SPABS+0.2*ABS(SP)
- IF (SPABS .GE. ACCA .OR. ACCA .GE. ACCB) SP=0.0
- IF (TOT .EQ. 0.0) THEN
- TOT=SP
- ELSE
- KP=K+1
- TEMP=SQRT(SP*SP+TOT*TOT)
- ALPHA=SP/TEMP
- BETA=TOT/TEMP
- TOT=TEMP
- DO 120 I=1,N
- TEMP=ALPHA*Z(I,K)+BETA*Z(I,KP)
- Z(I,KP)=ALPHA*Z(I,KP)-BETA*Z(I,K)
- 120 Z(I,K)=TEMP
- END IF
- K=K-1
- GOTO 100
- END IF
- C
- C Add the new constraint if this can be done without a deletion from the
- C active set.
- C
- IF (TOT .NE. 0.0) THEN
- NACT=NACT+1
- ZDOTA(NACT)=TOT
- VMULTC(ICON)=VMULTC(NACT)
- VMULTC(NACT)=0.0
- GOTO 210
- END IF
- C
- C The next instruction is reached if a deletion has to be made from the
- C active set in order to make room for the new active constraint, because
- C the new constraint gradient is a linear combination of the gradients of
- C the old active constraints. Set the elements of VMULTD to the multipliers
- C of the linear combination. Further, set IOUT to the index of the
- C constraint to be deleted, but branch if no suitable index can be found.
- C
- RATIO=-1.0
- K=NACT
- 130 ZDOTV=0.0
- ZDVABS=0.0
- DO 140 I=1,N
- TEMP=Z(I,K)*DXNEW(I)
- ZDOTV=ZDOTV+TEMP
- 140 ZDVABS=ZDVABS+ABS(TEMP)
- ACCA=ZDVABS+0.1*ABS(ZDOTV)
- ACCB=ZDVABS+0.2*ABS(ZDOTV)
- IF (ZDVABS .LT. ACCA .AND. ACCA .LT. ACCB) THEN
- TEMP=ZDOTV/ZDOTA(K)
- IF (TEMP .GT. 0.0 .AND. IACT(K) .LE. M) THEN
- TEMPA=VMULTC(K)/TEMP
- IF (RATIO .LT. 0.0 .OR. TEMPA .LT. RATIO) THEN
- RATIO=TEMPA
- IOUT=K
- END IF
- END IF
- IF (K .GE. 2) THEN
- KW=IACT(K)
- DO 150 I=1,N
- 150 DXNEW(I)=DXNEW(I)-TEMP*A(I,KW)
- END IF
- VMULTD(K)=TEMP
- ELSE
- VMULTD(K)=0.0
- END IF
- K=K-1
- IF (K .GT. 0) GOTO 130
- IF (RATIO .LT. 0.0) GOTO 490
- C
- C Revise the Lagrange multipliers and reorder the active constraints so
- C that the one to be replaced is at the end of the list. Also calculate the
- C new value of ZDOTA(NACT) and branch if it is not acceptable.
- C
- DO 160 K=1,NACT
- 160 VMULTC(K)=AMAX1(0.0,VMULTC(K)-RATIO*VMULTD(K))
- IF (ICON .LT. NACT) THEN
- ISAVE=IACT(ICON)
- VSAVE=VMULTC(ICON)
- K=ICON
- 170 KP=K+1
- KW=IACT(KP)
- SP=0.0
- DO 180 I=1,N
- 180 SP=SP+Z(I,K)*A(I,KW)
- TEMP=SQRT(SP*SP+ZDOTA(KP)**2)
- ALPHA=ZDOTA(KP)/TEMP
- BETA=SP/TEMP
- ZDOTA(KP)=ALPHA*ZDOTA(K)
- ZDOTA(K)=TEMP
- DO 190 I=1,N
- TEMP=ALPHA*Z(I,KP)+BETA*Z(I,K)
- Z(I,KP)=ALPHA*Z(I,K)-BETA*Z(I,KP)
- 190 Z(I,K)=TEMP
- IACT(K)=KW
- VMULTC(K)=VMULTC(KP)
- K=KP
- IF (K .LT. NACT) GOTO 170
- IACT(K)=ISAVE
- VMULTC(K)=VSAVE
- END IF
- TEMP=0.0
- DO 200 I=1,N
- 200 TEMP=TEMP+Z(I,NACT)*A(I,KK)
- IF (TEMP .EQ. 0.0) GOTO 490
- ZDOTA(NACT)=TEMP
- VMULTC(ICON)=0.0
- VMULTC(NACT)=RATIO
- C
- C Update IACT and ensure that the objective function continues to be
- C treated as the last active constraint when MCON>M.
- C
- 210 IACT(ICON)=IACT(NACT)
- IACT(NACT)=KK
- IF (MCON .GT. M .AND. KK .NE. MCON) THEN
- K=NACT-1
- SP=0.0
- DO 220 I=1,N
- 220 SP=SP+Z(I,K)*A(I,KK)
- TEMP=SQRT(SP*SP+ZDOTA(NACT)**2)
- ALPHA=ZDOTA(NACT)/TEMP
- BETA=SP/TEMP
- ZDOTA(NACT)=ALPHA*ZDOTA(K)
- ZDOTA(K)=TEMP
- DO 230 I=1,N
- TEMP=ALPHA*Z(I,NACT)+BETA*Z(I,K)
- Z(I,NACT)=ALPHA*Z(I,K)-BETA*Z(I,NACT)
- 230 Z(I,K)=TEMP
- IACT(NACT)=IACT(K)
- IACT(K)=KK
- TEMP=VMULTC(K)
- VMULTC(K)=VMULTC(NACT)
- VMULTC(NACT)=TEMP
- END IF
- C
- C If stage one is in progress, then set SDIRN to the direction of the next
- C change to the current vector of variables.
- C
- IF (MCON .GT. M) GOTO 320
- KK=IACT(NACT)
- TEMP=0.0
- DO 240 I=1,N
- 240 TEMP=TEMP+SDIRN(I)*A(I,KK)
- TEMP=TEMP-1.0
- TEMP=TEMP/ZDOTA(NACT)
- DO 250 I=1,N
- 250 SDIRN(I)=SDIRN(I)-TEMP*Z(I,NACT)
- GOTO 340
- C
- C Delete the constraint that has the index IACT(ICON) from the active set.
- C
- 260 IF (ICON .LT. NACT) THEN
- ISAVE=IACT(ICON)
- VSAVE=VMULTC(ICON)
- K=ICON
- 270 KP=K+1
- KK=IACT(KP)
- SP=0.0
- DO 280 I=1,N
- 280 SP=SP+Z(I,K)*A(I,KK)
- TEMP=SQRT(SP*SP+ZDOTA(KP)**2)
- ALPHA=ZDOTA(KP)/TEMP
- BETA=SP/TEMP
- ZDOTA(KP)=ALPHA*ZDOTA(K)
- ZDOTA(K)=TEMP
- DO 290 I=1,N
- TEMP=ALPHA*Z(I,KP)+BETA*Z(I,K)
- Z(I,KP)=ALPHA*Z(I,K)-BETA*Z(I,KP)
- 290 Z(I,K)=TEMP
- IACT(K)=KK
- VMULTC(K)=VMULTC(KP)
- K=KP
- IF (K .LT. NACT) GOTO 270
- IACT(K)=ISAVE
- VMULTC(K)=VSAVE
- END IF
- NACT=NACT-1
- C
- C If stage one is in progress, then set SDIRN to the direction of the next
- C change to the current vector of variables.
- C
- IF (MCON .GT. M) GOTO 320
- TEMP=0.0
- DO 300 I=1,N
- 300 TEMP=TEMP+SDIRN(I)*Z(I,NACT+1)
- DO 310 I=1,N
- 310 SDIRN(I)=SDIRN(I)-TEMP*Z(I,NACT+1)
- GO TO 340
- C
- C Pick the next search direction of stage two.
- C
- 320 TEMP=1.0/ZDOTA(NACT)
- DO 330 I=1,N
- 330 SDIRN(I)=TEMP*Z(I,NACT)
- C
- C Calculate the step to the boundary of the trust region or take the step
- C that reduces RESMAX to zero. The two statements below that include the
- C factor 1.0E-6 prevent some harmless underflows that occurred in a test
- C calculation. Further, we skip the step if it could be zero within a
- C reasonable tolerance for computer rounding errors.
- C
- 340 DD=RHO*RHO
- SD=0.0
- SS=0.0
- DO 350 I=1,N
- IF (ABS(DX(I)) .GE. 1.0E-6*RHO) DD=DD-DX(I)**2
- SD=SD+DX(I)*SDIRN(I)
- 350 SS=SS+SDIRN(I)**2
- IF (DD .LE. 0.0) GOTO 490
- TEMP=SQRT(SS*DD)
- IF (ABS(SD) .GE. 1.0E-6*TEMP) TEMP=SQRT(SS*DD+SD*SD)
- STPFUL=DD/(TEMP+SD)
- STEP=STPFUL
- IF (MCON .EQ. M) THEN
- ACCA=STEP+0.1*RESMAX
- ACCB=STEP+0.2*RESMAX
- IF (STEP .GE. ACCA .OR. ACCA .GE. ACCB) GOTO 480
- STEP=AMIN1(STEP,RESMAX)
- END IF
- C
- C Set DXNEW to the new variables if STEP is the steplength, and reduce
- C RESMAX to the corresponding maximum residual if stage one is being done.
- C Because DXNEW will be changed during the calculation of some Lagrange
- C multipliers, it will be restored to the following value later.
- C
- DO 360 I=1,N
- 360 DXNEW(I)=DX(I)+STEP*SDIRN(I)
- IF (MCON .EQ. M) THEN
- RESOLD=RESMAX
- RESMAX=0.0
- DO 380 K=1,NACT
- KK=IACT(K)
- TEMP=B(KK)
- DO 370 I=1,N
- 370 TEMP=TEMP-A(I,KK)*DXNEW(I)
- RESMAX=AMAX1(RESMAX,TEMP)
- 380 CONTINUE
- END IF
- C
- C Set VMULTD to the VMULTC vector that would occur if DX became DXNEW. A
- C device is included to force VMULTD(K)=0.0 if deviations from this value
- C can be attributed to computer rounding errors. First calculate the new
- C Lagrange multipliers.
- C
- K=NACT
- 390 ZDOTW=0.0
- ZDWABS=0.0
- DO 400 I=1,N
- TEMP=Z(I,K)*DXNEW(I)
- ZDOTW=ZDOTW+TEMP
- 400 ZDWABS=ZDWABS+ABS(TEMP)
- ACCA=ZDWABS+0.1*ABS(ZDOTW)
- ACCB=ZDWABS+0.2*ABS(ZDOTW)
- IF (ZDWABS .GE. ACCA .OR. ACCA .GE. ACCB) ZDOTW=0.0
- VMULTD(K)=ZDOTW/ZDOTA(K)
- IF (K .GE. 2) THEN
- KK=IACT(K)
- DO 410 I=1,N
- 410 DXNEW(I)=DXNEW(I)-VMULTD(K)*A(I,KK)
- K=K-1
- GOTO 390
- END IF
- IF (MCON .GT. M) VMULTD(NACT)=AMAX1(0.0,VMULTD(NACT))
- C
- C Complete VMULTC by finding the new constraint residuals.
- C
- DO 420 I=1,N
- 420 DXNEW(I)=DX(I)+STEP*SDIRN(I)
- IF (MCON .GT. NACT) THEN
- KL=NACT+1
- DO 440 K=KL,MCON
- KK=IACT(K)
- SUM=RESMAX-B(KK)
- SUMABS=RESMAX+ABS(B(KK))
- DO 430 I=1,N
- TEMP=A(I,KK)*DXNEW(I)
- SUM=SUM+TEMP
- 430 SUMABS=SUMABS+ABS(TEMP)
- ACCA=SUMABS+0.1*ABS(SUM)
- ACCB=SUMABS+0.2*ABS(SUM)
- IF (SUMABS .GE. ACCA .OR. ACCA .GE. ACCB) SUM=0.0
- 440 VMULTD(K)=SUM
- END IF
- C
- C Calculate the fraction of the step from DX to DXNEW that will be taken.
- C
- RATIO=1.0
- ICON=0
- DO 450 K=1,MCON
- IF (VMULTD(K) .LT. 0.0) THEN
- TEMP=VMULTC(K)/(VMULTC(K)-VMULTD(K))
- IF (TEMP .LT. RATIO) THEN
- RATIO=TEMP
- ICON=K
- END IF
- END IF
- 450 CONTINUE
- C
- C Update DX, VMULTC and RESMAX.
- C
- TEMP=1.0-RATIO
- DO 460 I=1,N
- 460 DX(I)=TEMP*DX(I)+RATIO*DXNEW(I)
- DO 470 K=1,MCON
- 470 VMULTC(K)=AMAX1(0.0,TEMP*VMULTC(K)+RATIO*VMULTD(K))
- IF (MCON .EQ. M) RESMAX=RESOLD+RATIO*(RESMAX-RESOLD)
- C
- C If the full step is not acceptable then begin another iteration.
- C Otherwise switch to stage two or end the calculation.
- C
- IF (ICON .GT. 0) GOTO 70
- IF (STEP .EQ. STPFUL) GOTO 500
- 480 MCON=M+1
- ICON=MCON
- IACT(MCON)=MCON
- VMULTC(MCON)=0.0
- GOTO 60
- C
- C We employ any freedom that may be available to reduce the objective
- C function before returning a DX whose length is less than RHO.
- C
- 490 IF (MCON .EQ. M) GOTO 480
- IFULL=0
- 500 RETURN
- END
- C------------------------------------------------------------------------------
- C Main program of test problems in Report DAMTP 1992/NA5.
- C------------------------------------------------------------------------------
- COMMON NPROB
- DIMENSION X(10),XOPT(10),W(3000),IACT(51)
- DO 180 NPROB=1,10
- IF (NPROB .EQ. 1) THEN
- C
- C Minimization of a simple quadratic function of two variables.
- C
- PRINT 10
- 10 FORMAT (/7X,'Output from test problem 1 (Simple quadratic)')
- N=2
- M=0
- XOPT(1)=-1.0
- XOPT(2)=0.0
- ELSE IF (NPROB .EQ. 2) THEN
- C
- C Easy two dimensional minimization in unit circle.
- C
- PRINT 20
- 20 FORMAT (/7X,'Output from test problem 2 (2D unit circle ',
- 1 'calculation)')
- N=2
- M=1
- XOPT(1)=SQRT(0.5)
- XOPT(2)=-XOPT(1)
- ELSE IF (NPROB .EQ. 3) THEN
- C
- C Easy three dimensional minimization in ellipsoid.
- C
- PRINT 30
- 30 FORMAT (/7X,'Output from test problem 3 (3D ellipsoid ',
- 1 'calculation)')
- N=3
- M=1
- XOPT(1)=1.0/SQRT(3.0)
- XOPT(2)=1.0/SQRT(6.0)
- XOPT(3)=-1.0/3.0
- ELSE IF (NPROB .EQ. 4) THEN
- C
- C Weak version of Rosenbrock's problem.
- C
- PRINT 40
- 40 FORMAT (/7X,'Output from test problem 4 (Weak Rosenbrock)')
- N=2
- M=0
- XOPT(1)=-1.0
- XOPT(2)=1.0
- ELSE IF (NPROB .EQ. 5) THEN
- C
- C Intermediate version of Rosenbrock's problem.
- C
- PRINT 50
- 50 FORMAT (/7X,'Output from test problem 5 (Intermediate ',
- 1 'Rosenbrock)')
- N=2
- M=0
- XOPT(1)=-1.0
- XOPT(2)=1.0
- ELSE IF (NPROB .EQ. 6) THEN
- C
- C This problem is taken from Fletcher's book Practical Methods of
- C Optimization and has the equation number (9.1.15).
- C
- PRINT 60
- 60 FORMAT (/7X,'Output from test problem 6 (Equation ',
- 1 '(9.1.15) in Fletcher)')
- N=2
- M=2
- XOPT(1)=SQRT(0.5)
- XOPT(2)=XOPT(1)
- ELSE IF (NPROB .EQ. 7) THEN
- C
- C This problem is taken from Fletcher's book Practical Methods of
- C Optimization and has the equation number (14.4.2).
- C
- PRINT 70
- 70 FORMAT (/7X,'Output from test problem 7 (Equation ',
- 1 '(14.4.2) in Fletcher)')
- N=3
- M=3
- XOPT(1)=0.0
- XOPT(2)=-3.0
- XOPT(3)=-3.0
- ELSE IF (NPROB .EQ. 8) THEN
- C
- C This problem is taken from page 66 of Hock and Schittkowski's book Test
- C Examples for Nonlinear Programming Codes. It is their test problem Number
- C 43, and has the name Rosen-Suzuki.
- C
- PRINT 80
- 80 FORMAT (/7X,'Output from test problem 8 (Rosen-Suzuki)')
- N=4
- M=3
- XOPT(1)=0.0
- XOPT(2)=1.0
- XOPT(3)=2.0
- XOPT(4)=-1.0
- ELSE IF (NPROB .EQ. 9) THEN
- C
- C This problem is taken from page 111 of Hock and Schittkowski's
- C book Test Examples for Nonlinear Programming Codes. It is their
- C test problem Number 100.
- C
- PRINT 90
- 90 FORMAT (/7X,'Output from test problem 9 (Hock and ',
- 1 'Schittkowski 100)')
- N=7
- M=4
- XOPT(1)=2.330499
- XOPT(2)=1.951372
- XOPT(3)=-0.4775414
- XOPT(4)=4.365726
- XOPT(5)=-0.624487
- XOPT(6)=1.038131
- XOPT(7)=1.594227
- ELSE IF (NPROB .EQ. 10) THEN
- C
- C This problem is taken from page 415 of Luenberger's book Applied
- C Nonlinear Programming. It is to maximize the area of a hexagon of
- C unit diameter.
- C
- PRINT 100
- 100 FORMAT (/7X,'Output from test problem 10 (Hexagon area)')
- N=9
- M=14
- END IF
- DO 160 ICASE=1,2
- DO 120 I=1,N
- 120 X(I)=1.0
- RHOBEG=0.5
- RHOEND=0.001
- IF (ICASE .EQ. 2) RHOEND=0.0001
- IPRINT=1
- MAXFUN=2000
- CALL COBYLA (N,M,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W,IACT)
- IF (NPROB .EQ. 10) THEN
- TEMPA=X(1)+X(3)+X(5)+X(7)
- TEMPB=X(2)+X(4)+X(6)+X(8)
- TEMPC=0.5/SQRT(TEMPA*TEMPA+TEMPB*TEMPB)
- TEMPD=TEMPC*SQRT(3.0)
- XOPT(1)=TEMPD*TEMPA+TEMPC*TEMPB
- XOPT(2)=TEMPD*TEMPB-TEMPC*TEMPA
- XOPT(3)=TEMPD*TEMPA-TEMPC*TEMPB
- XOPT(4)=TEMPD*TEMPB+TEMPC*TEMPA
- DO 130 I=1,4
- 130 XOPT(I+4)=XOPT(I)
- END IF
- TEMP=0.0
- DO 140 I=1,N
- 140 TEMP=TEMP+(X(I)-XOPT(I))**2
- PRINT 150, SQRT(TEMP)
- 150 FORMAT (/5X,'Least squares error in variables =',1PE16.6)
- 160 CONTINUE
- PRINT 170
- 170 FORMAT (2X,'----------------------------------------------',
- 1 '--------------------')
- 180 CONTINUE
- STOP
- END
- C------------------------------------------------------------------------------
- SUBROUTINE CALCFC (N,M,X,F,CON)
- COMMON NPROB
- DIMENSION X(*),CON(*)
- IF (NPROB .EQ. 1) THEN
- C
- C Test problem 1 (Simple quadratic)
- C
- F=10.0*(X(1)+1.0)**2+X(2)**2
- ELSE IF (NPROB .EQ. 2) THEN
- C
- C Test problem 2 (2D unit circle calculation)
- C
- F=X(1)*X(2)
- CON(1)=1.0-X(1)**2-X(2)**2
- ELSE IF (NPROB .EQ. 3) THEN
- C
- C Test problem 3 (3D ellipsoid calculation)
- C
- F=X(1)*X(2)*X(3)
- CON(1)=1.0-X(1)**2-2.0*X(2)**2-3.0*X(3)**2
- ELSE IF (NPROB .EQ. 4) THEN
- C
- C Test problem 4 (Weak Rosenbrock)
- C
- F=(X(1)**2-X(2))**2+(1.0+X(1))**2
- ELSE IF (NPROB .EQ. 5) THEN
- C
- C Test problem 5 (Intermediate Rosenbrock)
- C
- F=10.0*(X(1)**2-X(2))**2+(1.0+X(1))**2
- ELSE IF (NPROB .EQ. 6) THEN
- C
- C Test problem 6 (Equation (9.1.15) in Fletcher's book)
- C
- F=-X(1)-X(2)
- CON(1)=X(2)-X(1)**2
- CON(2)=1.0-X(1)**2-X(2)**2
- ELSE IF (NPROB .EQ. 7) THEN
- C
- C Test problem 7 (Equation (14.4.2) in Fletcher's book)
- C
- F=X(3)
- CON(1)=5.0*X(1)-X(2)+X(3)
- CON(2)=X(3)-X(1)**2-X(2)**2-4.0*X(2)
- CON(3)=X(3)-5.0*X(1)-X(2)
- ELSE IF (NPROB .EQ. 8) THEN
- C
- C Test problem 8 (Rosen-Suzuki)
- C
- F=X(1)**2+X(2)**2+2.0*X(3)**2+X(4)**2-5.0*X(1)-5.0*X(2)
- 1 -21.0*X(3)+7.0*X(4)
- CON(1)=8.0-X(1)**2-X(2)**2-X(3)**2-X(4)**2-X(1)+X(2)
- 1 -X(3)+X(4)
- CON(2)=10.0-X(1)**2-2.0*X(2)**2-X(3)**2-2.0*X(4)**2+X(1)+X(4)
- CON(3)=5.0-2.0*X(1)**2-X(2)**2-X(3)**2-2.0*X(1)+X(2)+X(4)
- ELSE IF (NPROB .EQ. 9) THEN
- C
- C Test problem 9 (Hock and Schittkowski 100)
- C
- F=(X(1)-10.0)**2+5.0*(X(2)-12.0)**2+X(3)**4+3.0*(X(4)-11.0)**2
- 1 +10.0*X(5)**6+7.0*X(6)**2+X(7)**4-4.0*X(6)*X(7)-10.0*X(6)
- 2 -8.0*X(7)
- CON(1)=127.0-2.0*X(1)**2-3.0*X(2)**4-X(3)-4.0*X(4)**2-5.0*X(5)
- CON(2)=282.0-7.0*X(1)-3.0*X(2)-10.0*X(3)**2-X(4)+X(5)
- CON(3)=196.0-23.0*X(1)-X(2)**2-6.0*X(6)**2+8.0*X(7)
- CON(4)=-4.0*X(1)**2-X(2)**2+3.0*X(1)*X(2)-2.0*X(3)**2-5.0*X(6)
- 1 +11.0*X(7)
- ELSE IF (NPROB .EQ. 10) THEN
- C
- C Test problem 10 (Hexagon area)
- C
- F=-0.5*(X(1)*X(4)-X(2)*X(3)+X(3)*X(9)-X(5)*X(9)+X(5)*X(8)
- 1 -X(6)*X(7))
- CON(1)=1.0-X(3)**2-X(4)**2
- CON(2)=1.0-X(9)**2
- CON(3)=1.0-X(5)**2-X(6)**2
- CON(4)=1.0-X(1)**2-(X(2)-X(9))**2
- CON(5)=1.0-(X(1)-X(5))**2-(X(2)-X(6))**2
- CON(6)=1.0-(X(1)-X(7))**2-(X(2)-X(8))**2
- CON(7)=1.0-(X(3)-X(5))**2-(X(4)-X(6))**2
- CON(8)=1.0-(X(3)-X(7))**2-(X(4)-X(8))**2
- CON(9)=1.0-X(7)**2-(X(8)-X(9))**2
- CON(10)=X(1)*X(4)-X(2)*X(3)
- CON(11)=X(3)*X(9)
- CON(12)=-X(5)*X(9)
- CON(13)=X(5)*X(8)-X(6)*X(7)
- CON(14)=X(9)
- END IF
- RETURN
- END
- -------------------------------------------------------------------------------
- Output from test problem 1 (Simple quadratic)
- Normal return from subroutine COBYLA
- NFVALS = 37 F = 1.809750E-05 MAXCV = 0.000000E+00
- X =-1.000879E+00 3.220609E-03
- Least squares error in variables = 3.338389E-03
- Normal return from subroutine COBYLA
- NFVALS = 65 F = 1.153291E-07 MAXCV = 0.000000E+00
- X =-9.999341E-01 2.682342E-04
- Least squares error in variables = 2.762020E-04
- ------------------------------------------------------------------
- Output from test problem 2 (2D unit circle calculation)
- Normal return from subroutine COBYLA
- NFVALS = 37 F =-4.999994E-01 MAXCV = 2.026558E-06
- X = 7.062163E-01 -7.079976E-01
- Least squares error in variables = 1.259601E-03
- Normal return from subroutine COBYLA
- NFVALS = 44 F =-5.000000E-01 MAXCV = 5.960464E-08
- X = 7.071500E-01 -7.070636E-01
- Least squares error in variables = 6.107080E-05
- ------------------------------------------------------------------
- Output from test problem 3 (3D ellipsoid calculation)
- Normal return from subroutine COBYLA
- NFVALS = 52 F =-7.856688E-02 MAXCV = 6.288290E-06
- X = 5.780203E-01 4.069204E-01 -3.340311E-01
- Least squares error in variables = 1.642899E-03
- Normal return from subroutine COBYLA
- NFVALS = 65 F =-7.856742E-02 MAXCV = 8.940697E-08
- X = 5.773363E-01 4.082997E-01 -3.332995E-01
- Least squares error in variables = 6.312904E-05
- ------------------------------------------------------------------
- Output from test problem 4 (Weak Rosenbrock)
- Normal return from subroutine COBYLA
- NFVALS = 100 F = 3.125441E-05 MAXCV = 0.000000E+00
- X =-9.958720E-01 9.879909E-01
- Least squares error in variables = 1.269875E-02
- Normal return from subroutine COBYLA
- NFVALS = 173 F = 6.409362E-07 MAXCV = 0.000000E+00
- X =-9.994782E-01 9.983495E-01
- Least squares error in variables = 1.730967E-03
- ------------------------------------------------------------------
- Output from test problem 5 (Intermediate Rosenbrock)
- Normal return from subroutine COBYLA
- NFVALS = 347 F = 4.008353E-03 MAXCV = 0.000000E+00
- X =-9.366965E-01 8.777190E-01
- Least squares error in variables = 1.376952E-01
- Normal return from subroutine COBYLA
- NFVALS = 698 F = 9.516375E-05 MAXCV = 0.000000E+00
- X =-9.904447E-01 9.803594E-01
- Least squares error in variables = 2.184159E-02
- ------------------------------------------------------------------
- Output from test problem 6 (Equation (9.1.15) in Fletcher)
- Normal return from subroutine COBYLA
- NFVALS = 30 F =-1.414216E+00 MAXCV = 2.950430E-06
- X = 7.071948E-01 7.070208E-01
- Least squares error in variables = 1.230355E-04
- Normal return from subroutine COBYLA
- NFVALS = 40 F =-1.414214E+00 MAXCV = 0.000000E+00
- X = 7.071791E-01 7.070344E-01
- Least squares error in variables = 1.023325E-04
- ------------------------------------------------------------------
- Output from test problem 7 (Equation (14.4.2) in Fletcher)
- Normal return from subroutine COBYLA
- NFVALS = 28 F =-2.999881E+00 MAXCV = 0.000000E+00
- X = 1.362514E-08 -2.999881E+00 -2.999881E+00
- Least squares error in variables = 1.689246E-04
- Normal return from subroutine COBYLA
- NFVALS = 32 F =-3.000046E+00 MAXCV = 4.673004E-05
- X = 1.207445E-08 -3.000000E+00 -3.000046E+00
- Least squares error in variables = 4.649224E-05
- ------------------------------------------------------------------
- Output from test problem 8 (Rosen-Suzuki)
- Normal return from subroutine COBYLA
- NFVALS = 66 F =-4.400000E+01 MAXCV = 3.099442E-06
- X =-6.020486E-04 9.995968E-01 2.000546E+00 -9.994259E-01
- Least squares error in variables = 1.073541E-03
- Normal return from subroutine COBYLA
- NFVALS = 79 F =-4.400000E+01 MAXCV = 1.251698E-06
- X =-2.246869E-04 9.996516E-01 2.000260E+00 -9.997578E-01
- Least squares error in variables = 5.460466E-04
- ------------------------------------------------------------------
- Output from test problem 9 (Hock and Schittkowski 100)
- Normal return from subroutine COBYLA
- NFVALS = 237 F = 6.806309E+02 MAXCV = 0.000000E+00
- X = 2.332463E+00 1.951341E+00 -4.587620E-01 4.364742E+00 -6.244753E-01
- 1.038812E+00 1.593632E+00
- Least squares error in variables = 1.892897E-02
- Normal return from subroutine COBYLA
- NFVALS = 248 F = 6.806310E+02 MAXCV = 1.907349E-05
- X = 2.332446E+00 1.951307E+00 -4.577461E-01 4.364753E+00 -6.241184E-01
- 1.039491E+00 1.593760E+00
- Least squares error in variables = 1.996995E-02
- ------------------------------------------------------------------
- Output from test problem 10 (Hexagon area)
- Normal return from subroutine COBYLA
- NFVALS = 150 F =-8.660254E-01 MAXCV = 1.192093E-06
- X = 6.605685E-01 7.507660E-01 -3.188329E-01 9.478114E-01 6.614124E-01
- 7.500232E-01 -3.198982E-01 9.474520E-01 -6.671554E-11
- Least squares error in variables = 1.124314E-03
- Normal return from subroutine COBYLA
- NFVALS = 173 F =-8.660254E-01 MAXCV = 3.352761E-07
- X = 6.606672E-01 7.506790E-01 -3.195507E-01 9.475691E-01 6.608437E-01
- 7.505235E-01 -3.197733E-01 9.474941E-01 -3.822812E-11
- Least squares error in variables = 2.350494E-04
- ------------------------------------------------------------------
|