Binpacking -- Multiple Constraints: Weight+volume

I have a dataset with 50,000 orders. Each order has ~20 products. Product volume and weight are present (as well as x,y,z dimensions). I have shipping boxes of constant volume V_ma

Solution 1:

Here is a quick prototype using cvxpy (1.0 branch!) and CoinOR's Cbc MIP-solver through cylp. (everything is open-source)

I'm using cvxpy as it allows beautiful concise modelling (at some cost as cvxpy does more than modelling!). In a real-world implementation, one would feed those solvers directly (less nice code) which will also improve performance (no time taken by cvxpy; and we can make use of Simplex-features like dedicated variable-bounds). This also allows tuning the solver for your problem (e.g. cutting-plane setup of Cbc/Cgl). You also also use time-limits or MIPgaps then to get good approximations if your instances are too hard (NP-hardness!).

The first approach of improving performance (using cvxpy; no solver-options in this version) would be some kind of symmetry-reduction (use the first N boxes; don't scramble those N << M boxes around). Edit: most simple approach added -> see below!

As you seem to got one unlimited supply of equal boxes, optimizing orders are independent! This fact is used and the code tackles the problem of optimizing one single order! (this would change if you got different boxes and cardinalities and using some box for some order disallow it's usage in other orders). Independent-solving follows the theory. In practice, when the outer-language is python, there might be some merit in doing one big solve over all orders at once (solver will somewhat recognize independence; but it's hard to say if that's something to try).

This code:

• is quite minimal
• is (imho) in a nice mathematical-form
• should scale well for larger and more complex examples
  • (given this high-quality solver; the default one which is shipped will struggle early)
• will find a global-optimum in finite time (complete)

(Install might be troublesome on non-Linux systems; In this case: take this as an approach instead of ready-to-use code)

Code

import numpy as np
import cvxpy as cvx
from timeit import default_timer as time

# Data
order_vols = [8, 4, 12, 18, 5, 2, 1, 4]
order_weights = [5, 3, 2, 5, 3, 4, 5, 6]

box_vol = 20
box_weight = 12

N_ITEMS = len(order_vols)
max_n_boxes = len(order_vols) # real-world: heuristic?""" Optimization """
M = N_ITEMS + 1# VARIABLES
box_used = cvx.Variable(max_n_boxes, boolean=True)
box_vol_content = cvx.Variable(max_n_boxes)
box_weight_content = cvx.Variable(max_n_boxes)
box_item_map = cvx.Variable((max_n_boxes, N_ITEMS), boolean=True)

# CONSTRAINTS
cons = []

# each item is shipped once
cons.append(cvx.sum(box_item_map, axis=0) == 1)

# box is used when >=1 item is using it
cons.append(box_used * M >= cvx.sum(box_item_map, axis=1))

# box vol constraints
cons.append(box_item_map * order_vols <= box_vol)

# box weight constraints
cons.append(box_item_map * order_weights <= box_weight)

problem = cvx.Problem(cvx.Minimize(cvx.sum(box_used)), cons)
start_t = time()
problem.solve(solver='CBC', verbose=True)
end_t = time()
print('time used (cvxpys reductions & solving): ', end_t - start_t)
print(problem.status)
print(problem.value)
print(box_item_map.value)

""" Reconstruct solution """
n_boxes_used = int(np.round(problem.value))
box_inds_used = np.where(np.isclose(box_used.value, 1.0))[0]

print('N_BOXES USED: ', n_boxes_used)
for box inrange(n_boxes_used):
    print('Box ', box)
    raw = box_item_map[box_inds_used[box]]
    items = np.where(np.isclose(raw.value, 1.0))[0]
    vol_used = 0
    weight_used = 0for item in items:
        print('   item ', item)
        print('       vol: ', order_vols[item])
        print('       weight: ', order_weights[item])
        vol_used += order_vols[item]
        weight_used += order_weights[item]
    print(' total vol: ', vol_used)
    print(' total weight: ', weight_used)

Output

WelcometotheCBCMILPSolverVersion:2.9.9Build Date:Jan152018commandline-ICbcModel-solve-quit(defaultstrategy1)Continuousobjectivevalueis0.888889-0.00secondsCgl0006I8SOS(64membersoutof72)with0overlaps-toomuchoverlaportoomanyothersCgl0009I8elementschangedCgl0003I0fixed,0tightenedbounds,19strengthenedrows,0substitutionsCgl0004Iprocessedmodelhas32rows,72columns(72integer(72ofwhichbinary))and280elementsCutoffincrementincreasedfrom1e-05to0.9999Cbc0038IInitialstate-9integersunsatisfiedsum-2.75909Cbc0038I Pass   1:suminf.1.60000(5)obj.3iterations10Cbc0038I Pass   2:suminf.0.98824(5)obj.3iterations5Cbc0038I Pass   3:suminf.0.90889(5)obj.3.02iterations12Cbc0038I Pass   4:suminf.0.84444(3)obj.4iterations8Cbc0038ISolutionfoundof4Cbc0038IBeforeminibranchandbound,60integersatboundfixedand0continuousCbc0038IFullproblem32rows72columns,reducedto0rows0columnsCbc0038IMinibranchandbounddidnotimprovesolution(0.01seconds)Cbc0038IRoundagainwithcutoffof2.97509Cbc0038I Pass   5:suminf.1.62491(7)obj.2.97509iterations2Cbc0038I Pass   6:suminf.1.67224(8)obj.2.97509iterations7Cbc0038I Pass   7:suminf.1.24713(5)obj.2.97509iterations3Cbc0038I Pass   8:suminf.1.77491(5)obj.2.97509iterations9Cbc0038I Pass   9:suminf.1.08405(6)obj.2.97509iterations8Cbc0038I Pass  10:suminf.1.57481(7)obj.2.97509iterations12Cbc0038I Pass  11:suminf.1.15815(6)obj.2.97509iterations1Cbc0038I Pass  12:suminf.1.10425(7)obj.2.97509iterations17Cbc0038I Pass  13:suminf.1.05568(8)obj.2.97509iterations17Cbc0038I Pass  14:suminf.0.45188(6)obj.2.97509iterations15Cbc0038I Pass  15:suminf.1.67468(8)obj.2.97509iterations22Cbc0038I Pass  16:suminf.1.42023(8)obj.2.97509iterations2Cbc0038I Pass  17:suminf.1.92437(7)obj.2.97509iterations15Cbc0038I Pass  18:suminf.1.82742(7)obj.2.97509iterations8Cbc0038I Pass  19:suminf.1.31741(10)obj.2.97509iterations15Cbc0038I Pass  20:suminf.1.01947(6)obj.2.97509iterations12Cbc0038I Pass  21:suminf.1.57481(7)obj.2.97509iterations14Cbc0038I Pass  22:suminf.1.15815(6)obj.2.97509iterations1Cbc0038I Pass  23:suminf.1.10425(7)obj.2.97509iterations15Cbc0038I Pass  24:suminf.1.08405(6)obj.2.97509iterations1Cbc0038I Pass  25:suminf.3.06344(10)obj.2.97509iterations13Cbc0038I Pass  26:suminf.2.57488(8)obj.2.97509iterations10Cbc0038I Pass  27:suminf.2.43925(7)obj.2.97509iterations1Cbc0038I Pass  28:suminf.0.91380(3)obj.2.97509iterations6Cbc0038I Pass  29:suminf.0.46935(3)obj.2.97509iterations6Cbc0038I Pass  30:suminf.0.46935(3)obj.2.97509iterations0Cbc0038I Pass  31:suminf.0.91380(3)obj.2.97509iterations8Cbc0038I Pass  32:suminf.1.96865(12)obj.2.97509iterations23Cbc0038I Pass  33:suminf.1.40385(6)obj.2.97509iterations13Cbc0038I Pass  34:suminf.1.90833(7)obj.2.79621iterations16Cbc0038INosolutionfoundthismajorpassCbc0038IBeforeminibranchandbound,42integersatboundfixedand0continuousCbc0038IFullproblem32rows72columns,reducedto20rows27columnsCbc0038IMinibranchandboundimprovedsolutionfrom4to3(0.06seconds)Cbc0038IAfter0.06seconds-Feasibilitypumpexitingwithobjectiveof3-took0.06secondsCbc0012IIntegersolutionof3foundbyfeasibilitypumpafter0iterationsand0nodes(0.06seconds)Cbc0001ISearchcompleted-bestobjective3,took0iterationsand0nodes(0.06seconds)Cbc0035IMaximumdepth0,0variablesfixedonreducedcostCutsatrootnodechangedobjectivefrom2.75to2.75Probingwastried0timesandcreated0cutsofwhich0wereactiveafteraddingroundsofcuts(0.000seconds)Gomorywastried0timesandcreated0cutsofwhich0wereactiveafteraddingroundsofcuts(0.000seconds)Knapsackwastried0timesandcreated0cutsofwhich0wereactiveafteraddingroundsofcuts(0.000seconds)Cliquewastried0timesandcreated0cutsofwhich0wereactiveafteraddingroundsofcuts(0.000seconds)MixedIntegerRounding2wastried0timesandcreated0cutsofwhich0wereactiveafteraddingroundsofcuts(0.000seconds)FlowCoverwastried0timesandcreated0cutsofwhich0wereactiveafteraddingroundsofcuts(0.000seconds)TwoMirCutswastried0timesandcreated0cutsofwhich0wereactiveafteraddingroundsofcuts(0.000seconds)Result-OptimalsolutionfoundObjective value:3.00000000Enumerated nodes:0Total iterations:0Time(CPUseconds):0.07Time(Wallclockseconds):0.05Totaltime(CPUseconds):0.07(Wallclockseconds):0.05

and the more interesting part:

timeused(cvxpysreductions&solving):0.07740794896380976optimal3.0
[[0.0.0.1.0.0.1.0.]
 [0.0.0.0.0.0.0.0.]
 [0.0.0.0.0.0.0.0.]
 [1.0.0.0.1.1.0.0.]
 [0.0.0.0.0.0.0.0.]
 [0.0.0.0.0.0.0.0.]
 [0.1.1.0.0.0.0.1.]
 [0.0.0.0.0.0.0.0.]]
N_BOXES USED:3Box0item3vol:18weight:5item6vol:1weight:5total vol:19total weight:10Box1item0vol:8weight:5item4vol:5weight:3item5vol:2weight:4total vol:15total weight:12Box2item1vol:4weight:3item2vol:12weight:2item7vol:4weight:6total vol:20total weight:11

An instance following your dimensions like:

order_vols = [8, 4, 12, 18, 5, 2, 1, 4, 6, 5, 3, 2, 5, 11, 17, 15, 14, 14, 12, 20]
order_weights = [5, 3, 2, 5, 3, 4, 5, 6, 3, 11, 3, 8, 12, 3, 1, 5, 3, 5, 6, 7]

box_vol = 20box_weight = 12

will be more work of course:

Result-OptimalsolutionfoundObjective value:11.00000000Enumerated nodes:0Total iterations:2581Time(CPUseconds):0.78Time(Wallclockseconds):0.72Totaltime(CPUseconds):0.78(Wallclockseconds):0.72N_BOXES USED:11

Edit: Symmetry-reduction

Playing around with different formulations really shows that it's sometimes hard to say a-priori what helps and what does not!

But a cheap small symmetry-exploiting change should always work (if the size of an order is not too big: 20 is ok; at 30 it probably starts to be critical). The approach is called lexicographic perturbation (Symmetry in Integer Linear Programming | François Margot).

We can add one variable-fixing (if there is a solution, there will always be a solution of same costs with this fixing):

cons.append(box_item_map[0,0] == 1)

and we change the objective:

# lex-perturbationc = np.power(2, np.arange(1, max_n_boxes+1))
problem = cvx.Problem(cvx.Minimize(c * box_used), cons)

# small change in reconstruction due to new objectiven_boxes_used = int(np.round(np.sum(box_used.value)))

For the above N=20 problem, we now achieve:

Result- Optimal solution found

Objective value:                4094.00000000
Enumerated nodes:               0
Total iterations:               474Time (CPU seconds):             0.60Time (Wallclock seconds):       0.44

Total time (CPU seconds):       0.60   (Wallclock seconds):       0.44time used (cvxpys reductions & solving):  0.46845901099732146

N_BOXES USED:  11