Monday, February 5, 2018

Fun with a quadratic model


In [1] a type of assignment problem is described:

I'd like to match n ~= 100 people to m ~= 5 teams, subject to a couple of constraints. For each of the n people, I have a vector of k numbers that correspond to various skill ratings. Each person also has a preference rating for each of the 5 teams, as well as a list of people they'd like to work with and a list of people they would like not to work with. 
Given some criteria for a 'good' team (for example, some linear combination of skill levels of all its team members, everyone is happy with the team they're on (given their preferences), etc.) that can be used for evaluation, how should I design an algorithm to match the people to teams?
The foundation for an optimization model for this problem can be an assignment structure with:

\[x_{p,t} = \begin{cases} 1 & \text{ if person $p$ is assigned to team $t$}\\  0 & \text{ otherwise}\end{cases}\]

The assignment constraints can look like: \[\begin{align} &\sum_t x_{p,t} \le 1 & \forall p\\&\sum_p x_{p,t} = \mathit{demand}_{t}& \forall t \end{align} \] The problem can be viewed as a multi-criteria optimization problem with objectives: \[\begin{align} &\max\> f_{skill} = \sum_{p,s,t} \mathit{skill}_{p,s} x_{p,t} && \text{skill sets} \\ &\max\> f_{\mathit{preft}} = \sum_{p,t} \mathit{prefteam}_{p,t} x_{p,t} && \text{team preferences}\\&\max\>f_{\mathit{prefw}} = \sum_{p,p',t|p' >p} \mathit{prefworker}_{p,p'} x_{p,t} x_{p',t} && \text{worker preferences}\end{align}\] where \(s\) is the set of different skills. A value \( \mathit{prefworker}_{p,p'}\lt 0 \) indicates workers \(p\) and \(p'\) dislike each other. Note that in the quadratic expression we sum only over \(p' \gt p\) to avoid double counting. A typical way to handle competing objectives is to form a weighted sum: \[\max\> z = \sum_k w_k f_k\] So our first model M1 can look like: \[\bbox[lightcyan,10px,border:3px solid darkblue]{ \begin{align} \max\> &z = \sum_k w_k f_k\\ & f_{skill} = \sum_{p,s,t} \mathit{skill}_{p,s} x_{p,t}\\ & f_{\mathit{preft}} = \sum_{p,t} \mathit{prefteam}_{p,t} x_{p,t} \\ & f_{\mathit{prefw}} = \sum_{p,p',t|p' >p} \mathit{prefworker}_{p,p'} x_{p,t} x_{p',t} \\ &\sum_t x_{p,t} \le 1 \\ &\sum_p x_{p,t} = \mathit{demand}_{t}\\ &x_{p,t} \in \{0,1\} \end{align}}\]


  • M1: Not many solvers can handle this problem as stated: Cplex is the major exception.
  • M2: We can convince Gurobi and Xpress to solve this by changing the quadratic equality constraint into a \(\le\) constraint. (Interestingly Cplex has much more of a tough time with this version than with M1).
  • M3: This model is formed by substituting out the \(f\) variables. We end up with an MIQP model.
  • M4: As suggested by Bob Fourer below in the comments, we can simplify the quadratic term in the objective \[\sum_{p,p',t|p' >p} \mathit{prefworker}_{p,p'} x_{p,t} x_{p',t} \] by \[\begin{align} & \sum_{p,t} x_{p,t} y_{p,t}\\& y_{p,t} = \sum_{p' > p} \mathit{prefworker}_{p,p'} x_{p',t}\\& y_{p,t} \text{ free}   \end{align}\] Note that this is a MIQP model but with a simpler Q matrix than M3
  • M5: We can formulate this problem as a straight linear MIP by introducing binary variables \(y_{p,p',t}=x_{p,t} x_{p',t}\). This multiplication can be written as: \[\begin{align}&y_{p,p',t}\le x_{p,t}\\ &y_{p,p',t}\le x_{p',t}\\ &y_{p,p',t}\ge x_{p,t} + x_{p',t} -1\\ & f_{\mathit{prefw}} =  \sum_{p,p',t|p' >p} \mathit{prefworker}_{p,p'} y_{p,p',t}\\ & y_{p,p',t} \in \{0,1\} \end{align}\] This formulation removes any non-convexity from the problem.

In the table below I show how different solvers react on these different models. If you want your model to solve with as many solvers as possible, it is a good idea to reformulate this as a MIP model.

Model Solver Results
M1 (MIQCP)CplexOK
GurobiGurobi can only solve QCP with =L= or =G= quadratic rows. Could not load problem into Gurobi database (1)
MosekNonconvex model detected. Constraint 'objprefworker'(2) is a nonlinear equality. (*1290*)
XpressQuadratic equality constraint not supported: objprefworker. ERROR: Could not load problem into XPRESS optimizer.
M2 (MIQCP)CplexOK
GurobiOK
MosekConstraint 'objprefworker'(2) is not convex. Q should be positive semidefinite in a constraint with finite upper bound. (*1293*)
XpressOK with warning:
?900 Warning: The quadratic part of row 'R3' defines a nonconvex region. Please check your model or use Xpress-SLP. Could not solve final fixed LP. We report the unfixed solution.
M3 (MIQP)CplexOK
GurobiOK
MosekThe quadratic coefficient matrix in the objective is not negative semidefinite as expected for a maximization problem. (*1296*)
XpressOK
M4 (MIQP)CplexOK
GurobiOK
MosekThe quadratic coefficient matrix in the objective is not negative semidefinite as expected for a maximization problem. (*1296*)
Xpress?899 Warning: The quadratic objective is not convex
Please check model or use Xpress-SLP
Problem is integer infeasible
M5 (MIP)CplexOK
GurobiOK
MosekOK
XpressOK

Some of the messages are a bit over the top. If an LP and Q matrix  is called a 'database', what should we call a real database?

Often I make the argument that a MIP can perform better than a MIQP/MIQCP model. Here is another reason why such a reformulation may make sense: different solvers have different rules about what type of MIQCP/MIQP models they can solve. A formulation that works for one solver is no guarantee it will work for another solver. With a MIP formulation this is a little bit more predictable.

Performance


To show the unexpectedly large variation in performance of the models M1 through M5 we reproduce the Cplex results on a small randomly generated data set.

Model Rows Columns (discrete) Iterations Nodes Seconds Objective value
M1 109 504 (500) 3420 0 3.16 3178
M2 109 504 (500) > 10,000 (time interrupt)
M3 106 501 (500) 2977 0 4.13 3178
M4 606 1,001 (500) 25203 1504 6 3178
M5 74,359 25,254 (25,250) 19983 53 65.92 3178

Model M2 takes a very long time: I was unable to solve this to optimality in runs that took several hours. This is somewhat surprising: the difference between models M1 and M2 is just the sign of equation: \[f_{\mathit{prefw}}\left\{\begin{matrix}=\\ \le\end{matrix}\right\}\sum_{p,p',t|p' >p} \mathit{prefworker}_{p,p'} x_{p,t} x_{p',t}\] Some models are very sensitive to minor differences in the formulation.

References


  1. Matching people with skills and preferences to groups, https://stackoverflow.com/questions/48482550/matching-people-with-skills-and-preferences-to-groups