Dear community team,
I am using docplex to apply a branch and cut algorithm by taking advantage of LazyConstraintCallback package of cplex.callbacks.
The main problem is that after custom callback is called for 2 times, the code stop running eventhough the problem is not solved. What I mean is that, 2 cuts are added to the problem and the code stops running, but the problem is not solved yet, it should add more cuts to be able to reach optimality.
my code is as follows ;
# -*- coding: utf-8 -*-
"""
Created on Tue Mar 8 23:47:33 2022
@author: 90531
"""
import sys
from math import fabs
import docplex
from cplex.callbacks import LazyConstraintCallback
from cplex.exceptions import CplexError
from docplex.mp.model import Model
from docplex.mp.solution import *
from docplex.mp.callbacks.cb_mixin import *
import pandas as pd
import numpy as np
class nazli(ConstraintCallbackMixin, LazyConstraintCallback):
# Callback constructor. Model object is set after registration.
def __init__(self, env):
LazyConstraintCallback.__init__(self, env)
ConstraintCallbackMixin.__init__(self)
self.eps = 1e-6
self.sepcpx = Model()
def ali(self, ct):
self.register_constraint(ct)
def SepRoutine(self,w,d, s, t_indexx, timeuntil):
#Create single scenario optimization model model,
gama = 0.12
l = 1
sepcpx = self.sepcpx
#Add decision variables
Vs = sepcpx.continuous_var_matrix(timeuntil, w+1, lb=0, ub=1193637041 , name = "Vs")
#Add constraint
sepcpx.add_constraints(( Vs[t,w] >= d[t,w] for t in range(timeuntil) ))
#Add objective function as minimization
objExpr = sepcpx.linear_expr()
objExpr.add_term(Vs[t_indexx,w],1)
sepcpx.set_objective("min", objExpr)
#Solve scenario optimization problem for given alpha for each scenario
sol = sepcpx.solve()
soln = sol.get_objective_value()
print("soln for ",w, " :",soln)
sepcpx.clear()
return soln
def __call__(self):
print("defcall a girdi")
solv = self.make_solution_from_vars(self.Vs.values())
sol = self.make_solution()
solz = self.make_solution_from_vars(self.z)
z = self.z
Vs = self.Vs
d = self.d
s = self.s
gama = 0.12
l = 1
sepcpx = self.sepcpx
timeuntil = self.timeuntil
tolerance = self.tolerance
numberScenario = self.numberScenario
# Get the current Vs solution
VsT = solv.as_df()
VsT = VsT["value"].to_numpy()
VsSol = np.zeros((timeuntil, numberScenario))
for t in range(timeuntil):
for w in range(numberScenario):
VsSol[t][w] = VsT[t+w]
# Get the current z solution
zs = solz.as_df()
zs = zs["value"].to_numpy()
zSol = []
for w in range(numberScenario):
zSol.append(zs[w])
print(zSol)
p = int(tolerance * numberScenario)
# Cut Seperation
viol = False
z_violate = []
cikis = 0
for w in range(numberScenario):
for t in range(timeuntil):
if VsSol[t][w] < d[t][w] and zSol[w] < 1:
viol = True
z_violate.append(w)
cikis = 1
if cikis ==1:
cikis = 0
break
fark =[]
a = 0
for w in z_violate :
w = int(w)
tviol = []
t_index = []
a = 0
for t in range(timeuntil):
a = (d[t][w] - VsSol[t][w])
tviol.append(a)
t_index.append(t)
a = 0
t_h_viol = np.column_stack((t_index,tviol))
t_h_viol = pd.DataFrame(t_h_viol,columns=("t index","t viol"))
t_h_viol = t_h_viol.sort_values(by="t viol", ascending=False)
t_h_viol = t_h_viol.to_numpy(dtype=int)
t_indexx = t_h_viol[0][0]
if not viol:
print("violation yokkkk")
if viol :
print("viol:", viol)
h_alpha = []
h_alpha_indx = []
for omega in range(numberScenario):
solution = self.SepRoutine(omega,d,s, t_indexx, timeuntil)
h_alpha.append(solution)
h_alpha_indx.append(omega)
sigma = np.column_stack((h_alpha_indx,h_alpha))
sigma_x = pd.DataFrame(sigma,columns=("h index","solution"))
sigma_indx = sigma_x.sort_values(by="solution", ascending=False)
sigma_indx = sigma_indx.to_numpy()
h_index = sigma_indx[:,0]
h_value = sigma_indx[:,1]
sce_viol = VsSol[t_indexx][w]
for i in range(1,p+1): #since i = {2,...,p} for inequalities (11) on page 9
print("i= ", i)
sep = sce_viol + zSol[int(h_index[0])] * (h_value[0] - h_value[i]) + zSol[int(h_index[i])] * (h_value[i] - h_value[p])
cutRhs = h_value[0]
ncuts = 0
if(sep<cutRhs):
cutLhs = sepcpx.linear_expr()
cutLhs.add_term(z[int(h_index[0])] , float((h_value[0] - h_value[i])))
cutLhs.add_term(z[int(h_index[i])] , float((h_value[i] - h_value[p])))
print("h1 : ", float((h_value[i] - h_value[p])))
print("h2 : ", float((h_value[0] - h_value[i])))
ncuts = ncuts + 1
if ncuts > 0:
cutLhs.add_term(Vs[t_indexx,w],1.0)
cutRhs = cutRhs.item()
ct = (cutLhs >= cutRhs )
print("cut : ",ct)
self.ali(ct)
unsats = self.get_cpx_unsatisfied_cts(self.cts, sol, tolerance=1e-18)
for ct, cpx_lhs, sense, cpx_rhs in unsats:
print("rhs : ",cpx_rhs)
self.add(cpx_lhs, sense, cpx_rhs)
break
break
class Data:
def __init__(self, filename):
# read the data in filename
solar = pd.read_excel(filename,sheet_name="solar_baby", header=None)
self.solar = solar.to_numpy()
demand = pd.read_excel(filename,sheet_name="demand", header=None)
self.demand = demand.to_numpy()
try:
use_cuts = True
filename = 'C:\\Users\\90531\\Dropbox\\My PC (LAPTOP-6JH4R671)\\Desktop\\python codes\\data2.xlsx'
data = Data(filename)
tolerance = 0.3
d = data.demand
s = data.solar
#define sets
timeuntil =10
numberScenario = 7
#Define parameters
cS =200
l = 1
gama =0.12
BigM=250001
p = int(tolerance * numberScenario)
#Create MasterILP
cpx = Model("master problem")
Smax = cpx.continuous_var(name = "Smax", lb=0)
Vs = cpx.continuous_var_matrix(timeuntil, numberScenario, lb=0, name = "Vs")
Ks = cpx.continuous_var_matrix(timeuntil, numberScenario, lb=0, name = "Ks")
z = cpx.integer_var_list(numberScenario, lb=0, ub=1, name = "z")
const2 = cpx.add_constraints((Vs[t,w] + Ks[t,w] == Smax*gama*s[t,w]*l for t in range(timeuntil) for w in range(numberScenario)), names = "const11")
cpx.add_constraint(sum(z[w] for w in range(numberScenario)) <= p)
cpx.minimize(cS*Smax)
if use_cuts:
nazlibenders = cpx.register_callback(nazli)
nazlibenders.z = z
nazlibenders.Vs = Vs
nazlibenders.d = d
nazlibenders.s = s
nazlibenders.tolerance = tolerance
nazlibenders.timeuntil = timeuntil
nazlibenders.numberScenario = numberScenario
cpx.lazy_callback = nazlibenders
sol = cpx.solve()
except CplexError as exc:
print (exc)
TotCost_sa = sol.get_objective_value()
VS_sa= sol.get_value_df(Vs)
z_sa= sol.get_values(z)
Smax_sa = sol.get_value(Smax)
Regards,
Nazli Kalkan
------------------------------
Nazli Kalkan
------------------------------
#DecisionOptimization