From cc24fc87b5ffe4b74af4947f2c335d588fed08a0 Mon Sep 17 00:00:00 2001 From: Mohammed Ghannam Date: Tue, 2 Apr 2024 23:54:43 +0200 Subject: [PATCH] Fix formatting of examples & recipes (#835) --- examples/finished/atsp.py | 228 +++++++++---------- examples/finished/bpp.py | 69 +++--- examples/finished/diet.py | 104 ++++----- examples/finished/eoq_en.py | 52 ++--- examples/finished/even.py | 23 +- examples/finished/flp-benders.py | 66 +++--- examples/finished/flp.py | 65 +++--- examples/finished/gcp.py | 108 ++++----- examples/finished/gcp_fixed_k.py | 59 ++--- examples/finished/kmedian.py | 71 +++--- examples/finished/lo_wines.py | 34 +-- examples/finished/logical.py | 70 +++--- examples/finished/lotsizing_lazy.py | 119 +++++----- examples/finished/markowitz_soco.py | 38 ++-- examples/finished/mctransp.py | 130 +++++------ examples/finished/mkp.py | 25 +- examples/finished/pfs.py | 77 +++---- examples/finished/piecewise.py | 217 +++++++++--------- examples/finished/prodmix_soco.py | 45 ++-- examples/finished/rcs.py | 131 +++++------ examples/finished/read_tsplib.py | 134 ++++++----- examples/finished/ssa.py | 85 +++---- examples/finished/ssp.py | 27 ++- examples/finished/sudoku.py | 26 +-- examples/finished/transp.py | 53 ++--- examples/finished/transp_nofn.py | 32 +-- examples/finished/tsp.py | 73 +++--- examples/finished/weber_soco.py | 152 +++++++------ examples/tutorial/even.py | 58 ++--- examples/tutorial/logical.py | 41 ++-- examples/tutorial/puzzle.py | 6 +- examples/unfinished/cutstock.py | 2 +- examples/unfinished/eld.py | 4 +- examples/unfinished/flp_nonlinear.py | 4 +- examples/unfinished/flp_nonlinear_soco.py | 8 +- examples/unfinished/gpp.py | 3 +- examples/unfinished/kcenter.py | 3 +- examples/unfinished/kcenter_binary_search.py | 3 +- examples/unfinished/lotsizing.py | 2 + examples/unfinished/lotsizing_echelon.py | 2 + examples/unfinished/portfolio_soco.py | 4 +- examples/unfinished/staff_sched.py | 1 - examples/unfinished/staff_sched_mo.py | 6 +- examples/unfinished/tsp_flow.py | 5 +- examples/unfinished/tsp_lazy.py | 4 +- examples/unfinished/tsp_mo.py | 4 +- examples/unfinished/tsptw.py | 4 +- examples/unfinished/vrp.py | 67 +++--- examples/unfinished/vrp_lazy.py | 4 +- src/pyscipopt/recipes/README.md | 5 +- src/pyscipopt/recipes/nonlinear.py | 30 +-- 51 files changed, 1336 insertions(+), 1247 deletions(-) diff --git a/examples/finished/atsp.py b/examples/finished/atsp.py index 90b15616c..03597a498 100644 --- a/examples/finished/atsp.py +++ b/examples/finished/atsp.py @@ -1,5 +1,5 @@ ##@file atsp.py -#@brief solve the asymmetric traveling salesman problem +# @brief solve the asymmetric traveling salesman problem """ @@ -11,9 +11,10 @@ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ -from pyscipopt import Model, quicksum, multidict +from pyscipopt import Model, quicksum -def mtz(n,c): + +def mtz(n, c): """mtz: Miller-Tucker-Zemlin's model for the (asymmetric) traveling salesman problem (potential formulation) Parameters: @@ -24,30 +25,29 @@ def mtz(n,c): model = Model("atsp - mtz") - x,u = {},{} - for i in range(1,n+1): - u[i] = model.addVar(lb=0, ub=n-1, vtype="C", name="u(%s)"%i) - for j in range(1,n+1): + x, u = {}, {} + for i in range(1, n + 1): + u[i] = model.addVar(lb=0, ub=n - 1, vtype="C", name="u(%s)" % i) + for j in range(1, n + 1): if i != j: - x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j)) + x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j)) - for i in range(1,n+1): - model.addCons(quicksum(x[i,j] for j in range(1,n+1) if j != i) == 1, "Out(%s)"%i) - model.addCons(quicksum(x[j,i] for j in range(1,n+1) if j != i) == 1, "In(%s)"%i) + for i in range(1, n + 1): + model.addCons(quicksum(x[i, j] for j in range(1, n + 1) if j != i) == 1, "Out(%s)" % i) + model.addCons(quicksum(x[j, i] for j in range(1, n + 1) if j != i) == 1, "In(%s)" % i) - for i in range(1,n+1): - for j in range(2,n+1): + for i in range(1, n + 1): + for j in range(2, n + 1): if i != j: - model.addCons(u[i] - u[j] + (n-1)*x[i,j] <= n-2, "MTZ(%s,%s)"%(i,j)) + model.addCons(u[i] - u[j] + (n - 1) * x[i, j] <= n - 2, "MTZ(%s,%s)" % (i, j)) - model.setObjective(quicksum(c[i,j]*x[i,j] for (i,j) in x), "minimize") - - model.data = x,u - return model + model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), "minimize") + model.data = x, u + return model -def mtz_strong(n,c): +def mtz_strong(n, c): """mtz_strong: Miller-Tucker-Zemlin's model for the (asymmetric) traveling salesman problem (potential formulation, adding stronger constraints) Parameters: @@ -57,34 +57,34 @@ def mtz_strong(n,c): """ model = Model("atsp - mtz-strong") - - x,u = {},{} - for i in range(1,n+1): - u[i] = model.addVar(lb=0, ub=n-1, vtype="C", name="u(%s)"%i) - for j in range(1,n+1): + + x, u = {}, {} + for i in range(1, n + 1): + u[i] = model.addVar(lb=0, ub=n - 1, vtype="C", name="u(%s)" % i) + for j in range(1, n + 1): if i != j: - x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j)) + x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j)) - for i in range(1,n+1): - model.addCons(quicksum(x[i,j] for j in range(1,n+1) if j != i) == 1, "Out(%s)"%i) - model.addCons(quicksum(x[j,i] for j in range(1,n+1) if j != i) == 1, "In(%s)"%i) + for i in range(1, n + 1): + model.addCons(quicksum(x[i, j] for j in range(1, n + 1) if j != i) == 1, "Out(%s)" % i) + model.addCons(quicksum(x[j, i] for j in range(1, n + 1) if j != i) == 1, "In(%s)" % i) - for i in range(1,n+1): - for j in range(2,n+1): + for i in range(1, n + 1): + for j in range(2, n + 1): if i != j: - model.addCons(u[i] - u[j] + (n-1)*x[i,j] + (n-3)*x[j,i] <= n-2, "LiftedMTZ(%s,%s)"%(i,j)) + model.addCons(u[i] - u[j] + (n - 1) * x[i, j] + (n - 3) * x[j, i] <= n - 2, "LiftedMTZ(%s,%s)" % (i, j)) + + for i in range(2, n + 1): + model.addCons(-x[1, i] - u[i] + (n - 3) * x[i, 1] <= -2, name="LiftedLB(%s)" % i) + model.addCons(-x[i, 1] + u[i] + (n - 3) * x[1, i] <= n - 2, name="LiftedUB(%s)" % i) - for i in range(2,n+1): - model.addCons(-x[1,i] - u[i] + (n-3)*x[i,1] <= -2, name="LiftedLB(%s)"%i) - model.addCons(-x[i,1] + u[i] + (n-3)*x[1,i] <= n-2, name="LiftedUB(%s)"%i) + model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), "minimize") - model.setObjective(quicksum(c[i,j]*x[i,j] for (i,j) in x), "minimize") - - model.data = x,u + model.data = x, u return model -def scf(n,c): +def scf(n, c): """scf: single-commodity flow formulation for the (asymmetric) traveling salesman problem Parameters: - n: number of nodes @@ -93,40 +93,39 @@ def scf(n,c): """ model = Model("atsp - scf") - x,f = {},{} - for i in range(1,n+1): - for j in range(1,n+1): + x, f = {}, {} + for i in range(1, n + 1): + for j in range(1, n + 1): if i != j: - x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j)) - if i==1: - f[i,j] = model.addVar(lb=0, ub=n-1, vtype="C", name="f(%s,%s)"%(i,j)) + x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j)) + if i == 1: + f[i, j] = model.addVar(lb=0, ub=n - 1, vtype="C", name="f(%s,%s)" % (i, j)) else: - f[i,j] = model.addVar(lb=0, ub=n-2, vtype="C", name="f(%s,%s)"%(i,j)) + f[i, j] = model.addVar(lb=0, ub=n - 2, vtype="C", name="f(%s,%s)" % (i, j)) - for i in range(1,n+1): - model.addCons(quicksum(x[i,j] for j in range(1,n+1) if j != i) == 1, "Out(%s)"%i) - model.addCons(quicksum(x[j,i] for j in range(1,n+1) if j != i) == 1, "In(%s)"%i) + for i in range(1, n + 1): + model.addCons(quicksum(x[i, j] for j in range(1, n + 1) if j != i) == 1, "Out(%s)" % i) + model.addCons(quicksum(x[j, i] for j in range(1, n + 1) if j != i) == 1, "In(%s)" % i) - model.addCons(quicksum(f[1,j] for j in range(2,n+1)) == n-1, "FlowOut") + model.addCons(quicksum(f[1, j] for j in range(2, n + 1)) == n - 1, "FlowOut") - for i in range(2,n+1): - model.addCons(quicksum(f[j,i] for j in range(1,n+1) if j != i) - \ - quicksum(f[i,j] for j in range(1,n+1) if j != i) == 1, "FlowCons(%s)"%i) + for i in range(2, n + 1): + model.addCons(quicksum(f[j, i] for j in range(1, n + 1) if j != i) - \ + quicksum(f[i, j] for j in range(1, n + 1) if j != i) == 1, "FlowCons(%s)" % i) - for j in range(2,n+1): - model.addCons(f[1,j] <= (n-1)*x[1,j], "FlowUB(%s,%s)"%(1,j)) - for i in range(2,n+1): + for j in range(2, n + 1): + model.addCons(f[1, j] <= (n - 1) * x[1, j], "FlowUB(%s,%s)" % (1, j)) + for i in range(2, n + 1): if i != j: - model.addCons(f[i,j] <= (n-2)*x[i,j], "FlowUB(%s,%s)"%(i,j)) + model.addCons(f[i, j] <= (n - 2) * x[i, j], "FlowUB(%s,%s)" % (i, j)) - model.setObjective(quicksum(c[i,j]*x[i,j] for (i,j) in x), "minimize") + model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), "minimize") - model.data = x,f + model.data = x, f return model - -def mcf(n,c): +def mcf(n, c): """mcf: multi-commodity flow formulation for the (asymmetric) traveling salesman problem Parameters: - n: number of nodes @@ -135,48 +134,47 @@ def mcf(n,c): """ model = Model("mcf") - x,f = {},{} - for i in range(1,n+1): - for j in range(1,n+1): + x, f = {}, {} + for i in range(1, n + 1): + for j in range(1, n + 1): if i != j: - x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j)) + x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j)) if i != j and j != 1: - for k in range(2,n+1): + for k in range(2, n + 1): if i != k: - f[i,j,k] = model.addVar(ub=1, vtype="C", name="f(%s,%s,%s)"%(i,j,k)) + f[i, j, k] = model.addVar(ub=1, vtype="C", name="f(%s,%s,%s)" % (i, j, k)) - for i in range(1,n+1): - model.addCons(quicksum(x[i,j] for j in range(1,n+1) if j != i) == 1, "Out(%s)"%i) - model.addCons(quicksum(x[j,i] for j in range(1,n+1) if j != i) == 1, "In(%s)"%i) + for i in range(1, n + 1): + model.addCons(quicksum(x[i, j] for j in range(1, n + 1) if j != i) == 1, "Out(%s)" % i) + model.addCons(quicksum(x[j, i] for j in range(1, n + 1) if j != i) == 1, "In(%s)" % i) - for k in range(2,n+1): - model.addCons(quicksum(f[1,i,k] for i in range(2,n+1) if (1,i,k) in f) == 1, "FlowOut(%s)"%k) - model.addCons(quicksum(f[i,k,k] for i in range(1,n+1) if (i,k,k) in f) == 1, "FlowIn(%s)"%k) + for k in range(2, n + 1): + model.addCons(quicksum(f[1, i, k] for i in range(2, n + 1) if (1, i, k) in f) == 1, "FlowOut(%s)" % k) + model.addCons(quicksum(f[i, k, k] for i in range(1, n + 1) if (i, k, k) in f) == 1, "FlowIn(%s)" % k) - for i in range(2,n+1): + for i in range(2, n + 1): if i != k: - model.addCons(quicksum(f[j,i,k] for j in range(1,n+1) if (j,i,k) in f) == \ - quicksum(f[i,j,k] for j in range(1,n+1) if (i,j,k) in f), - "FlowCons(%s,%s)"%(i,k)) + model.addCons(quicksum(f[j, i, k] for j in range(1, n + 1) if (j, i, k) in f) == \ + quicksum(f[i, j, k] for j in range(1, n + 1) if (i, j, k) in f), + "FlowCons(%s,%s)" % (i, k)) - for (i,j,k) in f: - model.addCons(f[i,j,k] <= x[i,j], "FlowUB(%s,%s,%s)"%(i,j,k)) + for (i, j, k) in f: + model.addCons(f[i, j, k] <= x[i, j], "FlowUB(%s,%s,%s)" % (i, j, k)) - model.setObjective(quicksum(c[i,j]*x[i,j] for (i,j) in x), "minimize") + model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), "minimize") - model.data = x,f + model.data = x, f return model - def sequence(arcs): """sequence: make a list of cities to visit, from set of arcs""" succ = {} - for (i,j) in arcs: + for (i, j) in arcs: succ[i] = j - curr = 1 # first node being visited + curr = 1 # first node being visited sol = [curr] - for i in range(len(arcs)-2): + for i in range(len(arcs) - 2): curr = succ[curr] sol.append(curr) return sol @@ -184,85 +182,85 @@ def sequence(arcs): if __name__ == "__main__": n = 5 - c = { (1,1):0, (1,2):1989, (1,3):102, (1,4):102, (1,5):103, - (2,1):104, (2,2):0, (2,3):11, (2,4):104, (2,5):108, - (3,1):107, (3,2):108, (3,3):0, (3,4):19, (3,5):102, - (4,1):109, (4,2):102, (4,3):107, (4,4):0, (4,5):15, - (5,1):13, (5,2):103, (5,3):104, (5,4):101, (5,5):0, + c = {(1, 1): 0, (1, 2): 1989, (1, 3): 102, (1, 4): 102, (1, 5): 103, + (2, 1): 104, (2, 2): 0, (2, 3): 11, (2, 4): 104, (2, 5): 108, + (3, 1): 107, (3, 2): 108, (3, 3): 0, (3, 4): 19, (3, 5): 102, + (4, 1): 109, (4, 2): 102, (4, 3): 107, (4, 4): 0, (4, 5): 15, + (5, 1): 13, (5, 2): 103, (5, 3): 104, (5, 4): 101, (5, 5): 0, } - model = mtz(n,c) - model.hideOutput() # silent mode + model = mtz(n, c) + model.hideOutput() # silent mode model.optimize() cost = model.getObjVal() print() print("Miller-Tucker-Zemlin's model:") print("Optimal value:", cost) - #model.printAttr("X") + # model.printAttr("X") for v in model.getVars(): if model.getVal(v) > 0.001: print(v.name, "=", model.getVal(v)) - x,u = model.data - sol = [i for (p,i) in sorted([(int(model.getVal(u[i])+.5),i) for i in range(1,n+1)])] + x, u = model.data + sol = [i for (p, i) in sorted([(int(model.getVal(u[i]) + .5), i) for i in range(1, n + 1)])] print(sol) - arcs = [(i,j) for (i,j) in x if model.getVal(x[i,j]) > .5] + arcs = [(i, j) for (i, j) in x if model.getVal(x[i, j]) > .5] sol = sequence(arcs) print(sol) # assert cost == 5 - model = mtz_strong(n,c) - model.hideOutput() # silent mode + model = mtz_strong(n, c) + model.hideOutput() # silent mode model.optimize() cost = model.getObjVal() print() print("Miller-Tucker-Zemlin's model with stronger constraints:") - print("Optimal value:",cost) - #model.printAttr("X") + print("Optimal value:", cost) + # model.printAttr("X") for v in model.getVars(): if model.getVal(v) > 0.001: print(v.name, "=", model.getVal(v)) - x,u = model.data - sol = [i for (p,i) in sorted([(int(model.getVal(u[i])+.5),i) for i in range(1,n+1)])] + x, u = model.data + sol = [i for (p, i) in sorted([(int(model.getVal(u[i]) + .5), i) for i in range(1, n + 1)])] print(sol) - arcs = [(i,j) for (i,j) in x if model.getVal(x[i,j]) > .5] + arcs = [(i, j) for (i, j) in x if model.getVal(x[i, j]) > .5] sol = sequence(arcs) print(sol) # assert cost == 5 - model = scf(n,c) - model.hideOutput() # silent mode + model = scf(n, c) + model.hideOutput() # silent mode model.optimize() cost = model.getObjVal() print() print("Single-commodity flow formulation:") - print("Optimal value:",cost) - #model.printAttr("X") + print("Optimal value:", cost) + # model.printAttr("X") for v in model.getVars(): if model.getVal(v) > 0.001: print(v.name, "=", model.getVal(v)) - x,f = model.data - arcs = [(i,j) for (i,j) in x if model.getVal(x[i,j]) > .5] + x, f = model.data + arcs = [(i, j) for (i, j) in x if model.getVal(x[i, j]) > .5] sol = sequence(arcs) print(sol) # assert cost == 5 - model = mcf(n,c) - model.hideOutput() # silent mode + model = mcf(n, c) + model.hideOutput() # silent mode model.optimize() cost = model.getObjVal() print() print("Multi-commodity flow formulation:") - print("Optimal value:",cost) - #model.printAttr("X") + print("Optimal value:", cost) + # model.printAttr("X") for v in model.getVars(): - if model.getVal(v)>0.001: + if model.getVal(v) > 0.001: print(v.name, "=", model.getVal(v)) - x,f = model.data - arcs = [(i,j) for (i,j) in x if model.getVal(x[i,j]) > .5] + x, f = model.data + arcs = [(i, j) for (i, j) in x if model.getVal(x[i, j]) > .5] sol = sequence(arcs) print(sol) # assert cost == 5 diff --git a/examples/finished/bpp.py b/examples/finished/bpp.py index 26a438cf8..3c0333305 100644 --- a/examples/finished/bpp.py +++ b/examples/finished/bpp.py @@ -1,5 +1,5 @@ ##@file bpp.py -#@brief use SCIP for solving the bin packing problem. +# @brief use SCIP for solving the bin packing problem. """ The instance of the bin packing problem is represented by the two lists of n items of sizes and quantity s=(s_i). @@ -13,28 +13,29 @@ from pyscipopt import Model, quicksum -def FFD(s,B): + +def FFD(s, B): """First Fit Decreasing heuristics for the Bin Packing Problem. Parameters: - s: list with item widths - B: bin capacity Returns a list of lists with bin compositions. """ - remain = [B] # keep list of empty space per bin - sol = [[]] # a list ot items (i.e., sizes) on each used bin - for item in sorted(s,reverse=True): - for (j,free) in enumerate(remain): + remain = [B] # keep list of empty space per bin + sol = [[]] # a list ot items (i.e., sizes) on each used bin + for item in sorted(s, reverse=True): + for (j, free) in enumerate(remain): if free >= item: remain[j] -= item sol[j].append(item) break - else: #does not fit in any bin + else: # does not fit in any bin sol.append([item]) - remain.append(B-item) + remain.append(B - item) return sol -def bpp(s,B): +def bpp(s, B): """bpp: Martello and Toth's model to solve the bin packing problem. Parameters: - s: list with item widths @@ -42,44 +43,44 @@ def bpp(s,B): Returns a model, ready to be solved. """ n = len(s) - U = len(FFD(s,B)) # upper bound of the number of bins + U = len(FFD(s, B)) # upper bound of the number of bins model = Model("bpp") # setParam("MIPFocus",1) - x,y = {},{} + x, y = {}, {} for i in range(n): for j in range(U): - x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j)) + x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j)) for j in range(U): - y[j] = model.addVar(vtype="B", name="y(%s)"%j) + y[j] = model.addVar(vtype="B", name="y(%s)" % j) # assignment constraints for i in range(n): - model.addCons(quicksum(x[i,j] for j in range(U)) == 1, "Assign(%s)"%i) + model.addCons(quicksum(x[i, j] for j in range(U)) == 1, "Assign(%s)" % i) # bin capacity constraints for j in range(U): - model.addCons(quicksum(s[i]*x[i,j] for i in range(n)) <= B*y[j], "Capac(%s)"%j) + model.addCons(quicksum(s[i] * x[i, j] for i in range(n)) <= B * y[j], "Capac(%s)" % j) # tighten assignment constraints for j in range(U): for i in range(n): - model.addCons(x[i,j] <= y[j], "Strong(%s,%s)"%(i,j)) + model.addCons(x[i, j] <= y[j], "Strong(%s,%s)" % (i, j)) # tie breaking constraints - for j in range(U-1): - model.addCons(y[j] >= y[j+1],"TieBrk(%s)"%j) + for j in range(U - 1): + model.addCons(y[j] >= y[j + 1], "TieBrk(%s)" % j) # SOS constraints for i in range(n): - model.addConsSOS1([x[i,j] for j in range(U)]) + model.addConsSOS1([x[i, j] for j in range(U)]) model.setObjective(quicksum(y[j] for j in range(U)), "minimize") - model.data = x,y + model.data = x, y return model -def solveBinPacking(s,B): +def solveBinPacking(s, B): """solveBinPacking: use an IP model to solve the in Packing Problem. Parameters: @@ -89,15 +90,15 @@ def solveBinPacking(s,B): Returns a solution: list of lists, each of which with the items in a roll. """ n = len(s) - U = len(FFD(s,B)) # upper bound of the number of bins - model = bpp(s,B) - x,y = model.data + U = len(FFD(s, B)) # upper bound of the number of bins + model = bpp(s, B) + x, y = model.data model.optimize() bins = [[] for i in range(U)] - for (i,j) in x: - if model.getVal(x[i,j]) > .5: + for (i, j) in x: + if model.getVal(x[i, j]) > .5: bins[j].append(s[i]) for i in range(bins.count([])): bins.remove([]) @@ -108,29 +109,31 @@ def solveBinPacking(s,B): import random -def DiscreteUniform(n=10,LB=1,UB=99,B=100): + + +def DiscreteUniform(n=10, LB=1, UB=99, B=100): """DiscreteUniform: create random, uniform instance for the bin packing problem.""" B = 100 - s = [0]*n + s = [0] * n for i in range(n): - s[i] = random.randint(LB,UB) - return s,B + s[i] = random.randint(LB, UB) + return s, B if __name__ == "__main__": random.seed(256) - s,B = DiscreteUniform() + s, B = DiscreteUniform() print("items:", s) print("bin size:", B) - ffd = FFD(s,B) + ffd = FFD(s, B) print("\n\n\n FFD heuristic:") print("Solution:") print(ffd) print(len(ffd), "bins") print("\n\n\n IP formulation:") - bins = solveBinPacking(s,B) + bins = solveBinPacking(s, B) print("Solution:") print(bins) print(len(bins), "bins") diff --git a/examples/finished/diet.py b/examples/finished/diet.py index 478b640e8..53fbf83f4 100644 --- a/examples/finished/diet.py +++ b/examples/finished/diet.py @@ -1,12 +1,13 @@ ##@file diet.py -#@brief model for the modern diet problem +# @brief model for the modern diet problem """ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ # todo: can we make it work as "from pyscipopt import *"? from pyscipopt import Model, quicksum, multidict -def diet(F,N,a,b,c,d): + +def diet(F, N, a, b, c, d): """diet -- model for the modern diet problem Parameters: - F: set of foods @@ -21,84 +22,85 @@ def diet(F,N,a,b,c,d): model = Model("modern diet") # Create variables - x,y,z = {},{},{} + x, y, z = {}, {}, {} for j in F: - x[j] = model.addVar(vtype="I", name="x(%s)"%j) - y[j] = model.addVar(vtype="B", name="y(%s)"%j) + x[j] = model.addVar(vtype="I", name="x(%s)" % j) + y[j] = model.addVar(vtype="B", name="y(%s)" % j) for i in N: - z[i] = model.addVar(lb=a[i], ub=b[i], name="z(%s)"%j) + z[i] = model.addVar(lb=a[i], ub=b[i], name="z(%s)" % j) v = model.addVar(vtype="C", name="v") # Constraints: for i in N: - model.addCons(quicksum(d[j][i]*x[j] for j in F) == z[i], name="Nutr(%s)"%i) + model.addCons(quicksum(d[j][i] * x[j] for j in F) == z[i], name="Nutr(%s)" % i) - model.addCons(quicksum(c[j]*x[j] for j in F) == v, name="Cost") + model.addCons(quicksum(c[j] * x[j] for j in F) == v, name="Cost") for j in F: - model.addCons(y[j] <= x[j], name="Eat(%s)"%j) + model.addCons(y[j] <= x[j], name="Eat(%s)" % j) # Objective: model.setObjective(quicksum(y[j] for j in F), "maximize") - model.data = x,y,z,v + model.data = x, y, z, v return model def make_inst(): """make_inst: prepare data for the diet model""" - F,c,d = multidict({ # cost # composition - "QPounder" : [ 1.84, {"Cal":510, "Carbo":34, "Protein":28, - "VitA":15, "VitC": 6, "Calc":30, "Iron":20}], - "McLean" : [ 2.19, {"Cal":370, "Carbo":35, "Protein":24, "VitA":15, - "VitC": 10, "Calc":20, "Iron":20}], - "Big Mac" : [ 1.84, {"Cal":500, "Carbo":42, "Protein":25, - "VitA": 6, "VitC": 2, "Calc":25, "Iron":20}], - "FFilet" : [ 1.44, {"Cal":370, "Carbo":38, "Protein":14, - "VitA": 2, "VitC": 0, "Calc":15, "Iron":10}], - "Chicken" : [ 2.29, {"Cal":400, "Carbo":42, "Protein":31, - "VitA": 8, "VitC": 15, "Calc":15, "Iron": 8}], - "Fries" : [ .77, {"Cal":220, "Carbo":26, "Protein": 3, - "VitA": 0, "VitC": 15, "Calc": 0, "Iron": 2}], - "McMuffin" : [ 1.29, {"Cal":345, "Carbo":27, "Protein":15, - "VitA": 4, "VitC": 0, "Calc":20, "Iron":15}], - "1% LFMilk": [ .60, {"Cal":110, "Carbo":12, "Protein": 9, - "VitA":10, "VitC": 4, "Calc":30, "Iron": 0}], - "OrgJuice" : [ .72, {"Cal": 80, "Carbo":20, "Protein": 1, - "VitA": 2, "VitC":120, "Calc": 2, "Iron": 2}], - }) - - N,a,b = multidict({ # min,max intake - "Cal" : [ 2000, None ], - "Carbo" : [ 350, 375 ], - "Protein" : [ 55, None ], - "VitA" : [ 100, None ], - "VitC" : [ 100, None ], - "Calc" : [ 100, None ], - "Iron" : [ 100, None ], - }) - - return F,N,a,b,c,d + F, c, d = multidict({ # cost # composition + "QPounder": [1.84, {"Cal": 510, "Carbo": 34, "Protein": 28, + "VitA": 15, "VitC": 6, "Calc": 30, "Iron": 20}], + "McLean": [2.19, {"Cal": 370, "Carbo": 35, "Protein": 24, "VitA": 15, + "VitC": 10, "Calc": 20, "Iron": 20}], + "Big Mac": [1.84, {"Cal": 500, "Carbo": 42, "Protein": 25, + "VitA": 6, "VitC": 2, "Calc": 25, "Iron": 20}], + "FFilet": [1.44, {"Cal": 370, "Carbo": 38, "Protein": 14, + "VitA": 2, "VitC": 0, "Calc": 15, "Iron": 10}], + "Chicken": [2.29, {"Cal": 400, "Carbo": 42, "Protein": 31, + "VitA": 8, "VitC": 15, "Calc": 15, "Iron": 8}], + "Fries": [.77, {"Cal": 220, "Carbo": 26, "Protein": 3, + "VitA": 0, "VitC": 15, "Calc": 0, "Iron": 2}], + "McMuffin": [1.29, {"Cal": 345, "Carbo": 27, "Protein": 15, + "VitA": 4, "VitC": 0, "Calc": 20, "Iron": 15}], + "1% LFMilk": [.60, {"Cal": 110, "Carbo": 12, "Protein": 9, + "VitA": 10, "VitC": 4, "Calc": 30, "Iron": 0}], + "OrgJuice": [.72, {"Cal": 80, "Carbo": 20, "Protein": 1, + "VitA": 2, "VitC": 120, "Calc": 2, "Iron": 2}], + }) + + N, a, b = multidict({ # min,max intake + "Cal": [2000, None], + "Carbo": [350, 375], + "Protein": [55, None], + "VitA": [100, None], + "VitC": [100, None], + "Calc": [100, None], + "Iron": [100, None], + }) + + return F, N, a, b, c, d if __name__ == "__main__": - F,N,a,b,c,d = make_inst() + F, N, a, b, c, d = make_inst() - for b["Cal"] in [None,3500,3000,2500]: + for b["Cal"] in [None, 3500, 3000, 2500]: print("\nDiet for a maximum of {0} calories".format(b["Cal"] if b["Cal"] != None else "unlimited")) - model = diet(F,N,a,b,c,d) - model.hideOutput() # silent mode + model = diet(F, N, a, b, c, d) + model.hideOutput() # silent mode model.optimize() - print("Optimal value:",model.getObjVal()) - x,y,z,v = model.data + print("Optimal value:", model.getObjVal()) + x, y, z, v = model.data for j in x: if model.getVal(x[j]) > 0: - print("{0:30s}: {1:3.1f} dishes --> {2:4.2f} added to objective".format(j,model.getVal(x[j]),model.getVal(y[j]))) - print("amount spent:",model.getObjVal()) + print("{0:30s}: {1:3.1f} dishes --> {2:4.2f} added to objective".format(j, model.getVal(x[j]), + model.getVal(y[j]))) + print("amount spent:", model.getObjVal()) print("amount of nutrients:") for i in z: - print("{0:30s}: {1:4.2f}".format(i,model.getVal(z[i]))) + print("{0:30s}: {1:4.2f}".format(i, model.getVal(z[i]))) diff --git a/examples/finished/eoq_en.py b/examples/finished/eoq_en.py index e16483675..b52bae1d0 100644 --- a/examples/finished/eoq_en.py +++ b/examples/finished/eoq_en.py @@ -1,5 +1,5 @@ ##@file eoq_en.py -#@brief piecewise linear model to the multi-item economic ordering quantity problem. +# @brief piecewise linear model to the multi-item economic ordering quantity problem. """ Approach: use a convex combination formulation. @@ -7,7 +7,8 @@ """ from pyscipopt import Model, quicksum, multidict -def eoq(I,F,h,d,w,W,a0,aK,K): + +def eoq(I, F, h, d, w, W, a0, aK, K): """eoq -- multi-item capacitated economic ordering quantity model Parameters: - I: set of items @@ -23,54 +24,53 @@ def eoq(I,F,h,d,w,W,a0,aK,K): """ # construct points for piecewise-linear relation, store in a,b - a,b = {},{} - delta = float(aK-a0)/K + a, b = {}, {} + delta = float(aK - a0) / K for i in I: for k in range(K): - T = a0 + delta*k - a[i,k] = T # abscissa: cycle time - b[i,k] = F[i]/T + h[i]*d[i]*T/2. # ordinate: (convex) cost for this cycle time + T = a0 + delta * k + a[i, k] = T # abscissa: cycle time + b[i, k] = F[i] / T + h[i] * d[i] * T / 2. # ordinate: (convex) cost for this cycle time model = Model("multi-item, capacitated EOQ") - x,c,w_ = {},{},{} + x, c, w_ = {}, {}, {} for i in I: - x[i] = model.addVar(vtype="C", name="x(%s)"%i) # cycle time for item i - c[i] = model.addVar(vtype="C", name="c(%s)"%i) # total cost for item i + x[i] = model.addVar(vtype="C", name="x(%s)" % i) # cycle time for item i + c[i] = model.addVar(vtype="C", name="c(%s)" % i) # total cost for item i for k in range(K): - w_[i,k] = model.addVar(ub=1, vtype="C", name="w(%s,%s)"%(i,k)) #todo ?? + w_[i, k] = model.addVar(ub=1, vtype="C", name="w(%s,%s)" % (i, k)) # todo ?? for i in I: - model.addCons(quicksum(w_[i,k] for k in range(K)) == 1) - model.addCons(quicksum(a[i,k]*w_[i,k] for k in range(K)) == x[i]) - model.addCons(quicksum(b[i,k]*w_[i,k] for k in range(K)) == c[i]) + model.addCons(quicksum(w_[i, k] for k in range(K)) == 1) + model.addCons(quicksum(a[i, k] * w_[i, k] for k in range(K)) == x[i]) + model.addCons(quicksum(b[i, k] * w_[i, k] for k in range(K)) == c[i]) - model.addCons(quicksum(w[i]*d[i]*x[i] for i in I) <= W) + model.addCons(quicksum(w[i] * d[i] * x[i] for i in I) <= W) model.setObjective(quicksum(c[i] for i in I), "minimize") - model.data = x,w + model.data = x, w return model - if __name__ == "__main__": # multiple item EOQ - I,F,h,d,w = multidict( - {1:[300,10,10,20], - 2:[300,10,30,40], - 3:[300,10,50,10]} - ) + I, F, h, d, w = multidict( + {1: [300, 10, 10, 20], + 2: [300, 10, 30, 40], + 3: [300, 10, 50, 10]} + ) W = 2000 K = 1000 - a0,aK = 0.1,10 - model = eoq(I,F,h,d,w,W,a0,aK,K) + a0, aK = 0.1, 10 + model = eoq(I, F, h, d, w, W, a0, aK, K) model.optimize() - x,w = model.data + x, w = model.data EPS = 1.e-6 for v in x: if model.getVal(x[v]) >= EPS: - print(x[v].name,"=",model.getVal(x[v])) + print(x[v].name, "=", model.getVal(x[v])) print("Optimal value:", model.getObjVal()) diff --git a/examples/finished/even.py b/examples/finished/even.py index 037c11aae..657ec86f2 100644 --- a/examples/finished/even.py +++ b/examples/finished/even.py @@ -1,6 +1,5 @@ ##@file finished/even.py -#@brief model to decide whether argument is even or odd - +# @brief model to decide whether argument is even or odd ################################################################################ @@ -23,18 +22,19 @@ from pyscipopt import Model verbose = False -sdic = {0:"even",1:"odd"} +sdic = {0: "even", 1: "odd"} + def parity(number): try: assert number == int(round(number)) m = Model() m.hideOutput() - + ### variables are non-negative by default since 0 is the default lb. ### To allow for negative values, give None as lower bound ### (None means -infinity as lower bound and +infinity as upper bound) - x = m.addVar("x", vtype="I", lb=None, ub=None) #ub=None is default + x = m.addVar("x", vtype="I", lb=None, ub=None) # ub=None is default n = m.addVar("n", vtype="I", lb=None) s = m.addVar("s", vtype="B") @@ -42,18 +42,18 @@ def parity(number): ### if x is set by default as non-negative and number is negative: ### there is no feasible solution (trivial) but the program ### does not highlight which constraints conflict. - m.addCons(x==number) + m.addCons(x == number) - m.addCons(s == x-2*n) + m.addCons(s == x - 2 * n) m.setObjective(s) m.optimize() assert m.getStatus() == "optimal" if verbose: for v in m.getVars(): - print("%s %d" % (v,m.getVal(v))) + print("%s %d" % (v, m.getVal(v))) print("%d%%2 == %d?" % (m.getVal(x), m.getVal(s))) - print(m.getVal(s) == m.getVal(x)%2) + print(m.getVal(s) == m.getVal(x) % 2) xval = m.getVal(x) sval = m.getVal(s) @@ -62,10 +62,13 @@ def parity(number): except (AssertionError, TypeError): print("%s is neither even nor odd!" % number.__repr__()) + if __name__ == "__main__": import sys from ast import literal_eval as leval - example_values = [0, 1, 1.5, "hallo welt", 20, 25, -101, -15., -10, -int(2**31), int(2**31-1), int(2**63)-1] + + example_values = [0, 1, 1.5, "hallo welt", 20, 25, -101, -15., -10, -int(2 ** 31), int(2 ** 31 - 1), + int(2 ** 63) - 1] try: try: n = leval(sys.argv[1]) diff --git a/examples/finished/flp-benders.py b/examples/finished/flp-benders.py index c6a947f8a..8a2d95894 100644 --- a/examples/finished/flp-benders.py +++ b/examples/finished/flp-benders.py @@ -1,5 +1,5 @@ ##@file flp-benders.py -#@brief model for solving the capacitated facility location problem using Benders' decomposition +# @brief model for solving the capacitated facility location problem using Benders' decomposition """ minimize the total (weighted) travel cost from n customers to some facilities with fixed costs and capacities. @@ -7,9 +7,9 @@ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ from pyscipopt import Model, quicksum, multidict, SCIP_PARAMSETTING -import pdb -def flp(I,J,d,M,f,c): + +def flp(I, J, d, M, f, c): """flp -- model for the capacitated facility location problem Parameters: - I: set of customers @@ -27,54 +27,53 @@ def flp(I,J,d,M,f,c): # creating the problem y = {} for j in J: - y[j] = master.addVar(vtype="B", name="y(%s)"%j) + y[j] = master.addVar(vtype="B", name="y(%s)" % j) master.setObjective( - quicksum(f[j]*y[j] for j in J), + quicksum(f[j] * y[j] for j in J), "minimize") master.data = y # creating the subproblem - x,y = {},{} + x, y = {}, {} for j in J: - y[j] = subprob.addVar(vtype="B", name="y(%s)"%j) + y[j] = subprob.addVar(vtype="B", name="y(%s)" % j) for i in I: - x[i,j] = subprob.addVar(vtype="C", name="x(%s,%s)"%(i,j)) + x[i, j] = subprob.addVar(vtype="C", name="x(%s,%s)" % (i, j)) for i in I: - subprob.addCons(quicksum(x[i,j] for j in J) == d[i], "Demand(%s)"%i) + subprob.addCons(quicksum(x[i, j] for j in J) == d[i], "Demand(%s)" % i) for j in M: - subprob.addCons(quicksum(x[i,j] for i in I) <= M[j]*y[j], "Capacity(%s)"%i) + subprob.addCons(quicksum(x[i, j] for i in I) <= M[j] * y[j], "Capacity(%s)" % i) - for (i,j) in x: - subprob.addCons(x[i,j] <= d[i]*y[j], "Strong(%s,%s)"%(i,j)) + for (i, j) in x: + subprob.addCons(x[i, j] <= d[i] * y[j], "Strong(%s,%s)" % (i, j)) subprob.setObjective( - quicksum(c[i,j]*x[i,j] for i in I for j in J), + quicksum(c[i, j] * x[i, j] for i in I for j in J), "minimize") - subprob.data = x,y + subprob.data = x, y return master, subprob def make_data(): """creates example data set""" - I,d = multidict({1:80, 2:270, 3:250, 4:160, 5:180}) # demand - J,M,f = multidict({1:[500,1000], 2:[500,1000], 3:[500,1000]}) # capacity, fixed costs - c = {(1,1):4, (1,2):6, (1,3):9, # transportation costs - (2,1):5, (2,2):4, (2,3):7, - (3,1):6, (3,2):3, (3,3):4, - (4,1):8, (4,2):5, (4,3):3, - (5,1):10, (5,2):8, (5,3):4, + I, d = multidict({1: 80, 2: 270, 3: 250, 4: 160, 5: 180}) # demand + J, M, f = multidict({1: [500, 1000], 2: [500, 1000], 3: [500, 1000]}) # capacity, fixed costs + c = {(1, 1): 4, (1, 2): 6, (1, 3): 9, # transportation costs + (2, 1): 5, (2, 2): 4, (2, 3): 7, + (3, 1): 6, (3, 2): 3, (3, 3): 4, + (4, 1): 8, (4, 2): 5, (4, 3): 3, + (5, 1): 10, (5, 2): 8, (5, 3): 4, } - return I,J,d,M,f,c - + return I, J, d, M, f, c if __name__ == "__main__": - I,J,d,M,f,c = make_data() - master, subprob = flp(I,J,d,M,f,c) + I, J, d, M, f, c = make_data() + master, subprob = flp(I, J, d, M, f, c) # initializing the default Benders' decomposition with the subproblem master.setPresolve(SCIP_PARAMSETTING.OFF) master.setBoolParam("misc/allowstrongdualreds", False) @@ -92,7 +91,7 @@ def make_data(): facilities = [j for j in y if master.getVal(y[j]) > EPS] x, suby = subprob.data - edges = [(i,j) for (i,j) in x if subprob.getVal(x[i,j]) > EPS] + edges = [(i, j) for (i, j) in x if subprob.getVal(x[i, j]) > EPS] print("Optimal value:", master.getObjVal()) print("Facilities at nodes:", facilities) @@ -105,24 +104,25 @@ def make_data(): # the solution will be lost master.freeBendersSubproblems() - try: # plot the result using networkx and matplotlib + try: # plot the result using networkx and matplotlib import networkx as NX import matplotlib.pyplot as P + P.clf() G = NX.Graph() other = [j for j in y if j not in facilities] - customers = ["c%s"%i for i in d] + customers = ["c%s" % i for i in d] G.add_nodes_from(facilities) G.add_nodes_from(other) G.add_nodes_from(customers) - for (i,j) in edges: - G.add_edge("c%s"%i,j) + for (i, j) in edges: + G.add_edge("c%s" % i, j) position = NX.drawing.layout.spring_layout(G) - NX.draw(G,position,node_color="y",nodelist=facilities) - NX.draw(G,position,node_color="g",nodelist=other) - NX.draw(G,position,node_color="b",nodelist=customers) + NX.draw(G, position, node_color="y", nodelist=facilities) + NX.draw(G, position, node_color="g", nodelist=other) + NX.draw(G, position, node_color="b", nodelist=customers) P.show() except ImportError: print("install 'networkx' and 'matplotlib' for plotting") diff --git a/examples/finished/flp.py b/examples/finished/flp.py index 50b716079..2c5391db7 100644 --- a/examples/finished/flp.py +++ b/examples/finished/flp.py @@ -1,5 +1,5 @@ ##@file flp.py -#@brief model for solving the capacitated facility location problem +# @brief model for solving the capacitated facility location problem """ minimize the total (weighted) travel cost from n customers to some facilities with fixed costs and capacities. @@ -8,7 +8,8 @@ """ from pyscipopt import Model, quicksum, multidict -def flp(I,J,d,M,f,c): + +def flp(I, J, d, M, f, c): """flp -- model for the capacitated facility location problem Parameters: - I: set of customers @@ -22,76 +23,76 @@ def flp(I,J,d,M,f,c): model = Model("flp") - x,y = {},{} + x, y = {}, {} for j in J: - y[j] = model.addVar(vtype="B", name="y(%s)"%j) + y[j] = model.addVar(vtype="B", name="y(%s)" % j) for i in I: - x[i,j] = model.addVar(vtype="C", name="x(%s,%s)"%(i,j)) + x[i, j] = model.addVar(vtype="C", name="x(%s,%s)" % (i, j)) for i in I: - model.addCons(quicksum(x[i,j] for j in J) == d[i], "Demand(%s)"%i) + model.addCons(quicksum(x[i, j] for j in J) == d[i], "Demand(%s)" % i) for j in M: - model.addCons(quicksum(x[i,j] for i in I) <= M[j]*y[j], "Capacity(%s)"%i) + model.addCons(quicksum(x[i, j] for i in I) <= M[j] * y[j], "Capacity(%s)" % i) - for (i,j) in x: - model.addCons(x[i,j] <= d[i]*y[j], "Strong(%s,%s)"%(i,j)) + for (i, j) in x: + model.addCons(x[i, j] <= d[i] * y[j], "Strong(%s,%s)" % (i, j)) model.setObjective( - quicksum(f[j]*y[j] for j in J) + - quicksum(c[i,j]*x[i,j] for i in I for j in J), + quicksum(f[j] * y[j] for j in J) + + quicksum(c[i, j] * x[i, j] for i in I for j in J), "minimize") - model.data = x,y + model.data = x, y return model def make_data(): """creates example data set""" - I,d = multidict({1:80, 2:270, 3:250, 4:160, 5:180}) # demand - J,M,f = multidict({1:[500,1000], 2:[500,1000], 3:[500,1000]}) # capacity, fixed costs - c = {(1,1):4, (1,2):6, (1,3):9, # transportation costs - (2,1):5, (2,2):4, (2,3):7, - (3,1):6, (3,2):3, (3,3):4, - (4,1):8, (4,2):5, (4,3):3, - (5,1):10, (5,2):8, (5,3):4, + I, d = multidict({1: 80, 2: 270, 3: 250, 4: 160, 5: 180}) # demand + J, M, f = multidict({1: [500, 1000], 2: [500, 1000], 3: [500, 1000]}) # capacity, fixed costs + c = {(1, 1): 4, (1, 2): 6, (1, 3): 9, # transportation costs + (2, 1): 5, (2, 2): 4, (2, 3): 7, + (3, 1): 6, (3, 2): 3, (3, 3): 4, + (4, 1): 8, (4, 2): 5, (4, 3): 3, + (5, 1): 10, (5, 2): 8, (5, 3): 4, } - return I,J,d,M,f,c - + return I, J, d, M, f, c if __name__ == "__main__": - I,J,d,M,f,c = make_data() - model = flp(I,J,d,M,f,c) + I, J, d, M, f, c = make_data() + model = flp(I, J, d, M, f, c) model.optimize() EPS = 1.e-6 - x,y = model.data - edges = [(i,j) for (i,j) in x if model.getVal(x[i,j]) > EPS] + x, y = model.data + edges = [(i, j) for (i, j) in x if model.getVal(x[i, j]) > EPS] facilities = [j for j in y if model.getVal(y[j]) > EPS] print("Optimal value:", model.getObjVal()) print("Facilities at nodes:", facilities) print("Edges:", edges) - try: # plot the result using networkx and matplotlib + try: # plot the result using networkx and matplotlib import networkx as NX import matplotlib.pyplot as P + P.clf() G = NX.Graph() other = [j for j in y if j not in facilities] - customers = ["c%s"%i for i in d] + customers = ["c%s" % i for i in d] G.add_nodes_from(facilities) G.add_nodes_from(other) G.add_nodes_from(customers) - for (i,j) in edges: - G.add_edge("c%s"%i,j) + for (i, j) in edges: + G.add_edge("c%s" % i, j) position = NX.drawing.layout.spring_layout(G) - NX.draw(G,position,node_color="y",nodelist=facilities) - NX.draw(G,position,node_color="g",nodelist=other) - NX.draw(G,position,node_color="b",nodelist=customers) + NX.draw(G, position, node_color="y", nodelist=facilities) + NX.draw(G, position, node_color="g", nodelist=other) + NX.draw(G, position, node_color="b", nodelist=customers) P.show() except ImportError: print("install 'networkx' and 'matplotlib' for plotting") diff --git a/examples/finished/gcp.py b/examples/finished/gcp.py index 9b46f282a..a07d2c44b 100644 --- a/examples/finished/gcp.py +++ b/examples/finished/gcp.py @@ -1,12 +1,13 @@ ##@file gcp.py -#@brief model for the graph coloring problem +# @brief model for the graph coloring problem """ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ -from pyscipopt import Model, quicksum, multidict +from pyscipopt import Model, quicksum -def gcp(V,E,K): + +def gcp(V, E, K): """gcp -- model for minimizing the number of colors in a graph Parameters: - V: set/list of nodes in the graph @@ -15,18 +16,18 @@ def gcp(V,E,K): Returns a model, ready to be solved. """ model = Model("gcp") - x,y = {},{} + x, y = {}, {} for k in range(K): - y[k] = model.addVar(vtype="B", name="y(%s)"%k) + y[k] = model.addVar(vtype="B", name="y(%s)" % k) for i in V: - x[i,k] = model.addVar(vtype="B", name="x(%s,%s)"%(i,k)) + x[i, k] = model.addVar(vtype="B", name="x(%s,%s)" % (i, k)) for i in V: - model.addCons(quicksum(x[i,k] for k in range(K)) == 1, "AssignColor(%s)"%i) + model.addCons(quicksum(x[i, k] for k in range(K)) == 1, "AssignColor(%s)" % i) - for (i,j) in E: + for (i, j) in E: for k in range(K): - model.addCons(x[i,k] + x[j,k] <= y[k], "NotSameColor(%s,%s,%s)"%(i,j,k)) + model.addCons(x[i, k] + x[j, k] <= y[k], "NotSameColor(%s,%s,%s)" % (i, j, k)) model.setObjective(quicksum(y[k] for k in range(K)), "minimize") @@ -34,7 +35,7 @@ def gcp(V,E,K): return model -def gcp_low(V,E,K): +def gcp_low(V, E, K): """gcp_low -- model for minimizing the number of colors in a graph (use colors with low indices) Parameters: @@ -44,22 +45,21 @@ def gcp_low(V,E,K): Returns a model, ready to be solved. """ model = Model("gcp - low colors") - x,y = {},{} + x, y = {}, {} for k in range(K): - y[k] = model.addVar(vtype="B", name="y(%s)"%k) + y[k] = model.addVar(vtype="B", name="y(%s)" % k) for i in V: - x[i,k] = model.addVar(vtype="B", name="x(%s,%s)"%(i,k)) - + x[i, k] = model.addVar(vtype="B", name="x(%s,%s)" % (i, k)) for i in V: - model.addCons(quicksum(x[i,k] for k in range(K)) == 1, "AssignColor(%s)" % i) + model.addCons(quicksum(x[i, k] for k in range(K)) == 1, "AssignColor(%s)" % i) - for (i,j) in E: + for (i, j) in E: for k in range(K): - model.addCons(x[i,k] + x[j,k] <= y[k], "NotSameColor(%s,%s,%s)"%(i,j,k)) + model.addCons(x[i, k] + x[j, k] <= y[k], "NotSameColor(%s,%s,%s)" % (i, j, k)) - for k in range(K-1): - model.addCons(y[k] >= y[k+1], "LowColor(%s)"%k) + for k in range(K - 1): + model.addCons(y[k] >= y[k + 1], "LowColor(%s)" % k) model.setObjective(quicksum(y[k] for k in range(K)), "minimize") @@ -67,7 +67,7 @@ def gcp_low(V,E,K): return model -def gcp_sos(V,E,K): +def gcp_sos(V, E, K): """gcp_sos -- model for minimizing the number of colors in a graph (use sos type 1 constraints) Parameters: @@ -78,22 +78,22 @@ def gcp_sos(V,E,K): """ model = Model("gcp - sos constraints") - x,y = {},{} + x, y = {}, {} for k in range(K): - y[k] = model.addVar(vtype="B", name="y(%s)"%k) + y[k] = model.addVar(vtype="B", name="y(%s)" % k) for i in V: - x[i,k] = model.addVar(vtype="B", name="x(%s,%s)"%(i,k)) + x[i, k] = model.addVar(vtype="B", name="x(%s,%s)" % (i, k)) for i in V: - model.addCons(quicksum(x[i,k] for k in range(K)) == 1, "AssignColor(%s)" % i) - model.addConsSOS1([x[i,k] for k in range(K)]) + model.addCons(quicksum(x[i, k] for k in range(K)) == 1, "AssignColor(%s)" % i) + model.addConsSOS1([x[i, k] for k in range(K)]) - for (i,j) in E: + for (i, j) in E: for k in range(K): - model.addCons(x[i,k] + x[j,k] <= y[k], "NotSameColor(%s,%s,%s)"%(i,j,k)) + model.addCons(x[i, k] + x[j, k] <= y[k], "NotSameColor(%s,%s,%s)" % (i, j, k)) - for k in range(K-1): - model.addCons(y[k] >= y[k+1], "LowColor(%s)"%k) + for k in range(K - 1): + model.addCons(y[k] >= y[k + 1], "LowColor(%s)" % k) model.setObjective(quicksum(y[k] for k in range(K)), "minimize") @@ -102,25 +102,27 @@ def gcp_sos(V,E,K): import random -def make_data(n,prob): + + +def make_data(n, prob): """make_data: prepare data for a random graph Parameters: - n: number of vertices - prob: probability of existence of an edge, for each pair of vertices Returns a tuple with a list of vertices and a list edges. """ - V = range(1,n+1) - E = [(i,j) for i in V for j in V if i < j and random.random() < prob] - return V,E + V = range(1, n + 1) + E = [(i, j) for i in V for j in V if i < j and random.random() < prob] + return V, E if __name__ == "__main__": random.seed(1) - V,E = make_data(20,.5) - K = 10 # upper bound to the number of colors - print("n,K=",len(V),K) + V, E = make_data(20, .5) + K = 10 # upper bound to the number of colors + print("n,K=", len(V), K) - model = gcp_low(V,E,K) + model = gcp_low(V, E, K) model.optimize() print("Optimal value:", model.getObjVal()) x = model.data @@ -128,36 +130,38 @@ def make_data(n,prob): color = {} for i in V: for k in range(K): - if model.getVal(x[i,k]) > 0.5: + if model.getVal(x[i, k]) > 0.5: color[i] = k - print("colors:",color) + print("colors:", color) + + import time, sys - import time,sys - models = [gcp,gcp_low,gcp_sos] + models = [gcp, gcp_low, gcp_sos] cpu = {} - N = 25 # number of observations + N = 25 # number of observations print("#size\t%s\t%s\t%s" % tuple(m.__name__ for m in models)) for size in range(250): - print(size,"\t",) + print(size, "\t", ) K = size for prob in [0.1]: for m in models: name = m.__name__ - if not (name,size-1,prob) in cpu or cpu[name,size-1,prob] < 100: #cpu.has_key((name,size-1,prob)) - cpu[name,size,prob] = 0. + if not (name, size - 1, prob) in cpu or cpu[ + name, size - 1, prob] < 100: # cpu.has_key((name,size-1,prob)) + cpu[name, size, prob] = 0. for t in range(N): tinit = time.time() random.seed(t) - V,E = make_data(size,prob) - model = m(V,E,K) - model.hideOutput() # silent mode + V, E = make_data(size, prob) + model = m(V, E, K) + model.hideOutput() # silent mode model.optimize() assert model.getObjVal() >= 0 and model.getObjVal() <= K tend = time.time() - cpu[name,size,prob] += tend - tinit - cpu[name,size,prob] /= N + cpu[name, size, prob] += tend - tinit + cpu[name, size, prob] /= N else: - cpu[name,size,prob] = "-" - print(cpu[name,size,prob],"\t",) + cpu[name, size, prob] = "-" + print(cpu[name, size, prob], "\t", ) print() sys.stdout.flush() diff --git a/examples/finished/gcp_fixed_k.py b/examples/finished/gcp_fixed_k.py index e539f5638..df403e0a8 100644 --- a/examples/finished/gcp_fixed_k.py +++ b/examples/finished/gcp_fixed_k.py @@ -1,12 +1,13 @@ ##@file gcp_fixed_k.py -#@brief solve the graph coloring problem with fixed-k model +# @brief solve the graph coloring problem with fixed-k model """ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ -from pyscipopt import Model, quicksum, multidict +from pyscipopt import Model, quicksum -def gcp_fixed_k(V,E,K): + +def gcp_fixed_k(V, E, K): """gcp_fixed_k -- model for minimizing number of bad edges in coloring a graph Parameters: - V: set/list of nodes in the graph @@ -16,27 +17,27 @@ def gcp_fixed_k(V,E,K): """ model = Model("gcp - fixed k") - x,z = {},{} + x, z = {}, {} for i in V: for k in range(K): - x[i,k] = model.addVar(vtype="B", name="x(%s,%s)"%(i,k)) - for (i,j) in E: - z[i,j] = model.addVar(vtype="B", name="z(%s,%s)"%(i,j)) + x[i, k] = model.addVar(vtype="B", name="x(%s,%s)" % (i, k)) + for (i, j) in E: + z[i, j] = model.addVar(vtype="B", name="z(%s,%s)" % (i, j)) for i in V: - model.addCons(quicksum(x[i,k] for k in range(K)) == 1, "AssignColor(%s)" % i) + model.addCons(quicksum(x[i, k] for k in range(K)) == 1, "AssignColor(%s)" % i) - for (i,j) in E: + for (i, j) in E: for k in range(K): - model.addCons(x[i,k] + x[j,k] <= 1 + z[i,j], "BadEdge(%s,%s,%s)"%(i,j,k)) + model.addCons(x[i, k] + x[j, k] <= 1 + z[i, j], "BadEdge(%s,%s,%s)" % (i, j, k)) - model.setObjective(quicksum(z[i,j] for (i,j) in E), "minimize") + model.setObjective(quicksum(z[i, j] for (i, j) in E), "minimize") - model.data = x,z + model.data = x, z return model -def solve_gcp(V,E): +def solve_gcp(V, E): """solve_gcp -- solve the graph coloring problem with bisection and fixed-k model Parameters: - V: set/list of nodes in the graph @@ -46,19 +47,19 @@ def solve_gcp(V,E): LB = 0 UB = len(V) color = {} - while UB-LB > 1: - K = int((UB+LB) / 2) - gcp = gcp_fixed_k(V,E,K) + while UB - LB > 1: + K = int((UB + LB) / 2) + gcp = gcp_fixed_k(V, E, K) # gcp.Params.OutputFlag = 0 # silent mode - #gcp.Params.Cutoff = .1 + # gcp.Params.Cutoff = .1 gcp.setObjlimit(0.1) gcp.optimize() status = gcp.getStatus() if status == "optimal": - x,z = gcp.data + x, z = gcp.data for i in V: for k in range(K): - if gcp.getVal(x[i,k]) > 0.5: + if gcp.getVal(x[i, k]) > 0.5: color[i] = k break # else: @@ -67,26 +68,28 @@ def solve_gcp(V,E): else: LB = K - return UB,color + return UB, color import random -def make_data(n,prob): + + +def make_data(n, prob): """make_data: prepare data for a random graph Parameters: - n: number of vertices - prob: probability of existence of an edge, for each pair of vertices Returns a tuple with a list of vertices and a list edges. """ - V = range(1,n+1) - E = [(i,j) for i in V for j in V if i < j and random.random() < prob] - return V,E + V = range(1, n + 1) + E = [(i, j) for i in V for j in V if i < j and random.random() < prob] + return V, E if __name__ == "__main__": random.seed(1) - V,E = make_data(75,.25) + V, E = make_data(75, .25) - K,color = solve_gcp(V,E) - print("minimum number of colors:",K) - print("solution:",color) + K, color = solve_gcp(V, E) + print("minimum number of colors:", K) + print("solution:", color) diff --git a/examples/finished/kmedian.py b/examples/finished/kmedian.py index db1910810..d6c4a338b 100644 --- a/examples/finished/kmedian.py +++ b/examples/finished/kmedian.py @@ -1,5 +1,5 @@ ##@file kmedian.py -#@brief model for solving the k-median problem. +# @brief model for solving the k-median problem. """ minimize the total (weighted) travel cost for servicing a set of customers from k facilities. @@ -8,9 +8,11 @@ """ import math import random -from pyscipopt import Model, quicksum, multidict -def kmedian(I,J,c,k): +from pyscipopt import Model, quicksum + + +def kmedian(I, J, c, k): """kmedian -- minimize total cost of servicing customers from k facilities Parameters: - I: set of customers @@ -21,73 +23,74 @@ def kmedian(I,J,c,k): """ model = Model("k-median") - x,y = {},{} + x, y = {}, {} for j in J: - y[j] = model.addVar(vtype="B", name="y(%s)"%j) + y[j] = model.addVar(vtype="B", name="y(%s)" % j) for i in I: - x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j)) + x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j)) for i in I: - model.addCons(quicksum(x[i,j] for j in J) == 1, "Assign(%s)"%i) + model.addCons(quicksum(x[i, j] for j in J) == 1, "Assign(%s)" % i) for j in J: - model.addCons(x[i,j] <= y[j], "Strong(%s,%s)"%(i,j)) + model.addCons(x[i, j] <= y[j], "Strong(%s,%s)" % (i, j)) model.addCons(quicksum(y[j] for j in J) == k, "Facilities") - model.setObjective(quicksum(c[i,j]*x[i,j] for i in I for j in J), "minimize") - model.data = x,y + model.setObjective(quicksum(c[i, j] * x[i, j] for i in I for j in J), "minimize") + model.data = x, y return model -def distance(x1,y1,x2,y2): +def distance(x1, y1, x2, y2): """return distance of two points""" - return math.sqrt((x2-x1)**2 + (y2-y1)**2) + return math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2) -def make_data(n,m,same=True): +def make_data(n, m, same=True): """creates example data set""" if same == True: I = range(n) J = range(m) - x = [random.random() for i in range(max(m,n))] # positions of the points in the plane - y = [random.random() for i in range(max(m,n))] + x = [random.random() for i in range(max(m, n))] # positions of the points in the plane + y = [random.random() for i in range(max(m, n))] else: I = range(n) - J = range(n,n+m) - x = [random.random() for i in range(n+m)] # positions of the points in the plane - y = [random.random() for i in range(n+m)] + J = range(n, n + m) + x = [random.random() for i in range(n + m)] # positions of the points in the plane + y = [random.random() for i in range(n + m)] c = {} for i in I: for j in J: - c[i,j] = distance(x[i],y[i],x[j],y[j]) + c[i, j] = distance(x[i], y[i], x[j], y[j]) - return I,J,c,x,y + return I, J, c, x, y if __name__ == "__main__": - import sys + random.seed(67) n = 200 m = n - I,J,c,x_pos,y_pos = make_data(n,m,same=True) + I, J, c, x_pos, y_pos = make_data(n, m, same=True) k = 20 - model = kmedian(I,J,c,k) + model = kmedian(I, J, c, k) # model.Params.Threads = 1 model.optimize() EPS = 1.e-6 - x,y = model.data - edges = [(i,j) for (i,j) in x if model.getVal(x[i,j]) > EPS] + x, y = model.data + edges = [(i, j) for (i, j) in x if model.getVal(x[i, j]) > EPS] facilities = [j for j in y if model.getVal(y[j]) > EPS] - print("Optimal value:",model.getObjVal()) + print("Optimal value:", model.getObjVal()) print("Selected facilities:", facilities) print("Edges:", edges) - print("max c:", max([c[i,j] for (i,j) in edges])) + print("max c:", max([c[i, j] for (i, j) in edges])) - try: # plot the result using networkx and matplotlib + try: # plot the result using networkx and matplotlib import networkx as NX import matplotlib.pyplot as P + P.clf() G = NX.Graph() @@ -97,16 +100,16 @@ def make_data(n,m,same=True): G.add_nodes_from(facilities) G.add_nodes_from(client) G.add_nodes_from(other) - for (i,j) in edges: - G.add_edge(i,j) + for (i, j) in edges: + G.add_edge(i, j) position = {} for i in range(len(x_pos)): - position[i] = (x_pos[i],y_pos[i]) + position[i] = (x_pos[i], y_pos[i]) - NX.draw(G,position,with_labels=False,node_color="w",nodelist=facilities) - NX.draw(G,position,with_labels=False,node_color="c",nodelist=other,node_size=50) - NX.draw(G,position,with_labels=False,node_color="g",nodelist=client,node_size=50) + NX.draw(G, position, with_labels=False, node_color="w", nodelist=facilities) + NX.draw(G, position, with_labels=False, node_color="c", nodelist=other, node_size=50) + NX.draw(G, position, with_labels=False, node_color="g", nodelist=client, node_size=50) P.show() except ImportError: print("install 'networkx' and 'matplotlib' for plotting") diff --git a/examples/finished/lo_wines.py b/examples/finished/lo_wines.py index 9b9e3ff92..172301c84 100644 --- a/examples/finished/lo_wines.py +++ b/examples/finished/lo_wines.py @@ -1,5 +1,5 @@ ##@file lo_wines.py -#@brief Simple SCIP example of linear programming. +# @brief Simple SCIP example of linear programming. """ It solves the same instance as lo_wines_simple.py: @@ -16,40 +16,40 @@ """ from pyscipopt import Model, quicksum, SCIP_PARAMSETTING -#Initialize model +# Initialize model model = Model("Wine blending") model.setPresolve(SCIP_PARAMSETTING.OFF) -Inventory = {"Alfrocheiro":60, "Baga":60, "Castelao":30} +Inventory = {"Alfrocheiro": 60, "Baga": 60, "Castelao": 30} Grapes = Inventory.keys() -Profit = {"Dry":15, "Medium":18, "Sweet":30} +Profit = {"Dry": 15, "Medium": 18, "Sweet": 30} Blends = Profit.keys() Use = { - ("Alfrocheiro","Dry"):2, - ("Alfrocheiro","Medium"):1, - ("Alfrocheiro","Sweet"):1, - ("Baga","Dry"):1, - ("Baga","Medium"):2, - ("Baga","Sweet"):1, - ("Castelao","Dry"):0, - ("Castelao","Medium"):0, - ("Castelao","Sweet"):1 - } + ("Alfrocheiro", "Dry"): 2, + ("Alfrocheiro", "Medium"): 1, + ("Alfrocheiro", "Sweet"): 1, + ("Baga", "Dry"): 1, + ("Baga", "Medium"): 2, + ("Baga", "Sweet"): 1, + ("Castelao", "Dry"): 0, + ("Castelao", "Medium"): 0, + ("Castelao", "Sweet"): 1 +} # Create variables x = {} for j in Blends: - x[j] = model.addVar(vtype="C", name="x(%s)"%j) + x[j] = model.addVar(vtype="C", name="x(%s)" % j) # Create constraints c = {} for i in Grapes: - c[i] = model.addCons(quicksum(Use[i,j]*x[j] for j in Blends) <= Inventory[i], name="Use(%s)"%i) + c[i] = model.addCons(quicksum(Use[i, j] * x[j] for j in Blends) <= Inventory[i], name="Use(%s)" % i) # Objective -model.setObjective(quicksum(Profit[j]*x[j] for j in Blends), "maximize") +model.setObjective(quicksum(Profit[j] * x[j] for j in Blends), "maximize") model.optimize() diff --git a/examples/finished/logical.py b/examples/finished/logical.py index e84f37893..b28cd4123 100644 --- a/examples/finished/logical.py +++ b/examples/finished/logical.py @@ -1,5 +1,5 @@ ##@file finished/logical.py -#@brief Tutorial example on how to use AND/OR/XOR constraints +# @brief Tutorial example on how to use AND/OR/XOR constraints from pyscipopt import Model from pyscipopt import quicksum @@ -17,7 +17,8 @@ """ -def printFunc(name,m): + +def printFunc(name, m): """prints results""" print("* %s *" % name) objSet = bool(m.getObjective().terms.keys()) @@ -29,54 +30,55 @@ def printFunc(name,m): print("%s: %d" % (v, round(m.getVal(v)))) print("\n") -# AND + +# AND model = Model() model.hideOutput() -x = model.addVar("x","B") -y = model.addVar("y","B") -z = model.addVar("z","B") -r = model.addVar("r","B") -model.addConsAnd([x,y,z],r) -model.addCons(x==1) -model.setObjective(r,sense="minimize") +x = model.addVar("x", "B") +y = model.addVar("y", "B") +z = model.addVar("z", "B") +r = model.addVar("r", "B") +model.addConsAnd([x, y, z], r) +model.addCons(x == 1) +model.setObjective(r, sense="minimize") model.optimize() -printFunc("AND",model) +printFunc("AND", model) # OR model = Model() model.hideOutput() -x = model.addVar("x","B") -y = model.addVar("y","B") -z = model.addVar("z","B") -r = model.addVar("r","B") -model.addConsOr([x,y,z],r) -model.addCons(x==0) -model.setObjective(r,sense="maximize") +x = model.addVar("x", "B") +y = model.addVar("y", "B") +z = model.addVar("z", "B") +r = model.addVar("r", "B") +model.addConsOr([x, y, z], r) +model.addCons(x == 0) +model.setObjective(r, sense="maximize") model.optimize() -printFunc("OR",model) +printFunc("OR", model) # XOR (r as boolean, standard) model = Model() model.hideOutput() -x = model.addVar("x","B") -y = model.addVar("y","B") -z = model.addVar("z","B") +x = model.addVar("x", "B") +y = model.addVar("y", "B") +z = model.addVar("z", "B") r = True -model.addConsXor([x,y,z],r) -model.addCons(x==1) +model.addConsXor([x, y, z], r) +model.addCons(x == 1) model.optimize() -printFunc("Standard XOR (as boolean)",model) +printFunc("Standard XOR (as boolean)", model) # XOR (r as variable, custom) model = Model() model.hideOutput() -x = model.addVar("x","B") -y = model.addVar("y","B") -z = model.addVar("z","B") -r = model.addVar("r","B") -n = model.addVar("n","I") #auxiliary -model.addCons(r+quicksum([x,y,z]) == 2*n) -model.addCons(x==0) -model.setObjective(r,sense="maximize") +x = model.addVar("x", "B") +y = model.addVar("y", "B") +z = model.addVar("z", "B") +r = model.addVar("r", "B") +n = model.addVar("n", "I") # auxiliary +model.addCons(r + quicksum([x, y, z]) == 2 * n) +model.addCons(x == 0) +model.setObjective(r, sense="maximize") model.optimize() -printFunc("Custom XOR (as variable)",model) +printFunc("Custom XOR (as variable)", model) diff --git a/examples/finished/lotsizing_lazy.py b/examples/finished/lotsizing_lazy.py index 39ed7e22f..934b16348 100644 --- a/examples/finished/lotsizing_lazy.py +++ b/examples/finished/lotsizing_lazy.py @@ -1,5 +1,5 @@ ##@file lotsizing_lazy.py -#@brief solve the single-item lot-sizing problem. +# @brief solve the single-item lot-sizing problem. """ Approaches: - sils: solve the problem using the standard formulation @@ -9,44 +9,45 @@ """ from pyscipopt import Model, Conshdlr, quicksum, multidict, SCIP_RESULT, SCIP_PRESOLTIMING, SCIP_PROPTIMING + class Conshdlr_sils(Conshdlr): def addcut(self, checkonly, sol): - D,Ts = self.data - y,x,I = self.model.data + D, Ts = self.data + y, x, I = self.model.data cutsadded = False for ell in Ts: lhs = 0 - S,L = [],[] - for t in range(1,ell+1): + S, L = [], [] + for t in range(1, ell + 1): yt = self.model.getSolVal(sol, y[t]) xt = self.model.getSolVal(sol, x[t]) - if D[t,ell]*yt < xt: + if D[t, ell] * yt < xt: S.append(t) - lhs += D[t,ell]*yt + lhs += D[t, ell] * yt else: L.append(t) lhs += xt - if lhs < D[1,ell]: + if lhs < D[1, ell]: if checkonly: return True else: # add cutting plane constraint self.model.addCons(quicksum([x[t] for t in L]) + \ quicksum(D[t, ell] * y[t] for t in S) - >= D[1, ell], removable = True) + >= D[1, ell], removable=True) cutsadded = True return cutsadded def conscheck(self, constraints, solution, checkintegrality, checklprows, printreason, completely): - if not self.addcut(checkonly = True, sol = solution): + if not self.addcut(checkonly=True, sol=solution): return {"result": SCIP_RESULT.INFEASIBLE} else: return {"result": SCIP_RESULT.FEASIBLE} def consenfolp(self, constraints, nusefulconss, solinfeasible): - if self.addcut(checkonly = False): + if self.addcut(checkonly=False): return {"result": SCIP_RESULT.CONSADDED} else: return {"result": SCIP_RESULT.FEASIBLE} @@ -54,7 +55,8 @@ def consenfolp(self, constraints, nusefulconss, solinfeasible): def conslock(self, constraint, locktype, nlockspos, nlocksneg): pass -def sils(T,f,c,d,h): + +def sils(T, f, c, d, h): """sils -- LP lotsizing for the single item lot sizing problem Parameters: - T: number of periods @@ -66,27 +68,28 @@ def sils(T,f,c,d,h): Returns a model, ready to be solved. """ model = Model("single item lotsizing") - Ts = range(1,T+1) + Ts = range(1, T + 1) M = sum(d[t] for t in Ts) - y,x,I = {},{},{} + y, x, I = {}, {}, {} for t in Ts: - y[t] = model.addVar(vtype="I", ub=1, name="y(%s)"%t) - x[t] = model.addVar(vtype="C", ub=M, name="x(%s)"%t) - I[t] = model.addVar(vtype="C", name="I(%s)"%t) + y[t] = model.addVar(vtype="I", ub=1, name="y(%s)" % t) + x[t] = model.addVar(vtype="C", ub=M, name="x(%s)" % t) + I[t] = model.addVar(vtype="C", name="I(%s)" % t) I[0] = 0 for t in Ts: - model.addCons(x[t] <= M*y[t], "ConstrUB(%s)"%t) - model.addCons(I[t-1] + x[t] == I[t] + d[t], "FlowCons(%s)"%t) + model.addCons(x[t] <= M * y[t], "ConstrUB(%s)" % t) + model.addCons(I[t - 1] + x[t] == I[t] + d[t], "FlowCons(%s)" % t) - model.setObjective(\ - quicksum(f[t]*y[t] + c[t]*x[t] + h[t]*I[t] for t in Ts),\ + model.setObjective( \ + quicksum(f[t] * y[t] + c[t] * x[t] + h[t] * I[t] for t in Ts), \ "minimize") - model.data = y,x,I + model.data = y, x, I return model -def sils_cut(T,f,c,d,h, conshdlr): + +def sils_cut(T, f, c, d, h, conshdlr): """solve_sils -- solve the lot sizing problem with cutting planes - start with a relaxed model - used lazy constraints to elimitate fractional setup variables with cutting planes @@ -99,64 +102,68 @@ def sils_cut(T,f,c,d,h, conshdlr): - h[t]: holding costs Returns the final model solved, with all necessary cuts added. """ - Ts = range(1,T+1) + Ts = range(1, T + 1) - model = sils(T,f,c,d,h) - y,x,I = model.data + model = sils(T, f, c, d, h) + y, x, I = model.data # relax integer variables for t in Ts: model.chgVarType(y[t], "C") - model.addVar(vtype="B", name="fake") # for making the problem MIP + model.addVar(vtype="B", name="fake") # for making the problem MIP # compute D[i,j] = sum_{t=i}^j d[t] D = {} for t in Ts: s = 0 - for j in range(t,T+1): + for j in range(t, T + 1): s += d[j] - D[t,j] = s + D[t, j] = s - #include the lot sizing constraint handler + # include the lot sizing constraint handler model.includeConshdlr(conshdlr, "SILS", "Constraint handler for single item lot sizing", - sepapriority = 0, enfopriority = -1, chckpriority = -1, sepafreq = -1, propfreq = -1, - eagerfreq = -1, maxprerounds = 0, delaysepa = False, delayprop = False, needscons = False, - presoltiming = SCIP_PRESOLTIMING.FAST, proptiming = SCIP_PROPTIMING.BEFORELP) - conshdlr.data = D,Ts + sepapriority=0, enfopriority=-1, chckpriority=-1, sepafreq=-1, propfreq=-1, + eagerfreq=-1, maxprerounds=0, delaysepa=False, delayprop=False, needscons=False, + presoltiming=SCIP_PRESOLTIMING.FAST, proptiming=SCIP_PROPTIMING.BEFORELP) + conshdlr.data = D, Ts - model.data = y,x,I + model.data = y, x, I return model + def mk_example(): """mk_example: book example for the single item lot sizing""" T = 5 - _,f,c,d,h = multidict({ - 1 : [3,1,5,1], - 2 : [3,1,7,1], - 3 : [3,3,3,1], - 4 : [3,3,6,1], - 5 : [3,3,4,1], - }) + _, f, c, d, h = multidict({ + 1: [3, 1, 5, 1], + 2: [3, 1, 7, 1], + 3: [3, 3, 3, 1], + 4: [3, 3, 6, 1], + 5: [3, 3, 4, 1], + }) + + return T, f, c, d, h - return T,f,c,d,h if __name__ == "__main__": - T,f,c,d,h = mk_example() + T, f, c, d, h = mk_example() - model = sils(T,f,c,d,h) - y,x,I = model.data + model = sils(T, f, c, d, h) + y, x, I = model.data model.optimize() - print("\nOptimal value [standard]:",model.getObjVal()) - print("%8s%8s%8s%8s%8s%8s%12s%12s" % ("t","fix","var","h","dem","y","x","I")) - for t in range(1,T+1): - print("%8d%8d%8d%8d%8d%8.1f%12.1f%12.1f" % (t,f[t],c[t],h[t],d[t],model.getVal(y[t]),model.getVal(x[t]),model.getVal(I[t]))) + print("\nOptimal value [standard]:", model.getObjVal()) + print("%8s%8s%8s%8s%8s%8s%12s%12s" % ("t", "fix", "var", "h", "dem", "y", "x", "I")) + for t in range(1, T + 1): + print("%8d%8d%8d%8d%8d%8.1f%12.1f%12.1f" % ( + t, f[t], c[t], h[t], d[t], model.getVal(y[t]), model.getVal(x[t]), model.getVal(I[t]))) conshdlr = Conshdlr_sils() - model = sils_cut(T,f,c,d,h, conshdlr) + model = sils_cut(T, f, c, d, h, conshdlr) model.setBoolParam("misc/allowstrongdualreds", 0) model.optimize() - y,x,I = model.data - print("\nOptimal value [cutting planes]:",model.getObjVal()) - print("%8s%8s%8s%8s%8s%8s%12s%12s" % ("t","fix","var","h","dem","y","x","I")) - for t in range(1,T+1): - print("%8d%8d%8d%8d%8d%8.1f%12.1f%12.1f" % (t,f[t],c[t],h[t],d[t],model.getVal(y[t]),model.getVal(x[t]),model.getVal(I[t]))) + y, x, I = model.data + print("\nOptimal value [cutting planes]:", model.getObjVal()) + print("%8s%8s%8s%8s%8s%8s%12s%12s" % ("t", "fix", "var", "h", "dem", "y", "x", "I")) + for t in range(1, T + 1): + print("%8d%8d%8d%8d%8d%8.1f%12.1f%12.1f" % ( + t, f[t], c[t], h[t], d[t], model.getVal(y[t]), model.getVal(x[t]), model.getVal(I[t]))) diff --git a/examples/finished/markowitz_soco.py b/examples/finished/markowitz_soco.py index 5aaaf24e9..82adab953 100644 --- a/examples/finished/markowitz_soco.py +++ b/examples/finished/markowitz_soco.py @@ -1,5 +1,5 @@ ##@file markowitz_soco.py -#@brief simple markowitz model for portfolio optimization. +# @brief simple markowitz model for portfolio optimization. """ Approach: use second-order cone optimization. @@ -7,7 +7,8 @@ """ from pyscipopt import Model, quicksum, multidict -def markowitz(I,sigma,r,alpha): + +def markowitz(I, sigma, r, alpha): """markowitz -- simple markowitz model for portfolio optimization. Parameters: - I: set of items @@ -20,41 +21,40 @@ def markowitz(I,sigma,r,alpha): x = {} for i in I: - x[i] = model.addVar(vtype="C", name="x(%s)"%i) # quantity of i to buy + x[i] = model.addVar(vtype="C", name="x(%s)" % i) # quantity of i to buy - model.addCons(quicksum(r[i]*x[i] for i in I) >= alpha) + model.addCons(quicksum(r[i] * x[i] for i in I) >= alpha) model.addCons(quicksum(x[i] for i in I) == 1) # set nonlinear objective: SCIP only allow for linear objectives hence the following - obj = model.addVar(vtype="C", name="objective", lb = None, ub = None) # auxiliary variable to represent objective - model.addCons(quicksum(sigma[i]**2 * x[i] * x[i] for i in I) <= obj) + obj = model.addVar(vtype="C", name="objective", lb=None, ub=None) # auxiliary variable to represent objective + model.addCons(quicksum(sigma[i] ** 2 * x[i] * x[i] for i in I) <= obj) model.setObjective(obj, "minimize") model.data = x return model - if __name__ == "__main__": # portfolio - import math - I,sigma,r = multidict( - {1:[0.07,1.01], - 2:[0.09,1.05], - 3:[0.1,1.08], - 4:[0.2,1.10], - 5:[0.3,1.20]} - ) + + I, sigma, r = multidict( + {1: [0.07, 1.01], + 2: [0.09, 1.05], + 3: [0.1, 1.08], + 4: [0.2, 1.10], + 5: [0.3, 1.20]} + ) alpha = 1.05 - model = markowitz(I,sigma,r,alpha) + model = markowitz(I, sigma, r, alpha) model.optimize() x = model.data EPS = 1.e-6 - print("%5s\t%8s" % ("i","x[i]")) + print("%5s\t%8s" % ("i", "x[i]")) for i in I: - print("%5s\t%8g" % (i,model.getVal(x[i]))) - print("sum:",sum(model.getVal(x[i]) for i in I)) + print("%5s\t%8g" % (i, model.getVal(x[i]))) + print("sum:", sum(model.getVal(x[i]) for i in I)) print print("Optimal value:", model.getObjVal()) diff --git a/examples/finished/mctransp.py b/examples/finished/mctransp.py index 248b7c65d..679f8839a 100644 --- a/examples/finished/mctransp.py +++ b/examples/finished/mctransp.py @@ -1,5 +1,5 @@ ##@file mctransp.py -#@brief a model for the multi-commodity transportation problem +# @brief a model for the multi-commodity transportation problem """ Model for solving the multi-commodity transportation problem: minimize the total transportation cost for satisfying demand at @@ -10,7 +10,8 @@ from pyscipopt import Model, quicksum, multidict -def mctransp(I,J,K,c,d,M): + +def mctransp(I, J, K, c, d, M): """mctransp -- model for solving the Multi-commodity Transportation Problem Parameters: - I: set of customers @@ -27,20 +28,20 @@ def mctransp(I,J,K,c,d,M): # Create variables x = {} - for (i,j,k) in c: - x[i,j,k] = model.addVar(vtype="C", name="x(%s,%s,%s)" % (i,j,k)) + for (i, j, k) in c: + x[i, j, k] = model.addVar(vtype="C", name="x(%s,%s,%s)" % (i, j, k)) # Demand constraints for i in I: for k in K: - model.addCons(sum(x[i,j,k] for j in J if (i,j,k) in x) == d[i,k], "Demand(%s,%s)" % (i,k)) + model.addCons(sum(x[i, j, k] for j in J if (i, j, k) in x) == d[i, k], "Demand(%s,%s)" % (i, k)) # Capacity constraints for j in J: - model.addCons(sum(x[i,j,k] for (i,j2,k) in x if j2 == j) <= M[j], "Capacity(%s)" % j) + model.addCons(sum(x[i, j, k] for (i, j2, k) in x if j2 == j) <= M[j], "Capacity(%s)" % j) # Objective - model.setObjective(quicksum(c[i,j,k]*x[i,j,k] for (i,j,k) in x), "minimize") + model.setObjective(quicksum(c[i, j, k] * x[i, j, k] for (i, j, k) in x), "minimize") model.data = x @@ -49,98 +50,99 @@ def mctransp(I,J,K,c,d,M): def make_inst1(): """creates example data set 1""" - d = {(1,1):80, (1,2):85, (1,3):300, (1,4):6, # {(customer,commodity):demand}} - (2,1):270, (2,2):160, (2,3):400, (2,4):7, - (3,1):250, (3,2):130, (3,3):350, (3,4):4, - (4,1):160, (4,2):60, (4,3):200, (4,4):3, - (5,1):180, (5,2):40, (5,3):150, (5,4):5 + d = {(1, 1): 80, (1, 2): 85, (1, 3): 300, (1, 4): 6, # {(customer,commodity):demand}} + (2, 1): 270, (2, 2): 160, (2, 3): 400, (2, 4): 7, + (3, 1): 250, (3, 2): 130, (3, 3): 350, (3, 4): 4, + (4, 1): 160, (4, 2): 60, (4, 3): 200, (4, 4): 3, + (5, 1): 180, (5, 2): 40, (5, 3): 150, (5, 4): 5 } - I = set([i for (i,k) in d]) - K = set([k for (i,k) in d]) - J,M = multidict({1:3000, 2:3000, 3:3000}) # capacity - - produce = {1:[2,4], 2:[1,2,3], 3:[2,3,4]} # products that can be produced in each facility - weight = {1:5, 2:2, 3:3, 4:4} # {commodity: weight} - cost = {(1,1):4, (1,2):6, (1,3):9, # {(customer,factory): cost} - (2,1):5, (2,2):4, (2,3):7, - (3,1):6, (3,2):3, (3,3):4, - (4,1):8, (4,2):5, (4,3):3, - (5,1):10, (5,2):8, (5,3):4 + I = set([i for (i, k) in d]) + K = set([k for (i, k) in d]) + J, M = multidict({1: 3000, 2: 3000, 3: 3000}) # capacity + + produce = {1: [2, 4], 2: [1, 2, 3], 3: [2, 3, 4]} # products that can be produced in each facility + weight = {1: 5, 2: 2, 3: 3, 4: 4} # {commodity: weight} + cost = {(1, 1): 4, (1, 2): 6, (1, 3): 9, # {(customer,factory): cost} + (2, 1): 5, (2, 2): 4, (2, 3): 7, + (3, 1): 6, (3, 2): 3, (3, 3): 4, + (4, 1): 8, (4, 2): 5, (4, 3): 3, + (5, 1): 10, (5, 2): 8, (5, 3): 4 } c = {} for i in I: for j in J: for k in produce[j]: - c[i,j,k] = cost[i,j] * weight[k] + c[i, j, k] = cost[i, j] * weight[k] + + return I, J, K, c, d, M - return I,J,K,c,d,M def make_inst2(): """creates example data set 2""" - d = {(1,1):45, # {(customer,commodity):demand}} - (2,1):20, - (3,1):30, - (4,1):30, + d = {(1, 1): 45, # {(customer,commodity):demand}} + (2, 1): 20, + (3, 1): 30, + (4, 1): 30, } - I = set([i for (i,k) in d]) - K = set([k for (i,k) in d]) - J,M = multidict({1:35, 2:50, 3:40}) # {factory: capacity}} - produce = {1:[1], 2:[1], 3:[1]} # products that can be produced in each facility - weight = {1:1} # {commodity: weight} - cost = {(1,1):8, (1,2):9, (1,3):14, # {(customer,factory): cost} - (2,1):6, (2,2):12, (2,3):9 , - (3,1):10, (3,2):13, (3,3):16, - (4,1):9, (4,2):7, (4,3):5 , + I = set([i for (i, k) in d]) + K = set([k for (i, k) in d]) + J, M = multidict({1: 35, 2: 50, 3: 40}) # {factory: capacity}} + produce = {1: [1], 2: [1], 3: [1]} # products that can be produced in each facility + weight = {1: 1} # {commodity: weight} + cost = {(1, 1): 8, (1, 2): 9, (1, 3): 14, # {(customer,factory): cost} + (2, 1): 6, (2, 2): 12, (2, 3): 9, + (3, 1): 10, (3, 2): 13, (3, 3): 16, + (4, 1): 9, (4, 2): 7, (4, 3): 5, } c = {} for i in I: for j in J: for k in produce[j]: - c[i,j,k] = cost[i,j] * weight[k] + c[i, j, k] = cost[i, j] * weight[k] - return I,J,K,c,d,M + return I, J, K, c, d, M def make_inst3(): """creates example data set 3""" - d = {(1,1):40, (1,2):30, (1,3):10, # {(customer,commodity):demand}} - (2,1):70, (2,2):100, (2,3):100, - (3,1):0, (3,2):0, (3,3):250, - (4,1):60, (4,2):100, (4,3):0, - (5,1):180, (5,2):0, (5,3):0 + d = {(1, 1): 40, (1, 2): 30, (1, 3): 10, # {(customer,commodity):demand}} + (2, 1): 70, (2, 2): 100, (2, 3): 100, + (3, 1): 0, (3, 2): 0, (3, 3): 250, + (4, 1): 60, (4, 2): 100, (4, 3): 0, + (5, 1): 180, (5, 2): 0, (5, 3): 0 } - I = set([i for (i,k) in d]) - K = set([k for (i,k) in d]) - J,M = multidict({1:500, 2:500, 3:500}) # capacity - - produce = {1:[2,4], 2:[1,2,3], 3:[2,3,4]} # products that can be produced in each facility - weight = {1:5, 2:2, 3:3, 4:4} # {commodity: weight} - cost = {(1,1):4, (1,2):6, (1,3):9, # {(customer,factory): cost} - (2,1):5, (2,2):4, (2,3):7, - (3,1):6, (3,2):3, (3,3):4, - (4,1):8, (4,2):5, (4,3):3, - (5,1):10, (5,2):8, (5,3):4 + I = set([i for (i, k) in d]) + K = set([k for (i, k) in d]) + J, M = multidict({1: 500, 2: 500, 3: 500}) # capacity + + produce = {1: [2, 4], 2: [1, 2, 3], 3: [2, 3, 4]} # products that can be produced in each facility + weight = {1: 5, 2: 2, 3: 3, 4: 4} # {commodity: weight} + cost = {(1, 1): 4, (1, 2): 6, (1, 3): 9, # {(customer,factory): cost} + (2, 1): 5, (2, 2): 4, (2, 3): 7, + (3, 1): 6, (3, 2): 3, (3, 3): 4, + (4, 1): 8, (4, 2): 5, (4, 3): 3, + (5, 1): 10, (5, 2): 8, (5, 3): 4 } c = {} for i in I: for j in J: for k in produce[j]: - c[i,j,k] = cost[i,j] * weight[k] + c[i, j, k] = cost[i, j] * weight[k] - return I,J,K,c,d,M + return I, J, K, c, d, M if __name__ == "__main__": - I,J,K,c,d,M = make_inst3(); - model = mctransp(I,J,K,c,d,M) + I, J, K, c, d, M = make_inst3(); + model = mctransp(I, J, K, c, d, M) model.writeProblem("transp.lp") model.optimize() - print("Optimal value:",model.getObjVal()) + print("Optimal value:", model.getObjVal()) EPS = 1.e-6 x = model.data - for (i,j,k) in x: - if model.getVal(x[i,j,k]) > EPS: - print("sending %10s units of %3s from plant %3s to customer %3s" % (model.getVal(x[i,j,k]),k,j,i)) + for (i, j, k) in x: + if model.getVal(x[i, j, k]) > EPS: + print("sending %10s units of %3s from plant %3s to customer %3s" % (model.getVal(x[i, j, k]), k, j, i)) diff --git a/examples/finished/mkp.py b/examples/finished/mkp.py index 8015714a7..cdb2eef73 100644 --- a/examples/finished/mkp.py +++ b/examples/finished/mkp.py @@ -1,11 +1,12 @@ ##@file mkp.py -#@brief model for the multi-constrained knapsack problem +# @brief model for the multi-constrained knapsack problem """ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ from pyscipopt import Model, quicksum, multidict -def mkp(I,J,v,a,b): + +def mkp(I, J, v, a, b): """mkp -- model for solving the multi-constrained knapsack Parameters: - I: set of dimensions @@ -20,14 +21,14 @@ def mkp(I,J,v,a,b): # Create Variables x = {} for j in J: - x[j] = model.addVar(vtype="B", name="x(%s)"%j) + x[j] = model.addVar(vtype="B", name="x(%s)" % j) # Create constraints for i in I: - model.addCons(quicksum(a[i,j]*x[j] for j in J) <= b[i], "Capacity(%s)"%i) + model.addCons(quicksum(a[i, j] * x[j] for j in J) <= b[i], "Capacity(%s)" % i) # Objective - model.setObjective(quicksum(v[j]*x[j] for j in J), "maximize") + model.setObjective(quicksum(v[j] * x[j] for j in J), "maximize") model.data = x return model @@ -35,17 +36,17 @@ def mkp(I,J,v,a,b): def example(): """creates example data set""" - J,v = multidict({1:16, 2:19, 3:23, 4:28}) - a = {(1,1):2, (1,2):3, (1,3):4, (1,4):5, - (2,1):3000, (2,2):3500, (2,3):5100, (2,4):7200, + J, v = multidict({1: 16, 2: 19, 3: 23, 4: 28}) + a = {(1, 1): 2, (1, 2): 3, (1, 3): 4, (1, 4): 5, + (2, 1): 3000, (2, 2): 3500, (2, 3): 5100, (2, 4): 7200, } - I,b = multidict({1:7, 2:10000}) - return I,J,v,a,b + I, b = multidict({1: 7, 2: 10000}) + return I, J, v, a, b if __name__ == "__main__": - I,J,v,a,b = example() - model = mkp(I,J,v,a,b) + I, J, v, a, b = example() + model = mkp(I, J, v, a, b) x = model.data model.optimize() diff --git a/examples/finished/pfs.py b/examples/finished/pfs.py index 523ed1a70..20279ed15 100644 --- a/examples/finished/pfs.py +++ b/examples/finished/pfs.py @@ -1,5 +1,5 @@ ##@file pfs.py -#@brief model for the permutation flow shop problem +# @brief model for the permutation flow shop problem """ Use a position index formulation for modeling the permutation flow shop problem, with the objective of minimizing the makespan (maximum @@ -7,11 +7,12 @@ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ -import math import random -from pyscipopt import Model, quicksum, multidict -def permutation_flow_shop(n,m,p): +from pyscipopt import Model, quicksum + + +def permutation_flow_shop(n, m, p): """gpp -- model for the graph partitioning problem Parameters: - n: number of jobs @@ -21,52 +22,52 @@ def permutation_flow_shop(n,m,p): """ model = Model("permutation flow shop") - x,s,f = {},{},{} - for j in range(1,n+1): - for k in range(1,n+1): - x[j,k] = model.addVar(vtype="B", name="x(%s,%s)"%(j,k)) + x, s, f = {}, {}, {} + for j in range(1, n + 1): + for k in range(1, n + 1): + x[j, k] = model.addVar(vtype="B", name="x(%s,%s)" % (j, k)) - for i in range(1,m+1): - for k in range(1,n+1): - s[i,k] = model.addVar(vtype="C", name="start(%s,%s)"%(i,k)) - f[i,k] = model.addVar(vtype="C", name="finish(%s,%s)"%(i,k)) + for i in range(1, m + 1): + for k in range(1, n + 1): + s[i, k] = model.addVar(vtype="C", name="start(%s,%s)" % (i, k)) + f[i, k] = model.addVar(vtype="C", name="finish(%s,%s)" % (i, k)) - for j in range(1,n+1): - model.addCons(quicksum(x[j,k] for k in range(1,n+1)) == 1, "Assign1(%s)"%(j)) - model.addCons(quicksum(x[k,j] for k in range(1,n+1)) == 1, "Assign2(%s)"%(j)) + for j in range(1, n + 1): + model.addCons(quicksum(x[j, k] for k in range(1, n + 1)) == 1, "Assign1(%s)" % (j)) + model.addCons(quicksum(x[k, j] for k in range(1, n + 1)) == 1, "Assign2(%s)" % (j)) - for i in range(1,m+1): - for k in range(1,n+1): + for i in range(1, m + 1): + for k in range(1, n + 1): if k != n: - model.addCons(f[i,k] <= s[i,k+1], "FinishStart(%s,%s)"%(i,k)) + model.addCons(f[i, k] <= s[i, k + 1], "FinishStart(%s,%s)" % (i, k)) if i != m: - model.addCons(f[i,k] <= s[i+1,k], "Machine(%s,%s)"%(i,k)) + model.addCons(f[i, k] <= s[i + 1, k], "Machine(%s,%s)" % (i, k)) - model.addCons(s[i,k] + quicksum(p[i,j]*x[j,k] for j in range(1,n+1)) <= f[i,k], - "StartFinish(%s,%s)"%(i,k)) + model.addCons(s[i, k] + quicksum(p[i, j] * x[j, k] for j in range(1, n + 1)) <= f[i, k], + "StartFinish(%s,%s)" % (i, k)) - model.setObjective(f[m,n], "minimize") + model.setObjective(f[m, n], "minimize") - model.data = x,s,f + model.data = x, s, f return model -def make_data(n,m): +def make_data(n, m): """make_data: prepare matrix of m times n random processing times""" p = {} - for i in range(1,m+1): - for j in range(1,n+1): - p[i,j] = random.randint(1,10) + for i in range(1, m + 1): + for j in range(1, n + 1): + p[i, j] = random.randint(1, 10) return p def example(): """creates example data set""" - proc = [[2,3,1],[4,2,3],[1,4,1]] + proc = [[2, 3, 1], [4, 2, 3], [1, 4, 1]] p = {} for i in range(3): for j in range(3): - p[i+1,j+1] = proc[j][i] + p[i + 1, j + 1] = proc[j][i] return p @@ -74,21 +75,21 @@ def example(): random.seed(1) n = 15 m = 10 - p = make_data(n,m) + p = make_data(n, m) # n = 3 # m = 3 # p = example() - print("processing times (%s jobs, %s machines):" % (n,m)) - for i in range(1,m+1): - for j in range(1,n+1): - print(p[i,j],) + print("processing times (%s jobs, %s machines):" % (n, m)) + for i in range(1, m + 1): + for j in range(1, n + 1): + print(p[i, j], ) print - model = permutation_flow_shop(n,m,p) + model = permutation_flow_shop(n, m, p) # model.write("permflow.lp") model.optimize() - x,s,f = model.data + x, s, f = model.data print("Optimal value:", model.getObjVal()) ### for (j,k) in sorted(x): @@ -102,5 +103,5 @@ def example(): ### print(f[i].VarName,f[i].X # x[j,k] = 1 if j is the k-th job; extract job sequence: - seq = [j for (k,j) in sorted([(k,j) for (j,k) in x if model.getVal(x[j,k]) > 0.5])] - print("optimal job permutation:",seq) + seq = [j for (k, j) in sorted([(k, j) for (j, k) in x if model.getVal(x[j, k]) > 0.5])] + print("optimal job permutation:", seq) diff --git a/examples/finished/piecewise.py b/examples/finished/piecewise.py index 7e4006243..f1aeb57ce 100644 --- a/examples/finished/piecewise.py +++ b/examples/finished/piecewise.py @@ -1,5 +1,5 @@ ##@file piecewise.py -#@brief several approaches for solving problems with piecewise linear functions. +# @brief several approaches for solving problems with piecewise linear functions. """ Approaches: - mult_selection: multiple selection model @@ -12,10 +12,11 @@ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ import math -import random -from pyscipopt import Model, quicksum, multidict -def mult_selection(model,a,b): +from pyscipopt import Model, quicksum + + +def mult_selection(model, a, b): """mult_selection -- add piecewise relation with multiple selection formulation Parameters: - model: a model where to include the piecewise linear relation @@ -24,29 +25,29 @@ def mult_selection(model,a,b): Returns the model with the piecewise linear relation on added variables X, Y, and z. """ - K = len(a)-1 - w,z = {},{} + K = len(a) - 1 + w, z = {}, {} for k in range(K): - w[k] = model.addVar(lb=-model.infinity()) # do not name variables for avoiding clash + w[k] = model.addVar(lb=-model.infinity()) # do not name variables for avoiding clash z[k] = model.addVar(vtype="B") X = model.addVar(lb=a[0], ub=a[K], vtype="C") Y = model.addVar(lb=-model.infinity()) for k in range(K): - model.addCons(w[k] >= a[k]*z[k]) - model.addCons(w[k] <= a[k+1]*z[k]) + model.addCons(w[k] >= a[k] * z[k]) + model.addCons(w[k] <= a[k + 1] * z[k]) model.addCons(quicksum(z[k] for k in range(K)) == 1) model.addCons(X == quicksum(w[k] for k in range(K))) - c = [float(b[k+1]-b[k])/(a[k+1]-a[k]) for k in range(K)] - d = [b[k]-c[k]*a[k] for k in range(K)] - model.addCons(Y == quicksum(d[k]*z[k] + c[k]*w[k] for k in range(K))) + c = [float(b[k + 1] - b[k]) / (a[k + 1] - a[k]) for k in range(K)] + d = [b[k] - c[k] * a[k] for k in range(K)] + model.addCons(Y == quicksum(d[k] * z[k] + c[k] * w[k] for k in range(K))) - return X,Y,z + return X, Y, z -def convex_comb_sos(model,a,b): +def convex_comb_sos(model, a, b): """convex_comb_sos -- add piecewise relation with gurobi's SOS constraints Parameters: - model: a model where to include the piecewise linear relation @@ -54,23 +55,23 @@ def convex_comb_sos(model,a,b): - b[k]: y-coordinate of the k-th point in the piecewise linear relation Returns the model with the piecewise linear relation on added variables X, Y, and z. """ - K = len(a)-1 + K = len(a) - 1 z = {} - for k in range(K+1): + for k in range(K + 1): z[k] = model.addVar(lb=0, ub=1, vtype="C") X = model.addVar(lb=a[0], ub=a[K], vtype="C") Y = model.addVar(lb=-model.infinity(), vtype="C") - model.addCons(X == quicksum(a[k]*z[k] for k in range(K+1))) - model.addCons(Y == quicksum(b[k]*z[k] for k in range(K+1))) + model.addCons(X == quicksum(a[k] * z[k] for k in range(K + 1))) + model.addCons(Y == quicksum(b[k] * z[k] for k in range(K + 1))) - model.addCons(quicksum(z[k] for k in range(K+1)) == 1) - model.addConsSOS2([z[k] for k in range(K+1)]) + model.addCons(quicksum(z[k] for k in range(K + 1)) == 1) + model.addConsSOS2([z[k] for k in range(K + 1)]) - return X,Y,z + return X, Y, z -def convex_comb_dis(model,a,b): +def convex_comb_dis(model, a, b): """convex_comb_dis -- add piecewise relation with convex combination formulation Parameters: - model: a model where to include the piecewise linear relation @@ -78,8 +79,8 @@ def convex_comb_dis(model,a,b): - b[k]: y-coordinate of the k-th point in the piecewise linear relation Returns the model with the piecewise linear relation on added variables X, Y, and z. """ - K = len(a)-1 - wL,wR,z = {},{},{} + K = len(a) - 1 + wL, wR, z = {}, {}, {} for k in range(K): wL[k] = model.addVar(lb=0, ub=1, vtype="C") wR[k] = model.addVar(lb=0, ub=1, vtype="C") @@ -87,22 +88,22 @@ def convex_comb_dis(model,a,b): X = model.addVar(lb=a[0], ub=a[K], vtype="C") Y = model.addVar(lb=-model.infinity(), vtype="C") - model.addCons(X == quicksum(a[k]*wL[k] + a[k+1]*wR[k] for k in range(K))) - model.addCons(Y == quicksum(b[k]*wL[k] + b[k+1]*wR[k] for k in range(K))) + model.addCons(X == quicksum(a[k] * wL[k] + a[k + 1] * wR[k] for k in range(K))) + model.addCons(Y == quicksum(b[k] * wL[k] + b[k + 1] * wR[k] for k in range(K))) for k in range(K): model.addCons(wL[k] + wR[k] == z[k]) model.addCons(quicksum(z[k] for k in range(K)) == 1) - return X,Y,z + return X, Y, z def gray(i): """returns i^int(i/2)""" - return i^(int(i/2)) + return i ^ (int(i / 2)) -def convex_comb_dis_log(model,a,b): +def convex_comb_dis_log(model, a, b): """convex_comb_dis_log -- add piecewise relation with a logarithmic number of binary variables using the convex combination formulation. Parameters: @@ -111,11 +112,11 @@ def convex_comb_dis_log(model,a,b): - b[k]: y-coordinate of the k-th point in the piecewise linear relation Returns the model with the piecewise linear relation on added variables X, Y, and z. """ - K = len(a)-1 - G = int(math.ceil((math.log(K)/math.log(2)))) # number of required bits - N = 1<zeros:",zeros,"ones:",ones - if (1 & gray(k)>>j) == 1 and (1 & gray(k-1)>>j) == 1: + if (1 & gray(k) >> j) == 1 and (1 & gray(k - 1) >> j) == 1: ones.append(k) - if (1 & gray(k)>>j) == 0 and (1 & gray(k-1)>>j) == 0: + if (1 & gray(k) >> j) == 0 and (1 & gray(k - 1) >> j) == 0: zeros.append(k) # print(j,k,"\tzeros>:",zeros,"ones:",ones # print(j,"\tzeros:",zeros,"ones:",ones model.addCons(quicksum(w[k] for k in ones) <= g[j]) - model.addCons(quicksum(w[k] for k in zeros) <= 1-g[j]) + model.addCons(quicksum(w[k] for k in zeros) <= 1 - g[j]) - return X,Y,w + return X, Y, w if __name__ == "__main__": # random.seed(1) - a = [ -10, 10, 15, 25, 30, 35, 40, 45, 50, 55, 60, 70] - b = [ -20,-20, 15, -21, 0, 50, 18, 0, 15, 24, 10, 15] + a = [-10, 10, 15, 25, 30, 35, 40, 45, 50, 55, 60, 70] + b = [-20, -20, 15, -21, 0, 50, 18, 0, 15, 24, 10, 15] print("\n\n\npiecewise: multiple selection") model = Model("multiple selection") - X,Y,z = mult_selection(model,a,b) # X,Y --> piecewise linear replacement of x,f(x) based on points a,b + X, Y, z = mult_selection(model, a, b) # X,Y --> piecewise linear replacement of x,f(x) based on points a,b # model using X and Y (and possibly other variables) u = model.addVar(vtype="C", name="u") - A = model.addCons(3*X + 4*Y <= 250, "A") - B = model.addCons(7*X - 2*Y + 3*u == 170, "B") - model.setObjective(2*X + 15*Y + 5*u, "maximize") + A = model.addCons(3 * X + 4 * Y <= 250, "A") + B = model.addCons(7 * X - 2 * Y + 3 * u == 170, "B") + model.setObjective(2 * X + 15 * Y + 5 * u, "maximize") model.optimize() - print("X:",model.getVal(X)) - print("Y:",model.getVal(Y)) - print("u:",model.getVal(u)) + print("X:", model.getVal(X)) + print("Y:", model.getVal(Y)) + print("u:", model.getVal(u)) print("\n\n\npiecewise: disaggregated convex combination") model = Model("disaggregated convex combination") - X,Y,z = convex_comb_dis(model,a,b) + X, Y, z = convex_comb_dis(model, a, b) u = model.addVar(vtype="C", name="u") - A = model.addCons(3*X + 4*Y <= 250, "A") - B = model.addCons(7*X - 2*Y + 3*u == 170, "B") - model.setObjective(2*X + 15*Y + 5*u, "maximize") + A = model.addCons(3 * X + 4 * Y <= 250, "A") + B = model.addCons(7 * X - 2 * Y + 3 * u == 170, "B") + model.setObjective(2 * X + 15 * Y + 5 * u, "maximize") model.optimize() - print("X:",model.getVal(X)) - print("Y:",model.getVal(Y)) - print("u:",model.getVal(u)) + print("X:", model.getVal(X)) + print("Y:", model.getVal(Y)) + print("u:", model.getVal(u)) print("\n\n\npiecewise: disaggregated convex combination, logarithmic number of variables") model = Model("disaggregated convex combination (log)") - X,Y,z = convex_comb_dis(model,a,b) + X, Y, z = convex_comb_dis(model, a, b) u = model.addVar(vtype="C", name="u") - A = model.addCons(3*X + 4*Y <= 250, "A") - B = model.addCons(7*X - 2*Y + 3*u == 170, "B") - model.setObjective(2*X + 15*Y + 5*u, "maximize") + A = model.addCons(3 * X + 4 * Y <= 250, "A") + B = model.addCons(7 * X - 2 * Y + 3 * u == 170, "B") + model.setObjective(2 * X + 15 * Y + 5 * u, "maximize") model.optimize() - print("X:",model.getVal(X)) - print("Y:",model.getVal(Y)) - print("u:",model.getVal(u)) + print("X:", model.getVal(X)) + print("Y:", model.getVal(Y)) + print("u:", model.getVal(u)) print("\n\n\npiecewise: SOS2 constraint") model = Model("SOS2") - X,Y,w = convex_comb_sos(model,a,b) + X, Y, w = convex_comb_sos(model, a, b) u = model.addVar(vtype="C", name="u") - A = model.addCons(3*X + 4*Y <= 250, "A") - B = model.addCons(7*X - 2*Y + 3*u == 170, "B") - model.setObjective(2*X + 15*Y + 5*u, "maximize") + A = model.addCons(3 * X + 4 * Y <= 250, "A") + B = model.addCons(7 * X - 2 * Y + 3 * u == 170, "B") + model.setObjective(2 * X + 15 * Y + 5 * u, "maximize") model.optimize() - print("X:",model.getVal(X)) - print("Y:",model.getVal(Y)) - print("u:",model.getVal(u)) + print("X:", model.getVal(X)) + print("Y:", model.getVal(Y)) + print("u:", model.getVal(u)) print("\n\n\npiecewise: aggregated convex combination") model = Model("aggregated convex combination") - X,Y,z = convex_comb_agg(model,a,b) + X, Y, z = convex_comb_agg(model, a, b) u = model.addVar(vtype="C", name="u") - A = model.addCons(3*X + 4*Y <= 250, "A") - B = model.addCons(7*X - 2*Y + 3*u == 170, "B") - model.setObjective(2*X + 15*Y + 5*u, "maximize") + A = model.addCons(3 * X + 4 * Y <= 250, "A") + B = model.addCons(7 * X - 2 * Y + 3 * u == 170, "B") + model.setObjective(2 * X + 15 * Y + 5 * u, "maximize") model.optimize() - print("X:",model.getVal(X)) - print("Y:",model.getVal(Y)) - print("u:",model.getVal(u)) + print("X:", model.getVal(X)) + print("Y:", model.getVal(Y)) + print("u:", model.getVal(u)) print("\n\n\npiecewise: aggregated convex combination, logarithmic number of variables") model = Model("aggregated convex combination (log)") - X,Y,w = convex_comb_agg_log(model,a,b) + X, Y, w = convex_comb_agg_log(model, a, b) u = model.addVar(vtype="C", name="u") - A = model.addCons(3*X + 4*Y <= 250, "A") - B = model.addCons(7*X - 2*Y + 3*u == 170, "B") - model.setObjective(2*X + 15*Y + 5*u, "maximize") + A = model.addCons(3 * X + 4 * Y <= 250, "A") + B = model.addCons(7 * X - 2 * Y + 3 * u == 170, "B") + model.setObjective(2 * X + 15 * Y + 5 * u, "maximize") model.optimize() - print("X:",model.getVal(X)) - print("Y:",model.getVal(Y)) - print("u:",model.getVal(u)) + print("X:", model.getVal(X)) + print("Y:", model.getVal(Y)) + print("u:", model.getVal(u)) diff --git a/examples/finished/prodmix_soco.py b/examples/finished/prodmix_soco.py index a3eec4511..3fbb258dc 100644 --- a/examples/finished/prodmix_soco.py +++ b/examples/finished/prodmix_soco.py @@ -1,11 +1,12 @@ ##@file prodmix_soco.py -#@brief product mix model using soco. +# @brief product mix model using soco. """ Copyright (c) by Joao Pedro PEDROSO, Masahiro MURAMATSU and Mikio KUBO, 2012 """ from pyscipopt import Model, quicksum, multidict -def prodmix(I,K,a,p,epsilon,LB): + +def prodmix(I, K, a, p, epsilon, LB): """prodmix: robust production planning using soco Parameters: I - set of materials @@ -18,41 +19,41 @@ def prodmix(I,K,a,p,epsilon,LB): model = Model("robust product mix") - x,rhs = {},{} + x, rhs = {}, {} for i in I: - x[i] = model.addVar(vtype="C", name="x(%s)"%i) + x[i] = model.addVar(vtype="C", name="x(%s)" % i) for k in K: - rhs[k] = model.addVar(vtype="C", name="rhs(%s)"%k) + rhs[k] = model.addVar(vtype="C", name="rhs(%s)" % k) model.addCons(quicksum(x[i] for i in I) == 1) for k in K: - model.addCons(rhs[k] == -LB[k]+ quicksum(a[i,k]*x[i] for i in I) ) - model.addCons(quicksum(epsilon*epsilon*x[i]*x[i] for i in I) <= rhs[k]*rhs[k]) + model.addCons(rhs[k] == -LB[k] + quicksum(a[i, k] * x[i] for i in I)) + model.addCons(quicksum(epsilon * epsilon * x[i] * x[i] for i in I) <= rhs[k] * rhs[k]) - model.setObjective(quicksum(p[i]*x[i] for i in I), "minimize") + model.setObjective(quicksum(p[i] * x[i] for i in I), "minimize") - model.data = x,rhs + model.data = x, rhs return model def make_data(): """creates example data set""" - a = { (1,1):.25, (1,2):.15, (1,3):.2, - (2,1):.3, (2,2):.3, (2,3):.1, - (3,1):.15, (3,2):.65, (3,3):.05, - (4,1):.1, (4,2):.05, (4,3):.8 - } + a = {(1, 1): .25, (1, 2): .15, (1, 3): .2, + (2, 1): .3, (2, 2): .3, (2, 3): .1, + (3, 1): .15, (3, 2): .65, (3, 3): .05, + (4, 1): .1, (4, 2): .05, (4, 3): .8 + } epsilon = 0.01 - I,p = multidict({1:5, 2:6, 3:8, 4:20}) - K,LB = multidict({1:.2, 2:.3, 3:.2}) - return I,K,a,p,epsilon,LB + I, p = multidict({1: 5, 2: 6, 3: 8, 4: 20}) + K, LB = multidict({1: .2, 2: .3, 3: .2}) + return I, K, a, p, epsilon, LB if __name__ == "__main__": - I,K,a,p,epsilon,LB = make_data() - model = prodmix(I,K,a,p,epsilon,LB) + I, K, a, p, epsilon, LB = make_data() + model = prodmix(I, K, a, p, epsilon, LB) model.optimize() - print("Objective value:",model.getObjVal()) - x,rhs = model.data + print("Objective value:", model.getObjVal()) + x, rhs = model.data for i in I: - print(i,": ",model.getVal(x[i])) + print(i, ": ", model.getVal(x[i])) diff --git a/examples/finished/rcs.py b/examples/finished/rcs.py index c08b7ff4d..aad89a184 100644 --- a/examples/finished/rcs.py +++ b/examples/finished/rcs.py @@ -1,11 +1,12 @@ ##@file rcs.py -#@brief model for the resource constrained scheduling problem +# @brief model for the resource constrained scheduling problem """ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ from pyscipopt import Model, quicksum, multidict -def rcs(J,P,R,T,p,c,a,RUB): + +def rcs(J, P, R, T, p, c, a, RUB): """rcs -- model for the resource constrained scheduling problem Parameters: - J: set of jobs @@ -20,106 +21,108 @@ def rcs(J,P,R,T,p,c,a,RUB): """ model = Model("resource constrained scheduling") - s,x = {},{} # s - start time variable; x=1 if job j starts on period t + s, x = {}, {} # s - start time variable; x=1 if job j starts on period t for j in J: - s[j] = model.addVar(vtype="C", name="s(%s)"%j) - for t in range(1,T-p[j]+2): - x[j,t] = model.addVar(vtype="B", name="x(%s,%s)"%(j,t)) + s[j] = model.addVar(vtype="C", name="s(%s)" % j) + for t in range(1, T - p[j] + 2): + x[j, t] = model.addVar(vtype="B", name="x(%s,%s)" % (j, t)) for j in J: # job execution constraints - model.addCons(quicksum(x[j,t] for t in range(1,T-p[j]+2)) == 1, "ConstrJob(%s,%s)"%(j,t)) + model.addCons(quicksum(x[j, t] for t in range(1, T - p[j] + 2)) == 1, "ConstrJob(%s,%s)" % (j, t)) # start time constraints - model.addCons(quicksum((t-1)*x[j,t] for t in range(2,T-p[j]+2)) == s[j], "ConstrJob(%s,%s)"%(j,t)) + model.addCons(quicksum((t - 1) * x[j, t] for t in range(2, T - p[j] + 2)) == s[j], "ConstrJob(%s,%s)" % (j, t)) # resource upper bound constraints - for t in range(1,T-p[j]+2): + for t in range(1, T - p[j] + 2): for r in R: model.addCons( - quicksum(a[j,r,t-t_]*x[j,t_] for j in J for t_ in range(max(t-p[j]+1,1),min(t+1,T-p[j]+2))) \ - <= RUB[r,t], "ResourceUB(%s)"%t) + quicksum(a[j, r, t - t_] * x[j, t_] for j in J for t_ in + range(max(t - p[j] + 1, 1), min(t + 1, T - p[j] + 2))) \ + <= RUB[r, t], "ResourceUB(%s)" % t) # time (precedence) constraints, i.e., s[k]-s[j] >= p[j] - for (j,k) in P: - model.addCons(s[k] - s[j] >= p[j], "Precedence(%s,%s)"%(j,k)) + for (j, k) in P: + model.addCons(s[k] - s[j] >= p[j], "Precedence(%s,%s)" % (j, k)) - model.setObjective(quicksum(c[j,t]*x[j,t] for (j,t) in x), "minimize") + model.setObjective(quicksum(c[j, t] * x[j, t] for (j, t) in x), "minimize") - model.data = x,s + model.data = x, s return model def make_1r(): """creates example data set 1""" - J, p = multidict({ # jobs, processing times - 1 : 1, - 2 : 3, - 3 : 2, - 4 : 2, - }) - P = [(1,2), (1,3), (2,4)] + J, p = multidict({ # jobs, processing times + 1: 1, + 2: 3, + 3: 2, + 4: 2, + }) + P = [(1, 2), (1, 3), (2, 4)] R = [1] T = 6 c = {} for j in J: - for t in range(1,T-p[j]+2): - c[j,t] = 1*(t-1+p[j]) + for t in range(1, T - p[j] + 2): + c[j, t] = 1 * (t - 1 + p[j]) a = { - (1,1,0):2, - (2,1,0):2, (2,1,1):1, (2,1,2):1, - (3,1,0):1, (3,1,1):1, - (4,1,0):1, (4,1,1):2, - } - RUB = {(1,1):2, (1,2):2, (1,3):1, (1,4):2, (1,5):2, (1,6):2} - return (J,P,R,T,p,c,a,RUB) + (1, 1, 0): 2, + (2, 1, 0): 2, (2, 1, 1): 1, (2, 1, 2): 1, + (3, 1, 0): 1, (3, 1, 1): 1, + (4, 1, 0): 1, (4, 1, 1): 2, + } + RUB = {(1, 1): 2, (1, 2): 2, (1, 3): 1, (1, 4): 2, (1, 5): 2, (1, 6): 2} + return (J, P, R, T, p, c, a, RUB) + def make_2r(): """creates example data set 2""" - J, p = multidict({ # jobs, processing times - 1 : 2, - 2 : 2, - 3 : 3, - 4 : 2, - 5 : 5, - }) - P = [(1,2), (1,3), (2,4)] - R = [1,2] + J, p = multidict({ # jobs, processing times + 1: 2, + 2: 2, + 3: 3, + 4: 2, + 5: 5, + }) + P = [(1, 2), (1, 3), (2, 4)] + R = [1, 2] T = 6 c = {} for j in J: - for t in range(1,T-p[j]+2): - c[j,t] = 1*(t-1+p[j]) + for t in range(1, T - p[j] + 2): + c[j, t] = 1 * (t - 1 + p[j]) a = { # resource 1: - (1,1,0):2, (1,1,1):2, - (2,1,0):1, (2,1,1):1, - (3,1,0):1, (3,1,1):1, (3,1,2):1, - (4,1,0):1, (4,1,1):1, - (5,1,0):0, (5,1,1):0, (5,1,2):1, (5,1,3):0, (5,1,4):0, + (1, 1, 0): 2, (1, 1, 1): 2, + (2, 1, 0): 1, (2, 1, 1): 1, + (3, 1, 0): 1, (3, 1, 1): 1, (3, 1, 2): 1, + (4, 1, 0): 1, (4, 1, 1): 1, + (5, 1, 0): 0, (5, 1, 1): 0, (5, 1, 2): 1, (5, 1, 3): 0, (5, 1, 4): 0, # resource 2: - (1,2,0):1, (1,2,1):0, - (2,2,0):1, (2,2,1):1, - (3,2,0):0, (3,2,1):0, (3,2,2):0, - (4,2,0):1, (4,2,1):2, - (5,2,0):1, (5,2,1):2, (5,2,2):1, (5,2,3):1, (5,2,4):1, - } - RUB = {(1,1):2, (1,2):2, (1,3):2, (1,4):2, (1,5):2, (1,6):2, (1,7):2, - (2,1):2, (2,2):2, (2,3):2, (2,4):2, (2,5):2, (2,6):2, (2,7):2 } - return (J,P,R,T,p,c,a,RUB) + (1, 2, 0): 1, (1, 2, 1): 0, + (2, 2, 0): 1, (2, 2, 1): 1, + (3, 2, 0): 0, (3, 2, 1): 0, (3, 2, 2): 0, + (4, 2, 0): 1, (4, 2, 1): 2, + (5, 2, 0): 1, (5, 2, 1): 2, (5, 2, 2): 1, (5, 2, 3): 1, (5, 2, 4): 1, + } + RUB = {(1, 1): 2, (1, 2): 2, (1, 3): 2, (1, 4): 2, (1, 5): 2, (1, 6): 2, (1, 7): 2, + (2, 1): 2, (2, 2): 2, (2, 3): 2, (2, 4): 2, (2, 5): 2, (2, 6): 2, (2, 7): 2} + return (J, P, R, T, p, c, a, RUB) if __name__ == "__main__": - (J,P,R,T,p,c,a,RUB) = make_2r() - model = rcs(J,P,R,T,p,c,a,RUB) + (J, P, R, T, p, c, a, RUB) = make_2r() + model = rcs(J, P, R, T, p, c, a, RUB) model.optimize() - x,s = model.data + x, s = model.data - print("Optimal value:",model.getObjVal()) - for (j,t) in x: - if model.getVal(x[j,t]) > 0.5: - print(x[j,t].name,"=",model.getVal(x[j,t])) + print("Optimal value:", model.getObjVal()) + for (j, t) in x: + if model.getVal(x[j, t]) > 0.5: + print(x[j, t].name, "=", model.getVal(x[j, t])) for j in s: if model.getVal(s[j]) > 0.: - print(s[j].name,"=",model.getVal(s[j])) + print(s[j].name, "=", model.getVal(s[j])) diff --git a/examples/finished/read_tsplib.py b/examples/finished/read_tsplib.py index a2ec84307..770de663a 100644 --- a/examples/finished/read_tsplib.py +++ b/examples/finished/read_tsplib.py @@ -1,5 +1,5 @@ ##@file read_tsplib.py -#@brief read standard instances of the traveling salesman problem +# @brief read standard instances of the traveling salesman problem """ Functions provided: * read_tsplib - read a symmetric tsp instance @@ -11,7 +11,8 @@ import gzip import math -def distL2(x1,y1,x2,y2): + +def distL2(x1, y1, x2, y2): """Compute the L2-norm (Euclidean) distance between two points. The distance is rounded to the closest integer, for compatibility @@ -21,10 +22,10 @@ def distL2(x1,y1,x2,y2): sent as parameters""" xdiff = x2 - x1 ydiff = y2 - y1 - return int(math.sqrt(xdiff*xdiff + ydiff*ydiff) + .5) + return int(math.sqrt(xdiff * xdiff + ydiff * ydiff) + .5) -def distL1(x1,y1,x2,y2): +def distL1(x1, y1, x2, y2): """Compute the L1-norm (Manhattan) distance between two points. The distance is rounded to the closest integer, for compatibility @@ -32,110 +33,117 @@ def distL1(x1,y1,x2,y2): The two points are located on coordinates (x1,y1) and (x2,y2), sent as parameters""" - return int(abs(x2-x1) + abs(y2-y1)+.5) + return int(abs(x2 - x1) + abs(y2 - y1) + .5) -def distLinf(x1,y1,x2,y2): +def distLinf(x1, y1, x2, y2): """Compute the Linfty distance between two points (see TSPLIB documentation)""" - return int(max(abs(x2-x1),abs(y2-y1))) + return int(max(abs(x2 - x1), abs(y2 - y1))) + -def distATT(x1,y1,x2,y2): +def distATT(x1, y1, x2, y2): """Compute the ATT distance between two points (see TSPLIB documentation)""" xd = x2 - x1 yd = y2 - y1 - rij = math.sqrt((xd*xd + yd*yd) /10.) + rij = math.sqrt((xd * xd + yd * yd) / 10.) tij = int(rij + .5) if tij < rij: return tij + 1 else: return tij -def distCEIL2D(x1,y1,x2,y2): + +def distCEIL2D(x1, y1, x2, y2): """returns smallest integer not less than the distance of two points""" xdiff = x2 - x1 ydiff = y2 - y1 - return int(math.ceil(math.sqrt(xdiff*xdiff + ydiff*ydiff))) + return int(math.ceil(math.sqrt(xdiff * xdiff + ydiff * ydiff))) + -def distGEO(x1,y1,x2,y2): +def distGEO(x1, y1, x2, y2): print("Implementation is wrong") assert False PI = 3.141592 deg = int(x1 + .5) min_ = x1 - deg - lat1 = PI * (deg + 5.*min_/3)/180. + lat1 = PI * (deg + 5. * min_ / 3) / 180. deg = int(y1 + .5) min_ = y1 - deg - long1 = PI * (deg + 5.*min_/3)/180. + long1 = PI * (deg + 5. * min_ / 3) / 180. deg = int(x2 + .5) min_ = x2 - deg - lat2 = PI * (deg + 5.*min_/3)/180. + lat2 = PI * (deg + 5. * min_ / 3) / 180. deg = int(y2 + .5) min_ = y2 - deg - long2 = PI * (deg + 5.*min_/3)/180. - + long2 = PI * (deg + 5. * min_ / 3) / 180. RRR = 6378.388 - q1 = math.cos( long1 - long2 ); - q2 = math.cos( lat1 - lat2 ); - q3 = math.cos( lat1 + lat2 ); - return int(RRR * math.acos(.5*((1.+q1)*q2 - (1.-q1)*q3)) + 1.) + q1 = math.cos(long1 - long2); + q2 = math.cos(lat1 - lat2); + q3 = math.cos(lat1 + lat2); + return int(RRR * math.acos(.5 * ((1. + q1) * q2 - (1. - q1) * q3)) + 1.) + -def read_explicit_lowerdiag(f,n): +def read_explicit_lowerdiag(f, n): c = {} - i,j = 1,1 + i, j = 1, 1 while True: line = f.readline() for data in line.split(): - c[j,i] = int(data) + c[j, i] = int(data) j += 1 - if j>i: + if j > i: i += 1 j = 1 if i > n: - return range(1,n+1),c,None,None + return range(1, n + 1), c, None, None + -def read_explicit_upper(f,n): +def read_explicit_upper(f, n): c = {} - i,j = 1,2 + i, j = 1, 2 while True: line = f.readline() for data in line.split(): - c[i,j] = int(data) + c[i, j] = int(data) j += 1 - if j>n: + if j > n: i += 1 - j = i+1 + j = i + 1 if i == n: - return range(1,n+1),c,None,None + return range(1, n + 1), c, None, None -def read_explicit_upperdiag(f,n): + +def read_explicit_upperdiag(f, n): c = {} - i,j = 1,1 + i, j = 1, 1 while True: line = f.readline() for data in line.split(): - c[i,j] = int(data) + c[i, j] = int(data) j += 1 - if j>n: + if j > n: i += 1 j = i if i == n: - return range(1,n+1),c,None,None + return range(1, n + 1), c, None, None + -def read_explicit_matrix(f,n): +def read_explicit_matrix(f, n): c = {} - i,j = 1,1 + i, j = 1, 1 while True: line = f.readline() for data in line.split(): - if j>i: - c[i,j] = int(data) + if j > i: + c[i, j] = int(data) j += 1 - if j>n: + if j > n: i += 1 j = 1 if i == n: - return range(1,n+1),c,None,None + return range(1, n + 1), c, None, None + def read_tsplib(filename): "basic function for reading a symmetric problem in the TSPLIB format" @@ -173,21 +181,21 @@ def read_tsplib(filename): if line.find("LOWER_DIAG_ROW") != -1: while line.find("EDGE_WEIGHT_SECTION") == -1: line = f.readline() - return read_explicit_lowerdiag(f,n) + return read_explicit_lowerdiag(f, n) if line.find("UPPER_ROW") != -1: while line.find("EDGE_WEIGHT_SECTION") == -1: line = f.readline() - return read_explicit_upper(f,n) + return read_explicit_upper(f, n) if line.find("UPPER_DIAG_ROW") != -1: while line.find("EDGE_WEIGHT_SECTION") == -1: line = f.readline() - return read_explicit_upperdiag(f,n) + return read_explicit_upperdiag(f, n) if line.find("FULL_MATRIX") != -1: while line.find("EDGE_WEIGHT_SECTION") == -1: line = f.readline() - return read_explicit_matrix(f,n) + return read_explicit_matrix(f, n) print("error reading line " + line) - raise(Exception) + raise (Exception) else: print("cannot deal with '%s' distances" % line) raise Exception @@ -195,22 +203,21 @@ def read_tsplib(filename): while line.find("NODE_COORD_SECTION") == -1: line = f.readline() - x,y = {},{} + x, y = {}, {} while 1: line = f.readline() if line.find("EOF") != -1 or not line: break - (i,xi,yi) = line.split() + (i, xi, yi) = line.split() x[i] = float(xi) y[i] = float(yi) V = x.keys() - c = {} # dictionary to hold n times n matrix + c = {} # dictionary to hold n times n matrix for i in V: for j in V: - c[i,j] = dist(x[i],y[i],x[j],y[j]) - - return V,c,x,y + c[i, j] = dist(x[i], y[i], x[j], y[j]) + return V, c, x, y def read_atsplib(filename): @@ -238,7 +245,7 @@ def read_atsplib(filename): else: raise IOError("'EDGE_WEIGHT_TYPE' is not 'EXPLICIT' in file '%s'" % filename) - for k,line in enumerate(data): + for k, line in enumerate(data): if line.find("EDGE_WEIGHT_SECTION") >= 0: break else: @@ -247,7 +254,7 @@ def read_atsplib(filename): c = {} # flatten list of distances dist = [] - for line in data[k+1:]: + for line in data[k + 1:]: if line.find("EOF") >= 0: break for val in line.split(): @@ -256,31 +263,32 @@ def read_atsplib(filename): k = 0 for i in range(n): for j in range(n): - c[i+1,j+1] = dist[k] + c[i + 1, j + 1] = dist[k] k += 1 - return n,c - + return n, c if __name__ == "__main__": import sys + # Parse argument if len(sys.argv) < 2: print('Usage: %s instance' % sys.argv[0]) exit(1) from read_tsplib import read_tsplib - V,c,x,y = read_tsplib(sys.argv[1]) + + V, c, x, y = read_tsplib(sys.argv[1]) print(len(V), "vertices,", len(c), "arcs") print("distance matrix:") for i in V: for j in V: if j > i: - print(c[i,j],) + print(c[i, j], ) elif j < i: - print(c[j,i],) + print(c[j, i], ) else: - print(0,) + print(0, ) print print diff --git a/examples/finished/ssa.py b/examples/finished/ssa.py index 9e722a758..7da500ce6 100644 --- a/examples/finished/ssa.py +++ b/examples/finished/ssa.py @@ -1,17 +1,19 @@ ##@file ssa.py -#@brief multi-stage (serial) safety stock allocation model +# @brief multi-stage (serial) safety stock allocation model """ Approach: use SOS2 constraints for modeling non-linear functions. Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ -from pyscipopt import Model, quicksum, multidict import math import random +from pyscipopt import Model, quicksum + from piecewise import convex_comb_sos -def ssa(n,h,K,f,T): + +def ssa(n, h, K, f, T): """ssa -- multi-stage (serial) safety stock allocation model Parameters: - n: number of stages @@ -25,72 +27,73 @@ def ssa(n,h,K,f,T): model = Model("safety stock allocation") # calculate endpoints for linear segments - a,b = {},{} - for i in range(1,n+1): + a, b = {}, {} + for i in range(1, n + 1): a[i] = [k for k in range(K)] - b[i] = [f(i,k) for k in range(K)] + b[i] = [f(i, k) for k in range(K)] # x: net replenishment time for stage i # y: corresponding cost # s: piecewise linear segment of variable x - x,y,s = {},{},{} - L = {} # service time of stage i - for i in range(1,n+1): - x[i],y[i],s[i] = convex_comb_sos(model,a[i],b[i]) + x, y, s = {}, {}, {} + L = {} # service time of stage i + for i in range(1, n + 1): + x[i], y[i], s[i] = convex_comb_sos(model, a[i], b[i]) if i == 1: - L[i] = model.addVar(ub=0, vtype="C", name="L[%s]"%i) + L[i] = model.addVar(ub=0, vtype="C", name="L[%s]" % i) else: - L[i] = model.addVar(vtype="C", name="L[%s]"%i) - L[n+1] = model.addVar(ub=0, vtype="C", name="L[%s]"%(n+1)) + L[i] = model.addVar(vtype="C", name="L[%s]" % i) + L[n + 1] = model.addVar(ub=0, vtype="C", name="L[%s]" % (n + 1)) - for i in range(1,n+1): + for i in range(1, n + 1): # net replenishment time for each stage i - model.addCons(x[i] + L[i] == T[i] + L[i+1]) + model.addCons(x[i] + L[i] == T[i] + L[i + 1]) - model.setObjective(quicksum(h[i]*y[i] for i in range(1,n+1)), "minimize") + model.setObjective(quicksum(h[i] * y[i] for i in range(1, n + 1)), "minimize") - model.data = x,s,L + model.data = x, s, L return model - def make_data(): """creates example data set""" - n = 30 # number of stages - z = 1.65 # for 95% service level - sigma = 100 # demand's standard deviation - h = {} # inventory cost - T = {} # production lead time + n = 30 # number of stages + z = 1.65 # for 95% service level + sigma = 100 # demand's standard deviation + h = {} # inventory cost + T = {} # production lead time h[n] = 1 - for i in range(n-1,0,-1): - h[i] = h[i+1] + random.randint(30,50) - K = 0 # number of segments (=sum of processing times) - for i in range(1,n+1): - T[i] = random.randint(3,5) # production lead time at stage i + for i in range(n - 1, 0, -1): + h[i] = h[i + 1] + random.randint(30, 50) + K = 0 # number of segments (=sum of processing times) + for i in range(1, n + 1): + T[i] = random.randint(3, 5) # production lead time at stage i K += T[i] - return z,sigma,h,T,K,n - + return z, sigma, h, T, K, n if __name__ == "__main__": random.seed(1) - z,sigma,h,T,K,n = make_data() - def f(i,k): - return sigma*z*math.sqrt(k) + z, sigma, h, T, K, n = make_data() + + + def f(i, k): + return sigma * z * math.sqrt(k) + - model = ssa(n,h,K,f,T) + model = ssa(n, h, K, f, T) model.optimize() # model.write("ssa.lp") - x,s,L = model.data - for i in range(1,n+1): + x, s, L = model.data + for i in range(1, n + 1): for k in range(K): if model.getVal(s[i][k]) >= 0.001: - print(s[i][k].name,model.getVal(s[i][k])) + print(s[i][k].name, model.getVal(s[i][k])) print - print("%10s%10s%10s%10s" % ("Period","x","L","T")) - for i in range(1,n+1): - print("%10s%10s%10s%10s" % (i,model.getVal(x[i]), model.getVal(L[i]), T[i])) + print("%10s%10s%10s%10s" % ("Period", "x", "L", "T")) + for i in range(1, n + 1): + print("%10s%10s%10s%10s" % (i, model.getVal(x[i]), model.getVal(L[i]), T[i])) - print("Objective:",model.getObjVal()) + print("Objective:", model.getObjVal()) diff --git a/examples/finished/ssp.py b/examples/finished/ssp.py index 0b3a31d66..df4da4dab 100644 --- a/examples/finished/ssp.py +++ b/examples/finished/ssp.py @@ -1,12 +1,13 @@ ##@file ssp.py -#@brief model for the stable set problem +# @brief model for the stable set problem """ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ -from pyscipopt import Model, quicksum, multidict +from pyscipopt import Model, quicksum -def ssp(V,E): + +def ssp(V, E): """ssp -- model for the stable set problem Parameters: - V: set/list of nodes in the graph @@ -17,10 +18,10 @@ def ssp(V,E): x = {} for i in V: - x[i] = model.addVar(vtype="B", name="x(%s)"%i) + x[i] = model.addVar(vtype="B", name="x(%s)" % i) - for (i,j) in E: - model.addCons(x[i] + x[j] <= 1, "Edge(%s,%s)"%(i,j)) + for (i, j) in E: + model.addCons(x[i] + x[j] <= 1, "Edge(%s,%s)" % (i, j)) model.setObjective(quicksum(x[i] for i in V), "maximize") @@ -29,23 +30,25 @@ def ssp(V,E): import random -def make_data(n,prob): + + +def make_data(n, prob): """make_data: prepare data for a random graph Parameters: - n: number of vertices - prob: probability of existence of an edge, for each pair of vertices Returns a tuple with a list of vertices and a list edges. """ - V = range(1,n+1) - E = [(i,j) for i in V for j in V if i < j and random.random() < prob] - return V,E + V = range(1, n + 1) + E = [(i, j) for i in V for j in V if i < j and random.random() < prob] + return V, E if __name__ == "__main__": random.seed(1) - V,E = make_data(100,.5) + V, E = make_data(100, .5) - model = ssp(V,E) + model = ssp(V, E) model.optimize() print("Optimal value:", model.getObjVal()) diff --git a/examples/finished/sudoku.py b/examples/finished/sudoku.py index d0cf352d4..172b21585 100755 --- a/examples/finished/sudoku.py +++ b/examples/finished/sudoku.py @@ -1,7 +1,7 @@ ##@file sudoku.py -#@brief Simple example of modeling a Sudoku as a binary program +# @brief Simple example of modeling a Sudoku as a binary program -#!/usr/bin/env python +# !/usr/bin/env python from pyscipopt import Model, quicksum @@ -23,31 +23,31 @@ for i in range(9): for j in range(9): for k in range(9): - name = str(i)+','+str(j)+','+str(k) - x[i,j,k] = m.addVar(name, vtype='B') + name = str(i) + ',' + str(j) + ',' + str(k) + x[i, j, k] = m.addVar(name, vtype='B') # fill in initial values for i in range(9): for j in range(9): - if init[j + 9*i] != 0: - m.addCons(x[i,j,init[j + 9*i]-1] == 1) + if init[j + 9 * i] != 0: + m.addCons(x[i, j, init[j + 9 * i] - 1] == 1) # only one digit in every field for i in range(9): for j in range(9): - m.addCons(quicksum(x[i,j,k] for k in range(9)) == 1) + m.addCons(quicksum(x[i, j, k] for k in range(9)) == 1) # set up row and column constraints for ind in range(9): for k in range(9): - m.addCons(quicksum(x[ind,j,k] for j in range(9)) == 1) - m.addCons(quicksum(x[i,ind,k] for i in range(9)) == 1) + m.addCons(quicksum(x[ind, j, k] for j in range(9)) == 1) + m.addCons(quicksum(x[i, ind, k] for i in range(9)) == 1) # set up square constraints for row in range(3): for col in range(3): for k in range(9): - m.addCons(quicksum(x[i+3*row, j+3*col, k] for i in range(3) for j in range(3)) == 1) + m.addCons(quicksum(x[i + 3 * row, j + 3 * col, k] for i in range(3) for j in range(3)) == 1) m.hideOutput() m.optimize() @@ -61,7 +61,7 @@ out = '' for j in range(9): for k in range(9): - if m.getVal(x[i,j,k]) == 1: - sol[i,j] = k+1 - out += str(sol[i,j]) + ' ' + if m.getVal(x[i, j, k]) == 1: + sol[i, j] = k + 1 + out += str(sol[i, j]) + ' ' print(out) diff --git a/examples/finished/transp.py b/examples/finished/transp.py index ca36c38bd..858e033aa 100644 --- a/examples/finished/transp.py +++ b/examples/finished/transp.py @@ -1,5 +1,5 @@ ##@file transp.py -#@brief a model for the transportation problem +# @brief a model for the transportation problem """ Model for solving a transportation problem: minimize the total transportation cost for satisfying demand at @@ -9,7 +9,8 @@ """ from pyscipopt import Model, quicksum, multidict -def transp(I,J,c,d,M): + +def transp(I, J, c, d, M): """transp -- model for solving the transportation problem Parameters: I - set of customers @@ -27,18 +28,18 @@ def transp(I,J,c,d,M): for i in I: for j in J: - x[i,j] = model.addVar(vtype="C", name="x(%s,%s)" % (i, j)) + x[i, j] = model.addVar(vtype="C", name="x(%s,%s)" % (i, j)) # Demand constraints for i in I: - model.addCons(quicksum(x[i,j] for j in J if (i,j) in x) == d[i], name="Demand(%s)" % i) + model.addCons(quicksum(x[i, j] for j in J if (i, j) in x) == d[i], name="Demand(%s)" % i) # Capacity constraints for j in J: - model.addCons(quicksum(x[i,j] for i in I if (i,j) in x) <= M[j], name="Capacity(%s)" % j) + model.addCons(quicksum(x[i, j] for i in I if (i, j) in x) <= M[j], name="Capacity(%s)" % j) # Objective - model.setObjective(quicksum(c[i,j]*x[i,j] for (i,j) in x), "minimize") + model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), "minimize") model.optimize() @@ -48,33 +49,33 @@ def transp(I,J,c,d,M): def make_inst1(): """creates example data set 1""" - I,d = multidict({1:80, 2:270, 3:250 , 4:160, 5:180}) # demand - J,M = multidict({1:500, 2:500, 3:500}) # capacity - c = {(1,1):4, (1,2):6, (1,3):9, # cost - (2,1):5, (2,2):4, (2,3):7, - (3,1):6, (3,2):3, (3,3):4, - (4,1):8, (4,2):5, (4,3):3, - (5,1):10, (5,2):8, (5,3):4, + I, d = multidict({1: 80, 2: 270, 3: 250, 4: 160, 5: 180}) # demand + J, M = multidict({1: 500, 2: 500, 3: 500}) # capacity + c = {(1, 1): 4, (1, 2): 6, (1, 3): 9, # cost + (2, 1): 5, (2, 2): 4, (2, 3): 7, + (3, 1): 6, (3, 2): 3, (3, 3): 4, + (4, 1): 8, (4, 2): 5, (4, 3): 3, + (5, 1): 10, (5, 2): 8, (5, 3): 4, } - return I,J,c,d,M + return I, J, c, d, M def make_inst2(): """creates example data set 2""" - I,d = multidict({1:45, 2:20, 3:30 , 4:30}) # demand - J,M = multidict({1:35, 2:50, 3:40}) # capacity - c = {(1,1):8, (1,2):9, (1,3):14 , # {(customer,factory) : cost} - (2,1):6, (2,2):12, (2,3):9 , - (3,1):10, (3,2):13, (3,3):16 , - (4,1):9, (4,2):7, (4,3):5 , + I, d = multidict({1: 45, 2: 20, 3: 30, 4: 30}) # demand + J, M = multidict({1: 35, 2: 50, 3: 40}) # capacity + c = {(1, 1): 8, (1, 2): 9, (1, 3): 14, # {(customer,factory) : cost} + (2, 1): 6, (2, 2): 12, (2, 3): 9, + (3, 1): 10, (3, 2): 13, (3, 3): 16, + (4, 1): 9, (4, 2): 7, (4, 3): 5, } - return I,J,c,d,M + return I, J, c, d, M if __name__ == "__main__": - I,J,c,d,M = make_inst1(); + I, J, c, d, M = make_inst1(); # I,J,c,d,M = make_inst2(); - model = transp(I,J,c,d,M) + model = transp(I, J, c, d, M) model.optimize() print("Optimal value:", model.getObjVal()) @@ -82,6 +83,6 @@ def make_inst2(): EPS = 1.e-6 x = model.data - for (i,j) in x: - if model.getVal(x[i,j]) > EPS: - print("sending quantity %10s from factory %3s to customer %3s" % (model.getVal(x[i,j]),j,i)) + for (i, j) in x: + if model.getVal(x[i, j]) > EPS: + print("sending quantity %10s from factory %3s to customer %3s" % (model.getVal(x[i, j]), j, i)) diff --git a/examples/finished/transp_nofn.py b/examples/finished/transp_nofn.py index d6ac409e1..b54c4cfd2 100644 --- a/examples/finished/transp_nofn.py +++ b/examples/finished/transp_nofn.py @@ -1,5 +1,5 @@ ##@file transp_nofn.py -#@brief a model for the transportation problem +# @brief a model for the transportation problem """ Model for solving a transportation problem: minimize the total transportation cost for satisfying demand at @@ -15,19 +15,19 @@ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ -from pyscipopt import Model, quicksum, multidict +from pyscipopt import Model, quicksum -d = {1:80, 2:270, 3:250 , 4:160, 5:180} # demand +d = {1: 80, 2: 270, 3: 250, 4: 160, 5: 180} # demand I = d.keys() -M = {1:500, 2:500, 3:500} # capacity +M = {1: 500, 2: 500, 3: 500} # capacity J = M.keys() -c = {(1,1):4, (1,2):6, (1,3):9, # cost - (2,1):5, (2,2):4, (2,3):7, - (3,1):6, (3,2):3, (3,3):4, - (4,1):8, (4,2):5, (4,3):3, - (5,1):10, (5,2):8, (5,3):4, +c = {(1, 1): 4, (1, 2): 6, (1, 3): 9, # cost + (2, 1): 5, (2, 2): 4, (2, 3): 7, + (3, 1): 6, (3, 2): 3, (3, 3): 4, + (4, 1): 8, (4, 2): 5, (4, 3): 3, + (5, 1): 10, (5, 2): 8, (5, 3): 4, } model = Model("transportation") @@ -37,18 +37,18 @@ for i in I: for j in J: - x[i,j] = model.addVar(vtype="C", name="x(%s,%s)" % (i,j)) + x[i, j] = model.addVar(vtype="C", name="x(%s,%s)" % (i, j)) # Demand constraints for i in I: - model.addCons(sum(x[i,j] for j in J if (i,j) in x) == d[i], name="Demand(%s)" % i) + model.addCons(sum(x[i, j] for j in J if (i, j) in x) == d[i], name="Demand(%s)" % i) # Capacity constraints for j in J: - model.addCons(sum(x[i,j] for i in I if (i,j) in x) <= M[j], name="Capacity(%s)" % j) + model.addCons(sum(x[i, j] for i in I if (i, j) in x) <= M[j], name="Capacity(%s)" % j) # Objective -model.setObjective(quicksum(c[i,j]*x[i,j] for (i,j) in x), "minimize") +model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), "minimize") model.optimize() @@ -56,6 +56,6 @@ EPS = 1.e-6 -for (i,j) in x: - if model.getVal(x[i,j]) > EPS: - print("sending quantity %10s from factory %3s to customer %3s" % (model.getVal(x[i,j]),j,i)) +for (i, j) in x: + if model.getVal(x[i, j]) > EPS: + print("sending quantity %10s from factory %3s to customer %3s" % (model.getVal(x[i, j]), j, i)) diff --git a/examples/finished/tsp.py b/examples/finished/tsp.py index 2d45d975b..3c9ca46d5 100644 --- a/examples/finished/tsp.py +++ b/examples/finished/tsp.py @@ -1,5 +1,5 @@ ##@file tsp.py -#@brief solve the traveling salesman problem +# @brief solve the traveling salesman problem """ minimize the travel cost for visiting n customers exactly once approach: @@ -13,10 +13,12 @@ """ import math import random + import networkx from pyscipopt import Model, quicksum -def solve_tsp(V,c): + +def solve_tsp(V, c): """solve_tsp -- solve the traveling salesman problem - start with assignment model - add cuts until there are no sub-cycles @@ -34,11 +36,10 @@ def addcut(cut_edges): return False model.freeTransform() for S in Components: - model.addCons(quicksum(x[i,j] for i in S for j in S if j>i) <= len(S)-1) - print("cut: len(%s) <= %s" % (S,len(S)-1)) + model.addCons(quicksum(x[i, j] for i in S for j in S if j > i) <= len(S) - 1) + print("cut: len(%s) <= %s" % (S, len(S) - 1)) return True - def addcut2(cut_edges): G = networkx.Graph() G.add_edges_from(cut_edges) @@ -49,64 +50,65 @@ def addcut2(cut_edges): model.freeTransform() for S in Components: T = set(V) - set(S) - print("S:",S) - print("T:",T) - model.addCons(quicksum(x[i,j] for i in S for j in T if j>i) + - quicksum(x[i,j] for i in T for j in S if j>i) >= 2) - print("cut: %s >= 2" % "+".join([("x[%s,%s]" % (i,j)) for i in S for j in T if j>i])) + print("S:", S) + print("T:", T) + model.addCons(quicksum(x[i, j] for i in S for j in T if j > i) + + quicksum(x[i, j] for i in T for j in S if j > i) >= 2) + print("cut: %s >= 2" % "+".join([("x[%s,%s]" % (i, j)) for i in S for j in T if j > i])) return True # main part of the solution process: model = Model("tsp") - model.hideOutput() # silent/verbose mode + model.hideOutput() # silent/verbose mode x = {} for i in V: for j in V: if j > i: - x[i,j] = model.addVar(ub=1, name="x(%s,%s)"%(i,j)) + x[i, j] = model.addVar(ub=1, name="x(%s,%s)" % (i, j)) for i in V: - model.addCons(quicksum(x[j,i] for j in V if j < i) + \ - quicksum(x[i,j] for j in V if j > i) == 2, "Degree(%s)"%i) + model.addCons(quicksum(x[j, i] for j in V if j < i) + \ + quicksum(x[i, j] for j in V if j > i) == 2, "Degree(%s)" % i) - model.setObjective(quicksum(c[i,j]*x[i,j] for i in V for j in V if j > i), "minimize") + model.setObjective(quicksum(c[i, j] * x[i, j] for i in V for j in V if j > i), "minimize") EPS = 1.e-6 isMIP = False while True: model.optimize() edges = [] - for (i,j) in x: - if model.getVal(x[i,j]) > EPS: - edges.append( (i,j) ) + for (i, j) in x: + if model.getVal(x[i, j]) > EPS: + edges.append((i, j)) if addcut(edges) == False: - if isMIP: # integer variables, components connected: solution found + if isMIP: # integer variables, components connected: solution found break model.freeTransform() - for (i,j) in x: # all components connected, switch to integer model - model.chgVarType(x[i,j], "B") + for (i, j) in x: # all components connected, switch to integer model + model.chgVarType(x[i, j], "B") isMIP = True - return model.getObjVal(),edges + return model.getObjVal(), edges -def distance(x1,y1,x2,y2): +def distance(x1, y1, x2, y2): """distance: euclidean distance between (x1,y1) and (x2,y2)""" - return math.sqrt((x2-x1)**2 + (y2-y1)**2) + return math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2) + def make_data(n): """make_data: compute matrix distance based on euclidean distance""" - V = range(1,n+1) - x = dict([(i,random.random()) for i in V]) - y = dict([(i,random.random()) for i in V]) + V = range(1, n + 1) + x = dict([(i, random.random()) for i in V]) + y = dict([(i, random.random()) for i in V]) c = {} for i in V: for j in V: if j > i: - c[i,j] = distance(x[i],y[i],x[j],y[j]) - return V,c + c[i, j] = distance(x[i], y[i], x[j], y[j]) + return V, c if __name__ == "__main__": @@ -119,16 +121,17 @@ def make_data(n): n = 200 seed = 1 random.seed(seed) - V,c = make_data(n) + V, c = make_data(n) else: from read_tsplib import read_tsplib + try: - V,c,x,y = read_tsplib(sys.argv[1]) + V, c, x, y = read_tsplib(sys.argv[1]) except: - print("Cannot read TSPLIB file",sys.argv[1]) + print("Cannot read TSPLIB file", sys.argv[1]) exit(1) - obj,edges = solve_tsp(V,c) + obj, edges = solve_tsp(V, c) - print("\nOptimal tour:",edges) - print("Optimal cost:",obj) + print("\nOptimal tour:", edges) + print("Optimal cost:", obj) diff --git a/examples/finished/weber_soco.py b/examples/finished/weber_soco.py index fa4ddccd2..29bb95a68 100644 --- a/examples/finished/weber_soco.py +++ b/examples/finished/weber_soco.py @@ -1,11 +1,12 @@ ##@file weber_soco.py -#@brief model for solving the weber problem using soco. +# @brief model for solving the weber problem using soco. """ Copyright (c) by Joao Pedro PEDROSO, Masahiro MURAMATSU and Mikio KUBO, 2012 """ -from pyscipopt import Model, quicksum, multidict +from pyscipopt import Model, quicksum -def weber(I,x,y,w): + +def weber(I, x, y, w): """weber: model for solving the single source weber problem using soco. Parameters: - I: set of customers @@ -17,65 +18,68 @@ def weber(I,x,y,w): model = Model("weber") - X,Y,z,xaux,yaux = {},{},{},{},{} + X, Y, z, xaux, yaux = {}, {}, {}, {}, {} X = model.addVar(lb=-model.infinity(), vtype="C", name="X") Y = model.addVar(lb=-model.infinity(), vtype="C", name="Y") for i in I: - z[i] = model.addVar(vtype="C", name="z(%s)"%(i)) - xaux[i] = model.addVar(lb=-model.infinity(), vtype="C", name="xaux(%s)"%(i)) - yaux[i] = model.addVar(lb=-model.infinity(), vtype="C", name="yaux(%s)"%(i)) + z[i] = model.addVar(vtype="C", name="z(%s)" % (i)) + xaux[i] = model.addVar(lb=-model.infinity(), vtype="C", name="xaux(%s)" % (i)) + yaux[i] = model.addVar(lb=-model.infinity(), vtype="C", name="yaux(%s)" % (i)) for i in I: - model.addCons(xaux[i]*xaux[i] + yaux[i]*yaux[i] <= z[i]*z[i], "MinDist(%s)"%(i)) - model.addCons(xaux[i] == (x[i]-X), "xAux(%s)"%(i)) - model.addCons(yaux[i] == (y[i]-Y), "yAux(%s)"%(i)) + model.addCons(xaux[i] * xaux[i] + yaux[i] * yaux[i] <= z[i] * z[i], "MinDist(%s)" % (i)) + model.addCons(xaux[i] == (x[i] - X), "xAux(%s)" % (i)) + model.addCons(yaux[i] == (y[i] - Y), "yAux(%s)" % (i)) - model.setObjective(quicksum(w[i]*z[i] for i in I), "minimize") + model.setObjective(quicksum(w[i] * z[i] for i in I), "minimize") - model.data = X,Y,z + model.data = X, Y, z return model import random -def make_data(n,m): + + +def make_data(n, m): """creates example data set""" - I = range(1,n+1) - J = range(1,m+1) - x,y,w = {},{},{} + I = range(1, n + 1) + J = range(1, m + 1) + x, y, w = {}, {}, {} for i in I: - x[i] = random.randint(0,100) - y[i] = random.randint(0,100) - w[i] = random.randint(1,5) - return I,J,x,y,w + x[i] = random.randint(0, 100) + y[i] = random.randint(0, 100) + w[i] = random.randint(1, 5) + return I, J, x, y, w if __name__ == "__main__": - import sys + random.seed(3) n = 7 m = 1 - I,J,x,y,w = make_data(n,m) + I, J, x, y, w = make_data(n, m) print("data:") - print("%s\t%8s\t%8s\t%8s" % ("i","x[i]","y[i]","w[i]")) + print("%s\t%8s\t%8s\t%8s" % ("i", "x[i]", "y[i]", "w[i]")) for i in I: - print("%s\t%8g\t%8g\t%8g" % (i,x[i],y[i],w[i])) + print("%s\t%8g\t%8g\t%8g" % (i, x[i], y[i], w[i])) print - model = weber(I,x,y,w) + model = weber(I, x, y, w) model.optimize() - X,Y,z = model.data - print("Optimal value=",model.getObjVal()) - print("Selected position:",) - print("\t",(round(model.getVal(X)),round(model.getVal(Y)))) + X, Y, z = model.data + print("Optimal value=", model.getObjVal()) + print("Selected position:", ) + print("\t", (round(model.getVal(X)), round(model.getVal(Y)))) print print("Solution:") - print("%s\t%8s" % ("i","z[i]")) + print("%s\t%8s" % ("i", "z[i]")) for i in I: - print("%s\t%8g" % (i,model.getVal(z[i]))) + print("%s\t%8g" % (i, model.getVal(z[i]))) print - try: # plot the result using networkx and matplotlib + try: # plot the result using networkx and matplotlib import networkx as NX import matplotlib.pyplot as P + P.clf() G = NX.Graph() @@ -84,20 +88,19 @@ def make_data(n,m): position = {} for i in I: - position[i] = (x[i],y[i]) - position["D"] = (round(model.getVal(X)),round(model.getVal(Y))) + position[i] = (x[i], y[i]) + position["D"] = (round(model.getVal(X)), round(model.getVal(Y))) - NX.draw(G,pos=position,node_size=200,node_color="g",nodelist=I) - NX.draw(G,pos=position,node_size=400,node_color="w",nodelist=["D"],alpha=0.5) - #P.savefig("weber.pdf",format="pdf",dpi=300) + NX.draw(G, pos=position, node_size=200, node_color="g", nodelist=I) + NX.draw(G, pos=position, node_size=400, node_color="w", nodelist=["D"], alpha=0.5) + # P.savefig("weber.pdf",format="pdf",dpi=300) P.show() except ImportError: print("install 'networkx' and 'matplotlib' for plotting") - -def weber_MS(I,J,x,y,w): +def weber_MS(I, J, x, y, w): """weber -- model for solving the weber problem using soco (multiple source version). Parameters: - I: set of customers @@ -107,34 +110,32 @@ def weber_MS(I,J,x,y,w): - w[i]: weight of customer i Returns a model, ready to be solved. """ - M = max([((x[i]-x[j])**2 + (y[i]-y[j])**2) for i in I for j in I]) + M = max([((x[i] - x[j]) ** 2 + (y[i] - y[j]) ** 2) for i in I for j in I]) model = Model("weber - multiple source") - X,Y,v,u = {},{},{},{} - xaux,yaux,uaux = {},{},{} + X, Y, v, u = {}, {}, {}, {} + xaux, yaux, uaux = {}, {}, {} for j in J: - X[j] = model.addVar(lb=-model.infinity(), vtype="C", name="X(%s)"%j) - Y[j] = model.addVar(lb=-model.infinity(), vtype="C", name="Y(%s)"%j) + X[j] = model.addVar(lb=-model.infinity(), vtype="C", name="X(%s)" % j) + Y[j] = model.addVar(lb=-model.infinity(), vtype="C", name="Y(%s)" % j) for i in I: - v[i,j] = model.addVar(vtype="C", name="v(%s,%s)"%(i,j)) - u[i,j] = model.addVar(vtype="B", name="u(%s,%s)"%(i,j)) - xaux[i,j] = model.addVar(lb=-model.infinity(), vtype="C", name="xaux(%s,%s)"%(i,j)) - yaux[i,j] = model.addVar(lb=-model.infinity(), vtype="C", name="yaux(%s,%s)"%(i,j)) - uaux[i,j] = model.addVar(vtype="C", name="uaux(%s,%s)"%(i,j)) - - + v[i, j] = model.addVar(vtype="C", name="v(%s,%s)" % (i, j)) + u[i, j] = model.addVar(vtype="B", name="u(%s,%s)" % (i, j)) + xaux[i, j] = model.addVar(lb=-model.infinity(), vtype="C", name="xaux(%s,%s)" % (i, j)) + yaux[i, j] = model.addVar(lb=-model.infinity(), vtype="C", name="yaux(%s,%s)" % (i, j)) + uaux[i, j] = model.addVar(vtype="C", name="uaux(%s,%s)" % (i, j)) for i in I: - model.addCons(quicksum(u[i,j] for j in J) == 1, "Assign(%s)"%i) + model.addCons(quicksum(u[i, j] for j in J) == 1, "Assign(%s)" % i) for j in J: - model.addCons(xaux[i,j]*xaux[i,j] + yaux[i,j]*yaux[i,j] <= v[i,j]*v[i,j], "MinDist(%s,%s)"%(i,j)) - model.addCons(xaux[i,j] == (x[i]-X[j]), "xAux(%s,%s)"%(i,j)) - model.addCons(yaux[i,j] == (y[i]-Y[j]), "yAux(%s,%s)"%(i,j)) - model.addCons(uaux[i,j] >= v[i,j] - M*(1-u[i,j]), "uAux(%s,%s)"%(i,j)) + model.addCons(xaux[i, j] * xaux[i, j] + yaux[i, j] * yaux[i, j] <= v[i, j] * v[i, j], + "MinDist(%s,%s)" % (i, j)) + model.addCons(xaux[i, j] == (x[i] - X[j]), "xAux(%s,%s)" % (i, j)) + model.addCons(yaux[i, j] == (y[i] - Y[j]), "yAux(%s,%s)" % (i, j)) + model.addCons(uaux[i, j] >= v[i, j] - M * (1 - u[i, j]), "uAux(%s,%s)" % (i, j)) - model.setObjective(quicksum(w[i]*uaux[i,j] for i in I for j in J), "minimize") + model.setObjective(quicksum(w[i] * uaux[i, j] for i in I for j in J), "minimize") - - model.data = X,Y,v,u + model.data = X, Y, v, u return model @@ -142,40 +143,41 @@ def weber_MS(I,J,x,y,w): random.seed(3) n = 7 m = 1 - I,J,x,y,w = make_data(n,m) - model = weber_MS(I,J,x,y,w) + I, J, x, y, w = make_data(n, m) + model = weber_MS(I, J, x, y, w) model.optimize() - X,Y,w,z = model.data - print("Optimal value:",model.getObjVal()) + X, Y, w, z = model.data + print("Optimal value:", model.getObjVal()) print("Selected positions:") for j in J: - print("\t",(model.getVal(X[j]),model.getVal(Y[j]))) - for (i,j) in sorted(w.keys()): - print("\t",(i,j),model.getVal(w[i,j]),model.getVal(z[i,j])) + print("\t", (model.getVal(X[j]), model.getVal(Y[j]))) + for (i, j) in sorted(w.keys()): + print("\t", (i, j), model.getVal(w[i, j]), model.getVal(z[i, j])) EPS = 1.e-4 - edges = [(i,j) for (i,j) in z if model.getVal(z[i,j]) > EPS] + edges = [(i, j) for (i, j) in z if model.getVal(z[i, j]) > EPS] - try: # plot the result using networkx and matplotlib + try: # plot the result using networkx and matplotlib import networkx as NX import matplotlib.pyplot as P + P.clf() G = NX.Graph() G.add_nodes_from(I) - G.add_nodes_from("%s"%j for j in J) # for distinguishing J from I, make nodes as strings - for (i,j) in edges: - G.add_edge(i,"%s"%j) + G.add_nodes_from("%s" % j for j in J) # for distinguishing J from I, make nodes as strings + for (i, j) in edges: + G.add_edge(i, "%s" % j) position = {} for i in I: - position[i] = (x[i],y[i]) + position[i] = (x[i], y[i]) for j in J: - position["%s"%j] = (model.getVal(X[j]),model.getVal(Y[j])) + position["%s" % j] = (model.getVal(X[j]), model.getVal(Y[j])) print(position) - NX.draw(G,position,node_color="g",nodelist=I) - NX.draw(G,position,node_color="w",nodelist=["%s"%j for j in J]) + NX.draw(G, position, node_color="g", nodelist=I) + NX.draw(G, position, node_color="w", nodelist=["%s" % j for j in J]) P.show() except ImportError: print("install 'networkx' and 'matplotlib' for plotting") diff --git a/examples/tutorial/even.py b/examples/tutorial/even.py index cd09ce2db..5fd4467a4 100644 --- a/examples/tutorial/even.py +++ b/examples/tutorial/even.py @@ -1,31 +1,33 @@ ##@file tutorial/even.py -#@brief Tutorial example to check whether values are even or odd +# @brief Tutorial example to check whether values are even or odd """ Public Domain, WTFNMFPL Public Licence """ -from pyscipopt import Model from pprint import pformat as pfmt +from pyscipopt import Model + example_values = [ - 0, - 1, - 1.5, - "helloworld", - 20, - 25, - -101, - -15., - -10, - -2**31, - -int(2**31), - "2**31-1", - int(2**31-1), - int(2**63)-1 + 0, + 1, + 1.5, + "helloworld", + 20, + 25, + -101, + -15., + -10, + -2 ** 31, + -int(2 ** 31), + "2**31-1", + int(2 ** 31 - 1), + int(2 ** 63) - 1 ] verbose = False -#verbose = True # uncomment for additional info on variables! -sdic = {0: "even", 1: "odd"} # remainder to 2 +# verbose = True # uncomment for additional info on variables! +sdic = {0: "even", 1: "odd"} # remainder to 2 + def parity(number): """ @@ -44,7 +46,7 @@ def parity(number): """ sval = -1 if verbose: - print(80*"*") + print(80 * "*") try: assert number == int(round(number)) m = Model() @@ -55,7 +57,7 @@ def parity(number): # since 0 is the default lb. To allow for negative values, give None # as lower bound. # (None means -infinity as lower bound and +infinity as upper bound) - x = m.addVar("x", vtype="I", lb=None, ub=None) #ub=None is default + x = m.addVar("x", vtype="I", lb=None, ub=None) # ub=None is default n = m.addVar("n", vtype="I", lb=None) s = m.addVar("s", vtype="B") # CAVEAT: if number is negative, x's lower bound must be None @@ -63,18 +65,18 @@ def parity(number): # there is no feasible solution (trivial) but the program # does not highlight which constraints conflict. - m.addCons(x==number) + m.addCons(x == number) # minimize the difference between the number and twice a natural number - m.addCons(s == x-2*n) + m.addCons(s == x - 2 * n) m.setObjective(s) m.optimize() assert m.getStatus() == "optimal" - boolmod = m.getVal(s) == m.getVal(x)%2 + boolmod = m.getVal(s) == m.getVal(x) % 2 if verbose: for v in m.getVars(): - print("%*s: %d" % (fmtlen, v,m.getVal(v))) + print("%*s: %d" % (fmtlen, v, m.getVal(v))) print("%*d%%2 == %d?" % (fmtlen, m.getVal(x), m.getVal(s))) print("%*s" % (fmtlen, boolmod)) @@ -86,10 +88,11 @@ def parity(number): print("%*s is neither even nor odd!" % (fmtlen, number.__repr__())) finally: if verbose: - print(80*"*") + print(80 * "*") print("") return sval + if __name__ == "__main__": """ If positional arguments are given: @@ -99,6 +102,7 @@ def parity(number): """ import sys from ast import literal_eval as leval + try: # check parity for each positional arguments sys.argv[1] @@ -107,10 +111,10 @@ def parity(number): # check parity for each default example value values = example_values # format lenght, cosmetics - fmtlen = max([len(fmt) for fmt in pfmt(values,width=1).split('\n')]) + fmtlen = max([len(fmt) for fmt in pfmt(values, width=1).split('\n')]) for value in values: try: n = leval(value) - except (ValueError, SyntaxError): # for numbers or str w/ spaces + except (ValueError, SyntaxError): # for numbers or str w/ spaces n = value parity(n) diff --git a/examples/tutorial/logical.py b/examples/tutorial/logical.py index 9d437558d..1553ae181 100644 --- a/examples/tutorial/logical.py +++ b/examples/tutorial/logical.py @@ -1,5 +1,5 @@ ##@file tutorial/logical.py -#@brief Tutorial example on how to use AND/OR/XOR constraints. +# @brief Tutorial example on how to use AND/OR/XOR constraints. """ N.B.: standard SCIP XOR constraint works differently from AND/OR by design. The constraint is set with a boolean rhs instead of an integer resultant. @@ -11,14 +11,16 @@ from pyscipopt import Model from pyscipopt import quicksum + def _init(): model = Model() model.hideOutput() - x = model.addVar("x","B") - y = model.addVar("y","B") - z = model.addVar("z","B") + x = model.addVar("x", "B") + y = model.addVar("y", "B") + z = model.addVar("z", "B") return model, x, y, z + def _optimize(name, m): m.optimize() print("* %s constraint *" % name) @@ -36,50 +38,53 @@ def _optimize(name, m): print("* No variable is printed if model status is not optimal") print("") + def and_constraint(v=1, sense="minimize"): """ AND constraint """ - assert v in [0,1], "v must be 0 or 1 instead of %s" % v.__repr__() + assert v in [0, 1], "v must be 0 or 1 instead of %s" % v.__repr__() model, x, y, z = _init() r = model.addVar("r", "B") - model.addConsAnd([x,y,z], r) - model.addCons(x==v) + model.addConsAnd([x, y, z], r) + model.addCons(x == v) model.setObjective(r, sense=sense) _optimize("AND", model) def or_constraint(v=0, sense="maximize"): """ OR constraint""" - assert v in [0,1], "v must be 0 or 1 instead of %s" % v.__repr__() + assert v in [0, 1], "v must be 0 or 1 instead of %s" % v.__repr__() model, x, y, z = _init() r = model.addVar("r", "B") - model.addConsOr([x,y,z], r) - model.addCons(x==v) + model.addConsOr([x, y, z], r) + model.addCons(x == v) model.setObjective(r, sense=sense) _optimize("OR", model) + def xors_constraint(v=1): """ XOR (r as boolean) standard constraint""" - assert v in [0,1], "v must be 0 or 1 instead of %s" % v.__repr__() + assert v in [0, 1], "v must be 0 or 1 instead of %s" % v.__repr__() model, x, y, z = _init() r = True - model.addConsXor([x,y,z], r) - model.addCons(x==v) + model.addConsXor([x, y, z], r) + model.addCons(x == v) _optimize("Standard XOR (as boolean)", model) + def xorc_constraint(v=0, sense="maximize"): """ XOR (r as variable) custom constraint""" - assert v in [0,1], "v must be 0 or 1 instead of %s" % v.__repr__() + assert v in [0, 1], "v must be 0 or 1 instead of %s" % v.__repr__() model, x, y, z = _init() r = model.addVar("r", "B") - n = model.addVar("n", "I") # auxiliary - model.addCons(r+quicksum([x,y,z]) == 2*n) - model.addCons(x==v) + n = model.addVar("n", "I") # auxiliary + model.addCons(r + quicksum([x, y, z]) == 2 * n) + model.addCons(x == v) model.setObjective(r, sense=sense) _optimize("Custom XOR (as variable)", model) + if __name__ == "__main__": and_constraint() or_constraint() xors_constraint() xorc_constraint() - diff --git a/examples/tutorial/puzzle.py b/examples/tutorial/puzzle.py index c507b77fc..b386cec58 100644 --- a/examples/tutorial/puzzle.py +++ b/examples/tutorial/puzzle.py @@ -1,5 +1,5 @@ ##@file tutorial/puzzle.py -#@brief solve a simple puzzle using SCIP +# @brief solve a simple puzzle using SCIP """ On a beach there are octopuses, turtles and cranes. The total number of legs of all animals is 80, while the number of heads is 32. @@ -18,7 +18,7 @@ model.addCons(x + y + z == 32, name="Heads") # Set up constraint for number of legs -model.addCons(8*x + 4*y + 2*z == 80, name="Legs") +model.addCons(8 * x + 4 * y + 2 * z == 80, name="Legs") # Set objective function model.setObjective(x + y, "minimize") @@ -26,7 +26,7 @@ model.hideOutput() model.optimize() -#solution = model.getBestSol() +# solution = model.getBestSol() print("Optimal value:", model.getObjVal()) print((x.name, y.name, z.name), " = ", (model.getVal(x), model.getVal(y), model.getVal(z))) diff --git a/examples/unfinished/cutstock.py b/examples/unfinished/cutstock.py index 0be31a730..c011867db 100644 --- a/examples/unfinished/cutstock.py +++ b/examples/unfinished/cutstock.py @@ -21,7 +21,7 @@ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ -from pyscipopt import Model, quicksum, multidict +from pyscipopt import Model, quicksum LOG = True EPS = 1.e-6 diff --git a/examples/unfinished/eld.py b/examples/unfinished/eld.py index 959e89665..0ea3351d8 100644 --- a/examples/unfinished/eld.py +++ b/examples/unfinished/eld.py @@ -5,11 +5,11 @@ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ -from pyscipopt import Model, quicksum, multidict import math -import random from piecewise import convex_comb_sos +from pyscipopt import Model, quicksum, multidict + def cost(a,b,c,e,f,p_min,p): """cost: fuel cost based on "standard" parameters diff --git a/examples/unfinished/flp_nonlinear.py b/examples/unfinished/flp_nonlinear.py index 86f517d28..5602daf11 100644 --- a/examples/unfinished/flp_nonlinear.py +++ b/examples/unfinished/flp_nonlinear.py @@ -16,8 +16,10 @@ """ import math import random -from pyscipopt import Model, quicksum, multidict + from piecewise import * +from pyscipopt import Model, quicksum, multidict + def flp_nonlinear_mselect(I,J,d,M,f,c,K): """flp_nonlinear_mselect -- use multiple selection model diff --git a/examples/unfinished/flp_nonlinear_soco.py b/examples/unfinished/flp_nonlinear_soco.py index 33c721e84..e4ccd32b4 100644 --- a/examples/unfinished/flp_nonlinear_soco.py +++ b/examples/unfinished/flp_nonlinear_soco.py @@ -11,9 +11,10 @@ Copyright (c) by Joao Pedro PEDROSO, Masahiro MURAMATSU and Mikio KUBO, 2012 """ -import math import random -from pyscipopt import Model, quicksum, multidict + +from pyscipopt import Model, quicksum + def flp_nonlinear_soco(I,J,d,M,f,c): """flp_nonlinear_soco -- use @@ -52,7 +53,8 @@ def flp_nonlinear_soco(I,J,d,M,f,c): if __name__ == "__main__": - from flp_nonlinear import distance,make_data,example,read_orlib,read_cortinhal + from flp_nonlinear import make_data + # I,J,d,M,f,c,x_pos,y_pos = example() K = 100 random.seed(1) diff --git a/examples/unfinished/gpp.py b/examples/unfinished/gpp.py index cf8378454..466e181ad 100644 --- a/examples/unfinished/gpp.py +++ b/examples/unfinished/gpp.py @@ -4,7 +4,8 @@ Copyright (c) by Joao Pedro PEDROSO, Masahiro MURAMATSU and Mikio KUBO, 2012 """ -from pyscipopt import Model, quicksum, multidict +from pyscipopt import Model, quicksum + def gpp(V,E): """gpp -- model for the graph partitioning problem diff --git a/examples/unfinished/kcenter.py b/examples/unfinished/kcenter.py index 9d0419691..aef661167 100644 --- a/examples/unfinished/kcenter.py +++ b/examples/unfinished/kcenter.py @@ -6,7 +6,8 @@ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ -from pyscipopt import Model, quicksum, multidict +from pyscipopt import Model, quicksum + def kcenter(I,J,c,k): """kcenter -- minimize the maximum travel cost from customers to k facilities. diff --git a/examples/unfinished/kcenter_binary_search.py b/examples/unfinished/kcenter_binary_search.py index ab746d677..f20f1a744 100644 --- a/examples/unfinished/kcenter_binary_search.py +++ b/examples/unfinished/kcenter_binary_search.py @@ -7,7 +7,8 @@ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ -from pyscipopt import Model, quicksum, multidict +from pyscipopt import Model, quicksum + def kcover(I,J,c,k): """kcover -- minimize the number of uncovered customers from k facilities. diff --git a/examples/unfinished/lotsizing.py b/examples/unfinished/lotsizing.py index 210abb18a..b5b18fee5 100644 --- a/examples/unfinished/lotsizing.py +++ b/examples/unfinished/lotsizing.py @@ -8,8 +8,10 @@ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ import random + from pyscipopt import Model, quicksum + def mils(T,P,f,g,c,d,h,M): """ mils: standard formulation for the multi-item lot-sizing problem diff --git a/examples/unfinished/lotsizing_echelon.py b/examples/unfinished/lotsizing_echelon.py index dbcaad8fd..0cd0e3cfa 100644 --- a/examples/unfinished/lotsizing_echelon.py +++ b/examples/unfinished/lotsizing_echelon.py @@ -8,8 +8,10 @@ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ import random + from pyscipopt import Model, quicksum, multidict + def mils_standard(T,K,P,f,g,c,d,h,a,M,UB,phi): """ mils_standard: standard formulation for the multi-item, multi-stage lot-sizing problem diff --git a/examples/unfinished/portfolio_soco.py b/examples/unfinished/portfolio_soco.py index bcd25bbfb..bd70a70ba 100644 --- a/examples/unfinished/portfolio_soco.py +++ b/examples/unfinished/portfolio_soco.py @@ -5,9 +5,11 @@ Copyright (c) by Joao Pedro PEDROSO, Masahiro MURAMATSU and Mikio KUBO, 2012 """ +import math + from pyscipopt import Model, quicksum, multidict -import math + def phi_inv(p): """phi_inv: inverse of gaussian (normal) CDF Source: diff --git a/examples/unfinished/staff_sched.py b/examples/unfinished/staff_sched.py index a30f35463..962482095 100644 --- a/examples/unfinished/staff_sched.py +++ b/examples/unfinished/staff_sched.py @@ -3,7 +3,6 @@ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ -import random from pyscipopt import Model, quicksum, multidict def staff(I,T,N,J,S,c,b): diff --git a/examples/unfinished/staff_sched_mo.py b/examples/unfinished/staff_sched_mo.py index 623e0f118..fc7ad2e68 100644 --- a/examples/unfinished/staff_sched_mo.py +++ b/examples/unfinished/staff_sched_mo.py @@ -7,8 +7,8 @@ Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 """ -import random -from pyscipopt import Model, quicksum, multidict +from pyscipopt import Model, quicksum + def staff_mo(I,T,N,J,S,c,b): """ @@ -135,7 +135,7 @@ def solve_segment(I,T,N,J,S,c,b): if __name__ == "__main__": from pareto_front import pareto_front - from staff_sched import make_data,make_data_trick + from staff_sched import make_data_trick # I,T,N,J,S,c,b = make_data() I,T,N,J,S,c,b = make_data_trick() cand_seg = solve_segment(I,T,N,J,S,c,b) diff --git a/examples/unfinished/tsp_flow.py b/examples/unfinished/tsp_flow.py index d49dfecc6..b6e8e8054 100644 --- a/examples/unfinished/tsp_flow.py +++ b/examples/unfinished/tsp_flow.py @@ -12,8 +12,9 @@ """ import math import random -#import networkx todo -from pyscipopt import Model, quicksum, multidict + +# import networkx todo +from pyscipopt import Model, quicksum EPS = 1.e-6 diff --git a/examples/unfinished/tsp_lazy.py b/examples/unfinished/tsp_lazy.py index 07f3911f0..c6fa8fcaf 100644 --- a/examples/unfinished/tsp_lazy.py +++ b/examples/unfinished/tsp_lazy.py @@ -10,8 +10,10 @@ """ import math import random + import networkx -from pyscipopt import Model, Conshdlr, quicksum, SCIP_RESULT, SCIP_PRESOLTIMING, SCIP_PROPTIMING, SCIP_PARAMSETTING +from pyscipopt import Model, Conshdlr, quicksum, SCIP_RESULT, SCIP_PRESOLTIMING, SCIP_PROPTIMING + class TSPconshdlr(Conshdlr): diff --git a/examples/unfinished/tsp_mo.py b/examples/unfinished/tsp_mo.py index 1b0430264..5cd01693e 100644 --- a/examples/unfinished/tsp_mo.py +++ b/examples/unfinished/tsp_mo.py @@ -11,7 +11,9 @@ import math import random -from pyscipopt import Model, quicksum, multidict + +from pyscipopt import quicksum + def optimize(model,cand): """optimize: function for solving the model, updating candidate solutions' list diff --git a/examples/unfinished/tsptw.py b/examples/unfinished/tsptw.py index 05d7a04db..5ee4ff1eb 100644 --- a/examples/unfinished/tsptw.py +++ b/examples/unfinished/tsptw.py @@ -11,7 +11,9 @@ """ import math import random -from pyscipopt import Model, quicksum, multidict + +from pyscipopt import Model, quicksum + def mtztw(n,c,e,l): """mtzts: model for the traveling salesman problem with time windows diff --git a/examples/unfinished/vrp.py b/examples/unfinished/vrp.py index a3e5f6d76..b34d7ed9f 100644 --- a/examples/unfinished/vrp.py +++ b/examples/unfinished/vrp.py @@ -9,10 +9,12 @@ """ import math import random + import networkx -from pyscipopt import Model, quicksum, multidict +from pyscipopt import Model, quicksum + -def solve_vrp(V,c,m,q,Q): +def solve_vrp(V, c, m, q, Q): """solve_vrp -- solve the vehicle routing problem. - start with assignment model (depot has a special status) - add cuts until all components of the graph are connected @@ -38,10 +40,11 @@ def addcut(cut_edges): for S in Components: S_card = len(S) q_sum = sum(q[i] for i in S) - NS = int(math.ceil(float(q_sum)/Q)) - S_edges = [(i,j) for i in S for j in S if i= 3 and (len(S_edges) >= S_card or NS > 1): - add = model.addCons(quicksum(x[i,j] for i in S for j in S if j > i) <= S_card-NS) + print("Adding a cut") + add = model.addCons(quicksum(x[i, j] for i in S for j in S if j > i) <= S_card - NS) cut = True return cut @@ -50,17 +53,17 @@ def addcut(cut_edges): x = {} for i in V: for j in V: - if j > i and i == V[0]: # depot - x[i,j] = model.addVar(ub=2, vtype="I", name="x(%s,%s)"%(i,j)) + if j > i and i == V[0]: # depot + x[i, j] = model.addVar(ub=2, vtype="I", name="x(%s,%s)" % (i, j)) elif j > i: - x[i,j] = model.addVar(ub=1, vtype="I", name="x(%s,%s)"%(i,j)) - - model.addCons(quicksum(x[V[0],j] for j in V[1:]) == 2*m, "DegreeDepot") + x[i, j] = model.addVar(ub=1, vtype="I", name="x(%s,%s)" % (i, j)) + + model.addCons(quicksum(x[V[0], j] for j in V[1:]) == 2 * m, "DegreeDepot") for i in V[1:]: - model.addCons(quicksum(x[j,i] for j in V if j < i) + - quicksum(x[i,j] for j in V if j > i) == 2, "Degree(%s)"%i) + model.addCons(quicksum(x[j, i] for j in V if j < i) + + quicksum(x[i, j] for j in V if j > i) == 2, "Degree(%s)" % i) - model.setObjective(quicksum(c[i,j]*x[i,j] for i in V for j in V if j>i), "minimize") + model.setObjective(quicksum(c[i, j] * x[i, j] for i in V for j in V if j > i), "minimize") model.hideOutput() @@ -68,44 +71,44 @@ def addcut(cut_edges): while True: model.optimize() edges = [] - for (i,j) in x: - if model.getVal(x[i,j]) > EPS: + for (i, j) in x: + if model.getVal(x[i, j]) > EPS: if i != V[0] and j != V[0]: - edges.append((i,j)) + edges.append((i, j)) + model.freeTransform() if addcut(edges) == False: break - return model.getObjVal(),edges + return model.getObjVal(), edges -def distance(x1,y1,x2,y2): +def distance(x1, y1, x2, y2): """distance: euclidean distance between (x1,y1) and (x2,y2)""" - return math.sqrt((x2-x1)**2 + (y2-y1)**2) + return math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2) + def make_data(n): """make_data: compute matrix distance based on euclidean distance""" - V = range(1,n+1) - x = dict([(i,random.random()) for i in V]) - y = dict([(i,random.random()) for i in V]) - c,q = {},{} + V = range(1, n + 1) + x = dict([(i, random.random()) for i in V]) + y = dict([(i, random.random()) for i in V]) + c, q = {}, {} Q = 100 for i in V: - q[i] = random.randint(10,20) + q[i] = random.randint(10, 20) for j in V: if j > i: - c[i,j] = distance(x[i],y[i],x[j],y[j]) - return V,c,q,Q + c[i, j] = distance(x[i], y[i], x[j], y[j]) + return V, c, q, Q if __name__ == "__main__": - import sys - - n = 19 + n = 5 m = 3 seed = 1 random.seed(seed) - V,c,q,Q = make_data(n) - z,edges = solve_vrp(V,c,m,q,Q) - print("Optimal solution:",z) + V, c, q, Q = make_data(n) + z, edges = solve_vrp(V, c, m, q, Q) + print("Optimal solution:", z) print("Edges in the solution:") print(sorted(edges)) diff --git a/examples/unfinished/vrp_lazy.py b/examples/unfinished/vrp_lazy.py index 307706633..7b88249ea 100644 --- a/examples/unfinished/vrp_lazy.py +++ b/examples/unfinished/vrp_lazy.py @@ -9,8 +9,10 @@ """ import math import random + import networkx -from pyscipopt import Model, Conshdlr, quicksum, multidict, SCIP_RESULT, SCIP_PRESOLTIMING, SCIP_PROPTIMING +from pyscipopt import Model, Conshdlr, quicksum, SCIP_RESULT, SCIP_PRESOLTIMING, SCIP_PROPTIMING + class VRPconshdlr(Conshdlr): diff --git a/src/pyscipopt/recipes/README.md b/src/pyscipopt/recipes/README.md index 9942db0c3..91184b144 100644 --- a/src/pyscipopt/recipes/README.md +++ b/src/pyscipopt/recipes/README.md @@ -1,3 +1,6 @@ # Recipes sub-package -This sub-package provides a set of functions for common usecases for pyscipopt. This sub-package is for all functions that don't necessarily reflect the core functionality of SCIP, but are useful for working with the solver. The functions implemented in this sub-package might not be the most efficient way to solve/formulate a problem but would provide a good starting point. +This sub-package provides a set of functions for common usecases for pyscipopt. This sub-package is for all functions +that don't necessarily reflect the core functionality of SCIP, but are useful for working with the solver. The functions +implemented in this sub-package might not be the most efficient way to solve/formulate a problem but would provide a +good starting point. diff --git a/src/pyscipopt/recipes/nonlinear.py b/src/pyscipopt/recipes/nonlinear.py index f804e68a8..d4f485786 100644 --- a/src/pyscipopt/recipes/nonlinear.py +++ b/src/pyscipopt/recipes/nonlinear.py @@ -1,18 +1,18 @@ from pyscipopt import Model + def set_nonlinear_objective(model: Model, expr, sense="minimize"): - """ - Takes a nonlinear expression and performs an epigraph reformulation. - """ - - assert expr.degree() > 1, "For linear objectives, please use the setObjective method." - new_obj = model.addVar(lb=-float("inf"),obj=1) - if sense == "minimize": - model.addCons(expr <= new_obj) - model.setMinimize() - elif sense == "maximize": - model.addCons(expr >= new_obj) - model.setMaximize() - else: - raise Warning("unrecognized optimization sense: %s" % sense) - \ No newline at end of file + """ + Takes a nonlinear expression and performs an epigraph reformulation. + """ + + assert expr.degree() > 1, "For linear objectives, please use the setObjective method." + new_obj = model.addVar(lb=-float("inf"), obj=1) + if sense == "minimize": + model.addCons(expr <= new_obj) + model.setMinimize() + elif sense == "maximize": + model.addCons(expr >= new_obj) + model.setMaximize() + else: + raise Warning("unrecognized optimization sense: %s" % sense)