# source:projects/Dippy/use-api/examples/facility.py@584

Last change on this file since 584 was 584, checked in by mosu001, 2 years ago
File size: 8.2 KB
Line
1import coinor.pulp as pulp
2from pulp import *
3##import coinor.dippy as dippy
4import dippy
5
6from math import floor, ceil
7
8tol = pow(pow(2, -24), 2.0 / 3.0)
9
10from facility_ex1 import REQUIREMENT, PRODUCTS, LOCATIONS, CAPACITY
11##from facility_test1 import REQUIREMENT, PRODUCTS, LOCATIONS, CAPACITY
12##from facility_test2 import REQUIREMENT, PRODUCTS, LOCATIONS, CAPACITY
13##from facility_test6 import REQUIREMENT, PRODUCTS, LOCATIONS, CAPACITY
14
15prob = dippy.DipProblem("Facility_Location")
16
17assign_vars = LpVariable.dicts("AtLocation",
18              [(i, j) for i in LOCATIONS
19                      for j in PRODUCTS],
20              0, 1, LpBinary)
21use_vars    = LpVariable.dicts("UseLocation",
22              LOCATIONS, 0, 1, LpBinary)
23waste_vars  = LpVariable.dicts("Waste",
24              LOCATIONS, 0, CAPACITY)
25
26# objective: minimise waste
27prob += lpSum(waste_vars[i] for i in LOCATIONS), "min"
28
29# assignment constraints
30for j in PRODUCTS:
31    prob += lpSum(assign_vars[(i, j)] for i in LOCATIONS) == 1
32
33# aggregate CAPACITY constraints
34for i in LOCATIONS:
35    prob += lpSum(assign_vars[(i, j)] * REQUIREMENT[j]
36                  for j in PRODUCTS) + waste_vars[i] == \
37            CAPACITY * use_vars[i]
38
39# disaggregate CAPACITY constraints
40for i in LOCATIONS:
41    for j in PRODUCTS:
42        prob += assign_vars[(i, j)] <= use_vars[i]
43
44# Ordering constraints
45for index, location in enumerate(LOCATIONS):
46    if index > 0:
47        prob += use_vars[LOCATIONS[index-1]] >= use_vars[location]
48
49# Anti-symmetry branches
50def choose_antisymmetry_branch(prob, sol):
51    num_locations = sum(sol[use_vars[i]] for i in LOCATIONS)
52    up   = ceil(num_locations)  # Round up to next nearest integer
53    down = floor(num_locations) # Round down
54    if  (up - num_locations   > tol) \
55    and (num_locations - down > tol): # Is fractional?
56        # Down branch: provide upper bounds, lower bounds are default
57        down_branch_ub = dict([(use_vars[LOCATIONS[n]], 0)
58                               for n in range(int(down), len(LOCATIONS))])
59        # Up branch: provide lower bounds, upper bounds are default
60        up_branch_lb = dict([(use_vars[LOCATIONS[n]], 1)
61                             for n in range(0, int(up))])
62        # Return the advanced branch to DIP
63        return {}, down_branch_ub, up_branch_lb, {}
64
65prob.branch_method = choose_antisymmetry_branch
66
67def generate_weight_cuts(prob, sol):
68    print "In generate_weight_cuts, sol = ", sol
69
70    # Define mu and T for each knapsack
71    mu = {}
72    S = {}
73    for i in LOCATIONS:
74        mu[i] = CAPACITY
75        S[i] = []
76
77    # Use current assign_var values to assign items to locations
78    assigning = True
79    while assigning:
80        bestValue = 0
81        bestAssign = None
82        for i in LOCATIONS:
83            for j in PRODUCTS:
84                if j not in S[i]: # If this product is not in the subset
85                    if (sol[assign_vars[(i, j)]] > bestValue) \
86                    and (REQUIREMENT[j] <= mu[i]):
87                        # The assignment variable for this product is closer
88                        # to 1 than any other product checked, and "fits" in
89                        # this location's remaining space
90                        bestValue = sol[assign_vars[(i, j)]]
91                        bestAssign = (i, j)
92        # Make the best assignment found across all products and locactions
93        if bestAssign:
94            (i,j) = bestAssign
95            mu[i] -= REQUIREMENT[j] # Decrease spare CAPACITY at this location
96            S[i].append(j) # Assign this product to this location's set
97        else:
98            assigning = False # Didn't find anything to assign - stop
99
100    # Generate the weight cuts from the sets found above
101    new_cuts = []
102    for i in LOCATIONS:
103        if len(S[i]) > 0: # If an item assigned to this location
104            con = LpAffineExpression() # Start a new constraint
105            con += sum(REQUIREMENT[j] * assign_vars[(i, j)]
106                            for j in S[i])
107            con += sum(max(0, REQUIREMENT[j] - mu[i]) *
108                            assign_vars[(i, j)] for j in PRODUCTS
109                            if j not in S[i])
110            new_cuts.append(con <= CAPACITY - mu[i])
111
112    print new_cuts
113    # Return the set of cuts we created to DIP
114    return new_cuts
115
116prob.generate_cuts = generate_weight_cuts
117
118def first_fit_heuristic():
119    # Sort the items in descending weight order
120    productReqs = [(REQUIREMENT[j],j) for j in PRODUCTS]
121    productReqs.sort(reverse=True)
122
123    # Add items to locations, fitting in as much
124    # as possible at each location.
125    allLocations = []
126    while len(productReqs) > 0:
127        waste = CAPACITY
128        currentLocation = []
129        j = 0
130        while j < len(productReqs):
131            # Can we fit this product?
132            if productReqs[j][0] <= waste:
133                currentLocation.append(productReqs[j][1]) # index
134                waste -= productReqs[j][0] # requirement
135                productReqs.pop(j)
136            else:
137                # Try to fit next item
138                j += 1
139        allLocations.append((currentLocation, waste))
140
141    # Return a list of tuples: ([products],waste)
142    return allLocations
143
144def first_fit():
145    # Use generic first-fit heuristic that is shared
146    # between heuristics and initial variable generation
147    allLocations = first_fit_heuristic()
148
149    # Convert to decision variable values
150    sol = {}
151    for i in LOCATIONS:
152        for j in PRODUCTS:
153            sol[assign_vars[(i, j)]] = 0
154        sol[use_vars[i]] = 0
155        sol[waste_vars[i]] = 0
156
157    index = 0
158    for loc in allLocations:
159        i = LOCATIONS[index]
160        sol[use_vars[i]] = 1
161        sol[waste_vars[i]] = loc[1]
162        for j in loc[0]:
163            sol[assign_vars[(i, j)]] = 1
164        index += 1
165
166    return sol
167
168def frac_fit(xhat):
169    # Initialise solution
170    sol = {}
171    for i in LOCATIONS:
172        for j in PRODUCTS: sol[assign_vars[(i, j)]] = 0
173        sol[use_vars[i]] = 0
174        sol[waste_vars[i]] = 0
175
176    # Get the list of non-zero fractional assignments
177    fracAssigns = [ (xhat[assign_vars[(i, j)]], (i, j))
178                     for i in LOCATIONS for j in PRODUCTS
179                     if xhat[assign_vars[(i, j)]] > tol ]
180    fracAssigns.sort()
181
182    # Track which products and locations have been used
183    notAllocated = dict((j,True) for j in PRODUCTS)
184    notUsed      = dict((i,True) for i in LOCATIONS)
185    while len(fracAssigns) > 0:
186        fracX = fracAssigns.pop() # Get best frac. assignment left
187        (i,j) = fracX[1]
188        if notAllocated[j]:
189            if notUsed[i]: # Create a new location if needed
190                notUsed[i] = False
191                sol[use_vars[i]] = 1
192                sol[waste_vars[i]] = CAPACITY
193            if REQUIREMENT[j] <= sol[waste_vars[i]]: # Space left?
194                sol[assign_vars[(i, j)]] = 1
195                notAllocated[j] = False
196                sol[waste_vars[i]] -= REQUIREMENT[j]
197
198    # Allocate the remaining products
199    unallocated = [(REQUIREMENT[j],j) for j in PRODUCTS
200                                      if notAllocated[j]]
201    unallocated.sort(reverse=True)
202    unused = [i for i in LOCATIONS if notUsed[i]]
203    while len(unallocated) > 0:
204        waste = CAPACITY
205        index = 0
206        loc = unused.pop()
207        while index < len(unallocated):
208            (j_req, j) = unallocated[index]
209            if j_req <= waste:
210                unallocated.pop(index)
211                sol[assign_vars[(loc, j)]] = 1
212                waste -= j_req
213            else: index += 1
214        sol[use_vars[loc]] = 1
215        sol[waste_vars[loc]] = waste
216
217    return sol
218
219def heuristics(prob, xhat, cost):
220    sols = []
221    if prob.root_heuristic:
222        prob.root_heuristic = False # Don't run twice
223        sol = first_fit()
224        sols.append(sol)
225    if prob.node_heuristic:
226        sol = frac_fit(xhat)
227        sols.append(sol)
228
229    return sols
230
231prob.heuristics = heuristics
232prob.root_heuristic = True
233prob.node_heuristic = True
234
235prob.writeLP("facility.lp")
236
237dippy.Solve(prob, {
238    'TolZero': '%s' % tol,
239    'LogDebugLevel': 5,
240    'LogDumpModel': 5,
241})
242
243# print solution
244for i in LOCATIONS:
245    if use_vars[i].varValue > tol:
246        print "Location ", i, " produces ", \
247              [j for j in PRODUCTS
248               if assign_vars[(i, j)].varValue > tol]
Note: See TracBrowser for help on using the repository browser.