Initializing help system before first use

Cutstk - Column generation for a cutting stock problem


Type: Cutting stock
Rating: 4 (medium-difficult)
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.
File(s): xbcutstk.java


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-2023 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 */
        }
    }
}

© 2001-2023 Fair Isaac Corporation. All rights reserved. This documentation is the property of Fair Isaac Corporation (“FICO”). Receipt or possession of this documentation does not convey rights to disclose, reproduce, make derivative works, use, or allow others to use it except solely for internal evaluation purposes to determine whether to purchase a license to the software described in this documentation, or as otherwise set forth in a written software license agreement between you and FICO (or a FICO affiliate). Use of this documentation and the software described in it must conform strictly to the foregoing permitted uses, and no other use is permitted.