Internet-shopping optimization problem

Let us take a deep dive into the so-called "internet-shopping optimization problem", also known as "supplier selection and order allocation". This is a very common problem that occurs basically whenever a large number of items needs to be ordered from a large number of suppliers. As it involves integer quantities, it is usually formulated as a mixed-integer linear problem (MILP). Unfortunately, MILPs are known to be NP-hard - which more or less says that they can be very hard to solve.

You can replicate all results by running this notebook. Note that this will execute the open-source solver CP-SAT on your machine while Quicopt runs on our API server. This means the hardware is not identical, which it should be for a strictly rigorous comparison.

Even if you're running CP-SAT on CPUs faster than our test server, you should still notice that Quicopt is a lot faster! On top of that, CP-SAT uses 8 cores while Quicopt uses just one! (CP-SAT's number of cores can be controlled by adjusting WORKERS, see below. Finally, be mindful of the fact that CP-SAT operates non-deterministically, i.e. runtimes will fluctuate quite a lot, as will solution quality.)

The problem as a MILP

Variables

The variables of our optimization problems are given in the following table:

Objective function

The problem has four "levels": items to buy, variants of each item, offers for each variant, and finally the suppliers offering the offers. The cost function has two main terms:

This means that for every item $i$ and every variant $v$ we pay $C^{\text{choice}}_v$ if the variant is chosen, and also for every used offer $o$ of that variant we pay $C^{\text{pack}}_o$ per pack bought as well as the fixed cost $C^{\text{use}}_o$. The second sum is over the suppliers: the charge $C^{\text{fee}}_k$ if the supplier is used, less the shipping $C^{\text{ship}}_k$ if shipping is free.

Constraints

This problem has many constraints that restrict the decision variables meaningfully. For $\boldsymbol{z}$, $\boldsymbol{x}$ and $\boldsymbol{n}$, we have the following five:

For the supplier variables, we have two more:

CP-SAT

Here you can load the MILP from its zipped .mps file and solve it with CP-SAT. Note that CP-SAT's performance can be quite sensitive to the number of WORKERS as well as the time limit.

from ortools.linear_solver.python import model_builder as mb
import gzip

# define the model and load the instance data (21538 variables, 35578 constraint rows)
model = mb.Model()
model.import_from_mps_string(gzip.open(f"{DATA}/{INSTANCE}.mps.gz", "rt").read())

solver = mb.Solver("CP-SAT")
solver.set_time_limit_in_seconds(60)
WORKERS = 8 # number of CP-SAT "workers" (roughly equal to number of CPUs)
solver.set_solver_specific_parameters(f"num_workers:{WORKERS}")
status = solver.solve(model)

print(f"status    {status.name}")
print(f"cost      {solver.objective_value:,.0f}")
print(f"time      {solver.wall_time:.1f} s on {WORKERS} workers")

status    FEASIBLE
cost      162,129,469,300
time      60.1 s on 8 workers

Quicopt

Below we rebuild the model with Pyomo (to illustrate how it's done), and then send it to Quicopt's API service for solving.

import pyomo.environ as pyo
from quicopt import Client
from quicopt.structured import active, block, semiinteger

# define a Pyomo model
model = pyo.ConcreteModel()

Next we introduce the decision variables $\boldsymbol{z}, \boldsymbol{n}$ and $\boldsymbol{x}$ that we already encountered above. Note that constraint type (1) is already covered in the Pyomo declaration of $\boldsymbol{n}$:

# introduce the binary variables z
model.z = pyo.Var(V, within=pyo.Binary)

# integer numbers of packs ordered, constraint type (1) via bounds and semiinteger
model.n = pyo.Var(O, within=pyo.NonNegativeIntegers, bounds=lambda model, o: (N_min[o], N_max[o]))
semiinteger(model.n)

# shortcut aliases for n and z
z, n = model.z, model.n

# the binary variables x become indicators of the pack counts n
x = {o: active(n[o]) for o in O}

# tell the solver how the variables are related
for i in I:
  block(model, f"item_{i}", [n[o] for v in V_i[i] for o in O_v[v]] + [z[v] for v in V_i[i]])

Next, we define the remaining constraints, followed by the objective function:

# constraint type (2)
model.covered = pyo.Constraint(V, rule=lambda model, v: sum(Q_pack[o] * n[o] - Q_lost[o] * x[o] for o in O_v[v]) >= Q_need[v] * z[v])

# constraint type (3)
model.one_variant = pyo.Constraint(I, rule=lambda model, i: sum(z[v] for v in V_i[i]) == 1)

# constraint type (4)
model.exclusive = pyo.Constraint(E, rule=lambda model, e: sum(x[o] for o in O_e[e]) <= 1)

# constraint type (5)
model.stock = pyo.Constraint(R, rule=lambda model, r: sum(Q_pack[o] * n[o] for o in O_r[r]) <= Q_stock[r])

# the objective function, including the final two constraint types
pyo_obj_expr = (sum(C_choice[v] * z[v] for v in V)
              + sum(C_pack[o] * n[o] + C_use[o] * x[o] for o in O)
              + sum(pyo.Expr_if(IF=sum(x[o] for o in O_k[k]) >= 1,
                                THEN=C_fee[k],
                                ELSE=0) for k in K)  # constraint (6)
              + sum(pyo.Expr_if(IF=sum(S_pack[o] * n[o] + S_use[o] * x[o] for o in O_k[k]) >= S_free[k],
                                THEN=-C_ship[k],
                                ELSE=0) for k in K)) # constraint (7)

# give the objective to Pyomo
model.cost = pyo.Objective(sense=pyo.minimize, expr=pyo_obj_expr)

And then we send it off:

# solve with our API service
result = Client().solve(model, time_limit=2.0, gzip=True)
print(result.display)

▄▀▀▀▄        ▄                     █
█   █ █   █ ▄▄  ▄▀▀▀▄ ▄▀▀▀▄ █▀▀▀▄ ▀█▀▀
█ ▀▄▀ █   █  █  █   ▄ █   █ █▄▄▄▀  █  ▄
 ▀▀ ▀  ▀▀▀▀ ▀▀▀  ▀▀▀   ▀▀▀  █       ▀▀

├── shots
│   └── 1 · primal   162044219300   1.98s  ◀ best
├── status:     heuristic
├── feasible:   true
├── objective:  1.620442193e11
├── x:          x1=1, x10=0, x100=0, x1000=0, x10000=0, x10001=0, …  (10933 variables)
├── time_limit: 2 s
└── solve_time: 1.9852 s

    Intel Xeon Gold 6338N @ 2.20GHz · 16 cores, 1 core per shot

Let's inspect the solution quality:

print(f"CP-SAT's total cost                              {solver.objective_value:,.0f}")
print(f"Total cost from Quicopt's API service            {result.objective:,.0f}")

CP-SAT's total cost                              162,129,469,300
Total cost from Quicopt's API service            162,044,219,300

Quicopt found a better solution in less than two seconds! CP-SAT invested 60s on eight CPUs (which actually sometimes suffices to find the same solution as Quicopt). The notebook on our GitHub has additional proofs that the solution of Quicopt's solver is correct in the original MILP formulation.