C++ OsiSolverInterface::isProvenPrimalInfeasible方法代码示例

// Override OsiSolver::solve() because Soplex doesn't have branch & bound
int OsiSpxSolver::solve() {
	OsiSolverInterface* si = OsiSolver::getSolver();
	if(!OsiSolver::prepareSolve()) {
		return UNSAT;
	}

	timer.reset();
	try {
		si->initialSolve();
	} catch (CoinError err) {
		err.print(true);
		return UNSAT;
	}
	time = timer.timeElapsed();

	if (si->isProvenOptimal())
		return SAT;
	else if (si->isProvenPrimalInfeasible())
		return UNSAT;
	else
		return UNKNOWN;
}

//.........这里部分代码省略.........
	int whichWay;
        if (bestColumn >= 0) {
	    bestColumnValue = newSolution[bestColumn];
            if (bestRound < 0) {
                originalBoundBestColumn = upper[bestColumn];
                solver->setColUpper(bestColumn, floor(bestColumnValue));
		whichWay=0;
            } else {
                originalBoundBestColumn = lower[bestColumn];
                solver->setColLower(bestColumn, ceil(bestColumnValue));
		whichWay=1;
            }
        } else {
            break;
        }
        int originalBestRound = bestRound;
        int saveModelOptions = model_->specialOptions();
	
        while (1) {

            model_->setSpecialOptions(saveModelOptions | 2048);
            solver->resolve();
            model_->setSpecialOptions(saveModelOptions);
            if (!solver->isAbandoned()&&!solver->isIterationLimitReached()) {
                numberSimplexIterations += solver->getIterationCount();
            } else {
                numberSimplexIterations = maxSimplexIterations + 1;
		reasonToStop += 100;
                break;
            }

            if (!solver->isProvenOptimal()) {
	        if (nodes) {
		  if (solver->isProvenPrimalInfeasible()) {
		    if (maxSimplexIterationsAtRoot_!=COIN_INT_MAX) {
		      // stop now
		      printf("stopping on first infeasibility\n");
		      break;
		    } else if (cuts) {
		      // can do conflict cut
		      printf("could do intermediate conflict cut\n");
		      bool localCut;
		      OsiRowCut * cut = model_->conflictCut(solver,localCut);
		      if (cut) {
			if (!localCut) {
			  model_->makePartialCut(cut,solver);
			  cuts[numberCuts++]=cut;
			} else {
			  delete cut;
			}
		      }
		    }
		  } else {
		    reasonToStop += 10;
		    break;
		  }
		}
                if (numberAtBoundFixed > 0) {
                    // Remove the bound fix for variables that were at bounds
                    for (int i = 0; i < numberAtBoundFixed; i++) {
                        int iColFixed = columnFixed[i];
                        if (fixedAtLowerBound[i])
                            solver->setColUpper(iColFixed, originalBound[i]);
                        else
                            solver->setColLower(iColFixed, originalBound[i]);
                    }

  /* After a CbcModel::resolve this can return a status
     -1 no effect
     0 treat as optimal
     1 as 0 but do not do any more resolves (i.e. no more cuts)
     2 treat as infeasible
  */
  int
  CbcNlpStrategy::status(CbcModel * model, CbcNodeInfo * parent,int whereFrom)
  {

    OsiSolverInterface * solver = model->solver();//get solver
    int feasible = 1;
    bool solved = true;
    int returnStatus = -1;
    BonCbcPartialNodeInfo * bmNodeInfo = dynamic_cast<BonCbcPartialNodeInfo *>(parent);
    if (!bmNodeInfo) return -1;

    int seqOfInfeasiblesSize = bmNodeInfo->getSequenceOfInfeasiblesSize();
    int seqOfUnsolvedSize = bmNodeInfo->getSequenceOfUnsolvedSize();


    if (solver->isAbandoned()) {
      solved = false;
      seqOfUnsolvedSize++;
      ;
    }
    else if (solver->isProvenPrimalInfeasible()) {
      feasible = 0;
      seqOfInfeasiblesSize++;
    }

    if ((seqOfUnsolvedSize==0) || (maxFailure_ == 0) &&
        (maxInfeasible_== 0) || (seqOfInfeasiblesSize==0))

      if (feasible && seqOfInfeasiblesSize > 1) {
        (*model->messageHandler())<<"Feasible node while father was infeasible."
        <<CoinMessageEol;
      }

    if (solved && seqOfUnsolvedSize > 1) {
      (*model->messageHandler())<<"Solved node while father was unsolved."
      <<CoinMessageEol;
    }

    if (seqOfInfeasiblesSize < maxInfeasible_ &&
        solved && !feasible) {
      (*model->messageHandler())<<"Branching on infeasible node, sequence of infeasibles size "
      <<seqOfInfeasiblesSize<<CoinMessageEol;
      // Have to make sure that we will branch
      OsiTMINLPInterface * ipopt = dynamic_cast<OsiTMINLPInterface *>(solver);
      ipopt->forceBranchable();
      //change objective value
      returnStatus = 0;

    }

    if (!solved && parent != NULL &&
        seqOfUnsolvedSize <= maxFailure_) {
      (*model->messageHandler())<<"Branching on unsolved node, sequence of unsolved size "<<seqOfUnsolvedSize<<CoinMessageEol;
      // Have to make sure that we will branch
      OsiTMINLPInterface * osiMinlp = dynamic_cast<OsiTMINLPInterface *>(solver);
      osiMinlp->forceBranchable();     //      feasible=1;
      returnStatus = 0;
    }

    if (solver->isAbandoned() && parent != NULL &&
        seqOfUnsolvedSize > maxFailure_) {
      hasFailed_ = true;
      OsiTMINLPInterface * osiMinlp =
        dynamic_cast<OsiTMINLPInterface *>(solver);
      if (pretendFailIsInfeasible_) {
        //force infeasible
        osiMinlp->forceInfeasible();
        returnStatus = 2;
      }
      else {
        std::string probName;
        osiMinlp->getStrParam(OsiProbName,probName);
        throw osiMinlp->newUnsolvedError(0, osiMinlp->problem(), probName);
      }
    }
    return returnStatus;
  }

  TNLPSolver::ReturnStatus LpBranchingSolver::
  solveFromHotStart(OsiTMINLPInterface* tminlp_interface)
  {
    TNLPSolver::ReturnStatus retstatus = TNLPSolver::solvedOptimal;

    // updated the bounds of the linear solver
    std::vector<int> diff_low_bnd_index;
    std::vector<double> diff_low_bnd_value;
    std::vector<int> diff_up_bnd_index;
    std::vector<double> diff_up_bnd_value;

    // Get the bounds.  We assume that the bounds in the linear solver
    // are always the original ones
    const int numCols = tminlp_interface->getNumCols();
    const double* colLow_orig = lin_->getColLower();
    const double* colUp_orig = lin_->getColUpper();
    const double* colLow = tminlp_interface->getColLower();
    const double* colUp = tminlp_interface->getColUpper();

    OsiSolverInterface * lin = lin_;
    // eventualy clone lin_
    if(warm_start_mode_ == Clone){
      lin = lin_->clone();
//      std::cout<<"Cloning it"<<std::endl;
    }
    // Set the bounds on the LP solver according to the changes in
    // tminlp_interface
    for (int i=0; i<numCols; i++) {
      const double& lo = colLow[i];
      if (colLow_orig[i] < lo) {
        if(warm_start_mode_ == Basis){
          diff_low_bnd_value.push_back(colLow_orig[i]);
          diff_low_bnd_index.push_back(i);
        }
        lin->setColLower(i,lo);
      }
      const double& up = colUp[i];
      if (colUp_orig[i] > up) {
        if(warm_start_mode_ == Basis){
          diff_up_bnd_index.push_back(i);
          diff_up_bnd_value.push_back(colUp_orig[i]);
        }
        lin->setColUpper(i,lo);
      }
    }

    if(warm_start_mode_ == Basis){
      lin->setWarmStart(warm_);
    }

    lin->resolve();

    double obj = lin->getObjValue();
    bool go_on = true;
    if (lin->isProvenPrimalInfeasible() || 
        lin->isDualObjectiveLimitReached()) {
      retstatus = TNLPSolver::provenInfeasible;
      go_on = false;
    }
    else if (lin->isIterationLimitReached()) {
      retstatus = TNLPSolver::iterationLimit;
      go_on = false;
    }
    else {
      if (maxCuttingPlaneIterations_ > 0 && go_on) {
        double violation;
        obj = ecp_->doEcpRounds(*lin, true, &violation);
        if (obj == COIN_DBL_MAX) {
          retstatus = TNLPSolver::provenInfeasible;
        }
        else if (violation <= 1e-8) {
          retstatus = TNLPSolver::solvedOptimal;
        }
      }
    }
    tminlp_interface->problem()->set_obj_value(obj);
    tminlp_interface->problem()->Set_x_sol(numCols, lin_->getColSolution());

    //restore the original bounds
    if(warm_start_mode_ == Basis){
      for (unsigned int i = 0; i < diff_low_bnd_index.size(); i++) {
        lin_->setColLower(diff_low_bnd_index[i],diff_low_bnd_value[i]);
      }
      for (unsigned int i = 0; i < diff_up_bnd_index.size(); i++) {
        lin_->setColUpper(diff_up_bnd_index[i],diff_up_bnd_value[i]);
      }
    }
    else {
      delete lin;
    }
    return retstatus;
  }

//.........这里部分代码省略.........
    std::vector<int>          A_row_index_coin(A_row_index,A_row_index+nnz);
    
    // declare the solver
    OsiSolverInterface* pSolver;
// 	initialize the solver
    if ( isMIP ) {
        pSolver = new OsiCbcSolverInterface;
    } else {
        pSolver = new OsiClpSolverInterface;
    }
   
//	OsiCbcSolverInterface is deprecated and CbcModel should be used instead but don't
//	know how to get that working with loadProblem. 
//     OsiCbcSolverInterface solver1;
//     CbcModel model(solver1);
//     CbcMain0(model);
//     OsiSolverInterface * pSolver = model.solver();
    
	if (nrhs>12) { // get stuff out of the options structure if provided
		// Finish me
	}
    
//     mexPrintf("Setting Log Level to 0.\n");
	// load the problem
	mexPrintf("Loading the problem.\n");
	pSolver->loadProblem( n_vars, n_cons, // problem size
						  &A_col_starts_coin[0], &A_row_index_coin[0], A_data, // the A matrix
						  x_lb,  x_ub, c, // the objective and bounds
						  Ax_lb, Ax_ub ); // the constraint bounds
    
//     pSolver->messageHandler()->setLogLevel(0); // This doesn't seem to work
    pSolver->setHintParam(OsiDoReducePrint,true,OsiHintTry);
    
	// deal with integer inputs
	if ( isMIP ) {
		for(int i=0;i<n_vars;i++) {
			if (isinteger[i]) pSolver->setInteger(i);
		}
	}
	if (isQP) {
		error("QP is not working yet");
		// need to call loadQuadraticObjective here ???
	}
    
//     CbcModel model(pSolver);
//     model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
    
    
	// solve the problem
	//mexPrintf("Trying to solve the problem.\n");
    if (isMIP) {
        pSolver->branchAndBound();
//         model.branchAndBound();
    } else {
//         model.initialSolve();
        pSolver->initialSolve();
    }
	
	// Allocate memory for return data
    plhs[0] = mxCreateDoubleMatrix(n_vars,1, mxREAL); // for the solution
    plhs[1] = mxCreateDoubleMatrix(n_cons,1, mxREAL); // for the constraint prices
    plhs[2] = mxCreateDoubleMatrix(1,1, mxREAL);      // for the return status
    double *x = mxGetPr(plhs[0]);
    double *y = mxGetPr(plhs[1]);
    double *returncode = mxGetPr(plhs[2]);
	
	// Copy solutions if available
	if ( pSolver->isProvenOptimal() ) {
//     if ( model.isProvenOptimal() ) {
//     if ( model.solver()->isProvenOptimal() ) {
        //mexPrintf("Solution found.\n");
		// extract the solutions
		const double * solution = pSolver->getColSolution();
		const double * dualvars = pSolver->getRowPrice();
		// copy the solution to the outpus
		memcpy(x,solution,n_vars*sizeof(double));
		memcpy(y,dualvars,n_cons*sizeof(double));
		*returncode = 1;
	} else {
		if ( pSolver->isProvenPrimalInfeasible() ) {
			mexPrintf("Primal problem is proven infeasible.\n");
			*returncode = 0;
		} else if ( pSolver->isProvenDualInfeasible() ) {
			mexPrintf("Dual problem is proven infeasible.\n");
			*returncode = -1;
		} else if ( pSolver->isPrimalObjectiveLimitReached() ) {
			mexPrintf("The primal objective limit was reached.\n");
			*returncode = -2;
		} else if ( pSolver->isDualObjectiveLimitReached() ) {
			mexPrintf("The dual objective limit was reached.\n");
			*returncode = -3;
		} else if ( pSolver->isIterationLimitReached() ) {
			mexPrintf("The iteration limit was reached\n");
			*returncode = -4;
		}
	}
	// clean up memory
	if ( mexprinter!= NULL) delete mexprinter;	
	delete pSolver;
}

//Solver function
int sci_rmps(char *fname) 
{
    //creating a problem pointer using base class of OsiSolverInterface and
    //instantiate the object using derived class of ClpSolverInterface
    OsiSolverInterface* si = new OsiClpSolverInterface();

    // Error management variable
	SciErr sciErr;

	//data declarations
	int *piAddressVarOne = NULL;                 //pointer used to access argument of the function
	char* ptr;                              	 //pointer to point to address of file name
    double* options_;                            //options to set maximum iterations 
	CheckInputArgument(pvApiCtx, 2,2 );          //Check we have exactly two arguments as input or not
	CheckOutputArgument(pvApiCtx, 6, 6);         //Check we have exactly six arguments on output side or not
    //Getting the input arguments from Scilab
    //Getting the MPS file path
	//Reading mps file
	getStringFromScilab(1,&ptr);

 	std::cout<<ptr;
	
    //get options from Scilab
    if(getFixedSizeDoubleMatrixInList(2 , 2 , 1 , 1 , &options_))
	{
		return 1;
	}

    //Read the MPS file
    si->readMps(ptr);

    //setting options for maximum iterations
    si->setIntParam(OsiMaxNumIteration,options_[0]);

    //Solve the problem
    si->initialSolve();
  
    //Quering about the problem
    //get number of variables
    double numVars_;
    numVars_ = si->getNumCols();
  
    //get number of constraint equations
    double numCons_;
    numCons_ = si->getNumRows();
   
    //Output the solution to Scilab
    //get solution for x
    const double* xValue = si->getColSolution();
   
    //get objective value
    double objValue = si->getObjValue();

    //get Status value
    double status;
    if(si->isProvenOptimal())
    	status=0;
    else if(si->isProvenPrimalInfeasible())
    	status=1;
    else if(si->isProvenDualInfeasible())
        status=2;
    else if(si->isIterationLimitReached())
        status=3;
   	else if(si->isAbandoned())
        status=4;
   	else if(si->isPrimalObjectiveLimitReached())
        status=5;
   	else if(si->isDualObjectiveLimitReached())
        status=6;

    //get number of iterations
    double iterations = si->getIterationCount();

    //get reduced cost 
    const double* reducedCost = si->getReducedCost();
   
    //get dual vector
    const double* dual = si->getRowPrice();
  
    returnDoubleMatrixToScilab(1 , 1 , numVars_ , xValue);
    returnDoubleMatrixToScilab(2 , 1 , 1 , &objValue);
    returnDoubleMatrixToScilab(3 , 1 , 1 , &status);
    returnDoubleMatrixToScilab(4 , 1 , 1 , &iterations);
    returnDoubleMatrixToScilab(5 , 1 , numVars_ , reducedCost);
    returnDoubleMatrixToScilab(6 , 1 , numCons_ , dual);
	
	free(xValue);
	free(dual);
	free(reducedCost);
}