In my master thesis, I chose a topic in combinatorial optimization. It introduces the Quadratic Maximum-Weight Independent Set (Q-MWIS) Problem as the quadratic extension of the Maximum-Weight Independent Set (MWIS) Problem, discusses its properties and suggests a linearization of the quadratic problem according to a scheme by Sherali and Adams. After optimizing the structure of the linearization, it is tested in practice and compared to alternative formulations, using the Python API of the optimization software Gurobi (external link). The testing showed that the constructed linearization outperforms both the trivial linearization and the default quadratic formulation in most cases.
To get a broad idea of the content of my thesis, its applications and the methodology used in the practical part, one may also skim through the slides of my master thesis presentation.
# Knowledgedump.org - qmwis_algorithm.py
# File containing the Gurobi model setup functions.
# >>>> To execute tests run qmwis.py <<<<
import gurobipy
import numpy as np
# Initializing Gurobi Model Q-MWIS1 with:
# - inequality constraints.
# - without linearization.
# Function that takes in costs and conflict sets. Returns Gurobi model object.
# -> costs_una is numpy array for unary cost terms.
# -> costs_pw is an upper triangle matrix, containing all pairwise cost terms.
# Implemented as a numpy matrix, with non-zero entries only for indices (i,j) with i<j.
# -> conflict sets in the form of numpy array, containing indices;
# each row corresponds to one conflict set.
def init_qmwis_inequality(costs_una, costs_pw, conflicts):
# Define Gurobi model object.
mod = gurobipy.Model()
# Set number of decision variables.
varNum = len(costs_una)
# Add variables to the model. Variables set to be binary and restricted to 0 <= x <= 1.
vars = mod.addMVar(varNum, name="x", lb = 0, ub = 1, vtype = gurobipy.GRB.BINARY)
# Updating is lazy, i.e. model has to be updated after changes are made.
mod.update()
# Set Objective function.
mod.setObjective(vars@costs_una + vars@costs_pw@vars, gurobipy.GRB.MAXIMIZE)
# Set conflict set constraints.
tempvars = mod.getVars()
for con in conflicts:
lhs = gurobipy.LinExpr()
for i in con:
lhs += tempvars[i]
mod.addConstr( lhs <= 1 )
mod.update()
return mod
# Initializing Gurobi Model Q-MWIS2 with:
# - equality constraints.
# - without linearization.
# Function that takes in costs and conflict sets. Returns Gurobi model object.
# Input same as init_qmwis_inequality.
def init_qmwis_equality(costs_una, costs_pw, conflicts):
mod = gurobipy.Model()
varNum = len(costs_una)
vars = mod.addMVar(varNum, name="x", lb = 0, ub = 1, vtype = gurobipy.GRB.BINARY)
# Add slack variables, to add in each conflict set.
vars_slack = mod.addVars(len(conflicts), name="x_s", lb = 0, ub = 1, vtype = gurobipy.GRB.BINARY)
# Updating is lazy, i.e. model has to be updated after changes are made.
mod.update()
# Set Objective function.
mod.setObjective(vars@costs_una + vars@costs_pw@vars, gurobipy.GRB.MAXIMIZE)
# Set conflict set constraints.
tempvars = mod.getVars()
j = 0
for con in conflicts:
lhs = gurobipy.LinExpr(vars_slack[j])
j += 1
for i in con:
lhs += tempvars[i]
mod.addConstr( lhs == 1 )
mod.update()
return mod
# Initializing Gurobi Model Q-MWIS3 with:
# - inequality constraints.
# - with trivial linearization.
# Function that takes in costs and conflict sets. Returns Gurobi model object.
# Input same as init_qmwis_inequality.
def init_qmwis_inequality_triv(costs_una, costs_pw, conflicts):
mod = gurobipy.Model()
varNum = len(costs_una)
# Create index list of integer tuples (i,j) with i<j, indexing non-zero pairwise cost entries.
index_pw = tuple(map(tuple, np.transpose(np.nonzero(costs_pw))))
# Add unary and pairwise variables to the model. Variables set to be binary restricted to 0 <= x <= 1.
vars_una = mod.addVars(varNum, name="x", lb = 0, ub = 1, vtype = gurobipy.GRB.BINARY)
vars_pw = mod.addVars(index_pw, name="y", lb = 0, ub = 1, vtype = gurobipy.GRB.BINARY)
mod.update()
# Set Objective function.
obj = gurobipy.LinExpr(costs_una, vars_una.select())
for ind in index_pw:
obj += costs_pw[ind] * vars_pw[ind]
mod.setObjective(obj, gurobipy.GRB.MAXIMIZE)
# Set conflict set constraints.
for con in conflicts:
lhs = gurobipy.LinExpr()
for i in con:
lhs += vars_una[i]
mod.addConstr( lhs <= 1 )
# Set constraints enforcing y[i,j] = x_i*x_j.
for ind in index_pw:
mod.addConstr( vars_una[ind[0]] - vars_pw[ind] >= 0 )
mod.addConstr( vars_una[ind[1]] - vars_pw[ind] >= 0 )
mod.addConstr( vars_pw[ind] - vars_una[ind[0]] - vars_una[ind[1]] + 1 >= 0 )
mod.update()
return mod
# Initializing Gurobi Model Q-MWIS4 with:
# - equality constraints.
# - with trivial linearization.
# Function that takes in costs and conflict sets. Returns Gurobi model object.
# Input same as init_qmwis_inequality.
def init_qmwis_equality_triv(costs_una, costs_pw, conflicts):
mod = gurobipy.Model()
varNum = len(costs_una)
# Create index list of integer tuples (i,j) with i<j, indexing non-zero pairwise cost entries.
index_pw = tuple(map(tuple, np.transpose(np.nonzero(costs_pw))))
# Add unary and pairwise variables to the model. Variables set to be binary restricted to 0 <= x <= 1.
vars_una = mod.addVars(varNum, name="x", lb = 0, ub = 1, vtype = gurobipy.GRB.BINARY)
vars_pw = mod.addVars(index_pw, name="y", lb = 0, ub = 1, vtype = gurobipy.GRB.BINARY)
# Add slack variables, to add in each conflict set.
vars_slack = mod.addVars(len(conflicts), name="x_s", lb = 0, ub = 1, vtype = gurobipy.GRB.BINARY)
mod.update()
# Set Objective function.
obj = gurobipy.LinExpr(costs_una, vars_una.select())
for ind in index_pw:
obj += costs_pw[ind] * vars_pw[ind]
mod.setObjective(obj, gurobipy.GRB.MAXIMIZE)
# Set conflict set constraints.
j = 0
for con in conflicts:
lhs = vars_slack[j]
j += 1
for i in con:
lhs += vars_una[i]
mod.addConstr( lhs == 1 )
# Set constraints enforcing y[i,j] = x_i*x_j.
for ind in index_pw:
mod.addConstr( vars_una[ind[0]] - vars_pw[ind] >= 0 )
mod.addConstr( vars_una[ind[1]] - vars_pw[ind] >= 0 )
mod.addConstr( vars_pw[ind] - vars_una[ind[0]] - vars_una[ind[1]] + 1 >= 0 )
mod.update()
return mod
#Initializing Gurobi Model Q-MWIS5 with:
#- equality constraints.
#- Sherali-Adams linearization for problems in standard form.
#Function that takes in costs and conflict sets. Returns Gurobi model object.
#Input same as init_qmwis_inequality.
def init_qmwis_inequality_sherali(costs_una, costs_pw, conflicts):
mod = gurobipy.Model()
# Set number of variables.
varNum = len(costs_una)
# Create index list of integer tuples (i,j) with i < j for pairwise terms w_ij.
index_pw = []
for i in range(varNum-1):
for j in range(i+1,varNum):
index_pw.append((i,j))
# Add unary and pairwise variables to the model. Variables set to be binary restricted to 0 <= x <= 1.
vars_una = mod.addVars(varNum, name="x", lb = 0, ub = 1, vtype = gurobipy.GRB.BINARY)
vars_pw = mod.addVars(index_pw, name="w", lb = 0, ub = 1, vtype = gurobipy.GRB.BINARY)
mod.update()
# Set Objective function.
obj = gurobipy.LinExpr(costs_una, vars_una.select())
index_nonzero = tuple(map(tuple, np.transpose(np.nonzero(costs_pw))))
for ind in index_nonzero:
obj += costs_pw[ind] * vars_pw[ind]
mod.setObjective(obj, gurobipy.GRB.MAXIMIZE)
# Set other constraints.
for j in range(len(conflicts)):
for k in range(varNum):
if k in conflicts[j]:
# Constraints multiplied by x_k:
lhs1 = gurobipy.LinExpr()
for i in np.delete(conflicts[j], np.argwhere(conflicts[j] == k)):
if i < k:
lhs1 += vars_pw[i,k]
else:
lhs1 += vars_pw[k,i]
mod.addConstr( lhs1 <= 0)
# Constraints multiplied by (1-x_k):
lhs2 = vars_una[k] - 1 - lhs1
for i in np.delete(conflicts[j], np.argwhere(conflicts[j] == k)):
lhs2 += vars_una[i]
mod.addConstr( lhs2 <= 0)
else:
# Constraints multiplied by x_k:
lhs1 = gurobipy.LinExpr()
for i in conflicts[j]:
if i < k:
lhs1 += vars_pw[i,k]
else:
lhs1 += vars_pw[k,i]
mod.addConstr( lhs1 - vars_una[k] <= 0)
# Constraints multiplied by (1-x_k):
lhs2 = vars_una[k] - 1 - lhs1
for i in conflicts[j]:
lhs2 += vars_una[i]
mod.addConstr( lhs2 <= 0)
# Check for labels that aren't present in any conflict set:
no_conf = sorted(set(range(varNum)) - set(np.unique(conflicts)))
# Set extra constraints, when problem is not in standard form:
if len(no_conf) > 0:
for i in no_conf:
for k in range(varNum):
if k < i:
mod.addConstr( vars_una[i] - vars_pw[k,i] >= 0 )
mod.addConstr( vars_una[k] - vars_pw[k,i] >= 0 )
mod.addConstr( vars_pw[k,i] - vars_una[i] - vars_una[k] + 1 >= 0 )
elif k > i:
mod.addConstr( vars_una[i] - vars_pw[i,k] >= 0 )
mod.addConstr( vars_una[k] - vars_pw[i,k] >= 0 )
mod.addConstr( vars_pw[i,k] - vars_una[i] - vars_una[k] + 1 >= 0 )
mod.update()
return mod
# Initializing Gurobi Model Q-MWIS6 with:
# - equality constraints.
# - Sherali-Adams linearization for problems in standard form.
# Function that takes in costs and conflict sets. Returns Gurobi model object.
# Input same as init_qmwis_inequality.
def init_qmwis_equality_sherali(costs_una, costs_pw, conflicts):
mod = gurobipy.Model()
# Set number of variables including and excluding slack variables for equality constraints.
varNum = len(costs_una) + len(conflicts)
varNum2 = len(costs_una)
# Define new cost vector for unary costs (slack variables have 0 cost).
n_costs_una = np.zeros(varNum)
n_costs_una[:costs_una.shape[0]] = costs_una
# Embed original pairwise cost matrix in new one with 0 cost for slack variables.
n_costs_pw = np.zeros((varNum, varNum))
n_costs_pw[:costs_pw.shape[0], :costs_pw.shape[1]] = costs_pw
# Add slack labels n+1, ..., n+m to the m conflict sets.
n_conflicts = np.append(conflicts, np.arange(len(costs_una), varNum, 1).reshape((len(conflicts),1)), axis = 1)
# Create index list of integer tuples (i,j) with i < j for pairwise terms w_ij.
index_pw = []
for i in range(varNum-1):
for j in range(i+1,varNum):
index_pw.append((i,j))
# Add unary and pairwise variables to the model. Variables set to be binary restricted to 0 <= x <= 1.
vars_una = mod.addVars(varNum, name="x", lb = 0, ub = 1, vtype = gurobipy.GRB.BINARY)
vars_pw = mod.addVars(index_pw, name="w", lb = 0, ub = 1, vtype = gurobipy.GRB.BINARY)
mod.update()
# Set Objective function.
obj = gurobipy.LinExpr(n_costs_una, vars_una.select())
index_nonzero = tuple(map(tuple, np.transpose(np.nonzero(n_costs_pw))))
for ind in index_nonzero:
obj += n_costs_pw[ind] * vars_pw[ind]
mod.setObjective(obj, gurobipy.GRB.MAXIMIZE)
# Set conflict set constraints.
for con in n_conflicts:
lhs = gurobipy.LinExpr()
for i in con:
lhs += vars_una[i]
mod.addConstr( lhs == 1 )
# Set other constraints.
for j in range(len(n_conflicts)):
for k in range(varNum2):
if k in n_conflicts[j]:
lhs = gurobipy.LinExpr()
for i in np.delete(n_conflicts[j], np.argwhere(n_conflicts[j] == k)):
if i < k:
lhs += vars_pw[i,k]
else:
lhs += vars_pw[k,i]
mod.addConstr( lhs == 0)
else:
lhs = gurobipy.LinExpr()
for i in n_conflicts[j]:
if i < k:
lhs += vars_pw[i,k]
else:
lhs += vars_pw[k,i]
mod.addConstr( lhs - vars_una[k] == 0)
# Set extra constraints, for problems containing negative costs.
if np.any(n_costs_una < 0) or np.any(n_costs_pw < 0):
for j in range(len(n_conflicts)):
for k in range(varNum2, varNum):
if k not in n_conflicts[j]:
lhs = gurobipy.LinExpr()
for i in n_conflicts[j]:
if i < k:
lhs += vars_pw[i,k]
else:
lhs += vars_pw[k,i]
mod.addConstr( lhs - vars_una[k] == 0)
# Check for labels that aren't present in any conflict set:
no_conf = sorted(set(range(varNum)) - set(np.unique(n_conflicts)))
# Set extra constraints, when problem is not in standard form:
if len(no_conf) > 0:
for i in no_conf:
for k in range(varNum2):
if k < i:
mod.addConstr( vars_una[i] - vars_pw[k,i] >= 0 )
mod.addConstr( vars_una[k] - vars_pw[k,i] >= 0 )
mod.addConstr( vars_pw[k,i] - vars_una[i] - vars_una[k] + 1 >= 0 )
elif i < k:
mod.addConstr( vars_una[i] - vars_pw[i,k] >= 0 )
mod.addConstr( vars_una[k] - vars_pw[i,k] >= 0 )
mod.addConstr( vars_pw[i,k] - vars_una[i] - vars_una[k] + 1 >= 0 )
mod.update()
return mod
# Knowledgedump.org - qmwis.py
# File to import/generate the data and run the solver.
import gurobipy
import numpy as np
import qmwis_algorithm as qmwis
import time
from itertools import repeat
import sys
import os
# Function to create example Q-MWIS problems.
# Input:
# -> varNum = n, number of decision variables.
# -> cost_range = [a,b], list of two integers defining cost interval; costs are integers here.
# -> conflict = [m,c], list of two integers defining number m and size c of conflict sets.
# -> seed = integer to fix random number generator seed for costs/conflict sets (optional).
# -> dens = density of pairwise cost matrix (default: 100%); If dens<1, sparse matrix is created.
# Note: This density factor only counts the values in the upper triangle of costs_pw,
# since lower triangle is unused (zero).
# Returns: Unary+pairwise costs and conflict sets of Q-MWIS problem, all stored as numpy arrays.
def create_example(varNum, cost_range, conflict, seed = 0.5, dens = 1):
# Set seed if parameter was changed to some integer.
if seed != 0.5:
np.random.seed(seed)
# costs_una is numpy array containing the varNum unary costs between costs[0] and costs[1].
costs_una = np.random.randint(cost_range[0], cost_range[1]+1, varNum)
# costs_pw is numpy array of shape varNum x varNum.
# Contains 1/2*varNum*(varNum-1) entries in the upper triangle, rest is 0.
# Pairwise cost values set randomly between costs[0] and costs[1].
# Initialize numpy array.
costs_pw = np.zeros((varNum, varNum))
# Fill upper triangle matrix with values. k=1 offset in order to not include diagonal.
pw_vals = np.random.randint(cost_range[0], cost_range[1]+1, int(varNum*(varNum-1)/2))
if dens == 1:
costs_pw[np.triu_indices(varNum, k=1)] = pw_vals
else:
# If density != 1, set a number of (1-dens)*len(pw_vals) values to 0 first.
zero_indices = np.random.choice(np.arange(len(pw_vals)), replace = False, size = int(len(pw_vals)*(1-dens)) )
pw_vals[zero_indices] = 0
costs_pw[np.triu_indices(varNum, k=1)] = pw_vals
# Conflict sets in the form of a numpy array.
# Each row corresponds to one conflict set, containing the indices of variables.
# Generate conflict[0] conflict sets, which contain conflict[1] unique elements each.
# Note: Conflict sets generated here may contain duplicate constraints.
conf = []
for i in repeat(None, conflict[0]):
conf.append(np.random.choice(varNum, size=conflict[1], replace=False))
conf = np.array(conf)
return [costs_una, costs_pw, conf]
# Creates different Gurobi models, solves them.
# Prints objective values, time taken to define and solve the model and the solution status (optimal/time-out).
# Input:
# -> inp = input list of unary costs, pairwise costs and conflict sets in the form of numpy arrays.
# -> time_lim = integer, setting time limit in seconds for solving each separate model. Default = 5*60 seconds.
# -> console = 0 (default) or 1, setting console output parameter of Gurobi solver.
# -> heur = 0 (default) or x time in (roughly) seconds, setting NoRelHeurWork time spent on heuristics without
# solving root relaxation.
# Returns:
# -> Gurobi model and the time taken to solve the problem.
# Included Gurobi model choices:
# -> 1: [] equality constr, [X] inequality constr, [] trivial linearization, [] Sherali-Adams linearization;
# -> 2: [X] equality constr, [] inequality constr, [] trivial linearization, [] Sherali-Adams linearization;
# -> 3: [] equality constr, [X] inequality constr, [X] trivial linearization, [] Sherali-Adams linearization;
# -> 4: [X] equality constr, [] inequality constr, [X] trivial linearization, [] Sherali-Adams linearization;
# -> 5: [] equality constr, [X] inequality constr, [] trivial linearization, [X] Sherali-Adams linearization;
# -> 6: [X] equality constr, [] inequality constr, [] trivial linearization, [X] Sherali-Adams linearization;
def solve_single(inp, model_num, time_lim = 5*60, console = 0, heur = 0):
# Set current time.
curr_time = time.time()
# Set Gurobi model, depending on model_num. Inputs inp = [costs_una, costs_pw, conf] are unpacked.
if model_num == 1:
mod = qmwis.init_qmwis_inequality(*inp)
elif model_num == 2:
mod = qmwis.init_qmwis_equality(*inp)
elif model_num == 3:
mod = qmwis.init_qmwis_inequality_triv(*inp)
elif model_num == 4:
mod = qmwis.init_qmwis_equality_triv(*inp)
elif model_num == 5:
mod = qmwis.init_qmwis_inequality_sherali(*inp)
elif model_num == 6:
mod = qmwis.init_qmwis_equality_sherali(*inp)
else:
raise ValueError("model_num has to be integer in interval [1,6].")
# Get time to set up the Gurobi model object.
time_setup = round(time.time() - curr_time, 2)
print("Model", model_num, " setup time:", time_setup, "seconds.")
# Set model parameters and call solver.
mod.Params.LogToConsole = console
mod.Params.TimeLimit = time_lim
# If heur > 0, set time spent on heuristics.
if heur > 0:
mod.Params.NoRelHeurWork = heur
mod.optimize()
# Get time elapsed.
time_taken = round(time.time() - curr_time, 2)
# Check optimization status; 2 = optimal solution found, 9 = solver timed out.
if mod.Status == 2:
status = "Optimal Solution found."
elif mod.Status == 9:
status = "Time-out."
else:
status = "ERROR."
print("Model", model_num, ": Obj:", round(mod.ObjVal, 2), "; Best bound:", round(mod.ObjBound, 2), "; Gap:", round(mod.MIPGap*100, 2),
"% ; Status:", status, "Time elapsed:", time_taken, "seconds.")
# Return model and the time it took to solve it (or until time-out).
return mod, time_taken
# Function to automatically build and solve all six Gurobi models for same inputs as in solve_single.
# The models are then ranked and compared according to solving speed or objective value (if optimal solution wasn't found yet).
# Ranking is printed to console.
def solve_all(inp, time_lim = 5*60, console = 0):
# Call solve_single for all models.
mod1, time_taken1 = solve_single(inp, 1, time_lim, console)
mod2, time_taken2 = solve_single(inp, 2, time_lim, console)
mod3, time_taken3 = solve_single(inp, 3, time_lim, console)
mod4, time_taken4 = solve_single(inp, 4, time_lim, console)
mod5, time_taken5 = solve_single(inp, 5, time_lim, console)
mod6, time_taken6 = solve_single(inp, 6, time_lim, console)
# Rank models according to solution speed / current best solution.
times_all = [time_taken1, time_taken2, time_taken3, time_taken4, time_taken5, time_taken6]
obj_all = [mod1.ObjVal, mod2.ObjVal, mod3.ObjVal, mod4.ObjVal, mod5.ObjVal, mod6.ObjVal]
status_all = [mod1.Status, mod2.Status, mod3.Status, mod4.Status, mod5.Status, mod6.Status]
gap_all = [mod1.MIPGap, mod2.MIPGap, mod3.MIPGap, mod4.MIPGap, mod5.MIPGap, mod6.MIPGap]
# Construct ranking:
ranking = []
# Get all model indices that were solved:
mods_solved = [ind for ind, stat in enumerate(status_all) if stat == 2]
# Add them to the ranking, ordered according to solution speed.
times_solved = [times_all[i] for i in mods_solved]
# Iterate over the sorted, unique elements of times_solved (conversion to set removes duplicates).
for i in sorted(list(set(times_solved))):
# Gets all indices in case two models take exactly the same amount of time to solve the problem.
indices_solved = [j+1 for j, x in enumerate(times_all) if x == i]
ranking.append(indices_solved)
# Get the model indices that timed out:
mods_timed = [ind for ind, stat in enumerate(status_all) if stat == 9]
# Add them to the ranking, ordered by current best objective value.
obj_timed = [obj_all[i] for i in mods_timed]
# Iterate over sorted, unique objective values (largest first).
for i in sorted(list(set(obj_timed)), reverse = True):
# If two timed out models have the same current best objective value, they're added simultaneously.
indices_timed = [j+1 for j, x in enumerate(obj_all) if x == i and status_all[j] == 9]
ranking.append(indices_timed)
print("Model ranking from fastest to slowest: ", ranking)
print("_______________________________\n")
mods = [mod1, mod2, mod3, mod4, mod5, mod6]
return mods, ranking, times_all, status_all, gap_all
if __name__ == "__main__":
print("Included Gurobi model setups:\n",
"-> 1: [] equality constr, [X] inequality constr, [] trivial linearization, [] Sherali-Adams linearization;\n",
"-> 2: [X] equality constr, [] inequality constr, [] trivial linearization, [] Sherali-Adams linearization;\n",
"-> 3: [] equality constr, [X] inequality constr, [X] trivial linearization, [] Sherali-Adams linearization;\n",
"-> 4: [X] equality constr, [] inequality constr, [X] trivial linearization, [] Sherali-Adams linearization;\n",
"-> 5: [] equality constr, [X] inequality constr, [] trivial linearization, [X] Sherali-Adams linearization;\n",
"-> 6: [X] equality constr, [] inequality constr, [] trivial linearization, [X] Sherali-Adams linearization;\n"
)
# Save path of python file directory.
current_dir_path = os.path.dirname(os.path.realpath(__file__))
# ************ PROBLEM SET 1 - Peformance of smaller problems with conflict sets of size 2 ************
# Write results to text-file.
sys.stdout = open(current_dir_path + "\\problem_set_1.txt", "w")
# Various smaller problems with n variables and n conflict sets of size 2, costs are integers in interval [1,10].
print("************ PROBLEM SET 1 - Peformance of smaller problems with conflict sets of size 2 ************")
for i in range(20,120,20):
print("_______________________________")
print("Problem size: %g variables and %g conflict sets of size 2. Results for Seeds 0 to 24, i.e. 25 Problems:" % (i,i))
print("_______________________________")
times_total = []
status_total = []
gap_total = []
for j in range(25):
# Save the time each model took, optimization status and MIP gap in a list of lists.
temp1, temp2, temp3 = solve_all(create_example(i, [1,10], [i,2], j))[2:5]
times_total.append(temp1)
temp2 = [1 if i == 2 else 0 for i in temp2]
status_total.append(temp2)
gap_total.append(temp3)
# Calculate average time and total number of models solved (within time limit), as well as avg MIP gap for each model.
times_avg = [round(sum(i)/len(times_total), 2) for i in zip(*times_total)]
gap_avg = [round(sum(i)/len(gap_total)*100, 2) for i in zip(*gap_total)]
status_sum = [sum(i) for i in zip(*status_total)]
print("Average completion time of the respective models for all problem instances (in seconds):")
print(times_avg)
print("Total number of problem instances solved by the respective models (within the time-limit of 5 minutes):")
print(status_sum)
print("Average MIP gap of the respective models for all problem instances (in %):")
print(gap_avg)
print("_______________________________")
sys.stdout.close()
#############################################################################################################################################
#############################################################################################################################################
# ************ PROBLEM SET 2 - Peformance of smaller problems with scaling conflict set size ************
# Write results to text-file.
sys.stdout = open(current_dir_path + "\\problem_set_2.txt", "w")
# Various smaller problems with n variables and n/4 conflict sets of size n/10, costs are integers in interval [1,3].
print("************ PROBLEM SET 2 - Peformance of smaller problems with scaling conflict set size ************")
for i in range(20,100,20):
print("_______________________________")
print("Problem size: %g variables and %g conflict sets of size %g. Results for Seeds 0 to 24, i.e. 25 Problems:" % (i,int(i/4),int(i/10)))
print("_______________________________")
times_total = []
status_total = []
gap_total = []
for j in range(25):
# Save the time, etc. each model took in a list of lists.
temp1, temp2, temp3 = solve_all(create_example(i, [1,3], [int(i/4),int(i/10)], j))[2:5]
times_total.append(temp1)
temp2 = [1 if i == 2 else 0 for i in temp2]
status_total.append(temp2)
gap_total.append(temp3)
# Calculate average time, etc. for each model.
times_avg = [round(sum(i)/len(times_total), 2) for i in zip(*times_total)]
status_sum = [sum(i) for i in zip(*status_total)]
gap_avg = [round(sum(i)/len(gap_total)*100, 2) for i in zip(*gap_total)]
print("Average completion time of the respective models for all problem instances (in seconds):")
print(times_avg)
print("Total number of problem instances solved by the respective models (within the time-limit of 5 minutes):")
print(status_sum)
print("Average MIP gap of the respective models for all problem instances (in %):")
print(gap_avg)
print("_______________________________")
sys.stdout.close()
#############################################################################################################################################
#############################################################################################################################################
# ************ PROBLEM SET 3 - Effect of conflict set size/number on mid-sized problems ************
# Write results to text-file.
sys.stdout = open(current_dir_path + "\\problem_set_3.txt", "w")
# n = 150, m = 150, |K_j| = 2, c \in [1,10]
# n = 150, m = 150, |K_j| = 4, c \in [1,10]
# n = 150, m = 150, |K_j| = 10, c \in [1,10]
# n = 150, m = 50, |K_j| = 2, c \in [1,10]
# n = 150, m = 50, |K_j| = 10, c \in [1,10]
# n = 150, m = 50, |K_j| = 20, c \in [1,10]
# n = 150, m = 10, |K_j| = 20, c \in [1,10]
# n = 150, m = 10, |K_j| = 30, c \in [1,10]
# n = 150, m = 10, |K_j| = 50, c \in [1,10]
print("************ PROBLEM SET 3 - Effect of conflict set size/number on mid-sized problems ************")
for m, kj in [[150,2],[150,4],[150,10],[50,2],[50,10],[50,20],[10,20],[10,30],[10,50]]:
print("_______________________________")
print("Problem size: %g variables and %g conflict sets of size %g. Results for Seeds 0 to 2, i.e. 3 Problems:" % (150,m,kj))
print("_______________________________")
times_total = []
status_total = []
gap_total = []
for j in range(3):
# Save the time, etc. each model took in a list of lists.
temp1, temp2, temp3 = solve_all(create_example(150, [1,10], [m,kj], j))[2:5]
times_total.append(temp1)
temp2 = [1 if i == 2 else 0 for i in temp2]
status_total.append(temp2)
gap_total.append(temp3)
# Calculate average time, etc. for each model.
times_avg = [round(sum(i)/len(times_total), 2) for i in zip(*times_total)]
status_sum = [sum(i) for i in zip(*status_total)]
gap_avg = [round(sum(i)/len(gap_total)*100, 2) for i in zip(*gap_total)]
print("Average completion time of the respective models for all problem instances (in seconds):")
print(times_avg)
print("Total number of problem instances solved by the respective models (within the time-limit of 5 minutes):")
print(status_sum)
print("Average MIP gap of the respective models for all problem instances (in %):")
print(gap_avg)
print("_______________________________")
sys.stdout.close()
#############################################################################################################################################
#############################################################################################################################################
# ************ PROBLEM SET 4 - Effect of cost range on mid-sized problems ************
# Write results to text-file.
sys.stdout = open(current_dir_path + "\\problem_set_4.txt", "w")
# n = 150, m = 50, |K_j| = 5, c \in [1,10]
# n = 150, m = 50, |K_j| = 5, c \in [-10,10]
# n = 150, m = 50, |K_j| = 5, c \in [-100,100]
# n = 150, m = 50, |K_j| = 5, c \in [1,1000]
print("************ PROBLEM SET 4 - Effect of cost range on mid-sized problems ************")
for c1,c2 in [[1,10],[-10,10],[-100,100],[1,1000]]:
print("_______________________________")
print("Problem size: %g variables and %g conflict sets of size %g. Results for Seeds 0 to 9, i.e. 10 Problems:" % (150,50,5))
print("_______________________________")
times_total = []
status_total = []
gap_total = []
for j in range(10):
# Save the time, etc. each model took in a list of lists.
temp1, temp2, temp3 = solve_all(create_example(150, [c1,c2], [50,5], j))[2:5]
times_total.append(temp1)
temp2 = [1 if i == 2 else 0 for i in temp2]
status_total.append(temp2)
gap_total.append(temp3)
# Calculate average time, etc. for each model.
times_avg = [round(sum(i)/len(times_total), 2) for i in zip(*times_total)]
status_sum = [sum(i) for i in zip(*status_total)]
gap_avg = [round(sum(i)/len(gap_total)*100, 2) for i in zip(*gap_total)]
print("Average completion time of the respective models for all problem instances (in seconds):")
print(times_avg)
print("Total number of problem instances solved by the respective models (within the time-limit of 5 minutes):")
print(status_sum)
print("Average MIP gap of the respective models for all problem instances (in %):")
print(gap_avg)
print("_______________________________")
sys.stdout.close()
#############################################################################################################################################
#############################################################################################################################################
# ************ PROBLEM SET 5 - large problem instances, Sherali with heuristic ************
# Write results to text-file.
sys.stdout = open(current_dir_path + "\\problem_set_5.txt", "w")
# n = 1500, m = 100, |K_j| = 50, c \in [1,100]
# n = 1500, m = 100, |K_j| = 50, c \in [-100,100]
print("************ PROBLEM SET 5 ************")
for c1,c2 in [[1,100],[-100,100]]:
print("_______________________________")
print("Problem size: %g variables and %g conflict sets of size %g. Results for Seeds 0 to 2, i.e. 3 Problems:" % (1500,100,50))
print("_______________________________")
gap_total = []
for j in range(3):
inp = create_example(1500, [c1,c2], [100,50], j)
mod1 = solve_single(inp, model_num = 1, time_lim = 15*60, console = 1)[0]
mod1h = solve_single(inp, model_num = 1, time_lim = 15*60, console = 1, heur = 15*60)[0]
mod6 = solve_single(inp, model_num = 6, time_lim = 15*60, console = 1, heur = 15*60)[0]
gap_total.append([mod1.MIPGap, mod1h.MIPGap, mod6.MIPGap])
# Calculate average MIP gap for each model.
gap_avg = [round(sum(i)/len(gap_total)*100, 2) for i in zip(*gap_total)]
print("Average MIP gap of the respective models 1 (without/with heuristics) and 6 for all problem instances (in %):")
print(gap_avg)
print("_______________________________")
sys.stdout.close()
# Supplementary: Checking, how Q-MWIS1 and 6 perform, when time limit is increased further:
sys.stdout = open(current_dir_path + "\\problem_set_5_single_trial_long.txt", "w")
solve_single(inp, model_num = 1, time_lim = 60*60, console = 1)[0]
solve_single(inp, model_num = 6, time_lim = 60*60, console = 1, heur = 60*60)[0]
sys.stdout.close()
#############################################################################################################################################
#############################################################################################################################################
# ************ PROBLEM SET 6 - Effect of sparsity on pairwise cost matrix ************
# Write results to text-file.
sys.stdout = open(current_dir_path + "\\problem_set_6.txt", "w")
# n = 200, m = 100, |K_j| = 10, c \in [1,20], density = 100%
# n = 200, m = 100, |K_j| = 10, c \in [1,20], density = 50%
# n = 200, m = 100, |K_j| = 10, c \in [1,20], density = 20%
# n = 200, m = 100, |K_j| = 10, c \in [1,20], density = 10%
# n = 200, m = 100, |K_j| = 10, c \in [1,20], density = 5%
# n = 200, m = 100, |K_j| = 10, c \in [1,20], density = 1%
print("************ PROBLEM SET 6 - Effect of sparsity on pairwise cost matrix ************")
for d in [1, 0.5, 0.2, 0.1, 0.05, 0.01]:
print("_______________________________")
print("Problem size: %g variables, %g conflict sets of size %g and pairwise cost density of %g%%. Results for Seeds 0 to 2, i.e. 3 Problems:" % (200,50,20,d*100))
print("_______________________________")
times_total = []
status_total = []
gap_total = []
for j in range(3):
# Save the time, etc. each model took in a list of lists.
temp1, temp2, temp3 = solve_all(create_example(200, [1,20], [50,20], j, dens = d))[2:5]
times_total.append(temp1)
temp2 = [1 if i == 2 else 0 for i in temp2]
status_total.append(temp2)
gap_total.append(temp3)
# Calculate average time, etc. for each model.
times_avg = [round(sum(i)/len(times_total), 2) for i in zip(*times_total)]
status_sum = [sum(i) for i in zip(*status_total)]
gap_avg = [round(sum(i)/len(gap_total)*100, 2) for i in zip(*gap_total)]
print("Average completion time of the respective models for all problem instances (in seconds):")
print(times_avg)
print("Total number of problem instances solved by the respective models (within the time-limit of 5 minutes):")
print(status_sum)
print("Average MIP gap of the respective models for all problem instances (in %):")
print(gap_avg)
print("_______________________________")
sys.stdout.close()