module MINLP use Constants use BranchAndBound_Options logical, parameter :: RANDOM_DIRECTION = .true. contains ! ************************************************************************ ! ************************************************************************ subroutine BranchAndBound(n,incumbent,lower,upper,m,lambda,equatn,linear, & rho,constf,epsfeas,epsopk,maxit,optSolverIter,objectiveValue,g,gpsupn, & cnorm,cnormu,geninfo,inform,numBinary,numInteger,iteration, & solutionFound,nodeSelection,branchingVariableStrategy, & multistartStrategy,maxSolverCalls,numNodesGenerated,numProblemsSolved, & iterationSolutionFound,variableTypes,minlpinfo,options) use Auxiliar use LinkedList use OptSolver use Problems use ProblemType use PseudoCosts use Random use Structures use VariableSelection implicit none ! SCALAR ARGUMENTS logical, intent(in) :: constf integer, intent(in) :: m,n,maxit,numBinary,numInteger integer, intent(in) :: maxSolverCalls integer, intent(in) :: nodeSelection integer, intent(in) :: branchingVariableStrategy,multistartStrategy integer, intent(out) :: iteration,geninfo,optSolverIter,minlpinfo integer, intent(in out) :: inform real(kind=8), intent(out) :: objectiveValue,gpsupn real(kind=8), intent(in) :: epsfeas,epsopk real(kind=8), intent(in out) :: cnorm,cnormu logical, intent(out) :: solutionFound integer, intent(out) :: iterationSolutionFound integer, intent(out) :: numNodesGenerated,numProblemsSolved ! ARRAY ARGUMENTS logical, intent(in) :: equatn(m),linear(m) integer, intent(in) :: variableTypes(n) real(kind=8), intent(in) :: lambda(m),rho(m) real(kind=8), intent(in out) :: lower(n),upper(n),incumbent(n) real(kind=8), intent(out) :: g(n) type(bb_options), intent(in) :: options ! Parameters of the subroutine: ! ============================= ! ! On Entry: ! ========= ! ! INCUMBENT real(kind=8) incumbent(n): Initial estimation of the solution. ! ----------------------------------- ! ! N integer: Number of variables. ! --------- ! ! LOWER real(kind=8) lower(n): Lower bounds of the variables. ! --------------------------- ! ! UPPER real(kind=8) upper(n): Upper bounds of the variables. ! --------------------------- ! ! M integer: Number of constraints (equalities plus inequalities) ! --------- without considering the bound constraints. ! ! LAMBDA real(kind=8) lambda(m): ! ----------------------------- ! ! Estimation of the Lagrange multipliers. ! ! EQUATN logical equatn(m): ! ------------------------ ! ! Logical array that, for each constraint, indicates whether the ! constraint is an equality constraint (TRUE) or an inequality ! constraint (FALSE). ! ! LINEAR logical linear(m) ! ------------------------ ! ! Logical array that, for each constraint, indicates whether the ! constraint is a linear constraint (TRUE) or a nonlinear constraint ! (FALSE). ! ! NUMBINARY integer: Number of binary variables. ! ----------------- ! ! NUMINTEGER integer: Number of integer variables. ! ------------------ ! ! NODESELECTION integer: ! --------------------- ! ! The rule for selecting the next node. ! ! BRANCHINGVARIABLESTRATEGY integer: ! --------------------------------- ! ! The strategy to be used to select the branching variable. The ! possible strategies are: ! ! 1) Choose the non-integer variable whose index is smaller ! (LOWEST_INDEX_FIRST). ! ! 2) Choose the non-integer variable which is most fractional, ! that is, the variable whose value is farther than an integer ! (MOST_FRACTIONAL). ! ! MULTISTARTSTRATEGY integer: ! -------------------------- ! ! The strategy to be used joint with multistart. It can be: ! ! 1) Among all feasible solutions found for a subproblem, choose ! the best one (BEST_FEASIBLE). ! ! 2) Solve a subproblem at most 'maxSolverCalls' times and stop ! with the first feasible solution (FIRST_FEASIBLE). ! ! If no feasible solution was found for a subproblem, this ! subproblem will be declared infeasible. ! ! MAXSOLVERCALLS integer: Maximum number of attempts to solve the ! ---------------------- same subproblem. ! ! On Return: ! ========== ! ! INCUMBENT real(kind=8) incumbent(n): The solution found for the problem. ! ----------------------------------- ! ! OBJECTIVEVALUE real(kind=8): Value of the solution found. ! --------------------------- ! ! ITERATION integer: Number of iterations of the algorithm. ! ----------------- ! ! SOLUTIONFOUND logical: Indicates whether a solution was ! --------------------- found (TRUE) or not (FALSE). ! ! ITERATIONSOLUTIONFOUND integer: Iteration at which the solution was ! ------------------------------ found. ! ! NUMNODESGENERATED integer: Number of generated nodes. ! ------------------------- ! ! NUMPROBLEMSSOLVED integer: Number of (sub)problems solved. ! ------------------------- ! LOCAL SCALARS integer :: branchingIndex,m_ logical :: satisfy_domain,feasibleSolutionFound real(kind=8) :: f,best_f,round_f,cutoff,boundTolerance integer :: listType real(kind=8), save :: seed = 101063.0d0 real(kind=8), save :: seed_prob = 7919.0d0 integer :: j,k,nextID,jj integer :: i integer :: multiStartIteration integer :: variableType integer :: maxSolverCallsInc,updatesBeforeTrust integer :: iter_,iter_2 real(kind=8) :: f_,f_2,lowerBound,frac,S0 ! Counters integer :: maxSize,pruningByBoundCount,pruningByStoredBound,& pruningByInfeas,pruningByOpt integer, save :: totalNumNodesGenerated = 0 integer, save :: totalNumProblemsSolved = 0 ! LOCAL ARRAYS real(kind=8) :: l(n),u(n),x(n),best_x(n),round_x(n),l_(n),u_(n) logical :: switched,circlePackingInitialPoint,solveRoundedSolution real :: time,startTime ! OTHER LOCAL VARIABLES type(data) :: data_ type(domain) :: mainDomain,subdomain type(list) :: openProblemsList type(problem) :: p,p1,p2 type(set), allocatable :: discrete_sets(:) ! ************** ! INITIALIZATION ! ************** boundTolerance = getBoundTolerance(options) solveRoundedSolution = getSolveRoundedSolution(options) circlePackingInitialPoint = .false. call cpu_time(startTime) cutoff = getCutOff(options) select case(nodeSelection) case(DEPTH_FIRST) listType = STACK case(SWITCH_DF_BBR) listType = STACK case(BREADTH_FIRST) listType = QUEUE case(BEST_BOUND_RULE) listType = PRIORITY end select switched = .false. lowerBound = LARGEST_VALUE optSolverIter = 0 iteration = 0 minlpinfo = 0 nextID = 2 solutionFound = .false. iterationSolutionFound = -1 maxSolverCallsInc = getNumNodeResolutionsIncrease(options) updatesBeforeTrust = getNumberOfUpdatesBeforeTrust(options) ! Counters. maxSize = 1 pruningByBoundCount = 0 pruningByStoredBound = 0 pruningByInfeas = 0 pruningByOpt = 0 numNodesGenerated = 0 numProblemsSolved = 0 ! Initialize the variables pseudo-costs. if(branchingVariableStrategy .eq. PSEUDO_COSTS) then call initializePseudoCosts(n) end if call initialize(n,numBinary,numInteger,mainDomain,discrete_sets,& variableTypes) ! Create the main problem. ! An empty array is given by (/ (0, i=1,0) /). p = createProblem(n,incumbent,createDomain(0,(/ (0, jj=1,0) /),& (/ (BINARY, jj=1,0) /),(/ (0.0d0, jj=1,0) /) ),1,0,cutoff,& -1,0.0d0,.false.,0) ! Create the list of open problems and add the main problem to it. data_%problem = p openProblemsList = createList(listType) call add(openProblemsList,data_) ! Verify whether the initial point is a feasible solution. if(satisfyDomain(n,incumbent,mainDomain,discrete_sets) .and. & satisfyConstraints(n,incumbent,lower,upper)) then ! The initial point is an integer feasible solution. ! Evaluate the objective function. call sevalal(n,incumbent,m,lambda,rho,equatn,linear,f,inform) if(inform .lt. 0) then call deallocateDomain(mainDomain) call destroyList(openProblemsList) deallocate(discrete_sets) if(branchingVariableStrategy .eq. PSEUDO_COSTS) then call deallocatePseudoCosts() end if return end if if(isNumber(f)) then ! Update the incumbent solution. objectiveValue = f solutionFound = .true. iterationSolutionFound = 0 end if else ! Initialize the objective value of the best solution with ! "infinity". objectiveValue = LARGEST_VALUE end if ! ********* ! MAIN LOOP ! ********* ! ! The problems are removed from the open problems list until it ! becomes empty. do while(.not. empty(openProblemsList)) call flush() if(nodeSelection .eq. SWITCH_DF_BBR .and. & iteration .eq. getNumIterationsForBestBoundRule(options)) then call changeListType(openProblemsList,PRIORITY) end if numNodesGenerated = numNodesGenerated + 1 iteration = iteration + 1 ! Get a subprolem. data_ = remove(openProblemsList) p = data_%problem ! Add the bounds on the variables. l = lower u = upper call addBounds(p%domain,l,u,n) ! Get the initial point for this subproblem. x = p%initialPoint ! PRUNING BY (STORED) BOUND ! ------------------------- ! ! Check whether the stored lower bound for this subproblem is ! greater than or equal to the value of the incumbent ! solution. If it occurs, the subproblem can be discarded, ! since it cannot provide a better solution. if(solutionFound .and. p%geninfoFather .eq. 0 .and. & boundTolerance * p%lowerBound .ge. objectiveValue) then pruningByStoredBound = pruningByStoredBound + 1 call deallocateProblem(p) ! Continue with the next subproblem from the list. cycle end if ! ******************* ! MULTISTART STRATEGY ! ******************* ! ! The relaxation of this subproblem is solved 'maxSolverCalls' ! times following the multistart strategy. Initial points are ! randomly generated according to the following rule: ! ! We have two procedures to apply joint with multistart: ! ! 1) The relaxation of the subproblem is solved at most ! 'maxSolverCalls' times and we accept the first feasible ! relaxed solution found. This strategy is activated by ! choosing FIRST_FEASIBLE as subproblem resolution strategy. ! ! 2) The relaxation of the subproblem is solved exactly ! 'maxSolverCalls' times and we get the best feasible relaxed ! solution found among all solutions found. This strategy is ! activated by choosing BEST_FEASIBLE as subproblem resolution ! strategy. ! Initially, the value of the best relaxed solution found is ! set to "infinity". best_f = LARGEST_VALUE ! This variable indicates whether a solution for the relaxed ! problem was found. feasibleSolutionFound = .false. ! Initialize the seed to generate random initial points for ! this subproblem. if(iteration .gt. 1) then call getSeed(seed,l,u,n) else seed = real(84543787,kind=8) end if do multiStartIteration = 1,max(1,maxSolverCalls + maxSolverCallsInc*p%depth) ! Get a random initial point for this subproblem. If i = 1, ! we consider the initial point previous defined (stored). if(multiStartIteration .gt. 1) then call randomInitialPoint(n,x,seed) end if if(circlePackingInitialPoint) then call defCirclePackingProblem() call optimizationSolver(n,x,l,u,m,lambda,equatn,& linear,rho,constf,epsfeas,epsopk,maxit,iter_2,f_2,g,& gpsupn,cnorm,cnormu,geninfo,inform) optSolverIter = optSolverIter + iter_2 call defRectanglePackingProblem() end if call solveProblem(p,n,x,l,u,m,lambda,equatn,linear,rho,constf,& epsfeas,epsopk,maxit,iter_,f,g,gpsupn,cnorm,cnormu,geninfo,& inform,numProblemsSolved) optSolverIter = optSolverIter + iter_ ! Update lower bound on the objective function. if(f < lowerBound) then lowerBound = f end if ! Updating the best solution. if(satisfyConstraints(n,x,l,u)) then ! A feasible solution for the relaxed problem was found. feasibleSolutionFound = .true. if(f < best_f) then ! Update the best relaxed solution found. best_f = f best_x = x end if ! Verify whether the relaxed solution found is feasible ! for the main problem, that is, if it also satisfies the ! integrality constraints. if(satisfyDomain(n,x,mainDomain,discrete_sets) .and. & isNumber(f)) then if(.not. solutionFound .or. f < objectiveValue) then if(PRINT_) then call cpu_time(time) write(*,*) ' Solution updated: ', f, ' time = ', time end if ! Update the incumbent solution. objectiveValue = f incumbent = x solutionFound = .true. iterationSolutionFound = iteration if(objectiveValue .le. cutoff) then exit end if end if else ! Round solution x. call round(n,x,round_x,mainDomain) if(solveRoundedSolution) then ! FIX THE INTEGER VARIABLES AND SOLVE THE PROBLEM ! ONLY FOR THE CONTINUOUS VARIABLES. do i = 1,n if(variableTypes(i) .eq. BINARY .or. & variableTypes(i) .eq. INTEGER) then l_(i) = round_x(i) u_(i) = l_(i) else l_(i) = l(i) u_(i) = u(i) end if end do call optimizationSolver(n,round_x,l_,u_,m,lambda,equatn,& linear,rho,constf,epsfeas,1.0D-8,maxit,iter_,f_,g,& gpsupn,cnorm,cnormu,geninfo,inform) numProblemsSolved = numProblemsSolved + 1 end if ! Check whether the rounded solution 'round_x' is a ! feasible (integer) solution for the main problem. if(satisfyDomain(n,round_x,mainDomain,discrete_sets) .and. & satisfyConstraints(n,round_x,l,u)) then ! Evaluate the objective function at round_x. call sevalal(n,round_x,m,lambda,rho,equatn,linear,& round_f,inform) if(inform .lt. 0) then call deallocateDomain(mainDomain) call destroyList(openProblemsList) deallocate(discrete_sets) if(branchingVariableStrategy .eq. PSEUDO_COSTS) then call deallocatePseudoCosts() end if return end if ! Rounded solution is feasible for the main problem. feasibleSolutionFound = .true. ! Check whether the rounded solution is better than ! the best relaxed solution found for this ! subproblem. if(round_f < best_f) then ! Update the best solution found. best_f = round_f best_x = round_x end if ! Check whether the rounded solution is better than ! the incumbent solution. if(inform .eq. 0 .and. isNumber(f) .and. & (.not. solutionFound .or. round_f < objectiveValue)) then ! The rounded solution is feasible for the main ! problem and it is better than the incumbent ! solution. if(PRINT_) then call cpu_time(time) write(*,*) ' Solution updated (rounded solution):',& round_f, ' time = ', time end if ! Update the incumbent solution. objectiveValue = round_f incumbent = round_x solutionFound = .true. iterationSolutionFound = iteration if(objectiveValue .le. cutoff) then exit end if end if end if ! End of testing rounded solution. end if end if ! Accept the first relaxed feasible solution if this ! strategy was chosen. if(feasibleSolutionFound .and. & multistartStrategy .eq. FIRST_FEASIBLE) then ! Exit from the multistart loop. exit end if end do ! best_x is the best feasible solution found fir the relaxed ! subproblem. f = best_f x = best_x if(iteration .eq. 1) then S0 = 0.0d0 do i = 1,n if(variableTypes(i) .eq. BINARY .or. & variableTypes(i) .eq. INTEGER) then frac = x(i) - real(floor(x(i)),kind=8) if(frac .le. 0.5d0) then S0 = S0 + frac else S0 = S0 + 1.0d0 - frac end if end if end do end if ! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ! EXIT ! exit ! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ! ************************** ! END OF MULTISTART STRATEGY ! ************************** p%solved = .true. ! PRUNING BY INFEASIBILITY ! ------------------------ ! ! If no solution was found for the relaxed problem, this ! subproblem is declared to be infeasible. if(.not. feasibleSolutionFound) then pruningByInfeas = pruningByInfeas + 1 call deallocateProblem(p) ! Continue with the next subproblem from the list. cycle end if ! UPDATE PSEUDO-COSTS ! ------------------- if(branchingVariableStrategy .eq. PSEUDO_COSTS) then if(p%indexBranchingVariable .ge. 1 .and. & p%indexBranchingVariable .le. n) then call updatePseudoCost(p%indexBranchingVariable,& p%branchingVariableValue,p%downBranch,& p%lowerBound,f) end if end if if(objectiveValue .le. cutoff) then pruningByOpt = pruningByOpt + 1 exit end if ! Check whether the solution found for the relaxed subproblem ! is a feasible solution for the main problem. satisfy_domain = satisfyDomain(n,x,mainDomain,discrete_sets) ! PRUNING BY OPTIMALITY ! --------------------- ! ! If the feasible solution of the relaxed subproblem also ! satisfies the integrality constraints then it is a feasible ! solution for the main problem. if(isNumber(f) .and. satisfy_domain) then pruningByOpt = pruningByOpt + 1 if(.not. solutionFound .or. f < objectiveValue) then if(PRINT_) then write(*,*) ' Iteration: ',iteration write(*,*) ' Solution updated: ',f end if ! Update the incumbent solution. objectiveValue = f incumbent = x solutionFound = .true. iterationSolutionFound = iteration if(objectiveValue .le. cutoff) then exit end if end if call deallocateProblem(p) ! Continue with the next subproblem from the list. cycle end if ! PRUNING BY BOUND ! ---------------- ! ! If the value of the solution of the relaxed subproblem is not ! less than the value of the incumbent solution, this node can ! be pruned. if(solutionFound .and. boundTolerance * f .ge. objectiveValue) then pruningByBoundCount = pruningByBoundCount + 1 call deallocateProblem(p) ! Continue with the next subproblem from the list. cycle end if ! ******************* ! BRANCHING PROCEDURE ! ******************* ! Choose a branching variable following the strategy specified by ! 'branchingVariableStrategy' variable. call chooseBranchingVariable(x,mainDomain,discrete_sets,& variableType,branchingVariableStrategy,branchingIndex,options) m_ = size(p%domain%type) + 1 ! BINARY VARIABLE ! --------------- if(variableType .eq. BINARY) then call branch(p,p1,p2,n,x,branchingIndex,variableType,f,nextID) if(RANDOM_DIRECTION .and. drand(seed_prob) .lt. 0.5d0) then data_%problem = p2 call add(openProblemsList,data_) data_%problem = p1 call add(openProblemsList,data_) else data_%problem = p1 call add(openProblemsList,data_) data_%problem = p2 call add(openProblemsList,data_) end if ! INTEGER VARIABLE ! ---------------- else if(variableType .eq. INTEGER) then call branch(p,p1,p2,n,x,branchingIndex,variableType,f,nextID) if(RANDOM_DIRECTION .and. drand(seed_prob) .lt. 0.5d0) then data_%problem = p2 call add(openProblemsList,data_) data_%problem = p1 call add(openProblemsList,data_) else data_%problem = p1 call add(openProblemsList,data_) data_%problem = p2 call add(openProblemsList,data_) end if ! DISCRETE VARIABLE ! ----------------- else if(variableType .eq. DISCRETE) then do j = 1,size(discrete_sets) if(discrete_sets(j)%index .eq. branchingIndex) then ! For each possible value for the branching variable, a ! subproblem is created and a constraint is added in ! order to fix the value of the variable. do k = 1,size(discrete_sets(j)%elements) subdomain = createDomain(m_,& (/ p%domain%indices,branchingIndex /),& (/ p%domain%type,EQUAL /),& (/ p%domain%value,discrete_sets(j)%elements(k) /)) p1 = createProblem(n,x,subdomain,nextID,p%id,f,& branchingIndex,x(branchingIndex),.false.,p%depth + 1) p1%geninfoFather = p%geninfo data_%problem = p1 call add(openProblemsList,data_) nextID = nextID + 1 end do exit end if end do end if call deallocateProblem(p) maxSize = max(maxSize, getSize(openProblemsList)) end do ! END OF MAIN LOOP call deallocateDomain(mainDomain) deallocate(discrete_sets) if(branchingVariableStrategy .eq. PSEUDO_COSTS) then call deallocatePseudoCosts() end if ! FIX THE INTEGER VARIABLES AND SOLVE THE PROBLEM ONLY FOR THE ! CONTINUOUS VARIABLES. do i = 1,n if(variableTypes(i) .eq. BINARY .or. variableTypes(i) .eq. INTEGER) then l(i) = incumbent(i) u(i) = l(i) else l(i) = lower(i) u(i) = upper(i) end if end do call optimizationSolver(n,incumbent,l,u,m,lambda,equatn,linear,rho,constf,& epsfeas,1.0D-8,maxit,iter_,objectiveValue,g,gpsupn,cnorm,cnormu,& geninfo,inform) totalNumNodesGenerated = totalNumNodesGenerated + numNodesGenerated totalNumProblemsSolved = totalNumProblemsSolved + numProblemsSolved call cpu_time(time) if(PRINT_) then write(*,*) write(*,*) 'Time: ', time write(*,*) 'Solution found: ', objectiveValue write(*,*) 'Nodes generated: ', totalNumNodesGenerated write(*,*) 'Subproblems solved: ', totalNumProblemsSolved write(*,*) end if ! Print formats. !100 format('Iteration Objective') !200 format (I9,1X,1PE20.10) end subroutine BranchAndBound ! ************************************************************************ ! ************************************************************************ subroutine branch(p,p1,p2,n,x,branchingIndex,variableType,f,nextID) use Problems implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n,branchingIndex,variableType integer, intent(in out) :: nextID real(kind=8), intent(in) :: x(n),f ! OTHER ARGUMENTS type(problem), intent(in) :: p type(problem), intent(out) :: p1,p2 ! LOCAL SCALARS integer :: m_ real(kind=8) :: boundValue ! OTHER SCALARS type(domain) :: subdomain ! p1 must be the down branch and p2 must be the up branch. ! Number of constraints on the subproblems. m_ = size(p%domain%type) + 1 ! BINARY VARIABLE ! --------------- if(variableType .eq. BINARY) then ! First problem: we set x_i = 0. subdomain = createDomain(m_,(/ p%domain%indices,branchingIndex /),& (/ p%domain%type,EQUAL /),(/ p%domain%value,0.0d0 /)) p1 = createProblem(n,x,subdomain,nextID,p%id,f,branchingIndex,& x(branchingIndex),.true.,p%depth + 1) p1%geninfoFather = p%geninfo nextID = nextID + 1 ! Second problem: we set x_i = 1. subdomain = createDomain(m_,(/ p%domain%indices,branchingIndex /),& (/ p%domain%type,EQUAL /),(/ p%domain%value,1.0d0 /)) p2 = createProblem(n,x,subdomain,nextID,p%id,f,branchingIndex,& x(branchingIndex),.false.,p%depth + 1) p2%geninfoFather = p%geninfo nextID = nextID + 1 ! INTEGER VARIABLE ! ---------------- else if(variableType .eq. INTEGER) then boundValue = real(floor(x(branchingIndex)),kind=8) ! First problem. We add the following constraint: ! x_i <= floor(x_i). subdomain = createDomain(m_,(/ p%domain%indices,branchingIndex /),& (/ p%domain%type,LESS_EQUAL /),(/ p%domain%value,boundValue /)) p1 = createProblem(n,x,subdomain,nextID,p%id,f,branchingIndex,& x(branchingIndex),.true.,p%depth + 1) p1%geninfoFather = p%geninfo nextID = nextID + 1 ! Second problem. We add the following constraint: ! x_i >= floor(x_i) + 1. subdomain = createDomain(m_,(/ p%domain%indices,branchingIndex /),& (/ p%domain%type,GREATER_EQUAL /),& (/ p%domain%value,boundValue + 1.0d0 /)) p2 = createProblem(n,x,subdomain,nextID,p%id,f,branchingIndex,& x(branchingIndex),.false.,p%depth + 1) p2%geninfoFather = p%geninfo nextID = nextID + 1 ! DISCRETE VARIABLE ! ----------------- ! else if(variableType .eq. DISCRETE) then ! TODO end if end subroutine branch ! ************************************************************************ ! ************************************************************************ subroutine solveProblem(p,n,x,l,u,m,lambda,equatn,linear,rho,constf,& epsfeas,epsopk,maxit,iter,f,g,gpsupn,cnorm,cnormu,geninfo,inform,& numProblemsSolved) use OptSolver use Problems implicit none ! SCALAR ARGUMENTS logical, intent(in) :: constf integer, intent(in) :: m,n,maxit integer, intent(out) :: geninfo,iter integer, intent(in out) :: inform real(kind=8), intent(in) :: epsfeas,epsopk real(kind=8), intent(out) :: f,gpsupn real(kind=8), intent(in out) :: cnorm,cnormu integer, intent(in out) :: numProblemsSolved ! ARRAY ARGUMENTS logical, intent(in) :: equatn(m),linear(m) real(kind=8), intent(in) :: lambda(m),rho(m) real(kind=8), intent(in) :: l(n),u(n) real(kind=8), intent(out) :: g(n) real(kind=8), intent(in out) :: x(n) ! OTHER ARGUMENTS type(problem), intent(in out) :: p ! Call the solver for the relaxation of this subproblem. if(.not. p%solved) then call optimizationSolver(n,x,l,u,m,lambda,equatn,linear,rho,constf,& epsfeas,epsopk,maxit,iter,f,g,gpsupn,cnorm,cnormu,geninfo,inform) numProblemsSolved = numProblemsSolved + 1 else ! This subproblem has already been solved. x = p%initialPoint f = p%solutionValue geninfo = p%geninfo end if end subroutine solveProblem ! ************************************************************************ ! ************************************************************************ subroutine round(n,x,round_x,domain_) use Problems implicit none ! SCALAR ARGUMENT integer, intent(in) :: n ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n) real(kind=8), intent(out) :: round_x(n) type(domain), intent(in) :: domain_ ! LOCAL SCALARS integer :: i, index real(kind=8) :: floor_x, ceil_x round_x = x do i = 1, size(domain_%type) index = domain_%indices(i) if(domain_%type(i) .eq. BINARY .or. domain_%type(i) .eq. INTEGER) then floor_x = real(floor(x(index)), kind=8) ceil_x = real(ceiling(x(index)), kind=8) ! Round to the closest value. if(x(index) - floor_x .lt. ceil_x - x(index)) then round_x(index) = floor_x else round_x(index) = ceil_x end if end if end do end subroutine round ! ************************************************************************ ! ************************************************************************ subroutine printSolution(objectiveValue,x,n,iterations,time,solutionFound,& numNodesGenerated,numProblemsSolved) implicit none ! SCALAR ARGUMENTS integer, intent(in) :: iterations,n real(kind=8), intent(in) :: objectiveValue real, intent(in) :: time logical, intent(in) :: solutionFound integer, intent(in) :: numNodesGenerated,numProblemsSolved ! ARRAY ARGUMENT real(kind=8) :: x(n) ! LOCAL SCALAR integer :: i if(PRINT_) then if(solutionFound) then print *,'\nBest solution found:\nx =' do i = 1,n write(*,300) x(i) end do write(*,400) objectiveValue,iterations,time else print *,'\nSolution not found.' write(*,500) iterations,time end if write(*,600) numNodesGenerated,numProblemsSolved end if 300 format(' ', 1PE30.20) 400 format(/, 'Objective function value = ', 0PF30.20, & /, 'Number of iterations = ', I14, & /, 'Runtime = ', 0PF30.20) 500 format(/, 'Number of iterations = ', I14, & /, 'Runtime = ', 0PF30.20) 600 format(/, 'Nodes generated = ', I14, & /, 'Problems solved = ', I14) end subroutine printSolution ! ************************************************************************ ! ************************************************************************ subroutine initialize(n,numBinary,numInteger,domain_,discrete_sets,& variableTypes) use Problems use Structures implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n,numBinary,numInteger ! ARRAY ARGUMENTS integer, intent(in) :: variableTypes(n) type(set), allocatable, intent(out) :: discrete_sets(:) ! OTHER ARGUMENTS type(domain), intent(out) :: domain_ ! LOCAL SCALARS integer :: i,count ! Allocate memory for the domain structure. allocate(domain_%indices(numBinary + numInteger)) allocate(domain_%type (numBinary + numInteger)) allocate(domain_%value (numBinary + numInteger)) allocate(discrete_sets(0)) count = 1 do i = 1,n if(variableTypes(i) .eq. BINARY) then domain_%type(count) = BINARY domain_%indices(count) = i count = count + 1 else if(variableTypes(i) .eq. INTEGER) then domain_%type(count) = INTEGER domain_%indices(count) = i count = count + 1 end if end do end subroutine initialize ! ************************************************************************ ! ************************************************************************ logical function satisfyConstraints(n, x, lower, upper) implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n), lower(n), upper(n) satisfyConstraints = satisfyBounds(n, x, lower, upper) end function satisfyConstraints ! ************************************************************************ ! ************************************************************************ logical function satisfyBounds(n, x, lower, upper) implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n), lower(n), upper(n) ! LOCAL SCALARS integer :: i do i = 1, n if(x(i) .lt. lower(i) - EPSILON_FEAS .or. & x(i) .gt. upper(i) + EPSILON_FEAS) then satisfyBounds = .false. return end if end do satisfyBounds = .true. end function satisfyBounds ! ************************************************************************ ! ************************************************************************ logical function satisfyDomain(n, x, domain_, discrete_sets) use Problems use Structures use Auxiliar implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n) type(set), intent(in) :: discrete_sets(:) ! OTHER ARGUMENTS type(domain), intent(in) :: domain_ integer :: i, j, k, index logical :: satisfyDiscrete satisfyDomain = .true. do i = 1, size(domain_%type) index = domain_%indices(i) ! INTEGER variable. if(domain_%type(i) .eq. INTEGER) then if(.not. isInteger(x(index))) then satisfyDomain = .false. return end if ! BINARY variable. else if(domain_%type(i) .eq. BINARY) then if(.not. isBinary(x(index))) then satisfyDomain = .false. return end if ! DISCRETE variable. else if(domain_%type(i) .eq. DISCRETE) then satisfyDiscrete = .false. ! What is the set of variable x_index? do j = 1, size(discrete_sets) if(discrete_sets(j)%index .eq. index) then do k = 1, size(discrete_sets(j)%elements) if(abs(x(index) - discrete_sets(j)%elements(k)) & .le. EPSILON_INT) then satisfyDiscrete = .true. exit end if end do if(.not. satisfyDiscrete) then satisfyDomain = .false. return end if end if end do end if end do end function satisfyDomain ! ************************************************************************ ! ************************************************************************ subroutine addBounds(domain_, lower, upper, n) use Problems implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n ! ARRAY ARGUMENTS real(kind=8), intent(in out) :: lower(n), upper(n) ! OTHER ARGUMENTS type(domain), intent(in) :: domain_ ! LOCAL SCALARS integer :: i, index do i = 1, size(domain_%type) index = domain_%indices(i) if(domain_%type(i) .eq. EQUAL) then lower(index) = domain_%value(i) upper(index) = lower(index) else if(domain_%type(i) .eq. LESS_EQUAL) then upper(index) = domain_%value(i) else if(domain_%type(i) .eq. GREATER_EQUAL) then lower(index) = domain_%value(i) end if end do end subroutine addBounds ! ************************************************************************ ! ************************************************************************ subroutine getSeed(seed, lower, upper, n) use Problems implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n double precision, intent(out) :: seed ! ARRAY ARGUMENTS real(kind=8), intent(in) :: lower(n), upper(n) ! LOCAL SCALARS integer :: i seed = 0.0d0 do i = 1, n seed = seed + real(i,kind=8) * (abs(lower(i)) + abs(upper(i))) end do seed = real(floor(seed),kind=8) if(seed .lt. 0.0d0 .or. seed .gt. 2.0d0**31.0d0 - 1.0d0) then seed = 84543787.0d0 end if end subroutine getSeed ! ************************************************************************ ! ************************************************************************ logical function isNumber(number) implicit none ! This function verifies whether the specified number is valid, ! that is, if the number is not +infinity, -infinity and NaN. ! SCALAR ARGUMENT real(kind=8), intent(in) :: number ! Parameters of the subroutine: ! ============================= ! ! On Entry: ! ========= ! ! NUMBER real(kind=8): The number to be verified. ! ------------------- ! ! Return: ! ======= ! ! Return false if the number is either +infinity, -infinity or NaN. if(.not. (abs(number) .le. LARGEST_VALUE)) then isNumber = .false. else isNumber = .true. end if end function isNumber ! ************************************************************************ ! ************************************************************************ end module MINLP