![]() | |||||||||||
| |||||||||||
Cutstk - Column generation for a cutting stock problem Description This example features iteratively adding new variables,
basis in/output and working with subproblems. The column
generation algorithm is implemented as a loop over the
root node of the MIP problem.
Source Files By clicking on a file name, a preview is opened at the bottom of this page.
xbcutstk.c /******************************************************** BCL Example Problems ==================== file xbcutstk.c ``````````````` Cutting stock problem, solved by column (= cutting pattern) generation heuristic looping over the root node. (c) 2008-2023 Fair Isaac Corporation author: S.Heipcke, 2001, rev. Mar. 2014 ********************************************************/ #include <stdio.h> #include <stdlib.h> #include <math.h> #include "xprb.h" #include "xprs.h" #define NWIDTHS 5 #define MAXWIDTH 94 #define EPS 1e-6 #define MAXCOL 10 /****DATA****/ double WIDTH[] = {17, 21, 22.5, 24, 29.5}; /* Possible widths */ int DEMAND[] = {150, 96, 48, 108, 227}; /* Demand per width */ int PATTERNS[NWIDTHS][NWIDTHS]; /* (Basic) cutting patterns */ XPRBvar use[NWIDTHS+MAXCOL]; /* Rolls per pattern */ XPRBctr dem[NWIDTHS]; /* Demand constraints */ XPRBctr cobj; /* Objective function */ double knapsack(int N, double *c, double *a, double R, int *d, int *xbest); /***********************************************************************/ void modcutstock(XPRBprob prob) { int i,j; for(j=0;j<NWIDTHS;j++) PATTERNS[j][j]=(int)floor(MAXWIDTH/WIDTH[j]); /****VARIABLES****/ for(j=0;j<NWIDTHS;j++) use[j]=XPRBnewvar(prob,XPRB_UI, XPRBnewname("pat_%d",j+1),0, (int)ceil((double)DEMAND[j]/PATTERNS[j][j])); /****OBJECTIVE****/ cobj = XPRBnewctr(prob,"OBJ",XPRB_N); /* Minimize total number of rolls */ for(j=0;j<NWIDTHS;j++) XPRBaddterm(cobj, use[j], 1); XPRBsetobj(prob,cobj); /****CONSTRAINTS****/ for(i=0;i<NWIDTHS;i++) { dem[i] = XPRBnewctr(prob,"Demand",XPRB_G); /* Satisfy the demand per width */ for(j=0;j<NWIDTHS;j++) XPRBaddterm(dem[i], use[j], PATTERNS[i][j]); XPRBaddterm(dem[i], NULL, DEMAND[i]); } } /**************************************************************************/ /* Column generation loop at the top node: */ /* solve the LP and save the basis */ /* get the solution values */ /* generate new column(s) (=cutting pattern) */ /* load the modified problem and load the saved basis */ /**************************************************************************/ void solvecutstock(XPRBprob prob) { double objval; /* Objective value */ int i,j; int starttime; int npatt, npass; /* Counters for columns and passes */ double soluse[NWIDTHS+MAXCOL]; /* Solution values for variables pat */ double dualdem[NWIDTHS]; /* Dual values of demand constraints */ XPRBbasis basis; double dw,z; int x[NWIDTHS]; starttime=XPRBgettime(); npatt = NWIDTHS; for(npass=0;npass<MAXCOL;npass++) { XPRBlpoptimize(prob,""); /* Solve the LP */ basis=XPRBsavebasis(prob); /* Save the current basis */ objval = XPRBgetobjval(prob); /* Get the objective value */ /* Get the solution values: */ for(j=0;j<npatt;j++) soluse[j]=XPRBgetsol(use[j]); for(i=0;i<NWIDTHS;i++) dualdem[i]=XPRBgetdual(dem[i]); /* for(j=0;j<npatt;j++) printf("%g, ",soluse[j]); printf("\n"); */ /* Solve integer knapsack problem z = min{cx : ax<=r, x in Z^n} with r=MAXWIDTH, n=NWIDTHS */ z = knapsack(NWIDTHS, dualdem, WIDTH, (double)MAXWIDTH, DEMAND, x); printf("(%g sec) Pass %d: ",(XPRBgettime()-starttime)/1000.0, npass+1); if(z < 1+EPS) { printf("no profitable column found.\n\n"); XPRBdelbasis(basis); /* No need to keep the basis any longer */ break; } else { /* Print the new pattern: */ printf("new pattern found with marginal cost %g\n ", z-1); printf("Widths distribution: "); dw=0; for(i=0;i<NWIDTHS;i++) { printf ("%g:%d ", WIDTH[i], x[i]); dw += WIDTH[i]*x[i]; } printf("Total width: %g\n", dw); /* Create a new variable for this pattern: */ use[npatt]=XPRBnewvar(prob,XPRB_UI, XPRBnewname("pat_%d",npatt+1),0,XPRB_INFINITY); XPRBaddterm(cobj, use[npatt], 1); /* Add new var. to the objective */ dw=0; for(i=0; i<NWIDTHS; i++) /* Add new var. to demand constraints*/ if(x[i] > EPS) { XPRBaddterm(dem[i], use[npatt], x[i]); if((int)ceil((double)DEMAND[i]/x[i]) > dw) dw = (int)ceil((double)DEMAND[i]/x[i]); } XPRBsetub(use[npatt],dw); /* Change the upper bound on the new var.*/ npatt++; XPRBloadmat(prob); /* Reload the problem */ XPRBloadbasis(basis); /* Load the saved basis */ XPRBdelbasis(basis); /* No need to keep the basis any longer */ } } XPRBmipoptimize(prob,"c"); /* Solve the MIP */ printf("(%g sec) Optimal solution: %g rolls, %d patterns\n ", (XPRBgettime()-starttime)/1000.0, XPRBgetobjval(prob), npatt); printf("Rolls per pattern: "); for(i=0;i<npatt;i++) printf("%g, ", XPRBgetsol(use[i])); printf("\n"); } /**************************************************************************/ /* Integer Knapsack Algorithm for solving the integer knapsack problem */ /* z = max{cx : ax <= R, x <= d, x in Z^N} */ /* */ /* Input data: */ /* N: Number of item types */ /* c[i]: Unit profit of item type i, i=1..n */ /* a[i]: Unit resource use of item type i, i=1..n */ /* R: Total resource available */ /* d[i]: Demand for item type i, i=1..n */ /* Return values: */ /* xbest[i]: Number of items of type i used in optimal solution */ /* zbest: Value of optimal solution */ /**************************************************************************/ double knapsack(int N, double *c, double *a, double R, int *d, int *xbest) { int j; double zbest = 0.0; XPRBarrvar x; XPRBctr cobj; XPRBprob kprob; /* printf ("Solving z = max{cx : ax <= b; x in Z^n}\n c ="); for (j = 0; j < N; j++) printf (" %g,", c[j]); printf("\n a ="); for (j = 0; j < N; j++) printf (" %g,", a[j]); printf("\n c/a ="); for (j = 0; j < N; j++) printf (" %g,", c[j]/a[j]); printf("\n b = %f\n", R); */ kprob=XPRBnewprob("Knapsack"); x=XPRBnewarrvar(kprob,N, XPRB_UI, "x", 0, XPRB_INFINITY); for(j=0;j<N;j++) XPRBsetub(x[j], d[j]); cobj = XPRBnewarrsum(kprob,"OBJ", x, c, XPRB_N, 0); XPRBsetobj(kprob,cobj); XPRBnewarrsum(kprob,"KnapCtr", x, a, XPRB_L, R); XPRBsetsense(kprob,XPRB_MAXIM); XPRBmipoptimize(kprob,""); zbest = XPRBgetobjval(kprob); for(j=0;j<N;j++) xbest[j]=(int)floor(XPRBgetsol(x[j])+0.5); /* printf (" z = %f\n x =", zbest); for(j=0; j<N; j++) printf (" %d,", xbest[j]); printf("\n"); */ XPRBdelprob(kprob); return (zbest); } /***********************************************************************/ int main(int argc, char **argv) { XPRBprob prob; prob=XPRBnewprob("CutStock"); /* Initialize a new problem in BCL */ modcutstock(prob); /* Model the problem */ solvecutstock(prob); /* Solve the problem */ XPRBdelprob(prob); /* Delete the main problem */ XPRBfree(); /* Terminate BCL */ return 0; }
| |||||||||||
© Copyright 2023 Fair Isaac Corporation. |