@@ -77,8 +77,8 @@ SUBROUTINE OGOPTI (VARACC, NUMEQU, FINISH, DESNOR, CALVAL)
7777C WRITE (STR,'("NUMACT = ",I4)') NUMACT
7878C CALL OGWRIT (3,STR)
7979 NUMCOR = NUMACT
80- CONCOR = ACTCON
81- NUMACT = 0
80+ CONCOR = ACTCON ! ACTCON = list of indices of active constraints
81+ NUMACT = 0 ! number of active constraints
8282 CONACT(1 :DES) = 0
8383C DO COR = 1, NUMCOR
8484C CON = CONCOR(COR)
@@ -89,7 +89,8 @@ SUBROUTINE OGOPTI (VARACC, NUMEQU, FINISH, DESNOR, CALVAL)
8989C CALL OGINCL (CON)
9090C ENDDO
9191C ======================================================================
92- 1100 CONTINUE
92+ 1100 CONTINUE ! Begin optimisation iteration
93+
9394C ======================================================================
9495C VECTOR OF STEEPEST ASCENT
9596C ----------------------------------------------------------------------
@@ -103,23 +104,23 @@ SUBROUTINE OGOPTI (VARACC, NUMEQU, FINISH, DESNOR, CALVAL)
103104C REMOVE PASSIVE INEQUALITY CONSTRAINTS
104105C ----------------------------------------------------------------------
105106 DO ACT = NUMACT, 1 , - 1
106- CON = ACTCON(ACT)
107- IF (CONVAL(CON) .LE. 1D0 ) CYCLE
107+ CON = ACTCON(ACT) ! CON = index of active constraint in CONVAL
108+ IF (CONVAL(CON) .LE. 1D0 ) CYCLE
108109 NAM = CONSTR(CON)
109110 LEN = CONLEN(CON)
110111 WRITE (STR,' (I4,5X,1X,A)' ) CON, NAM(1 :LEN)
111112 CALL OGWRIT (2 ,STR)
112- CALL OGEXCL (ACT)
113- CONACT(CON) = - 1
113+ CALL OGEXCL (ACT) ! Exclude constraint ACT from active set
114+ CONACT(CON) = - 1 ! Mark CON as blocked from inclusion in active set
114115 ENDDO
115116C ----------------------------------------------------------------------
116117C INCLUDE VIOLATED INEQUALITY CONSTRAINTS AND SELECT PASSIVE ONES
117118C SELECT PASSIVE INEQUALITY CONSTRAINTS
118119C ----------------------------------------------------------------------
119120 DO CON = 1 , NUMCON
120121 IF (CONTYP(CON) .EQ. - 2 ) CYCLE
121- IF (CONACT(CON) .GT. 0 ) THEN
122- ELSEIF (CONTYP(CON) .EQ. 0 ) THEN
122+ IF (CONACT(CON) .GT. 0 ) THEN ! CON is active, CONACT(CON) is the position of CON in active set
123+ ELSEIF (CONTYP(CON) .EQ. 0 ) THEN ! CON in inactive
123124 CONACT(CON) = 0
124125 ELSEIF (CONVAL(CON) .LT. - 1D0 ) THEN
125126 CONACT(CON) = 0
@@ -131,9 +132,9 @@ SUBROUTINE OGOPTI (VARACC, NUMEQU, FINISH, DESNOR, CALVAL)
131132 ENDDO
132133C ======================================================================
133134 NNN = 1
134- 1110 CONTINUE
135+ 1110 CONTINUE ! Active set refinement loop
135136 NNN = NNN + 1
136- IF (NNN .GT. 999 ) THEN
137+ IF (NNN .GT. 999 ) THEN ! saveguard active set refinement loop to continue forever
137138 FINISH = 0
138139 WRITE (STR,* ) " NNN=" ,NNN
139140 CALL OGWRIT (2 ,STR)
@@ -142,16 +143,16 @@ SUBROUTINE OGOPTI (VARACC, NUMEQU, FINISH, DESNOR, CALVAL)
142143C ======================================================================
143144C DERIVATIVES OF MERIT W.R.T. ACTIVE CONSTRAINTS
144145C ----------------------------------------------------------------------
145- CALL OGRIGT (- CONRED(COS,1 :NUMACT), COSACT)
146- DESNOR = DSQRT(SUM (CONRED(COS,NUMACT+1 :NUMVAR)** 2 ))
146+ CALL OGRIGT (- CONRED(COS,1 :NUMACT), COSACT) ! COSACT = Lagrange multipliers???
147+ DESNOR = DSQRT(SUM (CONRED(COS,NUMACT+1 :NUMVAR)** 2 )) ! norm of cost function gradient
147148C ----------------------------------------------------------------------
148149C CONSTRAINT REMOVAL
149150C ----------------------------------------------------------------------
150151 IND = 0
151152 EXC = - 1D-12
152153 MAX = EXC
153154 DO ACT = 1 , NUMACT
154- CON = ACTCON(ACT)
155+ CON = ACTCON(ACT) ! ACT = index in active set, CON = index of constraint
155156 IF (CONTYP(CON) .EQ. 0 ) CYCLE
156157 VAL = COSACT(ACT)
157158 FAC = DOT_PRODUCT (CONRED(CON,1 :NUMVAR),
@@ -163,6 +164,7 @@ SUBROUTINE OGOPTI (VARACC, NUMEQU, FINISH, DESNOR, CALVAL)
163164 MAX = VAL
164165 IND = ACT
165166 ENDDO
167+
166168C ----------------------------------------------------------------------
167169 IF (IND .NE. 0 ) THEN
168170 CON = ACTCON(IND)
@@ -171,7 +173,7 @@ SUBROUTINE OGOPTI (VARACC, NUMEQU, FINISH, DESNOR, CALVAL)
171173 WRITE (STR,' (I4,5X,3(1X,D10.3),1X,A)' )
172174 & CON,DESNOR,MAX,VARMAX,NAM(1 :LEN)
173175 CALL OGWRIT (3 ,STR)
174- CALL OGEXCL (IND)
176+ CALL OGEXCL (IND) ! Exclude constraint IND from active set
175177 GOTO 1110
176178 ENDIF
177179C ----------------------------------------------------------------------
@@ -279,7 +281,7 @@ SUBROUTINE OGOPTI (VARACC, NUMEQU, FINISH, DESNOR, CALVAL)
279281 VAL = DOT_PRODUCT (CONRED(CON,NUMACT+1 :NUMVAR),
280282 & CONRED(COS,NUMACT+1 :NUMVAR))
281283 IF (VAL .EQ. 0D0 ) CYCLE
282- VAL = - CONVAL(CON) / VAL * DESNOR
284+ VAL = - CONVAL(CON) / VAL * DESNOR
283285 IF (VAL .LE. 0D0 ) CYCLE
284286 IF (VAL .GE. DIS) CYCLE
285287 DIS = VAL
@@ -472,6 +474,20 @@ SUBROUTINE OGOPTI (VARACC, NUMEQU, FINISH, DESNOR, CALVAL)
472474 VARDES = THT * DESPRV + BET * VARDIR
473475 DESNOR = DSQRT(SUM (VARDES** 2 ))
474476C ----------------------------------------------------------------------
477+ C SAFEGUARD: ZERO STEEPEST-ASCENT NORM => TREAT AS CONVERGENCE
478+ C ----------------------------------------------------------------------
479+ ! NOTE: otherwise, DESNOR used in denominator later on may cause
480+ ! NaN issues
481+ IF (DESNOR .EQ. 0D0 ) THEN
482+ WRITE (LOGLUP,' ("OGOPTI WARN: ZERO SECOND ORDER CORRECTION")' )
483+ WRITE (LOGLUP,' (" TREATING AS CONVERGENCE")' )
484+ FOLDIS = 0D0
485+ QUACOR = 0D0
486+ COSIMP = 0D0
487+ FINISH = 1
488+ GOTO 9999
489+ ENDIF
490+ C ----------------------------------------------------------------------
475491C MAXIMUM TRAVEL DISTANCE
476492C ----------------------------------------------------------------------
477493 DIS = VARMAX
@@ -509,6 +525,7 @@ SUBROUTINE OGOPTI (VARACC, NUMEQU, FINISH, DESNOR, CALVAL)
509525 IF (SENOPT .LE. 0 .OR. NNN .EQ. 1 ) THEN
510526 FAC = VARSTP / DESNOR
511527 VARVEC = VARREF + VARDES * VARSTP / DESNOR
528+ ! TODO: do we need to initialise CONVEC and CONDER with zeros here?
512529 CALL OGEVAL (VARVEC, CONVEC, 0 , CONDER, CALVAL, CALVAL)
513530 CONQUA = MATMUL (CONDER(1 :COS, 1 :NUMVAR),
514531 & VARDES( 1 :NUMVAR))
0 commit comments