# source:projects/Dippy/trunk/examples/facility_decomp.py@524

Last change on this file since 524 was 524, checked in by mosu001, 3 years ago
File size: 6.9 KB
Line
1import coinor.pulp as pulp
2from pulp import *
3import coinor.dippy as dippy
4
5from math import floor, ceil
6
7tol = pow(pow(2, -24), 2.0 / 3.0)
8
9from facility_ex1 import REQUIREMENT, PRODUCTS, LOCATIONS, CAPACITY
10##from facility_test1 import REQUIREMENT, PRODUCTS, LOCATIONS, CAPACITY
11##from facility_test2 import REQUIREMENT, PRODUCTS, LOCATIONS, CAPACITY
12##from facility_test6 import REQUIREMENT, PRODUCTS, LOCATIONS, CAPACITY
13
14prob = dippy.DipProblem("Facility_Location")
15
16assign_vars = LpVariable.dicts("AtLocation",
17              [(i, j) for i in LOCATIONS
18                      for j in PRODUCTS],
19              0, 1, LpBinary)
20use_vars    = LpVariable.dicts("UseLocation",
21              LOCATIONS, 0, 1, LpBinary)
22waste_vars  = LpVariable.dicts("Waste",
23              LOCATIONS, 0, CAPACITY)
24
25# objective: minimise waste
26prob += lpSum(waste_vars[i] for i in LOCATIONS), "min"
27
28# assignment constraints
29for j in PRODUCTS:
30    prob += lpSum(assign_vars[(i, j)] for i in LOCATIONS) == 1
31
32# Aggregate capacity constraints
33for i in LOCATIONS:
34    prob.relaxation[i] += lpSum(assign_vars[(i, j)] * REQUIREMENT[j]
35                                for j in PRODUCTS) + waste_vars[i] \
36                                    == CAPACITY * use_vars[i]
37
38# Disaggregate capacity constraints
39for i in LOCATIONS:
40    for j in PRODUCTS:
41        prob.relaxation[i] += assign_vars[(i, j)] <= use_vars[i]
42
43# Ordering constraints
44for index, location in enumerate(LOCATIONS):
45    if index > 0:
46        prob += use_vars[LOCATIONS[index-1]] >= use_vars[location]
47
48# Anti-symmetry branches
49def choose_antisymmetry_branch(prob, sol):
50    num_locations = sum(sol[use_vars[i]] for i in LOCATIONS)
51    up   = ceil(num_locations)  # Round up to next nearest integer
52    down = floor(num_locations) # Round down
53    if  (up - num_locations   > tol) \
54    and (num_locations - down > tol): # Is fractional?
55        # Down branch: provide upper bounds, lower bounds are default
56        down_branch_ub = [(use_vars[LOCATIONS[n]], 0) for
57                                n in range(int(down), len(LOCATIONS))]
58        # Up branch: provide lower bounds, upper bounds are default
59        up_branch_lb = [(use_vars[LOCATIONS[n]], 1) for
60                                n in range(0, int(up))]
61        # Return the advanced branch to DIP
62        return ([], down_branch_ub, up_branch_lb, [])
63
64prob.branch_method = choose_antisymmetry_branch
65
66def solve_subproblem(prob, index, redCosts, convexDual):
67   loc = index
68
69   # Calculate effective objective coefficient of products
70   effs = {}
71   for j in PRODUCTS:
72       effs[j] = redCosts[assign_vars[(loc, j)]] \
73               - redCosts[waste_vars[loc]] * REQUIREMENT[j]
74
75   avars = [assign_vars[(loc, j)] for j in PRODUCTS]
76   obj = [-effs[j] for j in PRODUCTS]
77   weights = [REQUIREMENT[j] for j in PRODUCTS]
78
79   # Use 0-1 KP to max. total effective value of products at location
80   z, solution = knapsack01(obj, weights, CAPACITY)
81
82   # Get the reduced cost of the knapsack solution and waste
83   rc = redCosts[use_vars[loc]] -z + \
84        redCosts[waste_vars[loc]] * CAPACITY
85   waste = CAPACITY - sum(weights[i] for i in solution)
86   rc += redCosts[waste_vars[loc]] * waste
87
88   # Return the solution if the reduced cost is low enough
89   if rc < -tol: # The location is used and "useful"
90       if rc - convexDual < -tol:
91           var_values = [(avars[i], 1) for i in solution]
92           var_values.append((use_vars[loc], 1))
93           var_values.append((waste_vars[loc], waste))
94
95           dv = dippy.DecompVar(var_values, rc - convexDual, waste)
96           return [dv]
97
98   elif -convexDual < -tol: # An empty location is "useful"
99           var_values = []
100
101           dv = dippy.DecompVar(var_values, -convexDual, 0.0)
102           return [dv]
103
104   return []
105
106prob.relaxed_solver = solve_subproblem
107
108def knapsack01(obj, weights, capacity):
109    """ 0/1 knapsack solver, maximizes profit. weights and capacity integer """
110    assert len(obj) == len(weights)
111    n = len(obj)
112
113    if n == 0:
114        return 0, []
115
116    c = [[0]*(capacity+1) for i in range(n)]
117    added = [[False]*(capacity+1) for i in range(n)]
118    # c [items, remaining capacity]
119    # important: this code assumes strictly positive objective values
120    for i in range(n):
121        for j in range(capacity+1):
122            if (weights[i] > j):
123                c[i][j] = c[i-1][j]
124            else:
125                c_add = obj[i] + c[i-1][j-weights[i]]
129                else:
130                    c[i][j] = c[i-1][j]
131
132    # backtrack to find solution
133    i = n-1
134    j = capacity
135
136    solution = []
137    while i >= 0 and j >= 0:
139            solution.append(i)
140            j -= weights[i]
141        i -= 1
142
143    return c[n-1][capacity], solution
144
145def first_fit_heuristic():
146    # Sort the items in descending weight order
147    productReqs = [(REQUIREMENT[j],j) for j in PRODUCTS]
148    productReqs.sort(reverse=True)
149
150    # Add items to locations, fitting in as much
151    # as possible at each location.
152    allLocations = []
153    while len(productReqs) > 0:
154        waste = CAPACITY
155        currentLocation = []
156        j = 0
157        while j < len(productReqs):
158            # Can we fit this product?
159            if productReqs[j][0] <= waste:
160                currentLocation.append(productReqs[j][1]) # index
161                waste -= productReqs[j][0] # requirement
162                productReqs.pop(j)
163            else:
164                # Try to fit next item
165                j += 1
166        allLocations.append((currentLocation, waste))
167    # Return a list of tuples: ([products],waste)
168    return allLocations
169
170
171def first_fit(prob):
172    locations = first_fit_heuristic()
173    bvs = []
174    index = 0
175    for loc in locations:
176        i = LOCATIONS[index]
177        var_values = [(assign_vars[(i, j)], 1) for j in loc[0]]
178        var_values.append((use_vars[i], 1))
179        var_values.append((waste_vars[i], loc[1]))
180        dv = dippy.DecompVar(var_values, None, loc[1])
181        bvs.append((i, dv))
182        index += 1
183    return bvs
184
185def one_each(prob):
186   bvs = []
187   for index, loc in enumerate(LOCATIONS):
188       lc = [PRODUCTS[index]]
189       waste = CAPACITY - REQUIREMENT[PRODUCTS[index]]
190       var_values = [(assign_vars[(loc, j)], 1) for j in lc]
191       var_values.append((use_vars[loc], 1))
192       var_values.append((waste_vars[loc], waste))
193
194       dv = dippy.DecompVar(var_values, None, waste)
195       bvs.append((loc, dv))
196   return bvs
197
198prob.init_vars = first_fit
199##prob.init_vars = one_each
200
201prob.writeLP('facility_main.lp')
202for n, i in enumerate(LOCATIONS):
203    prob.writeRelaxed(n, 'facility_relax%s.lp' % i);
204
205dippy.Solve(prob, {
206    'TolZero': '%s' % tol,
207    'doPriceCut': '1',
208    'generateInitVars': '1', })
209
210# print solution
211for i in LOCATIONS:
212    if use_vars[i].varValue > tol:
213        print "Location ", i, \
214              " produces ", \
215              [j for j in PRODUCTS
216               if assign_vars[(i, j)].varValue > tol]
Note: See TracBrowser for help on using the repository browser.