Big M

Introduction

The example shows how to implement a big M formulation using the GLPK C library.

GMPL model

set S;
param a{S};
param b{S};
param M;

var w{S}, >= 0, <=1.5;
var x{S}, binary;

maximize obj :
  sum{i in S} a[i] * w[i];

s.t. c1{i in S}:
  w[i] <= x[i] * M;

s.t. c2 :
  sum{i in S} b[i] * w[i] <= 1;

s.t. c3 :
  sum{i in S} x[i] <= 5;

display S;

solve;

printf "   %6s %6s %6s %1s\n", "a", "b", "w", "x";
for {i in S}
  printf "%1d: %6.3f %6.3f %6.3f %1.0f\n", i, a[i], b[i], w[i], x[i];

data;

set S := 0 1 2 3 4 5 6 7 8 9;

param a :=
  0 0.44
  1 0.47
  2 0.11
  3 0.92
  4 0.95
  5 0.56
  6 0.45
  7 0.95
  8 0.68
  9 0.59;

param b :=
  0 0.49
  1 0.96
  2 0.33
  3 0.09
  4 0.75
  5 0.72
  6 0.01
  7 0.46
  8 0.11
  9 0.85;

param M := 1;

end;

C program

#include <stdio.h>
#include "glpk.h"

int main(){
  const double  a[] = 
    {0.44, 0.47, 0.11, 0.92, 0.95, 0.56, 0.45, 0.95, 0.68, 0.59};
  const double  b[] =
    {0.49, 0.96, 0.33, 0.09, 0.75, 0.72, 0.01, 0.46, 0.11, 0.85};
  char          colName[256];  
  int           i;
  int           ind[11];
  glp_prob     *lp;
  const int     n = 10;
  glp_iocp      parm;
  int           ret;
  const int     S[] = 
    { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 };
  const double  M = 1.; // any value >= max(w[])
  char          rowName[256];  
  double        val[11];
  double        w[10];
  int           x[10];

  // create problem object
  lp = glp_create_prob();
  // set problem name
  glp_set_prob_name(lp, "BigMDemo");  
  // assign objective function name
  glp_set_obj_name(lp, "obj");
  // set optimization direction
  glp_set_obj_dir(lp, GLP_MAX);
  // add columns to problem object
  glp_add_cols(lp, 2 * n);
  for(i = 1; i<=n; i++){
    sprintf(colName, "w_%1d", i);
    // set column name
    glp_set_col_name(lp, i, colName);
    // set column bounds
    glp_set_col_bnds(lp, i, GLP_DB, 0., 1.);
    // set column kind
    glp_set_col_kind(lp, i+n, GLP_CV);
    // set objective coefficient
    glp_set_obj_coef(lp, i, a[i-1]);
    sprintf(colName, "x_%1d", i);
    // set column name
    glp_set_col_name(lp, i+n, colName);
    // set column kind
    glp_set_col_kind(lp, i+n, GLP_BV);
  }
  // add rows to problem object
  glp_add_rows(lp, 2+n);
  for(i = 1; i<=n; i++){
    sprintf(rowName, "c1_%1d", i);
    // set row name
    glp_set_row_name(lp, i, rowName);
    // set row bounds
    glp_set_row_bnds(lp, i, GLP_LO, 0., 0.);
    // w[i] <= x[i] * M;
    ind[1] = i;
    val[1] = -1.;
    ind[2] = i+n;
    val[2] = M;
    glp_set_mat_row(lp, i, 2, ind, val);
  }

  // set row name
  glp_set_row_name(lp, n+1, "c2");
  // set row bounds
  glp_set_row_bnds(lp, n+1, GLP_UP, 0., 1.);
  //sum{i in S} b[i] * w[i] <= 1;
  for(i = 1; i<=n; i++){
    ind[i] = i;
    val[i] = b[i-1];
  }
  glp_set_mat_row(lp, n+1, n, ind, val);

  // set row name
  glp_set_row_name(lp, n+2, "c3");
  // set row bounds
  glp_set_row_bnds(lp, n+2, GLP_UP, 0., 5.);
  // sum{i in S} x[i] <= 5;
  for(i = 1; i<=n; i++){
    ind[i] = i+n;
    val[i] = 1.;
  }
  glp_set_mat_row(lp, n+2, n, ind, val);

  // Uncomment the following line to output the problem in CPLEX format
  // glp_write_lp(lp, NULL, "BigMDemo.lp");

  // initialize integer optimizer parameters
  glp_init_iocp(&parm);
  parm.presolve = GLP_ON;
  // solve MIP
  glp_intopt(lp, &parm);  
  // determine status of MIP solution 
  ret = glp_mip_status(lp);
  switch (ret) {
    case GLP_OPT:
      printf("Success\n");
      break;
    default:
      printf("Failed\n");
      // delete problem
      glp_delete_prob(lp);
      return 1;
  }
  for(i=1; i<=n; i++){
    // get column value
    w[i-1] = glp_mip_col_val(lp, i);
    x[i-1] = glp_mip_col_val(lp, i+n);
  }

  // delete problem
  glp_delete_prob(lp);

  // output the solution
  printf("   %6s %6s %6s %1s\n", "a", "b", "w", "x");
  for(i = 0; i < n; i++){
    printf("%1d: %6.3f %6.3f %6.3f %1d\n",i, a[i], b[i], w[i], (int) x[i]);
  }
  return 0;
}