From 3c9388ffb2cc7d1e9ad630e3f753e11aca9cf213 Mon Sep 17 00:00:00 2001 From: andreaslundell Date: Fri, 24 Sep 2021 11:34:46 +0300 Subject: [PATCH 1/9] Fix format of variables in problem output --- src/Model/AuxiliaryVariables.cpp | 137 +++++++++++++++++++++++-------- src/Model/Problem.cpp | 5 +- src/Model/Variables.cpp | 2 +- 3 files changed, 110 insertions(+), 34 deletions(-) diff --git a/src/Model/AuxiliaryVariables.cpp b/src/Model/AuxiliaryVariables.cpp index 3896ee8b..8a325b04 100644 --- a/src/Model/AuxiliaryVariables.cpp +++ b/src/Model/AuxiliaryVariables.cpp @@ -10,6 +10,8 @@ #include "AuxiliaryVariables.h" +#include "spdlog/fmt/fmt.h" + namespace SHOT { double AuxiliaryVariable::calculate(const VectorDouble& point) const @@ -44,82 +46,153 @@ Interval AuxiliaryVariable::calculate(const IntervalVector& intervalVector) cons std::ostream& operator<<(std::ostream& stream, AuxiliaryVariablePtr var) { - stream << "[" << var->index << "]:\t"; + std::stringstream type; switch(var->properties.type) { - case E_VariableType::Real: - stream << var->lowerBound << " <= " << var->name << " <= " << var->upperBound; + case(E_VariableType::Real): + type << "C "; break; - case E_VariableType::Binary: - stream << var->name << " in {0,1}"; + case(E_VariableType::Binary): + type << "B "; break; - case E_VariableType::Integer: - if(var->lowerBound == 0.0 && var->upperBound == 1.0) - stream << var->name << " in {0,1}"; - else - stream << var->name << " in {" << var->lowerBound << ",...," << var->upperBound << "}"; + case(E_VariableType::Integer): + type << "I "; break; - case E_VariableType::Semicontinuous: - if( var->semiBound < 0.0 ) - stream << var->name << " in {0} or " << var->lowerBound << " <= " << var->name << " <= " << var->semiBound; - else - stream << var->name << " in {0} or " << var->semiBound << " <= " << var->name << " <= " << var->upperBound; + case(E_VariableType::Semicontinuous): + type << "SC"; break; - case E_VariableType::Semiinteger: - if( var->semiBound < 0.0 ) - stream << var->name << " in {" << var->lowerBound << ",...," << var->semiBound << ",0}"; - else - stream << var->name << " in {0," << var->semiBound << ",...," << var->upperBound << "}"; + case(E_VariableType::Semiinteger): + type << "SI"; break; default: - stream << var->lowerBound << " <= " << var->name << " <= " << var->upperBound; + type << "? "; break; } + std::stringstream contains; + + if(var->properties.inObjectiveFunction) + contains << "O"; + else + contains << " "; + + if(var->properties.inLinearConstraints) + contains << "L"; + else + contains << " "; + + if(var->properties.inQuadraticConstraints) + contains << "Q"; + else + contains << " "; + + if(var->properties.inNonlinearConstraints) + contains << "N"; + else + contains << " "; + + std::stringstream inTerms; + + if(var->properties.inLinearTerms) + inTerms << "L"; + else + inTerms << " "; + + if(var->properties.inQuadraticTerms) + inTerms << "Q"; + else + inTerms << " "; + + if(var->properties.inMonomialTerms) + inTerms << "M"; + else + inTerms << " "; + + if(var->properties.inSignomialTerms) + inTerms << "S"; + else + inTerms << " "; + + if(var->properties.inNonlinearExpression) + inTerms << "N"; + else + inTerms << " "; + + std::stringstream auxtype; + switch(var->properties.auxiliaryType) { case E_AuxiliaryVariableType::NonlinearObjectiveFunction: - stream << " (objective auxiliary variable)"; + auxtype << "nonlinear obj. aux. var."; break; case E_AuxiliaryVariableType::NonlinearExpressionPartitioning: - stream << " (partition reformulation for nonlinear sum)"; + auxtype << "nonlinear sum part."; break; case E_AuxiliaryVariableType::MonomialTermsPartitioning: - stream << " (partition reformulation for monomial sum)"; + auxtype << "monomial sum part."; break; case E_AuxiliaryVariableType::SignomialTermsPartitioning: - stream << " (partition reformulation for signomial sum)"; + auxtype << "signomial sum part."; + break; + + case E_AuxiliaryVariableType::SquareTermsPartitioning: + auxtype << "square terms part."; + break; + + case E_AuxiliaryVariableType::ContinuousBilinear: + auxtype << "cont. bilinear lin."; break; case E_AuxiliaryVariableType::BinaryBilinear: - stream << " (binary bilinear linearization)"; + auxtype << "bin bilinear lin."; break; + case E_AuxiliaryVariableType::BinaryContinuousBilinear: - stream << " (mixed binary-continuous bilinear linearization)"; + auxtype << "mixed bin./cont. bilinear lin."; break; - case E_AuxiliaryVariableType::ContinuousBilinear: - stream << " (continuous bilinear linearization)"; + case E_AuxiliaryVariableType::IntegerBilinear: + auxtype << "int. bilinear lin."; break; - case E_AuxiliaryVariableType::IntegerBilinear: - stream << " (integer bilinear linearization)"; + case E_AuxiliaryVariableType::BinaryMonomial: + auxtype << "bin. monomial lin."; + break; + + case E_AuxiliaryVariableType::AbsoluteValue: + auxtype << "abs. value ref."; + break; + + case E_AuxiliaryVariableType::AntiEpigraph: + auxtype << "anti epigraph ref."; + break; + + case E_AuxiliaryVariableType::EigenvalueDecomposition: + auxtype << "eigenval. decomp. ref."; break; default: - stream << " (unknown auxiliary variable)"; + auxtype << "unspecified aux. var."; break; } + stream << fmt::format("[{:>6d},{:<1s}] [{:<4s}] [{:<5s}]\t{:>12f} {:1s} <= {:^16s} <= {:1s} {:<12f} {:>30s}", + var->index, type.str(), contains.str(), inTerms.str(), + (var->properties.type == E_VariableType::Semicontinuous || var->properties.type == E_VariableType::Semiinteger) + ? var->semiBound + : var->lowerBound, + var->properties.hasLowerBoundBeenTightened ? "*" : " ", var->name, + var->properties.hasUpperBoundBeenTightened ? "*" : " ", var->upperBound, auxtype.str()); + return stream; } diff --git a/src/Model/Problem.cpp b/src/Model/Problem.cpp index 75267368..5a4eeab9 100644 --- a/src/Model/Problem.cpp +++ b/src/Model/Problem.cpp @@ -2223,7 +2223,10 @@ std::ostream& operator<<(std::ostream& stream, const Problem& problem) for(auto& V : problem.allVariables) { - stream << V << '\n'; + if(!V->properties.isAuxiliary) + stream << V << '\n'; + else + stream << std::static_pointer_cast(V) << '\n'; } switch(problem.properties.convexity) diff --git a/src/Model/Variables.cpp b/src/Model/Variables.cpp index c2d8c304..cbfea973 100644 --- a/src/Model/Variables.cpp +++ b/src/Model/Variables.cpp @@ -204,7 +204,7 @@ std::ostream& operator<<(std::ostream& stream, VariablePtr var) else inTerms << " "; - stream << fmt::format("[{:>6d},{:<1s}] [{:<4s}] [{:<5s}]\t{:>10} {:1s} <= {:^14s} <= {:1s} {:<10}", var->index, + stream << fmt::format("[{:>6d},{:<1s}] [{:<4s}] [{:<5s}]\t{:>12f} {:1s} <= {:^16s} <= {:1s} {:<12f}", var->index, type.str(), contains.str(), inTerms.str(), (var->properties.type == E_VariableType::Semicontinuous || var->properties.type == E_VariableType::Semiinteger) ? var->semiBound From 9e37916a9933c19613cf61f5077b2499bcf6ceba Mon Sep 17 00:00:00 2001 From: andreaslundell Date: Fri, 24 Sep 2021 12:16:47 +0300 Subject: [PATCH 2/9] Add assert that checks variables are added in correct order --- src/Model/Problem.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/Model/Problem.cpp b/src/Model/Problem.cpp index 5a4eeab9..0356678b 100644 --- a/src/Model/Problem.cpp +++ b/src/Model/Problem.cpp @@ -892,6 +892,8 @@ void Problem::add(VariablePtr variable) break; } + assert(variable->index + 1 == allVariables.size()); + variable->takeOwnership(shared_from_this()); variablesUpdated = false; @@ -934,6 +936,8 @@ void Problem::add(AuxiliaryVariablePtr variable) break; } + assert(variable->index + 1 == allVariables.size()); + variable->takeOwnership(shared_from_this()); variablesUpdated = false; From 0435e572025f7f9f468829989aff7555b586c724 Mon Sep 17 00:00:00 2001 From: andreaslundell Date: Fri, 24 Sep 2021 12:17:10 +0300 Subject: [PATCH 3/9] Fix eigenvalue decomposition, add option to move coefficient to nonlinear constraint --- src/Enums.h | 6 ++++++ src/Solver.cpp | 18 +++++++++++++++--- src/Tasks/TaskReformulateProblem.cpp | 26 ++++++++++++++++++++------ 3 files changed, 41 insertions(+), 9 deletions(-) diff --git a/src/Enums.h b/src/Enums.h index 19296e52..3f6854d1 100644 --- a/src/Enums.h +++ b/src/Enums.h @@ -272,6 +272,12 @@ enum class ES_AddPrimalPointAsInteriorPoint OnlyAverage }; +enum class ES_EigenValueDecompositionFormulation +{ + CoefficientReformulated, + CoefficientRemains +}; + enum class ES_HyperplaneCutStrategy { ESH, diff --git a/src/Solver.cpp b/src/Solver.cpp index 939c04b2..c6304d7b 100644 --- a/src/Solver.cpp +++ b/src/Solver.cpp @@ -1086,7 +1086,19 @@ void Solver::initializeSettings() "How to treat quadratic functions", enumQPStrategy, 0); enumQPStrategy.clear(); - env->settings->createSetting("Reformulation.Quadratics.UseEigenValueDecomposition", "Model", false, + env->settings->createSetting("Reformulation.Quadratics.EigenValueDecomposition.Use", "Model", false, + "Whether to use the eigenvalue decomposition of convex quadratic functions"); + + VectorString enumEigenValueStrategy; + enumEigenValueStrategy.push_back("Term coefficient is included in reformulation"); + enumEigenValueStrategy.push_back("Term coefficient remains"); + + env->settings->createSetting("Reformulation.Quadratics.EigenValueDecomposition.Formulation", "Model", + static_cast(ES_EigenValueDecompositionFormulation::CoefficientReformulated), + "Which formulation to use in eigenvalue decomposition", enumEigenValueStrategy, 0); + enumEigenValueStrategy.clear(); + + env->settings->createSetting("Reformulation.Quadratics.EigenValueDecomposition.Method", "Model", false, "Whether to use the eigen value decomposition of convex quadratic functions"); // Modeling system settings @@ -1825,7 +1837,7 @@ void Solver::setConvexityBasedSettingsPreReformulation() #ifdef HAS_CBC if(static_cast(env->settings->getSetting("MIP.Solver", "Dual")) == ES_MIPSolver::Cbc) { - env->settings->updateSetting("Reformulation.Quadratics.UseEigenValueDecomposition", "Model", false); + env->settings->updateSetting("Reformulation.Quadratics.EigenValueDecomposition.Use", "Model", false); } #endif } @@ -1834,7 +1846,7 @@ void Solver::setConvexityBasedSettingsPreReformulation() #ifdef HAS_CBC if(static_cast(env->settings->getSetting("MIP.Solver", "Dual")) == ES_MIPSolver::Cbc) { - env->settings->updateSetting("Reformulation.Quadratics.UseEigenValueDecomposition", "Model", true); + env->settings->updateSetting("Reformulation.Quadratics.EigenValueDecomposition.Use", "Model", true); } #endif } diff --git a/src/Tasks/TaskReformulateProblem.cpp b/src/Tasks/TaskReformulateProblem.cpp index 0b60a6fd..7169fba1 100644 --- a/src/Tasks/TaskReformulateProblem.cpp +++ b/src/Tasks/TaskReformulateProblem.cpp @@ -782,6 +782,7 @@ void TaskReformulateProblem::createEpigraphConstraint() auto objectiveVariable = std::make_shared( "shot_objvar", auxVariableCounter, E_VariableType::Real, objectiveBound.l(), objectiveBound.u()); + auxVariableCounter++; objectiveVariable->properties.auxiliaryType = E_AuxiliaryVariableType::NonlinearObjectiveFunction; env->results->increaseAuxiliaryVariableCounter(E_AuxiliaryVariableType::NonlinearObjectiveFunction); @@ -1765,7 +1766,7 @@ std::tuple TaskReformulateProblem::reformulateAndPa } } - if(env->settings->getSetting("Reformulation.Quadratics.UseEigenValueDecomposition", "Model") + if(env->settings->getSetting("Reformulation.Quadratics.EigenValueDecomposition.Use", "Model") && partitionStrategy <= ES_PartitionNonlinearSums::IfConvex && quadraticSumConvex && !quadraticTerms.allSquares) // Use the eigenvalue decomposition reformulation { @@ -2253,15 +2254,28 @@ LinearTerms TaskReformulateProblem::doEigenvalueDecomposition(QuadraticTerms qua auto auxQuadVariable = std::make_shared("q_evd_" + std::to_string(auxVariableCounter), auxVariableCounter, E_VariableType::Real, bounds.l(), bounds.u()); auxVariableCounter++; + auxQuadVariable->properties.auxiliaryType = E_AuxiliaryVariableType::EigenvalueDecomposition; + reformulatedProblem->add(auxQuadVariable); + env->results->increaseAuxiliaryVariableCounter(E_AuxiliaryVariableType::EigenvalueDecomposition); - auto [auxVariable, newVariable] = getSquareAuxiliaryVariable( - auxQuadVariable, quadraticTerms.eigenvalues[i].real(), E_AuxiliaryVariableType::EigenvalueDecomposition); - resultLinearTerms.add(std::make_shared(0.5, auxVariable)); + if(env->settings->getSetting("Reformulation.Quadratics.EigenValueDecomposition.Formulation", "Model") + == static_cast(ES_EigenValueDecompositionFormulation::CoefficientReformulated)) + { + auto [auxVariable, newVariable] = getSquareAuxiliaryVariable(auxQuadVariable, + quadraticTerms.eigenvalues[i].real(), E_AuxiliaryVariableType::EigenvalueDecomposition); + resultLinearTerms.add(std::make_shared(0.5, auxVariable)); + } + else + { + auto [auxVariable, newVariable] + = getSquareAuxiliaryVariable(auxQuadVariable, 1.0, E_AuxiliaryVariableType::EigenvalueDecomposition); + resultLinearTerms.add( + std::make_shared(0.5 * quadraticTerms.eigenvalues[i].real(), auxVariable)); + } auxConstraint->add(std::make_shared(-1.0, auxQuadVariable)); - reformulatedProblem->add(auxConstraint); - reformulatedProblem->add(std::move(auxQuadVariable)); + reformulatedProblem->add(std::move(auxConstraint)); } return (resultLinearTerms); From 417522e4258eea7eee9f3dfb7fe292dcae157f97 Mon Sep 17 00:00:00 2001 From: andreaslundell Date: Fri, 24 Sep 2021 12:30:47 +0300 Subject: [PATCH 4/9] Fix missing output of number of square refs --- src/Report.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/Report.cpp b/src/Report.cpp index fc4a4cc5..b37b69e5 100644 --- a/src/Report.cpp +++ b/src/Report.cpp @@ -878,6 +878,10 @@ void Report::outputProblemInstanceReport() value > 0) env->output->outputInfo(fmt::format(" {:56s}{:d}", " - signomial terms partitioning:", value)); + if(auto value = env->results->getAuxiliaryVariableCounter(E_AuxiliaryVariableType::SquareTermsPartitioning); + value > 0) + env->output->outputInfo(fmt::format(" {:56s}{:d}", " - square terms partitioning:", value)); + if(auto value = env->results->getAuxiliaryVariableCounter(E_AuxiliaryVariableType::ContinuousBilinear); value > 0) env->output->outputInfo(fmt::format(" {:56s}{:d}", " - continuous bilinear term extraction:", value)); From c17079af6e48035d3983184b06d2cc2c7b9f31b1 Mon Sep 17 00:00:00 2001 From: andreaslundell Date: Wed, 29 Sep 2021 16:39:16 +0300 Subject: [PATCH 5/9] Improve fixed NLP calls in single-tree strategy --- src/MIPSolver/MIPSolverCplexSingleTree.cpp | 27 ++++--- .../MIPSolverCplexSingleTreeLegacy.cpp | 24 +++--- src/MIPSolver/MIPSolverGurobiSingleTree.cpp | 77 ++++++++++--------- 3 files changed, 71 insertions(+), 57 deletions(-) diff --git a/src/MIPSolver/MIPSolverCplexSingleTree.cpp b/src/MIPSolver/MIPSolverCplexSingleTree.cpp index fc7c327e..80aea014 100644 --- a/src/MIPSolver/MIPSolverCplexSingleTree.cpp +++ b/src/MIPSolver/MIPSolverCplexSingleTree.cpp @@ -315,21 +315,26 @@ void CplexCallback::invoke(const IloCplex::Callback::Context& context) if(checkFixedNLPStrategy(candidatePoints.at(0))) { - env->primalSolver->addFixedNLPCandidate(candidatePoints.at(0).point, E_PrimalNLPSource::FirstSolution, - context.getCandidateObjective(), env->results->getCurrentIteration()->iterationNumber, - candidatePoints.at(0).maxDeviation); if(taskSelectPrimNLPOriginal) - taskSelectPrimNLPOriginal->run(); + { + env->primalSolver->addFixedNLPCandidate(candidatePoints.at(0).point, + E_PrimalNLPSource::FirstSolution, context.getCandidateObjective(), + env->results->getCurrentIteration()->iterationNumber, candidatePoints.at(0).maxDeviation); - env->primalSolver->addFixedNLPCandidate(candidatePoints.at(0).point, E_PrimalNLPSource::FirstSolution, - context.getCandidateObjective(), env->results->getCurrentIteration()->iterationNumber, - candidatePoints.at(0).maxDeviation); + taskSelectPrimNLPOriginal->run(); + env->primalSolver->fixedPrimalNLPCandidates.clear(); + } if(taskSelectPrimNLPReformulated) - taskSelectPrimNLPReformulated->run(); + { + env->primalSolver->addFixedNLPCandidate(candidatePoints.at(0).point, + E_PrimalNLPSource::FirstSolution, context.getCandidateObjective(), + env->results->getCurrentIteration()->iterationNumber, candidatePoints.at(0).maxDeviation); - env->primalSolver->fixedPrimalNLPCandidates.clear(); + taskSelectPrimNLPReformulated->run(); + env->primalSolver->fixedPrimalNLPCandidates.clear(); + } env->primalSolver->checkPrimalSolutionCandidates(); } @@ -554,8 +559,8 @@ void CplexCallback::addLazyConstraint( for(auto& hp : env->dualSolver->hyperplaneWaitingList) { - this->createHyperplane(hp, context); - this->lastNumAddedHyperplanes++; + if(this->createHyperplane(hp, context)) + this->lastNumAddedHyperplanes++; } env->dualSolver->hyperplaneWaitingList.clear(); diff --git a/src/MIPSolver/MIPSolverCplexSingleTreeLegacy.cpp b/src/MIPSolver/MIPSolverCplexSingleTreeLegacy.cpp index e4c00b26..0f5871bb 100644 --- a/src/MIPSolver/MIPSolverCplexSingleTreeLegacy.cpp +++ b/src/MIPSolver/MIPSolverCplexSingleTreeLegacy.cpp @@ -428,19 +428,23 @@ void CtCallbackI::main() if(checkFixedNLPStrategy(solutionCandidate)) { - env->primalSolver->addFixedNLPCandidate(solution, E_PrimalNLPSource::FirstSolution, this->getObjValue(), - env->results->getCurrentIteration()->iterationNumber, solutionCandidate.maxDeviation); - if(taskSelectPrimNLPOriginal) - taskSelectPrimNLPOriginal->run(); + { + env->primalSolver->addFixedNLPCandidate(solution, E_PrimalNLPSource::FirstSolution, this->getObjValue(), + env->results->getCurrentIteration()->iterationNumber, solutionCandidate.maxDeviation); - env->primalSolver->addFixedNLPCandidate(solution, E_PrimalNLPSource::FirstSolution, this->getObjValue(), - env->results->getCurrentIteration()->iterationNumber, solutionCandidate.maxDeviation); + taskSelectPrimNLPOriginal->run(); + env->primalSolver->fixedPrimalNLPCandidates.clear(); + } if(taskSelectPrimNLPReformulated) - taskSelectPrimNLPReformulated->run(); + { + env->primalSolver->addFixedNLPCandidate(solution, E_PrimalNLPSource::FirstSolution, this->getObjValue(), + env->results->getCurrentIteration()->iterationNumber, solutionCandidate.maxDeviation); - env->primalSolver->fixedPrimalNLPCandidates.clear(); + taskSelectPrimNLPReformulated->run(); + env->primalSolver->fixedPrimalNLPCandidates.clear(); + } env->primalSolver->checkPrimalSolutionCandidates(); @@ -485,8 +489,8 @@ void CtCallbackI::main() for(auto& hp : env->dualSolver->hyperplaneWaitingList) { - this->createHyperplane(hp); - this->lastNumAddedHyperplanes++; + if(this->createHyperplane(hp)) + this->lastNumAddedHyperplanes++; } env->dualSolver->hyperplaneWaitingList.clear(); diff --git a/src/MIPSolver/MIPSolverGurobiSingleTree.cpp b/src/MIPSolver/MIPSolverGurobiSingleTree.cpp index c5f042ad..f5c11d31 100644 --- a/src/MIPSolver/MIPSolverGurobiSingleTree.cpp +++ b/src/MIPSolver/MIPSolverGurobiSingleTree.cpp @@ -186,6 +186,31 @@ void GurobiCallbackSingleTree::callback() try { + + // Add current primal bound as new incumbent candidate + auto primalBound = env->results->getPrimalBound(); + + if(((isMinimization && lastUpdatedPrimal < primalBound) + || (!isMinimization && lastUpdatedPrimal > primalBound))) + { + auto primalSol = env->results->primalSolution; + + if((int)primalSol.size() < env->reformulatedProblem->properties.numberOfVariables) + env->reformulatedProblem->augmentAuxiliaryVariableValues(primalSol); + + assert(env->reformulatedProblem->properties.numberOfVariables == primalSol.size()); + + if(env->dualSolver->MIPSolver->hasDualAuxiliaryObjectiveVariable()) + primalSol.push_back(env->reformulatedProblem->objectiveFunction->calculateValue(primalSol)); + + for(size_t i = 0; i < primalSol.size(); i++) + { + setSolution(vars[i], primalSol.at(i)); + } + + lastUpdatedPrimal = env->results->getPrimalBound(); + } + // Check if better dual bound double tmpDualObjBound; @@ -394,21 +419,25 @@ void GurobiCallbackSingleTree::callback() if(checkFixedNLPStrategy(candidatePoints.at(0))) { - env->primalSolver->addFixedNLPCandidate(candidatePoints.at(0).point, E_PrimalNLPSource::FirstSolution, - getDoubleInfo(GRB_CB_MIPSOL_OBJ), env->results->getCurrentIteration()->iterationNumber, - candidatePoints.at(0).maxDeviation); - if(taskSelectPrimNLPOriginal) - taskSelectPrimNLPOriginal->run(); + { + env->primalSolver->addFixedNLPCandidate(candidatePoints.at(0).point, + E_PrimalNLPSource::FirstSolution, getDoubleInfo(GRB_CB_MIPSOL_OBJ), + env->results->getCurrentIteration()->iterationNumber, candidatePoints.at(0).maxDeviation); - env->primalSolver->addFixedNLPCandidate(candidatePoints.at(0).point, E_PrimalNLPSource::FirstSolution, - getDoubleInfo(GRB_CB_MIPSOL_OBJ), env->results->getCurrentIteration()->iterationNumber, - candidatePoints.at(0).maxDeviation); + taskSelectPrimNLPOriginal->run(); + env->primalSolver->fixedPrimalNLPCandidates.clear(); + } if(taskSelectPrimNLPReformulated) - taskSelectPrimNLPReformulated->run(); + { + env->primalSolver->addFixedNLPCandidate(candidatePoints.at(0).point, + E_PrimalNLPSource::FirstSolution, getDoubleInfo(GRB_CB_MIPSOL_OBJ), + env->results->getCurrentIteration()->iterationNumber, candidatePoints.at(0).maxDeviation); - env->primalSolver->fixedPrimalNLPCandidates.clear(); + taskSelectPrimNLPReformulated->run(); + env->primalSolver->fixedPrimalNLPCandidates.clear(); + } env->primalSolver->checkPrimalSolutionCandidates(); } @@ -443,30 +472,6 @@ void GurobiCallbackSingleTree::callback() lastExploredNodes = (int)getDoubleInfo(GRB_CB_MIP_NODCNT); lastOpenNodes = (int)getDoubleInfo(GRB_CB_MIP_NODLFT); } - - // Add current primal bound as new incumbent candidate - auto primalBound = env->results->getPrimalBound(); - - if(((isMinimization && lastUpdatedPrimal < primalBound) - || (!isMinimization && lastUpdatedPrimal > primalBound))) - { - auto primalSol = env->results->primalSolution; - - if((int)primalSol.size() < env->reformulatedProblem->properties.numberOfVariables) - env->reformulatedProblem->augmentAuxiliaryVariableValues(primalSol); - - assert(env->reformulatedProblem->properties.numberOfVariables == primalSol.size()); - - if(env->dualSolver->MIPSolver->hasDualAuxiliaryObjectiveVariable()) - primalSol.push_back(env->reformulatedProblem->objectiveFunction->calculateValue(primalSol)); - - for(size_t i = 0; i < primalSol.size(); i++) - { - setSolution(vars[i], primalSol.at(i)); - } - - lastUpdatedPrimal = env->results->getPrimalBound(); - } } catch(GRBException& e) { @@ -669,8 +674,8 @@ void GurobiCallbackSingleTree::addLazyConstraint(std::vector cand for(auto& hp : env->dualSolver->hyperplaneWaitingList) { - this->createHyperplane(hp); - this->lastNumAddedHyperplanes++; + if(this->createHyperplane(hp)) + this->lastNumAddedHyperplanes++; } env->dualSolver->hyperplaneWaitingList.clear(); From 0c2ae7eee298d0f6906fbe3990348be1cc4aba61 Mon Sep 17 00:00:00 2001 From: andreaslundell Date: Wed, 29 Sep 2021 18:29:04 +0300 Subject: [PATCH 6/9] Fix cutoff functionality --- src/Results.cpp | 21 ++++++++++++++++--- src/Solver.cpp | 13 ++++++++++++ src/Tasks/TaskPerformBoundTightening.cpp | 26 +++++++++--------------- src/Tasks/TaskSolveIteration.cpp | 12 ++--------- 4 files changed, 43 insertions(+), 29 deletions(-) diff --git a/src/Results.cpp b/src/Results.cpp index 07f63a91..aa764032 100644 --- a/src/Results.cpp +++ b/src/Results.cpp @@ -1082,7 +1082,7 @@ int Results::getNumberOfIterations() { return (iterations.size()); } double Results::getPrimalBound() { - if(this->currentPrimalBound != NAN) + if(!isnan(this->currentPrimalBound)) return (this->currentPrimalBound); else if(env->problem->objectiveFunction->direction == E_ObjectiveFunctionDirection::Minimize) return (SHOT_DBL_MAX); @@ -1112,8 +1112,23 @@ void Results::setPrimalBound(double value) this->currentDualBound = value; } - env->dualSolver->cutOffToUse = value; - env->dualSolver->useCutOff = true; + if(env->problem->objectiveFunction->properties.isMinimize) + { + if(value < env->dualSolver->cutOffToUse) + { + env->dualSolver->cutOffToUse = value; + env->dualSolver->useCutOff = true; + } + } + else + { + if(value > env->dualSolver->cutOffToUse) + { + env->dualSolver->cutOffToUse = value; + env->dualSolver->useCutOff = true; + } + } + env->solutionStatistics.numberOfIterationsWithPrimalStagnation = 0; env->solutionStatistics.lastIterationWithSignificantPrimalUpdate = getNumberOfIterations() - 1; env->solutionStatistics.numberOfPrimalReductionCutsUpdatesWithoutEffect = 0; diff --git a/src/Solver.cpp b/src/Solver.cpp index c6304d7b..9572ebee 100644 --- a/src/Solver.cpp +++ b/src/Solver.cpp @@ -356,6 +356,19 @@ bool Solver::setProblem(std::string fileName) "───────────────────────────────────────────────────────────────────────────────────────────────────╴"); #endif + if(env->settings->getSetting("MIP.CutOff.UseInitialValue", "Dual") + && std::abs(env->settings->getSetting("MIP.CutOff.InitialValue", "Dual")) < SHOT_DBL_MAX) + { + env->dualSolver->cutOffToUse = env->settings->getSetting("MIP.CutOff.InitialValue", "Dual"); + env->dualSolver->useCutOff = true; + env->output->outputDebug( + fmt::format(" Setting user defined cutoff value to {}.", env->dualSolver->cutOffToUse)); + } + else + { + env->dualSolver->cutOffToUse = env->results->getPrimalBound(); + } + auto taskPerformBoundTightening = std::make_unique(env, env->problem); taskPerformBoundTightening->run(); diff --git a/src/Tasks/TaskPerformBoundTightening.cpp b/src/Tasks/TaskPerformBoundTightening.cpp index bc7ef289..57915adc 100644 --- a/src/Tasks/TaskPerformBoundTightening.cpp +++ b/src/Tasks/TaskPerformBoundTightening.cpp @@ -142,30 +142,24 @@ void TaskPerformBoundTightening::run() if(sourceProblem->objectiveFunction->properties.isMinimize) { - if(objectiveBoundsAfter.l() > -SHOT_DBL_INF) // Update dual bound - { - DualSolution sol = { {}, E_DualSolutionSource::MIPSolverBound, objectiveBoundsAfter.l(), 0, false }; - env->dualSolver->addDualSolutionCandidate(sol); - } + DualSolution sol = { {}, E_DualSolutionSource::MIPSolverBound, objectiveBoundsAfter.l(), 0, false }; + env->dualSolver->addDualSolutionCandidate(sol); - if(objectiveBoundsAfter.u() < SHOT_DBL_INF) // Update MIP cutoff + if(objectiveBoundsAfter.u() < env->dualSolver->cutOffToUse) // Update MIP cutoff { - env->settings->updateSetting("MIP.CutOff.InitialValue", "Dual", objectiveBoundsAfter.u()); - env->settings->updateSetting("MIP.CutOff.UseInitialValue", "Dual", true); + env->dualSolver->cutOffToUse = objectiveBoundsAfter.u(); + env->dualSolver->useCutOff = true; } } else if(sourceProblem->objectiveFunction->properties.isMaximize) { - if(objectiveBoundsAfter.u() < SHOT_DBL_INF) // Update dual bound - { - DualSolution sol = { {}, E_DualSolutionSource::MIPSolverBound, objectiveBoundsAfter.u(), 0, false }; - env->dualSolver->addDualSolutionCandidate(sol); - } + DualSolution sol = { {}, E_DualSolutionSource::MIPSolverBound, objectiveBoundsAfter.u(), 0, false }; + env->dualSolver->addDualSolutionCandidate(sol); - if(objectiveBoundsAfter.l() > -SHOT_DBL_INF) // Update MIP cutoff + if(objectiveBoundsAfter.l() > env->dualSolver->cutOffToUse) // Update MIP cutoff { - env->settings->updateSetting("MIP.CutOff.InitialValue", "Dual", objectiveBoundsAfter.l()); - env->settings->updateSetting("MIP.CutOff.UseInitialValue", "Dual", true); + env->dualSolver->cutOffToUse = objectiveBoundsAfter.l(); + env->dualSolver->useCutOff = true; } } } diff --git a/src/Tasks/TaskSolveIteration.cpp b/src/Tasks/TaskSolveIteration.cpp index 0c4b2ed7..9ab96a7e 100644 --- a/src/Tasks/TaskSolveIteration.cpp +++ b/src/Tasks/TaskSolveIteration.cpp @@ -62,16 +62,6 @@ void TaskSolveIteration::run() double cutOffValue; double cutOffValueConstraint; - // Setting a user provided cutoff value - if(env->dualSolver->cutOffToUse == SHOT_DBL_MAX - && env->settings->getSetting("MIP.CutOff.UseInitialValue", "Dual")) - { - env->dualSolver->useCutOff = true; - env->dualSolver->cutOffToUse = env->settings->getSetting("MIP.CutOff.InitialValue", "Dual"); - env->output->outputDebug( - fmt::format(" Setting user-provided cutoff value to {}.", env->dualSolver->cutOffToUse)); - } - if(isMinimization) { cutOffValue @@ -85,6 +75,8 @@ void TaskSolveIteration::run() cutOffValueConstraint = env->dualSolver->cutOffToUse; } + env->output->outputDebug(fmt::format(" Setting cutoff value to {}.", env->dualSolver->cutOffToUse)); + env->dualSolver->MIPSolver->setCutOff(cutOffValue); if(env->reformulatedProblem->objectiveFunction->properties.classification From 291560245e276961c7eff3e252e6b149dd6518e5 Mon Sep 17 00:00:00 2001 From: andreaslundell Date: Wed, 29 Sep 2021 18:39:24 +0300 Subject: [PATCH 7/9] Gurobi does now like large cutoffs for some problems, e.g. fo7_2 --- src/MIPSolver/MIPSolverGurobi.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/MIPSolver/MIPSolverGurobi.cpp b/src/MIPSolver/MIPSolverGurobi.cpp index 70846f58..f27718ce 100644 --- a/src/MIPSolver/MIPSolverGurobi.cpp +++ b/src/MIPSolver/MIPSolverGurobi.cpp @@ -1096,6 +1096,9 @@ void MIPSolverGurobi::setTimeLimit(double seconds) void MIPSolverGurobi::setCutOff(double cutOff) { + if(std::abs(cutOff) > 1e20) + return; + try { // Gurobi has problems if not an epsilon value is added to the cutoff... From 0697744a9e38ab5fae0534a55c08cee98d3009bd Mon Sep 17 00:00:00 2001 From: andreaslundell Date: Thu, 30 Sep 2021 15:36:49 +0300 Subject: [PATCH 8/9] Fix seg fault --- src/DualSolver.cpp | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/DualSolver.cpp b/src/DualSolver.cpp index 7a02d32b..f9696e9d 100644 --- a/src/DualSolver.cpp +++ b/src/DualSolver.cpp @@ -72,7 +72,13 @@ void DualSolver::checkDualSolutionCandidates() // New dual solution env->results->setDualBound(C.objValue); currDualBound = C.objValue; - env->solutionStatistics.iterationLastDualBoundUpdate = env->results->getCurrentIteration()->iterationNumber; + + if(env->results->getNumberOfIterations() > 0) + env->solutionStatistics.iterationLastDualBoundUpdate + = env->results->getCurrentIteration()->iterationNumber; + else + env->solutionStatistics.iterationLastDualBoundUpdate = 0; + env->solutionStatistics.iterationLastDualBoundUpdate = env->timing->getElapsedTime("Total"); if(C.sourceType == E_DualSolutionSource::MIPSolutionOptimal @@ -257,7 +263,8 @@ bool DualSolver::hasHyperplaneBeenAdded(double hash, int constraintIndex) void DualSolver::addIntegerCut(IntegerCut integerCut) { - if(env->reformulatedProblem->properties.numberOfIntegerVariables > 0 || env->reformulatedProblem->properties.numberOfSemiintegerVariables > 0) + if(env->reformulatedProblem->properties.numberOfIntegerVariables > 0 + || env->reformulatedProblem->properties.numberOfSemiintegerVariables > 0) integerCut.areAllVariablesBinary = false; else { From 48ed866a95e6b9705759331e86d9c3a11d908b1d Mon Sep 17 00:00:00 2001 From: andreaslundell Date: Wed, 6 Oct 2021 19:14:04 +0300 Subject: [PATCH 9/9] Change default value for numeric focus with Gurobi --- src/Solver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Solver.cpp b/src/Solver.cpp index 9572ebee..d028b2e9 100644 --- a/src/Solver.cpp +++ b/src/Solver.cpp @@ -1415,7 +1415,7 @@ void Solver::initializeSettings() enumGurobiNumericFocus.push_back("Mild"); enumGurobiNumericFocus.push_back("Moderate"); enumGurobiNumericFocus.push_back("Aggressive"); - env->settings->createSetting("Gurobi.NumericFocus", "Subsolver", 2, "MIP focus", enumGurobiNumericFocus, 0); + env->settings->createSetting("Gurobi.NumericFocus", "Subsolver", 1, "MIP focus", enumGurobiNumericFocus, 0); enumGurobiNumericFocus.clear(); VectorString enumGurobiPoolSearchMode;