UFL::create_model()

Here, we explain each step in this procedure. First, we must initialize the OSI solver interface. The solver choice is determined by the -DCOIN_USE_$(SOLVER) in the Makefile.

  si = initialize_solver();

Variables

Next, we must define the columns (variables) in the problem. The functions xindex(i,j) and yindex(j) define the order of the columns. First we list the x indices i = 0, ..., m-1, j = 0, ..., n-1, then the y indices j = 0, ..., n-1. We will also store the indices of those columns which are integral (the y variables).

  int n_cols = (M * N) + N;
  double * objective    = new double[n_cols];//the objective coefficients
  double * col_lb       = new double[n_cols];//the column lower bounds
  double * col_ub       = new double[n_cols];//the column upper bounds
  int    * integer_vars = new int[N];        //the indices of the integral vars

Here we use some COIN helper functions which are defined in: Coin/include/CoinHelperFunctions.hpp. The function CoinIotaN() copies the range: mn, mn + 1, ..., mn + n - 1, into the array integer_vars. These represent the column indices that correspond to those variables that are integral (i.e., all of the y variables).

  CoinIotaN(integer_vars, N, M * N);

The function CoinFillN() makes n_cols copies of the value 0.0 into the array col_lb starting at the beginning. This represents the lower bounds for all variables x and y. That is, $x_{ij} \geq 0, \forall i \in M, j \in N$ and $y_j \geq 0, \forall j \in N$.

  CoinFillN(col_lb, n_cols, 0.0);

Next, we set the objective coefficient $c_{ij} / d_i$ and the upper bound $d_i$, for each $x_{ij}$ variable. That is, $x_{ij} \leq d_i, \forall i \in M, j \in N$.

  int i, j, index;
  for(i = 0; i < M; i++){
    for(j = 0; j < N; j++){
      index            = xindex(i,j);
      objective[index] = trans_cost[index] / demand[i];
      col_ub[index]    = demand[i];
    }
  }

Then, we set the objective coefficient $f_j$ and the upper bound 1.0 for each $y_j$ variable. That is, $y_j \leq 1.0, \forall j \in N$. The function CoinDisjointCopyN() copies the n members of the array fixed_cost to the array objective starting at index $mn$.

  CoinFillN(col_ub + (M*N), N, 1.0);
  CoinDisjointCopyN(fixed_cost, N, objective + (M * N));

Constraints

Next we must define the set of rows (constraints). In order to do this, we will store the information in a CoinPackedMatrix.

There are several ways to construct a CoinPackedMatrix. In this implementation, we will define the matrix as row-ordered. That means, that we will construct the matrix with rows of coefficients. Alternatively, you can construct the matrix column-ordered. The choice typically depends on the application.

For our model, there are m demand constraints

$$f_i(x,y) = \sum_{j \in N} x_{ij} = d_i, \forall i \in M$$ (1)

and n linking constraints

$$g_j(x,y) = \sum_{i \in M} x_{ij} - \sum_{i \in M} d_i y_j \leq 0, \forall j \in N .$$ (2)

In order to explicitly define the row bounds, we rewrite the demand constraints as $d_i \leq f_i(x,y) \leq d_i$ and the linking constraints $- \infty \leq g_j(x,y) \leq 0$.

Since different solvers define infinity differently, we use the OSI access function getInfinity() to fetch the proper value.

  int n_rows = M + N;
  double * row_lb = new double[n_rows]; //the row lower bounds
  double * row_ub = new double[n_rows]; //the row upper bounds
  CoinDisjointCopyN(demand, M, row_lb);
  CoinDisjointCopyN(demand, M, row_ub);
  CoinFillN(row_lb + M, N, -1 * si->getInfinity());
  CoinFillN(row_ub + M, N, 0.0);

Next, we define an empty row-ordered CoinPackedMatrix and set its column dimension to n_cols. See the constructor definitions in Coin/include/CoinPackedMatrix.hpp.

  CoinPackedMatrix * matrix = new CoinPackedMatrix(false,0,0);
  matrix->setDimensions(0, n_cols);

Then, we construct the demand constraints (1) using a CoinPackedVector. Each constraint entry consists of the column index and its coefficient value. Once the row is constructed, we then append it to the matrix using appendRow().

  for(i = 0; i < M; i++){
    CoinPackedVector row;
    for(j = 0; j < N; j++)
      row.insert(xindex(i,j), 1.0);
    matrix->appendRow(row);
  }

And, similarly for the linking constraints (2).

  for(j = 0; j < N; j++){
    CoinPackedVector row;
    row.insert(yindex(j), -1.0 * total_demand);
    for(i = 0; i < M; i++)
      row.insert(xindex(i,j), 1.0);
    matrix->appendRow(row);
  }

Now, we can simply load the problem data into the OsiSolver interface. See Osi/include/OsiSolverInterface.hpp for the different methods for loading (or assigning) a problem to the interface.

  si->loadProblem(*matrix, col_lb, col_ub, objective, row_lb, row_ub);

Finally, we tell the interface which variables are integral.

  si->setInteger(integer_vars, N);

Next, we will show how to solve this model using explicit branch and bound UFL::solve_model().