Different Optimal Solutions from Equivalent Quadratic and Direct Formulations
AnsweredHi all,
First of all, I hope I create this post in the correct community section.
I’m working on a mixed-integer quadratic program (MIQP) to minimize an objective function representing the squared difference between cluster sizes and assigned sensing values. The decision variables are binary, and I’ve implemented the objective in two mathematically equivalent forms: a quadratic form (using precomputed matrices) and a direct sum-of-squares form. However, I’m observing different optimal objective values (1253.0 vs. 1009.0) despite identical constraints and binary variables, and I’d like to understand why this happens and how to ensure consistency.
Here is the setup:
Objective:\[ f(x) = \sum_{j=1}^M (c_j - \sum_{i=1}^N s_i x_{ij})^2\]
\[ x_{ij} \in \{0, 1\} \] : Binary decision variable, N=10 agents, M=3 tasks.
\[ s_i \]: Sensing ranges (e.g., [5, 11, 10, 13, 14, 9, 12, 8, 7, 6]).
Direct Form: \[ f(x) = \sum_{j=1}^M (c_j - \sum_{i=1}^N s_i x_{ij})^2\]
and here the code:
import numpy as np
import gurobipy as gp
from gurobipy import GRB
def qp_objective_1(values, clstr_size):
sensing_range = np.array(values).reshape(1, -1)
cluster_size = clstr_size
task_num = len(cluster_size)
q = np.kron(np.eye(task_num), sensing_range)
Q = 2 * (q.T @ q)
p = -2 * (cluster_size.T @ q)
r = cluster_size @ cluster_size
return Q, p, r
def qp_allocation(values, cluster_size, Q1=None, p1=None, r1=None, use_quadratic=True):
n = len(values) * len(cluster_size)
model = gp.Model("qp_allocation")
model.Params.LogToConsole = 0
model.Params.Presolve = 2
model.params.MIPFocus = 3
model.Params.Method = 1
model.params.NonConvex = 2
model.setParam('MIPGap', 0.0001)
model.setParam('TimeLimit', 1800)
x = model.addMVar((n), vtype=GRB.BINARY, name="x")
A1 = np.kron(np.eye(len(values)), np.ones((1, len(cluster_size))))
A2 = np.kron(np.ones((1, len(values))), np.eye(len(cluster_size)))
b1 = np.ones((1, len(values))).T
b2 = np.ones((1, len(cluster_size))).T
model.addConstr(A1 @ x == b1.flatten(), name="c1")
model.addConstr(A2 @ x >= b2.flatten(), name="c2")
if use_quadratic:
obj = 0.5 * x @ Q1 @ x + p1 @ x + r1
model.setObjective(obj, GRB.MINIMIZE)
else:
obj = 0
n_tasks = len(cluster_size)
n_agents = len(values)
for j in range(n_tasks):
sum_sx = gp.quicksum(values[i] * x[i * n_tasks + j] for i in range(n_agents))
obj += (cluster_size[j] - sum_sx) * (cluster_size[j] - sum_sx)
model.setObjective(obj, GRB.MINIMIZE)
model.optimize()
if model.status == GRB.OPTIMAL or model.status == GRB.TIME_LIMIT:
x_val = x.X
print(f"Gurobi objective: {model.ObjVal}")
x_reshaped = x_val.reshape(len(values), len(cluster_size))
print(f"Solution:\n{x_reshaped}")
metric_1 = 0.0
for j in range(len(cluster_size)):
sum_sx = sum(values[i] * x_reshaped[i, j] for i in range(len(values)))
diff = cluster_size[j] - sum_sx
metric_1 += diff ** 2
print(f"Metric 1: {metric_1}")
return x_val, metric_1
return None, None
if __name__ == "__main__":
np.random.seed(123)
values = np.random.choice(np.arange(5, 15), size=10, replace=False)
cluster_size = np.array([50, 50, 50])
Q1, p1, r1 = qp_objective_1(values, cluster_size)
print("Quadratic Form:")
x_quad, metric_quad = qp_allocation(values, cluster_size, Q1, p1, r1, True)
print("\nDirect Form:")
x_direct, metric_direct = qp_allocation(values, cluster_size, use_quadratic=False)
It would be helpful for me to understand why the different results can happen. Or is it perhaps my modeling is not good or any different setting that I missed, etc., and how can I ensure Gurobi converges to the same (global) optimum, possibly the better one (1009.0)
I’d appreciate your insights on the issue that I have.
Thank you
-
-
Gurobi Staff
Hi Muhamad,
It is possible that the optimal solution point is the same but the objective value is different. A trivial example would be an objective with and without a constant.
What you should check first is whether the optimal solution point is the same for both formulations. This means that you should solve formulation 1 and save the optimal solution values and then solve formulation 2 and check the solution point values against the ones of the first optimization. This can be done in a similar manner to the following pseudo code.
# formulation 1
# ...
model.optimize()
# save optimal solution for formulation 1
sol_formulation_1 = model.getAttr("X", model.getVars())
model.reset()
# formulation 2
# ...
model.optimize()
# save optimal solution for formulation 2
sol_formulation_2 = model.getAttr("X", model.getVars())
for i in range(0, len(sol_formulation_1)):
if (math.abs(sol_formulation_1[i] - sol_formulation_2[i]) > 2e-5):
print("SOLUTION NOT EQUAL")Since you did not change the contraints, the optimal solution point of each formulation has to be feasible in the other one. You could use the solution of the first formulation as a MIP Start for the second one to check it's value.
Just to be sure, please always use the latest Gurobi version if possible.
Best regards,
Jaromił0 -
-
Hi Jaromil,
Thank you for your insight and suggestion!
I managed to get the quadratic form to have the same optimal solution as the direct form by ignoring the Kronecker product.
It is hilarious to know, as I thought using Kronecker product could simplify the expression and has no significant difference in the optimal solution.
But it really does, even though I can confirm the version with Kronecker product is equivalent to the direct form.
Thank you once again!
0 -
Please sign in to leave a comment.
Comments
2 comments