[Symphony] Problem adding custom cutting planes

[Richard Russell]

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