""" A Hybrid Tabu Search for Lot Sizing Problems Copyright (C) Joao Pedro Pedroso, LIACC/UP, 2005 these routines are dependent on: * the model specified in "fri-jt.mod" - variables named x[*,*,*] and y[*,*,*] are read from glpk's set of variables * data/instance format specied in "fridata.py" - requires .Prod, .Compat[p], .Period - uses the data class to produce the (ampl) model's ".dat" in a temporary file """ from glpkpi import * from glpk import * import random LOG = 0 EPS = 1.e-9 lpx_no_print() import resource def clock(): """clock() -> floating point number Return the CPU time in seconds (user + system) since the start of the process. This is done via a call to resource.getrusage, so it avoids the wraparound problems in time.clock().""" res = resource.getrusage(resource.RUSAGE_SELF) return res[0]+res[1] class LotSizing(lpx): def __init__(self, mod_file, data, CUT = {}): self.data = data # "produce ampl data file from data information" tmpdat = 'test.dat' # make a temporary ampl data file with this name self.data.mkampl(tmpdat, CUT) lpx.__init__(self, mod_file, tmpdat) def getsol(self): """copy lot sizing model's x, y, obj solution into dictionaries x[p,b,t] - (continuous) production y[p,b,t] - (0-1) setup z - objective *** NOTE *** supposing a (heuristic or exact) solution in available in self.sol """ x = {} y = {} for p in self.data.Prod: for b in self.data.Compat[p]: for t in self.data.Period: xname = "x[%s,%s,%d]" % (p,b,t) yname = "y[%s,%s,%d]" % (p,b,t) x[p,b,t] = self.sol[xname] y[p,b,t] = self.sol[yname] z = self.x[0] return x,y,z def printsol(self): """print solution in a easy-to-read format""" print "%10s%10s%10s%16s%16s" % ("product", "block", "period", "y", "x") for p in self.data.Prod: for b in self.data.Compat[p]: for t in self.data.Period: yname = "y[%s,%s,%d]" % (p,b,t) xname = "x[%s,%s,%d]" % (p,b,t) if self.sol[yname] != 0: print "%10s%10s%10i%16g%16g" % (p, b, t, self.sol[yname], self.sol[xname]) print "(remaining variables are 0)" print "objective:", self.x[0] def writesol(self,filename): """write x and y solution as text in 'filename'""" f = open(filename,'w') keys = self.sol.keys() keys.sort() for i in keys: print >>f, "%s = %r" % (i, self.sol[i]) f.close() def printstat(self): "print the current status of all the integer variables" print "bnds,\tkind,\t[lb,val,ub]" for i in self.ivars: print lpx_get_col_name(self.lp,i),\ lpx_get_col_type(self.lp,i), "\t", lpx_get_col_kind(self.lp,i),\ "[", lpx_get_col_lb(self.lp,i), lpx_mip_col_val(self.lp,i),\ lpx_get_col_ub(self.lp,i),"]" def update_sol(self, y): """update the current lp to reflect the choices in the y variable y is a dictionary: y[ """ for p in self.data.Prod: for b in self.data.Compat[p]: for t in self.data.Period: yname = "y[%s,%s,%d]" % (p,b,t) i = self.indices[yname] value = y[p,b,t] lpx_set_col_bnds(self.lp, i, LPX_FX, value, 0.) def destruct(self, KIND): """partial destruct a solution KIND is the glpk kind of variable to set (LPX_CV or LPX_IV) """ unfix = self.ivars[:] random.shuffle(unfix) alpha = random.random() nfix = int(len(self.ivars)*alpha) del unfix[nfix:] if LOG: print "alpha:", alpha print "unfix:", unfix for ivar in range(len(unfix)): i = unfix[ivar] lpx_set_col_bnds(self.lp, i, LPX_DB, self.lb[ivar], self.ub[ivar]) lpx_set_col_kind(self.lp, i, KIND) if LOG: print "made", unfix, " variables", KIND, "..." return unfix def rnd_fix(self): """fix a y variable at a random value; return its index on the lp""" ### this is rather inefficient, but later might be required... !!! # t = random.choice(self.data.Period) # p = random.choice(self.data.Prod) # b = random.choice(self.data.Compat[p]) # i = self.indices["y[%s,%s,%d]" % (p,b,t)] # if LOG: # print "changing setup for y[%s,%s,%d]" % (p,b,t) i = random.choice(self.ivars) ival = 1 - self.x[i] lpx_set_col_bnds(self.lp, i, LPX_FX, ival, 0.) if LOG: print "changed setup for %s to %d" % (lpx_get_col_name(self.lp,i), ival) return i def limited_solve(self, tmlim): """solve with branch & bound, with a timeout limit (note: might return before finding a feasible solution!) """ lpx_set_real_parm(self.lp, LPX_K_TMLIM, tmlim) lpx.solve(self) stat = lpx_mip_status(self.lp) if stat == LPX_I_OPT or stat == LPX_I_FEAS: # fix variables at the solution found for ivar in range(len(self.ivars)): i = self.ivars[ivar] value = self.x[i] lpx_set_col_bnds(self.lp, i, LPX_FX, value, 0.) else: print "failed, B&B could not find integer solution" raise Exception lpx_set_real_parm(self.lp, LPX_K_TMLIM, -1) def relax_and_fix(self): """relax-and-fix heuristic: make the variables: * for the current period: integer * for the previous periods: fixed * for the subsequent periods: continuous solve for all periods t=1...T """ for ivar in range(len(self.ivars)): i = self.ivars[ivar] lpx_set_col_bnds(self.lp, i, LPX_DB, self.lb[ivar], self.ub[ivar]) lpx_set_col_kind(self.lp, i, LPX_CV) if LOG: print "made all integer variables continuous..." # order vars for t in self.data.Period: var = [] per = ",%d]" % t for i in self.ivars: nm = self.names[i] if nm.find(per) != -1: var.append(i) for i in var: lpx_set_col_kind(self.lp, i, LPX_IV) if LOG: print 'made', var, 'integer' lpx.solve(self) for i in var: value = round(self.x[i]) if LOG: print "fixing %s (%f) to %f" % (self.names[i], self.x[i], value) lpx_set_col_bnds(self.lp, i, LPX_FX, value, 0.) ### lpx_set_col_kind(self.lp, i, LPX_CV) ??? def relax_and_fix_one_prod(self): """a more restricted relax-and-fix: make integer only the current periods'variables that concern ONE product (this seems to get many times the same sol...) """ for ivar in range(len(self.ivars)): i = self.ivars[ivar] lpx_set_col_bnds(self.lp, i, LPX_DB, self.lb[ivar], self.ub[ivar]) lpx_set_col_kind(self.lp, i, LPX_CV) if LOG: print "made all integer variables continuous..." # order vars for t in self.data.Period: rndprods = self.data.Prod[:] # make a product list in random order random.shuffle(rndprods) for p in rndprods: if LOG: print "period", t, "product", p, var = [] per = ",%d]" % t prd = "[%s," % p for i in self.ivars: nm = self.names[i] if nm.find(per) != -1 and nm.find(p) != -1: var.append(i) if LOG: print '; made', var, 'integer' for i in var: lpx_set_col_kind(self.lp, i, LPX_IV) if LOG: print self.names[i] lpx.solve(self) for i in var: value = round(self.x[i]) if LOG: print "fixing %s (%f) to %f" % (self.names[i], self.x[i], value) lpx_set_col_bnds(self.lp, i, LPX_FX, value, 0.) ### lpx_set_col_kind(self.lp, i, LPX_CV) ### // ??? required ??? def reconstruct_one_prod(self, unfix): """reconstruct, starting from a given solutiona more restricted relax-and-fix: make integer only the current periods'variables that concern ONE product (this seems to get many times the same sol...) """ # order vars for t in self.data.Period: rndprods = self.data.Prod[:] # make a product list in random order random.shuffle(rndprods) for p in rndprods: if LOG: print "period", t, "product", p, var = [] per = ",%d]" % t prd = "[%s," % p for i in unfix: nm = self.names[i] if nm.find(per) != -1 and nm.find(p) != -1: var.append(i) if var == []: continue if LOG: print '; made', var, 'integer' for i in var: lpx_set_col_kind(self.lp, i, LPX_IV) if LOG: print self.names[i] lpx.solve(self) for i in var: value = round(self.x[i]) if LOG: print "fixing %s (%f) to %f" % (self.names[i], self.x[i], value) lpx_set_col_bnds(self.lp, i, LPX_FX, value, 0.) ### lpx_set_col_kind(self.lp, i, LPX_CV) ### // ??? required ??? def relax_and_fix_one_var(self): """an even more restricted relax-and-fix: for each period, select ONLY ONE variable at a time for fixing; variable selection is done in a random order """ for ivar in range(len(self.ivars)): i = self.ivars[ivar] lpx_set_col_bnds(self.lp, i, LPX_DB, self.lb[ivar], self.ub[ivar]) lpx_set_col_kind(self.lp, i, LPX_CV) if LOG: print "made all integer variables continuous..." # order vars for t in self.data.Period: var = [] per = ",%d]" % t for i in self.ivars: nm = self.names[i] if nm.find(per) != -1: var.append(i) random.shuffle(var) for i in var: lpx_set_col_kind(self.lp, i, LPX_IV) lpx.solve(self) value = round(self.x[i]) if LOG: print "fixing %s (%f) to %f" % (self.names[i], self.x[i], value) lpx_set_col_bnds(self.lp, i, LPX_FX, value, 0.) ### lpx_set_col_kind(self.lp, i, LPX_CV) ??? class Stats: def __init__(self): self.niter = 0 self.ndestr = 0 self.nalltabu = 0 self.nreconstr = 0 self.niunset = 0 self.niset = 0 self.nidestr = 0 def __del__(self): print "\n\n\nRUN STATISTTICS" print "number of iterations:", self.niter print "number of destructions:", self.ndestr print "number of reconstructions:", self.nreconstr print "number of blocked situations:", self.nalltabu print "number of un-setup improvements on best:", self.niunset print "number of setup improvements on best:", self.niset print "number of destruct/reconstruct improvements on best:", self.nidestr stats = Stats() # quick and dirty def tabumove(prob, par, curr, zcurr, zbest, tabu, it): """check current neighbors for an improving solution arguments: prob -> a LotSizing instance (hence holding a lpx structure par -> ????? curr -> current solution (y[ in the MIP, 0-1 vector) zbest -> best found solution's objective value tabu -> list of variables which are tabu """ stats.niter +=1 cand = None # get the set of y=1 on the current solution # required info: # * the product corresponding to each variable # * the variables (machines) corresponding to each product (i.e., variable) # * if LOG: print "TC: starting candidate search; initial solution:" prob.printsol() move = None for t in prob.data.Period: if LOG: print "TC: PERIOD", t rndprods = prob.data.Prod[:] # make a product list in random order random.shuffle(rndprods) for p in rndprods: if LOG: print "TC: PRODUCT", p tabulen = int(len(prob.ivars) * random.random()) # !!! tabu length !!! # if p is being produced in this period setup = [] nosetup = [] for b in prob.data.Compat[p]: index = prob.indices["y[%s,%s,%d]" % (p,b,t)] if curr[p,b,t] == 1: setup.append(index) else: nosetup.append(index) for i in setup: lpx_set_col_bnds(prob.lp, i, LPX_FX, 0., 0.) prob.solve() x,y,z = prob.getsol() if LOG: print "unsetting up", prob.names[i], "obj:", z if z < zbest-EPS: # aspiration criterion stats.niunset += 1 if LOG: print "returning (best updated)" return y, z if (cand == None or z < zcand-EPS) and \ it - tabu[t,p] > tabulen: # p is not tabu cand = y.copy() zcand = z move = (i,0,None,None) tmove = t,p # tabu information if LOG: print "updated candidate: unsetting", i, "(%s)" % prob.names[i], "z=", zcand if zcand < zcurr: if LOG: print "returning (improved current solution)" return cand, zcand for j in nosetup: lpx_set_col_bnds(prob.lp, j, LPX_FX, 1., 0.) prob.solve() x,y,z = prob.getsol() if LOG: print "setting up", prob.names[j], "obj:", z if z < zbest-EPS: # aspiration criterion stats.niset += 1 if LOG: print "returning (-- best updated)" return y, z if (cand == None or z < zcand-EPS) and \ it - tabu[t,p] > tabulen: # p is not tabu cand = y.copy() zcand = z move = (i,0,j,1) tmove = t,p # tabu information if LOG: print "updated candidate: setting", j, "(%s)" % prob.names[j], "z=", zcand if zcand < zcurr: if LOG: print "returning (improved current solution)" return cand, zcand lpx_set_col_bnds(prob.lp, j, LPX_FX, 0., 0.) lpx_set_col_bnds(prob.lp, i, LPX_FX, 1., 0.) if move == None: "all moves were tabu, make a random move" stats.nalltabu += 1 index = prob.rnd_fix() p = random.choice(prob.data.Prod) b = random.choice(prob.data.Compat[p]) t = random.choice(prob.data.Period) i = prob.indices["y[%s,%s,%d]" % (p,b,t)] print "blocked; changing setup for y[%s,%s,%d]" % (p,b,t) ival = 1 - curr[p,b,t] tabu[(t,p)] = it lpx_set_col_bnds(prob.lp, i, LPX_FX, ival, 0.) x,y,z = prob.getsol() print "returnning; random move:", prob.names[i], "->", ival return y,z else: stats.ndestr += 1 i,ival,j,jval = move unfix = prob.destruct(LPX_CV) lpx_set_col_bnds(prob.lp, i, LPX_FX, ival, 0.) try: unfix.remove(i) except Exception: pass if j != None: lpx_set_col_bnds(prob.lp, j, LPX_FX, jval, 0.) try: unfix.remove(j) except Exception: pass # to do: remove all the tabu variables from unfix !!! prob.reconstruct_one_prod(unfix) stats.nreconstr += 1 if LOG: prob.printsol() print "candidate objective:", zcand x,cand,zcand = prob.getsol() if zcand < zbest-EPS: stats.nidestr += 1 if LOG: print "candidate objective after reconstruction:", zcand print "returning: reconstructed move" tabu[tmove] = it return cand,zcand #### i,ival,j,jval = move #### print "candidate is", prob.names[i], "->", ival #### lpx_set_col_bnds(prob.lp, i, LPX_FX, ival, 0.) #### if j != None: #### print "candidate is", prob.names[j], "->", jval #### lpx_set_col_bnds(prob.lp, j, LPX_FX, jval, 0.) #### #### tabu[tmove] = it #### print "returning: normal tabu move" #### return cand,zcand def tabusearch(prob, sol, zsol, zbest, tmlim): "tabu search main procedure" # initialize tabu list: tabu = {} for t in prob.data.Period: for p in prob.data.Prod: tabu[t,p] = -len(prob.ivars) print "initial solution:", zsol, "/", zbest best = sol.copy() it = 0 while clock() < tmlim: sol, zsol = tabumove(prob, None, sol, zsol, zbest, tabu, it) it += 1 print "iter", it, "current solution:", zsol, "/", zbest, "CPU", clock() if zsol < zbest: "improved, updated best solution" best = sol.copy() zbest = zsol print "best found solution:", best, zbest return best, zbest if __name__ == "__main__": import sys args = sys.argv args = ["", 1, 9, 600] if len(args) != 4: print "usage:", args[0], "seed instance-number tmlim(s)" sys.exit(-1) seed = int(args[1]) num = int(args[2]) tmlim = float(args[3]) print seed, num, tmlim random.seed(seed) import fridata # data = fridata.toy # data = fridata.jt # data = fridata.rnd data = fridata.inst[num] prob = LotSizing("fri-jt.mod", data) # prob.limited_solve(10) ### prob.printstat() ### raise Exception ### print "EXACT SOLUTION:" ### prob.solve() ### prob.printsol() ### ### print "relax_and_fix SOLUTION:" ### prob.relax_and_fix() ### prob.printsol() print "relax_and_fix_one_prod SOLUTION:" prob.relax_and_fix_one_prod() x,y,z = prob.getsol() ### prob.printstat() prob.printsol() print clock() ### sys.exit(0) ### raise Exception print "\n\n\n\n\ntabu search + reconstruct_one_prod:" sol,zsol = tabusearch(prob, y, z, z, tmlim) print "final solution:", sol print "final objective:", zsol sys.exit(0) ### LOG=0 ### print prob.indices ### for i in range(100): ### print "\n\n\n\n\ndestroy and B&B:" ### prob.destruct(LPX_IV) ### prob.rnd_fix() ### prob.limited_solve(20) ### x,y,z = prob.getsol() ### prob.printsol() ### raise Exception LOG=0 for i in range(1000): unfix = prob.destruct(LPX_CV) fixed = prob.rnd_fix() try: unfix.remove(fixed) except Exception: pass prob.reconstruct_one_prod(unfix) ### x,y,z = prob.getsol() prob.printsol() raise Exception print "reconstruct_one_prod SOLUTION:" prob.reconstruct_one_prod() x,y,z = prob.getsol() print "x=", x print "y=", y prob.printsol() print "reconstruct_one_prod SOLUTION:" prob.reconstruct_one_prod() x,y,z = prob.getsol() print "x=", x print "y=", y prob.printsol() ### tabusearch(prob, y, 10000, z)