## Monday, June 25, 2012

### Birthing a Litter in CPLEX

In a previous post, I discussed (at a conceptual level) how to work around CPLEX's limitation of at most two children for any node in the search tree of a mixed-integer program. Following up on a question on one of their support forums, I'll add some implementation details and an example (using the Java API).

The example that I will use is a team assignment problem.  Given $N$ students (indexed $1\ldots N$), with student $i$ having a grade point average (GPA) equal to $g_i$ (on a scale from 0.0 to 4.0), form teams containing between $S_{min}$ and $S_{max}\ge S_{min}$ students on each team, so that the mean GPAs of the teams are as close as possible (minimal spread).

The use of mean GPA, combined with a variable team size (unless $S_{max}=S_{min}$), actually makes this a nonlinear model.  I'll finesse that by binary variables, in a Type 1 Special Ordered Set (SOS1), to control the size of each team.

### Mathematical Model

#### Parameters

• $N$ = number of students to assign
• $g_i\in [0.0,4.0]$ = GPA of student $i$
• $S_{min}$ = minimum allowable size of any team
• $S_{max}$ = maximum allowable size of any team
• $T_{max}=\left\lfloor \frac{N}{S_{min}}\right\rfloor$ = maximum number of teams needed
• $T_{min}=\left\lceil \frac{N}{S_{max}}\right\rceil$ = minimum number of teams needed

#### Variables

• $x_{ij}$ = 1 if student $i\in \{1,\dots,N\}$ is assigned to team $j\in \{1,\dots,T_{max}\}$, else 0
• $y_j$ = 1 if team $j\in \{1,\dots,T_{max}\}$ is used, else 0
• $s_{jk}$ = 1 if team $j\in \{1,\dots,T_{max}\}$ has size $k\in \{S_{min},\dots,S_{max}\}$, else 0
• $G_{max}$ = highest average (mean) GPA of any team created
• $G_{min}$ = lowest average (mean) GPA of any team created
• $G_j$ = sum of the GPAs of students assigned to team  $j\in \{1,\dots,T_{max}\}$

#### Objective

$\min G_{max} - G_{min}$

#### Constraints

• $\sum_{j=1}^{T_{max}} x_{ij}=1\quad\forall i\in \{1,\dots,N\}$ (1)
• $\sum_{i=1}^Nx_{ij}=\sum_{k=S_{min}}^{S_{max}}k*s_{jk}\quad\forall j\in \{1,\dots,T_{max}\}$ (2)
• $\sum_{k=S_{min}}^{S_{max}}s_{jk}=y_j\quad \forall j\in \{1,\dots,T_{max}\}$ (3)
• $y_j\le y_{j-1}\quad \forall j\in \{2,\dots,T_{max}\}$ (4)
• $G_j=\sum_{i=1}^Ng_i x_{ij}\quad\forall j\in \{1,\dots,T_{max}\}$ (5)
• $G_j \le k G_{max}+4S_{max}\left(1-s_{jk} \right) \quad \forall j\in \{1,\dots,T_{max}\}, k\in \{S_{min},\dots,S_{max}\}$ (6)
• $G_j \ge k G_{min}-4k\left(1-s_{jk} \right) \quad \forall j\in \{1,\dots,T_{max}\}, k\in \{S_{min},\dots,S_{max}\}$ (7)

Constraint (1) ensures that each student is assigned to exactly one team. Constraint (2) simultaneously defines the $s$ variables (team sizes) and ensures that every team is either empty or has a size between $S_{min}$ and $S_{max}$. Constraint (3) simultaneously makes the $s$ variables for each team a SOS1 and ensures that the size is nonzero if and only if the team is used ($y_j=1$). Constraint (4) eliminates some symmetry in the model by "packing" teams (so that, if fewer than $T_{max}$ teams are used, the ones used are those with lowest indices). Constraint (5) computes the sum of the GPAs of students assigned to each team.

Finally, constraints (6) and (7) define the smallest ($G_{min}$) and largest ($G_{max}$) team GPAs in a linear manner. If $s_{jk}=1$, then team $j$ has size $k$, and (6) and (7) resolve to $kG_{min}\le G_j \le kG_{max},$ or equivalently $G_{min}\le G_j/k \le G_{max},$ where $G_j/k$ is the mean GPA of team $j$. If $s_{jk}=0$, (6) and (7) are both nonbinding; given the 0-4 grading scale, the right side of (6) becomes $k G_{max}+4S_{max}\ge 4S_{max}$ (the highest possible sum of GPAs on any valid team), while the right side of (7) becomes $k\left(G_{min} - 4\right)\le 0$.

### Branching Strategy

To illustrate how to simulate more than two children at a node, I will adopt the following branching strategy. At the root node, I will branch on the number of teams being utilized (which must be between $T_{min}$ and $T_{max}$ inclusive), creating one child for each possible team count. At any other nodes, I will let CPLEX branch as it pleases. Note that I am not advocating this strategy as an effective method for problems of this nature; it is strictly a coding example.

Of course, I cannot literally create $T_{max}-T_{min}+1$ children at the root node, since CPLEX limits a node to at most children. So my branch callback will attach to each node a data object containing the fewest ($a$) and most ($b$) teams allowed at that node and its descendants.  When I branch on that node, I will bisect $\{a,\dots, b\}$ into $\{a,\dots,m\}$ and $\{m+1,\dots,b\}$, where $\m=\left\lfloor (a+b)/2 \right\rfloor$. If $a=m$ or $b=m+1$, I will create a "normal" child corresponding to that (with no attached data). If $a<m$ or $m+1 < b$, the corresponding child (which I will refer to as a "composite" node) will have the new range ($[a,m]$ or $[m+1,b]$) attached as node data, and will receive the same special treatment when the time comes to branch on it.

Suppose that $T_{min}=3$ and $T_{max}=7$. This is the top of the search tree we want:
The root node would be "composite", whereas all the children would be "normal". This, on the other hand, is the top of the tree we will generate (given CPLEX's two-child limit):
Here there are three "composite" nodes in addition to the root node; the five leaf nodes are "normal".

### Java Code

Shown next is the main program. (My usual disclaimer applies: I'm an amateur coder, and I make no claims that this is optimized for anything.) The entire program can be downloaded here (Java source code, ~5 KB). [Note: This link was updated on 10/02/2014 to point to a new repository.]

There is nothing particularly special about the main program. The GPA is generated randomly; model parameters are set at the top of the main method. I set the absolute tolerance (EpAGap parameter) to 0.01; there is no point in splitting hairs when working with GPA data.

/*
* Example of how to generate > 2 child nodes in CPLEX.
* Context: assigning students to teams to balance average GPA.
* Team sizes are constrained, and the head counts of the largest and smallest
* teams should differ by at most one.
*/
package multichild;

import ilog.concert.IloException;
import ilog.concert.IloLinearNumExpr;
import ilog.concert.IloNumVar;
import ilog.cplex.IloCplex;
import java.util.Random;

/**
* @author Paul A. Rubin (http://about.me/paul.a.rubin)
*/
public class MultiChild {

/**
* Main program
*/
public static void main(String[] args) {
// problem parameters
int nStudents = 60;  // number of students to form into teams
int minSize = 4;     // minimum team size
int maxSize = 6;    // maximum team size
int maxTeams = (int) Math.floor((1.0*nStudents)/minSize);   // max # of teams
int minTeams = (int) Math.ceil((1.0*nStudents)/maxSize);  // max # of teams
double[] gpa = new double[nStudents];
// generate random GPAs
Random gen = new Random(666L);
for (int i = 0; i < nStudents; i++) {
gpa[i] = Math.rint(200.0 * (1.0 + gen.nextDouble()))/100;
// using Math.rint to get GPAs with exactly two decimal places
}
/*
* The MIP model will look like the following, where
*   x[i][j] = 1 if student i is assigned to team j, else 0
*   y[j] = 1 if team j is used, else 0
*   s[j][k] = 1 if team j has size k (k = minSize, ..., maxSize), else 0
*   gmax = highest composite GPA of any team
*   gmin = lowest composite GPA of any team
*   gsum[j] = sum of GPAs of students on team j
*
* minimize gmax - gmin
* s.t.
*   sum_j x[i][j] = 1 for all i (assign every student once)
*   sum_k k*s[j][k] = sum_i x[i][j] for all j (size of each team)
*   sum_k s[j][k] = y[j] for all j  (size is SOS1; unused teams have size 0)
*   y[j] <= y[j - 1] for all j > 0 (unused teams come at end of list)
*   gsum[j] = sum_i gpa[i]*x[i][j] for all j (GPA sums)
*   gsum[j] <= k*gmax + 4*maxSize*(1 - s[j][k]) for all j, k (defines gmax)
*   gsum[j] >= k*gmin - 4*k*(1 - s[j][k]) for all j, k (defines gmin)
*/
try {
IloCplex cplex = new IloCplex();
// create the variables
IloNumVar[][] x = new IloNumVar[nStudents][maxTeams];
IloNumVar[] y = new IloNumVar[maxTeams];
IloNumVar[][] s = new IloNumVar[maxTeams][maxSize];
IloNumVar[] gsum = new IloNumVar[maxTeams];
IloNumVar gmax = cplex.numVar(0.0, 4.0, "gmax");
IloNumVar gmin = cplex.numVar(0.0, 4.0, "gmin");
for (int i = 0; i < nStudents; i++) {
for (int j = 0; j < maxTeams; j++) {
x[i][j] = cplex.boolVar("x_" + i + "_" + j);
}
}
for (int j = 0; j < maxTeams; j++) {
y[j] = cplex.boolVar("y_" + j);
gsum[j] = cplex.numVar(0.0, Double.MAX_VALUE, "gsum_" + j);
for (int k = minSize - 1; k < maxSize; k++) {
s[j][k] = cplex.boolVar("s_" + j + "_" + (k + 1) );
}
}
// add the objective function
// add the constraints
for (int i = 0; i < nStudents; i++) {
cplex.addEq(cplex.sum(x[i]), 1.0); // single assignment, student i
}
for (int j = 0; j < maxTeams; j++) {
IloLinearNumExpr expr1 = cplex.linearNumExpr();
IloLinearNumExpr expr2 = cplex.linearNumExpr();
for (int k = minSize - 1; k < maxSize; k++) {
expr1.addTerm(s[j][k], k + 1);
}
for (int i = 0; i < nStudents; i++) {
}
cplex.addEq(expr1, expr2); // size definition for team j
expr1.clear();
for (int k = minSize - 1; k < maxSize; k++) {
}
cplex.addEq(expr1, y[j]); // size is SOS1; unused teams have size 0
if (j > 0) {
cplex.addLe(y[j], y[j - 1]);  // used teams come first
}
expr1.clear();
for (int i = 0; i < nStudents; i++) {
}
cplex.addEq(expr1, gsum[j]);  // define GPA sums
for (int k = minSize - 1; k < maxSize; k++) {
int kk = k + 1;
expr1.clear();
expr2.clear();
expr1.setConstant(4*maxSize);
expr2.setConstant(-4*kk);
cplex.addLe(gsum[j], expr1); // upper bound on GPA sum
cplex.addGe(gsum[j], expr2); // lower bound on GPA sum
}
}
// set a tolerance so that we don't split hairs
cplex.setParam(IloCplex.DoubleParam.EpAGap, 0.01);
// attach a branch callback -- comment this out to run without
// special branching
cplex.use(new Brancher(minTeams, maxTeams, y));
// solve the model
cplex.solve();
if (cplex.getStatus() == IloCplex.Status.Optimal) {
// fetch the solution
double[][] xx = new double[nStudents][];
for (int i = 0; i < nStudents; i++) {
xx[i] = cplex.getValues(x[i]);
}
double[] yy = cplex.getValues(y);
double[] ggsum = cplex.getValues(gsum);
// display the optimal solution
System.out.println("\nAssign " + nStudents + " students into teams "
+ "of size " + minSize + " to " + maxSize + ".");
System.out.println("GPA spread is from " + cplex.getValue(gmin)
+ " to " + cplex.getValue(gmax)
+ " (a spread of " + cplex.getObjValue() + ")");
System.out.println("CPLEX explored " + cplex.getNnodes() + " nodes.");
System.out.println("\nTeam\tStudent\tGPA");
for (int j = 0; j < maxTeams; j++) {
if (yy[j] > 0.5) {
int size = 0;
for (int i = 0; i < nStudents; i++) {
if (xx[i][j] > 0.5) {
size++;
System.out.println(j + "\t" + i + "\t" + gpa[i]);
}
}
System.out.println("\tOVERALL\t" + ggsum[j]/size + "\n");
}
}
}
else {
System.err.println("CPLEX returned status " + cplex.getStatus());
System.exit(1);
}
} catch (IloException ex) {
System.err.println("CPLEX pitched a hissy fit:\n" + ex.getMessage());
System.exit(1);
} catch (Exception ex) {
System.err.println("Could not construct the callback!\n" + ex.getMessage());
System.exit(1);
}
}
}


The information object attached to composite nodes is a simple container class that provides just the minimum and maximum number of teams at that node:
/*
* BranchInfo defines an object that can be attached to nodes and guides
* branching decisions.
*/
package multichild;

/**
* @author Paul A. Rubin (http://about.me/paul.a.rubin)
*/
public class BranchInfo {
private int minTeams;  // minimum number of teams to use in all child nodes
private int maxTeams;  // maximum number of teams to use in all child nodes

/**
* Constructor.
* @param min minimum number of teams to be used
* @param max maximum number of teams to be used
*/
public BranchInfo(int min, int max) {
this.minTeams = min;
this.maxTeams = max;
}

/**
* Get the maximum number of teams to create among child nodes.
* @return maximum number of teams to create among child nodes
*/
public int getMaxTeams() {
return maxTeams;
}

/**
* Get the minimum number of teams to create among child nodes.
* @return minimum number of teams to create among child nodes
*/
public int getMinTeams() {
return minTeams;
}

}


Finally, I present below the branch callback implementation. It contains a bunch of extra print statements just to illustrate how it is behaving. The callback is called at every node. At a "normal" node, where the number of teams is already fixed, it returns without doing anything (allowing CPLEX to branch as usual). At a "composite" node, it produces two children as described above (other than in the pathological case of $T_{min}=T_{max}$, in which case it produces a single "normal" child).

The callback detects a "composite" node by the presence of an information object attached (by an earlier call to the callback) to that node. The one exception is the root node, which we know is composite but to which we cannot attach an information object. I handle that by setting a flag in the callback (at the time of construction) saying that the next node encountered will be the root. The first time branching occurs (which must be at the root), the callback recognizes it as a "composite" node and clears the flag.
/*
* Custom branch callback creating 3+ children.
*
* When branching at the root node, this will create one child node for
* each possible count on the number of teams. At other nodes, it will
* let CPLEX branch as usual.
*/
package multichild;

import ilog.concert.IloException;
import ilog.concert.IloNumVar;
import ilog.cplex.IloCplex;
import java.util.ArrayList;

/**
* @author Paul A. Rubin (http://about.me/paul.a.rubin)
*/
public class Brancher extends IloCplex.BranchCallback {
private int minTeams;  // minimum number of teams to create
private int maxTeams;  // maximum number of teams to create
private IloNumVar[] teamVars;  // boolean variables for team use
private boolean atRoot;  // are we at the root node?
private ArrayList<Double> bounds; // bounds for variables when branching
private ArrayList<IloCplex.BranchDirection> dirs; // branch directions
private ArrayList<IloNumVar> vars;  // variables involved in branching
private double[] aBds;  // bounds as doubles
private IloCplex.BranchDirection[] aDirs;  // directions as array
private IloNumVar[] aVars;  // variables as array

/**
* Constructor.
* @param minTeams  minimum number of teams to create
* @param maxTeams  maximum number of teams to create
* @param v         boolean variables for team use
* @throws an exception if the variable array has the wrong dimension
*/
public Brancher(int minTeams, int maxTeams, IloNumVar[] v) throws Exception {
super();
this.minTeams = minTeams;
this.maxTeams = maxTeams;
this.teamVars = v;
// make sure the variable vector has the correct dimension
if (v.length != maxTeams) {
throw new Exception("Variable array dimension " + v.length
+ " does not match maxTeams = " + maxTeams);
}
this.atRoot = true;
}

/**
* The actual callback operation.
* @throws IloException if CPLEX gets miffed
*/
@Override
protected void main() throws IloException {
BranchInfo info = null;
Object raw;
String msg = ">>> At node " + getNodeId() + ", ";
IloCplex.NodeId id;
// check whether we are at the root node (which will have no attachment)
if (atRoot) {
raw = new BranchInfo(minTeams, maxTeams);  // create an attachment
atRoot = false;  // clear the root flag
}
else {
raw = getNodeData();  // get the current node's attachment, if any
}
if (raw == null) {
// no data object -- let CPLEX handle branching
return;
}
else if (raw instanceof BranchInfo) {
// convert the object to an instance of BranchInfo
info = (BranchInfo) raw;
}
else {
// unknown node data type -- should never happen
System.err.println("Encountered unknown node data type!");
abort();
}
// the presence of node data means we are branching on team count
// if the range of team sizes at this node is one or two, just create
// normal children; otherwise, bisect the range of counts and create
// two "special" children
int min = info.getMinTeams();
int max = info.getMaxTeams();
msg += "team range is " + min + " to " + max;
double obj = getObjValue();  // current node bound
switch (max - min) {
case 0: // create a single child with exactly this many teams
// this case should only occur (at the outset) if minTeams = maxTeams
cut(min, max);
id = makeBranch(aVars, aBds, aDirs, obj); // attach no data
msg += " -- adding one normal child ("
+ id + ") with exactly " + min + " teams";
System.out.println(msg);
return;
case 1: // create two normal children
// first child contains exactly min teams, no attached data
cut(min, min);
id = makeBranch(aVars, aBds, aDirs, obj); // attach no data
msg += " -- adding two normal children (" + id + " using exactly "
+ min + " teams and ";
// second child contains exactly max = min + 1 teams, no attachment
cut(max, max);
id = makeBranch(aVars, aBds, aDirs, obj); // attach no data
msg += id + " using exactly " + max + " teams)";
System.out.println(msg);
return;
case 2: // create one normal child and one composite child
// first child uses min to min+1 teams
cut(min, min + 1);
// attach node info to this child
info = new BranchInfo(min, min + 1);
id = makeBranch(aVars, aBds, aDirs, obj, info);
msg += " -- adding one composite child (" + id + ") with range "
+ min + " to " + (min + 1);
// second child uses exactly max = min + 2 teams, no attached data
cut(max, max);
id = makeBranch(aVars, aBds, aDirs, obj);
msg += " and one normal child (" + id + ") using exactly "
+ max + " teams";
System.out.println(msg);
return;
default: // max > min + 2: create two composite children
int mid = (max + min)/2;  // midpoint of size range (rounding down)
// first child uses min to mid team count
cut(min, mid);
info = new BranchInfo(min, mid);
id = makeBranch(aVars, aBds, aDirs, obj, info);
msg += " -- adding two composite children (" + id
+ " with range " + min + " to " + mid;
// second child uses mid + 1 to max team count
cut(mid + 1, max);
info = new BranchInfo(mid + 1, max);
id = makeBranch(aVars, aBds, aDirs, obj, info);
msg += " and " + id + " with range " + (mid + 1) + " to " + max + ")";
System.out.println(msg);
return;
}
}

/**
* Sets up the arguments for a cut that ensures somewhere between min and
* max teams are used.
* @param min minimum number of teams to use
* @param max maximum number of teams to use
*/
private void cut(int min, int max) {
this.bounds = new ArrayList<Double>();
this.dirs = new ArrayList<IloCplex.BranchDirection>();
this.vars = new ArrayList<IloNumVar>();
for (int j = 0; j < maxTeams; j++) {
if (j < min) {  // force the first min variables to be 1
}
else if (j >= max) {  // force variables above max to be 0
}
}
// convert the lists to arrays
aVars = vars.toArray(new IloNumVar[0]);
aDirs = dirs.toArray(new IloCplex.BranchDirection[0]);
// unavoidable PITA -- bounds must be converted from Double[] to double[]
// with a loop
Double[] temp = bounds.toArray(new Double[0]);
aBds = new double[temp.length];
for (int i = 0; i < temp.length; i++) {
aBds[i] = temp[i];
}
}

} 


Related post: "Multiple Children - Again!" (updated code, changing the composite nodes to clones of their parents)