/********************************************************
BCL Example Problems
====================
file xbcutstk.c
```````````````
Cutting stock problem, solved by column (= cutting
pattern) generation heuristic looping over the root
node.
(c) 2008 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;
}
|
