商品のトラック積載率を最大化！Python×Pyomoで実現する数理最適化モデル

はじめまして、花王株式会社でデータサイエンティストをしている@m_jです。
本記事では、物流の効率化という現場課題に対して数理最適化を使ったアプローチを紹介します。特に「トラックへの商品の積載率」を最大化するために、制約充足問題（CSP）をどのように解くかについて解説します。

背景

花王では、商品を段ボール箱に入れて梱包した後、複数の段ボール箱を「パレット」と呼ばれる荷物を載せるための板の上に積み重ねます。そして、パレットごとトラックに積み込み、販売店へ運びます。販売店に到着した商品は、最終的にお客様のもとへ届けられます。

この物流の過程で重要なのがパレット上への段ボール箱の積み方です。段ボール箱を隙間なく配置できれば積載効率が向上し、物流や保管コストの削減、配送時の安定性確保につながります。逆に隙間が多いとトラックの積載効率が悪化し、コストや環境負荷にも影響します。

箱の積み方を決めるアルゴリズム

パレットへの箱の配置は、長方形詰め込み問題（Rectangle Packing Problem）として定式化できます。
「長方形詰め込み問題」とは、与えられた大きな長方形（コンテナやパレット）に、複数の小さな長方形（箱やアイテム）を重ならないように配置する組み合わせ問題です。花王で扱う箱の寸法は多岐に渡り、手作業での配置検討は困難なため、数理モデルによる自動配置が有効です。Pythonの最適化ライブラリ Pyomo を用いて数理最適化モデルを構築しました。

問題の定式化

前提条件として、パレットに積載する箱はすべて同じ種類です。パレットの寸法、箱の寸法、積載する箱の個数はあらかじめ与えられています。

目的関数

今回は制約充足問題とするため、目的関数は設定しません。

制約条件

1. 全ての箱がパレット内に収まる制約

トラックに収まるように、箱がパレットをはみ出さないことが求められます。

0 \le x_i + size\_x_i \le PL_x

0 \le y_i + size\_y_i \le PL_y

$x_i, y_i$ : 箱$i$の左下$x,y$座標
$PL_x, PL_y$ : パレットの$x$方向, $y$方向の寸法
$size\_x_i$, $size\_y_i$ : 箱$i$の$x$方向, $y$方向の寸法

2. 箱の90度回転を許容する制約

箱はパレットになるべく多く敷き詰めたいので箱の90度回転は可能とします。

size\_x_i = l_i(1 - r_i) + w_i r_i

size\_y_i = w_i(1 - r_i) + l_i r_i

$r_i \in {0, 1}$：箱$i$が$90^\circ$ 回転している時1、そうでない時0となる変数。
$l_i$ : 箱$i$の長辺
$w_i$ : 箱$i$の短辺

3. 箱同士が平面方向で重ならない制約

箱が貫通することを防ぐ制約です。

x_i + size\_x_i \le x_j + PL_x(1 - v_{x(i,j)})

y_i + size\_y_i \le y_j + PL_y(1 - v_{y(i,j)})

$$
v_{x(i,j)} + v_{y(i,j)} + v_{x(j,i)} + v_{y(j,i)} \ge 1
$$

$v_{x(i,j)}$, $v_{y(i,j)}\in {0, 1}$： $x$方向, $y$方向において箱$i$が箱$j$と重なっていない時1、そうでない時0となる変数。

ソースコード

本記事ではMINLPを扱えるソルバーとしてCouenneを使用しています。実行環境に応じて、Couenneが利用可能な状態に設定してください。

import matplotlib.pyplot as plt
import pandas as pd
import math
import matplotlib.patches as patches
from pyomo.environ import *

plt.rcParams['font.family'] = "MS Gothic"

import warnings
warnings.simplefilter('ignore')

# 空間、箱の寸法情報
area_x = 1000
area_y = 1000

l = 400
w = 300
box_num = 8

grid_unit = 1

# 最大公約数で割って探索空間を狭める
grid_unit = math.gcd(area_x, area_y, l, w) # mm

area_x = int(area_x / grid_unit)
area_y = int(area_y / grid_unit)
l = int(l / grid_unit)
w = int(w / grid_unit)
print(f"空間を{grid_unit}で割る")

# main
col_name = ["コード", "名称", "L", "W", "積載箱数"]
l = [[111111, "a", l, w, box_num]]

df = pd.DataFrame(data=l, columns=col_name)
box_meisai = df[["コード", 'L', 'W','積載箱数']].values.tolist()

ids = []
l = []
w = []
s = []

for id in range(len(box_meisai)):
    for i in range(box_meisai[id][3]):
        ids += [box_meisai[id][0]]
        l += [box_meisai[id][1]]
        w += [box_meisai[id][2]]
        s += [id]
        
# モデル作成
model = ConcreteModel()

# 箱の座標と寸法
model.x = Var(range(box_num), domain=NonNegativeIntegers) 
model.y = Var(range(box_num), domain=NonNegativeIntegers)

model.size_x = Var(range(box_num), domain=Integers)  
model.size_y = Var(range(box_num), domain=Integers) 

# 箱の回転フラグ
model.r = Var(range(box_num), within=Binary)

# 重なり判定用　iがjより原点方向にある時1
model.vx = Var(range(box_num), range(box_num), within=Binary) 
model.vy = Var(range(box_num), range(box_num), within=Binary)

# 目的関数
model.OBJ = Objective(expr=0) 

# PLをはみ出さない制約
def capacity_const3(model, i):
    return model.x[i] + model.size_x[i]<=area_x
def capacity_const4(model, i):
    return model.y[i] + model.size_y[i]<=area_y

# 箱が重ならない制約
def overlap_const1(model, i, j):
    return model.x[i] + model.size_x[i] <= model.x[j] + area_x*(1-model.vx[i,j])
def overlap_const2(model, i, j):
    return model.y[i] + model.size_y[i] <= model.y[j] + area_y*(1-model.vy[i,j])
def overlap_const3(model, i, j):
    return model.vx[i, j] + model.vy[i, j] + model.vx[j, i] + model.vy[j, i] >= 1

# 寸法制約
def size_const1(model, i):
    return model.size_x[i] == l[i]*(1-model.r[i]) + w[i]*model.r[i]
def size_const2(model, i):
    return model.size_y[i] == w[i]*(1-model.r[i]) + l[i]*model.r[i]

model.Const_capacity1 = ConstraintList()
model.Const_capacity2 = ConstraintList()
model.Const_capacity3 = ConstraintList()
model.Const_capacity4 = ConstraintList()

model.Const_size1 = ConstraintList()
model.Const_size2 = ConstraintList()

model.Const_overlap1 = ConstraintList()
model.Const_overlap2 = ConstraintList()
model.Const_overlap3 = ConstraintList()

# 物理制約
for i in range(box_num):

    model.Const_capacity3.add(capacity_const3(model, i))
    model.Const_capacity4.add(capacity_const4(model, i))

    model.Const_size1.add(size_const1(model, i))
    model.Const_size2.add(size_const2(model, i))

    for j in range(box_num):
        if i != j:
            model.Const_overlap1.add(overlap_const1(model, i, j))
            model.Const_overlap2.add(overlap_const2(model, i, j))
            model.Const_overlap3.add(overlap_const3(model, i, j))     


# ソルバー指定
opt = SolverFactory("couenne", executable="couenne")

# 最適化計算を実行
res = opt.solve(model, tee=True)

# 結果
print(model.pprint())

# グラフ描画
fig = plt.figure()
ax = plt.axes()

for i in range(box_num):

    x = int(value(model.x[i])*grid_unit)
    y = int(value(model.y[i])*grid_unit)
    size_x = int(value(model.size_x[i])*grid_unit)
    size_y = int(value(model.size_y[i])*grid_unit)

    r = patches.Rectangle((x, y), size_x, size_y, edgecolor='#000000', facecolor='lightblue')
    ax.add_patch(r)

plt.axis('scaled')
ax.set_aspect('equal')

ax.set_title("平面への敷き詰め方")
ax.set_xlim(0, area_x*grid_unit)
ax.set_ylim(0, area_y*grid_unit)

plt.show()

結果

1000 × 1000 のパレットに 400 × 300 の箱を 8 箱積載するケースで最適化計算を行ったところ、図のように効率的な配置が可能であることが確認できました。

ここでの工夫は、パレットと箱の寸法を最大公約数で割ることで探索空間を縮小した点です。
例えば、パレットが 1000 × 1000、箱が 400 × 300 の場合、最大公約数は 100 なので、パレットを 10 × 10、箱を 4 × 3 として扱っても同じ意味になります。この変換により、探索空間が大幅に減り、計算効率が向上します。ただし、寸法が互いに素の場合は、この方法は適用できません。

まとめ

限られた領域に箱を敷き詰める問題を、制約充足問題として解きました。今回は目的関数を設定していないため「積み方の美しさ」は保証されていませんが、パレットと箱の寸法の関係上、結果的に整った配置になっています。積載の見た目や安定性を重視する場合は、目的関数の追加が必要です。