Routine to solve constrained or unconstrained function minimization.
JUNE, 1979 VERSION
STORAGE REQUIREMENTS: PROGRAM - 7000 DECIMAL WORDS (CDC COMPUTER) ARRAYS - APPROX. 2(NDV2)+26NDV+4*NCON, WHERE N3 = NDV+2. RE-SCALE VARIABLES IF REQUIRED.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(conmin_class), | intent(inout) | :: | me | |||
| real(kind=wp), | intent(inout) | :: | x(n1) |
Vector of decision variables, |
||
| real(kind=wp), | intent(inout) | :: | vlb(n1) |
Used only if |
||
| real(kind=wp), | intent(inout) | :: | vub(n1) |
Used only if |
||
| real(kind=wp), | intent(inout) | :: | g(n2) |
Not used if |
||
| real(kind=wp), | intent(inout) | :: | scal(n1) |
Not used if |
||
| real(kind=wp), | intent(inout) | :: | df(n1) |
Analytic gradient of the objective function for the current
decision variables, |
||
| real(kind=wp), | intent(inout) | :: | a(n1,n3) |
Not used if |
||
| real(kind=wp), | intent(inout) | :: | s(n1) |
Move direction in the NDV-dimensional optimization space. |
||
| real(kind=wp), | intent(inout) | :: | g1(n2) |
Not used if |
||
| real(kind=wp), | intent(inout) | :: | g2(n2) |
Not used if |
||
| real(kind=wp), | intent(inout) | :: | b(n3,n3) |
Not used if |
||
| real(kind=wp), | intent(inout) | :: | c(n4) |
Not used in |
||
| integer, | intent(inout) | :: | isc(n2) |
Not used if |
||
| integer, | intent(inout) | :: | ic(n3) |
Identifies which constraints are active or violated. |
||
| integer, | intent(inout) | :: | ms1(n5) |
Not used if |
||
| integer, | intent(in) | :: | n1 |
|
||
| integer, | intent(in) | :: | n2 |
|
||
| integer, | intent(in) | :: | n3 |
|
||
| integer, | intent(in) | :: | n4 |
|
||
| integer, | intent(in) | :: | n5 |
|
subroutine conmin(me, x, vlb, vub, g, scal, df, a, s, g1, g2, b, c, isc, ic, ms1, n1, n2, n3, n4, n5) !! Routine to solve constrained or unconstrained function minimization. !! !! JUNE, 1979 VERSION !! !! STORAGE REQUIREMENTS: !! PROGRAM - 7000 DECIMAL WORDS (CDC COMPUTER) !! ARRAYS - APPROX. 2*(NDV**2)+26*NDV+4*NCON, WHERE N3 = NDV+2. !! RE-SCALE VARIABLES IF REQUIRED. class(conmin_class), intent(inout) :: me integer, intent(in) :: n1 !! `N1 = NDV + 2` integer, intent(in) :: n2 !! `N2 = NCON + 2*NDV` integer, intent(in) :: n3 !! `N3 = NACMX1` integer, intent(in) :: n4 !! `N4 = MAX (N3,NDV)` integer, intent(in) :: n5 !! `N5 = 2*N4` real(wp), intent(inout) :: x(n1) !! Vector of decision variables, `X(I), I = 1, NDV`. The initial !! X-vector contains the user's best estimate of the set of optimum !! design variables. real(wp), intent(inout) :: vlb(n1) !! Used only if `NSIDE/=0`. `VLB(I)` is the lower allowable value !! (lower bound) of variable `X(I)`. If one or more variables, `X(I)`, !! do not have lower bounds, the corresponding `VLB(I)` must be !! initialized to a very large negative number (say `-1.0E+10`). real(wp), intent(inout) :: vub(n1) !! Used only if `NSIDE/=0`. `VUB(I)` is the maximum allowable value !! (upper bound) of `X(I)`. If one or more variables, `X(I)`, do not !! have upper bounds, the corresponding `VUB(I)` must be initialized !! to a very large positive number (say `1.0E+10`). real(wp), intent(inout) :: g(n2) !! Not used if `NCON = NSIDE = 0`. Vector containing all constraint !! functions, `G(J), J = 1, NCON` for current decision variables, `X`. !! Calculate `G(J), J = 1, NCON` if `INFO = 2`. real(wp), intent(inout) :: scal(n1) !! Not used if `NSCAL = 0`. Vector of scaling parameters. If !! `NSCAL>0` vector SCAL need not be initialized since SCAL will !! be defined in CONMIN and its associated routines. If `NSCAL<0`, !! vector SCAL is initialized in the main program, and the scaled !! variables `X(I) = X(I)/SCAL(I)`. Efficiency of the optimization !! process can sometimes be improved if the variables are either !! normalized or are scaled in such a way that the partial deri- !! vative of the objective function, OBJ, with respect to variable !! `X(I)` is of the same order of magnitude for all `X(I)`. `SCAL(I)` !! must be greater than zero because a negative value of `SCAL(I)` !! will result in a change of sign of `X(I)` and possibly yield !! erroneous optimization results. The decision of if, and how, the !! variables should be scaled is highly problem dependent, and some !! experimentation is desirable for any given class of problems. real(wp), intent(inout) :: df(n1) !! Analytic gradient of the objective function for the current !! decision variables, `X(I)`. `DF(I)` contains the partial derivative !! of `OBJ` with respect to `X(I)`. Calculate `DF(I), I = 1, !! NDV` if `INFO = 3` or `INFO = 4` and if `NFDG = 0` or `NFDG = 2`. real(wp), intent(inout) :: a(n1, n3) !! Not used if `NCON = NSIDE = 0`. Gradients of active or violated !! constraints, for current decision variables, `X(I)`. !! `A(J,I)` contains the gradient of the Jth active or violated !! constraint, `G(J)`, with respect to the Ith decision variable, !! `X(I)` for `J = 1, NAC` and `I = 1, NDV`. Calculate if `INFO = 4` !! and `NFDG = 0`. real(wp), intent(inout) :: s(n1) !! Move direction in the NDV-dimensional optimization space. `S(I)` !! gives the rate at which variable `X(I)` changes with respect to !! `ALPHA`. real(wp), intent(inout) :: g1(n2) !! Not used if `NCON = NSIDE = NSCAL = 0`. Used for temporary !! storage of constraint values `G(J), J = 1, NCON` and decision !! variables `X(I), I = 1, NDV`. real(wp), intent(inout) :: g2(n2) !! Not used if `NCON = NSIDE = 0`. Used for temporary storage of !! constraint values `G(J), J = 1, NCON`. real(wp), intent(inout) :: b(n3, n3) !! Not used if `NCON = NSIDE = 0`. Used in determining direction !! vector S for constrained minimization problems. Array `B` may !! be used for temporary storage in external routine SUB1. real(wp), intent(inout) :: c(n4) !! Not used in `NCON = NSIDE = 0`. Used with array B in determining !! direction vector S for constrained minimization problems. Used !! for temporary storage of vector X if `NSCAL/=0`. routine SUB1. integer, intent(inout) :: isc(n2) !! Not used if `NCON = 0`. Linear constraint identification vector. !! If constraint `G(J)` is known to be a linear function of the !! decision variables, `X(I)`, `ISC(I)` should be initialized to !! `ISC(I) = 1`. If constraint `G(J)` is nonlinear `ISC(I)` is initialized !! to `ISC(I) = 0`. Identification of linear constraints may improve !! efficiency of the optimization process and is therefore desirable, !! but is not essential. If `G(J)` is not specifically known to be !! linear, set `ISC(I) = 0`. integer, intent(inout) :: ic(n3) !! Identifies which constraints are active or violated. `IC(J)` !! contains the number of the Jth active or violated constraint !! for `J = 1, NAC`. For example, if `G(10)` is the first active !! or violated constraint (`G(J)<CT, J = 1,9`), set `IC(1) = 10`. !! Calculate if `INFO = 4` and `NFDG = 0`. integer, intent(inout) :: ms1(n5) !! Not used if `NCON = NSIDE = 0`. Used with array `B` in determining !! direction vector `S` for constrained minimization problems. Array !! MS1 may be used for temporary storage in external routine SUB1. real(wp) :: alp1, alp11, alp12, c1, ct1, ctc, ff1, gi, objb, objd, & scj, si, sib, x1, x12, xid, xx integer :: i, ii, j, k, m1, m2, m3, mcn1, nci, ndv1, ndv2, nfeasct, nic, nnac if (me%nscal /= 0 .and. me%igoto /= 0) then do i = 1, me%ndv x(i) = c(i) end do end if ! CONSTANTS. ndv1 = me%ndv + 1 ndv2 = me%ndv + 2 if (me%igoto /= 0) then ! ------------------------------------------------------------------ ! CHECK FOR UNBOUNDED SOLUTION ! ------------------------------------------------------------------ ! STOP IF OBJ IS LESS THAN -1.0D+40 !if (me%obj <= -1.0e+40_wp) then ! original code if (me%obj <= -huge(1.0_real32)) then ! JW : replaced with this. note this is single precision huge (~3.0e38) ! This is so it will not overflow if using single precision real numbers. if (me%iprint > 0) write (me%iunit, 5100) go to 520 end if select case (me%igoto) case (1); go to 60 case (2); go to 210 case (3); go to 200 case (4); go to 410 case (5); go to 430 end select end if ! ------------------------------------------------------------------ ! SAVE INPUT CONTROL PARAMETERS ! ------------------------------------------------------------------ if (me%iprint > 0) write (me%iunit, 7500) if (.not. (me%linobj == 0 .or. (me%ncon > 0 .or. me%nside > 0))) then ! TOTALLY UNCONSTRAINED FUNCTION WITH LINEAR OBJECTIVE. ! SOLUTION IS UNBOUNDED. if (me%iprint > 0) write (me%iunit, 5000) me%linobj, me%ncon, me%nside return end if ! ------------------------------------------------------------------ ! DEFAULTS ! ------------------------------------------------------------------ if (me%itrm <= 0) me%itrm = 3 if (me%itmax <= 0) me%itmax = 20 ndv1 = me%ndv + 1 if (me%icndir == 0) me%icndir = ndv1 if (me%delfun <= 0.0_wp) me%delfun = 0.0001_wp me%ct = -abs(me%ct) if (me%ct >= 0.0_wp) me%ct = -0.1_wp me%ctmin = abs(me%ctmin) if (me%ctmin <= 0.0_wp) me%ctmin = 0.004_wp me%ctl = -abs(me%ctl) if (me%ctl >= 0.0_wp) me%ctl = -0.01_wp me%ctlmin = abs(me%ctlmin) if (me%ctlmin <= 0.0_wp) me%ctlmin = 0.001_wp if (me%theta <= 0.0_wp) me%theta = 1.0_wp if (me%abobj1 <= 0.0_wp) me%abobj1 = 0.1_wp if (me%alphax <= 0.0_wp) me%alphax = 0.1_wp if (me%fdch <= 0.0_wp) me%fdch = 0.01_wp if (me%fdchm <= 0.0_wp) me%fdchm = 0.01_wp ! ------------------------------------------------------------------ ! INITIALIZE INTERNAL PARAMETERS ! ------------------------------------------------------------------ me%infog = 0 me%iter = 0 me%jdir = 0 me%iobj = 0 me%kobj = 0 ndv2 = me%ndv + 2 me%kcount = 0 me%ncal(1) = 0 me%ncal(2) = 0 me%nac = 0 me%nfeas = 0 me%mscal = me%nscal ct1 = me%itrm ct1 = 1.0_wp/ct1 me%dct = (me%ctmin/abs(me%ct))**ct1 me%dctl = (me%ctlmin/abs(me%ctl))**ct1 me%phi = 5.0_wp me%abobj = me%abobj1 me%ncobj = 0 me%ctam = abs(me%ctmin) me%ctbm = abs(me%ctlmin) ! CALCULATE NUMBER OF LINEAR CONSTRAINTS, NLNC. me%nlnc = 0 if (me%ncon /= 0) then do i = 1, me%ncon if (isc(i) > 0) me%nlnc = me%nlnc + 1 end do end if ! ------------------------------------------------------------------ ! CHECK TO BE SURE THAT SIDE CONSTRAINTS ARE SATISFIED ! ------------------------------------------------------------------ if (me%nside /= 0) then do i = 1, me%ndv if (vlb(i) > vub(i)) then xx = 0.5_wp*(vlb(i) + vub(i)) x(i) = xx vlb(i) = xx vub(i) = xx if (me%iprint > 0) write (me%iunit, 6500) i end if xx = x(i) - vlb(i) if (xx < 0.0_wp) then ! LOWER BOUND VIOLATED. if (me%iprint > 0) write (me%iunit, 6600) x(i), vlb(i), i x(i) = vlb(i) else xx = vub(i) - x(i) if (xx < 0.0_wp) then if (me%iprint > 0) write (me%iunit, 6700) x(i), vub(i), i x(i) = vub(i) end if end if end do end if ! ------------------------------------------------------------------ ! INITIALIZE SCALING VECTOR, SCAL ! ------------------------------------------------------------------ if (me%nscal /= 0) then if (me%nscal >= 0) then scal(1:me%ndv) = 1.0_wp else do i = 1, me%ndv si = abs(scal(i)) if (si < me%small) si = 1.0e-5_wp scal(i) = si si = 1.0_wp/si x(i) = x(i)*si if (me%nside /= 0) then vlb(i) = vlb(i)*si vub(i) = vub(i)*si end if end do end if end if ! ------------------------------------------------------------------ ! ***** CALCULATE INITIAL FUNCTION AND CONSTRAINT VALUES ***** ! ------------------------------------------------------------------ me%info = 1 me%ncal(1) = 1 me%igoto = 1 go to 580 60 me%obj1 = me%obj if (me%dabfun <= 0.0_wp) me%dabfun = 0.001_wp*abs(me%obj) if (me%dabfun < 1.0e-10_wp) me%dabfun = 1.0e-10_wp if (me%iprint > 0) then ! ------------------------------------------------------------------ ! PRINT INITIAL DESIGN INFORMATION ! ------------------------------------------------------------------ if (me%iprint > 1) then if (me%nside == 0 .and. me%ncon == 0) write (me%iunit, 8200) if (me%nside /= 0 .or. me%ncon > 0) write (me%iunit, 7600) write (me%iunit, 7700) me%iprint, me%ndv, me%itmax, me%ncon, me%nside, me%icndir, me%nscal, & me%nfdg, me%linobj, me%itrm, n1, n2, n3, n4, n5 write (me%iunit, 7900) me%ct, me%ctmin, me%ctl, me%ctlmin, me%theta, me%phi, me%delfun, me%dabfun write (me%iunit, 7800) me%fdch, me%fdchm, me%alphax, me%abobj1 if (me%nside /= 0) then write (me%iunit, 8000) do i = 1, me%ndv, 6 m1 = min(me%ndv, i + 5) write (me%iunit, 5400) i, vlb(i:m1) end do write (me%iunit, 8100) do i = 1, me%ndv, 6 m1 = min(me%ndv, i + 5) write (me%iunit, 5400) i, vub(i:m1) end do end if if (me%nscal < 0) then write (me%iunit, 8300) write (me%iunit, 9900) scal(1:me%ndv) end if if (me%ncon /= 0) then if (me%nlnc /= 0 .and. me%nlnc /= me%ncon) then write (me%iunit, 5500) do i = 1, me%ncon, 15 m1 = min(me%ncon, i + 14) write (me%iunit, 5600) i, isc(i:m1) end do else if (me%nlnc == me%ncon) write (me%iunit, 5700) if (me%nlnc == 0) write (me%iunit, 5800) end if end if end if write (me%iunit, 9700) me%obj write (me%iunit, 9800) do i = 1, me%ndv x1 = 1.0_wp if (me%nscal /= 0) x1 = scal(i) g1(i) = x(i)*x1 end do do i = 1, me%ndv, 6 m1 = min(me%ndv, i + 5) write (me%iunit, 5400) i, g1(i:m1) end do if (associated(me%report)) call me%report(me%iter, g1(1:me%ndv), me%obj, g(1:me%ncon)) if (me%ncon /= 0) then write (me%iunit, 10000) do i = 1, me%ncon, 6 m1 = min(me%ncon, i + 5) write (me%iunit, 5400) i, g(i:m1) end do end if end if if (me%iprint > 1) write (me%iunit, 8900) ! ------------------------------------------------------------------ ! ******************** BEGIN MINIMIZATION ************************ ! ------------------------------------------------------------------ 130 me%iter = me%iter + 1 if (me%abobj1 < 0.0001_wp) me%abobj1 = 0.0001_wp if (me%abobj1 > 0.2_wp) me%abobj1 = 0.2_wp if (me%alphax > 1.0_wp) me%alphax = 1.0_wp if (me%alphax < 0.001_wp) me%alphax = 0.001_wp ! THE FOLLOWING TWO LINES OF CODE WERE COMMENTED OUT ON 3/5/81 ! NFEAS=NFEAS+1 ! IF (NFEAS>10) GO TO 810 if (me%iprint > 2) write (me%iunit, 8400) me%iter if (me%iprint > 3 .and. me%ncon > 0) write (me%iunit, 8500) me%ct, me%ctl, me%phi me%cta = abs(me%ct) if (me%ncobj /= 0) then ! ------------------------------------------------------------------ ! NO MOVE ON LAST ITERATION. DELETE CONSTRAINTS THAT ARE NO LONGER ACTIVE. ! ------------------------------------------------------------------ nnac = me%nac do i = 1, nnac if (ic(i) > me%ncon) me%nac = me%nac - 1 end do if (me%nac <= 0) go to 250 nnac = me%nac do i = 1, nnac 150 nic = ic(i) ct1 = me%ct if (isc(nic) > 0) ct1 = me%ctl if (g(nic) <= ct1) then me%nac = me%nac - 1 if (i > me%nac) go to 250 do k = i, me%nac ii = k + 1 do j = 1, ndv2 a(j, k) = a(j, ii) end do ic(k) = ic(ii) end do go to 150 end if end do go to 250 end if if (me%mscal >= me%nscal .and. me%nscal /= 0) then if (me%nscal >= 0 .or. me%kcount >= me%icndir) then me%mscal = 0 me%kcount = 0 ! ------------------------------------------------------------------ ! SCALE VARIABLES ! ------------------------------------------------------------------ do i = 1, me%ndv si = scal(i) me%xi = si*x(i) sib = si if (me%nscal > 0) si = abs(me%xi) if (si >= 1.0e-10_wp) then scal(i) = si si = 1.0_wp/si x(i) = me%xi*si if (me%nside /= 0) then vlb(i) = sib*si*vlb(i) vub(i) = sib*si*vub(i) end if end if end do if (.not. (me%iprint < 4 .or. (me%nscal < 0 .and. me%iter > 1))) then write (me%iunit, 8600) write (me%iunit, 9900) scal(1:me%ndv) end if end if end if me%mscal = me%mscal + 1 me%nac = 0 ! ------------------------------------------------------------------ ! OBTAIN GRADIENTS OF OBJECTIVE AND ACTIVE CONSTRAINTS ! ------------------------------------------------------------------ me%info = 2 me%ncal(2) = me%ncal(2) + 1 if (me%nfdg == 1) then me%igoto = 2 go to 580 end if me%jgoto = 0 200 call me%cnmn01(me%jgoto, x, df, g, isc, ic, a, g1, vub, scal, me%ncal, me%dx, me%dx1, & me%fi, me%xi, me%iii, n1, n2, n3) me%igoto = 3 if (me%jgoto > 0) go to 580 210 me%info = 1 if (me%nac >= n3) go to 520 if (me%nscal /= 0 .and. me%nfdg /= 0) then ! ------------------------------------------------------------------ ! SCALE GRADIENTS ! ------------------------------------------------------------------ ! SCALE GRADIENT OF OBJECTIVE FUNCTION. df(1:me%ndv) = df(1:me%ndv)*scal(1:me%ndv) if (me%nfdg /= 2 .and. me%nac /= 0) then ! SCALE GRADIENTS OF ACTIVE CONSTRAINTS. do j = 1, me%ndv scj = scal(j) a(j, 1:me%nac) = a(j, 1:me%nac)*scj end do end if end if 250 if (me%iprint >= 3 .and. me%ncon /= 0) then ! ------------------------------------------------------------------ ! PRINT ! ------------------------------------------------------------------ ! PRINT ACTIVE AND VIOLATED CONSTRAINT NUMBERS. m1 = 0 m2 = n3 if (me%nac /= 0) then do i = 1, me%nac j = ic(i) if (j <= me%ncon) then gi = g(j) c1 = me%ctam if (isc(j) > 0) c1 = me%ctbm gi = gi - c1 if (gi <= 0.0_wp) then ! ACTIVE CONSTRAINT. m1 = m1 + 1 ms1(m1) = j else m2 = m2 + 1 ! VIOLATED CONSTRAINT. ms1(m2) = j end if end if end do end if m3 = m2 - n3 write (me%iunit, 5900) m1 if (m1 /= 0) then write (me%iunit, 6000) write (me%iunit, 10100) ms1(1:m1) end if write (me%iunit, 6100) m3 if (m3 /= 0) then write (me%iunit, 6000) m3 = n3 + 1 write (me%iunit, 10100) ms1(m3:m2) end if end if ! ------------------------------------------------------------------ ! CALCULATE GRADIENTS OF ACTIVE SIDE CONSTRAINTS ! ------------------------------------------------------------------ if (me%nside /= 0) then mcn1 = me%ncon m1 = 0 do i = 1, me%ndv ! LOWER BOUND. me%xi = x(i) xid = vlb(i) x12 = abs(xid) if (x12 < 1.0_wp) x12 = 1.0_wp gi = (xid - me%xi)/x12 if (gi >= -1.0e-6_wp) then m1 = m1 + 1 ms1(m1) = -i me%nac = me%nac + 1 if (me%nac >= n3) go to 520 mcn1 = mcn1 + 1 do j = 1, me%ndv a(j, me%nac) = 0.0_wp end do a(i, me%nac) = -1.0_wp ic(me%nac) = mcn1 g(mcn1) = gi isc(mcn1) = 1 end if ! UPPER BOUND. xid = vub(i) x12 = abs(xid) if (x12 < 1.0_wp) x12 = 1.0_wp gi = (me%xi - xid)/x12 if (gi >= -1.0e-6_wp) then m1 = m1 + 1 ms1(m1) = i me%nac = me%nac + 1 if (me%nac >= n3) go to 520 mcn1 = mcn1 + 1 do j = 1, me%ndv a(j, me%nac) = 0.0_wp end do a(i, me%nac) = 1.0_wp ic(me%nac) = mcn1 g(mcn1) = gi isc(mcn1) = 1 end if end do ! ------------------------------------------------------------------ ! PRINT ! ------------------------------------------------------------------ ! PRINT ACTIVE SIDE CONSTRAINT NUMBERS. if (me%iprint >= 3) then write (me%iunit, 6200) m1 if (m1 /= 0) then write (me%iunit, 6300) write (me%iunit, 10100) (ms1(j), j=1, m1) end if end if end if ! PRINT GRADIENTS OF ACTIVE AND VIOLATED CONSTRAINTS. if (me%iprint >= 4) then write (me%iunit, 8700) do i = 1, me%ndv, 6 m1 = min(me%ndv, i + 5) write (me%iunit, 5400) i, df(i:m1) end do if (me%nac /= 0) then write (me%iunit, 8800) do i = 1, me%nac m1 = ic(i) m2 = m1 - me%ncon m3 = 0 if (m2 > 0) m3 = abs(ms1(m2)) if (m2 <= 0) write (me%iunit, 5200) m1 if (m2 > 0) write (me%iunit, 5300) m3 do k = 1, me%ndv, 6 m1 = min(me%ndv, k + 5) write (me%iunit, 5400) k, (a(j, i), j=k, m1) end do write (me%iunit, 8900) end do end if end if ! ------------------------------------------------------------------ ! ****************** DETERMINE SEARCH DIRECTION ******************* ! ------------------------------------------------------------------ me%alp = me%large if (me%nac > 0) go to 340 ! ------------------------------------------------------------------ ! UNCONSTRAINED FUNCTION ! ------------------------------------------------------------------ ! FIND DIRECTION OF STEEPEST DESCENT OR CONJUGATE DIRECTION. ! ! S. N. 575 ADDED ON 2/25/81 ! 330 me%nvc = 0 me%nfeas = 0 me%kcount = me%kcount + 1 ! IF KCOUNT>ICNDIR RESTART CONJUGATE DIRECTION ALGORITHM. if (me%kcount > me%icndir .or. me%iobj == 2) me%kcount = 1 if (me%kcount == 1) me%jdir = 0 ! IF JDIR = 0 FIND DIRECTION OF STEEPEST DESCENT. call me%cnmn02(me%jdir, me%slope, me%dftdf1, df, s) go to 380 ! ------------------------------------------------------------------ ! CONSTRAINED FUNCTION ! ------------------------------------------------------------------ ! FIND USABLE-FEASIBLE DIRECTION. 340 me%kcount = 0 me%jdir = 0 me%phi = 10.0_wp*me%phi if (me%phi > 1000.0_wp) me%phi = 1000.0_wp ! THE FOLLOWING LINE OF CODE WAS COMMENTED OUT ON 3/5/81 ! ! IF (NFEAS==1) PHI=5. ! CALCULATE DIRECTION, S. call me%cnmn05(g, df, a, s, b, c, me%slope, me%phi, isc, ic, ms1, me%nvc, n1, n2, n3, n4, n5) ! THE FOLLOWING LINE WAS ADDED ON 2/25/81 if (me%nac == 0) go to 330 ! THE FOLLOWING FIVE LINES WERE COMMENTED OUT ON 3/5/81 ! REASON : THEY WERE NOT IN G. VANDERPLAATS LISTING ! ! IF THIS DESIGN IS FEASIBLE AND LAST ITERATION WAS INFEASIBLE, ! SET ABOBJ1=.05 (5 PERCENT). ! IF (NVC==0 .AND. NFEAS>1) ABOBJ1=.05 ! IF (NVC==0) NFEAS=0 if (me%iprint >= 3) then write (me%iunit, 9000) do i = 1, me%nac, 6 m1 = min(me%nac, i + 5) write (me%iunit, 5400) i, (a(ndv1, j), j=i, m1) end do write (me%iunit, 7400) s(ndv1) end if ! ------------------------------------------------------------------ ! ****************** ONE-DIMENSIONAL SEARCH ************************ ! ------------------------------------------------------------------ if (s(ndv1) < 1.0e-6_wp .and. me%nvc == 0) go to 450 ! ------------------------------------------------------------------ ! FIND ALPHA TO OBTAIN A FEASIBLE DESIGN ! ------------------------------------------------------------------ if (me%nvc /= 0) then me%alp = -1.0_wp do i = 1, me%nac nci = ic(i) c1 = g(nci) ctc = me%ctam if (isc(nci) > 0) ctc = me%ctbm if (c1 > ctc) then alp1 = dot_product(s(1:me%ndv), a(1:me%ndv, i)) alp1 = alp1*a(ndv2, i) if (abs(alp1) >= me%small) then alp1 = -c1/alp1 if (alp1 > me%alp) me%alp = alp1 end if end if end do end if ! ------------------------------------------------------------------ ! LIMIT CHANCE TO ABOBJ1*OBJ ! ------------------------------------------------------------------ 380 alp1 = me%large si = abs(me%obj) if (si < 0.01_wp) si = 0.01_wp if (abs(me%slope) > me%small) alp1 = me%abobj1*si/me%slope alp1 = abs(alp1) if (me%nvc > 0) alp1 = 10.0_wp*alp1 if (alp1 < me%alp) me%alp = alp1 ! ------------------------------------------------------------------ ! LIMIT CHANGE IN VARIABLE TO ALPHAX ! ------------------------------------------------------------------ alp11 = me%large do i = 1, me%ndv si = abs(s(i)) me%xi = abs(x(i)) if (si >= 1.0e-10_wp .and. me%xi >= 0.1_wp) then alp1 = me%alphax*me%xi/si if (alp1 < alp11) alp11 = alp1 end if end do if (me%nvc > 0) alp11 = 10.0_wp*alp11 if (alp11 < me%alp) me%alp = alp11 if (me%alp > me%large) me%alp = me%large if (me%alp <= me%small) me%alp = me%small if (me%iprint >= 3) then write (me%iunit, 9100) do i = 1, me%ndv, 6 m1 = min(me%ndv, i + 5) write (me%iunit, 5400) i, (s(j), j=i, m1) end do write (me%iunit, 6400) me%slope, me%alp end if if (me%ncon > 0 .or. me%nside > 0) go to 420 ! ------------------------------------------------------------------ ! DO ONE-DIMENSIONAL SEARCH FOR UNCONSTRAINED FUNCTION ! ------------------------------------------------------------------ me%jgoto = 0 410 call me%cnmn03(x, s, me%slope, me%alp, me%fff, me%a1, me%a2, me%a3, me%a4, & me%f1, me%f2, me%f3, me%f4, me%app, me%ncal, me%kount, me%jgoto) me%igoto = 4 if (me%jgoto > 0) go to 580 me%jdir = 1 ! PROCEED TO CONVERGENCE CHECK. go to 450 ! ------------------------------------------------------------------ ! SOLVE ONE-DIMENSIONAL SEARCH PROBLEM FOR CONSTRAINED FUNCTION ! ------------------------------------------------------------------ 420 me%jgoto = 0 430 call me%cnmn06(x, vlb, vub, g, scal, df, s, g1, g2, me%ctam, me%ctbm, me%slope, me%alp, me%a2, me%a3, me%a4, & me%f1, me%f2, me%f3, me%cv1, me%cv2, me%cv3, me%cv4, me%alpca, me%alpfes, me%alpln, me%alpmin, me%alpnc, & me%alpsav, me%alpsid, me%alptot, isc, me%ncal, me%nvc, me%icount, me%igood1, & me%igood2, me%igood3, me%igood4, me%ibest, me%iii, me%nlnc, me%jgoto) me%igoto = 5 if (me%jgoto > 0) go to 580 if (me%nac == 0) me%jdir = 1 ! ------------------------------------------------------------------ ! ******************* UPDATE ALPHAX ************************** ! ------------------------------------------------------------------ 450 if (me%alp > 1.0e+19_wp) me%alp = 0.0_wp ! UPDATE ALPHAX TO BE AVERAGE OF MAXIMUM CHANGE IN X(I) AND ALHPAX. alp11 = 0.0_wp do i = 1, me%ndv si = abs(s(i)) me%xi = abs(x(i)) if (me%xi >= 1.0e-10_wp) then alp1 = me%alp*si/me%xi if (alp1 > alp11) alp11 = alp1 end if end do alp11 = 0.5_wp*(alp11 + me%alphax) alp12 = 5.0_wp*me%alphax if (alp11 > alp12) alp11 = alp12 me%alphax = alp11 me%ncobj = me%ncobj + 1 ! ABSOLUTE CHANGE IN OBJECTIVE. objd = me%obj1 - me%obj objb = abs(objd) if (objb < 1.0e-10_wp) objb = 0.0_wp if (me%nac == 0 .or. objb > 0.0_wp) me%ncobj = 0 if (me%ncobj > 1) me%ncobj = 0 ! ------------------------------------------------------------------ ! PRINT ! ------------------------------------------------------------------ ! PRINT MOVE PARAMETER, NEW X-VECTOR AND CONSTRAINTS. if (me%iprint >= 3) then write (me%iunit, 9200) me%alp end if if (me%iprint >= 2) then if (objb <= 0.0_wp) then if (me%iprint == 2) write (me%iunit, 9300) me%iter, me%obj if (me%iprint > 2) write (me%iunit, 9400) me%obj else if (me%iprint /= 2) then write (me%iunit, 9500) me%obj else write (me%iunit, 9600) me%iter, me%obj end if end if write (me%iunit, 9800) do i = 1, me%ndv ff1 = 1.0_wp if (me%nscal /= 0) ff1 = scal(i) g1(i) = ff1*x(i) end do do i = 1, me%ndv, 6 m1 = min(me%ndv, i + 5) write (me%iunit, 5400) i, (g1(j), j=i, m1) end do if (associated(me%report)) call me%report(me%iter, g1(1:me%ndv), me%obj, g(1:me%ncon)) if (me%ncon /= 0) then write (me%iunit, 10000) do i = 1, me%ncon, 6 m1 = min(me%ncon, i + 5) write (me%iunit, 5400) i, (g(j), j=i, m1) end do end if end if ! THE FOLLOWING CODE WAS ADDED ON 3/5/81 ! ! IT HAD NOT BEEN REPORTED AS A FIX TO MAOB ! BUT WAS SENT TO JEFF STROUD A YEAR AGO ! SEE OTHER COMMENTS IN CONMIN SUBROUTINE FOR DELETIONS OF CODE ! ON 3/5/81 PERTAINING TO THIS FIX ! CHECK FEASIBILITY if (me%ncon > 0) then nfeasct = 10 ! added by slp 11/17/94 do i = 1, me%ncon c1 = me%ctam if (isc(i) > 0) c1 = me%ctbm if (g(i) > c1) then me%nfeas = me%nfeas + 1 go to 510 end if end do if (me%nfeas > 0) me%abobj1 = 0.05_wp !cc me%nfeas = 0 me%phi = 5.0_wp 510 if (me%nfeas >= nfeasct) go to 520 end if ! END OF INSERTED FIX ! ------------------------------------------------------------------ ! CHECK CONVERGENCE ! ------------------------------------------------------------------ ! STOP IF ITER EQUALS ITMAX. if (me%iter < me%itmax) then ! ------------------------------------------------------------------ ! ABSOLUTE CHANGE IN OBJECTIVE ! ------------------------------------------------------------------ objb = abs(objd) me%kobj = me%kobj + 1 if (objb >= me%dabfun .or. me%nfeas > 0) me%kobj = 0 ! ------------------------------------------------------------------ ! RELATIVE CHANGE IN OBJECTIVE ! ------------------------------------------------------------------ if (abs(me%obj1) > 1.0e-10_wp) objd = objd/abs(me%obj1) me%abobj1 = .5*(abs(me%abobj) + abs(objd)) me%abobj = abs(objd) me%iobj = me%iobj + 1 if (me%nvc > 0 .or. objd >= me%delfun) me%iobj = 0 if (me%iobj < me%itrm .and. me%kobj < me%itrm) then me%obj1 = me%obj ! ------------------------------------------------------------------ ! REDUCE CT IF OBJECTIVE FUNCTION IS CHANGING SLOWLY ! ------------------------------------------------------------------ if (me%iobj < 1 .or. me%nac == 0) go to 130 me%ct = me%dct*me%ct me%ctl = me%ctl*me%dctl if (abs(me%ct) < me%ctmin) me%ct = -me%ctmin if (abs(me%ctl) < me%ctlmin) me%ctl = -me%ctlmin go to 130 end if end if 520 if (me%nac >= n3 .and. me%iprint > 0) write (me%iunit, 10200) ! ------------------------------------------------------------------ ! **************** FINAL FUNCTION INFORMATION ******************** ! ------------------------------------------------------------------ if (me%nscal /= 0) then ! UN-SCALE THE DESIGN VARIABLES. do i = 1, me%ndv me%xi = scal(i) if (me%nside /= 0) then vlb(i) = me%xi*vlb(i) vub(i) = me%xi*vub(i) end if x(i) = me%xi*x(i) end do end if ! ------------------------------------------------------------------ ! PRINT FINAL RESULTS ! ------------------------------------------------------------------ if (me%iprint /= 0 .and. me%nac < n3) then write (me%iunit, 10300) write (me%iunit, 9500) me%obj write (me%iunit, 9800) do i = 1, me%ndv, 6 m1 = min(me%ndv, i + 5) write (me%iunit, 5400) i, (x(j), j=i, m1) end do if (associated(me%report)) call me%report(me%iter, x(1:me%ndv), me%obj, g(1:me%ncon)) if (me%ncon /= 0) then write (me%iunit, 10000) do i = 1, me%ncon, 6 m1 = min(me%ncon, i + 5) write (me%iunit, 5400) i, (g(j), j=i, m1) end do ! DETERMINE WHICH CONSTRAINTS ARE ACTIVE AND PRINT. me%nac = 0 me%nvc = 0 do i = 1, me%ncon me%cta = me%ctam if (isc(i) > 0) me%cta = me%ctbm gi = g(i) if (gi <= me%cta) then if (gi < me%ct .and. isc(i) == 0) cycle if (gi < me%ctl .and. isc(i) > 0) cycle me%nac = me%nac + 1 ic(me%nac) = i else me%nvc = me%nvc + 1 ms1(me%nvc) = i end if end do write (me%iunit, 5900) me%nac if (me%nac /= 0) then write (me%iunit, 6000) write (me%iunit, 10100) ic(1:me%nac) end if write (me%iunit, 6100) me%nvc if (me%nvc /= 0) then write (me%iunit, 6000) write (me%iunit, 10100) ms1(1:me%nvc) end if end if if (me%nside /= 0) then ! DETERMINE WHICH SIDE CONSTRAINTS ARE ACTIVE AND PRINT. me%nac = 0 do i = 1, me%ndv me%xi = x(i) xid = vlb(i) x12 = abs(xid) if (x12 < 1.0_wp) x12 = 1.0_wp gi = (xid - me%xi)/x12 if (gi >= -1.0e-6_wp) then me%nac = me%nac + 1 ms1(me%nac) = -i end if xid = vub(i) x12 = abs(xid) if (x12 < 1.0_wp) x12 = 1.0_wp gi = (me%xi - xid)/x12 if (gi >= -1.0e-6_wp) then me%nac = me%nac + 1 ms1(me%nac) = i end if end do write (me%iunit, 6200) me%nac if (me%nac /= 0) then write (me%iunit, 6300) write (me%iunit, 10100) (ms1(j), j=1, me%nac) end if end if write (me%iunit, 6800) if (me%iter >= me%itmax) write (me%iunit, 6900) if (me%nfeas >= nfeasct) write (me%iunit, 7000) if (me%iobj >= me%itrm) write (me%iunit, 7100) me%itrm if (me%kobj >= me%itrm) write (me%iunit, 7200) me%itrm write (me%iunit, 7300) me%iter write (me%iunit, 10400) me%ncal(1) if (me%ncon > 0) write (me%iunit, 10500) me%ncal(1) if (me%nfdg /= 0) write (me%iunit, 10600) me%ncal(2) if (me%ncon > 0 .and. me%nfdg == 1) write (me%iunit, 10700) me%ncal(2) end if me%igoto = 0 580 if (me%nscal == 0 .or. me%igoto == 0) return ! UN-SCALE VARIABLES. do i = 1, me%ndv c(i) = x(i) x(i) = x(i)*scal(i) end do return ! ------------------------------------------------------------------ ! FORMATS ! ------------------------------------------------------------------ 5000 format(//t6, & 'A COMPLETELY UNCONSTRAINED FUNCTION WITH A LINEAR OBJECTIVE IS SPECIFIED'// & t11, 'LINOBJ =', i5/t11, 'NCON =', i5/ & t11, 'NSIDE =', i5//t6, 'CONTROL RETURNED TO CALLING PROGRAM') 5100 format(//t6, & 'CONMIN HAS ACHIEVED A SOLUTION OF OBJ LESS THAN -1.0E+40'/ & t6, 'SOLUTION APPEARS TO BE UNBOUNDED'/t6, 'OPTIMIZATION IS TERMINATED') 5200 format(t6, 'CONSTRAINT NUMBER', i5) 5300 format(t6, 'SIDE CONSTRAINT ON VARIABLE', i5) 5400 format(t4, i5, ') ', 6e13.5) 5500 format(/t6, 'LINEAR CONSTRAINT IDENTIFIERS (ISC)'/ & t6, 'NON-ZERO INDICATES LINEAR CONSTRAINT') 5600 format(t4, i5, ') ', 15i5) 5700 format(/t6, 'ALL CONSTRAINTS ARE LINEAR') 5800 format(/t6, 'ALL CONSTRAINTS ARE NON-LINEAR') 5900 format(/t6, 'THERE ARE', i5, ' ACTIVE CONSTRAINTS') 6000 format(t6, 'CONSTRAINT NUMBERS ARE') 6100 format(/t6, 'THERE ARE', i5, ' VIOLATED CONSTRAINTS') 6200 format(/t6, 'THERE ARE', i5, ' ACTIVE SIDE CONSTRAINTS') 6300 format(t6, 'DECISION VARIABLES AT LOWER OR UPPER BOUNDS', & ' (MINUS INDICATES LOWER BOUND)') 6400 format(/t6, 'ONE-DIMENSIONAL SEARCH'/t6, 'INITIAL SLOPE =', e12.4, & ' PROPOSED ALPHA =', e12.4) 6500 format(//t6, '* * CONMIN DETECTS VLB(I)>VUB(I)'/ & t6, 'FIX IS SET X(I)=VLB(I)=VUB(I) = .5*(VLB(I)+VUB(I) FOR I =', i5) 6600 format(//t6, '* * CONMIN DETECTS INITIAL X(I)<VLB(I)'/t6, & 'X(I) =', e12.4, ' VLB(I) =', e12.4/t6, & 'X(I) IS SET EQUAL TO VLB(I) FOR I =', i5) 6700 format(//t6, '* * CONMIN DETECTS INITIAL X(I)>VUB(I)'/t6, & 'X(I) =', e12.4, ' VUB(I) =', e12.4/t6, & 'X(I) IS SET EQUAL TO VUB(I) FOR I =', i5) 6800 format(/t6, 'TERMINATION CRITERION') 6900 format(t11, 'ITER EQUALS ITMAX') 7000 format(t11, & 'NFEASCT CONSECUTIVE ITERATIONS FAILED TO PRODUCE A FEASIBLE DESIGN') 7100 format(t11, 'ABS(1-OBJ(I-1)/OBJ(I)) LESS THAN DELFUN FOR', i3, & ' ITERATIONS') 7200 format(t11, 'ABS(OBJ(I)-OBJ(I-1)) LESS THAN DABFUN FOR', i3, & ' ITERATIONS') 7300 format(/t6, 'NUMBER OF ITERATIONS =', i5) 7400 format(/t6, 'CONSTRAINT PARAMETER, BETA =', e14.5) 7500 format(//t13, '* * * * * * * * * * * * * * * * * * * * * * * * * * * '/t13,& '*', t65, '*'/t13, '*', t34, & 'C O N M I N', t65, '*'/t13, '*', t65, '*'/t13, '*', t29, & ' FORTRAN PROGRAM FOR ', t65, '*'/t13, '*', t65, '*'/t13, '*', t23, & 'CONSTRAINED FUNCTION MINIMIZATION', t65, '*'/t13, '*', t65, '*'/ & t13, '* * * * * * * * * * * * * * * * * * * * * * * * * * * ') 7600 format(//t6, 'CONSTRAINED FUNCTION MINIMIZATION'// & t6, 'CONTROL PARAMETERS') 7700 format(/t6, 'IPRINT NDV ITMAX NCON NSIDE ICNDIR NSCAL NFDG'/ & 8i8//t6, 'LINOBJ ITRM N1 N2 N3 N4 N5'/8i8) 7800 format(/t10, 'FDCH', t26, 'FDCHM', t42, 'ALPHAX', t58, 'ABOBJ1'/ & ' ', *(' ', e14.5)) 7900 format(/t10, 'CT', t26, 'CTMIN', t42, 'CTL', t58, 'CTLMIN'/ & ' ', *(' ', e14.5)// & t10, 'THETA', t26, 'PHI', t42, 'DELFUN', t58, 'DABFUN'/ & ' ', *(' ', e14.5)) 8000 format(/t6, 'LOWER BOUNDS ON DECISION VARIABLES (VLB)') 8100 format(/t6, 'UPPER BOUNDS ON DECISION VARIABLES (VUB)') 8200 format(//t6, 'UNCONSTRAINED FUNCTION MINIMIZATION'//t6, & 'CONTROL PARAMETERS') 8300 format(/t6, 'SCALING VECTOR (SCAL)') 8400 format(//t6, 'BEGIN ITERATION NUMBER', i5) 8500 format(/t6, 'CT =', e14.5, ' CTL =', e14.5, ' PHI =', e14.5) 8600 format(/t6, 'NEW SCALING VECTOR (SCAL)') 8700 format(/t6, 'GRADIENT OF OBJ') 8800 format(/t6, 'GRADIENTS OF ACTIVE AND VIOLATED CONSTRAINTS') 8900 format(' ') 9000 format(/t6, 'PUSH-OFF FACTORS, (THETA(I), I=1,NAC)') 9100 format(/t6, 'SEARCH DIRECTION (S-VECTOR)') 9200 format(/t6, 'CALCULATED ALPHA =', e14.5) 9300 format(//t6, 'ITER =', i5, ' OBJ =', e14.5, ' NO CHANGE IN OBJ') 9400 format(/t6, 'OBJ =', e15.6, ' NO CHANGE ON OBJ') 9500 format(/t6, 'OBJ =', e15.6) 9600 format(//t6, 'ITER =', i5, ' OBJ =', e14.5) 9700 format(//t6, 'INITIAL FUNCTION INFORMATION'//t6, 'OBJ =', e15.6) 9800 format(/t6, 'DECISION VARIABLES (X-VECTOR)') 9900 format(t4, 7e13.4) 10000 format(/t6, 'CONSTRAINT VALUES (G-VECTOR)') 10100 format(t6, 15i5) 10200 format(/t6, 'THE NUMBER OF ACTIVE AND VIOLATED CONSTRAINTS EXCEEDS N3-1.'/ & t6, 'DIMENSIONED SIZE OF MATRICES A AND B AND VECTOR IC IS INSUFFICIENT'/ & t6, 'OPTIMIZATION TERMINATED AND CONTROL RETURNED TO MAIN PROGRAM.') 10300 format(///' FINAL OPTIMIZATION INFORMATION') 10400 format(/t6, 'OBJECTIVE FUNCTION WAS EVALUATED ', i5, ' TIMES') 10500 format(/t6, 'CONSTRAINT FUNCTIONS WERE EVALUATED', i10, ' TIMES') 10600 format(/t6, 'GRADIENT OF OBJECTIVE WAS CALCULATED', i9, ' TIMES') 10700 format(/t6, 'GRADIENTS OF CONSTRAINTS WERE CALCULATED', i5, ' TIMES') end subroutine conmin