From fe6276284374797c4593b4c7b75110cd374fcdc3 Mon Sep 17 00:00:00 2001 From: andreaslundell Date: Wed, 9 Oct 2024 13:51:41 +0300 Subject: [PATCH] Verify dual bound --- src/DualSolver.cpp | 27 ++- src/Enums.h | 6 +- src/MIPSolver/MIPSolverCplex.cpp | 4 +- src/PrimalSolver.cpp | 5 +- src/Results.cpp | 18 +- src/Results.h | 2 +- .../SolutionStrategyMIQCQP.cpp | 2 + .../SolutionStrategyMultiTree.cpp | 8 + src/Solver.cpp | 3 + src/Structs.h | 12 +- src/Tasks/TaskCreateDualProblem.cpp | 191 +-------------- src/Tasks/TaskCreateDualProblem.h | 6 +- src/Tasks/TaskCreateMIPProblem.cpp | 220 ++++++++++++++++++ src/Tasks/TaskCreateMIPProblem.h | 32 +++ src/Tasks/TaskPerformDualBounding.cpp | 190 +++++++++++++++ src/Tasks/TaskPerformDualBounding.h | 36 +++ src/Tasks/TaskSolveIteration.cpp | 4 + 17 files changed, 553 insertions(+), 213 deletions(-) create mode 100644 src/Tasks/TaskCreateMIPProblem.cpp create mode 100644 src/Tasks/TaskCreateMIPProblem.h create mode 100644 src/Tasks/TaskPerformDualBounding.cpp create mode 100644 src/Tasks/TaskPerformDualBounding.h diff --git a/src/DualSolver.cpp b/src/DualSolver.cpp index 3732635a..97b4befc 100644 --- a/src/DualSolver.cpp +++ b/src/DualSolver.cpp @@ -69,8 +69,11 @@ void DualSolver::checkDualSolutionCandidates() if(updateDual) { - // New dual solution - env->results->setDualBound(C.objValue); + if(C.sourceType == E_DualSolutionSource::ConvexBounding) + env->results->setDualBound(C.objValue, true); // Force valid dual bound update + else + env->results->setDualBound(C.objValue); + currDualBound = C.objValue; if(env->results->getNumberOfIterations() > 0) @@ -83,7 +86,8 @@ void DualSolver::checkDualSolutionCandidates() if(C.sourceType == E_DualSolutionSource::MIPSolutionOptimal || C.sourceType == E_DualSolutionSource::LPSolution - || C.sourceType == E_DualSolutionSource::MIPSolverBound) + || C.sourceType == E_DualSolutionSource::MIPSolverBound + || C.sourceType == E_DualSolutionSource::ConvexBounding) { env->results->addDualSolution(C); } @@ -104,11 +108,21 @@ void DualSolver::checkDualSolutionCandidates() case E_DualSolutionSource::MIPSolverBound: sourceDesc = "MIP solver bound"; break; + case E_DualSolutionSource::ConvexBounding: + sourceDesc = "Convex MIP bounding"; + break; default: break; } - env->output->outputDebug(fmt::format(" New dual bound {}, source: {}", C.objValue, sourceDesc)); + if(C.sourceType == E_DualSolutionSource::ConvexBounding) + { + env->output->outputInfo(fmt::format(" New dual bound {}, source: {}", C.objValue, sourceDesc)); + } + else + { + env->output->outputDebug(fmt::format(" New dual bound {}, source: {}", C.objValue, sourceDesc)); + } } } @@ -232,6 +246,11 @@ void DualSolver::addGeneratedHyperplane(const Hyperplane& hyperplane) currentIteration->totNumHyperplanes++; env->solutionStatistics.iterationLastDualCutAdded = currentIteration->iterationNumber; + if(hyperplane.isSourceConvex) + env->solutionStatistics.numberOfHyperplanesWithConvexSource++; + else + env->solutionStatistics.numberOfHyperplanesWithNonconvexSource++; + env->output->outputTrace(" Hyperplane generated from: " + source); } diff --git a/src/Enums.h b/src/Enums.h index 3f6854d1..0b4caac8 100644 --- a/src/Enums.h +++ b/src/Enums.h @@ -46,7 +46,8 @@ enum class E_DualSolutionSource LPSolution, MIPSolutionOptimal, ObjectiveConstraint, - MIPSolverBound + MIPSolverBound, + ConvexBounding }; enum class E_EventType @@ -152,7 +153,8 @@ enum class E_PrimalSolutionSource MIPSolutionPool, LPFixedIntegers, MIPCallback, - InteriorPointSearch + InteriorPointSearch, + ConvexBounding }; enum class E_ProblemConvexity diff --git a/src/MIPSolver/MIPSolverCplex.cpp b/src/MIPSolver/MIPSolverCplex.cpp index 16a0574c..5d2028e2 100644 --- a/src/MIPSolver/MIPSolverCplex.cpp +++ b/src/MIPSolver/MIPSolverCplex.cpp @@ -1139,9 +1139,7 @@ void MIPSolverCplex::setTimeLimit(double seconds) { try { - if(seconds > 1e+75) - { - } + if(seconds > 1e+75) { } else if(seconds > 0) { cplexInstance.setParam(IloCplex::Param::TimeLimit, seconds); diff --git a/src/PrimalSolver.cpp b/src/PrimalSolver.cpp index a31c1aa4..1dc6e1ac 100644 --- a/src/PrimalSolver.cpp +++ b/src/PrimalSolver.cpp @@ -116,7 +116,7 @@ bool PrimalSolver::checkPrimalSolutionPoint(PrimalSolution primalSol) sourceDesc = "NLP fixed"; break; case E_PrimalSolutionSource::MIPSolutionPool: - sourceDesc = "MILP sol. pool"; + sourceDesc = "MIP sol. pool"; break; case E_PrimalSolutionSource::LPFixedIntegers: sourceDesc = "LP fixed"; @@ -127,6 +127,9 @@ bool PrimalSolver::checkPrimalSolutionPoint(PrimalSolution primalSol) case E_PrimalSolutionSource::InteriorPointSearch: sourceDesc = "Interior point search"; break; + case E_PrimalSolutionSource::ConvexBounding: + sourceDesc = "Convex MIP bounding"; + break; default: sourceDesc = "other"; break; diff --git a/src/Results.cpp b/src/Results.cpp index b9376457..8dee8ef2 100644 --- a/src/Results.cpp +++ b/src/Results.cpp @@ -121,14 +121,16 @@ void Results::addPrimalSolution(PrimalSolution solution) if(env->problem->objectiveFunction->properties.isMinimize) { std::sort(this->primalSolutions.begin(), this->primalSolutions.end(), - [](const PrimalSolution& firstSolution, const PrimalSolution& secondSolution) - { return (firstSolution.objValue < secondSolution.objValue); }); + [](const PrimalSolution& firstSolution, const PrimalSolution& secondSolution) { + return (firstSolution.objValue < secondSolution.objValue); + }); } else { std::sort(this->primalSolutions.begin(), this->primalSolutions.end(), - [](const PrimalSolution& firstSolution, const PrimalSolution& secondSolution) - { return (firstSolution.objValue > secondSolution.objValue); }); + [](const PrimalSolution& firstSolution, const PrimalSolution& secondSolution) { + return (firstSolution.objValue > secondSolution.objValue); + }); } env->solutionStatistics.numberOfFoundPrimalSolutions++; @@ -1135,7 +1137,7 @@ double Results::getCurrentDualBound() { return (this->currentDualBound); } double Results::getGlobalDualBound() { return (globalDualBound); } -void Results::setDualBound(double value) +void Results::setDualBound(double value, bool forceGlobal) { double primalBound = this->getPrimalBound(); @@ -1152,8 +1154,12 @@ void Results::setDualBound(double value) this->currentDualBound = value; - if(this->solutionIsGlobal) + if(this->solutionIsGlobal || forceGlobal) + { + env->output->outputInfo(fmt::format(" New dual bound: {} ", this->currentDualBound)); this->globalDualBound = value; + this->solutionIsGlobal = true; + } env->solutionStatistics.numberOfIterationsWithDualStagnation = 0; diff --git a/src/Results.h b/src/Results.h index 26a9a108..ebe49552 100644 --- a/src/Results.h +++ b/src/Results.h @@ -46,7 +46,7 @@ class DllExport Results void addDualSolution(DualSolution solution); double getCurrentDualBound(); double getGlobalDualBound(); - void setDualBound(double value); + void setDualBound(double value, bool forceGlobal = false); // For minimization problems, the lower bound is the dual while the upper bound is the primal objective value // for maximization problems, the lower bound is the primal while the upper bound is the dual objective value diff --git a/src/SolutionStrategy/SolutionStrategyMIQCQP.cpp b/src/SolutionStrategy/SolutionStrategyMIQCQP.cpp index f04dc44f..76f0805b 100644 --- a/src/SolutionStrategy/SolutionStrategyMIQCQP.cpp +++ b/src/SolutionStrategy/SolutionStrategyMIQCQP.cpp @@ -63,6 +63,8 @@ #include "../Settings.h" #include "../Timing.h" +#include "../Model/Problem.h" + namespace SHOT { diff --git a/src/SolutionStrategy/SolutionStrategyMultiTree.cpp b/src/SolutionStrategy/SolutionStrategyMultiTree.cpp index 1e9452ad..60cfb489 100644 --- a/src/SolutionStrategy/SolutionStrategyMultiTree.cpp +++ b/src/SolutionStrategy/SolutionStrategyMultiTree.cpp @@ -66,6 +66,8 @@ #include "../Tasks/TaskAddIntegerCuts.h" +#include "../Tasks/TaskPerformDualBounding.h" + #include "../Output.h" #include "../Model/Problem.h" #include "../Model/ObjectiveFunction.h" @@ -137,6 +139,12 @@ SolutionStrategyMultiTree::SolutionStrategyMultiTree(EnvironmentPtr envPtr) auto tSolveIteration = std::make_shared(env); env->tasks->addTask(tSolveIteration, "SolveIter"); + // if(env->reformulatedProblem->properties.convexity != E_ProblemConvexity::Convex) + // { + auto tPerformDualBounding = std::make_shared(env); + env->tasks->addTask(tPerformDualBounding, "DualBounding"); + //} + auto tSelectPrimSolPool = std::make_shared(env); env->tasks->addTask(tSelectPrimSolPool, "SelectPrimSolPool"); std::dynamic_pointer_cast(tFinalizeSolution)->addTask(tSelectPrimSolPool); diff --git a/src/Solver.cpp b/src/Solver.cpp index f4d8e192..9c5073a0 100644 --- a/src/Solver.cpp +++ b/src/Solver.cpp @@ -1939,6 +1939,9 @@ void Solver::setConvexityBasedSettings() env->settings->updateSetting("BoundTightening.FeasibilityBased.TimeLimit", "Model", 5.0); + // Need to save these to perform dual bound updates + // env->settings->updateSetting("HyperplaneCuts.SaveHyperplanePoints", "Dual", true); + #ifdef HAS_CPLEX if(static_cast(env->settings->getSetting("MIP.Solver", "Dual")) == ES_MIPSolver::Cplex) diff --git a/src/Structs.h b/src/Structs.h index 9b1a47b1..a5d1e502 100644 --- a/src/Structs.h +++ b/src/Structs.h @@ -179,23 +179,17 @@ struct Hyperplane int sourceConstraintIndex; // -1 if objective function VectorDouble generatedPoint; double objectiveFunctionValue; // Used for the objective cuts only - E_HyperplaneSource source; + E_HyperplaneSource source = E_HyperplaneSource::None; bool isObjectiveHyperplane = false; bool isSourceConvex = false; double pointHash; }; -struct GeneratedHyperplane +struct GeneratedHyperplane : Hyperplane { - NumericConstraintPtr sourceConstraint; - int sourceConstraintIndex; // -1 if objective function - VectorDouble generatedPoint; - E_HyperplaneSource source = E_HyperplaneSource::None; bool isLazy = false; bool isRemoved = false; - bool isSourceConvex = false; int iterationGenerated = -1; - double pointHash; }; struct IntegerCut @@ -231,6 +225,8 @@ struct SolutionStatistics int numberOfConstraintsRemovedInPresolve = 0; int numberOfVariableBoundsTightenedInPresolve = 0; + int numberOfHyperplanesWithConvexSource = 0; + int numberOfHyperplanesWithNonconvexSource = 0; int numberOfIntegerCuts = 0; int numberOfIterationsWithDualStagnation = 0; diff --git a/src/Tasks/TaskCreateDualProblem.cpp b/src/Tasks/TaskCreateDualProblem.cpp index fb1cafc1..52e5af97 100644 --- a/src/Tasks/TaskCreateDualProblem.cpp +++ b/src/Tasks/TaskCreateDualProblem.cpp @@ -16,6 +16,9 @@ #include "../Timing.h" #include "../MIPSolver/IMIPSolver.h" +#include "../Model/Problem.h" + +#include "TaskCreateMIPProblem.h" namespace SHOT { @@ -26,18 +29,12 @@ TaskCreateDualProblem::TaskCreateDualProblem(EnvironmentPtr envPtr) : TaskBase(e env->output->outputDebug(" Creating dual problem"); - createProblem(env->dualSolver->MIPSolver, env->reformulatedProblem); - - env->dualSolver->MIPSolver->finalizeProblem(); + taskCreateMIPProblem + = std::make_shared(env, env->dualSolver->MIPSolver, env->reformulatedProblem); + taskCreateMIPProblem->run(); env->dualSolver->MIPSolver->initializeSolverSettings(); - if(env->settings->getSetting("Debug.Enable", "Output")) - { - env->dualSolver->MIPSolver->writeProblemToFile( - env->settings->getSetting("Debug.Path", "Output") + "/dualiter0_problem.lp"); - } - env->output->outputDebug(" Dual problem created"); env->timing->stopTimer("DualStrategy"); } @@ -53,11 +50,7 @@ void TaskCreateDualProblem::run() env->output->outputDebug(" Recreating dual problem"); - createProblem(env->dualSolver->MIPSolver, env->reformulatedProblem); - - env->dualSolver->MIPSolver->finalizeProblem(); - - env->dualSolver->MIPSolver->initializeSolverSettings(); + taskCreateMIPProblem->run(); if(env->settings->getSetting("Debug.Enable", "Output")) { @@ -70,176 +63,6 @@ void TaskCreateDualProblem::run() } } -bool TaskCreateDualProblem::createProblem(MIPSolverPtr destination, ProblemPtr sourceProblem) -{ - // Now creating the variables - - bool variablesInitialized = true; - - for(auto& V : sourceProblem->allVariables) - { - variablesInitialized = variablesInitialized - && destination->addVariable( - V->name.c_str(), V->properties.type, V->lowerBound, V->upperBound, V->semiBound); - } - - if(!variablesInitialized) - return false; - - if(sourceProblem->auxiliaryObjectiveVariable) // The source problem already has a nonlinear objective variable - { - destination->setDualAuxiliaryObjectiveVariableIndex(sourceProblem->auxiliaryObjectiveVariable->index); - } - else if(sourceProblem->objectiveFunction->properties.classification > E_ObjectiveFunctionClassification::Quadratic) - { - double objVarBound = env->settings->getSetting("Variables.NonlinearObjectiveVariable.Bound", "Model"); - - Interval objectiveBound; - - try - { - objectiveBound = sourceProblem->objectiveFunction->getBounds(); - } - catch(mc::Interval::Exceptions&) - { - objectiveBound = Interval(-objVarBound, objVarBound); - } - - destination->setDualAuxiliaryObjectiveVariableIndex(sourceProblem->properties.numberOfVariables); - - destination->addVariable("shot_dual_objvar", E_VariableType::Real, objectiveBound.l(), objectiveBound.u(), 0.0); - - env->output->outputDebug(fmt::format( - " SHOT internal dual objective variable created with index {} and bounds [{},{}] created.", - sourceProblem->properties.numberOfVariables, objectiveBound.l(), objectiveBound.u())); - } - - // Now creating the objective function - - bool objectiveInitialized = true; - - objectiveInitialized = objectiveInitialized && destination->initializeObjective(); - - if(destination->hasDualAuxiliaryObjectiveVariable()) - { - objectiveInitialized = objectiveInitialized - && destination->addLinearTermToObjective(1.0, destination->getDualAuxiliaryObjectiveVariableIndex()); - - objectiveInitialized = objectiveInitialized - && destination->finalizeObjective(sourceProblem->objectiveFunction->properties.isMinimize); - } - else - { - // Linear terms - for(auto& T : std::dynamic_pointer_cast(sourceProblem->objectiveFunction)->linearTerms) - { - objectiveInitialized - = objectiveInitialized && destination->addLinearTermToObjective(T->coefficient, T->variable->index); - } - - // Quadratic terms - if(sourceProblem->objectiveFunction->properties.hasQuadraticTerms) - { - for(auto& T : - std::dynamic_pointer_cast(sourceProblem->objectiveFunction)->quadraticTerms) - { - objectiveInitialized = objectiveInitialized - && destination->addQuadraticTermToObjective( - T->coefficient, T->firstVariable->index, T->secondVariable->index); - } - } - - objectiveInitialized = objectiveInitialized - && destination->finalizeObjective( - sourceProblem->objectiveFunction->properties.isMinimize, sourceProblem->objectiveFunction->constant); - } - - if(!objectiveInitialized) - return false; - - // Now creating the constraints - - bool constraintsInitialized = true; - - for(auto& C : sourceProblem->linearConstraints) - { - constraintsInitialized = constraintsInitialized && destination->initializeConstraint(); - - if(C->properties.hasLinearTerms) - { - for(auto& T : C->linearTerms) - { - constraintsInitialized = constraintsInitialized - && destination->addLinearTermToConstraint(T->coefficient, T->variable->index); - } - } - - constraintsInitialized - = constraintsInitialized && destination->finalizeConstraint(C->name, C->valueLHS, C->valueRHS, C->constant); - } - - for(auto& C : sourceProblem->quadraticConstraints) - { - constraintsInitialized = constraintsInitialized && destination->initializeConstraint(); - - if(C->properties.hasLinearTerms) - { - for(auto& T : C->linearTerms) - { - constraintsInitialized = constraintsInitialized - && destination->addLinearTermToConstraint(T->coefficient, T->variable->index); - } - } - - if(C->properties.hasQuadraticTerms) - { - for(auto& T : C->quadraticTerms) - { - constraintsInitialized = constraintsInitialized - && destination->addQuadraticTermToConstraint( - T->coefficient, T->firstVariable->index, T->secondVariable->index); - } - } - - constraintsInitialized - = constraintsInitialized && destination->finalizeConstraint(C->name, C->valueLHS, C->valueRHS, C->constant); - } - - if(!constraintsInitialized) - return false; - - bool SOSInitialized = true; - - for(auto& S : sourceProblem->specialOrderedSets) - { - std::vector variableIndexes; - std::vector variableWeights; - - for(auto& V : S->variables) - variableIndexes.push_back(V->index); - - if(S->weights.size() > 0) - { - for(auto& W : S->weights) - variableWeights.push_back(W); - - SOSInitialized - = SOSInitialized && destination->addSpecialOrderedSet(S->type, variableIndexes, variableWeights); - } - else - { - SOSInitialized = SOSInitialized && destination->addSpecialOrderedSet(S->type, variableIndexes); - } - } - - if(!SOSInitialized) - return false; - - bool problemFinalized = destination->finalizeProblem(); - - return (problemFinalized); -} - std::string TaskCreateDualProblem::getType() { std::string type = typeid(this).name(); diff --git a/src/Tasks/TaskCreateDualProblem.h b/src/Tasks/TaskCreateDualProblem.h index c25f0407..b856eb29 100644 --- a/src/Tasks/TaskCreateDualProblem.h +++ b/src/Tasks/TaskCreateDualProblem.h @@ -11,12 +11,10 @@ #pragma once #include "TaskBase.h" -#include "../MIPSolver/IMIPSolver.h" -#include "../Model/Problem.h" +#include "TaskCreateMIPProblem.h" namespace SHOT { - class TaskCreateDualProblem : public TaskBase { public: @@ -27,6 +25,6 @@ class TaskCreateDualProblem : public TaskBase std::string getType() override; private: - bool createProblem(MIPSolverPtr destinationProblem, ProblemPtr sourceProblem); + std::shared_ptr taskCreateMIPProblem; }; } // namespace SHOT \ No newline at end of file diff --git a/src/Tasks/TaskCreateMIPProblem.cpp b/src/Tasks/TaskCreateMIPProblem.cpp new file mode 100644 index 00000000..f2b675c8 --- /dev/null +++ b/src/Tasks/TaskCreateMIPProblem.cpp @@ -0,0 +1,220 @@ +/** + The Supporting Hyperplane Optimization Toolkit (SHOT). + + @author Andreas Lundell, Åbo Akademi University + + @section LICENSE + This software is licensed under the Eclipse Public License 2.0. + Please see the README and LICENSE files for more information. +*/ + +#include "TaskCreateMIPProblem.h" + +#include "../DualSolver.h" +#include "../Output.h" +#include "../Settings.h" +#include "../Timing.h" + +#include "../MIPSolver/IMIPSolver.h" +#include "../Model/Problem.h" + +namespace SHOT +{ + +TaskCreateMIPProblem::TaskCreateMIPProblem(EnvironmentPtr envPtr, MIPSolverPtr MIPSolver, ProblemPtr sourceProblem) + : TaskBase(envPtr) +{ + this->MIPSolver = MIPSolver; + this->sourceProblem = sourceProblem; +} + +TaskCreateMIPProblem::~TaskCreateMIPProblem() = default; + +void TaskCreateMIPProblem::run() +{ + env->output->outputDebug(" Creating MIP problem"); + + createProblem(this->MIPSolver, this->sourceProblem); + + this->MIPSolver->finalizeProblem(); + this->MIPSolver->initializeSolverSettings(); + + env->output->outputDebug(" MIP problem created"); +} + +bool TaskCreateMIPProblem::createProblem(MIPSolverPtr destination, ProblemPtr sourceProblem) +{ + // Now creating the variables + + bool variablesInitialized = true; + + for(auto& V : sourceProblem->allVariables) + { + variablesInitialized = variablesInitialized + && destination->addVariable( + V->name.c_str(), V->properties.type, V->lowerBound, V->upperBound, V->semiBound); + } + + if(!variablesInitialized) + return false; + + if(sourceProblem->auxiliaryObjectiveVariable) // The source problem already has a nonlinear objective variable + { + destination->setDualAuxiliaryObjectiveVariableIndex(sourceProblem->auxiliaryObjectiveVariable->index); + } + else if(sourceProblem->objectiveFunction->properties.classification > E_ObjectiveFunctionClassification::Quadratic) + { + double objVarBound = env->settings->getSetting("Variables.NonlinearObjectiveVariable.Bound", "Model"); + + Interval objectiveBound; + + try + { + objectiveBound = sourceProblem->objectiveFunction->getBounds(); + } + catch(mc::Interval::Exceptions&) + { + objectiveBound = Interval(-objVarBound, objVarBound); + } + + destination->setDualAuxiliaryObjectiveVariableIndex(sourceProblem->properties.numberOfVariables); + + destination->addVariable("shot_dual_objvar", E_VariableType::Real, objectiveBound.l(), objectiveBound.u(), 0.0); + + env->output->outputDebug(fmt::format( + " SHOT internal dual objective variable created with index {} and bounds [{},{}] created.", + sourceProblem->properties.numberOfVariables, objectiveBound.l(), objectiveBound.u())); + } + + // Now creating the objective function + + bool objectiveInitialized = true; + + objectiveInitialized = objectiveInitialized && destination->initializeObjective(); + + if(destination->hasDualAuxiliaryObjectiveVariable()) + { + objectiveInitialized = objectiveInitialized + && destination->addLinearTermToObjective(1.0, destination->getDualAuxiliaryObjectiveVariableIndex()); + + objectiveInitialized = objectiveInitialized + && destination->finalizeObjective(sourceProblem->objectiveFunction->properties.isMinimize); + } + else + { + // Linear terms + for(auto& T : std::dynamic_pointer_cast(sourceProblem->objectiveFunction)->linearTerms) + { + objectiveInitialized + = objectiveInitialized && destination->addLinearTermToObjective(T->coefficient, T->variable->index); + } + + // Quadratic terms + if(sourceProblem->objectiveFunction->properties.hasQuadraticTerms) + { + for(auto& T : + std::dynamic_pointer_cast(sourceProblem->objectiveFunction)->quadraticTerms) + { + objectiveInitialized = objectiveInitialized + && destination->addQuadraticTermToObjective( + T->coefficient, T->firstVariable->index, T->secondVariable->index); + } + } + + objectiveInitialized = objectiveInitialized + && destination->finalizeObjective( + sourceProblem->objectiveFunction->properties.isMinimize, sourceProblem->objectiveFunction->constant); + } + + if(!objectiveInitialized) + return false; + + // Now creating the constraints + + bool constraintsInitialized = true; + + for(auto& C : sourceProblem->linearConstraints) + { + constraintsInitialized = constraintsInitialized && destination->initializeConstraint(); + + if(C->properties.hasLinearTerms) + { + for(auto& T : C->linearTerms) + { + constraintsInitialized = constraintsInitialized + && destination->addLinearTermToConstraint(T->coefficient, T->variable->index); + } + } + + constraintsInitialized + = constraintsInitialized && destination->finalizeConstraint(C->name, C->valueLHS, C->valueRHS, C->constant); + } + + for(auto& C : sourceProblem->quadraticConstraints) + { + constraintsInitialized = constraintsInitialized && destination->initializeConstraint(); + + if(C->properties.hasLinearTerms) + { + for(auto& T : C->linearTerms) + { + constraintsInitialized = constraintsInitialized + && destination->addLinearTermToConstraint(T->coefficient, T->variable->index); + } + } + + if(C->properties.hasQuadraticTerms) + { + for(auto& T : C->quadraticTerms) + { + constraintsInitialized = constraintsInitialized + && destination->addQuadraticTermToConstraint( + T->coefficient, T->firstVariable->index, T->secondVariable->index); + } + } + + constraintsInitialized + = constraintsInitialized && destination->finalizeConstraint(C->name, C->valueLHS, C->valueRHS, C->constant); + } + + if(!constraintsInitialized) + return false; + + bool SOSInitialized = true; + + for(auto& S : sourceProblem->specialOrderedSets) + { + std::vector variableIndexes; + std::vector variableWeights; + + for(auto& V : S->variables) + variableIndexes.push_back(V->index); + + if(S->weights.size() > 0) + { + for(auto& W : S->weights) + variableWeights.push_back(W); + + SOSInitialized + = SOSInitialized && destination->addSpecialOrderedSet(S->type, variableIndexes, variableWeights); + } + else + { + SOSInitialized = SOSInitialized && destination->addSpecialOrderedSet(S->type, variableIndexes); + } + } + + if(!SOSInitialized) + return false; + + bool problemFinalized = destination->finalizeProblem(); + + return (problemFinalized); +} + +std::string TaskCreateMIPProblem::getType() +{ + std::string type = typeid(this).name(); + return (type); +} +} // namespace SHOT \ No newline at end of file diff --git a/src/Tasks/TaskCreateMIPProblem.h b/src/Tasks/TaskCreateMIPProblem.h new file mode 100644 index 00000000..93df4af3 --- /dev/null +++ b/src/Tasks/TaskCreateMIPProblem.h @@ -0,0 +1,32 @@ +/** + The Supporting Hyperplane Optimization Toolkit (SHOT). + + @author Andreas Lundell, Åbo Akademi University + + @section LICENSE + This software is licensed under the Eclipse Public License 2.0. + Please see the README and LICENSE files for more information. +*/ + +#pragma once +#include "TaskBase.h" + +namespace SHOT +{ + +class TaskCreateMIPProblem : public TaskBase +{ +public: + TaskCreateMIPProblem(EnvironmentPtr envPtr, MIPSolverPtr MIPSolver, ProblemPtr sourceProblem); + ~TaskCreateMIPProblem() override; + + void run() override; + std::string getType() override; + +private: + bool createProblem(MIPSolverPtr destinationProblem, ProblemPtr sourceProblem); + + MIPSolverPtr MIPSolver; + ProblemPtr sourceProblem; +}; +} // namespace SHOT \ No newline at end of file diff --git a/src/Tasks/TaskPerformDualBounding.cpp b/src/Tasks/TaskPerformDualBounding.cpp new file mode 100644 index 00000000..3631c3b3 --- /dev/null +++ b/src/Tasks/TaskPerformDualBounding.cpp @@ -0,0 +1,190 @@ +/** + The Supporting Hyperplane Optimization Toolkit (SHOT). + + @author Andreas Lundell, Åbo Akademi University + + @section LICENSE + This software is licensed under the Eclipse Public License 2.0. + Please see the README and LICENSE files for more information. +*/ + +#include "TaskPerformDualBounding.h" + +#include "../DualSolver.h" +#include "../Output.h" +#include "../Settings.h" +#include "../Timing.h" + +#include "../MIPSolver/IMIPSolver.h" + +#ifdef HAS_CPLEX +#include "../MIPSolver/MIPSolverCplex.h" +#include "../MIPSolver/MIPSolverCplexSingleTree.h" +#include "../MIPSolver/MIPSolverCplexSingleTreeLegacy.h" +#endif + +#ifdef HAS_GUROBI +#include "../MIPSolver/MIPSolverGurobi.h" +#include "../MIPSolver/MIPSolverGurobiSingleTree.h" +#endif + +#ifdef HAS_CBC +#include "../MIPSolver/MIPSolverCbc.h" +#endif + +#include "../Results.h" +#include "../Tasks/TaskCreateMIPProblem.h" + +#include "../PrimalSolver.h" + +namespace SHOT +{ + +TaskPerformDualBounding::TaskPerformDualBounding(EnvironmentPtr envPtr) : TaskBase(envPtr) { } + +TaskPerformDualBounding::~TaskPerformDualBounding() = default; + +void TaskPerformDualBounding::run() +{ + + if(env->solutionStatistics.numberOfHyperplanesWithConvexSource == 0) + { + env->output->outputInfo( + " Dual bounding not performed since no hyperplanes with convex source have been added."); + return; + } + + if(env->solutionStatistics.numberOfHyperplanesWithNonconvexSource == 0) + { + env->output->outputInfo( + " Dual bounding not performed since no hyperplanes with nonconvex source have been added."); + return; + } + + if(lastNumberOfHyperplanesWithNonconvexSource == env->solutionStatistics.numberOfHyperplanesWithNonconvexSource + && lastNumberOfHyperplanesWithConvexSource == env->solutionStatistics.numberOfHyperplanesWithConvexSource) + { + env->output->outputInfo( + " Dual bounding not performed since no hyperplanes with both convex and nonconvex source have been added."); + return; + } + + lastNumberOfHyperplanesWithConvexSource = env->solutionStatistics.numberOfHyperplanesWithConvexSource; + lastNumberOfHyperplanesWithNonconvexSource = env->solutionStatistics.numberOfHyperplanesWithNonconvexSource; + + env->output->outputInfo(" Performing dual bounding."); + + MIPSolverPtr MIPSolver; + +#ifdef HAS_CPLEX + if(env->results->usedMIPSolver == ES_MIPSolver::Cplex) + { + MIPSolver = MIPSolverPtr(std::make_shared(env)); + } +#endif + +#ifdef HAS_GUROBI + if(env->results->usedMIPSolver == ES_MIPSolver::Gurobi) + { + MIPSolver = MIPSolverPtr(std::make_shared(env)); + } +#endif + +#ifdef HAS_CBC + if(env->results->usedMIPSolver == ES_MIPSolver::Cbc) + { + MIPSolver = MIPSolverPtr(std::make_shared(env)); + } +#endif + + assert(MIPSolver); + + if(!MIPSolver->initializeProblem()) + throw Exception(" Cannot initialize selected MIP solver."); + + env->output->outputInfo(" Creating dual bounding problem"); + + taskCreateMIPProblem = std::make_shared(env, MIPSolver, env->reformulatedProblem); + taskCreateMIPProblem->run(); + + int numberHyperplanesAdded = 0; + + for(auto HP : env->dualSolver->generatedHyperplanes) + { + if(HP.isSourceConvex) + { + if(MIPSolver->createHyperplane((Hyperplane)HP)) + numberHyperplanesAdded++; + } + } + + env->output->outputInfo(fmt::format( + " Number of hyperplanes added: {}/{}", numberHyperplanesAdded, env->dualSolver->generatedHyperplanes.size())); + + if(env->settings->getSetting("Debug.Enable", "Output")) + { + MIPSolver->writeProblemToFile( + env->settings->getSetting("Debug.Path", "Output") + "/dualbounding_problem.lp"); + } + + auto timeLim = env->settings->getSetting("TimeLimit", "Termination") - env->timing->getElapsedTime("Total"); + MIPSolver->setTimeLimit(timeLim); + + if(MIPSolver->getDiscreteVariableStatus() && env->results->hasPrimalSolution()) + { + auto primalSol = env->results->primalSolution; + env->reformulatedProblem->augmentAuxiliaryVariableValues(primalSol); + MIPSolver->addMIPStart(primalSol); + } + + MIPSolver->setSolutionLimit(2100000000); + + env->output->outputInfo(" Dual bounding problem created"); + auto solutionStatus = MIPSolver->solveProblem(); + + env->output->outputInfo( + fmt::format(" Dual bounding problem solved with return code: {}", (int)solutionStatus)); + + auto solutionPoints = MIPSolver->getAllVariableSolutions(); + double objectiveBound = MIPSolver->getDualObjectiveValue(); + env->output->outputInfo(fmt::format(" Dual bounding problem objective bound: {}", objectiveBound)); + + int iterationNumber = env->results->getCurrentIteration()->iterationNumber; + + if(solutionPoints.size() > 0) + { + env->output->outputInfo( + fmt::format(" Number of solutions in solution pool: {} ", solutionPoints.size())); + + double objectiveValue = MIPSolver->getObjectiveValue(); + DualSolution sol = { solutionPoints.at(0).point, E_DualSolutionSource::ConvexBounding, objectiveValue, + iterationNumber, false }; + env->dualSolver->addDualSolutionCandidate(sol); + + if(env->reformulatedProblem->antiEpigraphObjectiveVariable) + { + for(auto& SOL : solutionPoints) + SOL.point.at(env->reformulatedProblem->antiEpigraphObjectiveVariable->index) = objectiveValue; + } + + env->primalSolver->addPrimalSolutionCandidates(solutionPoints, E_PrimalSolutionSource::ConvexBounding); + + for(auto& SOL : solutionPoints) + env->primalSolver->addFixedNLPCandidate(SOL.point, E_PrimalNLPSource::FirstSolutionNewDualBound, + SOL.objectiveValue, SOL.iterFound, SOL.maxDeviation); + } + else + { + DualSolution sol = { {}, E_DualSolutionSource::ConvexBounding, objectiveBound, iterationNumber, false }; + env->dualSolver->addDualSolutionCandidate(sol); + } + + env->output->outputInfo(" Dual bounding finished."); +} + +std::string TaskPerformDualBounding::getType() +{ + std::string type = typeid(this).name(); + return (type); +} +} // namespace SHOT \ No newline at end of file diff --git a/src/Tasks/TaskPerformDualBounding.h b/src/Tasks/TaskPerformDualBounding.h new file mode 100644 index 00000000..4ddfa8c8 --- /dev/null +++ b/src/Tasks/TaskPerformDualBounding.h @@ -0,0 +1,36 @@ +/** + The Supporting Hyperplane Optimization Toolkit (SHOT). + + @author Andreas Lundell, Åbo Akademi University + + @section LICENSE + This software is licensed under the Eclipse Public License 2.0. + Please see the README and LICENSE files for more information. +*/ + +#pragma once +#include "TaskBase.h" + +#include "../MIPSolver/IMIPSolver.h" +#include "../Model/Problem.h" + +#include "TaskCreateMIPProblem.h" + +namespace SHOT +{ + +class TaskPerformDualBounding : public TaskBase +{ +public: + TaskPerformDualBounding(EnvironmentPtr envPtr); + ~TaskPerformDualBounding() override; + + void run() override; + std::string getType() override; + +private: + std::shared_ptr taskCreateMIPProblem; + int lastNumberOfHyperplanesWithConvexSource = 0; + int lastNumberOfHyperplanesWithNonconvexSource = 0; +}; +} // namespace SHOT \ No newline at end of file diff --git a/src/Tasks/TaskSolveIteration.cpp b/src/Tasks/TaskSolveIteration.cpp index b5f2e6dc..92786b66 100644 --- a/src/Tasks/TaskSolveIteration.cpp +++ b/src/Tasks/TaskSolveIteration.cpp @@ -234,6 +234,10 @@ void TaskSolveIteration::run() else { env->output->outputDebug(" Dual solver reports no solutions found."); + + DualSolution sol = { {}, E_DualSolutionSource::MIPSolverBound, + env->dualSolver->MIPSolver->getDualObjectiveValue(), currIter->iterationNumber, false }; + env->dualSolver->addDualSolutionCandidate(sol); } currIter->usedMIPSolutionLimit = env->dualSolver->MIPSolver->getSolutionLimit();