From dfd4943e4b9ad2aa25f73f5c73cd60d0ae252c96 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Mon, 16 May 2022 21:57:17 -0400 Subject: [PATCH 01/26] fix typo checkFeasiblity --- src/MibSBilevel.cpp | 16 ++++++++++------ src/MibSBilevel.hpp | 2 +- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index 7189b322..e0621808 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -85,14 +85,18 @@ MibSBilevel::createBilevel(CoinPackedVector* sol, int numElements(sol->getNumElements()); // number of nonzero elements int * fixedInd = model_->fixedInd_; - if(!upperSolutionOrd_) + if(!upperSolutionOrd_){ upperSolutionOrd_ = new double[uN]; - if(!lowerSolutionOrd_) + } + if(!lowerSolutionOrd_){ lowerSolutionOrd_ = new double[lN]; - if(!optUpperSolutionOrd_) + } + if(!optUpperSolutionOrd_){ optUpperSolutionOrd_ = new double[uN]; - if(!optLowerSolutionOrd_) + } + if(!optLowerSolutionOrd_){ optLowerSolutionOrd_ = new double[lN]; + } CoinZeroN(upperSolutionOrd_, uN); CoinZeroN(lowerSolutionOrd_, lN); @@ -239,7 +243,7 @@ MibSBilevel::createBilevel(CoinPackedVector* sol, (solveSecondLevelWhenXVarsInt && isUpperIntegral_) || (solveSecondLevelWhenLVarsInt && isLinkVarsIntegral_) || (solveSecondLevelWhenLVarsFixed && isLinkVarsFixed_ )))){ - storeSol = checkBilevelFeasiblity(mibs->isRoot_); + storeSol = checkBilevelFeasibility(mibs->isRoot_); } } @@ -257,7 +261,7 @@ MibSBilevel::createBilevel(CoinPackedVector* sol, //############################################################################# MibSSolType -MibSBilevel::checkBilevelFeasiblity(bool isRoot) +MibSBilevel::checkBilevelFeasibility(bool isRoot) { bool warmStartLL(model_->MibSPar_->entry (MibSParams::warmStartLL)); diff --git a/src/MibSBilevel.hpp b/src/MibSBilevel.hpp index c24bcbd8..f017858f 100644 --- a/src/MibSBilevel.hpp +++ b/src/MibSBilevel.hpp @@ -104,7 +104,7 @@ class MibSBilevel { MibSSolType createBilevel(CoinPackedVector *sol, MibSModel *mibs=0); - MibSSolType checkBilevelFeasiblity(bool isRoot); + MibSSolType checkBilevelFeasibility(bool isRoot); void gutsOfDestructor(); private: From d1589912eaf04d874b5146d1de661774cec2f94d Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Mon, 16 May 2022 22:01:13 -0400 Subject: [PATCH 02/26] add slTargetGap to param --- src/MibSModel.cpp | 10 ++++++++++ src/MibSParams.cpp | 4 ++++ src/MibSParams.hpp | 1 + 3 files changed, 15 insertions(+) diff --git a/src/MibSModel.cpp b/src/MibSModel.cpp index 0003029b..81bb3e32 100644 --- a/src/MibSModel.cpp +++ b/src/MibSModel.cpp @@ -3745,5 +3745,15 @@ MibSModel::instanceStructure() std::cout << "Linking solution pool will not be used." << std::endl; } } + + //YX: Setting "slTargetGap" parameter + + if (printProblemInfo == true){ + if(MibSPar_->entry(MibSParams::slTargetGap) > -1){ + std::cout << "Second (lower) level optimality gap is set to "; + std::cout << MibSPar_->entry(MibSParams::slTargetGap) <<"%."<< std::endl; + } + } + std::cout << std::endl; } diff --git a/src/MibSParams.cpp b/src/MibSParams.cpp index 998dbb27..0aeeade0 100644 --- a/src/MibSParams.cpp +++ b/src/MibSParams.cpp @@ -244,6 +244,8 @@ MibSParams::createKeywordList() { keys_.push_back(make_pair(std::string("MibS_boundCutTimeLim"), AlpsParameter(AlpsDoublePar, boundCutTimeLim))); + keys_.push_back(make_pair(std::string("MibS_slTargetGap"), + AlpsParameter(AlpsDoublePar, slTargetGap))); //YX: param for setting SL gap } //############################################################################# @@ -386,6 +388,8 @@ MibSParams::setDefaultEntries() { //------------------------------------------------------------- setEntry(boundCutTimeLim, 3600); + + setEntry(slTargetGap, -1); //YX: SL gap default set to -1 //------------------------------------------------------------- // String Parameters diff --git a/src/MibSParams.hpp b/src/MibSParams.hpp index bd0c42be..c73efdba 100644 --- a/src/MibSParams.hpp +++ b/src/MibSParams.hpp @@ -99,6 +99,7 @@ class MibSParams : public AlpsParameterSet { /** Double parameters. */ enum dblParams{ boundCutTimeLim, + slTargetGap, //YX: lower/second level optimality gap endOfDblParams }; From 8e9a1fe4ba6c19ac95b1f0367516794de40cc8da Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Wed, 18 May 2022 21:57:06 -0400 Subject: [PATCH 03/26] add gap to feasibility check and intersection cut --- src/MibSBilevel.cpp | 56 +++++++++++++++++++------ src/MibSCutGenerator.cpp | 91 ++++++++++++++++++++++++++++++++++++---- src/MibSCutGenerator.hpp | 6 +-- 3 files changed, 129 insertions(+), 24 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index e0621808..00947d0d 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -283,6 +283,10 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) (MibSParams::computeBestUBWhenLVarsFixed)); int useLinkingSolutionPool(model_->MibSPar_->entry (MibSParams::useLinkingSolutionPool)); + + // YX: target optimality gap for bounded rationality + double targetGap(model_->MibSPar_->entry(MibSParams::slTargetGap)); + double timeLimit(model_->AlpsPar()->entry(AlpsParams::timeLimit)); double remainingTime(0.0), startTimeVF(0.0), startTimeUB(0.0); MibSSolType storeSol(MibSNoSol); @@ -513,17 +517,34 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) //step 15 /** Current solution is bilevel feasible **/ - if((fabs(objVal - lowerObj) < etol) && (isIntegral_)){ - LPSolStatus_ = MibSLPSolStatusFeasible; - useBilevelBranching_ = false; - shouldPrune_ = true; - storeSol = MibSRelaxationSol; - } - else{ - memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); - if(isUpperIntegral_){ - storeSol = MibSHeurSol; - } + // YX: can be combined by set gap=0; check later + if(targetGap < etol){ + if((fabs(objVal - lowerObj) < etol) && (isIntegral_)){ + LPSolStatus_ = MibSLPSolStatusFeasible; + useBilevelBranching_ = false; + shouldPrune_ = true; + storeSol = MibSRelaxationSol; + } + else{ + memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); + if(isUpperIntegral_){ + storeSol = MibSHeurSol; + } + } + }else{ + // YX: [d2^y - phi(A^x)]/abs(phi(A^x)) + if(((lowerObj - objVal) <= fabs(objVal) * targetGap/100) && (isIntegral_)){ + LPSolStatus_ = MibSLPSolStatusFeasible; + useBilevelBranching_ = false; + shouldPrune_ = true; + storeSol = MibSRelaxationSol; + } + else{ + memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); + if(isUpperIntegral_){ + storeSol = MibSHeurSol; + } + } } if(!shouldPrune_){ @@ -694,6 +715,7 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, std::string feasCheckSolver = model_->MibSPar_->entry(MibSParams::feasCheckSolver); + double targetGap(model_->MibSPar_->entry(MibSParams::slTargetGap)); OsiSolverInterface * nSolver; @@ -713,6 +735,14 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, int rowNum(model_->getNumOrigCons() + 1); int colNum(model_->getNumOrigVars()); int * uColIndices(model_->getUpperColInd()); + double gap = (targetGap < model_->etol_) ? 0.0 : targetGap; // YX: added SL gap + double objUb(0.0); + // YX: d2^y <= phi(A^x) + \delta * abs(phi(A^x)) + if(objValLL > 0){ + objUb = objValLL + objValLL* gap/100; + }else{ + objUb = objValLL - objValLL* gap/100; + } if(newOsi){ double objSense(model_->getLowerObjSense()); @@ -800,7 +830,7 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, newRow[i] = newRow[i] * objSense; } - rowUb[rowNum-1] = objValLL; + rowUb[rowNum-1] = objUb; // YX: gap added to phi(A^2x) rowLb[rowNum-1] = -1 * (oSolver->getInfinity()); CoinPackedMatrix * newMat = new CoinPackedMatrix(false, 0, 0); @@ -839,7 +869,7 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, } else{ nSolver = UBSolver_; - nSolver->setRowUpper(rowNum-1, objValLL); + nSolver->setRowUpper(rowNum-1, objUb); // YX: gap added to phi(A^2x) for(i = 0; i < uCols; i++){ index1 = uColIndices[i]; if(fixedInd[index1] == 1){ diff --git a/src/MibSCutGenerator.cpp b/src/MibSCutGenerator.cpp index f6750c76..dcf0d2d7 100644 --- a/src/MibSCutGenerator.cpp +++ b/src/MibSCutGenerator.cpp @@ -572,7 +572,8 @@ MibSCutGenerator::intersectionCuts(BcpsConstraintPool &conPool, CoinZeroN(uselessIneqs, lRows); CoinZeroN(lowerLevelSol, lCols); isTimeLimReached = false; - if (!findLowerLevelSol(uselessIneqs, lowerLevelSol, sol, isTimeLimReached)){ + if (!findLowerLevelSol(uselessIneqs, lowerLevelSol, optLowerSolution, + sol, isTimeLimReached)){ // YX: add y^* for nonzero gap case delete [] uselessIneqs; delete [] lowerLevelSol; goto TERM_INTERSECTIONCUT; @@ -609,8 +610,8 @@ MibSCutGenerator::intersectionCuts(BcpsConstraintPool &conPool, CoinZeroN(uselessIneqs, lRows); CoinZeroN(lowerLevelSol, lCols); isTimeLimReached = false; - if (!findLowerLevelSolWatermelonIC(uselessIneqs, lowerLevelSol, lpSol, - isTimeLimReached)){ + if (!findLowerLevelSolWatermelonIC(uselessIneqs, lowerLevelSol, + optLowerSolution, lpSol, isTimeLimReached)){ // YX: add y^* for nonzero gap case delete [] uselessIneqs; delete [] lowerLevelSol; goto TERM_INTERSECTIONCUT; @@ -811,7 +812,7 @@ MibSCutGenerator::intersectionCuts(BcpsConstraintPool &conPool, //############################################################################# bool MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, - const double *sol, bool &isTimeLimReached) + double *optLowerSol, const double *sol, bool &isTimeLimReached) { std::string feasCheckSolver(localModel_->MibSPar_->entry @@ -820,6 +821,10 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, (MibSParams::maxThreadsLL)); int whichCutsLL(localModel_->MibSPar_->entry (MibSParams::whichCutsLL)); + double targetGap(localModel_->MibSPar_->entry(MibSParams::slTargetGap)); + double etol(localModel_->etol_); + double gap = (targetGap < etol) ? 0.0 : targetGap; // YX: added SL gap + double templObj(0.0); // YX: for nonzero gap, track d^2y^* bool foundSolution(false); @@ -838,7 +843,7 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, int lRows(localModel_->getLowerRowNum()); int numBinCols(lRows); int newNumCols(lCols + numBinCols); - int newNumRows(lRows + 1); + int newNumRows = (targetGap > etol)? lRows + 2 : lRows + 1; // YX: nonzero gap add a constraint double lObjSense(localModel_->getLowerObjSense()); double *lObjCoeff(localModel_->getLowerObjCoeffs()); int *lRowInd(localModel_->getLowerRowInd()); @@ -935,6 +940,15 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, newMatrix->appendRow(addedRow); addedRow.clear(); } + + // YX: for nonzero gap add an constraint: d^2y >= d^2y^*(1+gap) + if(targetGap > etol){ + for(i = 0; i < lCols; i++){ + addedRow.insert(i, -lObjCoeff[i] * lObjSense); + } + newMatrix->appendRow(addedRow); + addedRow.clear(); + } //filling row bounds CoinFillN(newRowLb, newNumRows, -1 * infinity); @@ -958,6 +972,18 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, assert(0); } } + + //YX: set cosntraint UB: d^2y >= d^2y^*(1+gap); convert row sense to L + if(targetGap > etol){ + for(i = 0; i < lCols; i++){ + templObj += lObjSense * lObjCoeff[i] * optLowerSol[i]; // YX: track d^2y^* + } + if(templObj > 0){ + newRowUb[newNumRows-1] = -templObj - (templObj * gap/100); + }else{ + newRowUb[newNumRows-1] = -templObj + (templObj * gap/100); + } + } //filling col bounds for(i = 0; i < lCols; i++){ @@ -1084,7 +1110,9 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, //the optimal solution of relaxation which satisfies integrality requirements //throw CoinError("The MIP which gives the best lower-level sol, cannot be infeasible!", // "findLowerLevelSol", "MibSCutGenerator"); - } + if(targetGap > etol){ + std::cout << "Type2IC aux MILP with optimality gap is infeasible."<getLowerObjCoeffs()); double objSense(localModel_->getLowerObjSense()); bool getA2Matrix(false), getG2Matrix(false); - + + double targetGap(localModel_->MibSPar_->entry(MibSParams::slTargetGap)); + double gap = (targetGap < etol) ? 0.0 : targetGap; + double templObj(0.0); // YX: track SL optimal obj val + if(localModel_->getA2Matrix() == NULL){ getA2Matrix = true; } @@ -1181,6 +1213,16 @@ MibSCutGenerator::getAlphaIC(double** extRay, double* uselessIneqs, for(i = 0; i < lCols; i++){ colIndex = lColInd[i]; rhs[lRows] += objSense * lObjCoeffs[i] * (lpSol[colIndex] - lowerSolution[i]); + templObj += objSense * lObjCoeffs[i] * lowerSolution[i]; // YX: track d^2y^* + } + + // YX: type I intersection IC; -d^2y^* - gap*|-d^2y^*| + if ((targetGap > etol) && (!uselessIneqs)){ + if(templObj > 0){ + rhs[lRows] += -templObj * gap/100; + }else{ + rhs[lRows] += templObj * gap/100; + } } for (i = 0; i < numNonBasic; i++){ @@ -1284,6 +1326,11 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo (MibSParams::maxThreadsLL)); int whichCutsLL(localModel_->MibSPar_->entry (MibSParams::whichCutsLL)); + double targetGap(localModel_->MibSPar_->entry(MibSParams::slTargetGap)); + double etol(localModel_->etol_); + double gap = (targetGap < etol) ? 0.0 : targetGap; // YX: added SL gap + double templObj(0.0); // YX: for nonzero gap, track d^2y^* + double timeLimit(localModel_->AlpsPar()->entry(AlpsParams::timeLimit)); double remainingTime(0.0); bool foundSolution = false; @@ -1298,7 +1345,7 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo int lCols(localModel_->getLowerDim()); int lRows(localModel_->getLowerRowNum()); int numContCols(lRows + 2 * lCols); - int newNumCols(lCols + numContCols); + int newNumRows = (targetGap > etol)? (2 * lRows + 2 * lCols + 2) : (2 * lRows + 2 * lCols + 1); // YX: nonzero gap add a constraint int newNumRows(2 * lRows + 2 * lCols + 1); double lObjSense(localModel_->getLowerObjSense()); double *lObjCoeff(localModel_->getLowerObjCoeffs()); @@ -1383,6 +1430,15 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo addedRow.clear(); } + // YX: for nonzero gap add an constraint: d^2 \delta y >= -d^2\yhat + d^2y^*(1+gap) + if(targetGap > etol){ + for(i = 0; i < lCols; i++){ + addedRow.insert(i, -lObjCoeff[i] * lObjSense); + } + newMatrix->appendRow(addedRow); + addedRow.clear(); + } + //filling row bounds CoinFillN(newRowLb, newNumRows, -1 * infinity); @@ -1462,6 +1518,22 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo } nSolver->setRowUpper(2 * i + 1, rhs - lCoeffsTimesLpSol[i]); } + + //YX: set cosntraint UB: d^2 \delta y >= -d^2\yhat + d^2y^*(1+gap); convert row sense to L + if(targetGap > etol){ + rhs = 0; + for(i = 0; i < lCols; i++){ + colIndex = lColInd[i]; + rhs += lObjSense * lObjCoeffs[i] * (lpSol[colIndex] - lowerSolution[i]); + templObj += lObjSense * lObjCoeff[i] * lowerSolution[i]; // YX: track d^2y^* + } + if(templObj > 0){ + rhs += -templObj * gap/100; + }else{ + rhs += templObj * gap/100; + } + nSolver->setRowUpper(newNumRows-1, rhs); + } //modifying col bounds for(i = 0; i < lCols; i++){ @@ -1548,6 +1620,9 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo //std::cout << "current time = " << timeLimit - localModel_->broker_->subTreeTimer().getTime() << std::endl; //throw CoinError("The MIP which is solved for watermelon IC, cannot be infeasible!", // "findLowerLevelSolWatermelonIC", "MibSCutGenerator"); + if(targetGap > etol){ + std::cout << "Watermelon aux MILP with optimality gap is infeasible."< Date: Wed, 18 May 2022 22:56:43 -0400 Subject: [PATCH 04/26] fix bugs in the last commit... --- src/MibSCutGenerator.cpp | 21 +++++++++++---------- src/MibSModel.cpp | 2 +- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/src/MibSCutGenerator.cpp b/src/MibSCutGenerator.cpp index dcf0d2d7..14f77bcf 100644 --- a/src/MibSCutGenerator.cpp +++ b/src/MibSCutGenerator.cpp @@ -1110,9 +1110,10 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, //the optimal solution of relaxation which satisfies integrality requirements //throw CoinError("The MIP which gives the best lower-level sol, cannot be infeasible!", // "findLowerLevelSol", "MibSCutGenerator"); - if(targetGap > etol){ - std::cout << "Type2IC aux MILP with optimality gap is infeasible."< etol){ + // std::cout << "Type2IC aux MILP with optimality gap is infeasible."<MibSPar_->entry (MibSParams::feasCheckSolver)); @@ -1345,8 +1346,8 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo int lCols(localModel_->getLowerDim()); int lRows(localModel_->getLowerRowNum()); int numContCols(lRows + 2 * lCols); + int newNumCols(lCols + numContCols); int newNumRows = (targetGap > etol)? (2 * lRows + 2 * lCols + 2) : (2 * lRows + 2 * lCols + 1); // YX: nonzero gap add a constraint - int newNumRows(2 * lRows + 2 * lCols + 1); double lObjSense(localModel_->getLowerObjSense()); double *lObjCoeff(localModel_->getLowerObjCoeffs()); int *lColInd(localModel_->getLowerColInd()); @@ -1524,8 +1525,8 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo rhs = 0; for(i = 0; i < lCols; i++){ colIndex = lColInd[i]; - rhs += lObjSense * lObjCoeffs[i] * (lpSol[colIndex] - lowerSolution[i]); - templObj += lObjSense * lObjCoeff[i] * lowerSolution[i]; // YX: track d^2y^* + rhs += lObjSense * lObjCoeff[i] * (lpSol[colIndex] - optLowerSol[i]); + templObj += lObjSense * lObjCoeff[i] * optLowerSol[i]; // YX: track d^2y^* } if(templObj > 0){ rhs += -templObj * gap/100; @@ -1620,9 +1621,9 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo //std::cout << "current time = " << timeLimit - localModel_->broker_->subTreeTimer().getTime() << std::endl; //throw CoinError("The MIP which is solved for watermelon IC, cannot be infeasible!", // "findLowerLevelSolWatermelonIC", "MibSCutGenerator"); - if(targetGap > etol){ - std::cout << "Watermelon aux MILP with optimality gap is infeasible."< etol){ + // std::cout << "Watermelon aux MILP with optimality gap is infeasible."<entry(MibSParams::slTargetGap) > -1){ From 07db86e84d008134ff8faa7fc479da410a88d641 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Fri, 10 Jun 2022 23:15:10 -0400 Subject: [PATCH 05/26] add pessimistic related vars and params --- src/MibSConstants.hpp | 2 ++ src/MibSHelp.hpp | 2 ++ src/MibSParams.cpp | 4 ++++ src/MibSParams.hpp | 2 ++ 4 files changed, 10 insertions(+) diff --git a/src/MibSConstants.hpp b/src/MibSConstants.hpp index fb0ac47d..a9a494d9 100644 --- a/src/MibSConstants.hpp +++ b/src/MibSConstants.hpp @@ -58,6 +58,8 @@ enum MibSLinkingPoolTag{ MibSLinkingPoolTagIsNotSet = -4, MibSLinkingPoolTagLowerIsInfeasible, MibSLinkingPoolTagLowerIsFeasible, + MibSLinkingPoolTagPesIsInfeasible, // YX: pessimistic case + MibSLinkingPoolTagPesIsFeasible, // YX: pessimistic case MibSLinkingPoolTagUBIsSolved }; diff --git a/src/MibSHelp.hpp b/src/MibSHelp.hpp index ede7d102..6999105a 100644 --- a/src/MibSHelp.hpp +++ b/src/MibSHelp.hpp @@ -21,8 +21,10 @@ struct LINKING_SOLUTION{ int tag; double lowerObjValue; + double pesRFValue; // YX: pessimistic case double UBObjValue; std::vector lowerSolution; + std::vector pesSolution; // YX: pessimistic case std::vector UBSolution; }; diff --git a/src/MibSParams.cpp b/src/MibSParams.cpp index 0aeeade0..57ce8bef 100644 --- a/src/MibSParams.cpp +++ b/src/MibSParams.cpp @@ -66,6 +66,8 @@ MibSParams::createKeywordList() { keys_.push_back(make_pair(std::string("MibS_useNewPureIntCut"), AlpsParameter(AlpsBoolPar, useNewPureIntCut))); + keys_.push_back(make_pair(std::string("MibS_findPesSol"), + AlpsParameter(AlpsBoolPar, findPesSol))); // YX: pessimistic case //-------------------------------------------------------- // BoolArrayPar //-------------------------------------------------------- @@ -281,6 +283,8 @@ MibSParams::setDefaultEntries() { setEntry(useNewPureIntCut, false); + setEntry(findPesSol, false); + //------------------------------------------------------------- // Int Parameters. //------------------------------------------------------------- diff --git a/src/MibSParams.hpp b/src/MibSParams.hpp index c73efdba..435c83d6 100644 --- a/src/MibSParams.hpp +++ b/src/MibSParams.hpp @@ -39,6 +39,8 @@ class MibSParams : public AlpsParameterSet { printProblemInfo, allowRemoveCut, useNewPureIntCut, + // YX: pessimistic case + findPesSol, endOfBoolParams }; From 632da9266968d3f922323d30f29b694a6adb49e2 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Sun, 12 Jun 2022 04:32:58 -0400 Subject: [PATCH 06/26] mv mat coeff setup func to Model from CutGenerator --- src/MibSCutGenerator.cpp | 65 ++++++++++++-------- src/MibSCutGenerator.hpp | 2 +- src/MibSModel.cpp | 124 ++++++++++++++++++++++++++++++++++++++- src/MibSModel.hpp | 28 ++++++++- 4 files changed, 190 insertions(+), 29 deletions(-) diff --git a/src/MibSCutGenerator.cpp b/src/MibSCutGenerator.cpp index 14f77bcf..546ef7a1 100644 --- a/src/MibSCutGenerator.cpp +++ b/src/MibSCutGenerator.cpp @@ -891,16 +891,20 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, int *integerVars = new int[newNumCols]; CoinZeroN(integerVars, newNumCols); - if(!localModel_->getA2Matrix()){ - getA2Matrix = true; - } + // if(!localModel_->getA2Matrix()){ + // getA2Matrix = true; + // } - if(!localModel_->getG2Matrix()){ - getG2Matrix = true; - } + // if(!localModel_->getG2Matrix()){ + // getG2Matrix = true; + // } - if((getA2Matrix) || (getG2Matrix)){ - getLowerMatrices(false, getA2Matrix, getG2Matrix); + // if((getA2Matrix) || (getG2Matrix)){ + // getLowerMatrices(false, getA2Matrix, getG2Matrix); + // } + + if(!localModel_->getA2Matrix()){ + localModel_->setCoeffMatrices(); } //extracting optimal first-level solution of the relaxation problem and @@ -1165,18 +1169,22 @@ MibSCutGenerator::getAlphaIC(double** extRay, double* uselessIneqs, double gap = (targetGap < etol) ? 0.0 : targetGap; double templObj(0.0); // YX: track SL optimal obj val - if(localModel_->getA2Matrix() == NULL){ - getA2Matrix = true; - } + // if(localModel_->getA2Matrix() == NULL){ + // getA2Matrix = true; + // } - if(localModel_->getG2Matrix() == NULL){ - getG2Matrix = true; - } + // if(localModel_->getG2Matrix() == NULL){ + // getG2Matrix = true; + // } - if((getA2Matrix) || (getG2Matrix)){ - getLowerMatrices(false, getA2Matrix, getG2Matrix); + // if((getA2Matrix) || (getG2Matrix)){ + // getLowerMatrices(false, getA2Matrix, getG2Matrix); + // } + + if(!localModel_->getA2Matrix()){ + localModel_->setCoeffMatrices(); } - + CoinPackedMatrix *matrixA2 = localModel_->getA2Matrix(); CoinPackedMatrix *matrixG2 = localModel_->getG2Matrix(); @@ -1368,16 +1376,20 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo origMatrix.reverseOrdering(); - if(!localModel_->getLowerConstCoefMatrix()){ - getA2G2Matrix = true; - } + // if(!localModel_->getLowerConstCoefMatrix()){ + // getA2G2Matrix = true; + // } - if(!localModel_->getG2Matrix()){ - getG2Matrix = true; - } + // if(!localModel_->getG2Matrix()){ + // getG2Matrix = true; + // } + + // if((getA2G2Matrix) || (getG2Matrix)){ + // getLowerMatrices(getA2G2Matrix, false, getG2Matrix); + // } - if((getA2G2Matrix) || (getG2Matrix)){ - getLowerMatrices(getA2G2Matrix, false, getG2Matrix); + if(!localModel_->getLowerConstCoefMatrix()){ + localModel_->setCoeffMatrices(); } CoinPackedMatrix *matrixA2G2 = localModel_->getLowerConstCoefMatrix(); @@ -6580,6 +6592,7 @@ MibSCutGenerator::getBindingConsBasis() } //############################################################################# +/** void MibSCutGenerator::getLowerMatrices(bool getLowerConstCoefMatrix, bool getA2Matrix, bool getG2Matrix) @@ -6677,4 +6690,4 @@ MibSCutGenerator::getLowerMatrices(bool getLowerConstCoefMatrix, localModel_->setG2Matrix(matrixG2); } -} +}**/ diff --git a/src/MibSCutGenerator.hpp b/src/MibSCutGenerator.hpp index ff9f555e..39c0c4f3 100644 --- a/src/MibSCutGenerator.hpp +++ b/src/MibSCutGenerator.hpp @@ -186,7 +186,7 @@ class MibSCutGenerator : public BlisConGenerator { double * findDeepestLandPCut1();//stable (but maybe wrong) /** Store the matrices A2, G2 and lower-level coeffs (all constraints are 'L') **/ - void getLowerMatrices(bool getLowerConstCoefMatrix, bool getA2Matrix, bool getG2Matrix); + // void getLowerMatrices(bool getLowerConstCoefMatrix, bool getA2Matrix, bool getG2Matrix); private: diff --git a/src/MibSModel.cpp b/src/MibSModel.cpp index e9876173..e58d53ca 100644 --- a/src/MibSModel.cpp +++ b/src/MibSModel.cpp @@ -99,6 +99,8 @@ MibSModel::~MibSModel() if(rowName_) delete [] rowName_; if(MibSPar_) delete MibSPar_; if(lowerConstCoefMatrix_) delete lowerConstCoefMatrix_; + if(A1Matrix_) delete A1Matrix_; + if(G1Matrix_) delete G1Matrix_; if(A2Matrix_) delete A2Matrix_; if(G2Matrix_) delete G2Matrix_; if(bS_) delete bS_; @@ -134,6 +136,7 @@ MibSModel::initialize() counterVF_ = 0; counterUB_ = 0; timerVF_ = 0.0; + timerPES_ = 0.0; // YX: pessmistic case timerUB_ = 0.0; countIteration_ = 0; isInterdict_ = false; @@ -170,6 +173,8 @@ MibSModel::initialize() interdictCost_ = NULL; origConstCoefMatrix_ = NULL; lowerConstCoefMatrix_ = NULL; + A1Matrix_ = NULL; + G1Matrix_ = NULL; A2Matrix_ = NULL; G2Matrix_ = NULL; boundProbRoot_ = NULL; @@ -1048,6 +1053,11 @@ MibSModel::loadProblemData(const CoinPackedMatrix& matrix, //determine the list of first-stage variables participate in second-stage constraints setRequiredFixedList(); instanceStructure(); + + if(MibSPar_->entry(MibSParams::findPesSol)) + { + setCoeffMatrices(); + } } //############################################################################# @@ -3081,6 +3091,113 @@ MibSModel::decodeToSelf(AlpsEncoded& encoded) status = decodeMibS(encoded); } +//############################################################################# +void +MibSModel::setCoeffMatrices() +{ + /** Set coeff matrices using original input data (with 'L' sense) **/ + if(lowerConstCoefMatrix_){ + return; + } + + int i(0), j(0); + int index(0), numElements(0), pos(0); + const int *indices(0); + const double *elements(0); + + CoinPackedMatrix origMatrix(*origConstCoefMatrix_); + CoinShallowPackedVector origRow; + CoinPackedVector row, rowA, rowG; + CoinPackedMatrix *matrixA1(0); + CoinPackedMatrix *matrixG1(0); + CoinPackedMatrix *matrixA2G2(0); + CoinPackedMatrix *matrixA2(0); + CoinPackedMatrix *matrixG2(0); + + origMatrix.reverseOrdering(); + if(upperRowNum_ > 0){ + matrixA1 = new CoinPackedMatrix(false, 0, 0); + matrixA1->setDimensions(0, upperDim_); + matrixG1 = new CoinPackedMatrix(false, 0, 0); + matrixG1->setDimensions(0, lowerDim_); + } + matrixA2G2 = new CoinPackedMatrix(false, 0, 0); + matrixA2G2->setDimensions(0, numVars_); + assert(numVars_== upperDim_ + lowerDim_); // YX: debug only + matrixA2 = new CoinPackedMatrix(false, 0, 0); + matrixA2->setDimensions(0, upperDim_); + matrixG2 = new CoinPackedMatrix(false, 0, 0); + matrixG2->setDimensions(0, lowerDim_); + + // YX: set up A1 and G1 + // YX: need to check if it works for interdiction case with AUX rows + for(i = 0; i < upperRowNum_; i++){ + index = upperRowInd_[i]; + origRow = origMatrix.getVector(index); + if(origRowSense_[index] == 'G'){ + row = -1 * origRow; + }else{ + row = origRow; + } + numElements = row.getNumElements(); + indices = row.getIndices(); + elements = row.getElements(); + for(j = 0; j < numElements; j++){ + pos = binarySearch(0, upperDim_ -1, indices[j], upperColInd_); + if(pos >= 0){ + rowA.insert(pos, elements[j]); + }else{ + pos = binarySearch(0, lowerDim_ -1, indices[j], lowerColInd_); + rowG.insert(pos, elements[j]); + } + } + matrixA1->appendRow(rowA); + matrixG1->appendRow(rowG); + + rowA.clear(); + rowG.clear(); + } + + // YX: set up A2, G2, and A2G2 + for(i = 0; i < lowerRowNum_; i++){ + index = lowerRowInd_[i]; + origRow = origMatrix.getVector(index); + if(origRowSense_[index] == 'G'){ + row = -1 * origRow; + }else{ + row = origRow; + } + matrixA2G2->appendRow(row); + numElements = row.getNumElements(); + indices = row.getIndices(); + elements = row.getElements(); + for(j = 0; j < numElements; j++){ + pos = binarySearch(0, upperDim_ -1, indices[j], upperColInd_); + if(pos >= 0){ + rowA.insert(pos, elements[j]); + }else{ + pos = binarySearch(0, lowerDim_ -1, indices[j], lowerColInd_); + rowG.insert(pos, elements[j]); + } + } + matrixA2->appendRow(rowA); + matrixG2->appendRow(rowG); + + rowA.clear(); + rowG.clear(); + } + + if(matrixA1){ + A1Matrix_ = matrixA1; + } + if(matrixG1){ + G1Matrix_ = matrixG1; + } + lowerConstCoefMatrix_ = matrixA2G2; + A2Matrix_ = matrixA2; + G2Matrix_ = matrixG2; +} + //############################################################################# void MibSModel::setRequiredFixedList() @@ -3747,7 +3864,6 @@ MibSModel::instanceStructure() } //YX: Printing "slTargetGap" parameter; setting it as a model property? - if (printProblemInfo == true){ if(MibSPar_->entry(MibSParams::slTargetGap) > -1){ std::cout << "Second (lower) level optimality gap is set to "; @@ -3755,5 +3871,11 @@ MibSModel::instanceStructure() } } + //YX: Printing "findPesSol" parameter; + if (printProblemInfo == true){ + if(MibSPar_->entry(MibSParams::findPesSol) == PARAM_ON){ + std::cout << "Searching for a pessimistic solution." << std::endl; + } + } std::cout << std::endl; } diff --git a/src/MibSModel.hpp b/src/MibSModel.hpp index 4481ff99..7fbcc89b 100644 --- a/src/MibSModel.hpp +++ b/src/MibSModel.hpp @@ -105,6 +105,9 @@ class MibSModel : public BlisModel { /** Time for solving (VF) **/ double timerVF_; + /** Time for solving (PES) **/ + double timerPES_; // YX: pessimistic case; + /** Time for solving (UB) **/ double timerUB_; @@ -156,7 +159,7 @@ class MibSModel : public BlisModel { /** the right (positive) slope of the lower-level value function **/ double rightSlope_; - int countIteration_; + int countIteration_; // YX: HOW IS IT USED??? /** Indices of UL variables **/ int * upperColInd_; @@ -218,11 +221,18 @@ class MibSModel : public BlisModel { /** Original matrix of constraints coefficients **/ CoinPackedMatrix *origConstCoefMatrix_; + // YX: used in pessimistic feasibility check //we only fill next three matrices if they are required. //For now, we only need them for generating IC and watermelon IC. /** matrix of lower-level coefficients(all constraints are 'L') **/ CoinPackedMatrix *lowerConstCoefMatrix_; + /** matrix of upper-level vars in upper-level problem (with row sense 'L') **/ + CoinPackedMatrix *A1Matrix_; + + /** matrix of lower-level vars in upper-level problem (with row sense 'L') **/ + CoinPackedMatrix *G1Matrix_; + /** matrix of upper-level vars in lower-level problem(all constraints are 'L') **/ CoinPackedMatrix *A2Matrix_; @@ -403,6 +413,12 @@ class MibSModel : public BlisModel { /** Set pointer to the matrix of LL coefficients(all constraints are 'L') **/ void setLowerConstCoefMatrix(CoinPackedMatrix *ptr) {lowerConstCoefMatrix_ = ptr;} + /** YX: Set pointer to the matrix of UL vars in UL problem (with row sense 'L') **/ + void setA1Matrix(CoinPackedMatrix *ptr) {A1Matrix_ = ptr;} + + /** YX: Set pointer to the matrix of LL vars in UL problem (with row sense 'L') **/ + void setG1Matrix(CoinPackedMatrix *ptr) {G1Matrix_ = ptr;} + /** Set pointer to the matrix of UL vars in LL problem(all constraints are 'L') **/ void setA2Matrix(CoinPackedMatrix *ptr) {A2Matrix_ = ptr;} @@ -497,6 +513,12 @@ class MibSModel : public BlisModel { /** Get pointer to the matrix of LL coefficients(all constraints are 'L') **/ CoinPackedMatrix * getLowerConstCoefMatrix() const {return lowerConstCoefMatrix_;} + /** YX: Get pointer to the matrix of UL vars in UL problem (with row sense 'L') **/ + CoinPackedMatrix * getA1Matrix() const {return A1Matrix_;} + + /** YX: Get pointer to the matrix of LL vars in UL problem (with row sense 'L') **/ + CoinPackedMatrix * getG1Matrix() const {return G1Matrix_;} + /** Get pointer to the matrix of UL vars in LL problem(all constraints are 'L') **/ CoinPackedMatrix * getA2Matrix() const {return A2Matrix_;} @@ -605,6 +627,10 @@ class MibSModel : public BlisModel { /** The method that decodes the model from an encoded object. */ virtual void decodeToSelf(AlpsEncoded&); + /** YX: Separate and set A1, G1, A2, G2 matrices (with row sense 'L') **/ + // YX: this was moved from MibSCutGenerator + void setCoeffMatrices(); + /** Determine the list of first-stage variables participate in second-stage constraints */ void setRequiredFixedList(); From 6e36c9bc014880b21feccf9657b4b9f5e771777f Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Mon, 13 Jun 2022 01:26:19 -0400 Subject: [PATCH 07/26] add pessimistic MILP and related procedure --- src/MibSBilevel.cpp | 517 ++++++++++++++++++++++++++++++++++++++++---- src/MibSBilevel.hpp | 5 + 2 files changed, 477 insertions(+), 45 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index 00947d0d..4963b35b 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -226,6 +226,7 @@ MibSBilevel::createBilevel(CoinPackedVector* sol, //steps 5-6 if((isLinkVarsFixed_) && ((tagInSeenLinkingPool_ == MibSLinkingPoolTagLowerIsInfeasible) || + (tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsInfeasible) || // YX: pessimistic case (tagInSeenLinkingPool_ == MibSLinkingPoolTagUBIsSolved))){ useBilevelBranching_ = false; shouldPrune_ = true; @@ -234,6 +235,7 @@ MibSBilevel::createBilevel(CoinPackedVector* sol, //step 7 if(!shouldPrune_){ if((tagInSeenLinkingPool_ == MibSLinkingPoolTagLowerIsFeasible || + tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsFeasible || // YX: pessimistic case tagInSeenLinkingPool_ == MibSLinkingPoolTagUBIsSolved) || (!isContainedInLinkingPool_ && ((branchPar == MibSBranchingStrategyLinking && isIntegral_ && isLinkVarsFixed_) || @@ -265,6 +267,8 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) { bool warmStartLL(model_->MibSPar_->entry (MibSParams::warmStartLL)); + bool findPesSol(model_->MibSPar_->entry + (MibSParams::findPesSol)); // YX: pessimistic case int maxThreadsLL(model_->MibSPar_->entry (MibSParams::maxThreadsLL)); int whichCutsLL(model_->MibSPar_->entry @@ -286,22 +290,26 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) // YX: target optimality gap for bounded rationality double targetGap(model_->MibSPar_->entry(MibSParams::slTargetGap)); - double timeLimit(model_->AlpsPar()->entry(AlpsParams::timeLimit)); - double remainingTime(0.0), startTimeVF(0.0), startTimeUB(0.0); + double remainingTime(0.0), startTimeVF(0.0), startTimePES(0.0), startTimeUB(0.0); // YX: pessimistic case MibSSolType storeSol(MibSNoSol); int lN(model_->lowerDim_); // lower-level dimension int uN(model_->upperDim_); // upper-level dimension int i(0), index(0), length(0), pos(0); int sizeFixedInd(model_->sizeFixedInd_); double etol(model_->etol_), objVal(0.0), lowerObj(0.0); + double pesRF(0.0), lpPesVal(0.0), lowerPesVal(0.0); // YX: pessimistic case int * fixedInd = model_->fixedInd_; int * lowerColInd = model_->getLowerColInd(); int * upperColInd = model_->getUpperColInd(); + const double * valuesPes; double *lowerSol = new double[lN]; + double *pesSol = new double[lN + uN]; // YX: pessimistic case CoinFillN(lowerSol, lN, 0.0); + CoinFillN(pesSol, lN + uN, 0.0); std::vector shouldStoreValuesLowerSol(lN); + std::vector shouldStoreValuesPesSol(lN + uN); // YX: pessimistic case std::vector shouldStoreValuesUBSol(lN + uN); const double * sol = model_->solver()->getColSolution(); @@ -338,6 +346,9 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) OsiSolverInterface *lSolver = lSolver_; + lSolver->writeLp("lSolver"); + model_->getSolver()->writeLp("oSolver"); + remainingTime = timeLimit - model_->broker_->subTreeTimer().getTime(); if(remainingTime <= etol){ shouldPrune_ = true; @@ -486,6 +497,123 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) node->setLowerUB(objVal); node->setIsBoundSet(true); } + + // YX: pessimistic case; need to confirm repeated model setup? + if(findPesSol){ + if(pSolver_){ + pSolver_ = setUpPesModel(model_->getSolver(), objVal, false); + } + else{ + pSolver_ = setUpPesModel(model_->getSolver(), objVal, true); + } + + OsiSolverInterface *pSolver = pSolver_; + + pSolver->writeLp("pSolverLoaded"); + + // YX: add timer here + remainingTime = timeLimit - model_->broker_->subTreeTimer().getTime(); + + if(remainingTime <= etol){ + shouldPrune_ = true; + goto TERM_CHECKBILEVELFEAS; + } + + remainingTime = CoinMax(remainingTime, 0.00); + + if (feasCheckSolver == "Cbc"){ + dynamic_cast + (pSolver)->getModelPtr()->messageHandler()->setLogLevel(0); + }else if (feasCheckSolver == "SYMPHONY"){ +#if COIN_HAS_SYMPHONY + //dynamic_cast + // (pSolver)->setSymParam("prep_level", -1); + sym_environment *env = dynamic_cast + (pSolver)->getSymphonyEnvironment(); + //Always uncomment for debugging!! + sym_set_dbl_param(env, "time_limit", remainingTime); + sym_set_int_param(env, "do_primal_heuristic", FALSE); + sym_set_int_param(env, "verbosity", -2); + sym_set_int_param(env, "prep_level", -1); + sym_set_int_param(env, "max_active_nodes", maxThreadsLL); + sym_set_int_param(env, "tighten_root_bounds", FALSE); + sym_set_int_param(env, "max_sp_size", 100); + sym_set_int_param(env, "do_reduced_cost_fixing", FALSE); + if (whichCutsLL == 0){ + sym_set_int_param(env, "generate_cgl_cuts", FALSE); + }else{ + sym_set_int_param(env, "generate_cgl_gomory_cuts", GENERATE_DEFAULT); + } + if (whichCutsLL == 1){ + sym_set_int_param(env, "generate_cgl_knapsack_cuts", + DO_NOT_GENERATE); + sym_set_int_param(env, "generate_cgl_probing_cuts", + DO_NOT_GENERATE); + sym_set_int_param(env, "generate_cgl_clique_cuts", + DO_NOT_GENERATE); + sym_set_int_param(env, "generate_cgl_twomir_cuts", + DO_NOT_GENERATE); + sym_set_int_param(env, "generate_cgl_flowcover_cuts", + DO_NOT_GENERATE); + } +#endif + }else if (feasCheckSolver == "CPLEX"){ +#ifdef COIN_HAS_CPLEX + pSolver->setHintParam(OsiDoReducePrint); + pSolver->messageHandler()->setLogLevel(0); + CPXENVptr cpxEnv = + dynamic_cast(pSolver)->getEnvironmentPtr(); + assert(cpxEnv); + CPXsetintparam(cpxEnv, CPX_PARAM_SCRIND, CPX_OFF); + CPXsetintparam(cpxEnv, CPX_PARAM_THREADS, maxThreadsLL); +#endif + } + + startTimePES = model_->broker_->subTreeTimer().getTime(); + pSolver->branchAndBound(); + model_->timerPES_ += model_->broker_->subTreeTimer().getTime() - startTimePES; + // YX: use SL-MILP counter? + + if((feasCheckSolver == "SYMPHONY") && (sym_is_time_limit_reached + (dynamic_cast + (pSolver)->getSymphonyEnvironment()))){ + shouldPrune_ = true; + goto TERM_CHECKBILEVELFEAS; + } + + // isPESSolved_ = true; + if(!pSolver->isProvenOptimal()){ + isProvenOptimal_ = false; + if(useLinkingSolutionPool && isLinkVarsIntegral_){ + addSolutionToSeenLinkingSolutionPool + (MibSLinkingPoolTagPesIsInfeasible, shouldStoreValuesPesSol, 0.0); + } + if(isLinkVarsFixed_){ + useBilevelBranching_ = false; + shouldPrune_ = true; + } + }else{ + pesRF = pSolver->getObjValue(); + valuesPes = pSolver->getColSolution(); + + for (i = 0; i < uN + lN; i++){ + if ((pSolver->isInteger(i)) && + (((valuesPes[i] - floor(valuesPes[i])) < etol) || + ((ceil(valuesPes[i]) - valuesPes[i]) < etol))){ + pesSol[i] = (double) floor(valuesPes[i] + 0.5); + }else{ + pesSol[i] = (double) valuesPes[i]; + } + } + + if(useLinkingSolutionPool && isLinkVarsIntegral_){ + std::copy(pesSol, pesSol + lN + uN, shouldStoreValuesPesSol.begin()); + addSolutionToSeenLinkingSolutionPool + (MibSLinkingPoolTagPesIsFeasible, shouldStoreValuesPesSol, pesRF); + shouldStoreValuesPesSol.clear(); + } + } + } } /*if (!warmStartLL){ delete solver_; @@ -493,19 +621,23 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) } //step 13 - if(((!useLinkingSolutionPool || !isLinkVarsIntegral_) && (isProvenOptimal_)) || - ((tagInSeenLinkingPool_ == MibSLinkingPoolTagLowerIsFeasible) || - (tagInSeenLinkingPool_ == MibSLinkingPoolTagUBIsSolved))){ + if(((!useLinkingSolutionPool || !isLinkVarsIntegral_) && (isProvenOptimal_)) || + ((tagInSeenLinkingPool_ == MibSLinkingPoolTagLowerIsFeasible) || + (tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsFeasible) || + (tagInSeenLinkingPool_ == MibSLinkingPoolTagUBIsSolved))){ if(useLinkingSolutionPool && isLinkVarsIntegral_){ - //get optimal value of (VF) from solution pool - //model_->it = seenLinkingSolutions.find(linkSol); - //objVal = model_->it->second.lowerObjVal1; + //get optimal value of (VF) from solution pool objVal = model_->seenLinkingSolutions[linkSol].lowerObjValue; - //objVal = seenLinkingSolutions.find(linkSol). objVal_ = objVal; std::copy(model_->seenLinkingSolutions[linkSol].lowerSolution.begin(), model_->seenLinkingSolutions[linkSol].lowerSolution.end(), lowerSol); + + if(findPesSol){ // YX: pessimistic case + pesRF = model_->seenLinkingSolutions[linkSol].pesRFValue; + std::copy(model_->seenLinkingSolutions[linkSol].pesSolution.begin(), + model_->seenLinkingSolutions[linkSol].pesSolution.end(), pesSol); + } } lowerObj = getLowerObj(sol, model_->getLowerObjSense()); @@ -517,35 +649,44 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) //step 15 /** Current solution is bilevel feasible **/ - // YX: can be combined by set gap=0; check later - if(targetGap < etol){ - if((fabs(objVal - lowerObj) < etol) && (isIntegral_)){ - LPSolStatus_ = MibSLPSolStatusFeasible; - useBilevelBranching_ = false; - shouldPrune_ = true; - storeSol = MibSRelaxationSol; - } - else{ - memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); - if(isUpperIntegral_){ - storeSol = MibSHeurSol; - } - } - }else{ - // YX: [d2^y - phi(A^x)]/abs(phi(A^x)) - if(((lowerObj - objVal) <= fabs(objVal) * targetGap/100) && (isIntegral_)){ - LPSolStatus_ = MibSLPSolStatusFeasible; - useBilevelBranching_ = false; - shouldPrune_ = true; - storeSol = MibSRelaxationSol; - } - else{ - memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); - if(isUpperIntegral_){ - storeSol = MibSHeurSol; - } - } - } + // YX: [d2^y - phi(A^x)]/abs(phi(A^x)), combined check when gap is set to 0; + if(!findPesSol){ + if(((lowerObj - objVal) <= fabs(objVal) * targetGap/100) && (isIntegral_)){ + LPSolStatus_ = MibSLPSolStatusFeasible; + useBilevelBranching_ = false; + shouldPrune_ = true; + storeSol = MibSRelaxationSol; + }else{ + memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); + if(isUpperIntegral_){ + storeSol = MibSHeurSol; + } + } + }else{ + // YX: find pessimistic risk function value using (LR) and (SL-MILP) solutions + lpPesVal = getRiskFuncVal(model_->getSolver(), lowerSolutionOrd_, findPesSol); + // YX: whether (LR) solution is a pessimistic solution + if(abs(lpPesVal - pesRF) < etol){ + useBilevelBranching_ = false; + shouldPrune_ = true; + storeSol = MibSRelaxationSol; + }else{ + lowerPesVal = getRiskFuncVal(model_->getSolver(), lowerSol, findPesSol); + // YX: whether (SL-MILP) solution is a pessimistic solution + if(abs(lowerPesVal - pesRF) < etol){ + memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); + if(isUpperIntegral_){ + storeSol = MibSHeurSol; + } + }else{ + // copy pessimistic solution to optLower; + memcpy(optLowerSolutionOrd_, pesSol, sizeof(double) * lN); + if(isUpperIntegral_){ + storeSol = MibSHeurSol; + } + } + } + } if(!shouldPrune_){ //step 18 @@ -555,13 +696,28 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) ((computeBestUBWhenXVarsInt == PARAM_ON) && (isUpperIntegral_)) || ((computeBestUBWhenLVarsInt == PARAM_ON) && (isLinkVarsIntegral_)) || ((computeBestUBWhenLVarsFixed == PARAM_ON) && (isLinkVarsFixed_)))){ - if(UBSolver_){ - UBSolver_ = setUpUBModel(model_->getSolver(), objVal, false); - } - else{ - UBSolver_ =setUpUBModel(model_->getSolver(), objVal, true); - } - + if(!findPesSol){ + if(UBSolver_){ + UBSolver_ = setUpUBModel(model_->getSolver(), objVal, false); + } + else{ + UBSolver_ =setUpUBModel(model_->getSolver(), objVal, true); + } + }else{ // YX: set up later + if(UBSolver_){ + UBSolver_ = setUpUBModel(model_->getSolver(), objVal, false); + } + else{ + UBSolver_ =setUpUBModel(model_->getSolver(), objVal, true); + } + // if((sizeFixedInd - uN) < etol){ + // all variables are linking variables; solution fixed + // use and store existing bound, update tag, then jump to term; + // }else{ + // exists non-fixed upper vars, set up solver in a special way + // add another function to set up pes UB + // } + } OsiSolverInterface * UBSolver = UBSolver_; remainingTime = timeLimit - model_->broker_->subTreeTimer().getTime(); @@ -688,6 +844,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) TERM_CHECKBILEVELFEAS: delete [] lowerSol; + delete [] pesSol; return storeSol; } @@ -702,6 +859,7 @@ MibSBilevel::gutsOfDestructor() if(upperSolutionOrd_) delete [] upperSolutionOrd_; if(lowerSolutionOrd_) delete [] lowerSolutionOrd_; if(lSolver_) delete lSolver_; + if(pSolver_) delete pSolver_; // YX: pessimistic case if(UBSolver_) delete UBSolver_; //delete heuristic_; } @@ -884,6 +1042,230 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, } +//############################################################################# +OsiSolverInterface * +MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, + bool newOsi, const double *lpSol) +{ + /** Create pessimistic risk function model with fixed upper-level vars **/ + + std::string feasCheckSolver(model_->MibSPar()->entry + (MibSParams::feasCheckSolver)); + double targetGap(model_->MibSPar()->entry(MibSParams::slTargetGap)); + + OsiSolverInterface * nSolver; + + int i(0), j(0), idx(0); + double value(0.0), objUb(0.0); + + int uCols(model_->getUpperDim()); + int lCols(model_->getLowerDim()); + int uRows(model_->getOrigUpperRowNum()); + int lRows(model_->getLowerRowNum()); + int rowNum(model_->getNumOrigCons() + 1); + int colNum(model_->getNumOrigVars()); + int *uColIndices(model_->getUpperColInd()); + int *fixedInd(model_->getFixedInd()); + double gap = (targetGap < model_->etol_) ? 0.0 : targetGap; // YX: added SL gap + + // YX: d2^y <= phi(A^x) + \delta * abs(phi(A^x)) + if(objValLL > 0){ // YX: change to abs value later + objUb = objValLL + objValLL* gap/100; + }else{ + objUb = objValLL - objValLL* gap/100; + } + + if (!lpSol){ + lpSol = oSolver->getColSolution(); + } + + if(newOsi){ + double objSense(model_->getLowerObjSense()); + double *lObjCoeffs(model_->getLowerObjCoeffs()); + int *lColIndices(model_->getLowerColInd()); + int *uRowIndices(model_->getOrigUpperRowInd()); + int *lRowIndices(model_->getLowerRowInd()); + const double *origColLb(model_->getOrigColLb()); + const double *origColUb(model_->getOrigColUb()); + const double *origRowLb(model_->getOrigRowLb()); + const double *origRowUb(model_->getOrigRowUb()); + const double *uObjCoeffs(oSolver->getObjCoefficients()); + char *origRowSense(model_->getOrigRowSense()); + + int tmpRowNum(rowNum -1); + int intCnt(0); + double uObjSense(1); + + CoinPackedMatrix matrix = *model_->getOrigConstCoefMatrix(); + matrix.reverseOrdering(); + CoinPackedMatrix *matrixA1 = model_->getA1Matrix(); + CoinPackedMatrix *matrixG1 = model_->getG1Matrix(); + CoinPackedMatrix *matrixA2 = model_->getA2Matrix(); + CoinPackedMatrix *matrixG2 = model_->getG2Matrix(); + CoinShallowPackedVector origRow; + CoinPackedVector row; + CoinPackedVector addedRow; + CoinPackedMatrix *newMat = new CoinPackedMatrix(false, 0, 0); + newMat->setDimensions(0, lCols); + + double *rowUb = new double[rowNum]; + double *rowLb = new double[rowNum]; + double *colUb = new double[lCols]; + double *colLb = new double[lCols]; + int *integerVars = new int[lCols]; + double *objCoeffs = new double[lCols]; + double *newRow = new double[lCols]; + double *uLpSol = new double[uCols]; + double *multA1Xlp = new double[uRows]; + double *multA2Xlp = new double[lRows]; + + if(feasCheckSolver == "Cbc"){ + nSolver = new OsiCbcSolverInterface(); + }else if(feasCheckSolver == "SYMPHONY"){ +#ifdef COIN_HAS_SYMPHONY + nSolver = new OsiSymSolverInterface(); +#else + throw CoinError("SYMPHONY chosen as solver, but it has not been enabled", + "setUpPesModel", "MibsBilevel"); +#endif + }else if(feasCheckSolver == "CPLEX"){ +#ifdef COIN_HAS_CPLEX + nSolver = new OsiCpxSolverInterface(); +#else + throw CoinError("CPLEX chosen as solver, but it has not been enabled", + "setUpPesModel", "MibsBilevel"); +#endif + }else{ + throw CoinError("Unknown solver chosen", + "setUpPesModel", "MibsBilevel"); + } + + CoinZeroN(rowUb, rowNum); + CoinZeroN(rowLb, rowNum); + CoinZeroN(colUb, lCols); + CoinZeroN(colLb, lCols); + CoinZeroN(integerVars, lCols); + CoinZeroN(objCoeffs, lCols); + CoinZeroN(newRow, lCols); + CoinZeroN(uLpSol, uCols); + + /** Obtain A2x matrix **/ + for(i = 0; i < uCols; i++){ + idx = uRowIndices[i]; + uLpSol[i] = lpSol[idx]; + } + matrixA1->times(uLpSol, multA1Xlp); // YX: NEED TO RESET DIM + matrixA2->times(uLpSol, multA2Xlp); // YX: NEED TO RESET DIM + + /** Set row bounds **/ + for(i = 0; i < uRows; i++){ + idx = uRowIndices[i]; + if(origRowSense[i] == 'L'){ + rowLb[i] = origRowLb[idx]; + rowUb[i] = origRowUb[idx] - multA1Xlp[i]; + }else if(origRowSense[i] == 'G'){ + rowLb[i] = origRowLb[idx] - multA1Xlp[i]; + rowUb[i] = origRowUb[idx]; + } + } + for(i = 0; i < lRows; i++){ + idx = lRowIndices[i]; + if(origRowSense[i] == 'L'){ + rowLb[i+uRows] = origRowLb[idx]; + rowUb[i+uRows] = origRowUb[idx] - multA2Xlp[i]; + }else if(origRowSense[i] == 'G'){ + rowLb[i] = origRowLb[idx] - multA2Xlp[i]; + rowUb[i] = origRowUb[idx]; + } + } + + /** Set col bounds **/ + // YX: use lColLbInLProb/lColUbInLProb? + for(i = 0; i < lCols; i++){ + colLb[i] = origColLb[lColIndices[i]]; + colUb[i] = origColUb[lColIndices[i]]; + } + + /** Set up new coeffient matrix with G1 and G2 **/ + for(i = 0; i < rowNum-1; i++){ + if(i < uRows){ + row = matrixG1->getVector(i); + newMat->appendRow(row); + }else{ + row = matrixG1->getVector(i-uRows); + newMat->appendRow(row); + } + } + + /** Set up the value function constraint **/ + CoinDisjointCopyN(lObjCoeffs, lCols, newRow); + + for(i = 0; i < lCols; i++){ + newRow[i] = newRow[i] * objSense; + } + + rowUb[rowNum-1] = objUb; + rowLb[rowNum-1] = -1 * (oSolver->getInfinity()); + + for(i = 0; i < lCols; i++){ + // index1 = lColIndices[i]; + // row.insert(index1, newRow[i]); + row.insert(i, newRow[i]); + } + newMat->appendRow(row); + + + /** Fill in array of integer variables **/ + for(i = 0; i < lCols; i++){ + idx = lColIndices[i]; + if(oSolver->isInteger(idx)){ + integerVars[intCnt] = i; + intCnt ++; + } + // also set up obj: pessimistic risk function + objCoeffs[i] = uObjCoeffs[idx]; + } + + nSolver->loadProblem(*newMat, colLb, colUb, + objCoeffs, rowLb, rowUb); + + for(i = 0; i < intCnt; i++){ + nSolver->setInteger(integerVars[i]); + } + + nSolver->setObjSense(-1 * uObjSense); //1 min; -1 max // YX: pessimistic case + + // YX: still need this? + nSolver->setHintParam(OsiDoReducePrint, true, OsiHintDo); + + delete [] rowUb; + delete [] rowLb; + delete [] colUb; + delete [] colLb; + delete [] objCoeffs; + delete [] integerVars; + delete [] newRow; + delete [] uLpSol; + delete [] multA1Xlp; + delete [] multA2Xlp; + delete newMat; + + }else{ + nSolver = pSolver_; // YX: pessimistic case + nSolver->setRowUpper(rowNum-1, objUb); // YX: gap added to phi(A^2x) + for(i = 0; i < uCols; i++){ + idx = uColIndices[i]; + if(fixedInd[idx] == 1){ + value = floor(lpSol[idx] + 0.5); + nSolver->setColLower(idx, value); + nSolver->setColUpper(idx, value); + } + } + } + + return nSolver; +} + //############################################################################# OsiSolverInterface * MibSBilevel::setUpModel(OsiSolverInterface * oSolver, bool newOsi, @@ -1362,6 +1744,32 @@ MibSBilevel::getLowerObj(const double * sol, double objSense) return objVal * objSense; +} +//############################################################################# +double +MibSBilevel::getRiskFuncVal(OsiSolverInterface * oSolver, double * lowerSol, + bool pesType) +{ + int lN(model_->getLowerDim()); + int *lColIndices(model_->getLowerColInd()); + const double uObjSense(oSolver->getObjSense()); + const double * uObjCoeffs(oSolver->getObjCoefficients()); + + int i(0), idx(0); + double rfVal(0.0); + double rfSense = pesType? -1.0 : +1.0; + + for(i = 0; i < lN; i++){ + idx = lColIndices[i]; + if(1){ + std::cout << "sol[" << i << "]: " << lowerSol[i] << std::endl; + std::cout << "uObjCoeff[" << idx << "]: " << uObjCoeffs[idx] << std::endl; + } + rfVal += uObjCoeffs[idx] * lowerSol[i]; + } + + return rfVal * uObjSense * rfSense; + } //############################################################################# void @@ -1387,8 +1795,10 @@ void tagInSeenLinkingPool_ = solTag; linkingSolution.lowerSolution.push_back(0); + linkingSolution.pesSolution.push_back(0); // YX: pessimistic case linkingSolution.UBSolution.push_back(0); linkingSolution.lowerSolution.clear(); + linkingSolution.pesSolution.clear(); // YX: pessimistic case linkingSolution.UBSolution.clear(); switch(solTag){ @@ -1396,8 +1806,10 @@ void { linkingSolution.tag = solType; linkingSolution.lowerObjValue = 0.0; + linkingSolution.pesRFValue = 0.0; linkingSolution.UBObjValue = 0.0; linkingSolution.lowerSolution.push_back(0); + linkingSolution.pesSolution.push_back(0); linkingSolution.UBSolution.push_back(0); model_->seenLinkingSolutions[linkSol] = linkingSolution; break; @@ -1406,12 +1818,27 @@ void { linkingSolution.tag = solType; linkingSolution.lowerObjValue = objValue; + linkingSolution.pesRFValue = 0.0; linkingSolution.UBObjValue = 0.0; linkingSolution.lowerSolution = shouldStoreValues; + linkingSolution.pesSolution.push_back(0); linkingSolution.UBSolution.push_back(0); model_->seenLinkingSolutions[linkSol] = linkingSolution; break; } + case MibSLinkingPoolTagPesIsInfeasible: // YX: pessimistic case; combined to one using provenoptimal? + { + model_->seenLinkingSolutions[linkSol].tag = MibSLinkingPoolTagPesIsInfeasible; + break; + } + case MibSLinkingPoolTagPesIsFeasible: // YX: pessimistic case + { + model_->seenLinkingSolutions[linkSol].tag = MibSLinkingPoolTagPesIsFeasible; + model_->seenLinkingSolutions[linkSol].pesRFValue = objValue; + model_->seenLinkingSolutions[linkSol].pesSolution.clear(); + model_->seenLinkingSolutions[linkSol].pesSolution = shouldStoreValues; + break; + } case MibSLinkingPoolTagUBIsSolved: { model_->seenLinkingSolutions[linkSol].tag = MibSLinkingPoolTagUBIsSolved; diff --git a/src/MibSBilevel.hpp b/src/MibSBilevel.hpp index f017858f..cf035f41 100644 --- a/src/MibSBilevel.hpp +++ b/src/MibSBilevel.hpp @@ -73,6 +73,7 @@ class MibSBilevel { MibSModel *model_; MibSHeuristic *heuristic_; OsiSolverInterface * lSolver_; + OsiSolverInterface * pSolver_; // YX: pessimistic case OsiSolverInterface * UBSolver_; CoinWarmStart * ws_; @@ -94,6 +95,7 @@ class MibSBilevel { model_ = 0; heuristic_= 0; lSolver_ = 0; + pSolver_ = 0; // YX: pessimistic case UBSolver_ = 0; ws_ = 0; } @@ -112,9 +114,12 @@ class MibSBilevel { int findIndex(int index, int size, int * indices); OsiSolverInterface * setUpUBModel(OsiSolverInterface * solver, double objValLL, bool newOsi, const double *sol = NULL); + OsiSolverInterface * setUpPesModel(OsiSolverInterface * solver, double objValLL, + bool newOsi, const double *sol = NULL); // YX: pessimistic case OsiSolverInterface * setUpModel(OsiSolverInterface * solver, bool newOsi, const double *sol = NULL); double getLowerObj(const double * sol, double objSense); + double getRiskFuncVal(OsiSolverInterface * solver, double * lowerSol, bool pesType); // YX: pessimistic case int binarySearch(int index,int start, int stop, int * indexArray); CoinWarmStart * getWarmStart() {return ws_;} void setWarmStart(CoinWarmStart * ws) {ws_ = ws;} From 46a8f2403a5bbda28649251918b858e2fdeb8819 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Wed, 15 Jun 2022 01:19:06 -0400 Subject: [PATCH 08/26] modify pessimistic functions in MibSBilevel --- src/MibSBilevel.cpp | 106 ++++++++++++++++++++++++++++++-------------- src/MibSBilevel.hpp | 5 ++- 2 files changed, 77 insertions(+), 34 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index 4963b35b..cd63dcec 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -309,7 +309,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) CoinFillN(pesSol, lN + uN, 0.0); std::vector shouldStoreValuesLowerSol(lN); - std::vector shouldStoreValuesPesSol(lN + uN); // YX: pessimistic case + std::vector shouldStoreValuesPesSol(lN); // YX: pessimistic case std::vector shouldStoreValuesUBSol(lN + uN); const double * sol = model_->solver()->getColSolution(); @@ -596,7 +596,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) pesRF = pSolver->getObjValue(); valuesPes = pSolver->getColSolution(); - for (i = 0; i < uN + lN; i++){ + for (i = 0; i < lN; i++){ if ((pSolver->isInteger(i)) && (((valuesPes[i] - floor(valuesPes[i])) < etol) || ((ceil(valuesPes[i]) - valuesPes[i]) < etol))){ @@ -607,7 +607,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) } if(useLinkingSolutionPool && isLinkVarsIntegral_){ - std::copy(pesSol, pesSol + lN + uN, shouldStoreValuesPesSol.begin()); + std::copy(pesSol, pesSol + lN, shouldStoreValuesPesSol.begin()); addSolutionToSeenLinkingSolutionPool (MibSLinkingPoolTagPesIsFeasible, shouldStoreValuesPesSol, pesRF); shouldStoreValuesPesSol.clear(); @@ -636,7 +636,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) if(findPesSol){ // YX: pessimistic case pesRF = model_->seenLinkingSolutions[linkSol].pesRFValue; std::copy(model_->seenLinkingSolutions[linkSol].pesSolution.begin(), - model_->seenLinkingSolutions[linkSol].pesSolution.end(), pesSol); + model_->seenLinkingSolutions[linkSol].pesSolution.end(), pesSol); } } lowerObj = getLowerObj(sol, model_->getLowerObjSense()); @@ -664,22 +664,25 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) } }else{ // YX: find pessimistic risk function value using (LR) and (SL-MILP) solutions - lpPesVal = getRiskFuncVal(model_->getSolver(), lowerSolutionOrd_, findPesSol); + lpPesVal = getRiskFuncVal(lowerSolutionOrd_, findPesSol); // YX: whether (LR) solution is a pessimistic solution if(abs(lpPesVal - pesRF) < etol){ useBilevelBranching_ = false; shouldPrune_ = true; storeSol = MibSRelaxationSol; }else{ - lowerPesVal = getRiskFuncVal(model_->getSolver(), lowerSol, findPesSol); + lowerPesVal = getRiskFuncVal(lowerSol, findPesSol); // YX: whether (SL-MILP) solution is a pessimistic solution + // The upper level feasiblity will be verified later in MibSModel if(abs(lowerPesVal - pesRF) < etol){ memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); if(isUpperIntegral_){ storeSol = MibSHeurSol; } }else{ - // copy pessimistic solution to optLower; + // copy pessimistic solution to optLowerSolOrd; use in cut generation + // no need to verify upper-level feasibility + // may give it another tag; only used in feasibility check memcpy(optLowerSolutionOrd_, pesSol, sizeof(double) * lN); if(isUpperIntegral_){ storeSol = MibSHeurSol; @@ -696,28 +699,33 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) ((computeBestUBWhenXVarsInt == PARAM_ON) && (isUpperIntegral_)) || ((computeBestUBWhenLVarsInt == PARAM_ON) && (isLinkVarsIntegral_)) || ((computeBestUBWhenLVarsFixed == PARAM_ON) && (isLinkVarsFixed_)))){ - if(!findPesSol){ - if(UBSolver_){ - UBSolver_ = setUpUBModel(model_->getSolver(), objVal, false); - } - else{ - UBSolver_ =setUpUBModel(model_->getSolver(), objVal, true); - } - }else{ // YX: set up later - if(UBSolver_){ - UBSolver_ = setUpUBModel(model_->getSolver(), objVal, false); - } - else{ - UBSolver_ =setUpUBModel(model_->getSolver(), objVal, true); - } - // if((sizeFixedInd - uN) < etol){ + if(findPesSol){ + isUBSolved_ = true; // YX: and always feasible in either case; + if((sizeFixedInd - uN) < etol){ // all variables are linking variables; solution fixed - // use and store existing bound, update tag, then jump to term; - // }else{ + if(useLinkingSolutionPool && isLinkVarsIntegral_){ + // YX: use zero as a place holder, then jump to store solution; + // YX: UBSolution is not used anywhere; + objVal = getUpperObj(pesSol, findPesSol); + shouldStoreValuesUBSol.push_back(0); + } + }else{ // exists non-fixed upper vars, set up solver in a special way - // add another function to set up pes UB - // } + // add another function to set up pes UB + // if(UBSolver_){ + // UBSolver_ = setUpPesUBModel(objVal, false); + // }else{ + // UBSolver_ = setUpPesUBModel(objVal, true); + // } + } + }else{ + + if(UBSolver_){ + UBSolver_ = setUpUBModel(model_->getSolver(), objVal, false); + }else{ + UBSolver_ = setUpUBModel(model_->getSolver(), objVal, true); } + OsiSolverInterface * UBSolver = UBSolver_; remainingTime = timeLimit - model_->broker_->subTreeTimer().getTime(); @@ -823,14 +831,16 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) isProvenOptimal_ = false; storeSol = MibSNoSol; } + + } //step 22 //Adding x_L to set E if(useLinkingSolutionPool && isLinkVarsIntegral_){ addSolutionToSeenLinkingSolutionPool (MibSLinkingPoolTagUBIsSolved, shouldStoreValuesUBSol, objVal); + shouldStoreValuesUBSol.clear(); } - shouldStoreValuesUBSol.clear(); - + //step 23 if(isLinkVarsFixed_){ useBilevelBranching_ = false; @@ -1745,15 +1755,15 @@ MibSBilevel::getLowerObj(const double * sol, double objSense) return objVal * objSense; } + //############################################################################# double -MibSBilevel::getRiskFuncVal(OsiSolverInterface * oSolver, double * lowerSol, - bool pesType) +MibSBilevel::getRiskFuncVal(double *lowerSol, bool pesType) { int lN(model_->getLowerDim()); int *lColIndices(model_->getLowerColInd()); - const double uObjSense(oSolver->getObjSense()); - const double * uObjCoeffs(oSolver->getObjCoefficients()); + const double uObjSense(model_->solver()->getObjSense()); + const double *uObjCoeffs(model_->solver()->getObjCoefficients()); int i(0), idx(0); double rfVal(0.0); @@ -1762,7 +1772,7 @@ MibSBilevel::getRiskFuncVal(OsiSolverInterface * oSolver, double * lowerSol, for(i = 0; i < lN; i++){ idx = lColIndices[i]; if(1){ - std::cout << "sol[" << i << "]: " << lowerSol[i] << std::endl; + std::cout << "lowerSol[" << i << "]: " << lowerSol[i] << std::endl; std::cout << "uObjCoeff[" << idx << "]: " << uObjCoeffs[idx] << std::endl; } rfVal += uObjCoeffs[idx] * lowerSol[i]; @@ -1771,6 +1781,36 @@ MibSBilevel::getRiskFuncVal(OsiSolverInterface * oSolver, double * lowerSol, return rfVal * uObjSense * rfSense; } + +//############################################################################# +double +MibSBilevel::getUpperObj(double *lowerSol, bool pesType) +{ + /** calculate upper-level objective value for (x, hat{y}) or (x, y') **/ + int uN(model_->getUpperDim()); + int *uColIndices(model_->getUpperColInd()); + const double uObjSense(model_->solver()->getObjSense()); + const double *uObjCoeffs(model_->solver()->getObjCoefficients()); + const double *lrSol(model_->solver()->getColSolution()); + + int i(0), idx(0); + double rfuncVal(0.0), uObjVal(0.0); + + for(i = 0; i < uN; i++){ + idx = uColIndices[i]; + if(1){ + std::cout << "lrSol[" << i << "]: " << lrSol[i] << std::endl; + std::cout << "uObjCoeff[" << idx << "]: " << uObjCoeffs[idx] << std::endl; + } + uObjVal += uObjCoeffs[idx] * lrSol[i] * uObjSense; + } + + rfuncVal = getRiskFuncVal(lowerSol, pesType); + uObjVal += rfuncVal; + + return uObjVal; +} + //############################################################################# void MibSBilevel::addSolutionToSeenLinkingSolutionPool(MibSLinkingPoolTag solTag, std::vector diff --git a/src/MibSBilevel.hpp b/src/MibSBilevel.hpp index cf035f41..829c3388 100644 --- a/src/MibSBilevel.hpp +++ b/src/MibSBilevel.hpp @@ -114,12 +114,15 @@ class MibSBilevel { int findIndex(int index, int size, int * indices); OsiSolverInterface * setUpUBModel(OsiSolverInterface * solver, double objValLL, bool newOsi, const double *sol = NULL); + OsiSolverInterface * setUpPesUBModel(OsiSolverInterface * solver, double objValLL, + bool newOsi, const double *sol = NULL); // YX: pessimistic case OsiSolverInterface * setUpPesModel(OsiSolverInterface * solver, double objValLL, bool newOsi, const double *sol = NULL); // YX: pessimistic case OsiSolverInterface * setUpModel(OsiSolverInterface * solver, bool newOsi, const double *sol = NULL); double getLowerObj(const double * sol, double objSense); - double getRiskFuncVal(OsiSolverInterface * solver, double * lowerSol, bool pesType); // YX: pessimistic case + double getRiskFuncVal(double *lowerSol, bool pesType); // YX: pessimistic case + double getUpperObj(double *lowerSol, bool pesType); // YX: pessimistic case int binarySearch(int index,int start, int stop, int * indexArray); CoinWarmStart * getWarmStart() {return ws_;} void setWarmStart(CoinWarmStart * ws) {ws_ = ws;} From f8ce72ff1a8d5a13e571431b2e2e562026619f73 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Wed, 15 Jun 2022 22:33:26 -0400 Subject: [PATCH 09/26] add pess UB function --- src/MibSBilevel.cpp | 199 ++++++++++++++++++++++++++++++++++++++++---- src/MibSBilevel.hpp | 5 +- 2 files changed, 187 insertions(+), 17 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index cd63dcec..68aed60d 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -832,7 +832,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) storeSol = MibSNoSol; } - } + } // end of (findPesSol) //step 22 //Adding x_L to set E if(useLinkingSolutionPool && isLinkVarsIntegral_){ @@ -875,6 +875,177 @@ MibSBilevel::gutsOfDestructor() } +//############################################################################# +OsiSolverInterface * +MibSBilevel::setUpPesUBModel(double *lowerSol, bool newOsi) +{ + /** Set up upper bnding problem with x_L after pessimistic MILP is solved **/ + std::string feasCheckSolver = + model_->MibSPar()->entry(MibSParams::feasCheckSolver); + + OsiSolverInterface *nSolver; + + int i(0), idx(0); + double value(0.0); + + int uCols(model_->getUpperDim()); + int uRows(model_->getOrigUpperRowNum()); + int lCols(model_->getLowerDim()); + int *uRowIndices(model_->getOrigUpperRowInd()); + int *uColIndices(model_->getUpperColInd()); + int *lColIndices(model_->getLowerColInd()); + int *fixedInd(model_->getFixedInd()); + char *origRowSense(model_->getOrigRowSense()); + + const double *origRowLb(model_->getOrigRowLb()); + const double *origRowUb(model_->getOrigRowUb()); + const double *lpSol(model_->getSolver()->getColSolution()); + + CoinPackedMatrix *matrixG1 = model_->getG1Matrix(); + + double *lwSol = new double[lCols]; + double *multG1y = new double[uRows]; + + if(newOsi){ + double uObjSense(1); + const double *origColLb(model_->getOrigColLb()); + const double *origColUb(model_->getOrigColUb()); + const double *uObjCoeffs(model_->getSolver()->getObjCoefficients()); + CoinPackedMatrix *matrixA1 = model_->getA1Matrix(); + + double *rowUb = new double[uRows]; + double *rowLb = new double[uRows]; + double *colUb = new double[uCols]; + double *colLb = new double[uCols]; + + double *multG1y = new double[uRows]; + double *objCoeffs = new double[uCols]; + int *integerVars = new int[uCols]; + + CoinZeroN(colUb, uCols); + CoinZeroN(colLb, uCols); + CoinZeroN(integerVars, uCols); + CoinZeroN(objCoeffs, uCols); + + int intCnt(0); + + if (feasCheckSolver == "Cbc"){ + nSolver = new OsiCbcSolverInterface(); + }else if (feasCheckSolver == "SYMPHONY"){ +#ifdef COIN_HAS_SYMPHONY + nSolver = new OsiSymSolverInterface(); +#else + throw CoinError("SYMPHONY chosen as solver, but it has not been enabled", + "setUpUBModel", "MibsBilevel"); +#endif + }else if (feasCheckSolver == "CPLEX"){ +#ifdef COIN_HAS_CPLEX + nSolver = new OsiCpxSolverInterface(); +#else + throw CoinError("CPLEX chosen as solver, but it has not been enabled", + "setUpUBModel", "MibsBilevel"); +#endif + }else{ + throw CoinError("Unknown solver chosen", + "setUpUBModel", "MibsBilevel"); + } + + /** Obtain G1\hat{y} matrix, where \hat{y} is a pess (heur) sol **/ + for(i = 0; i < lCols; i++){ + idx = lColIndices[i]; + lwSol[i] = lowerSol[idx]; + } + matrixG1->times(lwSol, multG1y); + + /** Set row bounds **/ + for(i = 0; i < uRows; i++){ + idx = uRowIndices[i]; + if(origRowSense[idx] == 'L'){ + rowLb[i] = origRowLb[idx]; + rowUb[i] = origRowUb[idx] - multG1y[i]; + }else if(origRowSense[idx] == 'G'){ + rowLb[i] = origRowLb[idx] - multG1y[i]; + rowUb[i] = origRowUb[idx]; + } + } + + /** Set col bounds **/ + for(i = 0; i < uCols; i++){ + idx = uColIndices[i]; + if(fixedInd[idx] == 1){ + colLb[i] = floor(lpSol[idx] + 0.5); + colUb[i] = colLb[idx]; + }else{ + colLb[i] = origColLb[idx]; + colUb[i] = origColUb[idx]; + } + } + + /** Fill in array of integer variables and set up obj coeffs**/ + for(i = 0; i < uCols; i++){ + if(model_->getSolver()->isInteger(i)){ + integerVars[intCnt] = i; + intCnt ++; + } + // also set up obj: pessimistic risk function + objCoeffs[i] = uObjCoeffs[idx]; + } + + nSolver->loadProblem(*matrixA1, colLb, colUb, + objCoeffs, rowLb, rowUb); + + for(i = 0; i < intCnt; i++){ + nSolver->setInteger(integerVars[i]); + } + + nSolver->setObjSense(uObjSense); //1 min; -1 max + + nSolver->setHintParam(OsiDoReducePrint, true, OsiHintDo); + + delete [] integerVars; + delete [] rowUb; + delete [] rowLb; + delete [] colUb; + delete [] colLb; + delete [] objCoeffs; + + }else{ + nSolver = UBSolver_; + + /** Obtain G1\hat{y} matrix, where \hat{y} is a pess (heur) sol **/ + for(i = 0; i < lCols; i++){ + idx = lColIndices[i]; + lwSol[i] = lowerSol[idx]; + } + matrixG1->times(lwSol, multG1y); + + /** Update row bounds **/ + for(i = 0; i < uRows; i++){ + idx = uRowIndices[i]; + if(origRowSense[idx] == 'L'){ + nSolver->setRowUpper(i, origRowUb[idx] - multG1y[i]); + }else if(origRowSense[idx] == 'G'){ + nSolver->setRowLower(i, origRowLb[idx] - multG1y[i]); + } + } + + /** Update column bounds **/ + for(i = 0; i < uCols; i++){ + idx = uColIndices[i]; + if(fixedInd[idx] == 1){ + value = floor(lpSol[idx] + 0.5); + nSolver->setColLower(idx, value); + nSolver->setColUpper(idx, value); + } + } + } + + delete [] lwSol; + delete [] multG1y; + + return nSolver; +} + //############################################################################# OsiSolverInterface * MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, @@ -1161,29 +1332,29 @@ MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, /** Obtain A2x matrix **/ for(i = 0; i < uCols; i++){ - idx = uRowIndices[i]; + idx = uColIndices[i]; uLpSol[i] = lpSol[idx]; } - matrixA1->times(uLpSol, multA1Xlp); // YX: NEED TO RESET DIM - matrixA2->times(uLpSol, multA2Xlp); // YX: NEED TO RESET DIM + matrixA1->times(uLpSol, multA1Xlp); + matrixA2->times(uLpSol, multA2Xlp); /** Set row bounds **/ for(i = 0; i < uRows; i++){ idx = uRowIndices[i]; - if(origRowSense[i] == 'L'){ + if(origRowSense[idx] == 'L'){ rowLb[i] = origRowLb[idx]; rowUb[i] = origRowUb[idx] - multA1Xlp[i]; - }else if(origRowSense[i] == 'G'){ + }else if(origRowSense[idx] == 'G'){ rowLb[i] = origRowLb[idx] - multA1Xlp[i]; rowUb[i] = origRowUb[idx]; } } for(i = 0; i < lRows; i++){ idx = lRowIndices[i]; - if(origRowSense[i] == 'L'){ + if(origRowSense[idx] == 'L'){ rowLb[i+uRows] = origRowLb[idx]; rowUb[i+uRows] = origRowUb[idx] - multA2Xlp[i]; - }else if(origRowSense[i] == 'G'){ + }else if(origRowSense[idx] == 'G'){ rowLb[i] = origRowLb[idx] - multA2Xlp[i]; rowUb[i] = origRowUb[idx]; } @@ -1202,7 +1373,7 @@ MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, row = matrixG1->getVector(i); newMat->appendRow(row); }else{ - row = matrixG1->getVector(i-uRows); + row = matrixG2->getVector(i-uRows); newMat->appendRow(row); } } @@ -1218,9 +1389,9 @@ MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, rowLb[rowNum-1] = -1 * (oSolver->getInfinity()); for(i = 0; i < lCols; i++){ - // index1 = lColIndices[i]; - // row.insert(index1, newRow[i]); - row.insert(i, newRow[i]); + idx = lColIndices[i]; + row.insert(idx, newRow[i]); + // row.insert(i, newRow[i]); } newMat->appendRow(row); @@ -1229,8 +1400,8 @@ MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, for(i = 0; i < lCols; i++){ idx = lColIndices[i]; if(oSolver->isInteger(idx)){ - integerVars[intCnt] = i; - intCnt ++; + integerVars[intCnt] = i; + intCnt ++; } // also set up obj: pessimistic risk function objCoeffs[i] = uObjCoeffs[idx]; diff --git a/src/MibSBilevel.hpp b/src/MibSBilevel.hpp index 829c3388..893b98e2 100644 --- a/src/MibSBilevel.hpp +++ b/src/MibSBilevel.hpp @@ -111,15 +111,14 @@ class MibSBilevel { private: - int findIndex(int index, int size, int * indices); + OsiSolverInterface * setUpPesUBModel(double *lowerSol, bool newOsi); // YX: pessimistic case OsiSolverInterface * setUpUBModel(OsiSolverInterface * solver, double objValLL, bool newOsi, const double *sol = NULL); - OsiSolverInterface * setUpPesUBModel(OsiSolverInterface * solver, double objValLL, - bool newOsi, const double *sol = NULL); // YX: pessimistic case OsiSolverInterface * setUpPesModel(OsiSolverInterface * solver, double objValLL, bool newOsi, const double *sol = NULL); // YX: pessimistic case OsiSolverInterface * setUpModel(OsiSolverInterface * solver, bool newOsi, const double *sol = NULL); + int findIndex(int index, int size, int * indices); // YX: not used? double getLowerObj(const double * sol, double objSense); double getRiskFuncVal(double *lowerSol, bool pesType); // YX: pessimistic case double getUpperObj(double *lowerSol, bool pesType); // YX: pessimistic case From d909bc90f4101fa6501ec8364ff11eaa99e7f0d1 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Wed, 15 Jun 2022 22:57:20 -0400 Subject: [PATCH 10/26] change gap expressions using fabs() --- src/MibSBilevel.cpp | 20 ++++++-------------- src/MibSCutGenerator.cpp | 40 +++++++++++----------------------------- 2 files changed, 17 insertions(+), 43 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index 68aed60d..746bd481 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -1065,7 +1065,7 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, int * fixedInd = model_->fixedInd_; int i(0), j(0), index1(0); - double value(0.0); + double value(0.0), objUb(0.0); int uCols(model_->getUpperDim()); int uRows(model_->getOrigUpperRowNum()); @@ -1074,14 +1074,10 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, int rowNum(model_->getNumOrigCons() + 1); int colNum(model_->getNumOrigVars()); int * uColIndices(model_->getUpperColInd()); - double gap = (targetGap < model_->etol_) ? 0.0 : targetGap; // YX: added SL gap - double objUb(0.0); + double gap = (targetGap < model_->getTolerance()) ? 0.0 : targetGap; // YX: added SL gap + // YX: d2^y <= phi(A^x) + \delta * abs(phi(A^x)) - if(objValLL > 0){ - objUb = objValLL + objValLL* gap/100; - }else{ - objUb = objValLL - objValLL* gap/100; - } + objUb = objValLL + fabs(objValLL) * gap/100; if(newOsi){ double objSense(model_->getLowerObjSense()); @@ -1247,14 +1243,10 @@ MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, int colNum(model_->getNumOrigVars()); int *uColIndices(model_->getUpperColInd()); int *fixedInd(model_->getFixedInd()); - double gap = (targetGap < model_->etol_) ? 0.0 : targetGap; // YX: added SL gap + double gap = (targetGap < model_->getTolerance()) ? 0.0 : targetGap; // YX: added SL gap // YX: d2^y <= phi(A^x) + \delta * abs(phi(A^x)) - if(objValLL > 0){ // YX: change to abs value later - objUb = objValLL + objValLL* gap/100; - }else{ - objUb = objValLL - objValLL* gap/100; - } + objUb = objValLL + fabs(objValLL) * gap/100; if (!lpSol){ lpSol = oSolver->getColSolution(); diff --git a/src/MibSCutGenerator.cpp b/src/MibSCutGenerator.cpp index 546ef7a1..aa86c41f 100644 --- a/src/MibSCutGenerator.cpp +++ b/src/MibSCutGenerator.cpp @@ -686,7 +686,7 @@ MibSCutGenerator::intersectionCuts(BcpsConstraintPool &conPool, mult = -1; break; case 'E': - std::cout + std::couta << "MibS cannot currently handle equality constraints." << std::endl; abort(); break; @@ -822,21 +822,21 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, int whichCutsLL(localModel_->MibSPar_->entry (MibSParams::whichCutsLL)); double targetGap(localModel_->MibSPar_->entry(MibSParams::slTargetGap)); - double etol(localModel_->etol_); - double gap = (targetGap < etol) ? 0.0 : targetGap; // YX: added SL gap - double templObj(0.0); // YX: for nonzero gap, track d^2y^* bool foundSolution(false); double timeLimit(localModel_->AlpsPar()->entry(AlpsParams::timeLimit)); double remainingTime(0.0); - bool getA2Matrix(false), getG2Matrix(false); + // bool getA2Matrix(false), getG2Matrix(false); OsiSolverInterface * oSolver = localModel_->solver(); double infinity(oSolver->getInfinity()); int i, j; int index(0), cntA2(0), cntG2(0), cntInt(0); double coef(0.0), lObjVal(0.0), value(0.0); + double templObj(0.0); // YX: for nonzero gap, track d^2y^* + double etol(localModel_->getTolerance()); + double gap = (targetGap < etol) ? 0.0 : targetGap; // YX: added SL gap int numCols(localModel_->getNumCols()); int uCols(localModel_->getUpperDim()); int lCols(localModel_->getLowerDim()); @@ -982,11 +982,7 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, for(i = 0; i < lCols; i++){ templObj += lObjSense * lObjCoeff[i] * optLowerSol[i]; // YX: track d^2y^* } - if(templObj > 0){ - newRowUb[newNumRows-1] = -templObj - (templObj * gap/100); - }else{ - newRowUb[newNumRows-1] = -templObj + (templObj * gap/100); - } + newRowUb[newNumRows-1] = -templObj - (fabs(templObj) * gap/100); } //filling col bounds @@ -1114,9 +1110,6 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, //the optimal solution of relaxation which satisfies integrality requirements //throw CoinError("The MIP which gives the best lower-level sol, cannot be infeasible!", // "findLowerLevelSol", "MibSCutGenerator"); - // if(targetGap > etol){ - // std::cout << "Type2IC aux MILP with optimality gap is infeasible."<getOrigRowSense(); double *lObjCoeffs(localModel_->getLowerObjCoeffs()); double objSense(localModel_->getLowerObjSense()); - bool getA2Matrix(false), getG2Matrix(false); + // bool getA2Matrix(false), getG2Matrix(false); double targetGap(localModel_->MibSPar_->entry(MibSParams::slTargetGap)); double gap = (targetGap < etol) ? 0.0 : targetGap; @@ -1225,13 +1218,9 @@ MibSCutGenerator::getAlphaIC(double** extRay, double* uselessIneqs, templObj += objSense * lObjCoeffs[i] * lowerSolution[i]; // YX: track d^2y^* } - // YX: type I intersection IC; -d^2y^* - gap*|-d^2y^*| + // YX: type I intersection IC; -d^2y^* - gap*|d^2y^*| if ((targetGap > etol) && (!uselessIneqs)){ - if(templObj > 0){ - rhs[lRows] += -templObj * gap/100; - }else{ - rhs[lRows] += templObj * gap/100; - } + rhs[lRows] += -fabs(templObj) * gap/100; } for (i = 0; i < numNonBasic; i++){ @@ -1346,7 +1335,7 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo OsiSolverInterface *oSolver = localModel_->solver(); double infinity(oSolver->getInfinity()); - bool getA2G2Matrix(false), getG2Matrix(false); + // bool getA2G2Matrix(false), getG2Matrix(false); int i, j; int rowIndex(0), colIndex(0), numElements(0), pos(0), cntInt(0); double rhs(0.0), value(0.0); @@ -1540,11 +1529,7 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo rhs += lObjSense * lObjCoeff[i] * (lpSol[colIndex] - optLowerSol[i]); templObj += lObjSense * lObjCoeff[i] * optLowerSol[i]; // YX: track d^2y^* } - if(templObj > 0){ - rhs += -templObj * gap/100; - }else{ - rhs += templObj * gap/100; - } + rhs += -fabs(templObj) * gap/100; nSolver->setRowUpper(newNumRows-1, rhs); } @@ -1633,9 +1618,6 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo //std::cout << "current time = " << timeLimit - localModel_->broker_->subTreeTimer().getTime() << std::endl; //throw CoinError("The MIP which is solved for watermelon IC, cannot be infeasible!", // "findLowerLevelSolWatermelonIC", "MibSCutGenerator"); - // if(targetGap > etol){ - // std::cout << "Watermelon aux MILP with optimality gap is infeasible."< Date: Thu, 16 Jun 2022 00:05:55 -0400 Subject: [PATCH 11/26] update pess feasibility check (UB related) --- src/MibSBilevel.cpp | 81 ++++++++++++++++++++++++++++----------------- 1 file changed, 51 insertions(+), 30 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index 746bd481..3203ed65 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -666,7 +666,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) // YX: find pessimistic risk function value using (LR) and (SL-MILP) solutions lpPesVal = getRiskFuncVal(lowerSolutionOrd_, findPesSol); // YX: whether (LR) solution is a pessimistic solution - if(abs(lpPesVal - pesRF) < etol){ + if(fabs(lpPesVal - pesRF) < etol){ useBilevelBranching_ = false; shouldPrune_ = true; storeSol = MibSRelaxationSol; @@ -674,7 +674,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) lowerPesVal = getRiskFuncVal(lowerSol, findPesSol); // YX: whether (SL-MILP) solution is a pessimistic solution // The upper level feasiblity will be verified later in MibSModel - if(abs(lowerPesVal - pesRF) < etol){ + if(fabs(lowerPesVal - pesRF) < etol){ memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); if(isUpperIntegral_){ storeSol = MibSHeurSol; @@ -699,33 +699,32 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) ((computeBestUBWhenXVarsInt == PARAM_ON) && (isUpperIntegral_)) || ((computeBestUBWhenLVarsInt == PARAM_ON) && (isLinkVarsIntegral_)) || ((computeBestUBWhenLVarsFixed == PARAM_ON) && (isLinkVarsFixed_)))){ - if(findPesSol){ - isUBSolved_ = true; // YX: and always feasible in either case; - if((sizeFixedInd - uN) < etol){ - // all variables are linking variables; solution fixed - if(useLinkingSolutionPool && isLinkVarsIntegral_){ - // YX: use zero as a place holder, then jump to store solution; - // YX: UBSolution is not used anywhere; - objVal = getUpperObj(pesSol, findPesSol); - shouldStoreValuesUBSol.push_back(0); - } - }else{ + + if((findPesSol)&&((sizeFixedInd - uN) < etol)){ + // YX: all variables are linking variables; solution fixed + isUBSolved_ = true; + + // YX: UBSolution is not used anywhere? use zero as a place holder; + objVal = getUpperObj(optLowerSolutionOrd_, findPesSol); + shouldStoreValuesUBSol.push_back(0); + }else{ + // YX: need to solve an upper bounding problem + if(findPesSol){ // exists non-fixed upper vars, set up solver in a special way // add another function to set up pes UB - // if(UBSolver_){ - // UBSolver_ = setUpPesUBModel(objVal, false); - // }else{ - // UBSolver_ = setUpPesUBModel(objVal, true); - // } + if(UBSolver_){ + UBSolver_ = setUpPesUBModel(optLowerSolutionOrd_, false); + }else{ + UBSolver_ = setUpPesUBModel(optLowerSolutionOrd_, true); + } + }else{ + // regular optimistic case + if(UBSolver_){ + UBSolver_ = setUpUBModel(model_->getSolver(), objVal, false); + }else{ + UBSolver_ = setUpUBModel(model_->getSolver(), objVal, true); + } } - }else{ - - if(UBSolver_){ - UBSolver_ = setUpUBModel(model_->getSolver(), objVal, false); - }else{ - UBSolver_ = setUpUBModel(model_->getSolver(), objVal, true); - } - OsiSolverInterface * UBSolver = UBSolver_; remainingTime = timeLimit - model_->broker_->subTreeTimer().getTime(); @@ -800,7 +799,8 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) isUBSolved_ = true; - if (UBSolver->isProvenOptimal()){ + if(findPesSol){ + if (UBSolver->isProvenOptimal()){ isProvenOptimal_ = true; const double * valuesUB = UBSolver->getColSolution(); std::copy(valuesUB, valuesUB + uN + lN, shouldStoreValuesUBSol.begin()); @@ -831,9 +831,30 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) isProvenOptimal_ = false; storeSol = MibSNoSol; } - - } // end of (findPesSol) - //step 22 + }else{ + // YX: pess UB has to have an optimal solution b/c feasible (PES-MILP) + assert(UBSolver->isProvenOptimal()); + isProvenOptimal_ = true; + // YX: copy optimal upper-level solution to containers + const double * valuesUB = UBSolver->getColSolution(); + // DEBUG INDEX!! ARE BOTH UPPER AND LOWER SORTED ALREADY? + // NEED TO RECOMBINE THEM USING INDEX? + // std::copy(valuesUB, valuesUB + uN, shouldStoreValuesUBSol.begin()); + // std::copy(lowerSol, lowerSol + lN, shouldStoreValuesUBSol.end()); + for (i = 0; i < uN; i++){ + if ((UBSolver->isInteger(i)) && + (((valuesUB[i] - floor(valuesUB[i])) < etol) || + ((ceil(valuesUB[i]) - valuesUB[i]) < etol))){ + optUpperSolutionOrd_[i] = (double) floor(valuesUB[i] + 0.5); + }else{ + optUpperSolutionOrd_[i] = (double) valuesUB[i]; + } + } + } // end of (processing an UB solution) + + } // end of (solving an UB problem in either optimistic or pessimistic case) + + //step 22 //Adding x_L to set E if(useLinkingSolutionPool && isLinkVarsIntegral_){ addSolutionToSeenLinkingSolutionPool From f45078b7311ed3dfe16c8c79aeeaffd26492d249 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Wed, 22 Jun 2022 22:16:42 -0400 Subject: [PATCH 12/26] fix a bug in checkBilevelFeas; still has seg fault --- src/MibSBilevel.cpp | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index 3203ed65..8a2781f0 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -289,12 +289,13 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) (MibSParams::useLinkingSolutionPool)); // YX: target optimality gap for bounded rationality - double targetGap(model_->MibSPar_->entry(MibSParams::slTargetGap)); + double targetGap(model_->MibSPar()->entry(MibSParams::slTargetGap)); double timeLimit(model_->AlpsPar()->entry(AlpsParams::timeLimit)); double remainingTime(0.0), startTimeVF(0.0), startTimePES(0.0), startTimeUB(0.0); // YX: pessimistic case MibSSolType storeSol(MibSNoSol); int lN(model_->lowerDim_); // lower-level dimension int uN(model_->upperDim_); // upper-level dimension + int uRows(model_->getOrigUpperRowNum()); int i(0), index(0), length(0), pos(0); int sizeFixedInd(model_->sizeFixedInd_); double etol(model_->etol_), objVal(0.0), lowerObj(0.0); @@ -313,7 +314,8 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) std::vector shouldStoreValuesUBSol(lN + uN); const double * sol = model_->solver()->getColSolution(); - + double gap = (targetGap < etol) ? 0.0 : targetGap; // YX: added SL gap + std::vector linkSol; for(i = 0; i < uN; i++){ index = upperColInd[i]; @@ -502,8 +504,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) if(findPesSol){ if(pSolver_){ pSolver_ = setUpPesModel(model_->getSolver(), objVal, false); - } - else{ + }else{ pSolver_ = setUpPesModel(model_->getSolver(), objVal, true); } @@ -651,7 +652,8 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) /** Current solution is bilevel feasible **/ // YX: [d2^y - phi(A^x)]/abs(phi(A^x)), combined check when gap is set to 0; if(!findPesSol){ - if(((lowerObj - objVal) <= fabs(objVal) * targetGap/100) && (isIntegral_)){ + if(((lowerObj - objVal) <= fabs(objVal) * gap/100) && (isIntegral_)){ + // if(((lowerObj - objVal) < etol) && (isIntegral_)){ LPSolStatus_ = MibSLPSolStatusFeasible; useBilevelBranching_ = false; shouldPrune_ = true; @@ -675,6 +677,8 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) // YX: whether (SL-MILP) solution is a pessimistic solution // The upper level feasiblity will be verified later in MibSModel if(fabs(lowerPesVal - pesRF) < etol){ + // YX: need to verify upper-level feasibility here!!! + // if passed, mark as heurSol; otherwise, always use pesSol; memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); if(isUpperIntegral_){ storeSol = MibSHeurSol; @@ -700,7 +704,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) ((computeBestUBWhenLVarsInt == PARAM_ON) && (isLinkVarsIntegral_)) || ((computeBestUBWhenLVarsFixed == PARAM_ON) && (isLinkVarsFixed_)))){ - if((findPesSol)&&((sizeFixedInd - uN) < etol)){ + if((findPesSol)&&(((sizeFixedInd - uN) < 0)||(uRows > 0))){ // YX: all variables are linking variables; solution fixed isUBSolved_ = true; @@ -799,7 +803,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) isUBSolved_ = true; - if(findPesSol){ + if(!findPesSol){ if (UBSolver->isProvenOptimal()){ isProvenOptimal_ = true; const double * valuesUB = UBSolver->getColSolution(); @@ -837,10 +841,10 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) isProvenOptimal_ = true; // YX: copy optimal upper-level solution to containers const double * valuesUB = UBSolver->getColSolution(); - // DEBUG INDEX!! ARE BOTH UPPER AND LOWER SORTED ALREADY? // NEED TO RECOMBINE THEM USING INDEX? // std::copy(valuesUB, valuesUB + uN, shouldStoreValuesUBSol.begin()); // std::copy(lowerSol, lowerSol + lN, shouldStoreValuesUBSol.end()); + shouldStoreValuesUBSol.push_back(0); // YX: place holder for (i = 0; i < uN; i++){ if ((UBSolver->isInteger(i)) && (((valuesUB[i] - floor(valuesUB[i])) < etol) || @@ -1241,6 +1245,7 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, } //############################################################################# +// YX: CLEAN UP LATER OsiSolverInterface * MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, bool newOsi, const double *lpSol) @@ -1286,19 +1291,19 @@ MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, const double *uObjCoeffs(oSolver->getObjCoefficients()); char *origRowSense(model_->getOrigRowSense()); - int tmpRowNum(rowNum -1); + // int tmpRowNum(rowNum -1); int intCnt(0); double uObjSense(1); - CoinPackedMatrix matrix = *model_->getOrigConstCoefMatrix(); - matrix.reverseOrdering(); + // CoinPackedMatrix matrix = *model_->getOrigConstCoefMatrix(); + // matrix.reverseOrdering(); CoinPackedMatrix *matrixA1 = model_->getA1Matrix(); CoinPackedMatrix *matrixG1 = model_->getG1Matrix(); CoinPackedMatrix *matrixA2 = model_->getA2Matrix(); CoinPackedMatrix *matrixG2 = model_->getG2Matrix(); - CoinShallowPackedVector origRow; + // CoinShallowPackedVector origRow; CoinPackedVector row; - CoinPackedVector addedRow; + // CoinPackedVector addedRow; CoinPackedMatrix *newMat = new CoinPackedMatrix(false, 0, 0); newMat->setDimensions(0, lCols); From aa59b5fe663bdb0a71ad0cc0915531b1efa0dc63 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Thu, 23 Jun 2022 19:00:42 -0400 Subject: [PATCH 13/26] update pessimistic case functions --- src/MibSBilevel.cpp | 181 +++++++++++++++++++++++--------------------- src/MibSBilevel.hpp | 8 +- 2 files changed, 99 insertions(+), 90 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index 8a2781f0..ffe211ee 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -503,9 +503,9 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) // YX: pessimistic case; need to confirm repeated model setup? if(findPesSol){ if(pSolver_){ - pSolver_ = setUpPesModel(model_->getSolver(), objVal, false); + pSolver_ = setUpPesModel(objVal, false); }else{ - pSolver_ = setUpPesModel(model_->getSolver(), objVal, true); + pSolver_ = setUpPesModel(objVal, true); } OsiSolverInterface *pSolver = pSolver_; @@ -666,27 +666,27 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) } }else{ // YX: find pessimistic risk function value using (LR) and (SL-MILP) solutions - lpPesVal = getRiskFuncVal(lowerSolutionOrd_, findPesSol); + lpPesVal = getRiskFuncVal(lowerSolutionOrd_); // YX: whether (LR) solution is a pessimistic solution - if(fabs(lpPesVal - pesRF) < etol){ + if(((lowerObj - objVal) <= fabs(objVal) * gap/100) && (isIntegral_) && + (fabs(lpPesVal - pesRF) < etol)){ useBilevelBranching_ = false; shouldPrune_ = true; storeSol = MibSRelaxationSol; }else{ - lowerPesVal = getRiskFuncVal(lowerSol, findPesSol); + lowerPesVal = getRiskFuncVal(lowerSol); // YX: whether (SL-MILP) solution is a pessimistic solution - // The upper level feasiblity will be verified later in MibSModel if(fabs(lowerPesVal - pesRF) < etol){ - // YX: need to verify upper-level feasibility here!!! + // TODO: VERIFY UPPER LEVEL FEASIBILITY HERE? // if passed, mark as heurSol; otherwise, always use pesSol; - memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); + // memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); + memcpy(optLowerSolutionOrd_, pesSol, sizeof(double) * lN); if(isUpperIntegral_){ storeSol = MibSHeurSol; } }else{ // copy pessimistic solution to optLowerSolOrd; use in cut generation // no need to verify upper-level feasibility - // may give it another tag; only used in feasibility check memcpy(optLowerSolutionOrd_, pesSol, sizeof(double) * lN); if(isUpperIntegral_){ storeSol = MibSHeurSol; @@ -704,25 +704,23 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) ((computeBestUBWhenLVarsInt == PARAM_ON) && (isLinkVarsIntegral_)) || ((computeBestUBWhenLVarsFixed == PARAM_ON) && (isLinkVarsFixed_)))){ - if((findPesSol)&&(((sizeFixedInd - uN) < 0)||(uRows > 0))){ + if((findPesSol)&&(((sizeFixedInd - uN) >= 0)||(uRows == 0))){ // YX: all variables are linking variables; solution fixed isUBSolved_ = true; // YX: UBSolution is not used anywhere? use zero as a place holder; - objVal = getUpperObj(optLowerSolutionOrd_, findPesSol); + objVal = getUpperObj(optLowerSolutionOrd_); shouldStoreValuesUBSol.push_back(0); }else{ - // YX: need to solve an upper bounding problem if(findPesSol){ - // exists non-fixed upper vars, set up solver in a special way - // add another function to set up pes UB + // YX: when existing non-determined upper vars if(UBSolver_){ UBSolver_ = setUpPesUBModel(optLowerSolutionOrd_, false); }else{ UBSolver_ = setUpPesUBModel(optLowerSolutionOrd_, true); } }else{ - // regular optimistic case + // YX: regular optimistic case if(UBSolver_){ UBSolver_ = setUpUBModel(model_->getSolver(), objVal, false); }else{ @@ -926,7 +924,7 @@ MibSBilevel::setUpPesUBModel(double *lowerSol, bool newOsi) const double *origRowUb(model_->getOrigRowUb()); const double *lpSol(model_->getSolver()->getColSolution()); - CoinPackedMatrix *matrixG1 = model_->getG1Matrix(); + CoinPackedMatrix *matrixG1(model_->getG1Matrix()); double *lwSol = new double[lCols]; double *multG1y = new double[uRows]; @@ -936,7 +934,7 @@ MibSBilevel::setUpPesUBModel(double *lowerSol, bool newOsi) const double *origColLb(model_->getOrigColLb()); const double *origColUb(model_->getOrigColUb()); const double *uObjCoeffs(model_->getSolver()->getObjCoefficients()); - CoinPackedMatrix *matrixA1 = model_->getA1Matrix(); + CoinPackedMatrix *matrixA1(model_->getA1Matrix()); double *rowUb = new double[uRows]; double *rowLb = new double[uRows]; @@ -975,7 +973,7 @@ MibSBilevel::setUpPesUBModel(double *lowerSol, bool newOsi) "setUpUBModel", "MibsBilevel"); } - /** Obtain G1\hat{y} matrix, where \hat{y} is a pess (heur) sol **/ + /** Obtain G1\hat{y}: \hat{y} is a (heuristic) pessimistic sol **/ for(i = 0; i < lCols; i++){ idx = lColIndices[i]; lwSol[i] = lowerSol[idx]; @@ -989,8 +987,8 @@ MibSBilevel::setUpPesUBModel(double *lowerSol, bool newOsi) rowLb[i] = origRowLb[idx]; rowUb[i] = origRowUb[idx] - multG1y[i]; }else if(origRowSense[idx] == 'G'){ - rowLb[i] = origRowLb[idx] - multG1y[i]; - rowUb[i] = origRowUb[idx]; + rowLb[i] = -1*origRowUb[idx]; + rowUb[i] = -1*origRowLb[idx] - multG1y[i]; } } @@ -1037,7 +1035,7 @@ MibSBilevel::setUpPesUBModel(double *lowerSol, bool newOsi) }else{ nSolver = UBSolver_; - /** Obtain G1\hat{y} matrix, where \hat{y} is a pess (heur) sol **/ + /** Update G1\hat{y} **/ for(i = 0; i < lCols; i++){ idx = lColIndices[i]; lwSol[i] = lowerSol[idx]; @@ -1050,7 +1048,7 @@ MibSBilevel::setUpPesUBModel(double *lowerSol, bool newOsi) if(origRowSense[idx] == 'L'){ nSolver->setRowUpper(i, origRowUb[idx] - multG1y[i]); }else if(origRowSense[idx] == 'G'){ - nSolver->setRowLower(i, origRowLb[idx] - multG1y[i]); + nSolver->setRowUpper(i, -1*origRowLb[idx] - multG1y[i]); } } @@ -1245,10 +1243,8 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, } //############################################################################# -// YX: CLEAN UP LATER OsiSolverInterface * -MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, - bool newOsi, const double *lpSol) +MibSBilevel::setUpPesModel(double objValLL, bool newOsi, const double *lpSol) { /** Create pessimistic risk function model with fixed upper-level vars **/ @@ -1259,7 +1255,7 @@ MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, OsiSolverInterface * nSolver; int i(0), j(0), idx(0); - double value(0.0), objUb(0.0); + double objUb(0.0); int uCols(model_->getUpperDim()); int lCols(model_->getLowerDim()); @@ -1267,56 +1263,52 @@ MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, int lRows(model_->getLowerRowNum()); int rowNum(model_->getNumOrigCons() + 1); int colNum(model_->getNumOrigVars()); + int *uRowIndices(model_->getOrigUpperRowInd()); + int *lRowIndices(model_->getLowerRowInd()); int *uColIndices(model_->getUpperColInd()); int *fixedInd(model_->getFixedInd()); - double gap = (targetGap < model_->getTolerance()) ? 0.0 : targetGap; // YX: added SL gap + const double *origRowLb(model_->getOrigRowLb()); + const double *origRowUb(model_->getOrigRowUb()); + char *origRowSense(model_->getOrigRowSense()); + double gap = (targetGap < model_->getTolerance()) ? 0.0 : targetGap; + CoinPackedMatrix *matrixA1(model_->getA1Matrix()); + CoinPackedMatrix *matrixA2(model_->getA2Matrix()); + + double *rowUb = new double[rowNum]; + double *rowLb = new double[rowNum]; + double *uLpSol = new double[uCols]; + double *multA1Xlp = new double[uRows]; + double *multA2Xlp = new double[lRows]; // YX: d2^y <= phi(A^x) + \delta * abs(phi(A^x)) objUb = objValLL + fabs(objValLL) * gap/100; if (!lpSol){ - lpSol = oSolver->getColSolution(); + lpSol = model_->getSolver()->getColSolution(); } if(newOsi){ - double objSense(model_->getLowerObjSense()); + double lObjSense(model_->getLowerObjSense()); double *lObjCoeffs(model_->getLowerObjCoeffs()); int *lColIndices(model_->getLowerColInd()); - int *uRowIndices(model_->getOrigUpperRowInd()); - int *lRowIndices(model_->getLowerRowInd()); const double *origColLb(model_->getOrigColLb()); const double *origColUb(model_->getOrigColUb()); - const double *origRowLb(model_->getOrigRowLb()); - const double *origRowUb(model_->getOrigRowUb()); - const double *uObjCoeffs(oSolver->getObjCoefficients()); - char *origRowSense(model_->getOrigRowSense()); + const double *uObjCoeffs(model_->getSolver()->getObjCoefficients()); - // int tmpRowNum(rowNum -1); int intCnt(0); double uObjSense(1); - // CoinPackedMatrix matrix = *model_->getOrigConstCoefMatrix(); - // matrix.reverseOrdering(); - CoinPackedMatrix *matrixA1 = model_->getA1Matrix(); - CoinPackedMatrix *matrixG1 = model_->getG1Matrix(); - CoinPackedMatrix *matrixA2 = model_->getA2Matrix(); - CoinPackedMatrix *matrixG2 = model_->getG2Matrix(); - // CoinShallowPackedVector origRow; + CoinPackedMatrix *matrixG1(model_->getG1Matrix()); + CoinPackedMatrix *matrixG2(model_->getG2Matrix()); CoinPackedVector row; - // CoinPackedVector addedRow; CoinPackedMatrix *newMat = new CoinPackedMatrix(false, 0, 0); newMat->setDimensions(0, lCols); - double *rowUb = new double[rowNum]; - double *rowLb = new double[rowNum]; double *colUb = new double[lCols]; double *colLb = new double[lCols]; int *integerVars = new int[lCols]; double *objCoeffs = new double[lCols]; double *newRow = new double[lCols]; - double *uLpSol = new double[uCols]; - double *multA1Xlp = new double[uRows]; - double *multA2Xlp = new double[lRows]; if(feasCheckSolver == "Cbc"){ nSolver = new OsiCbcSolverInterface(); @@ -1346,14 +1338,15 @@ MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, CoinZeroN(integerVars, lCols); CoinZeroN(objCoeffs, lCols); CoinZeroN(newRow, lCols); - CoinZeroN(uLpSol, uCols); - /** Obtain A2x matrix **/ + /** Obtain A1x and A2x matrices **/ for(i = 0; i < uCols; i++){ idx = uColIndices[i]; uLpSol[i] = lpSol[idx]; } - matrixA1->times(uLpSol, multA1Xlp); + if(uRows > 0){ + matrixA1->times(uLpSol, multA1Xlp); + } matrixA2->times(uLpSol, multA2Xlp); /** Set row bounds **/ @@ -1363,8 +1356,8 @@ MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, rowLb[i] = origRowLb[idx]; rowUb[i] = origRowUb[idx] - multA1Xlp[i]; }else if(origRowSense[idx] == 'G'){ - rowLb[i] = origRowLb[idx] - multA1Xlp[i]; - rowUb[i] = origRowUb[idx]; + rowLb[i] = -1 * origRowUb[idx]; + rowUb[i] = -1 * origRowLb[idx] - multA1Xlp[i]; } } for(i = 0; i < lRows; i++){ @@ -1373,8 +1366,8 @@ MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, rowLb[i+uRows] = origRowLb[idx]; rowUb[i+uRows] = origRowUb[idx] - multA2Xlp[i]; }else if(origRowSense[idx] == 'G'){ - rowLb[i] = origRowLb[idx] - multA2Xlp[i]; - rowUb[i] = origRowUb[idx]; + rowLb[i+uRows] = -1 * origRowUb[idx]; + rowUb[i+uRows] = -1 * origRowLb[idx] - multA2Xlp[i]; } } @@ -1394,30 +1387,24 @@ MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, row = matrixG2->getVector(i-uRows); newMat->appendRow(row); } + row.clear(); } /** Set up the value function constraint **/ CoinDisjointCopyN(lObjCoeffs, lCols, newRow); for(i = 0; i < lCols; i++){ - newRow[i] = newRow[i] * objSense; + row.insert(i, newRow[i] * lObjSense); } + newMat->appendRow(row); rowUb[rowNum-1] = objUb; - rowLb[rowNum-1] = -1 * (oSolver->getInfinity()); - - for(i = 0; i < lCols; i++){ - idx = lColIndices[i]; - row.insert(idx, newRow[i]); - // row.insert(i, newRow[i]); - } - newMat->appendRow(row); + rowLb[rowNum-1] = -1 * (model_->getSolver()->getInfinity()); - /** Fill in array of integer variables **/ for(i = 0; i < lCols; i++){ idx = lColIndices[i]; - if(oSolver->isInteger(idx)){ + if(model_->getSolver()->isInteger(idx)){ integerVars[intCnt] = i; intCnt ++; } @@ -1431,37 +1418,60 @@ MibSBilevel::setUpPesModel(OsiSolverInterface * oSolver, double objValLL, for(i = 0; i < intCnt; i++){ nSolver->setInteger(integerVars[i]); } - - nSolver->setObjSense(-1 * uObjSense); //1 min; -1 max // YX: pessimistic case + + /** Pessimistic risk function **/ + nSolver->setObjSense(-1 * uObjSense); //1 min; -1 max // YX: still need this? nSolver->setHintParam(OsiDoReducePrint, true, OsiHintDo); - delete [] rowUb; - delete [] rowLb; delete [] colUb; delete [] colLb; delete [] objCoeffs; delete [] integerVars; delete [] newRow; - delete [] uLpSol; - delete [] multA1Xlp; - delete [] multA2Xlp; delete newMat; }else{ - nSolver = pSolver_; // YX: pessimistic case + + nSolver = pSolver_; + + /** Update value function constraint bound **/ nSolver->setRowUpper(rowNum-1, objUb); // YX: gap added to phi(A^2x) + + /** Update row bounds b-A1x and b-A2x **/ for(i = 0; i < uCols; i++){ idx = uColIndices[i]; - if(fixedInd[idx] == 1){ - value = floor(lpSol[idx] + 0.5); - nSolver->setColLower(idx, value); - nSolver->setColUpper(idx, value); + uLpSol[i] = lpSol[idx]; + } + if(uRows > 0){ + matrixA1->times(uLpSol, multA1Xlp); + } + matrixA2->times(uLpSol, multA2Xlp); + + for(i = 0; i < uRows; i++){ + idx = uRowIndices[i]; + if(origRowSense[idx] == 'L'){ + nSolver->setRowUpper(i, origRowUb[idx] - multA1Xlp[i]); + }else if(origRowSense[idx] == 'G'){ + nSolver->setRowUpper(i, -1*origRowLb[idx] - multA1Xlp[i]); + } + } + for(i = 0; i < lRows; i++){ + idx = lRowIndices[i]; + if(origRowSense[idx] == 'L'){ + nSolver->setRowUpper(i+uRows, origRowUb[idx] - multA2Xlp[i]); + }else if(origRowSense[idx] == 'G'){ + nSolver->setRowUpper(i+uRows, -1*origRowLb[idx] - multA2Xlp[i]); } } } + delete [] rowUb; + delete [] rowLb; + delete [] uLpSol; + delete [] multA1Xlp; + delete [] multA2Xlp; return nSolver; } @@ -1947,7 +1957,7 @@ MibSBilevel::getLowerObj(const double * sol, double objSense) //############################################################################# double -MibSBilevel::getRiskFuncVal(double *lowerSol, bool pesType) +MibSBilevel::getRiskFuncVal(double *lowerSol) { int lN(model_->getLowerDim()); int *lColIndices(model_->getLowerColInd()); @@ -1956,24 +1966,23 @@ MibSBilevel::getRiskFuncVal(double *lowerSol, bool pesType) int i(0), idx(0); double rfVal(0.0); - double rfSense = pesType? -1.0 : +1.0; for(i = 0; i < lN; i++){ idx = lColIndices[i]; - if(1){ + if(0){ std::cout << "lowerSol[" << i << "]: " << lowerSol[i] << std::endl; std::cout << "uObjCoeff[" << idx << "]: " << uObjCoeffs[idx] << std::endl; } rfVal += uObjCoeffs[idx] * lowerSol[i]; } - return rfVal * uObjSense * rfSense; + return rfVal * uObjSense; } //############################################################################# double -MibSBilevel::getUpperObj(double *lowerSol, bool pesType) +MibSBilevel::getUpperObj(double *lowerSol) { /** calculate upper-level objective value for (x, hat{y}) or (x, y') **/ int uN(model_->getUpperDim()); @@ -1987,14 +1996,14 @@ MibSBilevel::getUpperObj(double *lowerSol, bool pesType) for(i = 0; i < uN; i++){ idx = uColIndices[i]; - if(1){ + if(0){ std::cout << "lrSol[" << i << "]: " << lrSol[i] << std::endl; std::cout << "uObjCoeff[" << idx << "]: " << uObjCoeffs[idx] << std::endl; } uObjVal += uObjCoeffs[idx] * lrSol[i] * uObjSense; } - rfuncVal = getRiskFuncVal(lowerSol, pesType); + rfuncVal = getRiskFuncVal(lowerSol); uObjVal += rfuncVal; return uObjVal; diff --git a/src/MibSBilevel.hpp b/src/MibSBilevel.hpp index 893b98e2..f9b10436 100644 --- a/src/MibSBilevel.hpp +++ b/src/MibSBilevel.hpp @@ -114,14 +114,14 @@ class MibSBilevel { OsiSolverInterface * setUpPesUBModel(double *lowerSol, bool newOsi); // YX: pessimistic case OsiSolverInterface * setUpUBModel(OsiSolverInterface * solver, double objValLL, bool newOsi, const double *sol = NULL); - OsiSolverInterface * setUpPesModel(OsiSolverInterface * solver, double objValLL, - bool newOsi, const double *sol = NULL); // YX: pessimistic case + OsiSolverInterface * setUpPesModel(double objValLL, bool newOsi, + const double *sol = NULL); // YX: pessimistic case OsiSolverInterface * setUpModel(OsiSolverInterface * solver, bool newOsi, const double *sol = NULL); int findIndex(int index, int size, int * indices); // YX: not used? double getLowerObj(const double * sol, double objSense); - double getRiskFuncVal(double *lowerSol, bool pesType); // YX: pessimistic case - double getUpperObj(double *lowerSol, bool pesType); // YX: pessimistic case + double getRiskFuncVal(double *lowerSol); // YX: pessimistic case + double getUpperObj(double *lowerSol); // YX: pessimistic case int binarySearch(int index,int start, int stop, int * indexArray); CoinWarmStart * getWarmStart() {return ws_;} void setWarmStart(CoinWarmStart * ws) {ws_ = ws;} From e98920c8cb8a7ca043749176a05b62ebf77da125 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Thu, 23 Jun 2022 19:25:04 -0400 Subject: [PATCH 14/26] rm getLowerMatrices() and related lines; replaced by setCoeffMatrices() in MibSModel --- src/MibSCutGenerator.cpp | 48 +++++----------------------------------- 1 file changed, 5 insertions(+), 43 deletions(-) diff --git a/src/MibSCutGenerator.cpp b/src/MibSCutGenerator.cpp index aa86c41f..a836ae80 100644 --- a/src/MibSCutGenerator.cpp +++ b/src/MibSCutGenerator.cpp @@ -686,7 +686,7 @@ MibSCutGenerator::intersectionCuts(BcpsConstraintPool &conPool, mult = -1; break; case 'E': - std::couta + std::cout << "MibS cannot currently handle equality constraints." << std::endl; abort(); break; @@ -828,7 +828,6 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, double timeLimit(localModel_->AlpsPar()->entry(AlpsParams::timeLimit)); double remainingTime(0.0); - // bool getA2Matrix(false), getG2Matrix(false); OsiSolverInterface * oSolver = localModel_->solver(); double infinity(oSolver->getInfinity()); int i, j; @@ -891,18 +890,6 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, int *integerVars = new int[newNumCols]; CoinZeroN(integerVars, newNumCols); - // if(!localModel_->getA2Matrix()){ - // getA2Matrix = true; - // } - - // if(!localModel_->getG2Matrix()){ - // getG2Matrix = true; - // } - - // if((getA2Matrix) || (getG2Matrix)){ - // getLowerMatrices(false, getA2Matrix, getG2Matrix); - // } - if(!localModel_->getA2Matrix()){ localModel_->setCoeffMatrices(); } @@ -1156,24 +1143,11 @@ MibSCutGenerator::getAlphaIC(double** extRay, double* uselessIneqs, char *rowSense = localModel_->getOrigRowSense(); double *lObjCoeffs(localModel_->getLowerObjCoeffs()); double objSense(localModel_->getLowerObjSense()); - // bool getA2Matrix(false), getG2Matrix(false); double targetGap(localModel_->MibSPar_->entry(MibSParams::slTargetGap)); double gap = (targetGap < etol) ? 0.0 : targetGap; double templObj(0.0); // YX: track SL optimal obj val - // if(localModel_->getA2Matrix() == NULL){ - // getA2Matrix = true; - // } - - // if(localModel_->getG2Matrix() == NULL){ - // getG2Matrix = true; - // } - - // if((getA2Matrix) || (getG2Matrix)){ - // getLowerMatrices(false, getA2Matrix, getG2Matrix); - // } - if(!localModel_->getA2Matrix()){ localModel_->setCoeffMatrices(); } @@ -1325,9 +1299,6 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo int whichCutsLL(localModel_->MibSPar_->entry (MibSParams::whichCutsLL)); double targetGap(localModel_->MibSPar_->entry(MibSParams::slTargetGap)); - double etol(localModel_->etol_); - double gap = (targetGap < etol) ? 0.0 : targetGap; // YX: added SL gap - double templObj(0.0); // YX: for nonzero gap, track d^2y^* double timeLimit(localModel_->AlpsPar()->entry(AlpsParams::timeLimit)); double remainingTime(0.0); @@ -1335,10 +1306,12 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo OsiSolverInterface *oSolver = localModel_->solver(); double infinity(oSolver->getInfinity()); - // bool getA2G2Matrix(false), getG2Matrix(false); + double etol(localModel_->getTolerance()); + int i, j; int rowIndex(0), colIndex(0), numElements(0), pos(0), cntInt(0); double rhs(0.0), value(0.0); + double templObj(0.0); // YX: for nonzero gap, track d^2y^* int numCols(localModel_->getNumCols()); int lCols(localModel_->getLowerDim()); int lRows(localModel_->getLowerRowNum()); @@ -1354,6 +1327,7 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo double *origColLb(localModel_->getOrigColLb()); double *origColUb(localModel_->getOrigColUb()); char *origRowSense(localModel_->getOrigRowSense()); + double gap = (targetGap < etol) ? 0.0 : targetGap; // YX: added SL gap CoinPackedMatrix origMatrix = *localModel_->getOrigConstCoefMatrix(); CoinShallowPackedVector origRow; @@ -1365,18 +1339,6 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo origMatrix.reverseOrdering(); - // if(!localModel_->getLowerConstCoefMatrix()){ - // getA2G2Matrix = true; - // } - - // if(!localModel_->getG2Matrix()){ - // getG2Matrix = true; - // } - - // if((getA2G2Matrix) || (getG2Matrix)){ - // getLowerMatrices(getA2G2Matrix, false, getG2Matrix); - // } - if(!localModel_->getLowerConstCoefMatrix()){ localModel_->setCoeffMatrices(); } From 57197e14b762bb22e6091566846e0ce52d2eb560 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Fri, 24 Jun 2022 01:18:17 -0400 Subject: [PATCH 15/26] a index bug in getAlphaWatermelonIC()? --- src/MibSCutGenerator.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MibSCutGenerator.cpp b/src/MibSCutGenerator.cpp index a836ae80..99fe1edf 100644 --- a/src/MibSCutGenerator.cpp +++ b/src/MibSCutGenerator.cpp @@ -1635,11 +1635,11 @@ MibSCutGenerator::getAlphaWatermelonIC(double** extRay, double *uselessIneqs, for(i = 0; i < lRows; i++){ rowIndex = lRowInd[i]; - if(origRowSense[i] == 'L'){ + if(origRowSense[rowIndex] == 'L'){ rhs[i] = origRowUb[rowIndex] + 1 - G2TimeslSol[i] - lCoeffsTimesLpSol[i]; } - else if(origRowSense[i] == 'G'){ + else if(origRowSense[rowIndex] == 'G'){ rhs[i] = -1 * origRowLb[rowIndex] + 1 - G2TimeslSol[i] - lCoeffsTimesLpSol[i]; } From 5f69bf06e9238f8d6c96e431ce108bc57f6c9197 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Tue, 12 Jul 2022 16:34:43 -0400 Subject: [PATCH 16/26] minor format changes --- src/MibSBilevel.cpp | 8 ++++---- src/MibSConstants.hpp | 2 +- src/MibSModel.cpp | 14 +++++++------- src/MibSParams.hpp | 3 +-- 4 files changed, 13 insertions(+), 14 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index ffe211ee..4ae16ffa 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -1970,8 +1970,8 @@ MibSBilevel::getRiskFuncVal(double *lowerSol) for(i = 0; i < lN; i++){ idx = lColIndices[i]; if(0){ - std::cout << "lowerSol[" << i << "]: " << lowerSol[i] << std::endl; - std::cout << "uObjCoeff[" << idx << "]: " << uObjCoeffs[idx] << std::endl; + std::cout << "lowerSol[" << i << "]: " << lowerSol[i] << std::endl; + std::cout << "uObjCoeff[" << idx << "]: " << uObjCoeffs[idx] << std::endl; } rfVal += uObjCoeffs[idx] * lowerSol[i]; } @@ -1997,8 +1997,8 @@ MibSBilevel::getUpperObj(double *lowerSol) for(i = 0; i < uN; i++){ idx = uColIndices[i]; if(0){ - std::cout << "lrSol[" << i << "]: " << lrSol[i] << std::endl; - std::cout << "uObjCoeff[" << idx << "]: " << uObjCoeffs[idx] << std::endl; + std::cout << "lrSol[" << i << "]: " << lrSol[i] << std::endl; + std::cout << "uObjCoeff[" << idx << "]: " << uObjCoeffs[idx] << std::endl; } uObjVal += uObjCoeffs[idx] * lrSol[i] * uObjSense; } diff --git a/src/MibSConstants.hpp b/src/MibSConstants.hpp index a9a494d9..a77295f4 100644 --- a/src/MibSConstants.hpp +++ b/src/MibSConstants.hpp @@ -55,7 +55,7 @@ enum MibSLPSolStatus{ //############################################################################# enum MibSLinkingPoolTag{ - MibSLinkingPoolTagIsNotSet = -4, + MibSLinkingPoolTagIsNotSet = -6, MibSLinkingPoolTagLowerIsInfeasible, MibSLinkingPoolTagLowerIsFeasible, MibSLinkingPoolTagPesIsInfeasible, // YX: pessimistic case diff --git a/src/MibSModel.cpp b/src/MibSModel.cpp index e58d53ca..4bed5662 100644 --- a/src/MibSModel.cpp +++ b/src/MibSModel.cpp @@ -1054,10 +1054,10 @@ MibSModel::loadProblemData(const CoinPackedMatrix& matrix, setRequiredFixedList(); instanceStructure(); - if(MibSPar_->entry(MibSParams::findPesSol)) - { + if(MibSPar_->entry(MibSParams::findPesSol)){ setCoeffMatrices(); } + } //############################################################################# @@ -3121,13 +3121,13 @@ MibSModel::setCoeffMatrices() matrixG1 = new CoinPackedMatrix(false, 0, 0); matrixG1->setDimensions(0, lowerDim_); } - matrixA2G2 = new CoinPackedMatrix(false, 0, 0); - matrixA2G2->setDimensions(0, numVars_); + matrixA2G2 = new CoinPackedMatrix(false, 0, 0); + matrixA2G2->setDimensions(0, numVars_); assert(numVars_== upperDim_ + lowerDim_); // YX: debug only - matrixA2 = new CoinPackedMatrix(false, 0, 0); + matrixA2 = new CoinPackedMatrix(false, 0, 0); matrixA2->setDimensions(0, upperDim_); - matrixG2 = new CoinPackedMatrix(false, 0, 0); - matrixG2->setDimensions(0, lowerDim_); + matrixG2 = new CoinPackedMatrix(false, 0, 0); + matrixG2->setDimensions(0, lowerDim_); // YX: set up A1 and G1 // YX: need to check if it works for interdiction case with AUX rows diff --git a/src/MibSParams.hpp b/src/MibSParams.hpp index 435c83d6..0cfd90c3 100644 --- a/src/MibSParams.hpp +++ b/src/MibSParams.hpp @@ -39,8 +39,7 @@ class MibSParams : public AlpsParameterSet { printProblemInfo, allowRemoveCut, useNewPureIntCut, - // YX: pessimistic case - findPesSol, + findPesSol, // YX: pessimistic case endOfBoolParams }; From 821cd5750caa38ba93002f30317f3b3c39c12fb6 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Wed, 14 Sep 2022 23:16:30 -0400 Subject: [PATCH 17/26] update pess MILP/UB model setup and procedure --- src/MibSBilevel.cpp | 548 +++++++++++++------------------------------- src/MibSBilevel.hpp | 6 +- 2 files changed, 166 insertions(+), 388 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index 4ae16ffa..8bf73c3a 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -224,10 +224,10 @@ MibSBilevel::createBilevel(CoinPackedVector* sol, } //steps 5-6 - if((isLinkVarsFixed_) && ((tagInSeenLinkingPool_ == - MibSLinkingPoolTagLowerIsInfeasible) || - (tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsInfeasible) || // YX: pessimistic case - (tagInSeenLinkingPool_ == MibSLinkingPoolTagUBIsSolved))){ + if((isLinkVarsFixed_) && + ((tagInSeenLinkingPool_ == MibSLinkingPoolTagLowerIsInfeasible) || + (tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsInfeasible) || // YX: pessimistic case + (tagInSeenLinkingPool_ == MibSLinkingPoolTagUBIsSolved))){ useBilevelBranching_ = false; shouldPrune_ = true; } @@ -291,7 +291,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) // YX: target optimality gap for bounded rationality double targetGap(model_->MibSPar()->entry(MibSParams::slTargetGap)); double timeLimit(model_->AlpsPar()->entry(AlpsParams::timeLimit)); - double remainingTime(0.0), startTimeVF(0.0), startTimePES(0.0), startTimeUB(0.0); // YX: pessimistic case + double remainingTime(0.0), startTimeVF(0.0), startTimePES(0.0), startTimeUB(0.0); // YX: (PES-MILP) MibSSolType storeSol(MibSNoSol); int lN(model_->lowerDim_); // lower-level dimension int uN(model_->upperDim_); // upper-level dimension @@ -299,18 +299,18 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) int i(0), index(0), length(0), pos(0); int sizeFixedInd(model_->sizeFixedInd_); double etol(model_->etol_), objVal(0.0), lowerObj(0.0); - double pesRF(0.0), lpPesVal(0.0), lowerPesVal(0.0); // YX: pessimistic case + double pesRF(0.0), lpPesVal(0.0); // YX: pes rf values int * fixedInd = model_->fixedInd_; int * lowerColInd = model_->getLowerColInd(); int * upperColInd = model_->getUpperColInd(); - const double * valuesPes; + const double * valuesPes; // YX: (PES-MILP) solution (x,y) ptr double *lowerSol = new double[lN]; - double *pesSol = new double[lN + uN]; // YX: pessimistic case + double *pesSol = new double[lN + uN]; // YX: processed pes solution CoinFillN(lowerSol, lN, 0.0); CoinFillN(pesSol, lN + uN, 0.0); std::vector shouldStoreValuesLowerSol(lN); - std::vector shouldStoreValuesPesSol(lN); // YX: pessimistic case + std::vector shouldStoreValuesPesSol(lN + uN); // YX: pessimistic case std::vector shouldStoreValuesUBSol(lN + uN); const double * sol = model_->solver()->getColSolution(); @@ -500,7 +500,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) node->setIsBoundSet(true); } - // YX: pessimistic case; need to confirm repeated model setup? + // YX: pessimistic case; enter only when isUpperIntegral_? if(findPesSol){ if(pSolver_){ pSolver_ = setUpPesModel(objVal, false); @@ -510,7 +510,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) OsiSolverInterface *pSolver = pSolver_; - pSolver->writeLp("pSolverLoaded"); + pSolver->writeLp("pSolverLoaded"); // YX: debug only // YX: add timer here remainingTime = timeLimit - model_->broker_->subTreeTimer().getTime(); @@ -573,7 +573,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) startTimePES = model_->broker_->subTreeTimer().getTime(); pSolver->branchAndBound(); model_->timerPES_ += model_->broker_->subTreeTimer().getTime() - startTimePES; - // YX: use SL-MILP counter? + // YX: add tracker later? solved times to match SL-MILP counter? if((feasCheckSolver == "SYMPHONY") && (sym_is_time_limit_reached (dynamic_cast @@ -582,7 +582,6 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) goto TERM_CHECKBILEVELFEAS; } - // isPESSolved_ = true; if(!pSolver->isProvenOptimal()){ isProvenOptimal_ = false; if(useLinkingSolutionPool && isLinkVarsIntegral_){ @@ -594,21 +593,36 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) shouldPrune_ = true; } }else{ - pesRF = pSolver->getObjValue(); + pesRF = pSolver->getObjValue(); // YX: obtain max d^1y valuesPes = pSolver->getColSolution(); - for (i = 0; i < lN; i++){ - if ((pSolver->isInteger(i)) && - (((valuesPes[i] - floor(valuesPes[i])) < etol) || - ((ceil(valuesPes[i]) - valuesPes[i]) < etol))){ - pesSol[i] = (double) floor(valuesPes[i] + 0.5); + // YX: sort pes solutions here with position shift; + // Then saved linking pool pes solutions are sorted and rounded + // (！！which is different from UB solutions.) + for(i = 0; i < uN + lN; i++){ + pos = binarySearch(0, uN - 1, i, upperColInd); + if(pos >= 0){ + if((pSolver->isInteger(i)) && + (((valuesPes[i] - floor(valuesPes[i])) < etol) || + ((ceil(valuesPes[i]) - valuesPes[i]) < etol))){ + pesSol[pos] = (double) floor(valuesPes[i] + 0.5); + }else{ + pesSol[pos] = (double) valuesPes[i]; + } }else{ - pesSol[i] = (double) valuesPes[i]; + pos = binarySearch(0, lN - 1, i, lowerColInd); + if((pSolver->isInteger(i)) && + (((valuesPes[i] - floor(valuesPes[i])) < etol) || + ((ceil(valuesPes[i]) - valuesPes[i]) < etol))){ + pesSol[pos+uN] = (double) floor(valuesPes[i] + 0.5); + }else{ + pesSol[pos+uN] = (double) valuesPes[i]; + } } } if(useLinkingSolutionPool && isLinkVarsIntegral_){ - std::copy(pesSol, pesSol + lN, shouldStoreValuesPesSol.begin()); + std::copy(pesSol, pesSol + uN + lN, shouldStoreValuesPesSol.begin()); addSolutionToSeenLinkingSolutionPool (MibSLinkingPoolTagPesIsFeasible, shouldStoreValuesPesSol, pesRF); shouldStoreValuesPesSol.clear(); @@ -619,7 +633,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) /*if (!warmStartLL){ delete solver_; }*/ - } + } // end of (solving SL-MILP and PES-MILP if not solved before) //step 13 if(((!useLinkingSolutionPool || !isLinkVarsIntegral_) && (isProvenOptimal_)) || @@ -634,7 +648,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) std::copy(model_->seenLinkingSolutions[linkSol].lowerSolution.begin(), model_->seenLinkingSolutions[linkSol].lowerSolution.end(), lowerSol); - if(findPesSol){ // YX: pessimistic case + if(findPesSol){ // YX: pessimistic case; retrieve max d^1y pesRF = model_->seenLinkingSolutions[linkSol].pesRFValue; std::copy(model_->seenLinkingSolutions[linkSol].pesSolution.begin(), model_->seenLinkingSolutions[linkSol].pesSolution.end(), pesSol); @@ -649,7 +663,6 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) LPSolStatus_ = MibSLPSolStatusInfeasible; //step 15 - /** Current solution is bilevel feasible **/ // YX: [d2^y - phi(A^x)]/abs(phi(A^x)), combined check when gap is set to 0; if(!findPesSol){ if(((lowerObj - objVal) <= fabs(objVal) * gap/100) && (isIntegral_)){ @@ -665,33 +678,18 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) } } }else{ - // YX: find pessimistic risk function value using (LR) and (SL-MILP) solutions lpPesVal = getRiskFuncVal(lowerSolutionOrd_); - // YX: whether (LR) solution is a pessimistic solution + // YX: check whether (LR) solution is a pessimistic solution if(((lowerObj - objVal) <= fabs(objVal) * gap/100) && (isIntegral_) && (fabs(lpPesVal - pesRF) < etol)){ useBilevelBranching_ = false; shouldPrune_ = true; storeSol = MibSRelaxationSol; }else{ - lowerPesVal = getRiskFuncVal(lowerSol); - // YX: whether (SL-MILP) solution is a pessimistic solution - if(fabs(lowerPesVal - pesRF) < etol){ - // TODO: VERIFY UPPER LEVEL FEASIBILITY HERE? - // if passed, mark as heurSol; otherwise, always use pesSol; - // memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); - memcpy(optLowerSolutionOrd_, pesSol, sizeof(double) * lN); - if(isUpperIntegral_){ - storeSol = MibSHeurSol; - } - }else{ - // copy pessimistic solution to optLowerSolOrd; use in cut generation - // no need to verify upper-level feasibility - memcpy(optLowerSolutionOrd_, pesSol, sizeof(double) * lN); - if(isUpperIntegral_){ - storeSol = MibSHeurSol; - } - } + // YX: copy (PES-MILP) solution to optimal solution container + memcpy(optUpperSolutionOrd_, pesSol, sizeof(double) * uN); + memcpy(optLowerSolutionOrd_, pesSol + uN, sizeof(double) * lN); + storeSol = MibSHeurSol; } } @@ -704,30 +702,22 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) ((computeBestUBWhenLVarsInt == PARAM_ON) && (isLinkVarsIntegral_)) || ((computeBestUBWhenLVarsFixed == PARAM_ON) && (isLinkVarsFixed_)))){ + // YX: assume (PES-MILP) is solved and a solution is found; if((findPesSol)&&(((sizeFixedInd - uN) >= 0)||(uRows == 0))){ // YX: all variables are linking variables; solution fixed isUBSolved_ = true; - - // YX: UBSolution is not used anywhere? use zero as a place holder; - objVal = getUpperObj(optLowerSolutionOrd_); - shouldStoreValuesUBSol.push_back(0); + + // YX: update Upper bound value; UBSolution default is 0; + objVal = getUpperObj(optLowerSolutionOrd_, optUpperSolutionOrd_); + }else{ - if(findPesSol){ - // YX: when existing non-determined upper vars - if(UBSolver_){ - UBSolver_ = setUpPesUBModel(optLowerSolutionOrd_, false); - }else{ - UBSolver_ = setUpPesUBModel(optLowerSolutionOrd_, true); - } + if(UBSolver_){ + UBSolver_ = setUpUBModel(model_->getSolver(), objVal, false); }else{ - // YX: regular optimistic case - if(UBSolver_){ - UBSolver_ = setUpUBModel(model_->getSolver(), objVal, false); - }else{ - UBSolver_ = setUpUBModel(model_->getSolver(), objVal, true); - } + UBSolver_ = setUpUBModel(model_->getSolver(), objVal, true); } - OsiSolverInterface * UBSolver = UBSolver_; + + OsiSolverInterface * UBSolver = UBSolver_; remainingTime = timeLimit - model_->broker_->subTreeTimer().getTime(); @@ -801,7 +791,11 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) isUBSolved_ = true; - if(!findPesSol){ + // YX: pess UB has to have an optimal solution in case of feasible (PES-MILP) + if(findPesSol){ + assert(UBSolver->isProvenOptimal()); + } + if (UBSolver->isProvenOptimal()){ isProvenOptimal_ = true; const double * valuesUB = UBSolver->getColSolution(); @@ -827,32 +821,16 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) } } } - objVal = UBSolver->getObjValue() * model_->solver()->getObjSense(); + if(!findPesSol){ + objVal = UBSolver->getObjValue() * model_->solver()->getObjSense(); + }else{ + objVal = getUpperObj(optLowerSolutionOrd_, optUpperSolutionOrd_); + } storeSol = MibSHeurSol; }else{ isProvenOptimal_ = false; storeSol = MibSNoSol; } - }else{ - // YX: pess UB has to have an optimal solution b/c feasible (PES-MILP) - assert(UBSolver->isProvenOptimal()); - isProvenOptimal_ = true; - // YX: copy optimal upper-level solution to containers - const double * valuesUB = UBSolver->getColSolution(); - // NEED TO RECOMBINE THEM USING INDEX? - // std::copy(valuesUB, valuesUB + uN, shouldStoreValuesUBSol.begin()); - // std::copy(lowerSol, lowerSol + lN, shouldStoreValuesUBSol.end()); - shouldStoreValuesUBSol.push_back(0); // YX: place holder - for (i = 0; i < uN; i++){ - if ((UBSolver->isInteger(i)) && - (((valuesUB[i] - floor(valuesUB[i])) < etol) || - ((ceil(valuesUB[i]) - valuesUB[i]) < etol))){ - optUpperSolutionOrd_[i] = (double) floor(valuesUB[i] + 0.5); - }else{ - optUpperSolutionOrd_[i] = (double) valuesUB[i]; - } - } - } // end of (processing an UB solution) } // end of (solving an UB problem in either optimistic or pessimistic case) @@ -871,8 +849,8 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) shouldPrune_ = true; } } - } - } + } // end of (if not prune) + } // end of (bilevel feasibility check and UB finding) TERM_CHECKBILEVELFEAS: @@ -897,178 +875,6 @@ MibSBilevel::gutsOfDestructor() //delete heuristic_; } - -//############################################################################# -OsiSolverInterface * -MibSBilevel::setUpPesUBModel(double *lowerSol, bool newOsi) -{ - /** Set up upper bnding problem with x_L after pessimistic MILP is solved **/ - std::string feasCheckSolver = - model_->MibSPar()->entry(MibSParams::feasCheckSolver); - - OsiSolverInterface *nSolver; - - int i(0), idx(0); - double value(0.0); - - int uCols(model_->getUpperDim()); - int uRows(model_->getOrigUpperRowNum()); - int lCols(model_->getLowerDim()); - int *uRowIndices(model_->getOrigUpperRowInd()); - int *uColIndices(model_->getUpperColInd()); - int *lColIndices(model_->getLowerColInd()); - int *fixedInd(model_->getFixedInd()); - char *origRowSense(model_->getOrigRowSense()); - - const double *origRowLb(model_->getOrigRowLb()); - const double *origRowUb(model_->getOrigRowUb()); - const double *lpSol(model_->getSolver()->getColSolution()); - - CoinPackedMatrix *matrixG1(model_->getG1Matrix()); - - double *lwSol = new double[lCols]; - double *multG1y = new double[uRows]; - - if(newOsi){ - double uObjSense(1); - const double *origColLb(model_->getOrigColLb()); - const double *origColUb(model_->getOrigColUb()); - const double *uObjCoeffs(model_->getSolver()->getObjCoefficients()); - CoinPackedMatrix *matrixA1(model_->getA1Matrix()); - - double *rowUb = new double[uRows]; - double *rowLb = new double[uRows]; - double *colUb = new double[uCols]; - double *colLb = new double[uCols]; - - double *multG1y = new double[uRows]; - double *objCoeffs = new double[uCols]; - int *integerVars = new int[uCols]; - - CoinZeroN(colUb, uCols); - CoinZeroN(colLb, uCols); - CoinZeroN(integerVars, uCols); - CoinZeroN(objCoeffs, uCols); - - int intCnt(0); - - if (feasCheckSolver == "Cbc"){ - nSolver = new OsiCbcSolverInterface(); - }else if (feasCheckSolver == "SYMPHONY"){ -#ifdef COIN_HAS_SYMPHONY - nSolver = new OsiSymSolverInterface(); -#else - throw CoinError("SYMPHONY chosen as solver, but it has not been enabled", - "setUpUBModel", "MibsBilevel"); -#endif - }else if (feasCheckSolver == "CPLEX"){ -#ifdef COIN_HAS_CPLEX - nSolver = new OsiCpxSolverInterface(); -#else - throw CoinError("CPLEX chosen as solver, but it has not been enabled", - "setUpUBModel", "MibsBilevel"); -#endif - }else{ - throw CoinError("Unknown solver chosen", - "setUpUBModel", "MibsBilevel"); - } - - /** Obtain G1\hat{y}: \hat{y} is a (heuristic) pessimistic sol **/ - for(i = 0; i < lCols; i++){ - idx = lColIndices[i]; - lwSol[i] = lowerSol[idx]; - } - matrixG1->times(lwSol, multG1y); - - /** Set row bounds **/ - for(i = 0; i < uRows; i++){ - idx = uRowIndices[i]; - if(origRowSense[idx] == 'L'){ - rowLb[i] = origRowLb[idx]; - rowUb[i] = origRowUb[idx] - multG1y[i]; - }else if(origRowSense[idx] == 'G'){ - rowLb[i] = -1*origRowUb[idx]; - rowUb[i] = -1*origRowLb[idx] - multG1y[i]; - } - } - - /** Set col bounds **/ - for(i = 0; i < uCols; i++){ - idx = uColIndices[i]; - if(fixedInd[idx] == 1){ - colLb[i] = floor(lpSol[idx] + 0.5); - colUb[i] = colLb[idx]; - }else{ - colLb[i] = origColLb[idx]; - colUb[i] = origColUb[idx]; - } - } - - /** Fill in array of integer variables and set up obj coeffs**/ - for(i = 0; i < uCols; i++){ - if(model_->getSolver()->isInteger(i)){ - integerVars[intCnt] = i; - intCnt ++; - } - // also set up obj: pessimistic risk function - objCoeffs[i] = uObjCoeffs[idx]; - } - - nSolver->loadProblem(*matrixA1, colLb, colUb, - objCoeffs, rowLb, rowUb); - - for(i = 0; i < intCnt; i++){ - nSolver->setInteger(integerVars[i]); - } - - nSolver->setObjSense(uObjSense); //1 min; -1 max - - nSolver->setHintParam(OsiDoReducePrint, true, OsiHintDo); - - delete [] integerVars; - delete [] rowUb; - delete [] rowLb; - delete [] colUb; - delete [] colLb; - delete [] objCoeffs; - - }else{ - nSolver = UBSolver_; - - /** Update G1\hat{y} **/ - for(i = 0; i < lCols; i++){ - idx = lColIndices[i]; - lwSol[i] = lowerSol[idx]; - } - matrixG1->times(lwSol, multG1y); - - /** Update row bounds **/ - for(i = 0; i < uRows; i++){ - idx = uRowIndices[i]; - if(origRowSense[idx] == 'L'){ - nSolver->setRowUpper(i, origRowUb[idx] - multG1y[i]); - }else if(origRowSense[idx] == 'G'){ - nSolver->setRowUpper(i, -1*origRowLb[idx] - multG1y[i]); - } - } - - /** Update column bounds **/ - for(i = 0; i < uCols; i++){ - idx = uColIndices[i]; - if(fixedInd[idx] == 1){ - value = floor(lpSol[idx] + 0.5); - nSolver->setColLower(idx, value); - nSolver->setColUpper(idx, value); - } - } - } - - delete [] lwSol; - delete [] multG1y; - - return nSolver; -} - //############################################################################# OsiSolverInterface * MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, @@ -1078,6 +884,7 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, std::string feasCheckSolver = model_->MibSPar_->entry(MibSParams::feasCheckSolver); double targetGap(model_->MibSPar_->entry(MibSParams::slTargetGap)); + bool findPes(model_->MibSPar_->entry(MibSParams::findPesSol)); // YX: add pes option OsiSolverInterface * nSolver; @@ -1094,7 +901,7 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, int uRows(model_->getOrigUpperRowNum()); int lRows(model_->getLowerRowNum()); int lCols(model_->getLowerDim()); - int rowNum(model_->getNumOrigCons() + 1); + int rowNum = findPes ? model_->getNumOrigCons() + 2: model_->getNumOrigCons() + 1; int colNum(model_->getNumOrigVars()); int * uColIndices(model_->getUpperColInd()); double gap = (targetGap < model_->getTolerance()) ? 0.0 : targetGap; // YX: added SL gap @@ -1114,7 +921,7 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, const double * origRowUb(model_->getOrigRowUb()); double * lObjCoeffs(model_->getLowerObjCoeffs()); const double * uObjCoeffs(oSolver->getObjCoefficients()); - int tmpRowNum(rowNum -1); + int tmpRowNum = findPes ? rowNum - 2: rowNum - 1; CoinPackedMatrix matrix = *model_->origConstCoefMatrix_; matrix.reverseOrdering(); @@ -1182,14 +989,6 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, } CoinDisjointCopyN(uObjCoeffs, colNum, objCoeffs); - CoinDisjointCopyN(lObjCoeffs, lCols, newRow); - - for(i = 0; i < lCols; i++){ - newRow[i] = newRow[i] * objSense; - } - - rowUb[rowNum-1] = objUb; // YX: gap added to phi(A^2x) - rowLb[rowNum-1] = -1 * (oSolver->getInfinity()); CoinPackedMatrix * newMat = new CoinPackedMatrix(false, 0, 0); newMat->setDimensions(0, colNum); @@ -1197,7 +996,31 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, for(i = 0; i < rowNum - 1; i++){ newMat->appendRow(matrix.getVector(i)); } + /** Add pessimistic risk function constraint at (rowNum-2) **/ + if(findPes){ + CoinPackedVector row2; + for(i = 0; i < lCols; i++){ + index1 = lColIndices[i]; + row2.insert(index1, uObjCoeffs[index1] * uObjSense); + value += uObjCoeffs[index1] * optLowerSolutionOrd_[i]; + // YX: set up obj in pessimistic case cx+0y + objCoeffs[index1] = 0.0; // !! verify obj index + } + newMat->appendRow(row2); + rowUb[rowNum-2] = oSolver->getInfinity(); + rowLb[rowNum-2] = value * uObjSense; + } + + /** Add value function constraint at (rowNum-1) **/ + CoinDisjointCopyN(lObjCoeffs, lCols, newRow); + + for(i = 0; i < lCols; i++){ + newRow[i] = newRow[i] * objSense; + } + rowUb[rowNum-1] = objUb; // YX: gap added to phi(A^2x) + rowLb[rowNum-1] = -1 * (oSolver->getInfinity()); + CoinPackedVector row; for(i = 0; i < lCols; i++){ index1 = lColIndices[i]; @@ -1228,6 +1051,9 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, else{ nSolver = UBSolver_; nSolver->setRowUpper(rowNum-1, objUb); // YX: gap added to phi(A^2x) + if(findPes){ + nSolver->setRowLower(rowNum-2, getRiskFuncVal(optLowerSolutionOrd_)); + } for(i = 0; i < uCols; i++){ index1 = uColIndices[i]; if(fixedInd[index1] == 1){ @@ -1246,7 +1072,7 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, OsiSolverInterface * MibSBilevel::setUpPesModel(double objValLL, bool newOsi, const double *lpSol) { - /** Create pessimistic risk function model with fixed upper-level vars **/ + /** Setup pessimistic risk function model (PES-MILP) with fixed linking part **/ std::string feasCheckSolver(model_->MibSPar()->entry (MibSParams::feasCheckSolver)); @@ -1254,31 +1080,18 @@ MibSBilevel::setUpPesModel(double objValLL, bool newOsi, const double *lpSol) OsiSolverInterface * nSolver; - int i(0), j(0), idx(0); - double objUb(0.0); + int i(0), idx(0); + double value(0.0), objUb(0.0); int uCols(model_->getUpperDim()); int lCols(model_->getLowerDim()); - int uRows(model_->getOrigUpperRowNum()); - int lRows(model_->getLowerRowNum()); int rowNum(model_->getNumOrigCons() + 1); int colNum(model_->getNumOrigVars()); - int *uRowIndices(model_->getOrigUpperRowInd()); - int *lRowIndices(model_->getLowerRowInd()); int *uColIndices(model_->getUpperColInd()); int *fixedInd(model_->getFixedInd()); const double *origRowLb(model_->getOrigRowLb()); const double *origRowUb(model_->getOrigRowUb()); - char *origRowSense(model_->getOrigRowSense()); double gap = (targetGap < model_->getTolerance()) ? 0.0 : targetGap; - CoinPackedMatrix *matrixA1(model_->getA1Matrix()); - CoinPackedMatrix *matrixA2(model_->getA2Matrix()); - - double *rowUb = new double[rowNum]; - double *rowLb = new double[rowNum]; - double *uLpSol = new double[uCols]; - double *multA1Xlp = new double[uRows]; - double *multA2Xlp = new double[lRows]; // YX: d2^y <= phi(A^x) + \delta * abs(phi(A^x)) objUb = objValLL + fabs(objValLL) * gap/100; @@ -1298,16 +1111,19 @@ MibSBilevel::setUpPesModel(double objValLL, bool newOsi, const double *lpSol) int intCnt(0); double uObjSense(1); - CoinPackedMatrix *matrixG1(model_->getG1Matrix()); - CoinPackedMatrix *matrixG2(model_->getG2Matrix()); CoinPackedVector row; - CoinPackedMatrix *newMat = new CoinPackedMatrix(false, 0, 0); - newMat->setDimensions(0, lCols); + CoinPackedMatrix matrix = *model_->origConstCoefMatrix_; + matrix.reverseOrdering(); - double *colUb = new double[lCols]; - double *colLb = new double[lCols]; - int *integerVars = new int[lCols]; - double *objCoeffs = new double[lCols]; + CoinPackedMatrix *newMat = new CoinPackedMatrix(false, 0, 0); + newMat->setDimensions(0, colNum); + + double *colUb = new double[colNum]; + double *colLb = new double[colNum]; + double *rowUb = new double[rowNum]; + double *rowLb = new double[rowNum]; + int *integerVars = new int[colNum]; + double *objCoeffs = new double[colNum]; double *newRow = new double[lCols]; if(feasCheckSolver == "Cbc"){ @@ -1333,84 +1149,56 @@ MibSBilevel::setUpPesModel(double objValLL, bool newOsi, const double *lpSol) CoinZeroN(rowUb, rowNum); CoinZeroN(rowLb, rowNum); - CoinZeroN(colUb, lCols); - CoinZeroN(colLb, lCols); - CoinZeroN(integerVars, lCols); - CoinZeroN(objCoeffs, lCols); + CoinZeroN(colUb, colNum); + CoinZeroN(colLb, colNum); + CoinZeroN(integerVars, colNum); + CoinZeroN(objCoeffs, colNum); CoinZeroN(newRow, lCols); - /** Obtain A1x and A2x matrices **/ - for(i = 0; i < uCols; i++){ - idx = uColIndices[i]; - uLpSol[i] = lpSol[idx]; - } - if(uRows > 0){ - matrixA1->times(uLpSol, multA1Xlp); - } - matrixA2->times(uLpSol, multA2Xlp); - /** Set row bounds **/ - for(i = 0; i < uRows; i++){ - idx = uRowIndices[i]; - if(origRowSense[idx] == 'L'){ - rowLb[i] = origRowLb[idx]; - rowUb[i] = origRowUb[idx] - multA1Xlp[i]; - }else if(origRowSense[idx] == 'G'){ - rowLb[i] = -1 * origRowUb[idx]; - rowUb[i] = -1 * origRowLb[idx] - multA1Xlp[i]; - } - } - for(i = 0; i < lRows; i++){ - idx = lRowIndices[i]; - if(origRowSense[idx] == 'L'){ - rowLb[i+uRows] = origRowLb[idx]; - rowUb[i+uRows] = origRowUb[idx] - multA2Xlp[i]; - }else if(origRowSense[idx] == 'G'){ - rowLb[i+uRows] = -1 * origRowUb[idx]; - rowUb[i+uRows] = -1 * origRowLb[idx] - multA2Xlp[i]; - } - } + memcpy(rowLb, origRowLb, sizeof(double) * (rowNum - 1)); + memcpy(rowUb, origRowUb, sizeof(double) * (rowNum - 1)); /** Set col bounds **/ - // YX: use lColLbInLProb/lColUbInLProb? - for(i = 0; i < lCols; i++){ - colLb[i] = origColLb[lColIndices[i]]; - colUb[i] = origColUb[lColIndices[i]]; + memcpy(colLb, origColLb, sizeof(double) * colNum); + memcpy(colUb, origColUb, sizeof(double) * colNum); + + /** Fix linking variables **/ + for(i = 0; i < uCols; i++){ + idx = uColIndices[i]; + if(fixedInd[idx] == 1){ + colLb[idx] = floor(lpSol[idx] + 0.5); + colUb[idx] = colLb[idx]; + } } - /** Set up new coeffient matrix with G1 and G2 **/ - for(i = 0; i < rowNum-1; i++){ - if(i < uRows){ - row = matrixG1->getVector(i); - newMat->appendRow(row); - }else{ - row = matrixG2->getVector(i-uRows); - newMat->appendRow(row); + /** Fill in array of integer variables **/ + for(i = 0; i < colNum; i++){ + if(model_->getSolver()->isInteger(i)){ + integerVars[intCnt] = i; + intCnt ++; } - row.clear(); } - /** Set up the value function constraint **/ + /** Prepare for value function constraint **/ CoinDisjointCopyN(lObjCoeffs, lCols, newRow); - for(i = 0; i < lCols; i++){ - row.insert(i, newRow[i] * lObjSense); - } - newMat->appendRow(row); - rowUb[rowNum-1] = objUb; rowLb[rowNum-1] = -1 * (model_->getSolver()->getInfinity()); - /** Fill in array of integer variables **/ + /** Fill in new objective and constraint coeffs **/ for(i = 0; i < lCols; i++){ idx = lColIndices[i]; - if(model_->getSolver()->isInteger(idx)){ - integerVars[intCnt] = i; - intCnt ++; - } - // also set up obj: pessimistic risk function - objCoeffs[i] = uObjCoeffs[idx]; - } + objCoeffs[idx] = uObjCoeffs[idx]; // pess risk func 0x+d^1y + row.insert(idx, newRow[i] * lObjSense); // val func d^2y + } + + /** Set up constraints **/ + for(i = 0; i < rowNum - 1; i++){ + newMat->appendRow(matrix.getVector(i)); + } + + newMat->appendRow(row); nSolver->loadProblem(*newMat, colLb, colUb, objCoeffs, rowLb, rowUb); @@ -1427,6 +1215,8 @@ MibSBilevel::setUpPesModel(double objValLL, bool newOsi, const double *lpSol) delete [] colUb; delete [] colLb; + delete [] rowUb; + delete [] rowLb; delete [] objCoeffs; delete [] integerVars; delete [] newRow; @@ -1439,39 +1229,17 @@ MibSBilevel::setUpPesModel(double objValLL, bool newOsi, const double *lpSol) /** Update value function constraint bound **/ nSolver->setRowUpper(rowNum-1, objUb); // YX: gap added to phi(A^2x) - /** Update row bounds b-A1x and b-A2x **/ + /** Update fixed linking variable bounds **/ for(i = 0; i < uCols; i++){ idx = uColIndices[i]; - uLpSol[i] = lpSol[idx]; - } - if(uRows > 0){ - matrixA1->times(uLpSol, multA1Xlp); - } - matrixA2->times(uLpSol, multA2Xlp); - - for(i = 0; i < uRows; i++){ - idx = uRowIndices[i]; - if(origRowSense[idx] == 'L'){ - nSolver->setRowUpper(i, origRowUb[idx] - multA1Xlp[i]); - }else if(origRowSense[idx] == 'G'){ - nSolver->setRowUpper(i, -1*origRowLb[idx] - multA1Xlp[i]); - } - } - for(i = 0; i < lRows; i++){ - idx = lRowIndices[i]; - if(origRowSense[idx] == 'L'){ - nSolver->setRowUpper(i+uRows, origRowUb[idx] - multA2Xlp[i]); - }else if(origRowSense[idx] == 'G'){ - nSolver->setRowUpper(i+uRows, -1*origRowLb[idx] - multA2Xlp[i]); + if(fixedInd[idx] == 1){ + value = floor(lpSol[idx] + 0.5); + nSolver->setColLower(idx, value); + nSolver->setColUpper(idx, value); } } } - delete [] rowUb; - delete [] rowLb; - delete [] uLpSol; - delete [] multA1Xlp; - delete [] multA2Xlp; return nSolver; } @@ -1982,25 +1750,33 @@ MibSBilevel::getRiskFuncVal(double *lowerSol) //############################################################################# double -MibSBilevel::getUpperObj(double *lowerSol) +MibSBilevel::getUpperObj(double *lowerSol, double *upperSol) { /** calculate upper-level objective value for (x, hat{y}) or (x, y') **/ + // YX: assuming upperSol passed in are ordered; int uN(model_->getUpperDim()); int *uColIndices(model_->getUpperColInd()); const double uObjSense(model_->solver()->getObjSense()); const double *uObjCoeffs(model_->solver()->getObjCoefficients()); - const double *lrSol(model_->solver()->getColSolution()); + + if(!upperSol){ + const double *upperSol(model_->solver()->getColSolution()); + } int i(0), idx(0); double rfuncVal(0.0), uObjVal(0.0); for(i = 0; i < uN; i++){ - idx = uColIndices[i]; + if(!upperSol){ + idx = uColIndices[i]; + }else{ + idx = i; + } if(0){ - std::cout << "lrSol[" << i << "]: " << lrSol[i] << std::endl; + std::cout << "upperSol[" << i << "]: " << upperSol[i] << std::endl; std::cout << "uObjCoeff[" << idx << "]: " << uObjCoeffs[idx] << std::endl; } - uObjVal += uObjCoeffs[idx] * lrSol[i] * uObjSense; + uObjVal += uObjCoeffs[idx] * upperSol[i] * uObjSense; } rfuncVal = getRiskFuncVal(lowerSol); @@ -2064,7 +1840,7 @@ void model_->seenLinkingSolutions[linkSol] = linkingSolution; break; } - case MibSLinkingPoolTagPesIsInfeasible: // YX: pessimistic case; combined to one using provenoptimal? + case MibSLinkingPoolTagPesIsInfeasible: // YX: pessimistic case; following one tag above; { model_->seenLinkingSolutions[linkSol].tag = MibSLinkingPoolTagPesIsInfeasible; break; diff --git a/src/MibSBilevel.hpp b/src/MibSBilevel.hpp index f9b10436..c2094aa4 100644 --- a/src/MibSBilevel.hpp +++ b/src/MibSBilevel.hpp @@ -58,6 +58,9 @@ class MibSBilevel { bool isContainedInLinkingPool_; MibSLinkingPoolTag tagInSeenLinkingPool_; + // YX: this variable is not used anywhere other than MibSBilevel; + // (commented off in MibSTreeNode) + // looks like it is replaced by pool tag and can be removed; MibSLPSolStatus LPSolStatus_; /** Optimal value of LL objective **/ @@ -111,7 +114,6 @@ class MibSBilevel { private: - OsiSolverInterface * setUpPesUBModel(double *lowerSol, bool newOsi); // YX: pessimistic case OsiSolverInterface * setUpUBModel(OsiSolverInterface * solver, double objValLL, bool newOsi, const double *sol = NULL); OsiSolverInterface * setUpPesModel(double objValLL, bool newOsi, @@ -121,7 +123,7 @@ class MibSBilevel { int findIndex(int index, int size, int * indices); // YX: not used? double getLowerObj(const double * sol, double objSense); double getRiskFuncVal(double *lowerSol); // YX: pessimistic case - double getUpperObj(double *lowerSol); // YX: pessimistic case + double getUpperObj(double *lowerSol, double *upperSol = NULL); // YX: pessimistic case int binarySearch(int index,int start, int stop, int * indexArray); CoinWarmStart * getWarmStart() {return ws_;} void setWarmStart(CoinWarmStart * ws) {ws_ = ws;} From 5c00ef196750344f1f0364a545236337d4054ab8 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Thu, 15 Sep 2022 16:39:34 -0400 Subject: [PATCH 18/26] update pes conditions when NOT upperIntegral --- src/MibSBilevel.cpp | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index 8bf73c3a..75907c23 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -291,6 +291,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) // YX: target optimality gap for bounded rationality double targetGap(model_->MibSPar()->entry(MibSParams::slTargetGap)); double timeLimit(model_->AlpsPar()->entry(AlpsParams::timeLimit)); + bool isPesOptimal = false; // YX: for pes check when no linking pool double remainingTime(0.0), startTimeVF(0.0), startTimePES(0.0), startTimeUB(0.0); // YX: (PES-MILP) MibSSolType storeSol(MibSNoSol); int lN(model_->lowerDim_); // lower-level dimension @@ -501,7 +502,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) } // YX: pessimistic case; enter only when isUpperIntegral_? - if(findPesSol){ + if(findPesSol && isUpperIntegral_){ if(pSolver_){ pSolver_ = setUpPesModel(objVal, false); }else{ @@ -573,7 +574,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) startTimePES = model_->broker_->subTreeTimer().getTime(); pSolver->branchAndBound(); model_->timerPES_ += model_->broker_->subTreeTimer().getTime() - startTimePES; - // YX: add tracker later? solved times to match SL-MILP counter? + model_->counterPES_ ++; if((feasCheckSolver == "SYMPHONY") && (sym_is_time_limit_reached (dynamic_cast @@ -583,7 +584,6 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) } if(!pSolver->isProvenOptimal()){ - isProvenOptimal_ = false; if(useLinkingSolutionPool && isLinkVarsIntegral_){ addSolutionToSeenLinkingSolutionPool (MibSLinkingPoolTagPesIsInfeasible, shouldStoreValuesPesSol, 0.0); @@ -593,6 +593,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) shouldPrune_ = true; } }else{ + isPesOptimal = true; pesRF = pSolver->getObjValue(); // YX: obtain max d^1y valuesPes = pSolver->getColSolution(); @@ -639,6 +640,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) if(((!useLinkingSolutionPool || !isLinkVarsIntegral_) && (isProvenOptimal_)) || ((tagInSeenLinkingPool_ == MibSLinkingPoolTagLowerIsFeasible) || (tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsFeasible) || + (tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsInfeasible) || (tagInSeenLinkingPool_ == MibSLinkingPoolTagUBIsSolved))){ if(useLinkingSolutionPool && isLinkVarsIntegral_){ @@ -647,8 +649,9 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) objVal_ = objVal; std::copy(model_->seenLinkingSolutions[linkSol].lowerSolution.begin(), model_->seenLinkingSolutions[linkSol].lowerSolution.end(), lowerSol); - - if(findPesSol){ // YX: pessimistic case; retrieve max d^1y + + // YX: pessimistic case; retrieve max d^1y only when (PES-MILP) is feasible + if(findPesSol && (tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsFeasible)){ pesRF = model_->seenLinkingSolutions[linkSol].pesRFValue; std::copy(model_->seenLinkingSolutions[linkSol].pesSolution.begin(), model_->seenLinkingSolutions[linkSol].pesSolution.end(), pesSol); @@ -657,7 +660,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) lowerObj = getLowerObj(sol, model_->getLowerObjSense()); if(isIntegral_){ - assert((objVal - lowerObj) <= etol); + assert((objVal - lowerObj) <= etol); } LPSolStatus_ = MibSLPSolStatusInfeasible; @@ -677,7 +680,8 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) storeSol = MibSHeurSol; } } - }else{ + }else if((tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsFeasible) || + (isPesOptimal)){ lpPesVal = getRiskFuncVal(lowerSolutionOrd_); // YX: check whether (LR) solution is a pessimistic solution if(((lowerObj - objVal) <= fabs(objVal) * gap/100) && (isIntegral_) && @@ -691,11 +695,15 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) memcpy(optLowerSolutionOrd_, pesSol + uN, sizeof(double) * lN); storeSol = MibSHeurSol; } - } + }else{ + // YX: if (pes) and (PES-MILP) is unsolved or infeasible; for fractional cut + memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); + } if(!shouldPrune_){ //step 18 - if((tagInSeenLinkingPool_ != MibSLinkingPoolTagUBIsSolved) && + if(((tagInSeenLinkingPool_ != MibSLinkingPoolTagUBIsSolved) && + (tagInSeenLinkingPool_ != MibSLinkingPoolTagPesIsInfeasible)) && (((branchPar == MibSBranchingStrategyLinking) && (isIntegral_) && (isLinkVarsFixed_)) || ((computeBestUBWhenXVarsInt == PARAM_ON) && (isUpperIntegral_)) || @@ -796,7 +804,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) assert(UBSolver->isProvenOptimal()); } - if (UBSolver->isProvenOptimal()){ + if(UBSolver->isProvenOptimal()){ isProvenOptimal_ = true; const double * valuesUB = UBSolver->getColSolution(); std::copy(valuesUB, valuesUB + uN + lN, shouldStoreValuesUBSol.begin()); From 96e3161a936e8212867d699e80b0d76da352ece7 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Thu, 15 Sep 2022 16:53:40 -0400 Subject: [PATCH 19/26] add pessimistic related counter and output --- src/MibSModel.cpp | 1 + src/MibSModel.hpp | 3 +++ src/MibSSolution.cpp | 9 ++++++++- 3 files changed, 12 insertions(+), 1 deletion(-) diff --git a/src/MibSModel.cpp b/src/MibSModel.cpp index 4bed5662..dfb4919f 100644 --- a/src/MibSModel.cpp +++ b/src/MibSModel.cpp @@ -134,6 +134,7 @@ MibSModel::initialize() structRowNum_ = 0; sizeFixedInd_ = 0; counterVF_ = 0; + counterPES_ = 0; counterUB_ = 0; timerVF_ = 0.0; timerPES_ = 0.0; // YX: pessmistic case diff --git a/src/MibSModel.hpp b/src/MibSModel.hpp index 7fbcc89b..b777affa 100644 --- a/src/MibSModel.hpp +++ b/src/MibSModel.hpp @@ -99,6 +99,9 @@ class MibSModel : public BlisModel { /** Number of (VF) solved **/ int counterVF_; + /** Number of (PES) solved **/ + int counterPES_; // YX: counter for (PES-MILP); + /** Number of (UB) solved **/ int counterUB_; diff --git a/src/MibSSolution.cpp b/src/MibSSolution.cpp index c33c8743..f790cc16 100644 --- a/src/MibSSolution.cpp +++ b/src/MibSSolution.cpp @@ -82,7 +82,8 @@ void MibSSolution::print(std::ostream& os) const { - std::string inputFormat(localModel_->MibSPar_->entry + bool findPesSol(localModel_->MibSPar_->entry(MibSParams::findPesSol)); + std::string inputFormat(localModel_->MibSPar_->entry (MibSParams::inputFormat)); double nearInt = 0.0; //int size(localModel_->getNumOrigVars()); @@ -157,8 +158,14 @@ MibSSolution::print(std::ostream& os) const } std::cout << "Number of problems (VF) solved = " << localModel_->counterVF_ << std::endl; + if(findPesSol){ + std::cout << "Number of problems (PES) solved = " << localModel_->counterPES_ << std::endl; + } std::cout << "Number of problems (UB) solved = " << localModel_->counterUB_ << std::endl; std::cout << "Time for solving problem (VF) = " << localModel_->timerVF_ << std::endl; + if(findPesSol){ + std::cout << "Time for solving problem (PES) = " << localModel_->timerPES_ << std::endl; + } std::cout << "Time for solving problem (UB) = " << localModel_->timerUB_ << std::endl; } From ff75bed95802ad9ae2504b734bd87897e3fb1ee3 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Wed, 21 Sep 2022 01:23:04 -0400 Subject: [PATCH 20/26] fix pes bugs and update UB conds in cut generator --- src/MibSBilevel.cpp | 63 +++++++++++++++++++++++++--------------- src/MibSCutGenerator.cpp | 12 ++++---- src/MibSModel.cpp | 9 ++++-- src/MibSModel.hpp | 2 +- 4 files changed, 51 insertions(+), 35 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index 75907c23..fb8a93c3 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -694,6 +694,27 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) memcpy(optUpperSolutionOrd_, pesSol, sizeof(double) * uN); memcpy(optLowerSolutionOrd_, pesSol + uN, sizeof(double) * lN); storeSol = MibSHeurSol; + + // YX: no need to solve (UB) when all upper vars are linking + // or no upper-level constraints + if(((sizeFixedInd - uN) >= 0)||(uRows == 0)){ + isUBSolved_ = true; + + // YX: update upper bound value; UBSolution default is 0; + objVal = getUpperObj(optLowerSolutionOrd_, optUpperSolutionOrd_); + + if(useLinkingSolutionPool && isLinkVarsIntegral_){ + std::copy(pesSol, pesSol + uN + lN, shouldStoreValuesUBSol.begin()); + addSolutionToSeenLinkingSolutionPool + (MibSLinkingPoolTagUBIsSolved, shouldStoreValuesUBSol, objVal); + shouldStoreValuesUBSol.clear(); + } + + if(isLinkVarsFixed_){ + useBilevelBranching_ = false; + shouldPrune_ = true; + } + } } }else{ // YX: if (pes) and (PES-MILP) is unsolved or infeasible; for fractional cut @@ -710,23 +731,18 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) ((computeBestUBWhenLVarsInt == PARAM_ON) && (isLinkVarsIntegral_)) || ((computeBestUBWhenLVarsFixed == PARAM_ON) && (isLinkVarsFixed_)))){ - // YX: assume (PES-MILP) is solved and a solution is found; - if((findPesSol)&&(((sizeFixedInd - uN) >= 0)||(uRows == 0))){ - // YX: all variables are linking variables; solution fixed - isUBSolved_ = true; - - // YX: update Upper bound value; UBSolution default is 0; - objVal = getUpperObj(optLowerSolutionOrd_, optUpperSolutionOrd_); - + // YX: if (pes), then (PES-MILP) is solved and a solution is found; + + if(UBSolver_){ + UBSolver_ = setUpUBModel(model_->getSolver(), objVal, false); }else{ - if(UBSolver_){ - UBSolver_ = setUpUBModel(model_->getSolver(), objVal, false); - }else{ - UBSolver_ = setUpUBModel(model_->getSolver(), objVal, true); - } + UBSolver_ = setUpUBModel(model_->getSolver(), objVal, true); + } OsiSolverInterface * UBSolver = UBSolver_; + UBSolver->writeLp("UBSolverLoaded"); // YX: debug only + remainingTime = timeLimit - model_->broker_->subTreeTimer().getTime(); if(remainingTime <= etol){ @@ -800,7 +816,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) isUBSolved_ = true; // YX: pess UB has to have an optimal solution in case of feasible (PES-MILP) - if(findPesSol){ + if(findPesSol){ //YX: for debug assert(UBSolver->isProvenOptimal()); } @@ -839,9 +855,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) isProvenOptimal_ = false; storeSol = MibSNoSol; } - - } // end of (solving an UB problem in either optimistic or pessimistic case) - + //step 22 //Adding x_L to set E if(useLinkingSolutionPool && isLinkVarsIntegral_){ @@ -856,7 +870,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) //isProvenOptimal_ = false; shouldPrune_ = true; } - } + } // end of (solving optimistic or pessimistic UB and saving solution) } // end of (if not prune) } // end of (bilevel feasibility check and UB finding) @@ -929,7 +943,7 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, const double * origRowUb(model_->getOrigRowUb()); double * lObjCoeffs(model_->getLowerObjCoeffs()); const double * uObjCoeffs(oSolver->getObjCoefficients()); - int tmpRowNum = findPes ? rowNum - 2: rowNum - 1; + int origRowNum = findPes ? rowNum - 2: rowNum - 1; CoinPackedMatrix matrix = *model_->origConstCoefMatrix_; matrix.reverseOrdering(); @@ -943,8 +957,8 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, CoinZeroN(colLb, colNum); /** Set the row bounds **/ - memcpy(rowLb, origRowLb, sizeof(double) * tmpRowNum); - memcpy(rowUb, origRowUb, sizeof(double) * tmpRowNum); + memcpy(rowLb, origRowLb, sizeof(double) * origRowNum); + memcpy(rowUb, origRowUb, sizeof(double) * origRowNum); //Set the col bounds memcpy(colLb, origColLb, sizeof(double) * colNum); @@ -1001,18 +1015,19 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, CoinPackedMatrix * newMat = new CoinPackedMatrix(false, 0, 0); newMat->setDimensions(0, colNum); - for(i = 0; i < rowNum - 1; i++){ + for(i = 0; i < origRowNum; i++){ newMat->appendRow(matrix.getVector(i)); } /** Add pessimistic risk function constraint at (rowNum-2) **/ + // YX: assume optLowerSolutionOrd_ always contains the pes solution; may add to param later if(findPes){ CoinPackedVector row2; for(i = 0; i < lCols; i++){ index1 = lColIndices[i]; row2.insert(index1, uObjCoeffs[index1] * uObjSense); - value += uObjCoeffs[index1] * optLowerSolutionOrd_[i]; + value += uObjCoeffs[index1] * optLowerSolutionOrd_[i] * uObjSense; // YX: set up obj in pessimistic case cx+0y - objCoeffs[index1] = 0.0; // !! verify obj index + objCoeffs[index1] = 0.0; } newMat->appendRow(row2); rowUb[rowNum-2] = oSolver->getInfinity(); diff --git a/src/MibSCutGenerator.cpp b/src/MibSCutGenerator.cpp index 99fe1edf..b5f15dd6 100644 --- a/src/MibSCutGenerator.cpp +++ b/src/MibSCutGenerator.cpp @@ -558,6 +558,7 @@ MibSCutGenerator::intersectionCuts(BcpsConstraintPool &conPool, (bS->isProvenOptimal_ == false)))) || ((useLinkingSolutionPool == PARAM_ON) && ((bS->tagInSeenLinkingPool_ == MibSLinkingPoolTagLowerIsInfeasible) || + (bS->tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsInfeasible) || (bS->tagInSeenLinkingPool_ == MibSLinkingPoolTagUBIsSolved)))){ shouldFindBestSol = false; } @@ -1757,6 +1758,8 @@ MibSCutGenerator::storeBestSolHypercubeIC(const double* lpSol, double optLowerOb UBSolver = bS->UBSolver_; + UBSolver->writeLp("storeBestUBSolver"); // YX: debug only + remainingTime = timeLimit - localModel_->broker_->subTreeTimer().getTime(); if(remainingTime <= localModel_->etol_){ @@ -1845,18 +1848,13 @@ MibSCutGenerator::storeBestSolHypercubeIC(const double* lpSol, double optLowerOb if((useLinkingSolutionPool) && (isTimeLimReached == false)){ //Add to linking solution pool - //localModel_->it = localModel_->seenLinkingSolutions.find(linkSol); - //localModel_->it->second.tag = MibSSetETagUBIsSolved; - //localModel_->it->second.UBObjVal1 = objVal; localModel_->bS_->tagInSeenLinkingPool_ = MibSLinkingPoolTagUBIsSolved; localModel_->seenLinkingSolutions[linkSol].tag = MibSLinkingPoolTagUBIsSolved; localModel_->seenLinkingSolutions[linkSol].UBObjValue = objVal; if(UBSolver->isProvenOptimal()){ localModel_->seenLinkingSolutions[linkSol].UBSolution.clear(); - //localModel_->it->second.UBSol1.clear(); const double * valuesUB = UBSolver->getColSolution(); for(i = 0; i < uN + lN; i++){ - //localModel_->it->second.UBSol1.push_back(valuesUB[i]); localModel_->seenLinkingSolutions[linkSol].UBSolution.push_back(valuesUB[i]); } } @@ -3254,8 +3252,8 @@ MibSCutGenerator::generalNoGoodCut(BcpsConstraintPool &conPool) (bS->isProvenOptimal_ == false)))) || ((useLinkingSolutionPool == PARAM_ON) && ((bS->tagInSeenLinkingPool_ == MibSLinkingPoolTagLowerIsInfeasible) || - (bS->tagInSeenLinkingPool_ == - MibSLinkingPoolTagUBIsSolved)))){ + (bS->tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsInfeasible) || + (bS->tagInSeenLinkingPool_ == MibSLinkingPoolTagUBIsSolved)))){ shouldFindBestSol = false; } diff --git a/src/MibSModel.cpp b/src/MibSModel.cpp index dfb4919f..cb8bde4e 100644 --- a/src/MibSModel.cpp +++ b/src/MibSModel.cpp @@ -2006,9 +2006,12 @@ MibSModel::userFeasibleSolution(const double * solution, bool &userFeasible) this); } else if(solType == MibSHeurSol){ - if(!bS_->isUBSolved_){ - isHeurSolution = checkUpperFeasibility(lpSolution); - } + if(MibSPar_->entry(MibSParams::findPesSol)){ + // YX: skip upperFeasibility if PES-MILP is feasible + isHeurSolution = true; + }else if(!bS_->isUBSolved_){ + isHeurSolution = checkUpperFeasibility(lpSolution); + } if(isHeurSolution == true){ mibSol = new MibSSolution(getNumCols(), lpSolution, diff --git a/src/MibSModel.hpp b/src/MibSModel.hpp index b777affa..0e3ef986 100644 --- a/src/MibSModel.hpp +++ b/src/MibSModel.hpp @@ -224,7 +224,7 @@ class MibSModel : public BlisModel { /** Original matrix of constraints coefficients **/ CoinPackedMatrix *origConstCoefMatrix_; - // YX: used in pessimistic feasibility check + // YX: was used in pessimistic feasibility check //we only fill next three matrices if they are required. //For now, we only need them for generating IC and watermelon IC. /** matrix of lower-level coefficients(all constraints are 'L') **/ From 91fe7ee4cf00ed76d21200e915521667d7ea74b2 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Wed, 21 Sep 2022 11:31:18 -0400 Subject: [PATCH 21/26] debug writeLp off --- src/MibSBilevel.cpp | 8 ++++---- src/MibSCutGenerator.cpp | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index fb8a93c3..e159a8c1 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -349,8 +349,8 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) OsiSolverInterface *lSolver = lSolver_; - lSolver->writeLp("lSolver"); - model_->getSolver()->writeLp("oSolver"); + // lSolver->writeLp("lSolver"); // YX: debug only + // model_->getSolver()->writeLp("oSolver"); // YX: debug only remainingTime = timeLimit - model_->broker_->subTreeTimer().getTime(); if(remainingTime <= etol){ @@ -511,7 +511,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) OsiSolverInterface *pSolver = pSolver_; - pSolver->writeLp("pSolverLoaded"); // YX: debug only + // pSolver->writeLp("pSolverLoaded"); // YX: debug only // YX: add timer here remainingTime = timeLimit - model_->broker_->subTreeTimer().getTime(); @@ -741,7 +741,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) OsiSolverInterface * UBSolver = UBSolver_; - UBSolver->writeLp("UBSolverLoaded"); // YX: debug only + // UBSolver->writeLp("UBSolverLoaded"); // YX: debug only remainingTime = timeLimit - model_->broker_->subTreeTimer().getTime(); diff --git a/src/MibSCutGenerator.cpp b/src/MibSCutGenerator.cpp index b5f15dd6..072625a4 100644 --- a/src/MibSCutGenerator.cpp +++ b/src/MibSCutGenerator.cpp @@ -1758,7 +1758,7 @@ MibSCutGenerator::storeBestSolHypercubeIC(const double* lpSol, double optLowerOb UBSolver = bS->UBSolver_; - UBSolver->writeLp("storeBestUBSolver"); // YX: debug only + // UBSolver->writeLp("storeBestUBSolver"); // YX: debug only remainingTime = timeLimit - localModel_->broker_->subTreeTimer().getTime(); From 5be09a8665dd11f46e23faca7c1364208ee07e0a Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Tue, 27 Sep 2022 21:45:38 -0400 Subject: [PATCH 22/26] add gap to benders cuts; use vfLowerSolutionOrd_ to separate optimal SL-MILP solution from others --- src/MibSBilevel.cpp | 19 +++++++++------- src/MibSBilevel.hpp | 4 +++- src/MibSCutGenerator.cpp | 47 +++++++++++++++++++++------------------- 3 files changed, 39 insertions(+), 31 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index e159a8c1..8044a116 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -97,11 +97,15 @@ MibSBilevel::createBilevel(CoinPackedVector* sol, if(!optLowerSolutionOrd_){ optLowerSolutionOrd_ = new double[lN]; } + if(!vfLowerSolutionOrd_){ + vfLowerSolutionOrd_ = new double[lN]; // YX: SL-MILP solutions + } CoinZeroN(upperSolutionOrd_, uN); CoinZeroN(lowerSolutionOrd_, lN); CoinZeroN(optUpperSolutionOrd_, uN); CoinZeroN(optLowerSolutionOrd_, lN); + CoinZeroN(vfLowerSolutionOrd_, lN); // YX: SL-MILP solutions isIntegral_ = true; isUpperIntegral_ = true; @@ -475,7 +479,10 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) } } - if(useLinkingSolutionPool && isLinkVarsIntegral_){ + // YX: keep SL-MILP solutions separately for cut generation + memcpy(vfLowerSolutionOrd_, lowerSol, sizeof(double) * lN); + + if(useLinkingSolutionPool && isLinkVarsIntegral_){ std::copy(lowerSol, lowerSol + lN, shouldStoreValuesLowerSol.begin()); //step 12 @@ -640,7 +647,6 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) if(((!useLinkingSolutionPool || !isLinkVarsIntegral_) && (isProvenOptimal_)) || ((tagInSeenLinkingPool_ == MibSLinkingPoolTagLowerIsFeasible) || (tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsFeasible) || - (tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsInfeasible) || (tagInSeenLinkingPool_ == MibSLinkingPoolTagUBIsSolved))){ if(useLinkingSolutionPool && isLinkVarsIntegral_){ @@ -680,8 +686,7 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) storeSol = MibSHeurSol; } } - }else if((tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsFeasible) || - (isPesOptimal)){ + }else{ lpPesVal = getRiskFuncVal(lowerSolutionOrd_); // YX: check whether (LR) solution is a pessimistic solution if(((lowerObj - objVal) <= fabs(objVal) * gap/100) && (isIntegral_) && @@ -716,11 +721,8 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) } } } - }else{ - // YX: if (pes) and (PES-MILP) is unsolved or infeasible; for fractional cut - memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); } - + if(!shouldPrune_){ //step 18 if(((tagInSeenLinkingPool_ != MibSLinkingPoolTagUBIsSolved) && @@ -889,6 +891,7 @@ MibSBilevel::gutsOfDestructor() if(optUpperSolutionOrd_) delete [] optUpperSolutionOrd_; if(optLowerSolutionOrd_) delete [] optLowerSolutionOrd_; + if(vfLowerSolutionOrd_) delete [] vfLowerSolutionOrd_; // YX: SL-MILP solutions if(upperSolutionOrd_) delete [] upperSolutionOrd_; if(lowerSolutionOrd_) delete [] lowerSolutionOrd_; if(lSolver_) delete lSolver_; diff --git a/src/MibSBilevel.hpp b/src/MibSBilevel.hpp index c2094aa4..e35cb3fc 100644 --- a/src/MibSBilevel.hpp +++ b/src/MibSBilevel.hpp @@ -71,7 +71,8 @@ class MibSBilevel { double *lowerSolutionOrd_; //double *optLowerSolution_; double *optUpperSolutionOrd_;// result of solving (UB) - double *optLowerSolutionOrd_; + double *optLowerSolutionOrd_; // YX: optimal lower (heuristic) solutions + double *vfLowerSolutionOrd_; // YX: SL-MILP solutions; needed for gap > 0 MibSModel *model_; MibSHeuristic *heuristic_; @@ -95,6 +96,7 @@ class MibSBilevel { //optLowerSolution_ = 0; optUpperSolutionOrd_ = 0; optLowerSolutionOrd_ = 0; + vfLowerSolutionOrd_ = 0; // YX: SL-MILP solutions; for gap > 0 model_ = 0; heuristic_= 0; lSolver_ = 0; diff --git a/src/MibSCutGenerator.cpp b/src/MibSCutGenerator.cpp index 072625a4..0122cc4d 100644 --- a/src/MibSCutGenerator.cpp +++ b/src/MibSCutGenerator.cpp @@ -3324,7 +3324,7 @@ MibSCutGenerator::findDeepestLandPCut_ValFunc() int numLLCols = localModel_->getLowerDim(); double * upperSol = bS->upperSolutionOrd_; // UL portion from LR soln double * lowerSol = bS->lowerSolutionOrd_; // LL portion from LR soln - double * optLowerSol = bS->optLowerSolutionOrd_; // optimal LL solution + double * optLowerSol = bS->vfLowerSolutionOrd_; // optimal LL solution /* the negative and positive slopes of the LL value function */ double leftSlope = localModel_->leftSlope_; @@ -4286,7 +4286,7 @@ MibSCutGenerator::findDeepestLandPCut1() int numLLCols = localModel_->getLowerDim(); double * upperSol = bS->upperSolutionOrd_; // UL portion from LR soln double * lowerSol = bS->lowerSolutionOrd_; // LL portion from LR soln - double * optLowerSol = bS->optLowerSolutionOrd_; // optimal LL solution + double * optLowerSol = bS->vfLowerSolutionOrd_; // optimal LL solution /* the negative and positive slopes of the LL value function */ double leftSlope = localModel_->leftSlope_; @@ -4841,7 +4841,11 @@ MibSCutGenerator::generalWeakIncObjCutCurrent(BcpsConstraintPool &conPool) { bool allowRemoveCut(localModel_->MibSPar_->entry (MibSParams::allowRemoveCut)); - + int useLinkingSolutionPool(localModel_->MibSPar_->entry + (MibSParams::useLinkingSolutionPool)); + double targetGap(localModel_->MibSPar_->entry(MibSParams::slTargetGap)); + + MibSBilevel *bS = localModel_->bS_; OsiSolverInterface *solver = localModel_->solver(); int i; int index(0); @@ -4849,6 +4853,7 @@ MibSCutGenerator::generalWeakIncObjCutCurrent(BcpsConstraintPool &conPool) double cutub(0.0), cutlb(-solver->getInfinity()), bigM(0.0); double lObjVal(localModel_->bS_->objVal_); double etol(localModel_->etol_); + double gap = (targetGap < etol) ? 0.0 : targetGap; // YX: added SL gap int uN(localModel_->getUpperDim()); int lN(localModel_->getLowerDim()); int *fixedInd(localModel_->getFixedInd()); @@ -4858,11 +4863,6 @@ MibSCutGenerator::generalWeakIncObjCutCurrent(BcpsConstraintPool &conPool) double lowerObjSense(localModel_->getLowerObjSense()); const double *sol(solver->getColSolution()); - int useLinkingSolutionPool(localModel_->MibSPar_->entry - (MibSParams::useLinkingSolutionPool)); - - MibSBilevel *bS = localModel_->bS_; - bool shouldFindBestSol(true); std::vector indexList; @@ -4907,7 +4907,7 @@ MibSCutGenerator::generalWeakIncObjCutCurrent(BcpsConstraintPool &conPool) } } - cutub += lObjVal; + cutub += lObjVal * (1 + gap/100); // YX: add SL gap to bound assert(indexList.size() == valsList.size()); numCuts += addCut(conPool, cutlb, cutub, indexList, valsList, allowRemoveCut); @@ -5553,7 +5553,7 @@ MibSCutGenerator::incObjCutCurrent(BcpsConstraintPool &conPool) double * upperSol = bS->upperSolutionOrd_; // UL portion from LR soln double * lowerSol = bS->lowerSolutionOrd_; // LL portion from LR soln - double * optLowerSol = bS->optLowerSolutionOrd_; // optimal LL solution + double * optLowerSol = bS->vfLowerSolutionOrd_; // optimal LL solution /* returns a double with the values [alpha | beta | gamma] @@ -5639,6 +5639,7 @@ MibSCutGenerator::bendersInterdictionOneCut(BcpsConstraintPool &conPool, double (localModel_->MibSPar_->entry(MibSParams::branchStrategy)); bool allowRemoveCut(localModel_->MibSPar_->entry(MibSParams::allowRemoveCut)); + double targetGap(localModel_->MibSPar_->entry(MibSParams::slTargetGap)); //when the branching strategy is fractional and the optimal //solution of relaxation is integer, we are forced to generate cut. @@ -5658,6 +5659,7 @@ MibSCutGenerator::bendersInterdictionOneCut(BcpsConstraintPool &conPool, double double * lObjCoeffs = localModel_->getLowerObjCoeffs(); double cutlb(-localModel_->solver()->getInfinity()); double cutub(0.0); + double gap = (targetGap < etol) ? 0.0 : targetGap; // YX: added SL gap std::vector indexList; std::vector valsList; int bigM(10000); @@ -5698,6 +5700,7 @@ MibSCutGenerator::bendersInterdictionOneCut(BcpsConstraintPool &conPool, double } } assert(indexList.size() == valsList.size()); + cutub = cutub * (1 + gap/100); // YX: add SL gap to bound numCuts += addCut(conPool, cutlb, cutub, indexList, valsList, allowRemoveCut); indexList.clear(); @@ -6028,7 +6031,7 @@ MibSCutGenerator::generateConstraints(BcpsConstraintPool &conPool) localModel_->MibSPar_->entry(MibSParams::bendersCutType); if(bendersCutType == MibSBendersCutTypeJustOneCut){ numCuts += bendersInterdictionOneCut(conPool, - bS->optLowerSolutionOrd_); + bS->vfLowerSolutionOrd_); } else{ numCuts += bendersInterdictionMultipleCuts(conPool); @@ -6037,27 +6040,27 @@ MibSCutGenerator::generateConstraints(BcpsConstraintPool &conPool) if (useIntersectionCutTypeIC == PARAM_ON){ cutType = MibSIntersectionCutTypeIC; - numCuts += intersectionCuts(conPool, bS->optLowerSolutionOrd_, cutType); + numCuts += intersectionCuts(conPool, bS->vfLowerSolutionOrd_, cutType); } if (useIntersectionCutTypeWatermelon == PARAM_ON){ cutType = MibSIntersectionCutTypeWatermelon; - numCuts += intersectionCuts(conPool, bS->optLowerSolutionOrd_, cutType); + numCuts += intersectionCuts(conPool, bS->vfLowerSolutionOrd_, cutType); } if (useIntersectionCutTypeHypercubeIC == PARAM_ON){ cutType = MibSIntersectionCutTypeHypercubeIC; - numCuts += intersectionCuts(conPool, bS->optLowerSolutionOrd_, cutType); + numCuts += intersectionCuts(conPool, bS->vfLowerSolutionOrd_, cutType); } if (useIntersectionCutTypeTenderIC == PARAM_ON){ cutType = MibSIntersectionCutTypeTenderIC; - numCuts += intersectionCuts(conPool, bS->optLowerSolutionOrd_, cutType); + numCuts += intersectionCuts(conPool, bS->vfLowerSolutionOrd_, cutType); } if (useIntersectionCutTypeHybridIC == PARAM_ON){ cutType = MibSIntersectionCutTypeHybridIC; - numCuts += intersectionCuts(conPool, bS->optLowerSolutionOrd_, cutType); + numCuts += intersectionCuts(conPool, bS->vfLowerSolutionOrd_, cutType); } if (useGeneralNoGoodCut == PARAM_ON){ @@ -6080,7 +6083,7 @@ MibSCutGenerator::generateConstraints(BcpsConstraintPool &conPool) localModel_->MibSPar_->entry(MibSParams::bendersCutType); if(bendersCutType == MibSBendersCutTypeJustOneCut){ numCuts += bendersInterdictionOneCut(conPool, - bS->optLowerSolutionOrd_); + bS->vfLowerSolutionOrd_); } else{ numCuts += bendersInterdictionMultipleCuts(conPool); @@ -6090,17 +6093,17 @@ MibSCutGenerator::generateConstraints(BcpsConstraintPool &conPool) (relaxedObjVal > localModel_->bS_->objVal_ + localModel_->etol_ || localModel_->MibSPar_->entry(MibSParams::bilevelFreeSetTypeIC) == 1)){ cutType = MibSIntersectionCutTypeIC; - numCuts += intersectionCuts(conPool, bS->optLowerSolutionOrd_, cutType); + numCuts += intersectionCuts(conPool, bS->vfLowerSolutionOrd_, cutType); } if (useIntersectionCutTypeWatermelon == PARAM_ON){ cutType = MibSIntersectionCutTypeWatermelon; - numCuts += intersectionCuts(conPool, bS->optLowerSolutionOrd_, cutType); + numCuts += intersectionCuts(conPool, bS->vfLowerSolutionOrd_, cutType); } if (useIntersectionCutTypeHypercubeIC == PARAM_ON && haveSecondLevelSol){ cutType = MibSIntersectionCutTypeHypercubeIC; - numCuts += intersectionCuts(conPool, bS->optLowerSolutionOrd_, cutType); + numCuts += intersectionCuts(conPool, bS->vfLowerSolutionOrd_, cutType); } if (localModel_->allUpperBin_){ @@ -6122,13 +6125,13 @@ MibSCutGenerator::generateConstraints(BcpsConstraintPool &conPool) if (useIntersectionCutTypeWatermelon == PARAM_ON){ cutType = MibSIntersectionCutTypeWatermelon; - numCuts += intersectionCuts(conPool, bS->optLowerSolutionOrd_, cutType); + numCuts += intersectionCuts(conPool, bS->vfLowerSolutionOrd_, cutType); } if (useIntersectionCutTypeIC == PARAM_ON && ((haveSecondLevelSol && relaxedObjVal > localModel_->bS_->objVal_ + localModel_->etol_) || localModel_->MibSPar_->entry(MibSParams::bilevelFreeSetTypeIC) == 1)){ cutType = MibSIntersectionCutTypeIC; - numCuts += intersectionCuts(conPool, bS->optLowerSolutionOrd_, cutType); + numCuts += intersectionCuts(conPool, bS->vfLowerSolutionOrd_, cutType); } } From 8a3d7e1f27b4039b2d6164a1c623ef8050419117 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Tue, 27 Sep 2022 22:05:29 -0400 Subject: [PATCH 23/26] fix coeff matrix comp in interdiction --- src/MibSModel.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MibSModel.cpp b/src/MibSModel.cpp index cb8bde4e..3eeaffa7 100644 --- a/src/MibSModel.cpp +++ b/src/MibSModel.cpp @@ -3126,8 +3126,8 @@ MibSModel::setCoeffMatrices() matrixG1->setDimensions(0, lowerDim_); } matrixA2G2 = new CoinPackedMatrix(false, 0, 0); - matrixA2G2->setDimensions(0, numVars_); - assert(numVars_== upperDim_ + lowerDim_); // YX: debug only + matrixA2G2->setDimensions(0, upperDim_ + lowerDim_); + // assert(numVars_== upperDim_ + lowerDim_); // YX: debug only matrixA2 = new CoinPackedMatrix(false, 0, 0); matrixA2->setDimensions(0, upperDim_); matrixG2 = new CoinPackedMatrix(false, 0, 0); From 1a3cd4894c9c2969355d65e5f0e09f46ba1d8ee3 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Wed, 2 Nov 2022 23:50:38 -0400 Subject: [PATCH 24/26] update bound on Benders cuts --- src/MibSCutGenerator.cpp | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/MibSCutGenerator.cpp b/src/MibSCutGenerator.cpp index 0122cc4d..4e8d9052 100644 --- a/src/MibSCutGenerator.cpp +++ b/src/MibSCutGenerator.cpp @@ -965,7 +965,7 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, } } - //YX: set cosntraint UB: d^2y >= d^2y^*(1+gap); convert row sense to L + //YX: add cosntraint UB: -d^2y <= - d^2y* - gap|d^2y*| if(targetGap > etol){ for(i = 0; i < lCols; i++){ templObj += lObjSense * lObjCoeff[i] * optLowerSol[i]; // YX: track d^2y^* @@ -1395,7 +1395,7 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo addedRow.clear(); } - // YX: for nonzero gap add an constraint: d^2 \delta y >= -d^2\yhat + d^2y^*(1+gap) + // YX: for nonzero gap add an constraint: d^2\Dy >= -d^2\yhat + d^2y^*+gap|d^2y^*| if(targetGap > etol){ for(i = 0; i < lCols; i++){ addedRow.insert(i, -lObjCoeff[i] * lObjSense); @@ -1484,7 +1484,7 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo nSolver->setRowUpper(2 * i + 1, rhs - lCoeffsTimesLpSol[i]); } - //YX: set cosntraint UB: d^2 \delta y >= -d^2\yhat + d^2y^*(1+gap); convert row sense to L + //YX: add cosntraint UB: -d^2\Dy <= d^2\yhat - d^2y* - gap|d^2y*| if(targetGap > etol){ rhs = 0; for(i = 0; i < lCols; i++){ @@ -4879,7 +4879,8 @@ MibSCutGenerator::generalWeakIncObjCutCurrent(BcpsConstraintPool &conPool) isBigMIncObjSet_ = true; } - bigM = bigMIncObj_ - lObjVal + 1; + bigM = bigMIncObj_ + fabs(bigMIncObj_) * gap/100 + - lObjVal - fabs(lObjVal) * gap/100 + 1; for(i = 0; i < uN; i++){ index = upperColInd[i]; @@ -4907,7 +4908,8 @@ MibSCutGenerator::generalWeakIncObjCutCurrent(BcpsConstraintPool &conPool) } } - cutub += lObjVal * (1 + gap/100); // YX: add SL gap to bound + // YX: add SL gap to bound d^2y* + gap|d^2y*| + cutub += lObjVal + fabs(lObjVal) * gap/100; assert(indexList.size() == valsList.size()); numCuts += addCut(conPool, cutlb, cutub, indexList, valsList, allowRemoveCut); @@ -5658,7 +5660,7 @@ MibSCutGenerator::bendersInterdictionOneCut(BcpsConstraintPool &conPool, double int * lowerColInd = localModel_->getLowerColInd(); double * lObjCoeffs = localModel_->getLowerObjCoeffs(); double cutlb(-localModel_->solver()->getInfinity()); - double cutub(0.0); + double cutub(0.0), tempub(0.0); double gap = (targetGap < etol) ? 0.0 : targetGap; // YX: added SL gap std::vector indexList; std::vector valsList; @@ -5669,7 +5671,7 @@ MibSCutGenerator::bendersInterdictionOneCut(BcpsConstraintPool &conPool, double for(i = 0; i < uN; i++){ indexU = upperColInd[i]; indexL = lowerColInd[i]; - cutub += lObjCoeffs[i] * lSolution[i]; + tempub += lObjCoeffs[i] * lSolution[i]; valU = 0; valL = lObjCoeffs[i]; if(lSolution[i] > etol){ @@ -5700,7 +5702,9 @@ MibSCutGenerator::bendersInterdictionOneCut(BcpsConstraintPool &conPool, double } } assert(indexList.size() == valsList.size()); - cutub = cutub * (1 + gap/100); // YX: add SL gap to bound + + // YX: add SL gap to bound d^2y* + gap|d^2y*| + cutub = tempub + fabs(tempub) * gap/100; numCuts += addCut(conPool, cutlb, cutub, indexList, valsList, allowRemoveCut); indexList.clear(); From aab0d7104327d8f790900068d3de92856d57a299 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Fri, 4 Nov 2022 00:34:12 -0400 Subject: [PATCH 25/26] fix bug in pes feasiblity check --- src/MibSBilevel.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index 8044a116..74327474 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -657,7 +657,8 @@ MibSBilevel::checkBilevelFeasibility(bool isRoot) model_->seenLinkingSolutions[linkSol].lowerSolution.end(), lowerSol); // YX: pessimistic case; retrieve max d^1y only when (PES-MILP) is feasible - if(findPesSol && (tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsFeasible)){ + if(findPesSol){ + assert(tagInSeenLinkingPool_ != MibSLinkingPoolTagPesIsInfeasible); pesRF = model_->seenLinkingSolutions[linkSol].pesRFValue; std::copy(model_->seenLinkingSolutions[linkSol].pesSolution.begin(), model_->seenLinkingSolutions[linkSol].pesSolution.end(), pesSol); From 6cda486dbc9d92d9bd9153c7de2840a6d56b2ef6 Mon Sep 17 00:00:00 2001 From: Yu Date: Thu, 2 Feb 2023 01:37:45 -0500 Subject: [PATCH 26/26] fix obj coeff indices --- src/MibSBilevel.cpp | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index 74327474..8a21455a 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -1786,24 +1786,23 @@ MibSBilevel::getUpperObj(double *lowerSol, double *upperSol) const double uObjSense(model_->solver()->getObjSense()); const double *uObjCoeffs(model_->solver()->getObjCoefficients()); + bool useIdx(false); + int i(0), j(0), idx(0); + double rfuncVal(0.0), uObjVal(0.0); + if(!upperSol){ const double *upperSol(model_->solver()->getColSolution()); + useIdx = true; } - int i(0), idx(0); - double rfuncVal(0.0), uObjVal(0.0); - for(i = 0; i < uN; i++){ - if(!upperSol){ - idx = uColIndices[i]; - }else{ - idx = i; - } + j = (useIdx)? uColIndices[i] : i; + idx = uColIndices[i]; if(0){ - std::cout << "upperSol[" << i << "]: " << upperSol[i] << std::endl; + std::cout << "upperSol[" << j << "]: " << upperSol[j] << std::endl; std::cout << "uObjCoeff[" << idx << "]: " << uObjCoeffs[idx] << std::endl; } - uObjVal += uObjCoeffs[idx] * upperSol[i] * uObjSense; + uObjVal += uObjCoeffs[idx] * upperSol[j] * uObjSense; } rfuncVal = getRiskFuncVal(lowerSol);