Nonlinear integer (0/1) programming solver

Maybe I am missing something but it looks like there is no need for a library: begin{align} sum_i sum_j sum_k x_{ji} y_{kj} cost(i,k)&=sum_i sum_j x_{ji} sum_k y_{kj} cost(i,k) end{align} Now since $sum_k y_{kj}=1$, exactly one row is 1, the others zero. We pick the best one: $$ =sum_i sum_j x_{ji} max_k cost(i,k)$$ Since $sum_j x_{ji}=1$ we get $$=sum_i max_k cost(i,k)$$

So basically go through the cost matrix row by row and pick the largest entry. This must have been homework :)

Besides the traditional linearization suggested by @MarkL.Stone and @Richard, you might consider using the constraints to obtain a compact linearization. Explicitly, multiply both sides of your second constraint by $x_{j,i}$: $$sum_k x_{j,i} y_{k,j} = x_{j,i}$$ Now replace $x_{j,i} y_{kj}$ with $z_{i,j,k}$ and impose an additional constraint to enforce $y_{k,j} = 0 implies z_{i,j,k} = 0$. The resulting linear formulation is:

begin{align} &text{maximize} &sum_isum_jsum_k text{cost}(i,k) z_{i,j,k}\ &text{subject to} &sum_j x_{j,i} &= 1 &&text{for all $i$}\ &&sum_k z_{i,j,k} &= x_{j,i} &&text{for all $i$ and $j$} \ &&0 le z_{i,j,k} &le y_{k,j} &&text{for all $i$, $j$, and $k$} end{align}

Your optimization problem seems like a special variant of the quadratic assignment problem (qap). One difference is that you have only products of variables coming from two different sets (x and y). This structure is called separable or disjoint bilinear programming.

The standard qap is one of the simplest quadratic problems, there are often examples from the solvers for this problem (some official/some from third-party people). They can easily be changed for your problem.

Localsolver:

• https://www.localsolver.com/docs/last/exampletour/qap.html

Gurobi:

• https://rpubs.com/HeatWave2019/560446
• https://gurobi.github.io/modeling-examples/customer_assignment/customer_assignment.html

In my opinion your best bet is to define an auxiliary variable $z_{ijk}$ such that: begin{equation} z_{ijk} geq x_{ji} + y_{kj} -1 \ z_{ijk} leq x_{ji} \ z_{ijk} leq y_{kj} end{equation}

Now this may become a really huge problem depending on the dimensions of $i$, $j$ and $k$. However, you gain the linearity of the problem which is worth a lot in my experience.

Lastly, you may be able to prune some stuff if you know more about the problem, but nothing comes to my mind right now.

You can try our own Octeract Engine, it solves all classes of optimisation problems up to (and including) DMINLP to global optimality. If you are a student/academic, you can also use it for it free! On the free side, there's also SCIP (only free for academics) and Couenne (open source).

Option 1: Submit as is to a solver which can globally optimize MIQPs having non-convex objective, and which might reformulate to a linearized MILP model under the hood. Such solvers include CPLEX, Gurobi 9.x, and BARON, among others.

Option 2:

Step 1 Linearize the products of binary variables, per How to linearize the product of two binary variables? . <Edit: This step is written out explicitly in this thread in @Richard 's subsequent answer.>

Step 2: Submit your linearized model to a MILP solver, such as CPLEX, Gurobi, XPress, Mosek, SCIP, or many others.

Note:Some solvers, such as CPLEX, give you the option of specifying whether the solver should reformulate the Binary MIQP to a MILP. There may be a default to let the solver decide which way is better. See indefinite MIQP switch:Decides whether CPLEX will attempt to reformulate a MIQP or MIQCP model that contains only binary variables.