[Symphony] Problem adding custom cutting planes
Ashutosh Mahajan
asm4 at lehigh.edu
Tue Mar 3 10:14:39 EST 2009
hi richard,
what version of symphony are you using? i will take a look to see if symphony
is removing or not accepting the cuts for some reason.
regards
ashutosh
On Tue, 03 Mar 2009, Richard Russell wrote:
> Hi,
>
> I'm having problems generating custom cutting planes. I am identifying
> and attempting to add sensible cuts but they don't seem to be respected
> by the LP solver, i.e. they are violated in subsequent LP solutions.
> This has the result that the program just loops adding the same cut.
>
> The problem I'm working on has a collection of constraints that are
> represented explicitly in the constraint matrix. However there are some
> constraints that are too numerous to exhaustively list in the constraint
> matrix so I would like to add them dynamically as cuts.
>
> I am working with the USER function stubs in SYMPHONY/Applications and I
> have modified:
> (1) user_find_cuts in the file 'user_cg.c'
> (2) user_is_feasible in the file 'user_lp.c'
>
> A description of what each function is doing:
>
> (1) user_find_cuts first checks that the solution passed to it is
> integer. If it isn't it lets SYMPHONY return it's default behaviour
> without generating any cuts. If it is integer it then look for cycles in
> a graph constructed from the solution. The graph has a node for each
> variable that is set to 1.0 in the solution. These nodes are related in
> the problem and have directed arcs between them. If a cycle exists in
> the graph constructed from the solution then a cut needs to be added to
> prevent this cycle. This is done for each cycle by adding a cut of the form:
>
> var_1 + ... + var_n <= (n-1) for every cycle in the graph where var_1
> to var_n are the nodes in the cycle.
>
> So at least one of var_1 to var_n must be set to 0.0 in order for the
> constraint to hold. This way the cycle is prevented.
>
> (2) user_is_feasible checks that the solution is integer. If it isn't it
> returns infeasible. If it is integer it checks that there are no cycles
> in the graph constructed from the solution. If a cycle is found it
> returns infeasible. If no cycles are found it returns feasible.
>
> Here is most of the code, I have omitted some tedious initialisations of
> arrays and such.
>
> ///////////////////////////////////////////////////////////////////
> int user_find_cuts(void *user, int varnum, int iter_num, int level,
> int index, double objval, int *indices, double *values,
> double ub, double etol, int *num_cuts, int *alloc_cuts,
> cut_data ***cuts)
> {
> user_problem *prob = (user_problem *) user;
>
> int a,k,t,i;
> int horizon;
> int result;
> SASPtr sasptr;
> IPGraphPtr ipgraphptr;
> List<IPNodePtr>** cyclelist;
> int *nrcycles;
>
> // Cut related variables
> int nzcnt;
> int *cut_indices;
> double *cut_values;
> double rhs;
> bool integer = true;
>
> // Debug information
> printf("Entering cut generating routine.\n");
>
> printf("Node details:\n");
> printf("Index = %d\n", index);
> printf("Level = %d\n", level);
> printf("Iter = %d\n", iter_num);
>
> // Initialise cycle data structures (Code omitted)
>
> sasptr = prob->solverPtr->GetSASPtr();
> horizon = prob->solverPtr->GetHorizon();
>
> // Initialise vars[a][k][t] array (code omitted)
>
> // Construct full 0-1 solution from non-zero encoding
> for(t=0; t<horizon; t++) {
> for(a=0; a<sasptr->GetNrActions(); a++) {
> for(k=0; k<sasptr->GetNrActionProfiles(a); k++) {
> vars[a][k][t] = 0.0;
> for (i=0; i<varnum; i++) {
> if (indices[i] == prob->solverPtr->mX[a][k][t].getID()) {
> printf("%d = %e\n", indices[i], values[i]);
> vars[a][k][t] = values[i];
>
> if ( fabs(values[i]-1.0) > 0.001 ) {
> integer = false;
> }
> break;
> }
> }
> }
> }
> }
>
> // If not an integer solution then don't generate cuts
> if (!integer) {
> printf("LP solution not integer, do not generate cuts.\n");
> return (USER_DEFAULT);
> }
>
> // Count number of cycles for each time step
> printf("Checking for cycles in cut generation...\n");
> for(t=0; t<horizon; t++) {
> ipgraphptr = prob->solverPtr->GetIPGraph(t, vars);
>
> // Returns number of cycles and stores cycles in cyclelist[t] parameter
> nrcycles[t] = ipgraphptr->GetIntegralCycles(cyclelist[t], MAXCYCLES);
>
> printf("Time step = %d, Number of cycles = %d\n", t, nrcycles[t]);
> }
>
> int cut_count = 0;
>
> for(t=0; t<horizon; t++) {
> for(i=0; (i<MAXCYCLES) && (cyclelist[t][i].GetSize()>0); i++) {
> nzcnt = 0;
>
> cut_indices = new int[cyclelist[t][i].GetSize()];
> cut_values = new double[cyclelist[t][i].GetSize()];
>
> printf("Adding following cut:\n");
>
> ListIterator<IPNode*> iter(cyclelist[t][i]);
> for(iter.GotoBegin(); !iter.AtEnd(); iter.GotoNext()) {
> a = IPNodePtr(iter.GetElement())->GetA();
> k = IPNodePtr(iter.GetElement())->GetK();
> t = IPNodePtr(iter.GetElement())->GetT();
>
> cut_indices[nzcnt] = prob->solverPtr->mX[a][k][t].getID();
> cut_values[nzcnt] = 1.0;
>
> // Output debug information: the cut being generated
> printf("%e*%d + ",cut_values[nzcnt], cut_indices[nzcnt]);
>
> nzcnt++;
> }
>
> rhs = ((double)nzcnt) - 1.0;
> // Output debug information: rhs
> printf("<= %e\n", rhs);
>
> result = cg_add_explicit_cut(nzcnt, cut_indices, cut_values, rhs,
> 0.0,
> 'L', TRUE, num_cuts, alloc_cuts, cuts);
>
> if ( result != 1 ) {
> printf("WARNING: Cuts not being generated.\n");
> }
> else {
> cut_count++;
> }
> }
> }
>
> if (cut_count > 0)
> return (USER_SUCCESS);
>
> return (USER_DEFAULT);
> }
>
> //////////////////////////////////////////////////////////////////////////
> int user_is_feasible(void *user, double lpetol, int varnum, int *indices,
> double *values, int *feasible, double *objval,
> char branching, double *heur_solution)
> {
> user_problem *prob = (user_problem *) user;
>
> int a,k,t,i,j;
> int size, sum;
> int horizon;
> SASPtr sasptr;
> IPGraphPtr ipgraphptr;
> List<IPNodePtr>** cyclelist;
> int *nrcycles;
> // Cut related variables
> bool integer = true;
>
> printf("Entering user_is_feasible function.\n");
>
> // Initialise cycle data structures (code omitted)
>
> sasptr = prob->solverPtr->GetSASPtr();
> horizon = prob->solverPtr->GetHorizon();
>
> // Initialise vars[a][k][t] array (Code omitted)
>
> // Construct full 0-1 solution from non-zero encoding
> for(t=0; t<horizon; t++) {
> for(a=0; a<sasptr->GetNrActions(); a++) {
> for(k=0; k<sasptr->GetNrActionProfiles(a); k++) {
> vars[a][k][t] = 0.0;
> for (i=0; i<varnum; i++) {
> if (indices[i] == prob->solverPtr->mX[a][k][t].getID()) {
> printf("%d = %e\n", indices[i], values[i]);
> vars[a][k][t] = 1.0;
>
> if ( fabs(values[i] - 1.0) > lpetol) {
> integer = false;
> }
> break;
> }
> }
> }
> }
> }
>
> // First check that the solution is integer
> if (!integer) {
> printf("INFEASIBLE: Not integer.\n");
> *feasible = IP_INFEASIBLE;
> return (USER_SUCCESS);
> }
>
> printf("Checking for cycles...\n");
>
> bool violated_cycles = false;
>
> // Count the number of cycles for each time step
> for(t=0; t<horizon; t++) {
> ipgraphptr = prob->solverPtr->GetIPGraph(t, vars);
>
> nrcycles[t] = ipgraphptr->GetIntegralCycles(cyclelist[t], MAXCYCLES);
>
> if (nrcycles[t] > 0) {
> violated_cycles = true;
> }
>
> // Output debug information
> printf("Time step = %d, Number of cycles = %d\n", t, nrcycles[t]);
> }
>
> // If cycles have been found output the corresponding violated constraint
> if( violated_cycles ) {
> printf("INFEASIBLE: Violated cycle constraint.\n");
> printf("Found the following violated cycles:\n");
>
> for(t=0;t<horizon;t++) {
> for (i=0; (i<MAXCYCLES) && (cyclelist[t][i].GetSize()>0); i++) {
> sum = 0;
> size = 0;
>
> ListIterator<IPNode*> iter(cyclelist[t][i]);
> for (iter.GotoBegin(); !iter.AtEnd(); iter.GotoNext()) {
> a = IPNodePtr(iter.GetElement())->GetA();
> k = IPNodePtr(iter.GetElement())->GetK();
> t = IPNodePtr(iter.GetElement())->GetT();
>
> printf("%e*%d +", vars[a][k][t],
> prob->solverPtr->mX[a][k][t].getID());
>
> size++;
> }
>
> printf("<= %d\n", size-1);
> }
> }
>
> // Return infeasible solution
> *feasible = IP_INFEASIBLE;
> return (USER_SUCCESS);
> }
> else {
> // Integer and no cycles so solution is feasible
> *feasible = IP_FEASIBLE;
> return (USER_SUCCESS);
> }
> }
>
> Running the above code, here is a sample of the output I get:
>
> Entering user_is_feasible function.
> 74 = 1.000000e+00
> 116 = 1.000000e+00
> 99 = 1.000000e+00
> 131 = 1.000000e+00
> Checking for cycles...
> Time step = 0, Number of cycles = 1
> Time step = 1, Number of cycles = 0
> INFEASIBLE: Violated cycle constraint.
> Found the following violated cycles:
> 1.000000e+00*116 +1.000000e+00*74 +<= 1
> Entering cut generating routine.
> Node details:
> Index = 2
> Level = 1
> Iter = 142
> 74 = 1.000000e+00
> 116 = 1.000000e+00
> 99 = 1.000000e+00
> 131 = 1.000000e+00
> Checking for cycles in cut generation...
> Time step = 0, Number of cycles = 1
> Time step = 1, Number of cycles = 0
> Adding following cut:
> 1.000000e+00*116 + 1.000000e+00*74 + <= 1.000000e+00
> Entering user_is_feasible function.
> 74 = 1.000000e+00
> 116 = 1.000000e+00
> 99 = 1.000000e+00
> 131 = 1.000000e+00
> Checking for cycles...
> Time step = 0, Number of cycles = 1
> Time step = 1, Number of cycles = 0
> INFEASIBLE: Violated cycle constraint.
> Found the following violated cycles:
> 1.000000e+00*116 +1.000000e+00*74 +<= 1
>
> This just continues repeating with node index and level remaining
> unchanged and the iteration number increasing.
> Adding the cut should prohibit variable #74 and #116 being
> simultaneously set to 1.0 but the cut doesn't seem to have any effect.
>
> Does anybody notice anything I haven't done that I should or have any
> suggestions for things I could try?
>
> I would be very grateful for any help anyone could offer.
>
> Kind regards,
> Richard
>
> _______________________________________________
> Symphony mailing list
> Symphony at list.coin-or.org
> http://list.coin-or.org/mailman/listinfo/symphony
--
regards
Ashutosh Mahajan
http://coral.ie.lehigh.edu/~asm4
More information about the Symphony
mailing list