diff --git a/src/MibSBilevel.cpp b/src/MibSBilevel.cpp index 7189b322..8a21455a 100644 --- a/src/MibSBilevel.cpp +++ b/src/MibSBilevel.cpp @@ -85,19 +85,27 @@ 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]; + } + 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; @@ -220,9 +228,10 @@ MibSBilevel::createBilevel(CoinPackedVector* sol, } //steps 5-6 - if((isLinkVarsFixed_) && ((tagInSeenLinkingPool_ == - MibSLinkingPoolTagLowerIsInfeasible) || - (tagInSeenLinkingPool_ == MibSLinkingPoolTagUBIsSolved))){ + if((isLinkVarsFixed_) && + ((tagInSeenLinkingPool_ == MibSLinkingPoolTagLowerIsInfeasible) || + (tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsInfeasible) || // YX: pessimistic case + (tagInSeenLinkingPool_ == MibSLinkingPoolTagUBIsSolved))){ useBilevelBranching_ = false; shouldPrune_ = true; } @@ -230,6 +239,7 @@ MibSBilevel::createBilevel(CoinPackedVector* sol, //step 7 if(!shouldPrune_){ if((tagInSeenLinkingPool_ == MibSLinkingPoolTagLowerIsFeasible || + tagInSeenLinkingPool_ == MibSLinkingPoolTagPesIsFeasible || // YX: pessimistic case tagInSeenLinkingPool_ == MibSLinkingPoolTagUBIsSolved) || (!isContainedInLinkingPool_ && ((branchPar == MibSBranchingStrategyLinking && isIntegral_ && isLinkVarsFixed_) || @@ -239,7 +249,7 @@ MibSBilevel::createBilevel(CoinPackedVector* sol, (solveSecondLevelWhenXVarsInt && isUpperIntegral_) || (solveSecondLevelWhenLVarsInt && isLinkVarsIntegral_) || (solveSecondLevelWhenLVarsFixed && isLinkVarsFixed_ )))){ - storeSol = checkBilevelFeasiblity(mibs->isRoot_); + storeSol = checkBilevelFeasibility(mibs->isRoot_); } } @@ -257,10 +267,12 @@ MibSBilevel::createBilevel(CoinPackedVector* sol, //############################################################################# MibSSolType -MibSBilevel::checkBilevelFeasiblity(bool isRoot) +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 @@ -279,25 +291,36 @@ MibSBilevel::checkBilevelFeasiblity(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); + 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 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); + 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; // YX: (PES-MILP) solution (x,y) ptr double *lowerSol = new double[lN]; + 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 + uN); // YX: pessimistic case 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]; @@ -330,6 +353,9 @@ MibSBilevel::checkBilevelFeasiblity(bool isRoot) OsiSolverInterface *lSolver = lSolver_; + // lSolver->writeLp("lSolver"); // YX: debug only + // model_->getSolver()->writeLp("oSolver"); // YX: debug only + remainingTime = timeLimit - model_->broker_->subTreeTimer().getTime(); if(remainingTime <= etol){ shouldPrune_ = true; @@ -453,7 +479,10 @@ MibSBilevel::checkBilevelFeasiblity(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 @@ -478,66 +507,244 @@ MibSBilevel::checkBilevelFeasiblity(bool isRoot) node->setLowerUB(objVal); node->setIsBoundSet(true); } + + // YX: pessimistic case; enter only when isUpperIntegral_? + if(findPesSol && isUpperIntegral_){ + if(pSolver_){ + pSolver_ = setUpPesModel(objVal, false); + }else{ + pSolver_ = setUpPesModel(objVal, true); + } + + OsiSolverInterface *pSolver = pSolver_; + + // pSolver->writeLp("pSolverLoaded"); // YX: debug only + + // 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; + model_->counterPES_ ++; + + if((feasCheckSolver == "SYMPHONY") && (sym_is_time_limit_reached + (dynamic_cast + (pSolver)->getSymphonyEnvironment()))){ + shouldPrune_ = true; + goto TERM_CHECKBILEVELFEAS; + } + + if(!pSolver->isProvenOptimal()){ + if(useLinkingSolutionPool && isLinkVarsIntegral_){ + addSolutionToSeenLinkingSolutionPool + (MibSLinkingPoolTagPesIsInfeasible, shouldStoreValuesPesSol, 0.0); + } + if(isLinkVarsFixed_){ + useBilevelBranching_ = false; + shouldPrune_ = true; + } + }else{ + isPesOptimal = true; + pesRF = pSolver->getObjValue(); // YX: obtain max d^1y + valuesPes = pSolver->getColSolution(); + + // 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{ + 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 + uN + lN, shouldStoreValuesPesSol.begin()); + addSolutionToSeenLinkingSolutionPool + (MibSLinkingPoolTagPesIsFeasible, shouldStoreValuesPesSol, pesRF); + shouldStoreValuesPesSol.clear(); + } + } + } } /*if (!warmStartLL){ delete solver_; }*/ - } + } // end of (solving SL-MILP and PES-MILP if not solved before) //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); + + // YX: pessimistic case; retrieve max d^1y only when (PES-MILP) is feasible + if(findPesSol){ + assert(tagInSeenLinkingPool_ != MibSLinkingPoolTagPesIsInfeasible); + pesRF = model_->seenLinkingSolutions[linkSol].pesRFValue; + std::copy(model_->seenLinkingSolutions[linkSol].pesSolution.begin(), + model_->seenLinkingSolutions[linkSol].pesSolution.end(), pesSol); + } } lowerObj = getLowerObj(sol, model_->getLowerObjSense()); if(isIntegral_){ - assert((objVal - lowerObj) <= etol); + assert((objVal - lowerObj) <= etol); } LPSolStatus_ = MibSLPSolStatusInfeasible; //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: [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_)){ + // if(((lowerObj - objVal) < etol) && (isIntegral_)){ + LPSolStatus_ = MibSLPSolStatusFeasible; + useBilevelBranching_ = false; + shouldPrune_ = true; + storeSol = MibSRelaxationSol; + }else{ + memcpy(optLowerSolutionOrd_, lowerSol, sizeof(double) * lN); + if(isUpperIntegral_){ + storeSol = MibSHeurSol; + } + } + }else{ + lpPesVal = getRiskFuncVal(lowerSolutionOrd_); + // 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{ + // YX: copy (PES-MILP) solution to optimal solution container + 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; + } + } + } + } + if(!shouldPrune_){ //step 18 - if((tagInSeenLinkingPool_ != MibSLinkingPoolTagUBIsSolved) && + if(((tagInSeenLinkingPool_ != MibSLinkingPoolTagUBIsSolved) && + (tagInSeenLinkingPool_ != MibSLinkingPoolTagPesIsInfeasible)) && (((branchPar == MibSBranchingStrategyLinking) && (isIntegral_) && (isLinkVarsFixed_)) || ((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); - } + + // YX: if (pes), then (PES-MILP) is solved and a solution is found; + + if(UBSolver_){ + UBSolver_ = setUpUBModel(model_->getSolver(), objVal, false); + }else{ + UBSolver_ = setUpUBModel(model_->getSolver(), objVal, true); + } - OsiSolverInterface * UBSolver = UBSolver_; + OsiSolverInterface * UBSolver = UBSolver_; + + // UBSolver->writeLp("UBSolverLoaded"); // YX: debug only remainingTime = timeLimit - model_->broker_->subTreeTimer().getTime(); @@ -611,7 +818,12 @@ MibSBilevel::checkBilevelFeasiblity(bool isRoot) isUBSolved_ = true; - if (UBSolver->isProvenOptimal()){ + // YX: pess UB has to have an optimal solution in case of feasible (PES-MILP) + if(findPesSol){ //YX: for debug + assert(UBSolver->isProvenOptimal()); + } + + if(UBSolver->isProvenOptimal()){ isProvenOptimal_ = true; const double * valuesUB = UBSolver->getColSolution(); std::copy(valuesUB, valuesUB + uN + lN, shouldStoreValuesUBSol.begin()); @@ -636,33 +848,39 @@ MibSBilevel::checkBilevelFeasiblity(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; } - //step 22 + + //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; //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) TERM_CHECKBILEVELFEAS: delete [] lowerSol; + delete [] pesSol; return storeSol; } @@ -674,14 +892,15 @@ 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_; + if(pSolver_) delete pSolver_; // YX: pessimistic case if(UBSolver_) delete UBSolver_; //delete heuristic_; } - //############################################################################# OsiSolverInterface * MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, @@ -690,6 +909,8 @@ 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; @@ -700,15 +921,19 @@ 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()); 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 + + // YX: d2^y <= phi(A^x) + \delta * abs(phi(A^x)) + objUb = objValLL + fabs(objValLL) * gap/100; if(newOsi){ double objSense(model_->getLowerObjSense()); @@ -722,7 +947,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 origRowNum = findPes ? rowNum - 2: rowNum - 1; CoinPackedMatrix matrix = *model_->origConstCoefMatrix_; matrix.reverseOrdering(); @@ -736,8 +961,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); @@ -790,22 +1015,39 @@ 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] = objValLL; - rowLb[rowNum-1] = -1 * (oSolver->getInfinity()); 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] * uObjSense; + // YX: set up obj in pessimistic case cx+0y + objCoeffs[index1] = 0.0; + } + 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]; @@ -835,7 +1077,10 @@ 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) + if(findPes){ + nSolver->setRowLower(rowNum-2, getRiskFuncVal(optLowerSolutionOrd_)); + } for(i = 0; i < uCols; i++){ index1 = uColIndices[i]; if(fixedInd[index1] == 1){ @@ -850,6 +1095,181 @@ MibSBilevel::setUpUBModel(OsiSolverInterface * oSolver, double objValLL, } +//############################################################################# +OsiSolverInterface * +MibSBilevel::setUpPesModel(double objValLL, bool newOsi, const double *lpSol) +{ + /** Setup pessimistic risk function model (PES-MILP) with fixed linking part **/ + + std::string feasCheckSolver(model_->MibSPar()->entry + (MibSParams::feasCheckSolver)); + double targetGap(model_->MibSPar()->entry(MibSParams::slTargetGap)); + + OsiSolverInterface * nSolver; + + int i(0), idx(0); + double value(0.0), objUb(0.0); + + int uCols(model_->getUpperDim()); + int lCols(model_->getLowerDim()); + int rowNum(model_->getNumOrigCons() + 1); + int colNum(model_->getNumOrigVars()); + int *uColIndices(model_->getUpperColInd()); + int *fixedInd(model_->getFixedInd()); + const double *origRowLb(model_->getOrigRowLb()); + const double *origRowUb(model_->getOrigRowUb()); + double gap = (targetGap < model_->getTolerance()) ? 0.0 : targetGap; + + // YX: d2^y <= phi(A^x) + \delta * abs(phi(A^x)) + objUb = objValLL + fabs(objValLL) * gap/100; + + if (!lpSol){ + lpSol = model_->getSolver()->getColSolution(); + } + + if(newOsi){ + double lObjSense(model_->getLowerObjSense()); + double *lObjCoeffs(model_->getLowerObjCoeffs()); + int *lColIndices(model_->getLowerColInd()); + const double *origColLb(model_->getOrigColLb()); + const double *origColUb(model_->getOrigColUb()); + const double *uObjCoeffs(model_->getSolver()->getObjCoefficients()); + + int intCnt(0); + double uObjSense(1); + + CoinPackedVector row; + CoinPackedMatrix matrix = *model_->origConstCoefMatrix_; + matrix.reverseOrdering(); + + 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"){ + 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, colNum); + CoinZeroN(colLb, colNum); + CoinZeroN(integerVars, colNum); + CoinZeroN(objCoeffs, colNum); + CoinZeroN(newRow, lCols); + + /** Set row bounds **/ + memcpy(rowLb, origRowLb, sizeof(double) * (rowNum - 1)); + memcpy(rowUb, origRowUb, sizeof(double) * (rowNum - 1)); + + /** Set col bounds **/ + 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]; + } + } + + /** Fill in array of integer variables **/ + for(i = 0; i < colNum; i++){ + if(model_->getSolver()->isInteger(i)){ + integerVars[intCnt] = i; + intCnt ++; + } + } + + /** Prepare for value function constraint **/ + CoinDisjointCopyN(lObjCoeffs, lCols, newRow); + + rowUb[rowNum-1] = objUb; + rowLb[rowNum-1] = -1 * (model_->getSolver()->getInfinity()); + + /** Fill in new objective and constraint coeffs **/ + for(i = 0; i < lCols; i++){ + idx = lColIndices[i]; + 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); + + for(i = 0; i < intCnt; i++){ + nSolver->setInteger(integerVars[i]); + } + + /** Pessimistic risk function **/ + nSolver->setObjSense(-1 * uObjSense); //1 min; -1 max + + // YX: still need this? + nSolver->setHintParam(OsiDoReducePrint, true, OsiHintDo); + + delete [] colUb; + delete [] colLb; + delete [] rowUb; + delete [] rowLb; + delete [] objCoeffs; + delete [] integerVars; + delete [] newRow; + delete newMat; + + }else{ + + nSolver = pSolver_; + + /** Update value function constraint bound **/ + nSolver->setRowUpper(rowNum-1, objUb); // YX: gap added to phi(A^2x) + + /** Update fixed linking variable 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); + } + } + } + + return nSolver; +} + //############################################################################# OsiSolverInterface * MibSBilevel::setUpModel(OsiSolverInterface * oSolver, bool newOsi, @@ -1329,6 +1749,68 @@ MibSBilevel::getLowerObj(const double * sol, double objSense) return objVal * objSense; } + +//############################################################################# +double +MibSBilevel::getRiskFuncVal(double *lowerSol) +{ + int lN(model_->getLowerDim()); + int *lColIndices(model_->getLowerColInd()); + const double uObjSense(model_->solver()->getObjSense()); + const double *uObjCoeffs(model_->solver()->getObjCoefficients()); + + int i(0), idx(0); + double rfVal(0.0); + + 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; + } + rfVal += uObjCoeffs[idx] * lowerSol[i]; + } + + return rfVal * uObjSense; + +} + +//############################################################################# +double +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()); + + 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; + } + + for(i = 0; i < uN; i++){ + j = (useIdx)? uColIndices[i] : i; + idx = uColIndices[i]; + if(0){ + std::cout << "upperSol[" << j << "]: " << upperSol[j] << std::endl; + std::cout << "uObjCoeff[" << idx << "]: " << uObjCoeffs[idx] << std::endl; + } + uObjVal += uObjCoeffs[idx] * upperSol[j] * uObjSense; + } + + rfuncVal = getRiskFuncVal(lowerSol); + uObjVal += rfuncVal; + + return uObjVal; +} + //############################################################################# void MibSBilevel::addSolutionToSeenLinkingSolutionPool(MibSLinkingPoolTag solTag, std::vector @@ -1353,8 +1835,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){ @@ -1362,8 +1846,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; @@ -1372,12 +1858,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; following one tag above; + { + 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 c24bcbd8..e35cb3fc 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 **/ @@ -68,11 +71,13 @@ 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_; OsiSolverInterface * lSolver_; + OsiSolverInterface * pSolver_; // YX: pessimistic case OsiSolverInterface * UBSolver_; CoinWarmStart * ws_; @@ -91,9 +96,11 @@ class MibSBilevel { //optLowerSolution_ = 0; optUpperSolutionOrd_ = 0; optLowerSolutionOrd_ = 0; + vfLowerSolutionOrd_ = 0; // YX: SL-MILP solutions; for gap > 0 model_ = 0; heuristic_= 0; lSolver_ = 0; + pSolver_ = 0; // YX: pessimistic case UBSolver_ = 0; ws_ = 0; } @@ -104,17 +111,21 @@ class MibSBilevel { MibSSolType createBilevel(CoinPackedVector *sol, MibSModel *mibs=0); - MibSSolType checkBilevelFeasiblity(bool isRoot); + MibSSolType checkBilevelFeasibility(bool isRoot); void gutsOfDestructor(); private: - int findIndex(int index, int size, int * indices); OsiSolverInterface * setUpUBModel(OsiSolverInterface * solver, double objValLL, bool newOsi, const double *sol = NULL); + 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); // 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;} diff --git a/src/MibSConstants.hpp b/src/MibSConstants.hpp index fb0ac47d..a77295f4 100644 --- a/src/MibSConstants.hpp +++ b/src/MibSConstants.hpp @@ -55,9 +55,11 @@ enum MibSLPSolStatus{ //############################################################################# enum MibSLinkingPoolTag{ - MibSLinkingPoolTagIsNotSet = -4, + MibSLinkingPoolTagIsNotSet = -6, MibSLinkingPoolTagLowerIsInfeasible, MibSLinkingPoolTagLowerIsFeasible, + MibSLinkingPoolTagPesIsInfeasible, // YX: pessimistic case + MibSLinkingPoolTagPesIsFeasible, // YX: pessimistic case MibSLinkingPoolTagUBIsSolved }; diff --git a/src/MibSCutGenerator.cpp b/src/MibSCutGenerator.cpp index f6750c76..4e8d9052 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; } @@ -572,7 +573,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 +611,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 +813,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,25 +822,28 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, (MibSParams::maxThreadsLL)); int whichCutsLL(localModel_->MibSPar_->entry (MibSParams::whichCutsLL)); + double targetGap(localModel_->MibSPar_->entry(MibSParams::slTargetGap)); bool foundSolution(false); 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; 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()); 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()); @@ -887,15 +892,7 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, CoinZeroN(integerVars, newNumCols); if(!localModel_->getA2Matrix()){ - getA2Matrix = true; - } - - if(!localModel_->getG2Matrix()){ - getG2Matrix = true; - } - - if((getA2Matrix) || (getG2Matrix)){ - getLowerMatrices(false, getA2Matrix, getG2Matrix); + localModel_->setCoeffMatrices(); } //extracting optimal first-level solution of the relaxation problem and @@ -935,6 +932,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 +964,14 @@ MibSCutGenerator::findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, assert(0); } } + + //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^* + } + newRowUb[newNumRows-1] = -templObj - (fabs(templObj) * gap/100); + } //filling col bounds for(i = 0; i < lCols; i++){ @@ -1084,7 +1098,7 @@ 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"); - } + } } delete [] multA2XOpt; @@ -1130,20 +1144,15 @@ MibSCutGenerator::getAlphaIC(double** extRay, double* uselessIneqs, char *rowSense = localModel_->getOrigRowSense(); double *lObjCoeffs(localModel_->getLowerObjCoeffs()); double objSense(localModel_->getLowerObjSense()); - bool getA2Matrix(false), getG2Matrix(false); - - if(localModel_->getA2Matrix() == NULL){ - getA2Matrix = true; - } - if(localModel_->getG2Matrix() == NULL){ - getG2Matrix = true; - } + 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((getA2Matrix) || (getG2Matrix)){ - getLowerMatrices(false, getA2Matrix, getG2Matrix); + if(!localModel_->getA2Matrix()){ + localModel_->setCoeffMatrices(); } - + CoinPackedMatrix *matrixA2 = localModel_->getA2Matrix(); CoinPackedMatrix *matrixG2 = localModel_->getG2Matrix(); @@ -1181,6 +1190,12 @@ 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)){ + rhs[lRows] += -fabs(templObj) * gap/100; } for (i = 0; i < numNonBasic; i++){ @@ -1276,7 +1291,7 @@ MibSCutGenerator::solveModelIC(double *uselessIneqs, double *ray, double *rhs, //############################################################################# bool MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lowerLevelSol, - double* lpSol, bool &isTimeLimReached) + double *optLowerSol, double* lpSol, bool &isTimeLimReached) { std::string feasCheckSolver(localModel_->MibSPar_->entry (MibSParams::feasCheckSolver)); @@ -1284,22 +1299,26 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo (MibSParams::maxThreadsLL)); int whichCutsLL(localModel_->MibSPar_->entry (MibSParams::whichCutsLL)); + double targetGap(localModel_->MibSPar_->entry(MibSParams::slTargetGap)); + double timeLimit(localModel_->AlpsPar()->entry(AlpsParams::timeLimit)); double remainingTime(0.0); bool foundSolution = false; 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()); int numContCols(lRows + 2 * lCols); int newNumCols(lCols + numContCols); - int newNumRows(2 * lRows + 2 * lCols + 1); + int newNumRows = (targetGap > etol)? (2 * lRows + 2 * lCols + 2) : (2 * lRows + 2 * lCols + 1); // YX: nonzero gap add a constraint double lObjSense(localModel_->getLowerObjSense()); double *lObjCoeff(localModel_->getLowerObjCoeffs()); int *lColInd(localModel_->getLowerColInd()); @@ -1309,6 +1328,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; @@ -1321,15 +1341,7 @@ 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); + localModel_->setCoeffMatrices(); } CoinPackedMatrix *matrixA2G2 = localModel_->getLowerConstCoefMatrix(); @@ -1383,6 +1395,15 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo addedRow.clear(); } + // 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); + } + newMatrix->appendRow(addedRow); + addedRow.clear(); + } + //filling row bounds CoinFillN(newRowLb, newNumRows, -1 * infinity); @@ -1462,6 +1483,18 @@ MibSCutGenerator::findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lo } nSolver->setRowUpper(2 * i + 1, rhs - lCoeffsTimesLpSol[i]); } + + //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++){ + colIndex = lColInd[i]; + rhs += lObjSense * lObjCoeff[i] * (lpSol[colIndex] - optLowerSol[i]); + templObj += lObjSense * lObjCoeff[i] * optLowerSol[i]; // YX: track d^2y^* + } + rhs += -fabs(templObj) * gap/100; + nSolver->setRowUpper(newNumRows-1, rhs); + } //modifying col bounds for(i = 0; i < lCols; i++){ @@ -1603,11 +1636,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]; } @@ -1725,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_){ @@ -1813,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]); } } @@ -3222,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; } @@ -3294,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_; @@ -4256,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_; @@ -4811,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); @@ -4819,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()); @@ -4828,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; @@ -4849,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]; @@ -4877,7 +4908,8 @@ MibSCutGenerator::generalWeakIncObjCutCurrent(BcpsConstraintPool &conPool) } } - cutub += lObjVal; + // 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); @@ -5523,7 +5555,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] @@ -5609,6 +5641,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. @@ -5627,7 +5660,8 @@ 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; int bigM(10000); @@ -5637,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){ @@ -5668,6 +5702,9 @@ MibSCutGenerator::bendersInterdictionOneCut(BcpsConstraintPool &conPool, double } } assert(indexList.size() == valsList.size()); + + // 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(); @@ -5998,7 +6035,7 @@ MibSCutGenerator::generateConstraints(BcpsConstraintPool &conPool) localModel_->MibSPar_->entry(MibSParams::bendersCutType); if(bendersCutType == MibSBendersCutTypeJustOneCut){ numCuts += bendersInterdictionOneCut(conPool, - bS->optLowerSolutionOrd_); + bS->vfLowerSolutionOrd_); } else{ numCuts += bendersInterdictionMultipleCuts(conPool); @@ -6007,27 +6044,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){ @@ -6050,7 +6087,7 @@ MibSCutGenerator::generateConstraints(BcpsConstraintPool &conPool) localModel_->MibSPar_->entry(MibSParams::bendersCutType); if(bendersCutType == MibSBendersCutTypeJustOneCut){ numCuts += bendersInterdictionOneCut(conPool, - bS->optLowerSolutionOrd_); + bS->vfLowerSolutionOrd_); } else{ numCuts += bendersInterdictionMultipleCuts(conPool); @@ -6060,17 +6097,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_){ @@ -6092,13 +6129,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); } } @@ -6504,6 +6541,7 @@ MibSCutGenerator::getBindingConsBasis() } //############################################################################# +/** void MibSCutGenerator::getLowerMatrices(bool getLowerConstCoefMatrix, bool getA2Matrix, bool getG2Matrix) @@ -6601,4 +6639,4 @@ MibSCutGenerator::getLowerMatrices(bool getLowerConstCoefMatrix, localModel_->setG2Matrix(matrixG2); } -} +}**/ diff --git a/src/MibSCutGenerator.hpp b/src/MibSCutGenerator.hpp index 32a7284b..39c0c4f3 100644 --- a/src/MibSCutGenerator.hpp +++ b/src/MibSCutGenerator.hpp @@ -85,8 +85,8 @@ class MibSCutGenerator : public BlisConGenerator { int intersectionCuts(BcpsConstraintPool &conPool, double *optLowerSolution, MibSIntersectionCutType ICType); /** Helper function for IC*/ - bool findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, const double *sol, - bool &isTimeLimReached); + bool findLowerLevelSol(double *uselessIneqs, double *lowerLevelSol, double *optLowerSol, + const double *sol, bool &isTimeLimReached); // YX: add y^* for nonzero gap case /** Helper function for IC*/ bool getAlphaIC(double** extRay, double *uselessIneqs, double* lowerSolution, int numStruct, @@ -97,7 +97,7 @@ class MibSCutGenerator : public BlisConGenerator { /** Helper function for watermelon IC **/ bool findLowerLevelSolWatermelonIC(double *uselessIneqs, double *lowerLevelSol, - double* lpSol, bool &isTimeLimReached); + double *optLowerSol, double* lpSol, bool &isTimeLimReached); // YX: add y^* for nonzero gap case /** Helper function for watermelon IC*/ bool getAlphaWatermelonIC(double** extRay, double *uselessIneqs, double* lowerSolution, @@ -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/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/MibSModel.cpp b/src/MibSModel.cpp index 0003029b..3eeaffa7 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_; @@ -132,8 +134,10 @@ MibSModel::initialize() structRowNum_ = 0; sizeFixedInd_ = 0; counterVF_ = 0; + counterPES_ = 0; counterUB_ = 0; timerVF_ = 0.0; + timerPES_ = 0.0; // YX: pessmistic case timerUB_ = 0.0; countIteration_ = 0; isInterdict_ = false; @@ -170,6 +174,8 @@ MibSModel::initialize() interdictCost_ = NULL; origConstCoefMatrix_ = NULL; lowerConstCoefMatrix_ = NULL; + A1Matrix_ = NULL; + G1Matrix_ = NULL; A2Matrix_ = NULL; G2Matrix_ = NULL; boundProbRoot_ = NULL; @@ -1048,6 +1054,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(); + } + } //############################################################################# @@ -1995,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, @@ -3081,6 +3095,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, 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); + 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() @@ -3745,5 +3866,20 @@ MibSModel::instanceStructure() std::cout << "Linking solution pool will not be used." << std::endl; } } + + //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 "; + std::cout << MibSPar_->entry(MibSParams::slTargetGap) <<"%."<< std::endl; + } + } + + //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..0e3ef986 100644 --- a/src/MibSModel.hpp +++ b/src/MibSModel.hpp @@ -99,12 +99,18 @@ 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_; /** Time for solving (VF) **/ double timerVF_; + /** Time for solving (PES) **/ + double timerPES_; // YX: pessimistic case; + /** Time for solving (UB) **/ double timerUB_; @@ -156,7 +162,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 +224,18 @@ class MibSModel : public BlisModel { /** Original matrix of constraints coefficients **/ CoinPackedMatrix *origConstCoefMatrix_; + // 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') **/ 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 +416,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 +516,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 +630,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(); diff --git a/src/MibSParams.cpp b/src/MibSParams.cpp index 998dbb27..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 //-------------------------------------------------------- @@ -244,6 +246,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 } //############################################################################# @@ -279,6 +283,8 @@ MibSParams::setDefaultEntries() { setEntry(useNewPureIntCut, false); + setEntry(findPesSol, false); + //------------------------------------------------------------- // Int Parameters. //------------------------------------------------------------- @@ -386,6 +392,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..0cfd90c3 100644 --- a/src/MibSParams.hpp +++ b/src/MibSParams.hpp @@ -39,6 +39,7 @@ class MibSParams : public AlpsParameterSet { printProblemInfo, allowRemoveCut, useNewPureIntCut, + findPesSol, // YX: pessimistic case endOfBoolParams }; @@ -99,6 +100,7 @@ class MibSParams : public AlpsParameterSet { /** Double parameters. */ enum dblParams{ boundCutTimeLim, + slTargetGap, //YX: lower/second level optimality gap endOfDblParams }; 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; }