# Numerical Simulation of a Quadratic MIP with a highly rational term

I am interested in solving the following minimization problem: $$\begin{array}{cl} \displaystyle\min_{x,y}&\displaystyle\frac{1}{K}\sum_{i=1}^{K}\left(\frac{x_{i}}{y_{i}}-\frac{X}{Y}\right)^{2} \\[0.4cm] \mathrm{s.t.}&\sum_{i}x_{i}=X \\[0.2cm] &\sum_{i}y_{i}=Y \\[0.2cm] &x_{i}\leqslant y_{i},\quad i=1,2,\ldots,K \\[0.2cm] &x_{i},y_{i}\in\mathbb{Z}^{+},\quad i=1,2,\ldots, K \end{array} \tag{1}$$ with $$X$$ and $$Y$$ being fixed.

I attempted to simulate $$(1)$$ using $$\textsf{cvxpy}$$. However, it did not work as the quadratic term $$\frac{x_{i}^{2}}{y_{i}^{2}}$$ is neither convex, nor concave, and $$\textsf{cvxpy}$$ has no framework for such cases.

Question: How can we reformulate this problem so that it can be simulated by $$\textsf{cvx/cvxpy}$$, or what QMIP solver is able to handle an objective function of this form being highly rational and quadratic?

The solution to this minimization problem all boils down to $$x_{i}=\lceil y_{i}\frac{X}{Y}\rceil$$ or something close to that. For what its worth, I attempted to reformulate this problem by adding an envelop through $$f_{i}:=\frac{x_{i}}{y_{i}}$$. This leads to the following minimization problem: $$\begin{array}{cl} \displaystyle\min_{f}&\displaystyle\frac{1}{K}\sum_{i=1}^{K}\left(f_{i}-\frac{X}{Y}\right)^{2} \\[0.4cm] \mathrm{s.t.}&y_{i}f_{i}-x_{i}=0 \\[0.2cm] &\sum_{i}x_{i}=X \\[0.2cm] &\sum_{i}y_{i}=Y \\[0.2cm] &f_{i}\leqslant 1,\quad i=1,2,\ldots,K \\[0.2cm] &x_{i},y_{i}\in\mathbb{Z}^{+},\quad i=1,2,\ldots, K \end{array} \tag{2}$$ However this formulation in $$\textsf{cvxpy}$$ was to no avail as there is no guided way to handle cases with bilinear constraints.

• Didn't you post this elsewhere? Commented Jun 1, 2023 at 9:20
• I did cross-post it on cvxpy forum a couple of hours ago (ask.cvxr.com/t/…), and they recommended me to post it here. @RodrigodeAzevedo Commented Jun 1, 2023 at 9:23
• yes, that was a couple of days ago, I had problem configuring the optimizer solver, and it involved a different objective function (convex). @RodrigodeAzevedo Commented Jun 1, 2023 at 10:43
• What are representative values for your constants? Commented Jun 4, 2023 at 19:32
• you can think of them as a resource allocation problem. For instance, $x$ could be how much blue balls you have and $y$ could be how balls can fit in a basket. @Richard Commented Jun 4, 2023 at 19:50

# Getting an exact solution via brute force

I'll try to circle back later to formulate this as a MIP, but your problem as stated is small enough that you can just brute force the solution.

For instance, for K=3, X=941, Y=1223, the solution is

Best val = 8.22518e-12
xi =  307  317  317
yi =  399  412  412


Using 12 cores this takes 8.5 minutes to find using the code below.

This code finds said solution:

// Compile with: g++ -O3 42886-quadratic-mip.cpp -Wall -Wextra -fopenmp
#include <cassert>
#include <iostream>
#include <limits>
#include <iomanip>
#include <tuple>
#include <vector>

void generate_sum_sets_helper(
const int number_of_vars,
const int X,
const bool allow_zeros,
const int Xtotal,
std::vector<int>& running_sol,
const int running_sum,
std::vector<std::vector<int>>& solutions
){
if(X < 0){
return;
} else if(number_of_vars==0){
if(running_sum==Xtotal){
solutions.push_back(running_sol);
}
} else if(X==0){
if(allow_zeros && running_sum==Xtotal){
running_sol.push_back(0);
generate_sum_sets_helper(number_of_vars-1, X, allow_zeros, Xtotal, running_sol, running_sum, solutions);
running_sol.pop_back();
}
} else {
for(int i=(allow_zeros?0:1);i<=X;i++){
running_sol.push_back(i);
generate_sum_sets_helper(number_of_vars-1,X-i,allow_zeros,Xtotal,running_sol,running_sum+i, solutions);
running_sol.pop_back();
}
}
}

std::vector<std::vector<int>> generate_sum_sets(
const int number_of_vars,
const int X,
const bool allow_zeros
){
std::vector<int> running_sol;
std::vector<std::vector<int>> solutions;
generate_sum_sets_helper(number_of_vars, X, allow_zeros, X, running_sol, 0, solutions);
return solutions;
}

double objective(
const std::vector<int>& xvals,
const std::vector<int>& yvals,
const int X,
const int Y
){
assert(xvals.size()==yvals.size());
double sum = 0;
for(size_t i=0;i<xvals.size();i++){
const auto temp = xvals[i] / static_cast<double>(yvals[i]) - X / static_cast<double>(Y);
sum += temp * temp;
}
return sum / xvals.size();
}

bool all_less(
const std::vector<int>& xvals,
const std::vector<int>& yvals
){
assert(xvals.size()==yvals.size());
for(size_t i=0;i<xvals.size();i++){
if(xvals[i]>=yvals[i]){
return false;
}
}
return true;
}

int main(){
// Using prime numbers to prevent an exact solution
constexpr auto Xmax =  941;
constexpr auto Ymax = 1223;
constexpr auto K = 3;
static_assert(Xmax < Ymax);
std::cerr<<"Generating xvals..."<<std::endl;
const auto xvals = generate_sum_sets(K, Xmax, true);
std::cerr<<"Generated "<< xvals.size() << " xvals..."<<std::endl;
std::cerr<<"Generating yvals..."<<std::endl;
const auto yvals = generate_sum_sets(K, Ymax, false);
std::cerr<<"Generated "<< yvals.size() << " yvals..."<<std::endl;
std::tuple<double, std::vector<int>, std::vector<int>> bestval = {std::numeric_limits<double>::infinity(), {}, {}};
auto best_overall = bestval;
std::cerr<<"Searching for solutions..."<<std::endl;
#pragma omp parallel firstprivate(bestval)
{
#pragma omp for collapse(2)
for(size_t xi=0;xi<xvals.size();xi++){
for(size_t yi=0;yi<yvals.size();yi++){
if(!all_less(xvals.at(xi),yvals.at(yi))){
continue;
}
const auto oval = objective(xvals.at(xi),yvals.at(yi),Xmax, Ymax);
if(oval < std::get<0>(bestval)){
std::cerr<<"Best objective = "<<std::setprecision(10)<<oval<<std::endl;
bestval = std::make_tuple(oval, xvals.at(xi), yvals.at(yi));
}
}
}

#pragma omp critical
{
if(std::get<0>(bestval) < std::get<0>(best_overall)){
best_overall = bestval;
}
}
}
std::cout<<"Best val = "<<std::get<0>(best_overall)<<std::endl;
std::cout<<"xi = ";
for(const auto &x: std::get<1>(best_overall)){
std::cout<<std::setw(4)<<x<<" ";
}
std::cout<<std::endl;
std::cout<<"yi = ";
for(const auto &x: std::get<2>(best_overall)){
std::cout<<std::setw(4)<<x<<" ";
}
std::cout<<std::endl;

return 0;
}


# cvxpy doesn't work

In cvxpy you would write your envelope formulation as:

import cvxpy as cp

# Create variables
K = 3
X = 941
Y = 1223

xi = cp.Variable(K, integer=True)
yi = cp.Variable(K, integer=True)
fi = cp.Variable(K)

constraints = [
sum(xi) == X,
sum(yi) == Y,
xi >= 0,
yi >= 1,
xi <= X,
yi <= Y,
fi >= 0,
fi <= 1,
]
for i in range(K):
constraints.append(xi[i] <= yi[i]-1)
for i in range(K):
constraints.append(yi[i]*fi[i] - xi[i] == 0)

objective_expr = cp.Minimize(sum([(fi[i] - X / Y)**2 for i in range(K)]) / K)

# Optimize the model
model = cp.Problem(objective_expr, constraints)
optval = model.solve()
print(optval)


but this doesn't work because the problem isn't convex:

cvxpy.error.DCPError: Problem does not follow DCP rules. Specifically:
The following constraints are not DCP:
var2[0] @ var3[0] + -var1[0] == 0.0 , because the following subexpressions are not:
|--  var2[0] @ var3[0]
var2[1] @ var3[1] + -var1[1] == 0.0 , because the following subexpressions are not:
|--  var2[1] @ var3[1]
var2[2] @ var3[2] + -var1[2] == 0.0 , because the following subexpressions are not:
|--  var2[2] @ var3[2]


# Getting a maybe-good solution using MIPs in Gurobi

So, instead, we switch to using Gurobi, which is free with an academic license and otherwise will cost a blood sacrifice beneath a full moon. Note that you could try using alternatives like XPRESS.

In Gurobi the problem is formulated as follows:

import gurobipy as gp
from gurobipy import GRB

# Create a new model
model = gp.Model("qp")

#Set parameters
model.params.NonConvex = 2
model.setParam('FeasibilityTol', 1e-9)
model.setParam('OptimalityTol', 1e-9)
model.setParam('MIPGap', 1e-13)
model.setParam('IntFeasTol', 1e-9)
model.setParam('NumericFocus', 3) # VERY VERY IMPORTANT

# Create variables
K = 3
X = 941
Y = 1223

xi = model.addVars(K, lb=0, ub=X, vtype=GRB.INTEGER, name="x")
yi = model.addVars(K, lb=1, ub=Y, vtype=GRB.INTEGER, name="y")
fi = model.addVars(K, lb=0, ub=1, vtype=GRB.CONTINUOUS, name="f")

model.addConstr(sum(xi[i] for i in range(K)) == X, name="constraint_sum_x")
model.addConstr(sum(yi[i] for i in range(K)) == Y, name="constraint_sum_y")
for i in range(K):
for i in range(K):
model.addConstr(yi[i]*fi[i] - xi[i] == 0, name=f"constraint_f{i}")

objective_expr = sum((fi[i] - X / Y)**2 for i in range(K)) / K
model.setObjective(objective_expr, GRB.MINIMIZE)

# Optimize the model
model.optimize()

# Get the solution
if model.status == GRB.OPTIMAL:
print(f"min_value = {model.objVal}")
print(f"x_values = {[xi[i].x for i in range(K)]}")
print(f"y_values = {[yi[i].x for i in range(K)]}")
# Double check that we got the objective right
print(f"Objective: ", sum([(xi[i].x/yi[i].x - X/Y)**2 for i in range(K)])/K)
else:
print("Could not solve!")

model.printStats()
model.printQuality()


Note that there's a line I've labeled as VERY VERY IMPORTANT.

If we remove that line after 1.62s Gurobi confidently gives us:

Optimal solution found (tolerance 1.00e-13)
Best objective 5.844680295297e-11, best bound 5.844680295297e-11, gap 0.0000%
min_value = 5.844680295297167e-11
x_values = [327.0, 297.0, 317.0]
y_values = [425.0, 386.0, 412.0]
Objective:  5.844692423821425e-11

Solution quality statistics for model qp :
Maximum violation:
Bound       : 0.00000000e+00
Constraint  : 0.00000000e+00
Integrality : 0.00000000e+00


This is wrong! We know this because we brute forced the solution on a small problem.

We need to tighten every tolerance to its tightest value and set numeric focus to get the right answer:

model.setParam('FeasibilityTol', 1e-9)
model.setParam('OptimalityTol', 1e-9)
model.setParam('MIPGap', 1e-13)
model.setParam('IntFeasTol', 1e-9)
model.setParam('NumericFocus', 3) # VERY VERY IMPORTANT


and after 0.39s Gurobi hands that to us:

Optimal solution found (tolerance 1.00e-13)
Best objective 8.225198300238e-12, best bound 8.225198300238e-12, gap 0.0000%
min_value = 8.225198300237935e-12
x_values = [317.0, 307.0, 317.0]
y_values = [412.0, 399.0, 412.0]
Objective:  8.225179646366615e-12

Solution quality statistics for model qp :
Maximum violation:
Bound       : 0.00000000e+00
Constraint  : 0.00000000e+00
Integrality : 0.00000000e+00


# What have we learned?

What we've learned is that Gurobi is very fast in comparison to the brute force solution because it's able to use the problem's struture and its advanced heuristics to ignore enormous parts of the search space. We've also learned that it will confidently hand us a suboptimal answer because the problem you've formulated has many answers that are very close to optimality and Gurobi has difficulty distinguishing between them unless we really push it to be careful.

With large values of K I'd certainly become worried that even with the settings used above I might not be getting the very best answer. Of course, with large values of K a brute force solution will also take much longer.

• Thank you very much for your answer! I will also try to take this code and play with it to see how I can optimize speed. Typically, I am aiming to loop over $X=1,2,\ldots,Y-1$ and execute the code that will find the minimum for each $X$. That sounds pretty exhaustive!, but I would appreciate it if we can tackle the issue of speeding maybe brute force is a good plan B if eventually we dont find a way to represent it as a MIP Commented Jun 5, 2023 at 10:29
• @SPARSE - I've added a MIP-based solution. Commented Jun 5, 2023 at 16:53