#include <stdio.h>
#include <stdlib.h>
#include <glpk.h>
#include <float.h>
#include <string.h>
#include <assert.h>
#define MAX( a, b ) ( a > b ) ? ( a ) : ( b )
#define WEIGHT_SLACK_VARIABLE 1000
/** Input Data */
int L = 250; // bar size
int m = 4; // nr. of orders
int w[] = { 187, 119, 74, 90 }; // size of bar in order i
int b[] = { 1, 2, 2, 1 }; // number of bars in order i
/* vector of dual variables */
double *pi = NULL;
/* GLPK linear program */
glp_prob *lp = NULL;
/** solves the knapsack problem where repetition of items are allowed
cap: capacity n: number o items
weight[0,...n-1]: vector of weights
profit[0,...n-1]: vector of profits
selected[0,...n-1]: solution will be filled here, where
selected[i] contains an integer indicating how many times item i was included
work[] is a vector of size (cap+1)*2 where dynamic programming computation will occur **/
int solve_knapsack( int cap, int n, int weight[], int profit[],
int selected[], int work[] );
/** creates constraints and adds
artificial variables ***/
void fill_initial_LP();
int main( int argc, char **argv )
{
int exit_code = EXIT_SUCCESS;
lp = glp_create_prob();
pi = malloc( sizeof(double)*(m+1) );
int iteration = 0;
/// dynamic programming related memory structures
int *dpWork = NULL;
dpWork = malloc( sizeof(int)*((L+1)*2) );
int *selected = NULL;
selected = malloc( sizeof(int)*m );
int *ipi = NULL;
ipi = malloc( sizeof(int)*m );
int nColsLambda = 0;
int *colIdx = NULL;
double *colCoef = NULL;
colIdx = malloc( sizeof(int)*(m+1) );
colCoef = malloc( sizeof(double)*(m+1) );
printf("Cutting Stock Solver.\n\n");
fill_initial_LP();
int newCols;
SOLVE_RELAXATION:
++iteration;
printf("-> on iteration %d\n\n", iteration );
glp_simplex( lp, NULL );
if ( glp_get_status( lp ) != GLP_OPT )
{
fprintf( stderr, "ERROR: No optimal solution obtained, check problem data.\n" );
exit_code = EXIT_FAILURE;
goto TERMINATE;
}
/* getting dual values */
int i;
for ( i=1 ; (i<=m) ; ++i )
pi[i] = glp_get_row_dual( lp, i );
printf("\npi vector - dual variables for constraints: \n\t");
for ( i=1 ; (i<=m) ; ++i )
printf("[%d] %g ", i, pi[i] );
printf("\n\n");
/** solving the pricing problem ***/
for ( i=0 ; (i<m) ; ++i ) // knapsack solver needs
ipi[i] = ( pi[i+1] * 1000 ); // integer values
double rc = 1.0 - (((double)( solve_knapsack( L, m, w, ipi, selected, dpWork ) )) / 1000.0);
newCols = 0;
if ( rc < -0.000001 )
{
newCols++;
printf("\tattractive column rc:(%g)\n", rc );
printf("\t\tselected: ");
int k;
for ( k=0 ; (k<m) ; ++k )
if (selected[k])
printf("%d(%d) ", w[k], selected[k] );
printf("\n");
int nSizes = 0;
for ( k=0 ; (k<m) ; ++k )
if (selected[k])
{
nSizes++;
colIdx[ nSizes ] = k+1;
colCoef[ nSizes ] = selected[k];
}
/// inserting column
glp_add_cols( lp, 1 );
int col = glp_get_num_cols( lp );
nColsLambda++;
char colName[ 256 ];
sprintf( colName, "lbda[%d]", nColsLambda );
glp_set_col_name( lp, col, colName );
glp_set_col_bnds( lp, col, GLP_DB, 0.0, 1.0 );
glp_set_obj_coef( lp, col, 1 );
glp_set_mat_col( lp, col, nSizes, colIdx, colCoef );
}
char pNameIt[256];
sprintf( pNameIt, "cg_%d.lp", iteration );
glp_write_lp( lp, NULL, pNameIt );
if (newCols)
{
printf("\nPress Enter to proceed to the next iteration\n");
getchar();
printf("New columns added. Starting a new iteration.\n\n");
goto SOLVE_RELAXATION;
}
/// solution
double objVal = glp_get_obj_val( lp );
printf("solution found with cost %g\n", objVal );
int nPattern = 0, j;
for ( j=1 ; (j<=glp_get_num_cols(lp)) ; ++j )
if ( ( glp_get_col_prim( lp, j ) >= 0.01 ) && ( strstr( glp_get_col_name( lp, j ), "lbda" ) ) )
{
printf(" %s %g \n", glp_get_col_name(lp, j), glp_get_col_prim( lp, j ) );
int nzCol = glp_get_mat_col( lp, j, colIdx-1, colCoef-1 );
int k;
printf("\t");
for ( k=0 ; k<nzCol ; ++k )
printf( "%d(%d) ", w[colIdx[k]-1], (int)colCoef[k] );
printf("\n");
nPattern++;
}
printf("\n");
glp_free( lp );
glp_free_env();
TERMINATE:
if (pi)
free( pi );
if (dpWork)
free(dpWork);
if (selected)
free( selected );
if (ipi)
free( ipi );
if (colIdx)
free(colIdx);
if (colCoef)
free(colCoef);
return exit_code;
}
void fill_initial_LP()
{
glp_set_prob_name( lp, "cutting stock" );
int i;
/// artificial variables, one for each
/// constraint which will be created after
glp_add_cols( lp, m );
for ( i=1 ; (i<=m) ; ++i )
{
char name[64];
sprintf( name, "art_%d", i );
glp_set_col_name( lp, i, name );
glp_set_col_bnds( lp, i, GLP_DB, 0.0, 1.0 );
glp_set_obj_coef( lp, i, WEIGHT_SLACK_VARIABLE );
}
/// creating one constraint for each size
glp_add_rows( lp, m );
for ( i=1 ; (i<=m) ; ++i )
{
int idx[] = { 0, i };
double coef[] = { 0.0, b[i-1] };
glp_set_mat_row( lp, i, 1, idx, coef );
glp_set_row_bnds( lp, i, GLP_LO, b[i-1], DBL_MAX );
}
}
int solve_knapsack( int cap, int n, int weight[], int profit[],
int selected[], int work[] )
{
int bestProfit = -1;
/// best profit within a given capacity
int *bestProfitWC = work;
int *bestItemToIncludeWC = work + (cap+1);
memset( bestProfitWC, 0, sizeof(int)*(cap+1) );
memset( selected, 0, sizeof(int)*n );
int i;
for ( i=0 ; (i<=cap) ; ++i )
bestItemToIncludeWC[i] = -1;
int lastUsefullCapacity = -1;
int curCap;
for ( curCap=1 ; (curCap<=cap) ; ++curCap )
{
for ( i=0 ; (i<n) ; ++i )
{
int profitWith = -1;
int rCap = curCap - weight[ i ];
if ( rCap >= 0 )
profitWith = profit[i] + bestProfitWC[rCap];
if ( profitWith > bestProfitWC[curCap] )
{
bestProfitWC[curCap] = profitWith;
bestItemToIncludeWC[curCap] = i;
lastUsefullCapacity = curCap;
}
}
}
if ( lastUsefullCapacity == -1 )
return 0;
bestProfit = bestProfitWC[lastUsefullCapacity];
curCap = lastUsefullCapacity;
while ( curCap >0 )
{
int curItem = bestItemToIncludeWC[curCap];
if ( curItem == -1 )
return bestProfit;
selected[ curItem ]++;
curCap -= weight[ curItem ];
}
return bestProfit;
}