Automatized model tuning
Contents
Combinatorial Optimization course, FEE CTU in Prague. Created by Industrial Informatics Department.
%pip install -i https://pypi.gurobi.com gurobipy
import gurobipy as g
import matplotlib.pyplot as plt
import numpy as np
import networkx as nx
Looking in indexes: https://pypi.gurobi.com
Collecting gurobipy
Downloading gurobipy-9.5.1-cp310-cp310-manylinux2014_x86_64.whl (11.4 MB)
?25l ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/11.4 MB ? eta -:--:--
━━━━╺━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 1.2/11.4 MB 36.5 MB/s eta 0:00:01
━━━━━━━━━━━━━━━━━━━━━━━╸━━━━━━━━━━━━━━━━ 6.8/11.4 MB 100.4 MB/s eta 0:00:01
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 11.4/11.4 MB 174.4 MB/s eta 0:00:01
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 11.4/11.4 MB 135.3 MB/s eta 0:00:00
?25h
Installing collected packages: gurobipy
Successfully installed gurobipy-9.5.1
WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv
Note: you may need to restart the kernel to use updated packages.
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Input In [1], in <cell line: 3>()
1 get_ipython().run_line_magic('pip', 'install -i https://pypi.gurobi.com gurobipy')
2 import gurobipy as g
----> 3 import matplotlib.pyplot as plt
4 import numpy as np
5 import networkx as nx
ModuleNotFoundError: No module named 'matplotlib'
Automatized model tuning#
As an example, we will revisit the model for Game of Fiver which for larger board sizes took considerable time to solve:
def game_of_fivers(n, params=None):
m = g.Model()
x = m.addVars(n+2, n+2, vtype=g.GRB.BINARY, obj=1)
k = m.addVars(range(1,n+1), range(1,n+1), vtype=g.GRB.INTEGER)
for i in range(1,n+1):
for j in range(1,n+1):
m.addConstr(x[i,j] + x[i+1, j] + x[i-1,j] + x[i,j+1] + x[i,j-1] == 2*k[i,j] + 1)
m.addConstr(x.sum(0,"*") + x.sum(n+1,"*") + x.sum("*",0) + x.sum("*",n+1) == 0)
m.write('fivers_n{}.lp'.format(n))
if params is not None:
m.params.CutPasses = 10
m.params.PreDual = 1
m.params.outputflag = 0
m.optimize()
X = [[int(round(x[i,j].X)) for j in range(1,n+1)] for i in range(1,n+1)]
return m.runtime
# run with default parameters
def_time = game_of_fivers(21)
print('Time with default settings: {}s'.format(def_time))
# lets try to tune some of the params
model = g.read('fivers_n21.lp')
model.params.tuneResults = 1
model.params.TuneTimeLimit = 30 # how much time to invest into the tuning
model.tune()
if model.tuneResultCount > 0:
model.getTuneResult(0)
model.write('fivers_tuned_params.prm')
print('Time before tuning: {}s'.format(def_time))
print('Time after tuning: {}s'.format(game_of_fivers(21, {'CutPasses': 10, 'PreDual':1})))
Nonlinear constrains#
n = 200
t = np.linspace(0, 20, n)
y = 3*np.sin(t)+np.cos(6*t)+0.5*t+3
plt.plot(t, y)
m = g.Model()
u = m.addVar(vtype=g.GRB.CONTINUOUS)
v = m.addVar(vtype=g.GRB.CONTINUOUS)
m.addGenConstrPWL(u, v, t, y)
m.setObjective(v)
m.optimize()
plt.plot(t, y)
plt.plot(u.x, v.x, marker='o', markersize=8, color="red")
Solution pool#
dg = nx.DiGraph()
dg = nx.random_k_out_graph(40, 5, 0.4, seed=22)
pos = nx.spring_layout(dg)
nx.draw(dg, pos, node_color='k', node_size=3, edge_color='grey', with_labels=True)
s = 28
t = 4
E = dg.edges()
E = dict.fromkeys(E)
V = dg.nodes()
np.random.seed(69)
w = np.random.randint(0, 50, len(E))
m = g.Model()
x = m.addVars(E.keys(), vtype=g.GRB.BINARY, ub=1, obj=w)
m.addConstr(g.quicksum([x[s, j] for k, j in E if k == s]) == 1)
m.addConstr(g.quicksum([x[i, t] for i, k in E if k == t]) == 1)
for i in V:
if i not in [s, t]:
m.addConstr(g.quicksum([x[i, j] for k, j in E if k == i]) == g.quicksum([x[j, i] for j, k in E if k == i]))
m.setParam(g.GRB.Param.PoolSolutions, 3)
m.setParam(g.GRB.Param.PoolSearchMode, 2) # k-best solutions
m.optimize()
sols = [0]*m.solcount
colorlist = ['r', 'g', 'b']
print('Found {} solutions.'.format(m.solcount))
for sol_idx in range(m.solcount):
print('Sol no. {}'.format(sol_idx+1))
sols[sol_idx] = []
m.setParam(g.GRB.Param.SolutionNumber, sol_idx)
for i, j in E:
if x[i, j].xn > 0.5:
print(i, j)
sols[sol_idx] += [(i, j)]
nx.draw(dg, pos, node_color='k', node_size=3, edge_color='grey')
for k in range(m.solcount):
nx.draw_networkx_edges(dg, pos, edgelist=sols[k], edge_color=colorlist[k], width=3)
Bratley revisited: CP model#
# Visualization
import matplotlib.pyplot as plt
def plot_solution(s, p):
"""
s: solution vector
p: processing times
"""
fig = plt.figure(figsize=(10,2))
ax = plt.gca()
ax.set_xlabel('time')
ax.grid(True)
ax.set_yticks([2.5])
ax.set_yticklabels(["oven"])
eps = 0.25 # just to show spaces between the dishes
ax.broken_barh([(s[i], p[i]-eps) for i in range(len(s))], (0, 5),
facecolors=('tab:orange', 'tab:green', 'tab:red', 'tab:blue', 'tab:gray'))
# Load data
path = "./bratley_data/instances/public_3.txt"
with open(path, "r") as f_in:
lines = f_in.readlines()
n = int(lines[0].strip())
r,d,p = [],[],[]
for i in range(n):
(pi, ri, di) = list(map(int, lines[1+i].split()))
r.append(ri)
p.append(pi)
d.append(di)
print("r", r, "d", d, "p", p, sep="\n")
# Generate data
from numpy import random as rnd
import numpy as np
rnd.seed(15)
n = 200
p = [rnd.randint(1,100) for i in range(n)]
r = [0 for i in range(n)]
for i in range(1,n):
r[i] = int(round(r[i-1] + rnd.exponential(0.5* sum(p)/len(p))))
d = [ int(round(r[i] + p[i] + rnd.exponential(100*sum(p)/len(p)) )) for i in range(n)]
print("r", r, "d", d, "p", p, sep="\n")
CP model#
from docplex.cp.model import CpoModel
from docplex.cp.config import context
import sys
# Create model
m = CpoModel()
# - add variables
tasks = [m.interval_var(name="task{:d}".format(i), optional=False, size=p[i]) for i in range(n)]
seq = m.sequence_var(tasks, name='seq')
# - set objective
m.add(m.minimize(m.max([m.end_of(tasks[i]) for i in range(n)]) ) ) # minimize C_max
# - add constraints
for i in range(n):
m.add(m.start_of(tasks[i]) >= r[i]) # release time
m.add(m.end_of(tasks[i]) <= d[i]) # deadline
m.add(m.no_overlap(seq)) # one task executed at one time
# Solve the model
msol = m.solve(TimeLimit=10, LogVerbosity="Normal", LogPeriod=1, Workers=1)
# Print the solution
print()
if msol.is_solution():
starts = [msol.get_value(tasks[i])[0] for i in range(n)]
print(*starts, sep="\n")
else:
print("No solution found.")
print("Done")
plot_solution(starts, p)
ILP model#
import gurobipy as g
m = g.Model()
# - add variables
s = m.addVars(n, vtype=g.GRB.CONTINUOUS, lb=0)
x = {}
for i in range(n):
for j in range(i + 1, n):
x[i, j] = m.addVar(vtype=g.GRB.BINARY)
Cmax = m.addVar(vtype=g.GRB.CONTINUOUS, obj=1)
# - add constraints
for i in range(n):
m.addConstr(s[i] + p[i] <= Cmax)
m.addConstr(s[i] >= r[i])
m.addConstr(s[i] + p[i] <= d[i])
M = max(d)
for i in range(n):
for j in range(i + 1, n):
m.addConstr(s[i] + p[i] <= s[j] + M*(1-x[i, j]))
m.addConstr(s[j] + p[j] <= s[i] + M*x[i, j])
# call the solver -----------------------------------------------
m.optimize()
print()
if m.SolCount > 0:
starts = [s[i].X for i in range(n)]
else:
print("No solution was found.")
print("Done")
plot_solution(starts, p)
TSP#
import math
from collections import namedtuple
import gurobipy as g
Point = namedtuple("Point", ['x', 'y'])
def length(point1, point2):
return int(round(math.sqrt((point1.x - point2.x)**2 + (point1.y - point2.y)**2))) # CP works with int only
class TSP:
def __init__(self):
D = None # distance matrix
points = None # vertices
def load_instance(self, path):
input_data_file = open(path, 'r')
input_data = ''.join(input_data_file.readlines())
# parse the input
lines = input_data.split('\n')
nodeCount = int(lines[0])
points = []
for i in range(1, nodeCount+1):
parts = lines[i].split()
points.append(Point(float(parts[0]), float(parts[1])))
# distance matrix
D = [[0 for _ in range(nodeCount+1)] for _ in range(nodeCount+1)] # Add dummy vertex last
for i in range(nodeCount):
for j in range(nodeCount):
D[i][j] = length(points[i], points[j])
# the last vertex is the same as the first one
for i in range(nodeCount):
D[i][-1] = D[i][0]
D[-1][i] = D[0][i]
self.D = D
self.points = points
return self
instances = [
{"inst": TSP().load_instance("./tsp_data/tsp_5_1"),
"init": [0, 1, 2, 4, 3]},
{"inst": TSP().load_instance("./tsp_data/tsp_51_1"),
"init": [0, 5, 2, 28, 10, 9, 45, 3, 46, 8, 4, 35, 13, 7, 19, 40, 18, 11, 42, 37, 20, 25, 1, 31, 22, 48, 32, 17, 49, 39, 50, 38, 15, 44, 14, 16, 29, 43, 21, 30, 12, 23, 34, 24, 41, 27, 36, 6, 26, 47, 33]},
{"inst": TSP().load_instance("./tsp_data/tsp_70_1"),
"init": [0, 35, 50, 11, 57, 2, 56, 27, 21, 49, 58, 53, 41, 36, 38, 52, 6, 5, 7, 51, 55, 68, 46, 67, 24, 16, 44, 39, 22, 1, 14, 15, 20, 29, 28, 45, 12, 31, 18, 26, 3, 59, 9, 25, 4, 10, 61, 43, 32, 8, 64, 54, 48, 62, 13, 19, 60, 42, 37, 66, 40, 17, 30, 23, 69, 33, 65, 34, 47, 63]},
{"inst": TSP().load_instance("./tsp_data/tsp_100_1"),
"init": [0, 6, 69, 61, 76, 35, 84, 11, 9, 26, 72, 47, 40, 94, 81, 60, 64, 66, 8, 23, 70, 59, 33, 67, 43, 37, 65, 71, 19, 15, 75, 14, 53, 46, 5, 29, 80, 38, 91, 57, 41, 50, 12, 55, 98, 39, 24, 68, 2, 28, 73, 87, 48, 85, 21, 96, 42, 77, 16, 7, 10, 74, 30, 18, 17, 34, 22, 99, 93, 51, 3, 89, 13, 31, 44, 62, 25, 82, 86, 54, 1, 27, 45, 88, 79, 97, 49, 90, 20, 63, 52, 92, 95, 78, 83, 32, 4, 56, 58, 36]},
{"inst": TSP().load_instance("./tsp_data/tsp_200_1"),
"init": [0, 103, 62, 192, 5, 48, 89, 148, 117, 9, 128, 83, 136, 23, 37, 108, 177, 181, 98, 106, 35, 160, 125, 131, 123, 58, 73, 20, 145, 71, 111, 46, 97, 22, 114, 112, 178, 59, 61, 163, 119, 154, 141, 34, 85, 26, 11, 19, 146, 130, 166, 76, 164, 179, 60, 24, 80, 101, 134, 68, 167, 129, 188, 158, 102, 172, 88, 168, 41, 30, 79, 55, 199, 132, 144, 96, 180, 196, 3, 64, 65, 195, 25, 186, 151, 110, 183, 147, 69, 21, 15, 87, 143, 162, 93, 150, 115, 17, 78, 52, 165, 18, 191, 198, 118, 109, 74, 135, 156, 173, 7, 113, 91, 159, 57, 176, 50, 86, 56, 6, 8, 105, 153, 174, 82, 54, 107, 121, 33, 28, 45, 116, 124, 133, 189, 42, 2, 13, 197, 157, 40, 70, 99, 187, 47, 127, 138, 137, 170, 29, 171, 182, 161, 84, 67, 72, 122, 49, 43, 169, 175, 190, 193, 194, 149, 38, 185, 95, 155, 51, 77, 104, 4, 142, 36, 32, 75, 12, 94, 81, 1, 63, 39, 120, 53, 140, 66, 27, 92, 126, 90, 44, 184, 31, 100, 152, 14, 16, 10, 139]}
]
inst_id = 4
inst = instances[inst_id]["inst"]
init_order = instances[inst_id]["init"]
INITIALIZE = False
ILP#
nodeCount = len(inst.points)
points = inst.points
# Create model
m = g.Model("tsp")
# - add variables
x = m.addVars(nodeCount,nodeCount, vtype=g.GRB.BINARY, name="x")
u = m.addVars(nodeCount, vtype=g.GRB.INTEGER, lb=0, name="u")
# - set objective
obj = g.quicksum(g.quicksum(inst.D[i][j]*x[i,j] for j in range(nodeCount)) for i in range(nodeCount))
m.setObjective(obj, g.GRB.MINIMIZE)
# - add constraints
m.addConstrs((1 == g.quicksum(x[i,j] for j in range(nodeCount)) for i in range(nodeCount)))
m.addConstrs((1 == g.quicksum(x[j,i] for j in range(nodeCount)) for i in range(nodeCount)))
for i in range(1, nodeCount):
for j in range(1, nodeCount):
m.addConstr(u[i]-u[j]+1 <= nodeCount*(1-x[i,j]))
# Initialization
if INITIALIZE:
for i, order in enumerate(init_order):
u[order].start = i
m.Params.TimeLimit = 10
m.optimize()
# Print the solution
print()
if m.SolCount > 0:
obj = m.objVal
print("Objective {}".format(obj))
order = [u[i].X for i in range(nodeCount)]
indices = range(nodeCount)
s = sorted(zip(order,indices), key=lambda x: x[0])
print([x[1] for x in s])
else:
print("No solution was found.")
print("Done")
CP#
from docplex.cp.model import CpoModel
from docplex.cp.config import context
import sys
node_count = len(inst.points)
points = inst.points
# Create model
m = CpoModel()
# - add variables
cities = [m.interval_var(name="city{:d}".format(i), optional=False, size=1) for i in range(node_count+1)]
seq = m.sequence_var(cities, name='seq', types=([i for i in range(node_count)] + [0]))
# - set objective
m.add(m.minimize(m.max([m.end_of(cities[i]) for i in range(len(cities))]) - len(cities)))
# - add constraints
m.add(m.first(seq, cities[0])) # start from city 0
m.add(m.last(seq, cities[-1])) # repeat the same city last
m.add(m.no_overlap(seq, inst.D, True))
# Solve the model
msol = m.solve(TimeLimit=10, LogVerbosity="Verbose", LogPeriod=1, Workers=1)
# Print the solution
print()
if msol.is_solution():
ovals = msol.get_objective_values()
print("Objective {}".format(ovals[0]))
starts = [msol.get_value(cities[i])[0] for i in range(len(cities)-1)]
indices = range(len(cities)-1)
s = sorted(zip(starts,indices), key=lambda x: x[0])
print([x[1] for x in s])
else:
print("No solution found.")
print("Done")
Comparison#
Best objective found by ILP and CP model under 10s timelimit. (without initialization)
instance |
ILP |
CP |
|---|---|---|
5 |
3 |
3 |
50 |
496 |
441 |
70 |
994 |
710 |
100 |
- |
22247 |
200 |
- |
31962 |