[Symphony] Problem adding custom cutting planes
Ashutosh Mahajan
asm4 at lehigh.edu
Tue Mar 3 15:32:27 EST 2009
hi richard,
i can confirm that this is a bug in symphony. the way of packing cuts was
changed some time back but the appropriate changes were not made to
cg_add_explicit_cut() function. this bug is present in version
5.2 as well.
i think the quickest way out for you may be to just edit the source file
SYMPHONY/src/CutGen/cg_func.c yourself. you will have to compile symphony
again (and also your application code). here are the differences.
Index: ../../src/CutGen/cg_func.c
===================================================================
--- ../../src/CutGen/cg_func.c (revision 1588)
+++ ../../src/CutGen/cg_func.c (working copy)
@@ -188,11 +188,12 @@
cut->sense = sense;
cut->rhs = rhs;
cut->range = range;
- cut->size = ISIZE + nzcnt * (ISIZE + DSIZE);
+ cut->size = DSIZE + nzcnt * (ISIZE + DSIZE);
cut->coef = (char *) malloc (cut->size);
+ ((double *) cut->coef)[0] = 0; /* important to zero out */
((int *) cut->coef)[0] = nzcnt;
- memcpy(cut->coef + ISIZE, (char *)indices, nzcnt*ISIZE);
- memcpy(cut->coef + (nzcnt + 1) * ISIZE, (char *)values, nzcnt * DSIZE);
+ memcpy(cut->coef + DSIZE, (char *)values, nzcnt*DSIZE);
+ memcpy(cut->coef + (nzcnt + 1) * DSIZE, (char *)indices, nzcnt * ISIZE);
cut->branch = DO_NOT_BRANCH_ON_THIS_ROW;
cut->deletable = TRUE;
cut->name = send_to_cp ? CUT__SEND_TO_CP : CUT__DO_NOT_SEND_TO_CP;
let me know if you run into further problems. Thanks for catching this bug!
hopefully we will fix it soon.
ashutosh
On Tue, 03 Mar 2009, Richard Russell wrote:
> Hi Ashutosh,
> I'm using SYMPHONY 5.1.10. Thank you for having a look at this.
> Regards,
> Richard
>
> Ashutosh Mahajan wrote:
>> 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
>>
>> _______________________________________________
>> 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