Unbounded subproblems for general decomposition problem with dw

[kibaekkim]

The Dantzig-Wolfe decomposition does not address unbounded subproblems and terminates with segfault for general decomposition problem instance.

Example is: https://github.com/Argonne-National-Laboratory/DSP/blob/41c3586948ef5c1653823d5d73f3c780eb4ad3b2/examples/cpp/farmer.cpp

[phumthep]

Has this issue been solved? I used the DW method to solve an MILP and the algorithm does not seem to converge after 68 iterations (kept running for 10 hours+). The DW bound did not decrease.

[kibaekkim]

@phumthep Not sure.. I have to revisit this again.. Is your comment related to #253 ?

[hideakiv]

Example is: https://github.com/Argonne-National-Laboratory/DSP/blob/41c3586948ef5c1653823d5d73f3c780eb4ad3b2/examples/cpp/farmer.cpp

The subproblem model appears to be defined improperly

DSP/examples/cpp/farmer.cpp

Line 151 in 41c3586

0,

[hideakiv]

When some subproblem is dual infeasible, DSP tries to use an extreme ray. However, getPrimalRays is not implemented in most solvers in Osi (except for Cplex as far as I checked).

DSP/src/Solver/DantzigWolfe/DwWorker.cpp

Line 276 in 41c3586

std::vector<double*> rays = osi_[s]->si_->getPrimalRays(1);

If getPrimalRays is implemented, a bundle cut may be generated by using the extreme ray of the dual infeasible subproblem. Otherwise, we would better return a message and terminate DSP gracefully.

[hideakiv]

If getPrimalRays is implemented, a bundle cut may be generated by using the extreme ray of the dual infeasible subproblem.

However, when some subproblems are primal/dual infeasible, the generated columns are not added to the master.

DSP/src/Solver/DantzigWolfe/DwMaster.cpp

Lines 788 to 847 in 41c3586

	/** any subproblem primal/dual infeasible? */
	bool isInfeasible = false;
	for (auto status = status_subs_.begin(); status != status_subs_.end(); status++)
	if (*status == DSP_STAT_PRIM_INFEASIBLE ||
	*status == DSP_STAT_DUAL_INFEASIBLE) {
	isInfeasible = true;
	break;
	}
	if (!isInfeasible) {
	/** calculate lower bound */
	if (phase_ == 2) {
	DSP_RTN_CHECK_RTN_CODE(getLagrangianBound(subobjs));
	DSPdebugMessage("Current lower bound %e, best lower bound %e\n", dualobj_, bestdualobj_);
	}
	/** clear recent subproblem solutions */
	for (unsigned i = 0; i < recent_subsols_.size(); ++i)
	FREE_PTR(recent_subsols_[i]);
	recent_subsols_.clear();
	/** assign the current subproblem solutions */
	recent_subsols_.reserve(subsols.size());
	for (unsigned i = 0; i < subsols.size(); ++i)
	recent_subsols_.push_back(subsols[i]);
	/** create and add columns */
	DSP_RTN_CHECK_RTN_CODE(
	addCols(subinds, status_subs_, subcxs, subobjs, subsols));
	if (model_->isStochastic() &&
	par_->getIntParam("DW/MAX_EVAL_UB") > 0 &&
	par_->getIntParam("DW/EVAL_UB") >= 0) {
	/** maximum number of solutions to evaluate */
	int max_stores = par_->getIntParam("DW/MAX_EVAL_UB");
	/** store solutions to distribute */
	TssModel* tss = dynamic_cast<TssModel*>(model_);
	for (unsigned i = 0; i < subsols.size(); ++i) {
	if (max_stores == 0) break;
	/** assign first-stage solution value */
	CoinPackedVector* first_stage_solution = new CoinPackedVector;
	first_stage_solution->reserve(tss->getNumCols(0));
	for (int j = 0; j < subsols[i]->getNumElements(); ++j) {
	if (subsols[i]->getIndices()[j] >= tss->getNumScenarios() * tss->getNumCols(0))
	break;
	first_stage_solution->insert(subsols[i]->getIndices()[j] % tss->getNumCols(0), subsols[i]->getElements()[j]);
	// printf("i %d first_stage_solution[%d] = %e\n", i, subsols[i]->getIndices()[j] % tss->getNumCols(0), subsols[i]->getElements()[j]);
	}
	/** store it if not duplicated */
	if (!duplicateVector(first_stage_solution, stored_solutions_)) {
	stored_solutions_.push_back(first_stage_solution);
	/** count */
	max_stores--;
	}
	}
	}
	} else
	ngenerated_ = 0;

DwBundleDual needs to handle the case that no cut is added because all the generated cuts are from dual infeasible subproblems. This can possibly cause the master problem to be unbounded due to the bounds defined below, when all generated cuts are from dual infeasible solutions.

DSP/src/Solver/DantzigWolfe/DwBundleDual.cpp

Lines 240 to 242 in 41c3586

	std::fill(clbd.begin(), clbd.begin() + nrows_conv_, -COIN_DBL_MAX);
	std::fill(cubd.begin(), cubd.begin() + nrows_conv_, +COIN_DBL_MAX);
	std::fill(obj.begin(), obj.begin() + nrows_conv_, -1.0);

[hideakiv]

Similarly, when all generated cuts are from dual infeasible solutions, the primal form of the master problem becomes infeasible due to the following row bounds (0=1).

DSP/src/Solver/DantzigWolfe/DwBundleDual.cpp

Lines 167 to 168 in 41c3586

	std::fill(rlbd.begin(), rlbd.begin() + nrows_conv_, 1.0);
	std::fill(rubd.begin(), rubd.begin() + nrows_conv_, 1.0);

[hideakiv]

For the dual form of the master problem, I set the column bounds to be zero.

std::fill(clbd.begin(), clbd.begin() + nrows_conv_, 0.0); 
std::fill(cubd.begin(), cubd.begin() + nrows_conv_, 0.0); 

and modified it to COIN_DBL_MAX and -COIN_DBL_MAX when dual feasible subproblem solution is discovered.

Similarly, for the primal form of the master problem, I set the row bounds to be zero.

 std::fill(rlbd.begin(), rlbd.begin() + nrows_conv_, 0.0); 
 std::fill(rubd.begin(), rubd.begin() + nrows_conv_, 0.0); 

and modified it to 1 when dual feasible subproblem solution is discovered.

[hideakiv]

I am now getting optimal status for the master problem. However, in the second round of solving the subproblem, I got optimal status for one problem but the solution is extremely large, and gets stuck at the assertion below.

DSP/src/Solver/DantzigWolfe/DwWorker.cpp

Line 320 in 41c3586

assert(fabs(xval) < 1e+20);