-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathRun_Bilevel_KP.py
379 lines (343 loc) · 12.4 KB
/
Run_Bilevel_KP.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
"""
author: Margarida Carvalho
License: MIT
PYTHON 2
Run CCLW approach
Requires to install knapsack generator: http://hjemmesider.diku.dk/~pisinger/codes.html
In this code the knapsack generator is ./gen2
"""
from time import time
import os
from random import randint,seed
from gurobipy import *
# make xOpt from the MIP maximal (with respect to the MIP)
# ratio = Pi/wi by decreasing order
def MakeMaximal(ratio,F,c_F,xOpt,n,R,L):
i = 0
while c_F>0 and i<=n-1:
lixo, item = ratio[i]
if xOpt[item] <0.5: # xOpt=0
c_F = c_F - F[item]
i = i+1
while i<=n-1:
lixo, item = ratio[i]
if xOpt[item]==0 and R-L[item]>=0:
R = R - L[item]
xOpt[item] = 1
i = i+1
return xOpt
# ADD CUT: Follower does not improve constraint
# INPUT
# model - name of the model
# yOpt - we dont want solutions that let these items available
# n - number of items
# k - if k>2 then we can erase the last constraint added
# Profit - best (minor) value for the objective function found
# wmax - item with the largest weight
# OUTPUT
# xnew - return new solution of the leader
# m - new model
def NoImpFollower4(m,x,z,yOpt,n,Profit,wmax,k,P,upRow,oldProfit):
if k>2:
m.remove(m.getConstrs()[-1])
m.update()
# row generation cut
rowGen = LinExpr(0)
# add no improvement constraint for the follower
for i in range(1,n+1):
if yOpt[i-1]>0.5: # sum_ {i: yi=1} xi >=1
rowGen = rowGen + P[i-1]*(1-x[i])
m.addConstr(rowGen <= Profit-1)
m.update()
# it remains to update the rhs of rowGen constraints
if upRow == 1:
for j in range(k-2): # k-2 cutting plane constraints to update
m.getConstrs()[-j-2].setAttr("rhs",m.getConstrs()[-j-2].getAttr("rhs")+Profit-oldProfit) # remove the old profit to the new one
#m.getConstrs()[-j-2].setAttr("rhs",Profit-oldProfit) # remove the old profit to the new one
m.update()
m.update()
# add improve4 idea
m.addConstr( m.getObjective()-z[1]*wmax<= Profit-1)
m.update()
# # save problem - OPTIONAL
# m.write('current.lp')
# solve
m.optimize()
# show solution x and
xnew=[] # when it is empty we know that the problem is infeasible
value = 0
if m.Status == GRB.OPTIMAL:
value = m.ObjVal
for i in range(1,n+1):
xnew.append(x[i].x)
if m.status == 9:
xnew = 0 # time limit reached
return xnew,m,value
# Just stops when it has found a solution (i.e., the problem becames infeasible)
# using bounds for the followers profit (using Nogood function)
# INPUT
# m - model
# xOpt - current best solution
# OUTPUT
# iteraOpt - iteration in which the optimal solution is found
# iteraProve - number of iteration needed to prove optimality
# PLUS XOPT ALWAYS MAXIMAL WITH RESPECT TO THE MIP SOLUTION
# PLUS cut that avoid computing bi-level feasible solutions uninteresting (difference from improve2)
def EnumStopImproved4(m,x,z,xOpt,n,P,F,c_F,L,c_L,ti):
first_MIP_value = m.ObjVal
MIP_value_anterior = first_MIP_value
MIP_time = m.Runtime # running time of all solved MIPs
first_MIP_time = MIP_time
worst_MIP_time = MIP_time
first_MIP_nodes = m.NodeCount
worst_MIP_nodes = first_MIP_nodes
KP_time = 0 # running time of all solved KP
ratio =[((P[i]*1.)/F[i],i) for i in range(n)]
ratio.sort()
ratio.reverse()
# make xOpt maximal with respect to the mip solution
value = WeightInterdict(F,L,c_L,n)
xOpt = MakeMaximal(ratio,F,c_F,xOpt,n,m.getConstrs()[n*2+1].getAttr(GRB.Attr.Slack),L)
wmax,pmax = stopCritImproved(F,c_F,n,ratio,value,P)
#pmax,wmax = stopCrit(P,F,c_F,n,ratio)
#pmax = max(P)
#wmax = max(F)
yOpt,Profit,FollowerTime = ReactionFollower(n,P,F,c_F,xOpt,1) #maybe here we can opt for a differente algorithm
KP_time = KP_time+FollowerTime
x_new = xOpt
y_new = yOpt
s = 1
i=2
iteraOpt = 1 # step in which we find the best solution
upRow = 1 # do not update rowGen constraints rhs, if 1 update
oldProfit = Profit
while s==1:
x_new,m,LastMIPValue = NoImpFollower4(m,x,z,y_new,n,Profit,wmax,i,P,upRow,oldProfit) #to add the constraint sum_ {i: yi=1} xi >=1
upRow = 0
if m.Runtime>worst_MIP_time:
worst_MIP_time=m.Runtime
worst_MIP_nodes = m.NodeCount
MIP_time = MIP_time + m.Runtime
# fnovo = open('Evolucao4.txt','a')
# fnovo.write(str(m.ObjVal)+' ########## '+str(Profit)+'\n')
# fnovo.close()
if x_new ==[] or x_new==0 or time()-ti>3600: # the last model m is an infeasible problem
if x_new ==[]:
print '\nStopped in iteration ', i,'\nTHE LAST MODEL BUILT IS INFEASIBLE'
s = 0
LastMIPValue = MIP_value_anterior
else:
print '\nStopped in iteration ', i,'\nTIME LIMIT REACHED'
s = 0
LastMIPValue = MIP_value_anterior
xOpt = []
else:
if Profit + pmax <= LastMIPValue:
print '\nStopped in iteration ', i,'\nIT IS NOT POSSIBLE TO IMPROVE'
s = 0
else:
x_new = MakeMaximal(ratio,F,c_F,x_new,n,m.getConstrs()[n*2+1].getAttr(GRB.Attr.Slack),L)
y_new,profit_new,FollowerTime = ReactionFollower(n,P,F,c_F,x_new,1)
KP_time = KP_time+FollowerTime
if profit_new<Profit:
upRow = 1
oldProfit = Profit
Profit = profit_new
xOpt = x_new
yOpt = y_new
iteraOpt = i
MIP_value_anterior = LastMIPValue
i = i+1
# iteraProve = i-1
return xOpt,yOpt,Profit,iteraOpt,i-1,LastMIPValue,MIP_time,KP_time,first_MIP_time,first_MIP_nodes,worst_MIP_time,worst_MIP_nodes,first_MIP_value
# INPUT
# n - number of items
# P - list with the items profit
# F - list with the weights of the items
# c_F - follower capacity
# xOpt - strategy of the leader
# Bi - 1 if the variables are binary and 0 otherwise
# OUTPUT
# yOpt - best reaction of the follower
# profit
def ReactionFollower(n,P,F,c_F,xOpt,Bi):
Reaction = Model("ReactionF")
y = {}
for i in range(1,n+1):
# create decision variables and objective function
if Bi == 1:
y[i] = Reaction.addVar(obj=P[i-1],vtype="B",name="yf"+str(i))
else:
y[i] = Reaction.addVar(obj=P[i-1],vtype="C",name="yf"+str(i))
Reaction.update()
Reaction.addConstr(y[i]>=0)
Reaction.addConstr(xOpt[i-1]+y[i]<=1)
Reaction.update()
# constrain
Reaction.addConstr(quicksum(F[i]*y[i+1] for i in range(n))<=c_F)
Reaction.ModelSense = -1 # maximize
Reaction.update()
# solve
Reaction.optimize()
# view solution
profit = Reaction.ObjVal
yOpt=[]
for i in range(1,n+1):
yOpt.append(y[i].x)
return yOpt, profit,Reaction.Runtime
############## PRE-PROCESSING
def stopCrit(F,c_F,n,ratio,P):
C = 0
i = 0
while C <= c_F:
a,b = ratio[i]
C = C+ F[b]
i = i+1
wmax = 0
pmax = 0
for j in range(i-1,n): # items that may be critical
a,b = ratio[j]
if F[b]> wmax:
wmax = F[b]
if P[b] >pmax:
pmax = P[b]
return wmax,pmax
############## PRE-PROCESSING Improved
# ratio (of the Follower) - sorted by decreasing order
def stopCritImproved(F,c_F,n,ratio,value,P):
aux1,aux2 = ratio[-1]
if sum(F)>c_F+value+F[aux2]:# F[aux2]=wn - item with smaller ratio
C = 0
i = 0
while C <= c_F:
a,b = ratio[i]
C = C+ F[b]
i = i+1
wmax = F[b]
pmax = P[b]
a,b = ratio[i]
C = C+F[b]
while C-F[b]<c_F+value:
if F[b]> wmax:
wmax = F[b]
if P[b]>pmax:
pmax = P[b]
i = i+1
a,b = ratio[i]
C = C+F[b]
else:
wmax,pmax = stopCrit(F,c_F,n,ratio,P) #finalitem is the one above which the items are never critical
return wmax,pmax
### computation of the maximum weight of the follower that the leader can interdict (relaxed version)
def WeightInterdict(F,L,c_L,n):
value = 0
# sort the items by their value wi/vi
ratio1 = range(n)
ratio1.sort(key=lambda i: (-F[i]/float(L[i])))
R = c_L
aux3 = 0
while R>=0 and aux3<=n-1:
item = ratio1[aux3]
R = R-L[item]
value = value+F[item]
aux3 = aux3+1
if R < 0:
R = R +L[item]
value = value-F[item]+ F[item]*(c_L-R)/float(L[item])
return int(value)
def Run_CCLW(n,P,F,c_F,L,c_L):
ti = time()
####################################### FIRST ITERATIONS
setParam("TimeLimit", 3600)
# BUILT the first MIP
LeaderRelax = Model("Leader_Relax")
x,z = {},{}
z[1] = LeaderRelax.addVar(obj = c_F,vtype="C",name="z1")
LeaderRelax.update()
LeaderRelax.addConstr(z[1]>=0)
for j in range(1,n+1):
x[j] = LeaderRelax.addVar(vtype="B",name="x"+str(j))
z[j+1] = LeaderRelax.addVar(obj = 1,vtype="C",name="z"+str(j+1))
LeaderRelax.update()
LeaderRelax.addConstr(z[j+1]>=0)
# these dual constraints were updated to the new formulation of the followers primal problem
LeaderRelax.addConstr(F[j-1] * z[1] + z[j+1]>=P[j-1]*(1-x[j]))
LeaderRelax.update()
LeaderRelax.update()
# knapsasck constraint
LeaderRelax.addConstr(quicksum(L[i]*x[i+1] for i in range(n))<=c_L)
LeaderRelax.update()
# OBJECTIVE FUNCTION : add xi to the objective function in order to have a maximal solution
#LeaderRelax.setObjective(z[1]*c_F+sum(u[i] for i in range(1,n+1)),GRB.MINIMIZE)
LeaderRelax.ModelSense = 1 # to minimize
LeaderRelax.update()
# save problem - OPTIONAL
# LeaderRelax.write('Experiencia2.lp')
# solve
LeaderRelax.optimize()
# save solution file
# LeaderRelax.write('Experiencia2.sol')
if LeaderRelax.status != 9: # Optimization terminated because the time exceeded the value specified in the TimeLimit parameter
# x optimal:
xOpt=[]
for j in range(1,n+1):
xOpt.append(x[j].x)
################################# end of first iteration
if LeaderRelax.status != 9:
xOpt,yOpt,Profit,iteraOpt,iteraProve,LastMIPValue,MIP_time,KP_time,first_MIP_time,first_MIP_nodes,worst_MIP_time,worst_MIP_nodes,first_MIP_value = EnumStopImproved4(LeaderRelax,x,z,xOpt,n,P,F,c_F,L,c_L,ti)
else:
xOpt = []
timecpu = time()-ti
return xOpt,yOpt,Profit,iteraOpt,iteraProve,LastMIPValue,MIP_time,KP_time,first_MIP_time,first_MIP_nodes,worst_MIP_time,worst_MIP_nodes,first_MIP_value
if __name__ == "__main__":
seed(1)
################## Run a single instance #####################################
n = 35
P = [19, 87, 97, 22, 47, 22, 30, 5, 32, 54, 75, 70, 7, 76, 31, 61, 98, 27, 39, 91, 50, 10, 38, 25, 43, 37, 77, 18, 99, 9, 19, 97, 75, 26, 41]
F = [1, 96, 67, 90, 13, 74, 22, 86, 23, 63, 89, 25, 100, 76, 23, 56, 5, 47, 45, 39, 47, 11, 100, 27, 60, 76, 75, 37, 1, 81, 51, 30, 40, 48, 68]
L = [14, 85, 77, 26, 50, 45, 66, 79, 10, 3, 84, 44, 77, 1, 45, 73, 23, 95, 91, 4, 3, 55, 94, 39, 22, 43, 3, 23, 44, 50, 24, 24, 22, 46, 29]
c_L = 152
c_F = 162
xOpt,yOpt,Profit,iteraOpt,iteraProve,LastMIPValue,MIP_time,KP_time,first_MIP_time,first_MIP_nodes,worst_MIP_time,worst_MIP_nodes,first_MIP_value = Run_CCLW(n,P,F,c_F,L,c_L)
################## Run many instances #########################################
# fResults = open('Results.txt','a')
# fResults.write(' ins & FirstMIP.Obj & LastMIP.Obj & OPT & cpu_time\n')
# fResults.close()
#
# #for n in [10,15,20,25,30,35,40]:
# for n in [10,15,20]:
# for t in [1,2,3,4,5]:
# for i in range(1,101):
# # generate instance
# f = open('trash.txt','w')
# f.write(str(n)+'\n') # define n
# r = 100 # coeficients
# f.write(str(r)+'\n')
# f.write(str(t)+'\n')
# f.write(str(i)+'\n') # instance number
# f.write(str(1000)) # number of tests
# f.close()
# # Built instance using gen2.c for the follower
# os.system("./gen2<trash.txt")
# # read input file with the followers instances: test.in
# f = open('test.in','r')
# n = int(f.readline()) # number of items
# P_ord = [] # benefits - items values
# L_ord = [] # leader weights
# F_ord = [] # followers weights
# # formulate the knapsack problem for F
# for j in range(n):
# a,p,wf = map(int,f.readline().split())
# P_ord.append(p)
# F_ord.append(wf)
# wl = randint(1,r) # r is given in trash.txt
# L_ord.append(wl)
# # capacity
# c_F = int(f.readline())
# f.close()
# c_L = randint(c_F-10,c_F+10)
#
# xOpt,yOpt,Profit,iteraOpt,iteraProve,LastMIPValue,MIP_time,KP_time,first_MIP_time,first_MIP_nodes,worst_MIP_time,worst_MIP_nodes,first_MIP_value = Run_CCLW(n,P,F,c_F,L,c_L)
# fResults = open('Results.txt','a')
# fResults.write(str(i)+' & '+ str(first_MIP_value)+' & '+str(LastMIPValue)+' & '+str(Profit)+' & '+str(timecpu)+'\\\ \n\hline \n')
# fResults.close()