C++ IloModel::remove方法代码示例

/** Separation function.
 * This function is invoked whenever CPLEX finds an integer feasible
 * solution. It then separates either feasibility or optimality cuts
 * on this solution.
 */
void BendersOpt::LazyConstraintCallback::main() {
   std::cout << "Callback invoked. Separate Benders cuts." << std::endl;

   IloNumArray x(getEnv());
   IloNumArray rayVals(getEnv());
   IloNumVarArray rayVars(getEnv());
   IloNumArray cutVal(getEnv());
   IloNumVarArray cutVar(getEnv());

   getValues(x, master->vars);

   bool error = false;

   // Iterate over blocks and trigger a separation on each of them.
   // The separation is triggered asynchronously so that it can happen
   // on different remote objects simultaneously.
   for (BlockVector::size_type b = 0; b < blocks.size(); ++b) {
      Block *const block = blocks[b];

      // Remove current objective from the block's model.
      IloModel model = block->cplex.getModel();
      IloObjective obj = block->obj;
      model.remove(obj);
      IloExpr newObj = obj.getExpr();
         
      // Iterate over the fixed master variables in this block to update
      // the block's objective function.
      // Each fixed variable goes to the right-hand side and therefore
      // into the objective function.
      for (std::vector<Block::FixData>::const_iterator it = block->fixed.begin(); it != block->fixed.end(); ++it)
         newObj -= block->vars[it->row] * (it->val * x[it->col]);
      obj.setExpr(newObj);
      model.add(obj);
      newObj.end();

      // If the problem is unbounded we need to get an infinite ray in
      // order to be able to generate the respective Benders cut. If
      // CPLEX proves unboundedness in presolve then it will return
      // CPX_STAT_INForUNBD and no ray will be available. So we need to
      // disable presolve.
      block->cplex.setParam(IloCplex::Param::Preprocessing::Presolve, false);
      block->cplex.setParam(IloCplex::Param::Preprocessing::Reduce, 0);

      // Solve the updated problem to optimality.
      block->cplex.setParam(IloCplex::Param::RootAlgorithm, IloCplex::Primal);
      try {
         handles[b] = block->cplex.solve(true);
      } catch (...) {
         // If there is an exception then we need to kill and join
         // all remaining solves. Otherwise we may leak handles.
         while (--b > 0) {
            handles[b].kill();
            handles[b].join();
         }
         throw;
      }
   }

   // Wait for the various LP solves to complete.
   for (BlockVector::size_type b = 0; b < blocks.size(); ++b)
      handles[b].join();

   // See if we need to generate cuts.
   for (BlockVector::size_type b = 0; b < blocks.size(); ++b) {
      Block *const block = blocks[b];
      cutVal.clear();
      cutVar.clear();
      double cutlb = -IloInfinity;
      double cutub = IloInfinity;

      // We ust STL types here since they are exception safe.
      std::vector<double> tmp(master->vars.getSize(), 0);
      std::map<IloNumVar,double,ExtractableLess<IloNumVar> > rayMap;

      // Depending on the status either seperate a feasibility or an
      // optimality cut.
      switch (block->cplex.getStatus()) {
      case IloAlgorithm::Unbounded:
         {
            // The subproblem is unbounded. We need to extract a feasibility
            // cut from an unbounded ray of the problem (see also the comments
            // at the top of this file).
            std::cout << "unbounded ";
            block->cplex.getRay(rayVals, rayVars);
            cutub = 0.0;
            for (IloInt j = 0; j < rayVars.getSize(); ++j) {
               cutub -= rayVals[j] * block->objMap[rayVars[j]];
               rayMap[rayVars[j]] = rayVals[j];
            }
            for (std::vector<Block::FixData>::const_iterator it = block->fixed.begin(); it != block->fixed.end(); ++it)
               tmp[it->col] -= it->val * rayMap[block->vars[it->row]];
            for (IloInt j = 0; j < master->vars.getSize(); ++j) {
               if ( fabs (tmp[j]) > 1e-6 ) {
                  cutVar.add(master->vars[j]);
                  cutVal.add(tmp[j]);
//.........这里部分代码省略.........

// This routine separates Benders' cuts violated by the current x solution.
// Violated cuts are found by solving the worker LP
//
IloBool
separate(const Arcs x, const IloNumArray2 xSol, IloCplex cplex,
         const IloNumVarArray v, const IloNumVarArray u,
         IloObjective obj, IloExpr cutLhs, IloNum& cutRhs)
{
   IloBool violatedCutFound = IloFalse;

   IloEnv env = cplex.getEnv();
   IloModel mod = cplex.getModel();

   IloInt numNodes = xSol.getSize();
   IloInt numArcs  = numNodes * numNodes;
   IloInt i, j, k, h;

   // Update the objective function in the worker LP:
   // minimize sum(k in V0) sum((i,j) in A) x(i,j) * v(k,i,j)
   //          - sum(k in V0) u(k,0) + sum(k in V0) u(k,k)
   
   mod.remove(obj);
   IloExpr objExpr = obj.getExpr();
   objExpr.clear();
   for (k = 1; k < numNodes; ++k) {
      for (i = 0; i < numNodes; ++i) {
         for (j = 0; j < numNodes; ++j) {
               objExpr +=  xSol[i][j] * v[(k-1)*numArcs + i*numNodes + j];
         }
      }
   }
   for (k = 1; k < numNodes; ++k) {
      objExpr += u[(k-1)*numNodes + k];
      objExpr -= u[(k-1)*numNodes];
   }
   obj.setExpr(objExpr);
   mod.add(obj);
   objExpr.end(); 

   // Solve the worker LP

   cplex.solve();

   // A violated cut is available iff the solution status is Unbounded

   if ( cplex.getStatus() == IloAlgorithm::Unbounded ) {

      IloInt vNumVars = (numNodes-1) * numArcs;
      IloNumVarArray var(env);
      IloNumArray val(env);

      // Get the violated cut as an unbounded ray of the worker LP

      cplex.getRay(val, var);

      // Compute the cut from the unbounded ray. The cut is:
      // sum((i,j) in A) (sum(k in V0) v(k,i,j)) * x(i,j) >=
      // sum(k in V0) u(k,0) - u(k,k)

      cutLhs.clear();
      cutRhs = 0.;

      for (h = 0; h < val.getSize(); ++h) {

         IloInt *index_p = (IloInt*) var[h].getObject();
         IloInt index = *index_p;

         if ( index >= vNumVars ) {
            index -= vNumVars;
            k = index / numNodes + 1;
            i = index - (k-1)*numNodes;
            if ( i == 0 )
               cutRhs += val[h];
            else if ( i == k )
               cutRhs -= val[h];
         }
         else {
            k = index / numArcs + 1;
            i = (index - (k-1)*numArcs) / numNodes;
            j = index - (k-1)*numArcs - i*numNodes;
            cutLhs += val[h] * x[i][j];
         }
      }

      var.end();
      val.end();

      violatedCutFound = IloTrue;
   }

   return violatedCutFound;

} // END separate

//.........这里部分代码省略.........
		for(IloInt b=0;b<Node-1;++b)
		{
			detaP_node[b]=IloNumExpr(env);
			detaPw_node[b]=IloNumExpr(env);
			IloInt i=0;
			for(;i<NG;++i)
			{
				if(Unit[i][0]==b-1)break;
			}
			
			if(i<NG)
			{
				detaP_node[b]+=detaP[i];
			}	
			
			if(Sw[b]>=0)
			{
				detaPw_node[b]+=detaPw[ Sw[b] ];
			}
		}
		
		for(IloInt b=0;b<Node-1;++b)
		{
			Theta[b]=IloNumExpr(env);
			for(IloInt k=0;k<Node-1;++k)
			{			
				Theta[b]+=B0l[b][k]*(detaP_node[k]+detaPw_node[k]);	
			}
		}
		Theta[Node-1]=IloNumExpr(env);
		
		for(IloInt h=0;h<Branch;++h)
		{
			IloNumExpr exprTheta(env);//莫明其妙的错误
			exprTheta+=(Theta[(IloInt)Info_Branch[h][0]-1]-Theta[(IloInt)Info_Branch[h][1]-1]);
			
			Master_Model.add(exprTheta<=Info_Branch[h][3]*(Info_Branch[h][4]-PL[h]));			
			Master_Model.add(exprTheta>=Info_Branch[h][3]*(-Info_Branch[h][4]-PL[h]));
			exprTheta.end();
			//两个相减的节点顺序没有影响么？
		}
		Theta.end();
		detaP_node.end();
		detaPw_node.end();
		
		Master_Cplex.extract(Master_Model);
		Master_Cplex.solve();
		if (Master_Cplex.getStatus() == IloAlgorithm::Infeasible)//输出结果
		{
			output<<"Master Problem Have No Solution"<<endl;
			goto lable2;
		}
		
/************************************************************输出显示过程**************************************************/
		
		// output/*<<endl<<"Min:"*/<<Master_Cplex.getObjValue()<<endl;
		
		Master_Model.remove(min);
		Master_Model.add(max);
		Master_Cplex.extract(Master_Model);
		Master_Cplex.solve();
		if (Master_Cplex.getStatus() == IloAlgorithm::Infeasible)//输出结果
		{
			output<<"Master Problem Have No Solution"<<endl;
			goto lable2;
		}
		output/*<<endl<<"Max:"*/<<Master_Cplex.getObjValue()<<endl<<endl;;
		
		// output<<endl<<"常规机组出力调整量:"<<endl;
		// for(IloInt i=0;i<NG;++i)
		// {
			// output<<Master_Cplex.getValue(detaP[i])<<"  ";
		// }
		// output<<endl<<"风电机组出力调整量:"<<endl;
		// for(IloInt i=0;i<NW;++i)
		// {
			// output<<Master_Cplex.getValue(detaPw[i])<<"  ";
		// }
		// output<<endl;
		
		lable2:
		Master_Model.end();
		Master_Cplex.end();
		env.end();
	}
	catch(IloException& ex)//异常捕获
	{
		output<<"Error: "<<ex<<endl;
	}
	catch(...)
	{
		output<<"Error: Unknown exception caught!" << endl;
	}
	
	finish=clock();
	totaltime=(double)(finish-start)/CLOCKS_PER_SEC;
	output<<"totaltime: "<<totaltime<<"s"<<endl<<endl;
	output.close();	
	return 0;
}