Error occurred while defining a callback function in class
Answered
Hi,
When I implement Benders decomposition method with gurobipy, some bugs occur.
This 1st - 4th pictures show 1) the whole codes, 2)code implementing optimization with user-defined cut (EPI, callback function ) , 3) the details of callback function, and 4) the details of errors, respectively.
0
-
Official comment
This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum. Or why not try our AI Gurobot?. -
Gurobi expects that the callback function passed to Model.optimize() is a user-defined function, not a method. The callback function should take two arguments: a model and a where value. Can you try restructuring your code so \( \texttt{EPI} \) is a user-defined function of this form? For example:
def EPI(model, where):
...Then, invoke the callback function with:
self.model.optimize(EPI)
0 -
Hi Eli,
I have revised my code according to what you said, however, there is a new error occurred.
NameError: name 'EPI' is not defined
My code is uploaded as the attached file.
# Import Gurobi Library
# import csv
import gurobipy as gb
import numpy as np
from scipy.io import loadmat
import time
import sys
from csv import writer
# Class which can have attributes set.
class expando(object):
pass
# Master problem
class Benders_Master:
def __init__(self, benders_gap, run, node, vehicle, delta_val,gamma):
self.max_iters = 1e4
self.data = expando()
self.variables = expando()
self.constraints = expando()
self.results = expando()
self._load_data(benders_gap, run, node, vehicle, delta_val,gamma)
self.data.run = run
self.data.node = node
self.data.vehicle = vehicle
self.data.delta_val = delta_val
self.data.z_fea_sub = []
self._build_model()
def EPI(model, where):
# model = self.model
GRB = gb.GRB
if where == GRB.Callback.MIPNODE: #
nodecnt = int(model.cbGet(GRB.callback.MIPNODE_NODCNT))
nodecnt = nodecnt / model._numcheckedNode
status = model.cbGet(GRB.Callback.MIPNODE_STATUS) #
# print(nodecnt)
if status == GRB.OPTIMAL:
if (nodecnt.is_integer()) & (model._numEPI <= model._numEPImax): # model._numEPImax
# if nodecnt <= 5e3:
sol = model.cbGetNodeRel(model._vars) # A relaxed sol. obtained at the explored node
SubFun = model._Gfun
N = model._coeff[0]
K = model._coeff[1]
cut = np.zeros((K, N ** 2))
Z = np.zeros(K)
label_cut = np.zeros((K, N ** 2))
cut_reindex = np.zeros((K, N ** 2))
pi = np.zeros((K, N ** 2))
for i in range(K):
# for j in range(N ** 2):
cut[i, :] = sol[i * (N ** 2): (i + 1) * (N ** 2)]
for i in range(K):
Z[i] = sol[i + K * (N ** 2)]
# Indeices of relaxed sol in descending order (max-min).
for i in range(K):
label_cut[i, :] = np.argsort(-cut[i, :])
# label_cut.astype(int)
for i in range(K):
for j in range(N ** 2):
pi[i, j] = SubFun[int(label_cut[i, j])]
cut_reindex[i, j] = cut[i, int(label_cut[i, j])]
# EPIs construction
# for i in range(K):
# if np.dot(pi[i, :], cut[i, :]) > Z[i]:
# model.cbCut( sum((pi[i, j] * model._vars[j + i * (N ** 2)]) for j in range(N ** 2)) <= model._vars[ i + K * (N ** 2)])
# model._numEPI += 1
if sum(np.dot(pi[i, :], cut_reindex[i, :]) for i in range(K)) > sum((Z[i]) for i in range(K)):
model.cbCut(
sum(sum((pi[i, j] * model._vars[int(label_cut[i, j]) + i * (N ** 2)]) for i in range(K))
for j in range(N ** 2)) <= sum((model._vars[k + K * (N ** 2)]) for k in range(K)))
model._numEPI += 1
def optimize(self, simple_results=False):
data = self.data
N = data.N
K = data.K
T = data.T
c = data.c
maxEPI = 1e4
# tlimit = 1e4
# Initial solution
print('\n' + '#' * 50)
print('Master problem optimization')
# self.model.optimize() # Optimize MP
# Submodular set function (quadratic form)
GFun = []
for i in range(N):
for j in range(N):
GFun.append(T[i, j] + c[i])
self.model._Gfun = GFun
self.model._coeff = [N, K]
self.model._numEPI = 0;
self.model._numEPImax = maxEPI;
self.model._numcheckedNode = 1e0;
self.model._num_inte_sol = 0;
# self.model.setParam('TimeLimit', tlimit)
self.model.setParam('Heuristics', 0.00)
self.model.setParam('Threads', 1)
self.model.setParam('PreCrush', 1)
# self.model.optimize() # Optimize MP with EPIs
self.model.optimize(EPI) # Optimize MP with EPIs
# self.model.optimize(self.EPI) # Optimize MP with EPIs
self.submodel = Benders_Subproblem(self,data.run,data.node , data.vehicle,data.delta_val,data.gamma) # Build subproblem from solution
self.submodel.update_fixed_vars(self)
print('\n' + '#' * 50)
print('Subproblem optimization')
self.submodel.optimize()
self._add_cut()
self._update_bounds()
self._save_vars()
while (self.data.ub - self.data.lb)/self.data.lb > self.data.benders_gap and len(self.data.cutlist) < self.max_iters:
print('\n' + '#' * 50)
print('Master problem optimization_{}'.format(len(self.data.cutlist)) )
self.model.optimize() # Optimize master problem
self.submodel.update_fixed_vars(self) # Update the decision var. of MP
print('\n' + '#' * 50)
print('Subproblem optimization_{}'.format(len(self.data.cutlist)))
self.submodel.optimize() # Optimize subproblem
self._add_cut() # Add Generalized Benders cut
self._update_bounds()
self._save_vars()
print('\n' + '#' * 50)
print('The # of iteration = {}'.format(len(self.data.cutlist)))
print('LB, UB = {},{}'.format(self.data.lb, self.data.ub))
pass
print('\n' + '#' * 50)
print('The # of iteration = {}'.format(len(self.data.cutlist)))
print('LB, UB = {},{}'.format(self.data.lb, self.data.ub))
print('# of added cut={}'.format(len(self.data.cutlist)))
# print('Solution of x = {}'.format(self.data.xs))
# print('Solution of theta = {}'.format(self.data.thetas))
###
# Loading functions
###
def _load_data(self, benders_gap, run, node, vehicle, delta_val,gamma):
self.data.cutlist = []
self.data.upper_bounds = []
self.data.lower_bounds = []
self.data.lambdas = {}
self.data.benders_gap = benders_gap
self.data.gamma = gamma
self.data.ub = gb.GRB.INFINITY
self.data.lb = -gb.GRB.INFINITY
self.data.xs = []
self.data.ys = []
self.data.thetas = []
data = loadmat('DataMap/Vehicles_{}/RealMap_{}_test.mat'.format(run, node))
self.data.Emax = data['Emax'][0][0]
self.data.rmax = self.data.Emax
self.data.e = data['e']
self.data.T = 1e0 * data['T']
self.data.M = int(data['M'])
self.data.p = data['p'][0]
self.data.g = data['g'][0]
self.data.t = 1e0 * data['t'][0]
self.data.c = data['c'][0]
self.data.K = int(vehicle)
self.data.N = int(data['N'])
self.data.alpha = 0.05
self.data.Q = 200
self.data.deltabar = np.array(delta_val)
###
# Model Building
###
def _build_model(self):
self.model = gb.Model()
self._build_variables()
self._build_objective()
self._build_constraints()
self.model.update()
def _build_variables(self):
m = self.model
K = self.data.K
N = self.data.N
self.variables.theta = m.addVar(lb=0, ub=gb.GRB.INFINITY, name='theta')
self.variables.x = m.addVars(K, N, N, vtype=gb.GRB.BINARY, name="x")
self.variables.omega1 = m.addVars(N, vtype=gb.GRB.BINARY, name="omega1")
self.variables.omega2 = m.addVars(N, vtype=gb.GRB.BINARY, name="omega2")
self.variables.omega3 = m.addVars(N, vtype=gb.GRB.BINARY, name="omega3")
self.variables.omega4 = m.addVars(N, vtype=gb.GRB.BINARY, name="omega4")
if len(self.data.cutlist) ==0:
m = self.model
# K = self.data.K
N = self.data.N
# deltabar = self.data.deltabar
self.variables.eta = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="eta")
self.variables.epsilon = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="epsilon")
self.variables.q = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="q")
self.variables.u = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="u")
self.variables.v = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="v")
self.variables.zeta1 = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="zeta1")
self.variables.zeta2 = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="zeta2")
self.variables.delta = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="delta")
self.variables.tau = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="tau")
m.update()
def _build_objective(self):
K = self.data.K
N = self.data.N
c = self.data.c
T = self.data.T
qq = 0
for k in range(K):
for i in range(N):
for j in range(N):
qq += self.variables.x[k, i, j] * (c[i] + T[i, j])
if len(self.data.cutlist) ==0:
self.model.setObjective(self.variables.theta + qq, gb.GRB.MINIMIZE)
else:
self.model.setObjective(self.variables.theta + qq, gb.GRB.MINIMIZE)
def _build_constraints(self):
m = self.model
K = self.data.K
N = self.data.N
self.constraints.c1 = m.addConstrs(self.variables.x.sum('*', i, '*') <= 1 for i in range(N) if i != 0)
self.constraints.c2 = m.addConstrs(self.variables.x[k, i, i] == 0 for i in range(N) for k in range(K))
self.constraints.c3 = m.addConstrs(self.variables.x[k, N - 1, i] == 0 for i in range(N) for k in range(K))
self.constraints.c4 = m.addConstrs(self.variables.x[k, i, 0] == 0 for i in range(N) for k in range(K))
# Flow conservation
self.constraints.c5 = m.addConstrs(
self.variables.x.sum(k, 0, '*') - self.variables.x.sum(k, '*', 0) == 1 for k in range(K))
self.constraints.c6 = m.addConstrs(
self.variables.x.sum(k, N - 1, '*') - self.variables.x.sum(k, '*', N - 1) == -1 for k in range(K))
self.constraints.c7 = m.addConstrs(
self.variables.x.sum(k, i, '*') - self.variables.x.sum(k, '*', i) == 0 for i in range(N - 1) if i != 0 for k
in range(K))
if len(self.data.cutlist) ==0:
M = self.data.M
N = self.data.N
t = self.data.t
T = self.data.T
deltabar = self.data.deltabar
gamma = self.data.gamma
x = self.variables.x
omega1 = self.variables.omega1
omega2 = self.variables.omega2
omega3 = self.variables.omega3
omega4 = self.variables.omega4
tau = self.variables.tau
delta = self.variables.delta
u = self.variables.u
v = self.variables.v
eta = self.variables.eta
epsilon = self.variables.epsilon
zeta1 = self.variables.zeta1
zeta2 = self.variables.zeta2
q = self.variables.q
self.constraints.c8 = m.addConstrs(t[j] - tau[j] <= 0 for j in range(N))
self.constraints.c81 = m.addConstrs(tau[j] - t[j] - delta[j] <= 0 for j in range(N))
self.constraints.c9 = m.addConstrs(
tau[j] >= tau[i] + T[i, j] - M * (1 - x[k, i, j]) for j in range(N) for i in range(N) for k in range(K))
self.constraints.c10 = m.addConstrs(
eta[j] >= epsilon[j] + u[j] * deltabar - M * (1 - x.sum('*', '*', j)) for j in range(N))
# KKT condition of lower level problem
self.constraints.c11 = m.addConstrs(
-gamma * zeta1[j] + gamma * zeta2[j] - q[j] + u[j] - v[j] == 0 for j in range(N) if 0 < j < N - 1)
self.constraints.c12 = m.addConstrs(1 - zeta1[j] - zeta2[j] == 0 for j in range(N) if 0 < j < N - 1)
self.constraints.c13a = m.addConstrs( 0 <= epsilon[j] + gamma * delta[j] - gamma*deltabar/2
for j in range(N) if 0 < j < N - 1)
self.constraints.c13b = m.addConstrs(epsilon[j] + gamma * delta[j] - gamma*deltabar/2 <= M * omega1[j]
for j in range(N) if 0 < j < N - 1)
self.constraints.c14 = m.addConstrs(zeta1[j] <= M * (1 - omega1[j]) for j in range(N) if 0 < j < N - 1)
self.constraints.c15 = m.addConstrs( 0<= epsilon[j] - gamma * delta[j] + gamma*deltabar/2 for j in range(N) if 0 < j < N - 1)
self.constraints.c15 = m.addConstrs(epsilon[j] - gamma * delta[j] + gamma*deltabar/2 <= M * omega2[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c16 = m.addConstrs(zeta2[j] <= M * (1 - omega2[j]) for j in range(N) if 0 < j < N - 1)
self.constraints.c17 = m.addConstrs(deltabar - delta[j] <= M * omega3[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c18 = m.addConstrs(u[j] <= M * (1 - omega3[j]) for j in range(N) if 0 < j < N - 1)
self.constraints.c19 = m.addConstrs(delta[j] <= M * omega4[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c20 = m.addConstrs(v[j] <= M * (1 - omega4[j]) for j in range(N) if 0 < j < N - 1)
self.constraints.c8a = m.addConstrs( delta[j] >= 0 for j in range(N))
self.constraints.c8b = m.addConstrs( delta[j] <= deltabar for j in range(N))
self.constraints.cuts = {}
pass
###
# Cut adding
###
def _add_cut(self):
if self.submodel.model.Status == 2:
K = self.data.K
N = self.data.N
x = self.variables.x
omega1 = self.variables.omega1
omega2 = self.variables.omega2
omega3 = self.variables.omega3
omega4 = self.variables.omega4
cut = len(self.data.cutlist)
self.data.cutlist.append(cut)
# Get sensitivity from subproblem
sens_x = 0
sens_omega1 = 0
sens_omega2 = 0
sens_omega3 = 0
sens_omega4 = 0
for k in range(K):
for i in range(N):
for j in range(N):
sens_x += (x[k, i, j] - x[k, i, j].x) * self.submodel.constraints.fix_x[k, i, j].pi
sens_omega1 += (omega1[j] - omega1[j].x) * self.submodel.constraints.fix_omega1[j].pi
sens_omega2 += (omega2[j] - omega2[j].x) * self.submodel.constraints.fix_omega2[j].pi
sens_omega3 += (omega3[j] - omega3[j].x) * self.submodel.constraints.fix_omega3[j].pi
sens_omega4 += (omega4[j] - omega4[j].x) * self.submodel.constraints.fix_omega4[j].pi
z_sub = self.submodel.model.ObjVal
# Generate cut
self.constraints.cuts[cut] = self.model.addConstr(self.variables.theta >= z_sub + sens_x
+ sens_omega1 + sens_omega2+ sens_omega3 + sens_omega4)
elif self.submodel.model.Status == 3:
data = self.data
print('\n' + '#' * 50)
print('Feasibility Subproblem optimization')
self.feasibility_submodel = Benders_Feasibility_Subproblem(self,data.run,data.node , data.vehicle,data.delta_val,data.gamma) # Build feasibility subproblem from solution
self.feasibility_submodel.optimize()
K = self.data.K
N = self.data.N
x = self.variables.x
omega1 = self.variables.omega1
omega2 = self.variables.omega2
omega3 = self.variables.omega3
omega4 = self.variables.omega4
cut = len(self.data.cutlist)
self.data.cutlist.append(cut)
# Get sensitivity from subproblem
sens_x = 0
sens_omega1 = 0
sens_omega2 = 0
sens_omega3 = 0
sens_omega4 = 0
for k in range(K):
for i in range(N):
for j in range(N):
sens_x += (x[k, i, j] - x[k, i, j].x) * self.feasibility_submodel.constraints.fix_x[k, i, j].pi
sens_omega1 += (omega1[j] - omega1[j].x) * self.feasibility_submodel.constraints.fix_omega1[j].pi
sens_omega2 += (omega2[j] - omega2[j].x) * self.feasibility_submodel.constraints.fix_omega2[j].pi
sens_omega3 += (omega3[j] - omega3[j].x) * self.feasibility_submodel.constraints.fix_omega3[j].pi
sens_omega4 += (omega4[j] - omega4[j].x) * self.feasibility_submodel.constraints.fix_omega4[j].pi
z_sub = self.feasibility_submodel.model.ObjVal
self.data.z_fea_sub.append(z_sub)
# Generate feasibility cut
self.constraints.cuts[cut] = self.model.addConstr( 0 >= z_sub + sens_x + sens_omega1+ sens_omega2
+ sens_omega3 + sens_omega4)
###
# Update upper and lower bounds
###
def _update_bounds(self):
if self.submodel.model.Status == 2:
z_sub = self.submodel.model.ObjVal
z_master = self.model.ObjVal
self.data.ub = z_master - self.variables.theta.x + z_sub
# The best lower bound is the current bestbound,
# This will equal z_master at optimality
self.data.lb = self.model.ObjBound
self.data.upper_bounds.append(self.data.ub)
self.data.lower_bounds.append(self.data.lb)
def _save_vars(self):
K = self.data.K
N = self.data.N
self.data.thetas.append(self.variables.theta.x)
for k in range(K):
for i in range(N):
for j in range(N):
self.data.xs.append(self.variables.x[k,i,j].x)
# self.data.ys.append(self.submodel.variables.y.x)
# Subproblem
class Benders_Subproblem:
def __init__(self, MP, run, node, vehicle, delta_val,gamma):
self.data = expando()
self.variables = expando()
self.constraints = expando()
self.results = expando()
self._load_data(run, node, vehicle, delta_val,gamma)
self._build_model()
self.data.MP = MP
self.update_fixed_vars()
def optimize(self):
# self.model.Params.InfUnbdInfo = 1
self.model.optimize()
qq = 0
def _load_data(self, run, node, vehicle, delta_val,gamma):
self.data.cutlist = []
self.data.upper_bounds = []
self.data.lower_bounds = []
self.data.lambdas = {}
self.data.gamma = gamma
self.data.ub = gb.GRB.INFINITY
self.data.lb = -gb.GRB.INFINITY
self.data.xs = []
self.data.ys = []
self.data.thetas = []
data = loadmat('DataMap/Vehicles_{}/RealMap_{}_test.mat'.format(run, node))
self.data.Emax = data['Emax'][0][0]
self.data.rmax = self.data.Emax
self.data.e = data['e']
self.data.T = 1e0 * data['T']
self.data.M = int(data['M'])
self.data.p = data['p'][0]
self.data.g = data['g'][0]
self.data.t = 1e0 * data['t'][0]
self.data.c = data['c'][0]
self.data.K = int(vehicle)
self.data.N = int(data['N'])
# self.data.alpha = 0.05
# self.data.Q = 200
self.data.deltabar = delta_val
###
# Model Building
###
def _build_model(self):
self.model = gb.Model()
self._build_variables()
self._build_objective()
self._build_constraints()
self.model.update()
def _build_variables(self):
m = self.model
K = self.data.K
N = self.data.N
deltabar = self.data.deltabar
self.variables.theta = m.addVar(lb=-gb.GRB.INFINITY, ub=gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name='theta')
self.variables.x = m.addVars(K, N, N,lb = 0, ub = 1, vtype=gb.GRB.CONTINUOUS, name="x")
self.variables.omega1 = m.addVars(N,lb = 0, ub = 1, vtype=gb.GRB.CONTINUOUS, name="omega1")
self.variables.omega2 = m.addVars(N,lb = 0, ub = 1, vtype=gb.GRB.CONTINUOUS, name="omega2")
self.variables.omega3 = m.addVars(N,lb = 0, ub = 1, vtype=gb.GRB.CONTINUOUS, name="omega3")
self.variables.omega4 = m.addVars(N,lb = 0, ub = 1, vtype=gb.GRB.CONTINUOUS, name="omega4")
# self.variables.x = m.addVars(K, N, N,lb = 0, ub = 1, vtype=gb.GRB.BINARY, name="x")
# self.variables.omega1 = m.addVars(N,lb = 0, ub = 1, vtype=gb.GRB.BINARY, name="omega1")
# self.variables.omega2 = m.addVars(N,lb = 0, ub = 1, vtype=gb.GRB.BINARY, name="omega2")
# self.variables.omega3 = m.addVars(N,lb = 0, ub = 1, vtype=gb.GRB.BINARY, name="omega3")
# self.variables.omega4 = m.addVars(N,lb = 0, ub = 1, vtype=gb.GRB.BINARY, name="omega4")
self.variables.eta = m.addVars( N, vtype=gb.GRB.CONTINUOUS, name="eta")
self.variables.epsilon = m.addVars( N, vtype=gb.GRB.CONTINUOUS, name="epsilon")
self.variables.q = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="q")
# self.variables.r = m.addVars(N, K, vtype=gb.GRB.CONTINUOUS, name="r")
# self.variables.E = m.addVars(N, K, vtype=gb.GRB.CONTINUOUS, name="E")
self.variables.u = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="u")
self.variables.v = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="v")
self.variables.zeta1 = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="zeta1")
self.variables.zeta2 = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="zeta2")
self.variables.delta = m.addVars(N,lb = 0, ub = deltabar, vtype=gb.GRB.CONTINUOUS, name="delta")
self.variables.tau = m.addVars(N,vtype=gb.GRB.CONTINUOUS ,name="tau")
m.update()
def _build_objective(self):
m = self.model
N = self.data.N
qq=0
for i in range(N):
if i<N-1:
if i>0:
qq += self.variables.eta[i]
m.setObjective(qq, gb.GRB.MINIMIZE)
def _build_constraints(self):
m = self.model
K = self.data.K
M = self.data.M
N = self.data.N
t = self.data.t
T = self.data.T
deltabar = self.data.deltabar
gamma = self.data.gamma
x = self.variables.x
omega1 = self.variables.omega1
omega2 = self.variables.omega2
omega3 = self.variables.omega3
omega4 = self.variables.omega4
tau = self.variables.tau
delta = self.variables.delta
u = self.variables.u
v = self.variables.v
eta = self.variables.eta
epsilon = self.variables.epsilon
zeta1 = self.variables.zeta1
zeta2 = self.variables.zeta2
q = self.variables.q
self.constraints.c1 = m.addConstrs(x.sum('*', i, '*') <= 1 for i in range(N) if i != 0)
self.constraints.c2 = m.addConstrs(x[k, i, i] == 0 for i in range(N) for k in range(K))
self.constraints.c3 = m.addConstrs(x[k, N - 1, i] == 0 for i in range(N) for k in range(K))
self.constraints.c4 = m.addConstrs(x[k, i, 0] == 0 for i in range(N) for k in range(K))
# Flow conservation
self.constraints.c5 = m.addConstrs(x.sum(k, 0, '*') - x.sum(k, '*', 0) == 1 for k in range(K))
self.constraints.c6 = m.addConstrs(x.sum(k, N - 1, '*') - x.sum(k, '*', N - 1) == -1 for k in range(K))
self.constraints.c7 = m.addConstrs(x.sum(k, i, '*') - x.sum(k, '*', i) == 0 for i in range(N - 1) if i != 0 for k
in range(K))
self.constraints.c8a = m.addConstrs( t[j]- tau[j] <=0 for j in range(N) )
self.constraints.c8b = m.addConstrs( tau[j] - t[j]-delta[j] <= 0 for j in range(N) )
self.constraints.c9 = m.addConstrs( tau[j] >= tau[i] + T[i,j] - M*(1-x[k,i,j]) for j in range(N) for i in range(N) for k in range(K) )
self.constraints.c10 = m.addConstrs( eta[j] >= epsilon[j] + u[j]*deltabar - M*(1 - x.sum('*', '*',j) ) for j in range(N) )
# KKT condition of lower level problem
self.constraints.c11 = m.addConstrs(- gamma * zeta1[j] + gamma * zeta2[j] - q[j] + u[j] - v[j] == 0 for j in range(N) if 0 < j < N - 1)
self.constraints.c12 = m.addConstrs( 1- zeta1[j]- zeta2[j] == 0 for j in range(N) if 0 < j < N - 1)
self.constraints.c13a = m.addConstrs( 0<= epsilon[j] + gamma * delta[j] - gamma*deltabar/2 for j in range(N) if 0 < j < N - 1)
self.constraints.c13b = m.addConstrs( epsilon[j] + gamma * delta[j] - gamma*deltabar/2 <= M * omega1[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c14 = m.addConstrs(zeta1[j] <= M * (1 - omega1[j]) for j in range(N) if 0 < j < N - 1)
self.constraints.c15a = m.addConstrs( 0 <= epsilon[j] - gamma * delta[j] + gamma*deltabar/2 for j in range(N) if 0 < j < N - 1)
self.constraints.c15b = m.addConstrs( epsilon[j] - gamma * delta[j] + gamma*deltabar/2 <= M * omega2[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c16 = m.addConstrs(zeta2[j] <= M * (1 - omega2[j]) for j in range(N) if 0 < j < N - 1)
self.constraints.c17 = m.addConstrs(deltabar - delta[j] <= M * omega3[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c18 = m.addConstrs(u[j] <= M * (1 - omega3[j]) for j in range(N) if 0 < j < N - 1)
self.constraints.c19 = m.addConstrs(delta[j] <= M * omega4[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c20 = m.addConstrs(v[j] <= M * (1 - omega4[j]) for j in range(N) if 0 < j < N - 1)
self.constraints.fix_x = m.addConstrs(x[k, i, j] == 0 for k in range(K) for i in range(N) for j in range(N))
self.constraints.fix_omega1 = m.addConstrs(omega1[i] == 0 for i in range(N))
self.constraints.fix_omega2 = m.addConstrs(omega2[i] == 0 for i in range(N))
self.constraints.fix_omega3 = m.addConstrs(omega3[i] == 0 for i in range(N))
self.constraints.fix_omega4 = m.addConstrs(omega4[i] == 0 for i in range(N))
def update_fixed_vars(self, MP=None):
K = self.data.K
N = self.data.N
if MP is None:
MP = self.data.MP
for k in range(K):
for i in range(N):
for j in range(N):
self.constraints.fix_omega1[j].rhs = MP.variables.omega1[ j].x
self.constraints.fix_omega2[j].rhs = MP.variables.omega2[ j].x
self.constraints.fix_omega3[j].rhs = MP.variables.omega3[ j].x
self.constraints.fix_omega4[j].rhs = MP.variables.omega4[ j].x
self.constraints.fix_x[k, i, j].rhs = MP.variables.x[k, i, j].x
# Feasibility Subproblem
class Benders_Feasibility_Subproblem:
def __init__(self, MP, run, node, vehicle, delta_val,gamma):
self.data = expando()
self.variables = expando()
self.constraints = expando()
self.results = expando()
self._load_data(run, node, vehicle, delta_val,gamma)
self._build_model()
self.data.MP = MP
self.update_fixed_vars()
def optimize(self):
self.model.Params.InfUnbdInfo = 1
self.model.optimize()
qq = 0
def _load_data(self, run, node, vehicle, delta_val,gamma):
self.data.cutlist = []
self.data.upper_bounds = []
self.data.lower_bounds = []
self.data.lambdas = {}
self.data.gamma = gamma
self.data.ub = gb.GRB.INFINITY
self.data.lb = -gb.GRB.INFINITY
self.data.xs = []
self.data.ys = []
self.data.thetas = []
data = loadmat('DataMap/Vehicles_{}/RealMap_{}_test.mat'.format(run, node))
self.data.Emax = data['Emax'][0][0]
self.data.rmax = self.data.Emax
self.data.e = data['e']
self.data.T = 1e0 * data['T']
self.data.M = int(data['M'])
self.data.p = data['p'][0]
self.data.g = data['g'][0]
self.data.t = 1e0 * data['t'][0]
self.data.c = data['c'][0]
self.data.K = int(vehicle)
self.data.N = int(data['N'])
# self.data.alpha = 0.05
# self.data.Q = 200
self.data.deltabar = delta_val
###
# Model Building
###
def _build_model(self):
self.model = gb.Model()
self._build_variables()
self._build_objective()
self._build_constraints()
self.model.update()
def _build_variables(self):
m = self.model
K = self.data.K
N = self.data.N
deltabar = self.data.deltabar
self.variables.theta = m.addVar(lb=-gb.GRB.INFINITY, ub=gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name='theta')
self.variables.x = m.addVars(K, N, N,lb = 0, ub = 1, vtype=gb.GRB.CONTINUOUS, name="x")
self.variables.omega1 = m.addVars(N,lb = 0, ub = 1, vtype=gb.GRB.CONTINUOUS, name="omega1")
self.variables.omega2 = m.addVars(N,lb = 0, ub = 1, vtype=gb.GRB.CONTINUOUS, name="omega2")
self.variables.omega3 = m.addVars(N,lb = 0, ub = 1, vtype=gb.GRB.CONTINUOUS, name="omega3")
self.variables.omega4 = m.addVars(N,lb = 0, ub = 1, vtype=gb.GRB.CONTINUOUS, name="omega4")
self.variables.eta = m.addVars( N, vtype=gb.GRB.CONTINUOUS, name="eta")
self.variables.epsilon = m.addVars( N, vtype=gb.GRB.CONTINUOUS, name="epsilon")
self.variables.q = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="q")
# self.variables.r = m.addVars(N, K, vtype=gb.GRB.CONTINUOUS, name="r")
# self.variables.E = m.addVars(N, K, vtype=gb.GRB.CONTINUOUS, name="E")
self.variables.u = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="u")
self.variables.v = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="v")
self.variables.zeta1 = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="zeta1")
self.variables.zeta2 = m.addVars(N, vtype=gb.GRB.CONTINUOUS, name="zeta2")
self.variables.delta = m.addVars(N,lb = 0, ub = deltabar, vtype=gb.GRB.CONTINUOUS, name="delta")
self.variables.tau = m.addVars(N,vtype=gb.GRB.CONTINUOUS ,name="tau")
self.variables.Sthree = m.addVars(K, N, N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sthree") #Slack var.for feasibility
self.variables.Sone_1 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_1")
self.variables.Sone_2 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_2")
self.variables.Sone_3 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_3")
self.variables.Sone_4 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_4")
self.variables.Sone_5 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_5")
self.variables.Sone_6 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_6")
self.variables.Sone_7 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_7")
self.variables.Sone_8 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_8")
self.variables.Sone_9 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_9")
self.variables.Sone_10 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_10")
self.variables.Sone_11 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_11")
self.variables.Sone_12 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_12")
self.variables.Sone_13 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_13")
self.variables.Sone_14 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_14")
self.variables.Sone_15 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_15")
self.variables.Sone_m1 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_m1")
self.variables.Sone_m2 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_m2")
self.variables.Sone_m3 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_m3")
self.variables.Sone_m4 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_m4")
self.variables.Sone_m5 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_m5")
self.variables.Sone_m6 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_m6")
self.variables.Sone_m7 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_m7")
self.variables.Sone_m8 = m.addVars( N,lb = 0, ub = gb.GRB.INFINITY, vtype=gb.GRB.CONTINUOUS, name="Sone_m8")
m.update()
def _build_objective(self):
m = self.model
N = self.data.N
vars = self.variables
m.setObjective(vars.Sthree.sum('*', '*', '*')+ vars.Sone_1.sum('*')+ vars.Sone_2.sum('*')+ vars.Sone_3.sum('*')
+ vars.Sone_4.sum('*')+ vars.Sone_5.sum('*')+ vars.Sone_6.sum('*')+ vars.Sone_7.sum('*')
+ vars.Sone_8.sum('*')+ vars.Sone_9.sum('*')+ vars.Sone_10.sum('*')+ vars.Sone_11.sum('*')
+ vars.Sone_12.sum('*')+ vars.Sone_13.sum('*')+ vars.Sone_14.sum('*')+ vars.Sone_15.sum('*')
+ vars.Sone_m1.sum('*')+ vars.Sone_m2.sum('*')+ vars.Sone_m3.sum('*')+ vars.Sone_m4.sum('*')
+ vars.Sone_m5.sum('*')+ vars.Sone_m6.sum('*')+ vars.Sone_m7.sum('*')+ vars.Sone_m8.sum('*')
,
gb.GRB.MINIMIZE) #- vars.Sone_3[0] -vars.Sone_3[N-1]
def _build_constraints(self):
m = self.model
K = self.data.K
M = self.data.M
N = self.data.N
t = self.data.t
T = self.data.T
deltabar = self.data.deltabar
gamma = self.data.gamma
x = self.variables.x
omega1 = self.variables.omega1
omega2 = self.variables.omega2
omega3 = self.variables.omega3
omega4 = self.variables.omega4
tau = self.variables.tau
delta = self.variables.delta
u = self.variables.u
v = self.variables.v
eta = self.variables.eta
epsilon = self.variables.epsilon
zeta1 = self.variables.zeta1
zeta2 = self.variables.zeta2
q = self.variables.q
Sthree = self.variables.Sthree
Sone_1 = self.variables.Sone_1
Sone_2 = self.variables.Sone_2
Sone_3 = self.variables.Sone_3
Sone_4 = self.variables.Sone_4
Sone_5 = self.variables.Sone_5
Sone_6 = self.variables.Sone_6
Sone_7 = self.variables.Sone_7
Sone_8 = self.variables.Sone_8
Sone_9 = self.variables.Sone_9
Sone_10 = self.variables.Sone_10
Sone_11 = self.variables.Sone_11
Sone_12 = self.variables.Sone_12
Sone_13 = self.variables.Sone_13
Sone_14 = self.variables.Sone_14
Sone_15 = self.variables.Sone_15
Sone_m1 = self.variables.Sone_m1
Sone_m2 = self.variables.Sone_m2
Sone_m3 = self.variables.Sone_m3
Sone_m4 = self.variables.Sone_m4
Sone_m5 = self.variables.Sone_m5
Sone_m6 = self.variables.Sone_m6
Sone_m7 = self.variables.Sone_m7
Sone_m8 = self.variables.Sone_m8
self.constraints.c1 = m.addConstrs(x.sum('*', i, '*') <= 1 for i in range(N) if i != 0)
self.constraints.c2 = m.addConstrs(x[k, i, i] == 0 for i in range(N) for k in range(K))
self.constraints.c3 = m.addConstrs(x[k, N - 1, i] == 0 for i in range(N) for k in range(K))
self.constraints.c4 = m.addConstrs(x[k, i, 0] == 0 for i in range(N) for k in range(K))
# Flow conservation
self.constraints.c5 = m.addConstrs(x.sum(k, 0, '*') - x.sum(k, '*', 0) == 1 for k in range(K))
self.constraints.c6 = m.addConstrs(x.sum(k, N - 1, '*') - x.sum(k, '*', N - 1) == -1 for k in range(K))
self.constraints.c7 = m.addConstrs(x.sum(k, i, '*') - x.sum(k, '*', i) == 0 for i in range(N - 1) if i != 0 for k
in range(K))
self.constraints.c8a = m.addConstrs( t[j]- tau[j] - Sone_1[j] <=0 for j in range(N) )
self.constraints.c8b = m.addConstrs( tau[j] - t[j]-delta[j] - Sone_2[j] <= 0 for j in range(N) )
self.constraints.c9 = m.addConstrs( tau[j] + Sthree[k,i,j] >= tau[i] + T[i,j] - M*(1-x[k,i,j]) for j in range(N) for i in range(N) for k in range(K) )
self.constraints.c10 = m.addConstrs( eta[j] +Sone_3[j] >= epsilon[j] + u[j]*deltabar - M*(1 - x.sum('*', '*',j) ) for j in range(N) )
# KKT condition of lower level problem
self.constraints.c11a = m.addConstrs(- gamma * zeta1[j] + gamma * zeta2[j] - q[j] + u[j] - v[j]- Sone_4[j] <= 0 for j in range(N) if 0 < j < N - 1)
self.constraints.c11b = m.addConstrs(- gamma * zeta1[j] + gamma * zeta2[j] - q[j] + u[j] - v[j] +Sone_5[j] >= 0 for j in range(N) if 0 < j < N - 1)
self.constraints.c12a = m.addConstrs( 1- zeta1[j]- zeta2[j]- Sone_6[j] <= 0 for j in range(N) if 0 < j < N - 1)
self.constraints.c12b = m.addConstrs( 1- zeta1[j]- zeta2[j]+ Sone_7[j] >= 0 for j in range(N) if 0 < j < N - 1)
self.constraints.c13a = m.addConstrs( 0 <= epsilon[j] + gamma * delta[j] - gamma*deltabar/2 +Sone_m1[j] for j in range(N) if 0 < j < N - 1) # m=1, Slack var. for piecewise inconv. func.
self.constraints.c13b = m.addConstrs( epsilon[j] + gamma * delta[j]- gamma*deltabar/2 <= Sone_m2[j] + M * omega1[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c14a = m.addConstrs( 0<= zeta1[j] + Sone_m3[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c14b = m.addConstrs(zeta1[j] <= Sone_m4[j] + M * (1 - omega1[j]) for j in range(N) if 0 < j < N - 1)
self.constraints.c15a = m.addConstrs( 0 <= epsilon[j] - gamma * delta[j] + gamma*deltabar/2 + Sone_m5[j] for j in range(N) if 0 < j < N - 1) # m=2, Slack var. for piecewise inconv. func.
self.constraints.c15b = m.addConstrs( epsilon[j] - gamma * delta[j] + gamma*deltabar/2 <= Sone_m6[j]+ M * omega2[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c16a = m.addConstrs( 0<= zeta2[j] + Sone_m7[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c16b = m.addConstrs(zeta2[j] <= Sone_m8[j] + M * (1 - omega2[j]) for j in range(N) if 0 < j < N - 1)
self.constraints.c17a = m.addConstrs( 0<= deltabar - delta[j] + Sone_8[j] for j in range(N) if 0 < j < N - 1) # Linearizatio of comple. cons.
self.constraints.c17b = m.addConstrs(deltabar - delta[j] <= Sone_9[j] + M * omega3[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c18a = m.addConstrs( 0<= u[j] + Sone_10[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c18b = m.addConstrs(u[j] <= Sone_11[j] + M * (1 - omega3[j]) for j in range(N) if 0 < j < N - 1)
self.constraints.c19a = m.addConstrs( 0<= delta[j] + Sone_12[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c19b = m.addConstrs(delta[j] <= Sone_13[j] + M * omega4[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c20a = m.addConstrs( 0<= v[j] + Sone_14[j] for j in range(N) if 0 < j < N - 1)
self.constraints.c20b = m.addConstrs(v[j] <= Sone_15[j] + M * (1 - omega4[j]) for j in range(N) if 0 < j < N - 1)
self.constraints.fix_x = m.addConstrs(x[k, i, j] == 0 for k in range(K) for i in range(N) for j in range(N))
self.constraints.fix_omega1 = m.addConstrs(omega1[i] == 0 for i in range(N))
self.constraints.fix_omega2 = m.addConstrs(omega2[i] == 0 for i in range(N))
self.constraints.fix_omega3 = m.addConstrs(omega3[i] == 0 for i in range(N))
self.constraints.fix_omega4 = m.addConstrs(omega4[i] == 0 for i in range(N))
def update_fixed_vars(self, MP=None):
K = self.data.K
N = self.data.N
if MP is None:
MP = self.data.MP
for k in range(K):
for i in range(N):
for j in range(N):
self.constraints.fix_omega1[j].rhs = MP.variables.omega1[ j].x
self.constraints.fix_omega2[j].rhs = MP.variables.omega2[ j].x
self.constraints.fix_omega3[j].rhs = MP.variables.omega3[ j].x
self.constraints.fix_omega4[j].rhs = MP.variables.omega4[ j].x
self.constraints.fix_x[k, i, j].rhs = MP.variables.x[k, i, j].x
### To Run:
# m = Benders_Master(benders_gap=0.001, run=1, node=11, vehicle=3, delta_val=0.5)
# m.optimize()
##
if __name__ == '__main__':
gamma0 = np.array([0.05, 0.01])
deltabar0 = np.array([2, 1.5])
# tlimit = 3600 * 2
map = 100
vehicle = 5
delta_val = .5
index = 1
runmax = 1
benders_gap = 1e-4
gamma = 0.01
# m = Benders_Master(benders_gap, run + 1, node, vehicle, delta_val)
# for k in range(run):
for gam in range(1):
# if gam >0 :
gamma = gamma0[gam]
for del1 in range(1):
deltabar = deltabar0[del1]
for i in range(index):
for run in range(runmax):
node = 17 + 2 * i
ini_time = time.time()
m = Benders_Master(benders_gap, run + 1, node, vehicle, deltabar, gamma)
m.optimize()
end_time = time.time()
Runtime = end_time - ini_time
biterm = 0
for j in range(index):
if j < index - 1:
if j > 0:
biterm += m.submodel.variables.eta[j].x
with open(r'solution_Benders/BVRP_deltabar_{}_gamma_{}_gap_{}'
r'/Sol_Benders_BDD.csv'.format(delta_val, gamma, benders_gap), mode='a') as f:
Sol = writer(f, delimiter='\t', lineterminator='\n')
Sol.writerow(
['sol_map{}_N{}_K{}_run{}'.format(map, node, vehicle, run + 1), '%.2f' % (Runtime),
'%.8f' % (m.submodel.model.ObjVal), '%.2f ' % (m.data.lb), '%.2f ' % (m.data.ub),
'%i ' % (len(m.data.cutlist))])
f.close()0 -
Can you try moving the definition of \( \texttt{EPI} \) outside of the \( \texttt{Benders_Master} \) class? E.g.:
def EPI(model, where):
...
class Benders_Master:
...0 -
Hi Eli,
Thank you very much, it seems that I have fixed this issue.
Best,
Canqi
0
Post is closed for comments.
Comments
5 comments