| |||||||||||
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.java /******************************************************** Xpress-BCL Java Example Problems ================================ file xbcutstk.java `````````````````` Cutting stock problem, solved by column (= cutting pattern) generation heuristic looping over the root node. (c) 2008-2024 Fair Isaac Corporation author: S.Heipcke, 2001, rev. Mar. 2014 ********************************************************/ import java.lang.*; import com.dashoptimization.*; public class xbcutstk { static final int NWIDTHS = 5; static final int MAXWIDTH = 94; static final double EPS = 1e-6; static final int MAXCOL = 10; /****DATA****/ static final double[] WIDTH = {17, 21, 22.5, 24, 29.5}; /* Possible widths */ static final int[] DEMAND = {150, 96, 48, 108, 227}; /* Demand per width */ static int[][] PATTERNS; /* (Basic) cutting patterns */ static XPRBvar[] pat; /* Rolls per pattern */ static XPRBctr[] dem; /* Demand constraints */ static XPRBctr cobj; /* Objective function */ /***********************************************************************/ static void modCutStock(XPRBprob p) throws XPRSexception { int i,j; XPRBexpr le; PATTERNS = new int[NWIDTHS][NWIDTHS]; for(j=0;j<NWIDTHS;j++) PATTERNS[j][j] = (int)Math.floor((double)MAXWIDTH/WIDTH[j]); /****VARIABLES****/ pat = new XPRBvar[NWIDTHS+MAXCOL]; for(j=0;j<NWIDTHS;j++) pat[j]=p.newVar("pat_"+(j+1), XPRB.UI, 0, (int)Math.ceil((double)DEMAND[j]/PATTERNS[j][j])); /****OBJECTIVE****/ le = new XPRBexpr(); for(j=0;j<NWIDTHS;j++) le.add(pat[j]); cobj = p.newCtr("OBJ", le); p.setObj(cobj); /* Minimize total number of rolls */ /****CONSTRAINTS****/ dem = new XPRBctr[NWIDTHS]; for(i=0;i<NWIDTHS;i++) { /* Satisfy the demand per width */ le = new XPRBexpr(); for(j=0;j<NWIDTHS;j++) le .add(pat[j].mul(PATTERNS[i][j])); dem[i] = p.newCtr("Demand", le.gEql(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 */ /**************************************************************************/ static void solveCutStock(XPRB bcl, XPRBprob p) throws XPRSexception { double objval=0; /* Objective value */ int i,j; int starttime; int npatt, npass; /* Counters for columns and passes */ double[] solpat; /* Solution values for variables pat */ double[] dualdem; /* Dual values of demand constraints */ XPRBbasis basis; double dw,z; int[] x; x = new int[NWIDTHS]; solpat = new double[NWIDTHS+MAXCOL]; dualdem = new double[NWIDTHS]; starttime=XPRB.getTime(); npatt = NWIDTHS; for(npass=0;npass<MAXCOL;npass++) { p.lpOptimize(""); /* Solve the LP */ basis = p.saveBasis(); /* Save the current basis */ objval = p.getObjVal(); /* Get the objective value */ /* Get the solution values: */ for(j=0;j<npatt;j++) solpat[j]=pat[j].getSol(); for(i=0;i<NWIDTHS;i++) dualdem[i]=dem[i].getDual(); /* Solve integer knapsack problem z = min{cx : ax<=r, x in Z^n} with r=MAXWIDTH, n=NWIDTHS */ z = knapsack(bcl, NWIDTHS, dualdem, WIDTH, (double)MAXWIDTH, DEMAND, x); /* Force garbage collection to clean up subproblem: */ /* System.gc(); System.runFinalization(); */ System.out.print("(" + (XPRB.getTime()-starttime)/1000.0 + " sec) Pass " + (npass+1) + ": "); if(z < 1+EPS) { System.out.println("no profitable column found."); System.out.println(); basis=null; /* No need to keep the basis any longer */ break; } else { /* Print the new pattern: */ System.out.println("new pattern found with marginal cost " + (z-1) ); System.out.print(" Widths distribution: "); dw=0; for(i=0;i<NWIDTHS;i++) { System.out.print(WIDTH[i] +":"+ x[i] +" "); dw += WIDTH[i]*x[i]; } System.out.println("Total width: " + dw); /* Create a new variable for this pattern: */ pat[npatt]=p.newVar("pat_"+(npatt+1),XPRB.UI); cobj .add(pat[npatt]); /* Add new var. to the objective */ dw=0; for(i=0; i<NWIDTHS; i++) /* Add new var. to demand constraints*/ if(x[i] > EPS) { dem[i].add(pat[npatt].mul(x[i])); if((int)Math.ceil((double)DEMAND[i]/x[i]) > dw) dw = (int)Math.ceil((double)DEMAND[i]/x[i]); } pat[npatt].setUB(dw); /* Change the upper bound on the new var.*/ npatt++; p.loadMat(); /* Reload the problem */ p.loadBasis(basis); /* Load the saved basis */ basis=null; /* No need to keep the basis any longer */ } } p.mipOptimize(""); /* Solve the MIP */ System.out.println("(" + (XPRB.getTime()-starttime)/1000.0 + " sec) Optimal solution: " + p.getObjVal() + " rolls, " + npatt + " patterns"); System.out.print(" Rolls per pattern: "); for(i=0;i<npatt;i++) System.out.print(pat[i].getSol() + ", "); System.out.println(); } /**************************************************************************/ /* 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 */ /**************************************************************************/ static double knapsack(XPRB bcl, int N, double[] c, double[] a, double R, int[] d, int[] xbest) { /* System.out.println("Solving z = max{cx : ax <= b; x in Z^n}"); System.out.print(" c ="); for(j = 0; j < N; j++) System.out.print(" " + c[j] + ","); System.out.println(); System.out.print(" a ="); for(j = 0; j < N; j++) System.out.print(" " + a[j] + ","); System.out.println(); System.out.print(" c/a ="); for (j = 0; j < N; j++) System.out.print(" " + c[j]/a[j] + ","); System.out.println(); System.out.print(" b = " + R); */ try (XPRBprob pk = bcl.newProb("Knapsack")) { XPRBvar[] x = new XPRBvar[N]; for(int j=0;j<N;j++) x[j]=pk.newVar("x", XPRB.UI, 0, d[j]); XPRBexpr klobj = new XPRBexpr(); for(int j=0;j<N;j++) klobj.add(x[j].mul(c[j])); pk.setObj(klobj); XPRBexpr knap = new XPRBexpr(); for(int j=0;j<N;j++) knap.add(x[j].mul(a[j])); pk.newCtr("KnapXPRBctr", knap.lEql(R) ); pk.setSense(XPRB.MAXIM); pk.mipOptimize(""); double zbest = pk.getObjVal(); for(int j=0;j<N;j++) xbest[j]=(int)Math.floor(x[j].getSol() + 0.5); /* System.out.println(" z = " + zbest); System.out.print(" x ="); for(j=0; j<N; j++) System.out.print(" " + xbest[j] + ","); System.out.println(); */ return (zbest); } } /***********************************************************************/ public static void main(String[] args) throws XPRSexception { try (XPRB bcl = new XPRB(); /* Initialize BCL */ XPRBprob p = bcl.newProb("Els"); /* Create a new problem */ XPRS xprs = new XPRS()) { /* Initialize Xpress-Optimizer */ modCutStock(p); /* Model the problem */ solveCutStock(bcl, p); /* Solve the problem */ } } } | |||||||||||
© Copyright 2023 Fair Isaac Corporation. |