1"""Tests for Python-MIP""" 2from itertools import product 3import pytest 4import networkx as nx 5from mip import Model, xsum, OptimizationStatus, MAXIMIZE, BINARY, INTEGER 6from mip import ConstrsGenerator, CutPool, maximize, CBC, GUROBI, Column 7from os import environ 8import math 9 10TOL = 1e-4 11 12SOLVERS = [CBC] 13if "GUROBI_HOME" in environ: 14 SOLVERS += [GUROBI] 15 16 17@pytest.mark.parametrize("solver", SOLVERS) 18def test_column_generation(solver: str): 19 L = 250 # bar length 20 m = 4 # number of requests 21 w = [187, 119, 74, 90] # size of each item 22 b = [1, 2, 2, 1] # demand for each item 23 24 # creating master model 25 master = Model(solver_name=solver) 26 27 # creating an initial set of patterns which cut one item per bar 28 # to provide the restricted master problem with a feasible solution 29 lambdas = [master.add_var(obj=1, name="lambda_%d" % (j + 1)) for j in range(m)] 30 31 # creating constraints 32 constraints = [] 33 for i in range(m): 34 constraints.append(master.add_constr(lambdas[i] >= b[i], name="i_%d" % (i + 1))) 35 36 # creating the pricing problem 37 pricing = Model(solver_name=solver) 38 39 # creating pricing variables 40 a = [ 41 pricing.add_var(obj=0, var_type=INTEGER, name="a_%d" % (i + 1)) for i in range(m) 42 ] 43 44 # creating pricing constraint 45 pricing += xsum(w[i] * a[i] for i in range(m)) <= L, "bar_length" 46 47 new_vars = True 48 while new_vars: 49 ########## 50 # STEP 1: solving restricted master problem 51 ########## 52 master.optimize() 53 54 ########## 55 # STEP 2: updating pricing objective with dual values from master 56 ########## 57 pricing += 1 - xsum(constraints[i].pi * a[i] for i in range(m)) 58 59 # solving pricing problem 60 pricing.optimize() 61 62 ########## 63 # STEP 3: adding the new columns (if any is obtained with negative reduced cost) 64 ########## 65 # checking if columns with negative reduced cost were produced and 66 # adding them into the restricted master problem 67 if pricing.objective_value < -TOL: 68 pattern = [a[i].x for i in range(m)] 69 column = Column(constraints, pattern) 70 lambdas.append( 71 master.add_var( 72 obj=1, column=column, name="lambda_%d" % (len(lambdas) + 1) 73 ) 74 ) 75 76 # if no column with negative reduced cost was produced, then linear 77 # relaxation of the restricted master problem is solved 78 else: 79 new_vars = False 80 81 # printing the solution 82 assert len(lambdas) == 8 83 assert round(master.objective_value) == 3 84 85 86@pytest.mark.parametrize("solver", SOLVERS) 87def test_cutting_stock(solver: str): 88 n = 10 # maximum number of bars 89 L = 250 # bar length 90 m = 4 # number of requests 91 w = [187, 119, 74, 90] # size of each item 92 b = [1, 2, 2, 1] # demand for each item 93 94 # creating the model 95 model = Model(solver_name=solver) 96 x = { 97 (i, j): model.add_var(obj=0, var_type=INTEGER, name="x[%d,%d]" % (i, j)) 98 for i in range(m) 99 for j in range(n) 100 } 101 y = {j: model.add_var(obj=1, var_type=BINARY, name="y[%d]" % j) for j in range(n)} 102 103 # constraints 104 for i in range(m): 105 model.add_constr(xsum(x[i, j] for j in range(n)) >= b[i]) 106 for j in range(n): 107 model.add_constr(xsum(w[i] * x[i, j] for i in range(m)) <= L * y[j]) 108 109 # additional constraints to reduce symmetry 110 for j in range(1, n): 111 model.add_constr(y[j - 1] >= y[j]) 112 113 # optimizing the model 114 model.optimize() 115 116 # sanity tests 117 assert model.status == OptimizationStatus.OPTIMAL 118 assert abs(model.objective_value - 3) <= 1e-4 119 assert sum(x.x for x in model.vars) >= 5 120 121 122@pytest.mark.parametrize("solver", SOLVERS) 123def test_knapsack(solver: str): 124 p = [10, 13, 18, 31, 7, 15] 125 w = [11, 15, 20, 35, 10, 33] 126 c, I = 47, range(len(w)) 127 128 m = Model("knapsack", solver_name=solver) 129 130 x = [m.add_var(var_type=BINARY) for i in I] 131 132 m.objective = maximize(xsum(p[i] * x[i] for i in I)) 133 134 m += xsum(w[i] * x[i] for i in I) <= c, "cap" 135 136 m.optimize() 137 138 assert m.status == OptimizationStatus.OPTIMAL 139 assert round(m.objective_value) == 41 140 141 m.constr_by_name("cap").rhs = 60 142 m.optimize() 143 144 assert m.status == OptimizationStatus.OPTIMAL 145 assert round(m.objective_value) == 51 146 147 # modifying objective function 148 m.objective = m.objective + 10 * x[0] + 15 * x[1] 149 assert abs(m.objective.expr[x[0]] - 20) <= 1e-10 150 assert abs(m.objective.expr[x[1]] - 28) <= 1e-10 151 152 153@pytest.mark.parametrize("solver", SOLVERS) 154def test_queens(solver: str): 155 """MIP model n-queens""" 156 n = 50 157 queens = Model("queens", MAXIMIZE, solver_name=solver) 158 queens.verbose = 0 159 160 x = [ 161 [queens.add_var("x({},{})".format(i, j), var_type=BINARY) for j in range(n)] 162 for i in range(n) 163 ] 164 165 # one per row 166 for i in range(n): 167 queens += xsum(x[i][j] for j in range(n)) == 1, "row({})".format(i) 168 169 # one per column 170 for j in range(n): 171 queens += xsum(x[i][j] for i in range(n)) == 1, "col({})".format(j) 172 173 # diagonal \ 174 for p, k in enumerate(range(2 - n, n - 2 + 1)): 175 queens += ( 176 xsum(x[i][j] for i in range(n) for j in range(n) if i - j == k) <= 1, 177 "diag1({})".format(p), 178 ) 179 180 # diagonal / 181 for p, k in enumerate(range(3, n + n)): 182 queens += ( 183 xsum(x[i][j] for i in range(n) for j in range(n) if i + j == k) <= 1, 184 "diag2({})".format(p), 185 ) 186 187 queens.optimize() 188 assert queens.status == OptimizationStatus.OPTIMAL # "model status" 189 190 # querying problem variables and checking opt results 191 total_queens = 0 192 for v in queens.vars: 193 # basic integrality test 194 assert v.x <= TOL or v.x >= 1 - TOL 195 total_queens += v.x 196 197 # solution feasibility 198 rows_with_queens = 0 199 for i in range(n): 200 if abs(sum(x[i][j].x for j in range(n)) - 1) <= TOL: 201 rows_with_queens += 1 202 203 assert abs(total_queens - n) <= TOL # "feasible solution" 204 assert rows_with_queens == n # "feasible solution" 205 206 207@pytest.mark.parametrize("solver", SOLVERS) 208def test_tsp(solver: str): 209 """tsp related tests""" 210 N = ["a", "b", "c", "d", "e", "f", "g"] 211 n = len(N) 212 i0 = N[0] 213 214 A = { 215 ("a", "d"): 56, 216 ("d", "a"): 67, 217 ("a", "b"): 49, 218 ("b", "a"): 50, 219 ("d", "b"): 39, 220 ("b", "d"): 37, 221 ("c", "f"): 35, 222 ("f", "c"): 35, 223 ("g", "b"): 35, 224 ("b", "g"): 25, 225 ("a", "c"): 80, 226 ("c", "a"): 99, 227 ("e", "f"): 20, 228 ("f", "e"): 20, 229 ("g", "e"): 38, 230 ("e", "g"): 49, 231 ("g", "f"): 37, 232 ("f", "g"): 32, 233 ("b", "e"): 21, 234 ("e", "b"): 30, 235 ("a", "g"): 47, 236 ("g", "a"): 68, 237 ("d", "c"): 37, 238 ("c", "d"): 52, 239 ("d", "e"): 15, 240 ("e", "d"): 20, 241 } 242 243 # input and output arcs per node 244 Aout = {n: [a for a in A if a[0] == n] for n in N} 245 Ain = {n: [a for a in A if a[1] == n] for n in N} 246 m = Model(solver_name=solver) 247 m.verbose = 1 248 249 x = {a: m.add_var(name="x({},{})".format(a[0], a[1]), var_type=BINARY) for a in A} 250 251 m.objective = xsum(c * x[a] for a, c in A.items()) 252 253 for i in N: 254 m += xsum(x[a] for a in Aout[i]) == 1, "out({})".format(i) 255 m += xsum(x[a] for a in Ain[i]) == 1, "in({})".format(i) 256 257 # continuous variable to prevent subtours: each 258 # city will have a different "identifier" in the planned route 259 y = {i: m.add_var(name="y({})".format(i), lb=0.0) for i in N} 260 261 # subtour elimination 262 for (i, j) in A: 263 if i0 not in [i, j]: 264 m.add_constr(y[i] - (n + 1) * x[(i, j)] >= y[j] - n) 265 266 m.relax() 267 m.optimize() 268 269 assert m.status == OptimizationStatus.OPTIMAL # "lp model status" 270 assert abs(m.objective_value - 238.75) <= TOL # "lp model objective" 271 272 # setting all variables to integer now 273 for v in m.vars: 274 v.var_type = INTEGER 275 m.optimize() 276 277 assert m.status == OptimizationStatus.OPTIMAL # "mip model status" 278 assert abs(m.objective_value - 262) <= TOL # "mip model objective" 279 280 281class SubTourCutGenerator(ConstrsGenerator): 282 """Class to generate cutting planes for the TSP""" 283 284 def generate_constrs(self, model: Model): 285 G = nx.DiGraph() 286 r = [(v, v.x) for v in model.vars if v.name.startswith("x(")] 287 U = [v.name.split("(")[1].split(",")[0] for v, f in r] 288 V = [v.name.split(")")[0].split(",")[1] for v, f in r] 289 N = list(set(U + V)) 290 cp = CutPool() 291 for i in range(len(U)): 292 G.add_edge(U[i], V[i], capacity=r[i][1]) 293 for (u, v) in product(N, N): 294 if u == v: 295 continue 296 val, (S, NS) = nx.minimum_cut(G, u, v) 297 if val <= 0.99: 298 arcsInS = [ 299 (v, f) for i, (v, f) in enumerate(r) if U[i] in S and V[i] in S 300 ] 301 if sum(f for v, f in arcsInS) >= (len(S) - 1) + 1e-4: 302 cut = xsum(1.0 * v for v, fm in arcsInS) <= len(S) - 1 303 cp.add(cut) 304 if len(cp.cuts) > 256: 305 for cut in cp.cuts: 306 model.add_cut(cut) 307 return 308 for cut in cp.cuts: 309 model.add_cut(cut) 310 311 312@pytest.mark.parametrize("solver", SOLVERS) 313def test_tsp_cuts(solver: str): 314 """tsp related tests""" 315 N = ["a", "b", "c", "d", "e", "f", "g"] 316 n = len(N) 317 i0 = N[0] 318 319 A = { 320 ("a", "d"): 56, 321 ("d", "a"): 67, 322 ("a", "b"): 49, 323 ("b", "a"): 50, 324 ("d", "b"): 39, 325 ("b", "d"): 37, 326 ("c", "f"): 35, 327 ("f", "c"): 35, 328 ("g", "b"): 35, 329 ("b", "g"): 25, 330 ("a", "c"): 80, 331 ("c", "a"): 99, 332 ("e", "f"): 20, 333 ("f", "e"): 20, 334 ("g", "e"): 38, 335 ("e", "g"): 49, 336 ("g", "f"): 37, 337 ("f", "g"): 32, 338 ("b", "e"): 21, 339 ("e", "b"): 30, 340 ("a", "g"): 47, 341 ("g", "a"): 68, 342 ("d", "c"): 37, 343 ("c", "d"): 52, 344 ("d", "e"): 15, 345 ("e", "d"): 20, 346 } 347 348 # input and output arcs per node 349 Aout = {n: [a for a in A if a[0] == n] for n in N} 350 Ain = {n: [a for a in A if a[1] == n] for n in N} 351 m = Model(solver_name=solver) 352 m.verbose = 0 353 354 x = {a: m.add_var(name="x({},{})".format(a[0], a[1]), var_type=BINARY) for a in A} 355 356 m.objective = xsum(c * x[a] for a, c in A.items()) 357 358 for i in N: 359 m += xsum(x[a] for a in Aout[i]) == 1, "out({})".format(i) 360 m += xsum(x[a] for a in Ain[i]) == 1, "in({})".format(i) 361 362 # continuous variable to prevent subtours: each 363 # city will have a different "identifier" in the planned route 364 y = {i: m.add_var(name="y({})".format(i), lb=0.0) for i in N} 365 366 # subtour elimination 367 for (i, j) in A: 368 if i0 not in [i, j]: 369 m.add_constr(y[i] - (n + 1) * x[(i, j)] >= y[j] - n) 370 371 m.cuts_generator = SubTourCutGenerator() 372 373 # tiny model, should be enough to find the optimal 374 m.max_seconds = 10 375 m.max_nodes = 100 376 m.max_solutions = 1000 377 m.optimize() 378 379 assert m.status == OptimizationStatus.OPTIMAL # "mip model status" 380 assert abs(m.objective_value - 262) <= TOL # "mip model objective" 381 382 383@pytest.mark.parametrize("solver", SOLVERS) 384def test_tsp_mipstart(solver: str): 385 """tsp related tests""" 386 N = ["a", "b", "c", "d", "e", "f", "g"] 387 n = len(N) 388 i0 = N[0] 389 390 A = { 391 ("a", "d"): 56, 392 ("d", "a"): 67, 393 ("a", "b"): 49, 394 ("b", "a"): 50, 395 ("d", "b"): 39, 396 ("b", "d"): 37, 397 ("c", "f"): 35, 398 ("f", "c"): 35, 399 ("g", "b"): 35, 400 ("b", "g"): 25, 401 ("a", "c"): 80, 402 ("c", "a"): 99, 403 ("e", "f"): 20, 404 ("f", "e"): 20, 405 ("g", "e"): 38, 406 ("e", "g"): 49, 407 ("g", "f"): 37, 408 ("f", "g"): 32, 409 ("b", "e"): 21, 410 ("e", "b"): 30, 411 ("a", "g"): 47, 412 ("g", "a"): 68, 413 ("d", "c"): 37, 414 ("c", "d"): 52, 415 ("d", "e"): 15, 416 ("e", "d"): 20, 417 } 418 419 # input and output arcs per node 420 Aout = {n: [a for a in A if a[0] == n] for n in N} 421 Ain = {n: [a for a in A if a[1] == n] for n in N} 422 m = Model(solver_name=solver) 423 m.verbose = 0 424 425 x = {a: m.add_var(name="x({},{})".format(a[0], a[1]), var_type=BINARY) for a in A} 426 427 m.objective = xsum(c * x[a] for a, c in A.items()) 428 429 for i in N: 430 m += xsum(x[a] for a in Aout[i]) == 1, "out({})".format(i) 431 m += xsum(x[a] for a in Ain[i]) == 1, "in({})".format(i) 432 433 # continuous variable to prevent subtours: each 434 # city will have a different "identifier" in the planned route 435 y = {i: m.add_var(name="y({})".format(i), lb=0.0) for i in N} 436 437 # subtour elimination 438 for (i, j) in A: 439 if i0 not in [i, j]: 440 m.add_constr(y[i] - (n + 1) * x[(i, j)] >= y[j] - n) 441 442 route = ["a", "g", "f", "c", "d", "e", "b", "a"] 443 m.start = [(x[route[i - 1], route[i]], 1.0) for i in range(1, len(route))] 444 m.optimize() 445 446 assert m.status == OptimizationStatus.OPTIMAL 447 assert abs(m.objective_value - 262) <= TOL 448 449 450class TestAPI(object): 451 def build_model(self, solver): 452 """MIP model n-queens""" 453 n = 50 454 queens = Model("queens", MAXIMIZE, solver_name=solver) 455 queens.verbose = 0 456 457 x = [ 458 [queens.add_var("x({},{})".format(i, j), var_type=BINARY) for j in range(n)] 459 for i in range(n) 460 ] 461 462 # one per row 463 for i in range(n): 464 queens += xsum(x[i][j] for j in range(n)) == 1, "row{}".format(i) 465 466 # one per column 467 for j in range(n): 468 queens += xsum(x[i][j] for i in range(n)) == 1, "col{}".format(j) 469 470 # diagonal \ 471 for p, k in enumerate(range(2 - n, n - 2 + 1)): 472 queens += ( 473 xsum(x[i][j] for i in range(n) for j in range(n) if i - j == k) <= 1, 474 "diag1({})".format(p), 475 ) 476 477 # diagonal / 478 for p, k in enumerate(range(3, n + n)): 479 queens += ( 480 xsum(x[i][j] for i in range(n) for j in range(n) if i + j == k) <= 1, 481 "diag2({})".format(p), 482 ) 483 484 return n, queens 485 486 @pytest.mark.parametrize("solver", SOLVERS) 487 def test_constr_get_index(self, solver): 488 _, model = self.build_model(solver) 489 490 idx = model.solver.constr_get_index("row0") 491 assert idx >= 0 492 493 idx = model.solver.constr_get_index("col0") 494 assert idx >= 0 495 496 @pytest.mark.parametrize("solver", SOLVERS) 497 def test_remove_constrs(self, solver): 498 _, model = self.build_model(solver) 499 500 idx1 = model.solver.constr_get_index("row0") 501 assert idx1 >= 0 502 503 idx2 = model.solver.constr_get_index("col0") 504 assert idx2 >= 0 505 506 model.solver.remove_constrs([idx1, idx2]) 507 508 @pytest.mark.parametrize("solver", SOLVERS) 509 def test_constr_get_rhs(self, solver): 510 n, model = self.build_model(solver) 511 512 # test RHS of rows 513 for i in range(n): 514 idx1 = model.solver.constr_get_index("row{}".format(i)) 515 assert idx1 >= 0 516 assert model.solver.constr_get_rhs(idx1) == 1 517 518 # test RHS of columns 519 for i in range(n): 520 idx1 = model.solver.constr_get_index("col{}".format(i)) 521 assert idx1 >= 0 522 assert model.solver.constr_get_rhs(idx1) == 1 523 524 @pytest.mark.parametrize("solver", SOLVERS) 525 def test_constr_set_rhs(self, solver): 526 n, model = self.build_model(solver) 527 528 idx1 = model.solver.constr_get_index("row0") 529 assert idx1 >= 0 530 531 val = 10 532 model.solver.constr_set_rhs(idx1, val) 533 assert model.solver.constr_get_rhs(idx1) == val 534 535 @pytest.mark.parametrize("solver", SOLVERS) 536 def test_constr_by_name_rhs(self, solver): 537 n, model = self.build_model(solver) 538 539 val = 10 540 model.constr_by_name("row0").rhs = val 541 assert model.constr_by_name("row0").rhs == val 542 543 @pytest.mark.parametrize("solver", SOLVERS) 544 def test_var_by_name_rhs(self, solver): 545 n, model = self.build_model(solver) 546 547 v = model.var_by_name("x({},{})".format(0, 0)) 548 assert v is not None 549 550 551@pytest.mark.parametrize("val", range(1, 4)) 552@pytest.mark.parametrize("solver", SOLVERS) 553def test_variable_bounds(solver: str, val: int): 554 m = Model("bounds", solver_name=solver) 555 556 x = m.add_var(var_type=INTEGER, lb=0, ub=2 * val) 557 y = m.add_var(var_type=INTEGER, lb=val, ub=2 * val) 558 m.objective = maximize(x - y) 559 m.optimize() 560 assert m.status == OptimizationStatus.OPTIMAL 561 assert round(m.objective_value) == val 562 assert round(x.x) == 2 * val 563 assert round(y.x) == val 564 565 566@pytest.mark.parametrize("val", range(1, 4)) 567@pytest.mark.parametrize("solver", SOLVERS) 568def test_linexpr_x(solver: str, val: int): 569 m = Model("bounds", solver_name=solver) 570 571 x = m.add_var(lb=0, ub=2 * val) 572 y = m.add_var(lb=val, ub=2 * val) 573 obj = x - y 574 575 assert obj.x is None # No solution yet. 576 577 m.objective = maximize(obj) 578 m.optimize() 579 580 assert m.status == OptimizationStatus.OPTIMAL 581 assert round(m.objective_value) == val 582 assert round(x.x) == 2 * val 583 assert round(y.x) == val 584 585 # Check that the linear expression value is equal to the same expression 586 # calculated from the values of the variables. 587 assert abs((x + y).x - (x.x + y.x)) < TOL 588 assert abs((x + 2 * y).x - (x.x + 2 * y.x)) < TOL 589 assert abs((x + 2 * y + x).x - (x.x + 2 * y.x + x.x)) < TOL 590 assert abs((x + 2 * y + x + 1).x - (x.x + 2 * y.x + x.x + 1)) < TOL 591 assert abs((x + 2 * y + x + 1 + x / 2).x - (x.x + 2 * y.x + x.x + 1 + x.x / 2)) < TOL 592 593 594@pytest.mark.parametrize("solver", SOLVERS) 595def test_add_column(solver: str): 596 """Simple test which add columns in a specific way""" 597 m = Model() 598 x = m.add_var() 599 600 example_constr1 = m.add_constr(x >= 1, "constr1") 601 example_constr2 = m.add_constr(x <= 2, "constr2") 602 column1 = Column() 603 column1.constrs = [example_constr1] 604 column1.coeffs = [1] 605 second_var = m.add_var("second", column=column1) 606 column2 = Column() 607 column2.constrs = [example_constr2] 608 column2.coeffs = [2] 609 m.add_var("third", column=column2) 610 611 vthird = m.vars["third"] 612 assert vthird is not None 613 assert len(vthird.column.coeffs) == len(vthird.column.constrs) 614 assert len(vthird.column.coeffs) == 1 615 616 pconstr2 = m.constrs["constr2"] 617 assert vthird.column.constrs[0].name == pconstr2.name 618 assert len(example_constr1.expr.expr) == 2 619 assert second_var in example_constr1.expr.expr 620 assert x in example_constr1.expr.expr 621 622 623@pytest.mark.parametrize("val", range(1, 4)) 624@pytest.mark.parametrize("solver", SOLVERS) 625def test_float(solver: str, val: int): 626 m = Model("bounds", solver_name=solver) 627 x = m.add_var(lb=0, ub=2 * val) 628 y = m.add_var(lb=val, ub=2 * val) 629 obj = x - y 630 # No solution yet. __float__ MUST return a float type, so it returns nan. 631 assert obj.x is None 632 assert math.isnan(float(obj)) 633 m.objective = maximize(obj) 634 m.optimize() 635 assert m.status == OptimizationStatus.OPTIMAL 636 # test vars. 637 assert x.x == float(x) 638 assert y.x == float(y) 639 # test linear expressions. 640 assert float(x + y) == (x + y).x 641