LoginSignup
8
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@KTasknin株式会社LIGHTz

数理計画法(数理最適化)をはじめよう

Posted at

数理計画法(数理最適化)をはじめよう

数理計画法とはなんぞやっていうのは、詳しいところは詳しい人・資料に譲ります

日本オペレーションズリサーチ学会

しっかり学ぶ数理最適化 モデルからアルゴリズムまで (KS情報科学専門書)

組み合わせ最適化による問題解決の実践的アプローチ

数理最適化のここが良いと最近思っている

私はプログラマなので、何か問題を解こうと思ったら、長々とfor文をガリガリとコードを書いてしまうのだけど、コードを書くってやっぱり難しいし、いわゆる良いコードって曖昧だなと思う。
いわゆるフレームワークを忠実に使うことが一番最善なのだけど、何かを最適化するという意味で、フレームワークの究極系(というと雑すぎて怒られるかもしれない)のが、この数理最適化。問題を下記の数式にとにかく落とし込めば、あとはソルバーと呼ばれるものが解いてくれる。

\displaylines{
\text{min.} \quad \sum_{i=1}^{n} c_i x_i \\
\text{st.} \quad \sum_{j=1}^{n} a_{ij} x_j \leq b_i, \quad \quad i = 1, 2, \dots, m \\
x_j \geq 0, \quad  \quad j = 1, 2, \dots, n
}

実際に使ってみよう

といっても、数式とか考えたこともないって人にはなんのことやらさっぱりだと思うので、実際にプログラムとして使ってみる。まず、ソルバーを用意する。SCIPというのが無料で有名らしい(ソルバーは有料もいっぱいある、有料が主流かもしれない）、SCIPをpythonでラップしたpyscipoptを使うと楽。

長方形箱詰め問題

複数の長方形があるときにそれを四角い箱にみっちり詰め込みたいという問題を下記の資料の定式化をもとに解く。なんとこれを実行するだけで、四角い箱に長方形をみっちり詰めてくれる。
長方形箱詰め問題、長方形箱詰め問題とその定式化

import pyscipopt
class Box:
    w: int
    h: int
    x: pyscipopt.Variable
    y: pyscipopt.Variable
    topz: list[pyscipopt.Variable]
    bottomz: list[pyscipopt.Variable]
    leftz: list[pyscipopt.Variable]
    rightz: list[pyscipopt.Variable]

    def __init__(self, w, h, x, y, topz, bottomz, leftz, rightz):
        self.w = w
        self.h = h
        self.x = x
        self.y = y
        self.topz = topz
        self.bottomz = bottomz
        self.leftz = leftz
        self.rightz = rightz

    @classmethod
    def create(cls, x, y, topz, bottomz, leftz, rightz):
        return cls(randint(1, 10), randint(1, 10), x, y, topz, bottomz, leftz, rightz)

N = 10
BIG_M = 10000

def main():
    model = pyscipopt.Model("packing")
    boxes = [
        # ランダムなサイズの長方形を作成
        Box.create(
            model.addVar(name=f"x_{idx}", vtype="C"),
            model.addVar(name=f"y_{idx}", vtype="C"),
            [model.addVar(name=f"top_z_{idx}_{zidx}", vtype="B") for zidx in range(N)],
            [model.addVar(name=f"bottom_z_{idx}_{zidx}", vtype="B") for zidx in range(N)],
            [model.addVar(name=f"left_z_{idx}_{zidx}", vtype="B") for zidx in range(N)],
            [model.addVar(name=f"right_z_{idx}_{zidx}", vtype="B") for zidx in range(N)],
        ) for idx in range(N)
    ]
    # 箱のサイズを決定
    W, H = 20, model.addVar(name="H", vtype="C")
    
    # 高さを最小化する
    model.setObjective(H, sense="minimize")
    for idx0, target_box in enumerate(boxes):
        # 箱の中におさまるようにする
        model.addCons(0.0 <= target_box.y)
        model.addCons(target_box.y + target_box.h <= H)
        model.addCons(0 <= target_box.x)
        model.addCons(target_box.x + target_box.w <= W)
        for idx1, occupied_box in enumerate(boxes):
            if idx0 != idx1:
                # 長方形同士が重ならないようにする
                model.addCons(target_box.y + target_box.h <= occupied_box.y + BIG_M * (1 - target_box.bottomz[idx1]))
                model.addCons(occupied_box.y + occupied_box.h <= target_box.h + BIG_M * (1 - target_box.topz[idx1]))
                model.addCons(target_box.x + target_box.w <= occupied_box.x + BIG_M * (1 - target_box.rightz[idx1]))
                model.addCons(occupied_box.x + occupied_box.w <= target_box.x + BIG_M * (1 - target_box.leftz[idx1]))
                model.addCons(target_box.topz[idx1] + target_box.bottomz[idx1] + target_box.leftz[idx1] + target_box.rightz[idx1] == 1)

    model.optimize()

回転も考慮にいれてみよう

上記は回転を考慮に入れていないが、考慮に入れた場合どうすればよいか。
１つの長方形を、縦横が反転した２つの長方形をスイッチすることで表現することでうまいことできる。

import pyscipopt
class Box:
    w0: int
    h0: int
    w1: int
    h1: int
    x: pyscipopt.Variable
    y: pyscipopt.Variable
    topz: list[pyscipopt.Variable]
    bottomz: list[pyscipopt.Variable]
    leftz: list[pyscipopt.Variable]
    rightz: list[pyscipopt.Variable]

    def __init__(self, w, h, x, y, topz, bottomz, leftz, rightz, sw):
        self.w0, self.h0 = w, h
        self.w1, self.h1 = h, w
        self.sw = sw
        self.x = x
        self.y = y
        self.topz = topz
        self.bottomz = bottomz
        self.leftz = leftz
        self.rightz = rightz

    @classmethod
    def create(cls, x, y, topz, bottomz, leftz, rightz, sw):
        w, h = randint(1, 10), randint(1, 10)
        return cls(w, h, x, y, topz, bottomz, leftz, rightz, sw)

N = 10
BIG_M = 10000

def main():
    model = pyscipopt.Model("packing")
    boxes = [
        # ランダムなサイズの長方形を作成
        Box.create(
            model.addVar(name=f"x_{idx}", vtype="C"),
            model.addVar(name=f"y_{idx}", vtype="C"),
            [model.addVar(name=f"top_z_{idx}_{zidx}", vtype="B") for zidx in range(N)],
            [model.addVar(name=f"bottom_z_{idx}_{zidx}", vtype="B") for zidx in range(N)],
            [model.addVar(name=f"left_z_{idx}_{zidx}", vtype="B") for zidx in range(N)],
            [model.addVar(name=f"right_z_{idx}_{zidx}", vtype="B") for zidx in range(N)],
            model.addVar(name=f"sw_{idx}", vtype="B"),
        ) for idx in range(N)
    ]
    # 箱のサイズを決定
    W, H = 20, model.addVar(name="H", vtype="C")
    
    # 高さを最小化する
    model.setObjective(H, sense="minimize")
    for idx0, target_box in enumerate(boxes):
        # 箱の中におさまるようにする
        # 長方形の向きによって(w0, h0)/(w1, h1)のどちらかが選択される
        model.addCons(0.0 <= target_box.y)
        model.addCons(target_box.y + target_box.h0 <= H + BIG_M * (1 - target_box.sw))
        model.addCons(target_box.y + target_box.h1 <= H + BIG_M * target_box.sw)
        model.addCons(0 <= target_box.x)
        model.addCons(target_box.x + target_box.w0 <= W + BIG_M * (1 - target_box.sw))
        model.addCons(target_box.x + target_box.w1 <= W + BIG_M * target_box.sw)
        for idx1, occupied_box in enumerate(boxes):
            if idx0 != idx1:
                # 長方形同士が重ならないようにする
                # 長方形の向きによって(w0, h0)/(w1, h1)のどちらかが選択される
                model.addCons(target_box.y + target_box.h0 <= occupied_box.y + BIG_M * (1 - target_box.bottomz[idx1]) + BIG_M * (1 - target_box.sw))
                model.addCons(target_box.y + target_box.h1 <= occupied_box.y + BIG_M * (1 - target_box.bottomz[idx1]) + BIG_M * target_box.sw)
                model.addCons(occupied_box.y + occupied_box.h0 <= target_box.y + BIG_M * (1 - target_box.topz[idx1]) + BIG_M * (1 - occupied_box.sw))
                model.addCons(occupied_box.y + occupied_box.h1 <= target_box.y + BIG_M * (1 - target_box.topz[idx1]) + BIG_M * occupied_box.sw)
                model.addCons(target_box.x + target_box.w0 <= occupied_box.x + BIG_M * (1 - target_box.rightz[idx1]) + BIG_M * (1 - target_box.sw))
                model.addCons(target_box.x + target_box.w1 <= occupied_box.x + BIG_M * (1 - target_box.rightz[idx1]) + BIG_M * target_box.sw)
                model.addCons(occupied_box.x + occupied_box.w0 <= target_box.x + BIG_M * (1 - target_box.leftz[idx1]) + BIG_M * (1 - occupied_box.sw))
                model.addCons(occupied_box.x + occupied_box.w1 <= target_box.x + BIG_M * (1 - target_box.leftz[idx1]) + BIG_M * occupied_box.sw)
                
                model.addCons(target_box.topz[idx1] + target_box.bottomz[idx1] + target_box.leftz[idx1] + target_box.rightz[idx1] == 1)

    model.optimize()

BIG_Mがおもしろい

数理計画法でビッグMと呼ばれるテクニックがおもしろい。
下記の部分がそれにあたる。
sw=1のときは、条件AのBIG_MがBIG_M * (1 - 1)で0になる一方、条件BのBIG_MはBIG_M * 1で100000になる。つまりsw=1のとき条件Aは難しくなり、条件Bは簡単(無視できる)になる。sw=0のときはその逆。という仕組み

BIG_M = 100000
# 条件A
target_box.y + target_box.h0 <= H + BIG_M * (1 - target_box.sw)
# 条件B
target_box.y + target_box.h1 <= H + BIG_M * target_box.sw

実際使う上での問題

そもそも長方形箱詰め問題というのは個数が増えれば増えるほど計算がびっくりするぐらい遅くなる問題なので、数理計画法で解こうと思うと20個ぐらいでいっきに計算が辛い感じになっちゃう。
詳しい人からすると、ちゃんと適切なアルゴリズムを使いましょうと言われちゃうかもしれない。

でも、同じフレームワークに押し込むだけで、いろいろな問題が解けるって大事なことだと思う

なんか、うまいやり方ないかなぁ。

今年1年の学びをアウトプット！Advent Calendar 2024

Day 7
8
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
8
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?