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)	Cplex	OK
	Gurobi	Gurobi can only solve QCP with =L= or =G= quadratic rows. Could not load problem into Gurobi database (1)
	Mosek	Nonconvex model detected. Constraint 'objprefworker'(2) is a nonlinear equality. (*1290*)
	Xpress	Quadratic equality constraint not supported: objprefworker. ERROR: Could not load problem into XPRESS optimizer.
M2 (MIQCP)	Cplex	OK
	Gurobi	OK
	Mosek	Constraint 'objprefworker'(2) is not convex. Q should be positive semidefinite in a constraint with finite upper bound. (*1293*)
	Xpress	OK 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)	Cplex	OK
	Gurobi	OK
	Mosek	The quadratic coefficient matrix in the objective is not negative semidefinite as expected for a maximization problem. (*1296*)
	Xpress	OK
M4 (MIQP)	Cplex	OK
	Gurobi	OK
	Mosek	The 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)	Cplex	OK
	Gurobi	OK
	Mosek	OK
	Xpress	OK

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