|
| 1 | +""" |
| 2 | +## GAMSSOURCE: https://gams.com/latest/gamslib_ml/libhtml/gamslib_solmpool.html |
| 3 | +## LICENSETYPE: Demo |
| 4 | +## MODELTYPE: MIP |
| 5 | +## KEYWORDS: mixed integer linear programming, facility location problem, CPLEX, solution pool |
| 6 | +
|
| 7 | +A simple version of a facility location problem is used to show how the |
| 8 | +solution pool and the tools associated with it work. This example is taken |
| 9 | +from the Cplex 11 User's Manual (ILOG, Cplex 11 User's Manual, 2007) |
| 10 | +
|
| 11 | +A company is considering opening as many as four warehouses in order to serve |
| 12 | +nine different regions. The goal is to minimize the sum of fixed costs |
| 13 | +associated with opening warehouses as well as the various transportation |
| 14 | +costs incurred to ship goods from the warehouses to the regions. |
| 15 | +
|
| 16 | +Whether or not to open a warehouse is represented by binary variable ow. |
| 17 | +Whether or not to ship goods from warehouse i to region j is represented |
| 18 | +by binary variable oa. |
| 19 | +
|
| 20 | +Each region needs a specified amount of goods, and each warehouse can store |
| 21 | +only a limited quantity of goods. In addition, each region must be served |
| 22 | +by exactly one warehouse. |
| 23 | +
|
| 24 | +The following GAMSPy program demonstrates a number of different |
| 25 | +approaches to collecting solution pools. GAMSPy will store the solutions |
| 26 | +in a merged GDX files which can then be further used by other programs |
| 27 | +or the same GAMSPy run. GAMS/Cplex will have the variables with an extra |
| 28 | +index as parameters in the merged solution file. |
| 29 | +""" |
| 30 | + |
| 31 | +import numpy as np |
| 32 | + |
| 33 | +import gamspy as gp |
| 34 | + |
| 35 | + |
| 36 | +def load_solution(m: gp.Container) -> None: |
| 37 | + m.loadRecordsFromGdx( |
| 38 | + "solnpool.gdx", |
| 39 | + symbol_names={ |
| 40 | + "index": "solnpool", |
| 41 | + "oa": "oaX", |
| 42 | + "ow": "owX", |
| 43 | + "totcost": "totcostX", |
| 44 | + "tcost": "tcostX", |
| 45 | + "fcost": "fcostX", |
| 46 | + }, |
| 47 | + ) |
| 48 | + |
| 49 | + solnpool = m["solnpool"] |
| 50 | + num_solutions = len(solnpool) |
| 51 | + print(f"{num_solutions=}") |
| 52 | + |
| 53 | + xcostX = m["xcostX"] |
| 54 | + soln = m["soln"] |
| 55 | + u = m["u"] |
| 56 | + totcostX = m["totcostX"] |
| 57 | + tcostX = m["tcostX"] |
| 58 | + fcostX = m["fcostX"] |
| 59 | + |
| 60 | + xcostX[soln, u] = 0 |
| 61 | + xcostX[solnpool, "totcost"] = totcostX[solnpool] |
| 62 | + xcostX[solnpool, "tcost"] = tcostX[solnpool] |
| 63 | + xcostX[solnpool, "fcost"] = fcostX[solnpool] |
| 64 | + xcostX[solnpool, "fcost^0.96"] = fcostX[solnpool] ** 0.96 |
| 65 | + |
| 66 | + |
| 67 | +def calculate_diff(oaX: gp.Parameter) -> int: |
| 68 | + aggdiff = 0 |
| 69 | + oax_records = oaX.records.to_numpy() |
| 70 | + for row1 in oax_records: |
| 71 | + for row2 in oax_records: |
| 72 | + if row1[0] == row2[0] or row1[1] != row2[1] or row1[2] != row2[2]: |
| 73 | + continue |
| 74 | + |
| 75 | + aggdiff += np.logical_xor(row1[3], row2[3]) |
| 76 | + |
| 77 | + return aggdiff |
| 78 | + |
| 79 | + |
| 80 | +def main(): |
| 81 | + m = gp.Container() |
| 82 | + i = gp.Set(m, records=[f"w{idx}" for idx in range(1, 5)]) |
| 83 | + j = gp.Set(m, records=[f"r{idx}" for idx in range(1, 10)]) |
| 84 | + f = gp.Parameter(m, domain=i, records=np.array([130, 150, 170, 180])) |
| 85 | + c = gp.Parameter(m, domain=i, records=np.array([90, 110, 130, 150])) |
| 86 | + d = gp.Parameter( |
| 87 | + m, domain=j, records=np.array([10, 10, 12, 15, 15, 15, 20, 20, 30]) |
| 88 | + ) |
| 89 | + t = gp.Parameter( |
| 90 | + m, |
| 91 | + domain=[j, i], |
| 92 | + records=np.array( |
| 93 | + [ |
| 94 | + [10, 30, 25, 55], |
| 95 | + [10, 25, 25, 45], |
| 96 | + [20, 23, 30, 40], |
| 97 | + [25, 10, 26, 40], |
| 98 | + [28, 12, 20, 29], |
| 99 | + [36, 19, 16, 22], |
| 100 | + [40, 39, 22, 27], |
| 101 | + [75, 65, 55, 35], |
| 102 | + [34, 43, 41, 62], |
| 103 | + ] |
| 104 | + ), |
| 105 | + ) |
| 106 | + totcost = gp.Variable(m) |
| 107 | + fcost = gp.Variable(m) |
| 108 | + tcost = gp.Variable(m) |
| 109 | + ow = gp.Variable(m, domain=i, type=gp.VariableType.BINARY) |
| 110 | + oa = gp.Variable(m, domain=[i, j], type=gp.VariableType.BINARY) |
| 111 | + _ = gp.Equation(m, name="deftotcost", definition=totcost == fcost + tcost) |
| 112 | + _ = gp.Equation(m, name="deffcost", definition=fcost == gp.Sum(i, f[i] * ow[i])) |
| 113 | + _ = gp.Equation( |
| 114 | + m, name="deftcost", definition=tcost == gp.Sum((i, j), t[j, i] * oa[i, j]) |
| 115 | + ) |
| 116 | + _ = gp.Equation( |
| 117 | + m, name="defwcap", domain=i, definition=gp.Sum(j, d[j] * oa[i, j]) <= c[i] |
| 118 | + ) |
| 119 | + _ = gp.Equation(m, name="onew", domain=j, definition=gp.Sum(i, oa[i, j]) == 1) |
| 120 | + _ = gp.Equation(m, name="defow", domain=[i, j], definition=ow[i] >= oa[i, j]) |
| 121 | + loc = gp.Model(m, equations=m.getEquations(), sense=gp.Sense.MIN, objective=totcost) |
| 122 | + |
| 123 | + soln = gp.Set(m, records=[f"soln_loc_p{idx}" for idx in range(1, 1001)]) |
| 124 | + _ = gp.Set(m, name="solnpool", domain=soln) |
| 125 | + _ = gp.UniverseAlias(m, name="u") |
| 126 | + |
| 127 | + owX = gp.Parameter(m, domain=[soln, i]) |
| 128 | + oaX = gp.Parameter(m, domain=[soln, i, j]) |
| 129 | + _ = gp.Parameter(m, name="totcostX", domain=soln) |
| 130 | + _ = gp.Parameter(m, name="tcostX", domain=soln) |
| 131 | + _ = gp.Parameter(m, name="fcostX", domain=soln) |
| 132 | + xcostX = gp.Parameter(m, domain=[soln, "*"]) |
| 133 | + |
| 134 | + # 1. Collect the incumbents found during the regular optimize procedure |
| 135 | + # The Cplex option 'solnpool' triggers the collection of solutions in |
| 136 | + # the GDX container solnpool. |
| 137 | + loc.solve( |
| 138 | + solver="cplex", |
| 139 | + options=gp.Options(relative_optimality_gap=0), |
| 140 | + solver_options={"solnpoolmerge": "solnpool.gdx"}, |
| 141 | + ) |
| 142 | + |
| 143 | + load_solution(m) |
| 144 | + |
| 145 | + # 2. Use the populate procedure instead of regular optimize procedure (option |
| 146 | + # 'solnpoolpop 2'). By default we will generate 20 solutions determined by |
| 147 | + # the default of option populatelim (only valid for one thread). This is a |
| 148 | + # simple model which is quickly solved with heuristics and cuts, so we need |
| 149 | + # to let Cplex retain sufficient exploration space to find alternative solutions. |
| 150 | + # This is done with option 'solnpoolintensity 4'. With solutions where the optimal |
| 151 | + # solution can not so quickly be found, the default for this option should be suffict. |
| 152 | + loc.solve( |
| 153 | + solver="cplex", |
| 154 | + options=gp.Options(relative_optimality_gap=0), |
| 155 | + solver_options={ |
| 156 | + "solnpoolmerge": "solnpool.gdx", |
| 157 | + "solnpoolintensity": 4, |
| 158 | + "solnpoolpop": 2, |
| 159 | + }, |
| 160 | + ) |
| 161 | + load_solution(m) |
| 162 | + print(xcostX.records) |
| 163 | + |
| 164 | + # 3. Lets look at the diversity of the solution by counting the differences |
| 165 | + # between the shipment indicator variables. Lets limit the number of |
| 166 | + # solutions in the pool by 10 and require solution within 10% of the |
| 167 | + # optimum. |
| 168 | + loc.solve( |
| 169 | + solver="cplex", |
| 170 | + options=gp.Options(relative_optimality_gap=0), |
| 171 | + solver_options={ |
| 172 | + "solnpoolmerge": "solnpool.gdx", |
| 173 | + "threads": 1, |
| 174 | + "solnpoolintensity": 4, |
| 175 | + "solnpoolpop": 2, |
| 176 | + "solnpoolcapacity": 10, |
| 177 | + "solnpoolgap": 0.1, |
| 178 | + }, |
| 179 | + ) |
| 180 | + load_solution(m) |
| 181 | + |
| 182 | + aggdiff = calculate_diff(oaX) |
| 183 | + |
| 184 | + # 4. We repeat the experiment by now setting the solution pool replacement |
| 185 | + # strategy to 'diversity' and let the populate procedure find many more |
| 186 | + # solutions, we should see an increase in the aggregated difference |
| 187 | + loc.solve( |
| 188 | + solver="cplex", |
| 189 | + options=gp.Options(relative_optimality_gap=0), |
| 190 | + solver_options={ |
| 191 | + "solnpoolmerge": "solnpool.gdx", |
| 192 | + "threads": 1, |
| 193 | + "solnpoolintensity": 4, |
| 194 | + "solnpoolpop": 2, |
| 195 | + "solnpoolcapacity": 10, |
| 196 | + "solnpoolgap": 0.1, |
| 197 | + "populatelim": 10000, |
| 198 | + "solnpoolreplace": 2, |
| 199 | + }, |
| 200 | + ) |
| 201 | + load_solution(m) |
| 202 | + |
| 203 | + aggdiffX = calculate_diff(oaX) |
| 204 | + assert aggdiffX > aggdiff |
| 205 | + |
| 206 | + # 5. We can fine tune diversity by using a diversity filter. Suppose that |
| 207 | + # facilities w1 and w2 are open. Let a solution keeping those two facilities |
| 208 | + # open be the reference. We use a diversity filter to stipulate that any |
| 209 | + # solution added to the solution pool must differ from the reference by |
| 210 | + # decisions to open or close at least two other facilities. The following |
| 211 | + # filter enforces this diversity by specifying a minimum diversity of 2. |
| 212 | + # Note that the reference solution becomes the incumbent and is reported as |
| 213 | + # the first solution in the pool. |
| 214 | + loc.solve( |
| 215 | + solver="cplex", |
| 216 | + options=gp.Options(relative_optimality_gap=0), |
| 217 | + solver_options={ |
| 218 | + "solnpoolmerge": "solnpool.gdx", |
| 219 | + "threads": 1, |
| 220 | + "solnpoolintensity": 4, |
| 221 | + "solnpoolpop": 2, |
| 222 | + "divfltlo": 2, |
| 223 | + "writeflt": "ow2.flt", |
| 224 | + "ow.divflt('w1')": 1, |
| 225 | + "ow.divflt('w2')": 1, |
| 226 | + "ow.divflt('w3')": 0, |
| 227 | + "ow.divflt('w4')": 0, |
| 228 | + }, |
| 229 | + ) |
| 230 | + load_solution(m) |
| 231 | + print(owX.records) |
| 232 | + |
| 233 | + # 6. Solutions ordered best to worst. |
| 234 | + loc.solve( |
| 235 | + solver="cplex", |
| 236 | + options=gp.Options(relative_optimality_gap=0), |
| 237 | + solver_options={ |
| 238 | + "solnpoolmerge": "solnpool.gdx", |
| 239 | + "threads": 1, |
| 240 | + "solnpoolintensity": 4, |
| 241 | + "solnpoolpop": 2, |
| 242 | + }, |
| 243 | + ) |
| 244 | + print(xcostX[soln, "totcost"].records) |
| 245 | + |
| 246 | + |
| 247 | +if __name__ == "__main__": |
| 248 | + main() |
0 commit comments