C ****************************************************************** C ****************************************************************** subroutine auglag(n,x,l,u,m,lambda,equatn,linear,epsfeas,epsopt, +f,c,cnorm,snorm,nl,nlpsupn,fu,cnormu,outiter,totiter,nwcalls, +nwtotit,msqcalls,msqtotit,alinfo,inform,nbv,niv,variableTypes, +minlpinfo) implicit none C SCALAR ARGUMENTS integer alinfo,inform,m,msqcalls,msqtotit,n,nwcalls,nwtotit, + outiter,totiter,nbv,niv,minlpinfo double precision cnorm,cnormu,epsfeas,epsopt,f,fu,nlpsupn,snorm C ARRAY ARGUMENTS logical equatn(m),linear(m) double precision c(m),l(n),lambda(m),nl(n),rho(m),u(n),x(n) integer variableTypes(n) C Solves the nonlinear programming problem C C min f(x) C C subject to C C c_j(x) = 0, j in E, C c_j(x) <= 0, j in I, C l <= x <= u, C C where E is the set of indices of the equality constraints, I is C the set of indices of the inequality constraints, and there are C n variables and m constraints, using the method of multipliers C described in C C R. Andreani, E. G. Birgin, J. M. Martínez and M. L. Schuverdt, C "On Augmented Lagrangian methods with general lower-level C constraints", SIAM Journal on Optimization 18, pp. 1286-1309, 2007. C alinfo: C C 0: Feasibility, optimality and complementarity satisfied C 1: Maximum number of algencan iterations reached C 2: Infeasible problem? C 3: Exaggerated optimality tolerance required? #include "dim.par" #include "machconst.com" #include "algconst.par" #include "outtyp.com" #include "algparam.com" #include "scaling.com" C DATA BLOCKS character genstp(0:9) data genstp(0) /'C'/ data genstp(1) /'M'/ data genstp(2) /' '/ data genstp(3) /'F'/ data genstp(4) /'G'/ data genstp(5) /'P'/ data genstp(6) /'U'/ data genstp(7) /'S'/ data genstp(8) /' '/ data genstp(9) /' '/ C LOCAL SCALARS character innstp logical scaletmp integer geninfo,i,iter,maxit,msqiter,nwinfo,nwiter,outitnp, + rhoincr,xind double precision al,cnorma,epsfeas12,epsfeas14,epsopt12,epsopt14, + epsopk,fa,nlpi,nlpsupna,rsupn,snormprev,xsupn C LOCAL ARRAYS double precision lambar(mmax),lambdaa(mmax),sigma(mmax), + sigpre(mmax),xa(nmax) C FUNCTION logical satisfyIntegerConstraints C ================================================================== C Print initial information C ================================================================== if ( iprintout .ge. 1 ) then write(* ,1000) n,m write(10,1000) n,m if ( iprintout .ge. 4 .and. nprint .ne. 0 ) then write(* ,1010) nprint,(l(i),i=1,nprint) write(* ,1020) nprint,(u(i),i=1,nprint) write(10,1010) nprint,(l(i),i=1,nprint) write(10,1020) nprint,(u(i),i=1,nprint) end if end if C ================================================================== C Initialization C ================================================================== C Counters outiter = 0 outitnp = 0 totiter = 0 nwcalls = 0 nwtotit = 0 msqcalls = 0 msqtotit = 0 C Constants epsopt12 = sqrt( epsopt ) epsopt14 = sqrt( epsopt12 ) epsfeas12 = sqrt( epsfeas ) epsfeas14 = sqrt( epsfeas12 ) C Project initial point do i = 1,n x(i) = max( l(i), min( x(i), u(i) ) ) end do C Compute current point norm xsupn = 0.0d0 do i = 1,n if ( abs( x(i) ) .ge. xsupn ) then xsupn = abs( x(i) ) xind = i end if end do C Compute objective function and constraints call sevalobjc(n,x,f,m,c,inform) if ( inform .lt. 0 ) return C Set penalty parameters and compute maximum rho if ( rhoauto ) then call comprhoini(f,m,c,equatn,rho) end if rsupn = 0.0d0 do i = 1,m rsupn = max( rsupn, rho(i) ) end do C Compute safeguarded Lagrange multipliers do i = 1,m lambar(i) = max( lammin, min( lambda(i), lammax ) ) end do C Compute complementarity and feasibility violations snormprev = bignum cnorm = 0.0d0 snorm = 0.0d0 do i = 1,m if ( equatn(i) ) then cnorm = max( cnorm, abs( c(i) ) ) sigma(i) = c(i) else cnorm = max( cnorm, c(i) ) sigma(i) = max( c(i), - lambar(i) / rho(i) ) end if snorm = max( snorm, abs( sigma(i) ) ) end do C Compute unscaled objective function and constraints norm if ( scale ) then fu = f * sf cnormu = 0.0d0 do i = 1,m if ( equatn(i) ) then cnormu = max( cnormu, abs( sc(i) * c(i) ) ) else cnormu = max( cnormu, sc(i) * c(i) ) end if end do else fu = f cnormu = cnorm end if C Compute continuous projected Lagrangian gradient norm call sevalnl(n,x,m,lambda,equatn,linear,nl,inform) if ( inform .lt. 0 ) return nlpsupn = 0.0d0 do i = 1,n if( variableTypes(i) .ne. CONTINUOUS ) cycle nlpi = x(i) - nl(i) if ( l(i) .le. nlpi .and. nlpi .le. u(i) ) then nlpi = - nl(i) else nlpi = max( l(i), min( nlpi, u(i) ) ) - x(i) end if nlpsupn = max( nlpsupn, abs( nlpi ) ) end do C ================================================================== C Main loop C ================================================================== 100 continue C ================================================================== C Print information of this iteration C ================================================================== if ( iprintout .eq. 1 ) then if ( outiter .gt. 0 ) then innstp = genstp(geninfo) else innstp = ' ' end if if ( .not. scale ) then if ( mod(outiter,10) .eq. 0 ) then write(* ,1030) write(10,1030) end if write(* ,1040) outiter,f,cnorm,snorm,nlpsupn,xsupn,rsupn, + totiter,innstp,nwcalls,nwtotit write(10,1040) outiter,f,cnorm,snorm,nlpsupn,xsupn,rsupn, + totiter,innstp,nwcalls,nwtotit else if ( mod(outiter,10) .eq. 0 ) then write(* ,1050) write(10,1050) end if write(* ,1060) outiter,fu,cnormu,f,cnorm,snorm,nlpsupn, + rsupn,totiter,innstp,nwcalls,nwtotit write(10,1060) outiter,fu,cnormu,f,cnorm,snorm,nlpsupn, + rsupn,totiter,innstp,nwcalls,nwtotit end if else if ( iprintout .ge. 2 ) then write(* ,1070) outiter write(10,1070) outiter if ( iprintout .ge. 3 ) then write(* ,1080) totiter,nwcalls,nwtotit write(10,1080) totiter,nwcalls,nwtotit end if if ( .not. scale ) then write(* ,1090) f,cnorm write(10,1090) f,cnorm else write(* ,1100) f,fu,cnorm,cnormu write(10,1100) f,fu,cnorm,cnormu end if write(* ,1110) snorm,nlpsupn,xind,xsupn,rsupn write(10,1110) snorm,nlpsupn,xind,xsupn,rsupn if ( iprintout .ge. 4 ) then if ( nprint .ne. 0 ) then write(* ,1120) nprint,(x(i),i=1,nprint) write(10,1120) nprint,(x(i),i=1,nprint) end if if ( mprint .ne. 0 ) then write(* ,1130) mprint,(lambda(i),i=1,mprint) write(* ,1140) mprint,(rho(i),i=1,mprint) write(10,1130) mprint,(lambda(i),i=1,mprint) write(10,1140) mprint,(rho(i),i=1,mprint) end if end if end if C Save intermediate data for crash report open(20,file='algencan-tabline.out') write(20,1400) fu,cnormu,f,cnorm,nlpsupn,inform,9,n,m,outiter, + totiter,0,nwcalls,nwtotit,msqcalls,msqtotit,999.0d0 close(20) C ================================================================== C Test stopping criteria C ================================================================== C Test feasibility, optimality and complementarity if ( max( snorm, cnormu ) .le. epsfeas .and. + ( nlpsupn .le. epsopt .or. + ( outiter .gt. 0 .and. nbv + niv .gt. 0 ) ) .and. + ( nbv + niv .eq. 0 .or. + satisfyIntegerConstraints(n,x,l,u,variableTypes) ) ) then alinfo = 0 if ( iprintout .ge. 1 ) then write(*, 1300) write(10,1300) end if return end if C Test whether the number of iterations is exhausted if ( outiter .ge. maxoutit ) then alinfo = 1 if ( iprintout .ge. 1 ) then write(*, 1310) write(10,1310) end if return end if C Test whether the problem seems to be infeasible if ( snorm .gt. epsfeas .and. snorm .ge. snormprev ) then outitnp = outitnp + 1 else outitnp = 0 end if if ( outitnp .ge. maxoutitnp .or. rsupn .gt. rhomax ) then alinfo = 2 if ( iprintout .ge. 1 ) then write(*, 1320) write(10,1320) end if return end if C ================================================================== C Near the solution, try to solve the KKT system by Newton's method C ================================================================== if ( .not. skipacc .and. truehl .and. ( nwcalls .gt. 0 .or. + ( snorm .le. epsfeas12 .and. nlpsupn .le. epsopt12 ) .or. + ( snorm .le. epsfeas14 .and. nlpsupn .le. epsopt14 .and. + outiter .gt. 0 .and. geninfo .ne. 0 ) ) ) then nwcalls = nwcalls + 1 do i = 1,n xa(i) = x(i) end do do i = 1,m lambdaa(i) = lambar(i) end do scaletmp = scale if ( scale ) then do i = 1,m lambdaa(i) = lambdaa(i) * sf / sc(i) end do scale = .false. end if call newtonkkt(n,xa,l,u,m,lambdaa,equatn,linear,epsfeas, + epsopt,fa,cnorma,nlpsupna,nwiter,msqiter,nwinfo,inform) nwtotit = nwtotit + nwiter if ( nwinfo .eq. 5 .or. nwinfo .eq. 6 ) then skipacc = .true. end if if ( msqiter .gt. 0 ) then msqcalls = msqcalls + 1 msqtotit = msqtotit + msqiter end if scale = scaletmp if ( inform .lt. 0 ) return C Test feasibility, optimality and complementarity if ( cnorma .le. epsfeas .and. nlpsupna .le. epsopt ) then fu = fa cnormu = cnorma f = fa cnorm = cnorma snorm = 0.0d0 nlpsupn = nlpsupna do i = 1,n x(i) = xa(i) end do do i = 1,m lambar(i) = lambdaa(i) end do alinfo = 0 if ( iprintout .ge. 1 ) then write(*, 1300) write(10,1300) end if return end if end if C ================================================================== C Iteration C ================================================================== outiter = outiter + 1 C ================================================================== C Solve the augmented Lagrangian subproblem C ================================================================== C Set optimality requeriment for the subproblem if ( snorm .gt. epsfeas12 .or. nlpsupn .gt. epsopt12 .or. + ( nbv + niv ) .gt. 0 ) then epsopk = sqrt( epsopt ) else epsopk = max( epsopt, 0.1d0 * epsopk ) end if if ( outiter .eq. 1 ) then maxit = 10 else maxit = 1000 end if C Call the inner-solver if ( nbv + niv .eq. 0 ) then call gencan(n,x,l,u,m,lambar,equatn,linear,rho,.false., + epsfeas,epsopk,maxit,iter,al,nl,nlpsupn,cnorm,cnormu,geninfo, + inform) else call minlp_solver(n,x,l,u,m,lambar,equatn,linear,rho,.false., + epsfeas,epsopk,maxit,iter,al,nl,nlpsupn,cnorm,cnormu,geninfo, + inform,nbv,niv,variableTypes,minlpinfo) end if totiter = totiter + iter if ( inform .lt. 0 ) return C ================================================================== C Prepare for the next iteration C ================================================================== C Compute current point norm xsupn = 0.0d0 do i = 1,n if ( abs( x(i) ) .ge. xsupn ) then xsupn = abs( x(i) ) xind = i end if end do C Compute objective function and constraints call sevalobjc(n,x,f,m,c,inform) if ( inform .lt. 0 ) return fu = f * sf C Compute complementarity and feasibility violations snormprev = snorm cnorm = 0.0d0 snorm = 0.0d0 do i = 1,m if ( equatn(i) ) then cnorm = max( cnorm, abs( c(i) ) ) sigma(i) = c(i) else cnorm = max( cnorm, c(i) ) sigma(i) = max( c(i), - lambar(i) / rho(i) ) end if snorm = max( snorm, abs( sigma(i) ) ) end do C Compute unscaled objective function and constraints norm if ( scale ) then fu = f * sf cnormu = 0.0d0 do i = 1,m if ( equatn(i) ) then cnormu = max( cnormu, abs( sc(i) * c(i) ) ) else cnormu = max( cnormu, sc(i) * c(i) ) end if end do else fu = f cnormu = cnorm end if C Update Lagrange multipliers approximation do i = 1,m call evaldpdy(c(i),rho(i),lambar(i),equatn(i),lambda(i)) lambar(i) = max( lammin, min( lambda(i), lammax ) ) end do C Compute continuous projected Lagrangian gradient norm call sevalnl(n,x,m,lambda,equatn,linear,nl,inform) if ( inform .lt. 0 ) return nlpsupn = 0.0d0 do i = 1,n if( variableTypes(i) .ne. CONTINUOUS ) cycle nlpi = x(i) - nl(i) if ( l(i) .le. nlpi .and. nlpi .le. u(i) ) then nlpi = - nl(i) else nlpi = max( l(i), min( nlpi, u(i) ) ) - x(i) end if nlpsupn = max( nlpsupn, abs( nlpi ) ) end do C Update penalty parameters if ( outiter .eq. 1 ) then call comprhoini(f,m,c,equatn,rho) if ( iprintout .ge. 5 ) then write(*, 1200) write(10,1200) end if else C Note that |sigma(i)| <= eps implies c(i) <= eps. if ( rhoiden ) then if ( max( snorm, cnormu ) .gt. epsfeas .and. + snorm .gt. rhofrac * snormprev ) then do i = 1,m rho(i) = rhomult * rho(i) end do if ( iprintout .ge. 5 ) then write(*, 1210) rhomult write(10,1210) rhomult end if else if ( iprintout .ge. 5 ) then write(*, 1230) write(10,1230) end if end if else rhoincr = 0 do i = 1,m if ( max(abs(sigma(i)), sc(i)*c(i)) .gt. epsfeas .and. + abs(sigma(i)) .gt. rhofrac*abs(sigpre(i)) ) then rho(i) = rhomult * rho(i) rhoincr = rhoincr + 1 end if end do if ( iprintout .ge. 5 ) then if ( rhoincr .ne. 0 ) then write(*, 1220) rhoincr,m,rhomult write(10,1220) rhoincr,m,rhomult else write(*, 1230) write(10,1230) end if end if end if end if C Compute maximum rho rsupn = 0.0d0 do i = 1,m rsupn = max( rsupn, rho(i) ) end do C ================================================================== C Iterate C ================================================================== go to 100 C ================================================================== C End of main loop C ================================================================== C NON-EXECUTABLE STATEMENTS 1000 format(/,1X,'Entry to ALGENCAN.', + /,1X,'Number of variables : ',I7, + /,1X,'Number of constraints: ',I7) 1010 format(/,1X,'Lower bounds (first ',I7,' components):', + /,6(1X,1P,D11.4)) 1020 format(/,1X,'Upper bounds (first ',I7,' components):', + /,6(1X,1P,D11.4)) 1030 format(/,1X,'It',15X,'f',4X,'cnrm',4X,'snrm',2X,'nlpnrm',4X, + 'xnrm',1X,'penalty',2X,'innit',1X,'nwcalls',2X,'nwit') 1040 format(I3,1P,D16.8,5(1X,1P,D7.1),1X,I5,A1,2(1X,I6)) 1050 format(/,1X,'It',4X,'unscaled f & cnrm',6X,'scaled f & cnrm', + 3X,'snrm',1X,'nlpnrm',1X,'penalty',1X,'innit',2X,'nwkkt') 1060 format(I3,2(1P,D14.6,1X,1P,D6.0),3(1X,1P,D6.0),1X,I5,A1,1X,I2,1X, + I3) 1070 format(/,1X,'ALGENCAN OUTER ITERATION = ',I7) 1080 format(/,1X,'Up-to-now total number of iterations = ',I7, + /,1X,'Up-to-now acceleration trials = ',I7, + /,1X,'Up-to-now total number of Newton iterations = ',I7) 1090 format(/,1X,'Functional value = ', + 1P,D24.16, + /,1X,'Sup-norm of constraints = ', + 17X,1P,D7.1) 1100 format(/,1X,'Functional value (scaled = ',1P,D7.1, + ') = ',1P,D24.16, + /,1X,'Sup-norm of constraints (scaled = ',1P,D7.1, + ') = ',17X,1P,D7.1) 1110 format( 1X,'Sup-norm of complementarity-feasibility = ', + 17X,1P,D7.1, + /,1X,'Sup-norm of the Lagrangian projected gradient = ', + 17X,1P,D7.1, + /,1X,'Sup-norm of x (attained at x_{', I7,'}) = ', + 17X,1P,D7.1, + /,1X,'Largest penalty parameter = ', + 17X,1P,D7.1) 1120 format(/,1X,'Current point (first ',I7,' components):', + /,6(1X,1P,D11.4)) 1130 format(/,1X,'Updated Lagrange multipliers (first ',I7, + 1X,'components):', + /,6(1X,1P,D11.4)) 1140 format(/,1X,'Updated penalty parameters (first ',I7, + 1X,'components):', + /,6(1X,1P,D11.4)) 1200 format(/,1X,'Penalty parameters are re-initiated after the ', + 'resolution of the first',/,1X,'subproblem.') 1210 format(/,1X,'The desired infeasibility was not achieved. ', + 'The penalty',/,1X,'parameter will be ', + 'increased multiplying by rhofrac = ',1PD11.4,'.') 1220 format(/,1X,'The desired infeasibility was not achieved in ',I7, + 1X,'over ',I7,/,1X,'constraints.',/,1X,' Penalty ', + 'parameters will be increased multiplying by ', + 'rhofrac = ',1PD11.4,'.') 1230 format(/,1X,'Desired feasibility improvement was achieved, ', + 'so, penalty parameters will not',/,1X,'be ', + 'modified.') 1300 format(/,1X,'Flag of ALGENCAN: Solution was found.') 1310 format(/,1X,'Flag of ALGENCAN: Maximum of iterations reached.') 1320 format(/,1X,'Flag of ALGENCAN: The problem seems to be ', + 'infeasible.') 1330 format(/,1X,'Flag of ALGENCAN: Are you requiring an exaggerated ', + 'optimality tolerance?') 1400 format(1X,1P,D24.16,1X,1P,D7.1,1X,1P,D24.16,1X,1P,D7.1,1X,1P,D7.1, + 1X,I3,1X,I1,1X,I6,1X,I6,1X,I2,1X,I7,1X,I7,1X,I2,1X,I7,1X, + I7,1X,I7,0P,F8.2) end C ****************************************************************** C ****************************************************************** subroutine comprhoini(f,m,c,equatn,rho) C SCALAR ARGUMENTS integer m double precision f C ARRAY ARGUMENTS logical equatn(m) double precision c(m),rho(m) C Consider the Augmented Lagrangian function C C al = f(x) + \sum_i P(lambda_i,c_i(x),rho), C C where C C P(lambda,y,rho) = y ( lambda + 0.5 rho y ), C C If c_i(x) is an equality constraint or lambda + rho y > 0, and C C P(lambda,y,rho) = - 0.5 lambda^2 / rho, C C otherwise. C C Assuming that lambda_i = 0 for all i, it is clear that C C P(lambda_i,c_i(x),rho) = 0.5 rho c_i(x)^2 C C and that the value of rho that balances f(x) and C C \sum_i P(lambda_i,c_i(x),rho) is given by C C rho = f(x) / ( 0.5 \sum_i c_i(x)^2 ). #include "machconst.com" C LOCAL SCALARS integer i double precision rhoini,sumc sumc = 0.0d0 do i = 1,m if ( equatn(i) .or. c(i) .gt. 0.0d0 ) then sumc = sumc + 0.5d0 * c(i) ** 2 end if end do rhoini = 10.0d0 * max( 1.0d0, abs( f ) ) / max( 1.0d0, sumc ) rhoini = max( macheps12, min( rhoini, 1.0d0 / macheps12 ) ) do i = 1,m rho(i) = rhoini end do end C ****************************************************************** C ****************************************************************** subroutine sevalfeas(n,x,m,equatn,cnorm,cnormu,inform) implicit none C SCALAR ARGUMENTS double precision cnorm,cnormu integer inform,m,n C ARRAY ARGUMENTS logical equatn(m) double precision x(n) #include "dim.par" #include "graddat.com" #include "scaling.com" C LOCAL SCALARS integer i C Subroutine sevalnal is always called by GENCAN with the same point C before calling this subroutine. Then, constraints are computed and C saved. C Compute infeasibility cnorm = 0.0d0 do i = 1,m if ( equatn(i) ) then cnorm = max( cnorm, abs( c(i) ) ) else cnorm = max( cnorm, c(i) ) end if end do if ( scale ) then cnormu = 0.0d0 do i = 1,m if ( equatn(i) ) then cnormu = max( cnormu, abs( sc(i) * c(i) ) ) else cnormu = max( cnormu, sc(i) * c(i) ) end if end do else cnormu = cnorm end if end C ****************************************************************** C ****************************************************************** logical function satisfyIntegerConstraints(n,x,l,u,variableTypes) implicit none C SCALAR ARGUMENTS integer n C ARRAY ARGUMENTS integer variableTypes(n) double precision x(n),l(n),u(n) #include "algconst.par" C LOCAL SCALARS integer i C FUNCTION logical isInteger satisfyIntegerConstraints = .true. do i = 1,n if(variableTypes(i) .eq. BINARY .or. + variableTypes(i) .eq. INTEGER) then if(x(i) .lt. l(i) .or. x(i) .gt. u(i) .or. + .not. isInteger(x(i))) then satisfyIntegerConstraints = .false. return end if end if end do return end C ************************************************************************ C ************************************************************************ logical function isInteger(a) implicit none C This function verifies if the real argument 'a' is integer, i.e., C if "a = floor(a)". Return 'true' if the value of 'a' is integer, C and 'false' otherwise. C SCALAR ARGUMENT double precision a #include "algconst.par" isInteger = .false. if( abs(a) .le. EPSILON_INT ) then isInteger = .true. else if(min(abs(a - floor(a)), abs(a - ceiling(a))) / + max(1.0d0, abs(a)) .le. EPSILON_INT) then isInteger = .true. end if return end