数理最適化の概要とPythonでの実装

数理最適化とは

数理最適化は、現実の問題を数学的なモデルに変換し、そのモデルを解析して最適な解を見つける手法です。このアプローチは、多くの分野で利用されており、製造業、物流、金融、医療、エネルギー管理など、様々な産業や領域で問題解決に役立っています。数理最適化は、リソースの効率的な利用やコスト削減、意思決定のサポートにおいて非常に有力なツールとなっています。

数理最適化の基本的な構造は、最適な決定を行うための目的関数と、その際に守る必要がある制約条件から成り立っています。目的関数は最大化または最小化したい対象を表し、制約条件は問題に課される制約を表しています（後述）。これらの数学的な表現に基づいて、最適な決定を導き出すことが数理最適化の目的です。

実用例　

例えば、製造業においては、生産コストを最小化するために生産計画を最適化する問題があります。物流では、効率的な経路を見つけて輸送コストを削減する問題があります。これらの問題は、数理最適化を用いて、最適な生産計画や経路を見つけることができます。

数理最適化の手法の種類

一般的な数理最適化の手法には、線形計画法、整数計画法、非線形計画法、ダイナミック計画法、巡回セールスマン問題などがあるそう。これらの手法は、問題の性質や大きさに応じて適用されます。例えば、線形計画法は目的関数や制約条件が線形である場合に有効であり、整数計画法は変数が整数である問題に適しています。（まだあまり詳しくない...)

流れ

1、問題の定式化

実際に解決したい問題を数式を使って表現する。
何を最適化（最大化or最小化）したいかを表す「目的関数」、最適化によって値を決める決定変数、満たすべき制約条件を数式で表します。

２、アルゴリズム実装

選択した手法を用いて最適解を求める。これにはコンピュータを使用し、アルゴリズムによって最適解を見つけます。要はプログラムを書きます。

３、分析・検証・結果の解釈

簡単な問題例に対して想定通りの最適解が得られるかを確認してみることを推奨。
実行後、得られた解が制約条件を満たしているか、実行時間が許容範囲内かを確認。
得られた最適解を解釈し、実際の問題においてどのような意味を持つかを理解する。

４、修正

得られた解の精度や実行時間について課題がある場合、定式化やコードの修正作業を行う。

以下、例題を使って実際に最適化をやってみます。
なお、今回は1、問題の定式化と２、アルゴリズム実装までを行います。

1、問題の定式化

要は「何をどう解きたいか」を考えます。
例えば、製造業において、在庫最適化を考えてみます。
A社では、在庫の量を最小化し無駄を削減したいと考えています。
（トヨタ生産方式の、「必要な時に必要な分だけ在庫する」という考えですね）
しかし当然減らしすぎると、それだけ欠品のリスクがあります。

最適化問題として考えてみると、
ここで最適化したいものを「在庫の最少化」となります。
しかし、すでに述べた通り、在庫を減らすということは、その分だけ欠品のリスクが上がります。
A社としては、例えば欠品率は5%までなら許容できると考えたとします。これが制約条件です。

まとめると、
最適化したいこと：「在庫の最少化」
制約条件：欠品率5%まで
これらを満たすためには、常にどのぐらい在庫しておけば良いかを解きます。

もちろん、「在庫最適化」は在庫の最少化だけではありません。

例えばB社では、
欠品により納期遅れが多発しており、それを改善したいとします。
もちろん欠品をなくすためには大量に在庫しておけば解決できそうですが、保管コストや倉庫の容量に制限（上限）があります。これが制約条件となります。

つまりB社の場合、最適化の解くべき内容は
最適化したいこと：欠品率の最小化
制約条件：倉庫の容量や保管コストの上限
となります。

このように、同じ「在庫最適化」といっても、抱えている問題によって解くべき内容は変わってきます。もちろん「在庫最適化」以外のテーマでも同様です。

つまり、自社がどのような問題を抱えていて、「何を解くべきか」をはっきりさせることが重要です。
ここがずれていると、たとえ数理最適化で問題を解いたとしても、実際に活用して効果を得ることができません。

例題

実際に例題を使って解いてみます。
（オリジナルですので、その設定おかしいだろというのは心に留めておいてください。。。）

製品A、Bの２つがある。

	サイズ	1個あたり利益	1日あたり需要量
製品A	70	1,500円	100個
製品B	50	800円	150個

	サイズ
倉庫	10,000

必然的に欠品してしまうが、増設しろよ 倉庫を拡張したいが間に合わないので、なんとか機会損失を最小にしたい。

製品A、Bそれぞれいくつずつ在庫として置いておくべきか。

最適化したいこと：機会損失の最小化（＝利益の最大化）
制約条件：倉庫の容量が10,000

※製品と倉庫のサイズ単位は何にすべきか迷ったが、本筋ではないので一旦考慮していません。

数式化

以下、数式にまとめますが、$x^A$、$x^B$はそれぞれ製品Aと製品Bの在庫量を表します。

目的関数

Minimize $Z = 1500・max(0,100 - x^A)+800・max(0,120-x^B)$

需要量-在庫の分が不足分となり、利益を掛けることで機会損失を算出しています。
これを最小化したいのでMinimizeとしています。

制約条件

１、製品Aと製品Bのサイズ制約

$70・x^A+50・x^B \leq 10000$

製品Aと製品Bの合計が倉庫サイズである10,000を超えてはいけないので。

2、製品Aの需要制約
製品Aは需要量を超えて在庫する必要がないので、定式化しておきます。
$x^A\leq100 $

3、製品Bの需要制約
製品Bも同様
$x^B\leq120 $

4、在庫量は非負
$x^A\geq0, x^B \geq0$

以上をまとめると、以下のように定式化できます。

Minimize
$Z = 1500・max(0,100 - x^A)+800・max(0,120-x^B)$

Subject to
$70・x^A+50・x^B \leq 10000$
$x^A\leq100 $
$x^B\leq120 $
$x^A\geq0, x^B \geq0$

２、アルゴリズム実装

今回は、Pythonで実装してみます。

from scipy.optimize import linprog

# 製品A、Bのサイズ、利益、需要
size_A, size_B = 70, 50
profit_A, profit_B = 1500, 800
demand_A, demand_B = 100, 120

# 倉庫の容量
warehouse_capacity = 10000

# 目的関数の係数（機会損失を最小化するので、符号を反転）
#　inprogは最小化問題を解くために設計されているため、目的関数の符号を変更して最小化問題として扱う。
c = [-(profit_A * demand_A), -(profit_B * demand_B)]

# 制約条件の係数
A = [[size_A, size_B], [-1, 0], [0, -1]]
b = [warehouse_capacity, 0, 0]

# 在庫数は非負でなければならない
x_bounds = (0, demand_A)
y_bounds = (0, demand_B)

# 線形計画問題を解く
result = linprog(c, A_ub=A, b_ub=b, bounds=[x_bounds, y_bounds], method='highs')

# 結果を表示
inventory_A = round(result.x[0])
inventory_B = round(result.x[1])

# 製品A、Bの機会損失を計算
opportunity_loss_A = max(0, demand_A - inventory_A) * profit_A
opportunity_loss_B = max(0, demand_B - inventory_B) * profit_B

# 全体の機会損失
total_opportunity_loss = opportunity_loss_A + opportunity_loss_B

print("製品Aの在庫数:", inventory_A, "個")
print("製品Bの在庫数:", inventory_B, "個")
print("製品Aの機会損失:", opportunity_loss_A, "円")
print("製品Bの機会損失:", opportunity_loss_B, "円")
print("全体の機会損失:", total_opportunity_loss, "円")

実行結果

製品Aの在庫数: 100 個
製品Bの在庫数: 60 個
製品Aの機会損失: 0 円
製品Bの機会損失: 48000 円
全体の機会損失: 48000 円

結果として、製品Bの分だけ機会損失が出るという結果となりました。
（検証省略してしまってるんで、今度やります）

コード解説

まず、最適化に使えるライブラリであるlinprogを今回はインポートしています。
その下は、製品の情報、倉庫の容量を変数にセットしています。

from scipy.optimize import linprog

# 製品A、Bのサイズ、利益、需要
size_A, size_B = 70, 50
profit_A, profit_B = 1500, 800
demand_A, demand_B = 100, 120

# 倉庫の容量
warehouse_capacity = 10000

ここで目的関数でcという変数をセット。
コメントにも書いた通り、符号をマイナスにすることで、inprogライブラリの威力を発揮してもらう。

# 目的関数の係数（機会損失を最小化するので、符号を反転）
#　inprogは最小化問題を解くために設計されているため、目的関数の符号を変更して最小化問題として扱う。
c = [-(profit_A * demand_A), -(profit_B * demand_B)]

ここで制約条件を設定。

# 制約条件の係数
A = [[size_A, size_B], [-1, 0], [0, -1]]
b = [warehouse_capacity, 0, 0]

# 在庫数は非負でなければならない
x_bounds = (0, demand_A)
y_bounds = (0, demand_B)

ここで線形計画問題を解いています。
c: 目的関数の係数
A_ub: 不等式制約の係数行列。製品A、Bのサイズおよび在庫数に関する制約が含まれる。
b_ub: 不等式制約の右辺ベクトル。倉庫の容量制約および製品A、Bの需要に関する制約が含まれる。
bounds:変数の上下限。在庫数が非負あるとおいう制約
method: 線形計画問題を解くためのアルゴリズムを指定。今回はhighsを指定。
'highs'はHigh-Performance Interior-Point Optimizerの略で、線形計画問題の解法の一つ。
linprog関数はこれらのパラメータを受け取り、与えられた制約条件の下で目的関数を最小化（または最大化）する変数の値を見つけようとする。

# 線形計画問題を解く
result = linprog(c, A_ub=A, b_ub=b, bounds=[x_bounds, y_bounds], method='highs')

製品A、Bそれぞれの在庫量を変数に格納。

# 結果を表示
inventory_A = round(result.x[0])
inventory_B = round(result.x[1])

ちなみにresultの中身を見るとこのようになっている様。
詳しく読み解かないですが、x:[ 1.000e+02 6.000e+01]となっているところが、在庫数の解のようです。

message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
        success: True
         status: 0
            fun: -20760000.0
              x: [ 1.000e+02  6.000e+01]
            nit: 1
          lower:  residual: [ 1.000e+02  6.000e+01]
                 marginals: [ 0.000e+00  0.000e+00]
          upper:  residual: [ 0.000e+00  6.000e+01]
                 marginals: [-1.560e+04  0.000e+00]
          eqlin:  residual: []
                 marginals: []
        ineqlin:  residual: [ 0.000e+00  1.000e+02  6.000e+01]
                 marginals: [-1.920e+03 -0.000e+00 -0.000e+00]
 mip_node_count: 0
 mip_dual_bound: 0.0
        mip_gap: 0.0

ここで、機会損失を計算しています。
定式化した計算式をそのまま利用しているだけ。
下のコードは単純に製品AとBを合計しているだけ。

# 製品A、Bの機会損失を計算
opportunity_loss_A = max(0, demand_A - inventory_A) * profit_A
opportunity_loss_B = max(0, demand_B - inventory_B) * profit_B

# 全体の機会損失
total_opportunity_loss = opportunity_loss_A + opportunity_loss_B

最後に結果を表示

print("製品Aの在庫数:", inventory_A, "個")
print("製品Bの在庫数:", inventory_B, "個")
print("製品Aの機会損失:", opportunity_loss_A, "円")
print("製品Bの機会損失:", opportunity_loss_B, "円")
print("全体の機会損失:", total_opportunity_loss, "円")

最後に

今回、数理最適化をやってみましたが、実装に関しては、昔よりも飛躍的に計算速度が上がり、プログラムも書きやすくなったようですので、やはりどのような問題をどう解くかという定式化でセンスが問われそうな気がします。