Column Generation

Hello, 

I realize that each column generation formulation is specific to the problem at hand but I want to ensure that I'm on the right track.  A little bit of background, the problem itself is quite massive in it's full size so I wanted to explore something like multiple column insertion within each iteration.  You'll also notice that I delete the model within each iteration due to the sheer size of even the reduced model.  It actually runs much faster when I clear the memory and rebuild with modified/new column list each iteration.  Anyways, any feedback on the actual structure and any advice on further speeding up computations is very much appreciated.  I don't think there's a way to add attachments so I pasted below.  Thank you in advance - Taylor

import sys
import gc
from gurobipy import *
import pandas as pd
import numpy as np

T0 = np.load('T0.npy')
T1 = np.load('T1.npy')
T2 = np.load('T2.npy')
T3 = np.load('T3.npy')
T4 = np.load('T4.npy')
T6 = np.load('T6.npy')
T7 = np.load('T7.npy')
T8 = np.load('T8.npy')
T9 = np.load('T9.npy')
T10 = np.load('T10.npy')

EPS = .000001
#SAP = int(sys.argv[5])
SAP = 10
iter = 1
sets = 399
airports = 440
DA = airports*26
secmea = 26
secmea2 = 52
pclass = 2
riskt = 17
bsm = 4
csm = 9
asm = 12
psm = 26

E = np.array(T0)
K = np.array(T1[54:80,:440])
U = np.array(T2)
M = np.array(T10)
#U = np.array(T3.iloc[:26,0])
fr = np.array(T3[0,1:3])
SA = np.array(T4[:,:])
SL = np.array(T6)
qnew1 = np.transpose(np.array(T7))
qnew2 = np.transpose(np.array(T8))
Rx = np.transpose(np.array(T9))

#Initial Solution to run problem
T12 = pd.read_csv('uniquecombo.csv',header=None)
uniquecombo = np.array(T12)
uc = uniquecombo.shape
uc = int(uc[0])
numair = 440

p = Model("Pricing Problem")
p.Params.OutputFlag=0
p.read("p_param.prm")

dual = Model("Dual Problem")
dual.Params.OutputFlag=0
dual.read("dual_param.prm")

lambda0 = {}
lambda1 = {}
lambda2a = {}
lambda2b = {}
lambda3 = {}
lambda4 = {}
lambda5 = {}

lambda0 = dual.addVars(secmea,lb=0, vtype=GRB.INTEGER,name="lambda0")
lambda1 = dual.addVars(secmea,lb=0, vtype=GRB.INTEGER,name="lambda1")
lambda2a = dual.addVars(airports,lb=0, vtype=GRB.INTEGER, name="lambda2a")
lambda2b = dual.addVars(airports,lb=0, vtype=GRB.INTEGER, name="lambda2b")
lambda3 = dual.addVars(secmea,airports,lb=0,vtype=GRB.INTEGER,name='lambda3')
lambda4 = dual.addVars(secmea,airports,lb=0,vtype=GRB.INTEGER,name='lambda4')
lambda5 = dual.addVars(secmea,airports,lb=0,vtype=GRB.INTEGER,name='lambda5')

dual.update()
dobj1 = quicksum(lambda1[d]*U[d] - lambda0[d]*E[d] for d in range(0,secmea))
dobj2 = quicksum(lambda2a[k]-lambda2b[k] for k in range(0,airports))
dobj3 = quicksum(K[d,k]*lambda3[d,k] + K[d,k]*lambda4[d,k] for k in range(0,airports) for d in range(0,secmea))
dual.setObjectiveN(dobj1,index=0,priority=0,weight=1)
dual.setObjectiveN(dobj2,index=1,priority=0,weight=1)
dual.setObjectiveN(dobj3,index=2,priority=0,weight=1)

d0 = {}
d1 = {}
d2 = {}

for j in range(0,uc):
#For the r variable
d0[j] = dual.addConstr(quicksum(lambda2a[k]-lambda2b[k] + quicksum(K[d,k]*lambda4[d,k] for d in range(0,secmea)) for k in range(0,airports))
>= (Rx[int(uniquecombo[j])] - ((1-SA[SAP,1])*fr[0]*qnew1[int(uniquecombo[j])]) - ((1-SA[SAP,2])*fr[1]*qnew2[int(uniquecombo[j])])),"Dual_Constraint1_"+str(j))
for k in range(0,airports):
#For the z variable
d1[k,j] = dual.addConstr(quicksum(lambda1[d] - lambda0[d] + lambda3[d,k] - lambda4[d,k] + lambda5[d,k] for d in range(0,secmea))
>= SL[k,int(uniquecombo[j])],"Dual_Constraint3_"+str(k)+str(j))
#For the s variable
d2[0] = dual.addConstr(quicksum(lambda4[d,k] - lambda5[d,k] for d in range(0,secmea) for k in range(0,airports)) >= 0,"Dual_Constraint2")

dual.update()
dual.optimize()
print("Dual Complete")

dvar = np.asarray(dual.X)
sec2 = 2*secmea
sec2a = sec2+2*airports
smar = secmea*airports

lvar0 = dvar[0:secmea] #single dimension
lvar1 = dvar[secmea:sec2] #single dimension
lvar2a = dvar[sec2:(sec2+airports)]
lvar2b = dvar[(sec2+airports):(sec2+2*airports)]
lvar3 = dvar[sec2a:(sec2a+smar)] #security measure by airport, d,k
lvar4 = dvar[(sec2a+smar):(sec2a+2*smar)] #security measure by airport, d,k
lvar5 = dvar[(sec2a+2*smar):] #security measure by airport, d,k

lvar3 = np.reshape(lvar3,(airports,secmea)).transpose()
lvar4 = np.reshape(lvar4,(airports,secmea)).transpose()
lvar5 = np.reshape(lvar5,(airports,secmea)).transpose()

#Now set up initial pricing problem and add variables
p.update()
# This model help us find the right variables to consider from P
# Pricing problem has the same decision variables as the master problem
r = {}
s = {}
z = {}
r = p.addVars(airports,uc,lb=0,ub=1,vtype=GRB.BINARY,name='r')
s = p.addVars(secmea,airports,uc,lb=0,vtype=GRB.INTEGER,name='s')
z = p.addVars(secmea,airports,uc,lb=0,vtype=GRB.INTEGER,name='z')
p0 = {}
p1 = {}
p2 = {}
p3 = {}
p4 = {}
p5 = {}
for d in range(0,secmea):
p0[d] = p.addConstr(quicksum(z[d,k,j] for k in range(0,airports) for j in range(0,uc)) >= E[d],"Existing_Avail_"+str(d))
p1[d] = p.addConstr(quicksum(z[d,k,j] for k in range(0,airports) for j in range(0,uc)) <= U[d],"Resource_Avail_"+str(d))
for k in range(0,airports):
p2[k] = p.addConstr(quicksum(r[k,j] for j in range(0,uc)) == 1, 'Airport_'+str(k))
for d in range(0,secmea):
p3[d,k] = p.addConstr(quicksum(z[d,k,j] for j in range(0,uc)) <= quicksum(K[d,k] for j in range(0,uc)),"linearity1_"+str(d)+str(k))
p4[d,k] = p.addConstr(quicksum(s[d,k,j] - (1-r[k,j])*K[d,k] for j in range(0,uc)) <= quicksum(z[d,k,j] for j in range(0,uc)),"linearity2_"+str(d)+str(k))
p5[d,k] = p.addConstr(quicksum(z[d,k,j] for j in range(0,uc)) <= quicksum(s[d,k,j] for j in range(0,uc)),"linearity3_"+str(d)+str(k))
p.update()

Cost1 = np.zeros((uc))
RC1 = np.zeros((airports))
RC2 = np.zeros((secmea,airports))
RC3 = np.zeros((airports))

for j in range(0,uc):
Cost1[j] = (Rx[int(uniquecombo[j])] - ((1-SA[SAP,1])*fr[0]*qnew1[int(uniquecombo[j])]) -
((1-SA[SAP,2])*fr[1]*qnew2[int(uniquecombo[j])]))
for k in range(0,airports):
RC1[k] = sum(lvar2a[k]-lvar2b[k] + K[d,k]*lvar4[d,k] for d in range(0,secmea))
RC3[k] = sum(lvar1[d] - lvar0[d] + lvar3[d,k] - lvar4[d,k] + lvar5[d,k] for d in range(0,secmea))
for d in range(0,secmea):
RC2[d,k] = (lvar4[d,k] - lvar5[d,k]) #reduced cost for the s variable

objn1p = quicksum((Cost1[j] - r[k,j]*RC1[k]) for k in range(0,airports) for j in range(0,uc))
objn2p = quicksum((0 - s[d,k,j]*RC2[d,k]) for k in range(0,airports) for j in range(0,uc) for d in range(0,secmea))
objn3p = quicksum((SL[k,int(uniquecombo[j])] - z[d,k,j]*RC3[k]) for k in range(0,airports) for j in range(0,uc) for d in range(0,secmea))

p.setObjectiveN(objn1p,index=0,priority=0,weight=1)
p.setObjectiveN(objn2p,index=1,priority=0,weight=1)
p.setObjectiveN(objn3p,index=2,priority=0,weight=1)

p.update()
p.optimize()

for num in range (0,p.NumObj):
p.setParam(GRB.Param.ObjNumber,num)
print('Obj:',p.objNVal)
print("Pricing Complete")

#print('\n')

dvar = np.asarray(p.X)
xc = uc*numair
rvar = dvar[0:xc]
svar = dvar[xc:(xc+xc*secmea)]
zvar = dvar[(xc+xc*secmea):]
rvar = np.array(rvar)
svar = np.array(svar)
zvar = np.array(zvar)
RSOL = np.reshape(rvar,(uc,numair)).transpose()
SSOL = np.reshape(svar,(uc,numair,secmea)).transpose()
ZSOL = np.reshape(zvar,(uc,numair,secmea)).transpose()
#np.savetxt("RSOL.csv",RSOL,delimiter=",")

uniquecombo = np.reshape(uniquecombo,(uc)).transpose()
CSOL = RSOL*uniquecombo
CSOL = np.delete(CSOL,np.where(~CSOL.any(axis=0))[0], axis=1)
uniquecombo = np.unique(CSOL)

qty = 300
uc1 = uc
uc2 = uniquecombo.shape
uc2 = int(uc2[0])
if uc2 == uc1:
comborand = range(0,211909)
np.random.shuffle(comborand)
crsub = comborand[:qty]
crsub = np.array(crsub)
uniquecombo = np.concatenate((uniquecombo,crsub),axis=0)
uniquecombo = np.unique(uniquecombo)
uc2 = uniquecombo.shape
uc2 = int(uc2[0])

Iter = 0
if Iter % 20 == 0:
print('Iteration DualValue PricingValue Master Problem')
print('%8d %12.5g %12.5g' % (Iter, dual.ObjVal, p.ObjVal))

del dual
del p
gc.collect()

removalbig = 0
Iter += 1

while Iter < 410:
dual = Model("Dual Problem")
dual.Params.OutputFlag=0
dual.read("dual_param.prm")
dual.update()

lambda0 = {}
lambda1 = {}
lambda2a = {}
lambda2b = {}
lambda3 = {}
lambda4 = {}
lambda5 = {}

lambda0 = dual.addVars(secmea,lb=0, vtype=GRB.INTEGER,name="lambda0")
lambda1 = dual.addVars(secmea,lb=0, vtype=GRB.INTEGER,name="lambda1")
lambda2a = dual.addVars(airports,lb=0, vtype=GRB.INTEGER, name="lambda2a")
lambda2b = dual.addVars(airports,lb=0, vtype=GRB.INTEGER, name="lambda2b")
lambda3 = dual.addVars(secmea,airports,lb=0,vtype=GRB.INTEGER,name='lambda3')
lambda4 = dual.addVars(secmea,airports,lb=0,vtype=GRB.INTEGER,name='lambda4')
lambda5 = dual.addVars(secmea,airports,lb=0,vtype=GRB.INTEGER,name='lambda5')

dual.update()
dobj1 = quicksum(lambda1[d]*U[d] - lambda0[d]*E[d] for d in range(0,secmea))
dobj2 = quicksum(lambda2a[k]-lambda2b[k] for k in range(0,airports))
dobj3 = quicksum(K[d,k]*lambda3[d,k] + K[d,k]*lambda4[d,k] for k in range(0,airports) for d in range(0,secmea))
dual.setObjectiveN(dobj1,index=0,priority=0,weight=1)
dual.setObjectiveN(dobj2,index=1,priority=0,weight=1)
dual.setObjectiveN(dobj3,index=2,priority=0,weight=1)

## delete only if file exists ##
d0 = {}
d1 = {}
d2 = {}

for j in range(0,uc2):
#For the r variable
d0[j] = dual.addConstr(quicksum(lambda2a[k]-lambda2b[k] + quicksum(K[d,k]*lambda4[d,k] for d in range(0,secmea)) for k in range(0,airports))
>= (Rx[int(uniquecombo[j])] - ((1-SA[SAP,1])*fr[0]*qnew1[int(uniquecombo[j])]) - ((1-SA[SAP,2])*fr[1]*qnew2[int(uniquecombo[j])])),"Dual_Constraint1_"+str(j))
for k in range(0,airports):
#For the z variable
d1[k,j] = dual.addConstr(quicksum(lambda1[d] - lambda0[d] + lambda3[d,k] - lambda4[d,k] + lambda5[d,k] for d in range(0,secmea))
>= SL[k,int(uniquecombo[j])],"Dual_Constraint3_"+str(k)+str(j))
#For the s variable
d2[0] = dual.addConstr(quicksum(lambda4[d,k] - lambda5[d,k] for d in range(0,secmea) for k in range(0,airports)) >= 0,"Dual_Constraint2")

dual.update()
dual.optimize()
duval = dual.objVal
dvar = np.asarray(dual.X)

del dual
gc.collect()
#print('\n')
#solve the dual
#use dual solution in the pricing problem
#solve pricing problem

#Set up the pricing problem based off of the dual
sec2 = 2*secmea
sec2a = sec2+2*airports
smar = secmea*airports

lvar0 = dvar[0:secmea] #single dimension
lvar1 = dvar[secmea:sec2] #single dimension
lvar2a = dvar[sec2:(sec2+airports)]
lvar2b = dvar[(sec2+airports):(sec2+2*airports)]
lvar3 = dvar[sec2a:(sec2a+smar)] #security measure by airport, d,k
lvar4 = dvar[(sec2a+smar):(sec2a+2*smar)] #security measure by airport, d,k
lvar5 = dvar[(sec2a+2*smar):] #security measure by airport, d,k

lvar3 = np.reshape(lvar3,(airports,secmea)).transpose()
lvar4 = np.reshape(lvar4,(airports,secmea)).transpose()
lvar5 = np.reshape(lvar5,(airports,secmea)).transpose()

#Now add variables to pricing problem, should be able to just add columns to system
p = Model("Pricing Problem")
p.Params.OutputFlag=0
p.read("p_param.prm")
p.update()
p.reset(0)

r = {}
s = {}
z = {}
r = p.addVars(airports,uc2,lb=0,ub=1,vtype=GRB.BINARY,name='r')
s = p.addVars(secmea,airports,uc2,lb=0,vtype=GRB.INTEGER,name='s')
z = p.addVars(secmea,airports,uc2,lb=0,vtype=GRB.INTEGER,name='z')
p0 = {}
p1 = {}
p2 = {}
p3 = {}
p4 = {}
p5 = {}
for d in range(0,secmea):
p0[d] = p.addConstr(quicksum(z[d,k,j] for k in range(0,airports) for j in range(0,uc2)) >= E[d],"Existing_Avail_"+str(d))
p1[d] = p.addConstr(quicksum(z[d,k,j] for k in range(0,airports) for j in range(0,uc2)) <= U[d],"Resource_Avail_"+str(d))
for k in range(0,airports):
p2[k] = p.addConstr(quicksum(r[k,j] for j in range(0,uc2)) == 1, 'Airport_'+str(k))
for d in range(0,secmea):
p3[d,k] = p.addConstr(quicksum(z[d,k,j] for j in range(0,uc2)) <= quicksum(K[d,k] for j in range(0,uc2)),"linearity1_"+str(d)+str(k))
p4[d,k] = p.addConstr(quicksum(s[d,k,j] - (1-r[k,j])*K[d,k] for j in range(0,uc2)) <= quicksum(z[d,k,j] for j in range(0,uc2)),"linearity2_"+str(d)+str(k))
p5[d,k] = p.addConstr(quicksum(z[d,k,j] for j in range(0,uc2)) <= quicksum(s[d,k,j] for j in range(0,uc2)),"linearity3_"+str(d)+str(k))
p.update()

Cost1 = np.zeros((uc2))
RC1 = np.zeros((airports))
RC2 = np.zeros((secmea,airports))
RC3 = np.zeros((airports))

for j in range(0,uc2):
Cost1[j] = (Rx[int(uniquecombo[j])] - ((1-SA[SAP,1])*fr[0]*qnew1[int(uniquecombo[j])]) -
((1-SA[SAP,2])*fr[1]*qnew2[int(uniquecombo[j])]))
for k in range(0,airports):
RC1[k] = sum(lvar2a[k]-lvar2b[k] + K[d,k]*lvar4[d,k] for d in range(0,secmea))
RC3[k] = sum(lvar1[d] - lvar0[d] + lvar3[d,k] - lvar4[d,k] + lvar5[d,k] for d in range(0,secmea))
for d in range(0,secmea):
RC2[d,k] = (lvar4[d,k] - lvar5[d,k]) #reduced cost for the s variable

objn1p = quicksum((Cost1[j] - r[k,j]*RC1[k]) for k in range(0,airports) for j in range(0,uc2))
objn2p = quicksum((0 - s[d,k,j]*RC2[d,k]) for k in range(0,airports) for j in range(0,uc2) for d in range(0,secmea))
objn3p = quicksum((SL[k,int(uniquecombo[j])] - z[d,k,j]*RC3[k]) for k in range(0,airports) for j in range(0,uc2) for d in range(0,secmea))

p.setObjectiveN(objn1p,index=0,priority=0,weight=1)
p.setObjectiveN(objn2p,index=1,priority=0,weight=1)
p.setObjectiveN(objn3p,index=2,priority=0,weight=1)

p.update()
p.optimize()
#print('\n')

if p.Status != GRB.OPTIMAL:
raise('Unexpected optimization status')

# If improvement between master prblem and the pricing problem is too small, then stop iterations
# if p.ObjVal > -1.0001:
# break

pval = p.objVal
dvar = np.asarray(p.X)
del p
gc.collect()

xc = uc2*numair
rvar = dvar[0:xc]
svar = dvar[xc:(xc+xc*secmea)]
zvar = dvar[(xc+xc*secmea):]
rvar = np.array(rvar)
svar = np.array(svar)
zvar = np.array(zvar)
RSOL = np.reshape(rvar,(uc2,numair)).transpose()
SSOL = np.reshape(svar,(uc2,numair,secmea)).transpose()
ZSOL = np.reshape(zvar,(uc2,numair,secmea)).transpose()
#np.savetxt("RSOL.csv",RSOL,delimiter=",")

#uniquecombo = np.reshape(uniquecombo,(uc2)).transpose()
indr, indc = np.where(RSOL == 1)
uniqueopt = uniquecombo[indc]
uniqueopt = np.unique(uniqueopt)
# if uniqueopt.shape != indc.shape:
# removal, x_ind, y_ind = np.setdiff1d(uniquecombo,uniqueopt,return_indices=True)
# indr, indc = np.where(uniquecombo != uniqueopt)
# removal = uniquecombo[[indc]]
# removalbig = np.concatenate((removalbig,removal),axis=0)

#This tells me all of the unique combos that should be kept for the optimization
#x_ind is the index of the combo from the uniquecombo array
uniquecombo, x_ind, y_ind = np.intersect1d(uniquecombo,uniqueopt,return_indices=True)

uc1 = uc2
crsub = comborand[(Iter*qty):((Iter+1)*qty)]
crsub = np.array(crsub)
uniquecombo = np.concatenate((uniquecombo,crsub),axis=0)
uniquecombo = np.unique(uniquecombo)
# indr, indc = np.where(uniquecombo != removalbig)
# uniquecombo = uniquecombo[[indc]]
uc2 = uniquecombo.shape
uc2 = int(uc2[0])

if Iter % 20 == 0:
print('Iteration DualValue PricingValue Master Problem')
print('%8d %12.5g %12.5g %12.5g' % (Iter, dual.ObjVal, p.ObjVal, m.ObjVal))
Iter += 1
#Any variable that isn't selected needs to be priced out and stored in removal array

#Once here, weve solved the relaxation of m (to epsilon-optimality), now we solve the reduced MIP

#Set up the Master
uc = uniquecombo.shape
uc = int(uc[0])

#Set up the Master
m = Model("Master Problem")
m.Params.OutputFlag=0
m.read("m_param.prm")

m.update()
r = {}
s = {}
z = {}
r = m.addVars(airports,uc,lb=0,ub=1,vtype=GRB.BINARY,name='r')
s = m.addVars(secmea,airports,uc,lb=0,vtype=GRB.INTEGER,name='s')
z = m.addVars(secmea,airports,uc,lb=0,vtype=GRB.INTEGER,name='z')
objn1 = quicksum(r[k,j]*Rx[int(uniquecombo[j])] for k in range(0,airports) for j in range(0,uc))
objn2 = -1*quicksum(
(1-SA[SAP,1])*fr[0]*r[k,j]*qnew1[int(uniquecombo[j])] + (1-SA[SAP,2])*fr[1]*r[k,j]*qnew2[int(uniquecombo[j])]
for k in range(0,airports) for j in range(0,uc))
objn3 = quicksum(SL[k,int(uniquecombo[j])]*z[d,k,j] for d in range(0,secmea) for k in range(0,airports) for j in range(0,uc))

m.setObjectiveN(objn1,index=0,priority=0,weight=1)
m.setObjectiveN(objn2,index=1,priority=0,weight=1)
m.setObjectiveN(objn3,index=2,priority=0,weight=1)

c0 = {}
c1 = {}
c2 = {}
c3 = {}
c4 = {}
c5 = {}
for d in range(0,secmea):
c0[d] = m.addConstr(quicksum(z[d,k,j] for k in range(0,airports) for j in range(0,uc)) >= E[d],"Existing_Avail_"+str(d))
c1[d] = m.addConstr(quicksum(z[d,k,j] for k in range(0,airports) for j in range(0,uc)) <= U[d],"Resource_Avail_"+str(d))
for k in range(0,airports):
c2[k] = m.addConstr(quicksum(r[k,j] for j in range(0,uc)) == 1, 'Airport_'+str(k))
for d in range(0,secmea):
c3[d,k] = m.addConstr(quicksum(z[d,k,j] for j in range(0,uc)) <= quicksum(K[d,k] for j in range(0,uc)),"linearity1_"+str(d)+str(k))
c4[d,k] = m.addConstr(quicksum(s[d,k,j] - (1-r[k,j])*K[d,k] for j in range(0,uc)) <= quicksum(z[d,k,j] for j in range(0,uc)),"linearity2_"+str(d)+str(k))
c5[d,k] = m.addConstr(quicksum(z[d,k,j] for j in range(0,uc)) <= quicksum(s[d,k,j] for j in range(0,uc)),"linearity3_"+str(d)+str(k))

m.update
m.optimize()
m.write('final.lp')
rootbound = m.ObjVal

print('Proven Column Generation Gap %.3f%%' % (100*(m.ObjVal-rootbound)/m.ObjVal))
print("Columns Selected, ",uc)
print(m.ObjVal)
for num in range (0,m.NumObj):
m.setParam(GRB.Param.ObjNumber,num)
print('Obj:',m.objNVal)
print("Final Master Complete")
np.savetxt('ColGenUniqueCombos.csv',uniquecombo,delimiter=",")

#Note that we can only ensure that there is no solution better than rootbound, because we solved (M) on a restricted set of variables and not in the full set of variables. However, these gaps tend to be fairly small, but is something to double check.