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:
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 cplex.addMinimize(cplex.diff(gmax, gmin)); // 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++) { expr2.addTerm(x[i][j], 1.0); } cplex.addEq(expr1, expr2); // size definition for team j expr1.clear(); for (int k = minSize - 1; k < maxSize; k++) { expr1.addTerm(s[j][k], 1.0); } 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++) { expr1.addTerm(x[i][j], gpa[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); expr1.addTerm(gmax, kk); expr2.addTerm(gmin, kk); expr1.addTerm(s[j][k], -4*maxSize); expr2.addTerm(s[j][k], 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 vars.add(teamVars[j]); bounds.add(1.0); dirs.add(IloCplex.BranchDirection.Up); } else if (j >= max) { // force variables above max to be 0 vars.add(teamVars[j]); bounds.add(0.0); dirs.add(IloCplex.BranchDirection.Down); } } // 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)
No comments:
Post a Comment
Due to intermittent spamming, comments are being moderated. If this is your first time commenting on the blog, please read the Ground Rules for Comments. In particular, if you want to ask an operations research-related question not relevant to this post, consider asking it on Operations Research Stack Exchange.