線形最適化と混合整数最適化

線形最適化

  • 線形最適化
  • 簡単な例題と双対最適解
  • 輸送問題,多品種輸送問題
  • 栄養問題

混合整数最適化問題

  • 簡単な例題(パズル)
  • 多制約ナップサック問題

モデリングのコツ

  • 最大値の最小化
  • 絶対値
  • 離接制約
  • if A then B 条件

Google Colab.で動かす場合には、pipでmypulpモジュールをインストールしておく。

!pip install mypulp

線形最適化ソルバー(モデラー)


  • Gurobi(商用,アカデミックフリー)のソルバー
    • 独自のPythonインターフェイス(「あたらしい数理最適化」(近代科学社)で採用)
    • 凸2次(制約)整数,2次錐最適化,非凸2次
  • PuLP (MITライセンス)のモデラー
    • メインソルバーはCBC(COIN Branch & Cut; EPLライセンス)
    • Gurobiと同じインターフェイスをもつmypulpパッケージを使う
    • ソルバーにSCIP(https://www.scipopt.org/; ZIB Academic License; CBCより高速) を使いたい場合には, 引数solverにpulp.SCIP()を入れる.

例題(線形最適化問題)

整数や非線形関数を含まない線形最適化問題は、比較的簡単に解くことができるが、大規模な実際問題を解決する際には、以下のような注意が必要である。

  1. 数値の桁数を不用意に大きくしない。例えば、数百億円規模の最適化問題において、小数点以下の桁数が非常に多い問題例においては、線形最適化でも非常に時間がかかることがある。必要最小限の桁数に丸め、目的関数値は10万(6桁)程度に設定すると、計算が安定する。

  2. 問題の疎性を考慮する。例えば、複数製品を工場から倉庫、倉庫から顧客まで運ぶサプライ・チェインの最適化を考える。固定費用を考えない場合には、線形最適化問題に帰着されるが、問題例の規模によっては解けないこともある。そういう場合には、工場で生産できない製品や、顧客需要のない製品に対しては経路からあらかじめ除外してネットワークを生成すると劇的に速度が改善されることがある。また、輸送費用が変わらない顧客を集約したり、製品のABC分析を行い、需要の小さい製品はとりあえず除外するか、集約して扱うことも推奨される。

  3. 実行不可能にならないように注意する。実際問題では、ユーザーは無理な注文をしがちだ。それをそのまま真に受けて定式化すると実行不可能になることがある。数理最適化ソルバーは、制約を満たす解が存在しないしないことを示してくれるが、どの制約のせいで解がないのかまでは示してくれない。Gurobiには、Irreducible Inconsistent Subsystem (IIS)を計算してくれるメソッドが準備されているが、定式化を行う際に、主要な制約の逸脱を許してソフト制約にしておくことを推奨する。また、全ての費用を合計した1つの目的関数を設定するのではなく、輸送費用、配送費用、倉庫費用など小分けにして計算して変数に保存しておき、最後にそれを合計したものを最適化すると、あとでどの個別の費用を計算する手間が省ける。

  • トンコケ,コケトン,ミックスの丼を販売
  • 資源制約の下で,利益を最大化

変数

  • トンコケ丼 $x_1$,コケトン丼 $x_2$,ミックス丼 $x_3$

  • 定式化

$$ \begin{array}{l c c c c c} maximize & 15 x_1 & + 18 x_2 & +30 x_3 & & \\ s.t. & 2x_1 & + x_2 & + x_3 & \leq & 60 \\ & x_1 & + 2 x_2 & + x_3 &\leq & 60 \\ & & & x_3 &\leq & 30 \\ & & & x_1,x_2,x_3 & \geq & 0 \end{array} $$
from mypulp import GRB, Model, quicksum, multidict, tuplelist

model = Model("lo1")

x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
x3 = model.addVar(ub=30.0, name="x3")

model.update()  # Gurobiの怠惰な更新(lazy update)という仕様(忘れずに!)

model.addConstr(2 * x1 + x2 + x3 <= 60)
# 別の定義方法 1
# L1 = LinExpr([2,1,1],[x1,x2,x3]) #線形表現(式)
# 別の定義方法 2
# L1 = LinExpr()     #線形表現(式)
# L1.addTerms(2,x1)  #項を追加
# L1.addTerms(1,x2)
# L1.addTerms(1,x3)
# model.addConstr(lhs=L1,sense='<',rhs=60)  #制約を追加

model.addConstr(x1 + 2 * x2 + x3 <= 60)

model.setObjective(15 * x1 + 18 * x2 + 30 * x3, GRB.MAXIMIZE)

# SCIPを使用したい場合
# import pulp
# solver = pulp.SCIP()
# model.optimize(solver)

model.optimize()

if model.Status == GRB.Status.OPTIMAL:
    print("Opt. Value=", model.ObjVal)
    for v in model.getVars():
        print(v.VarName, v.X)
Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - /Users/mikiokubo/Library/Caches/pypoetry/virtualenvs/analytics-v-sH3Dza-py3.8/lib/python3.8/site-packages/pulp/apis/../solverdir/cbc/osx/64/cbc /var/folders/69/5y96sdc94jxf6khgc8mlmxrr0000gn/T/4475c5f229824915a8fe2d72d4485e3a-pulp.mps max timeMode elapsed branch printingOptions all solution /var/folders/69/5y96sdc94jxf6khgc8mlmxrr0000gn/T/4475c5f229824915a8fe2d72d4485e3a-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 7 COLUMNS
At line 17 RHS
At line 20 BOUNDS
At line 24 ENDATA
Problem MODEL has 2 rows, 3 columns and 6 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Presolve 2 (0) rows, 3 (0) columns and 6 (0) elements
0  Obj -0 Dual inf 63 (3)
0  Obj -0 Dual inf 63 (3)
3  Obj 1230
Optimal - objective value 1230
Optimal objective 1230 - 3 iterations time 0.002
Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

Opt. Value= 1230.0
x1 10.0
x2 10.0
x3 30.0

モデルファイルの出力


  • model.write('ファイル名.lp')で LPフォーマット(Linear Programming (LP) format)で保存

  • model.write('ファイル名.mps')でMPSフォーマット (Mathematical Programming System (MPS) format) で保存

    • 可読性はないが,ほとんどの最適化ソルバーが対応している古典的な書式
  • 注意:Gurobiの場合にはmodel.update()を直前にするのを忘れずに

  • PuLPだと print(model) でも画面にLPフォーマットを出力

model.write("lo1.lp")
model.write("lo1.mps")
print("LP file =========================")
!cat lo1.lp
print("MPS file =========================")
!cat lo1.mps
LP file =========================
\* lo1 *\
Maximize
OBJ: 15 x1 + 18 x2 + 30 x3
Subject To
c_1: 2 x1 + x2 + x3 <= 60
c_2: x1 + 2 x2 + x3 <= 60
Bounds
 x1 <= 1e+100
 x2 <= 1e+100
 x3 <= 30
End
MPS file =========================
*SENSE:Maximize
NAME          lo1
ROWS
 N  OBJ
 L  c_1
 L  c_2
COLUMNS
    x1        c_1        2.000000000000e+00
    x1        c_2        1.000000000000e+00
    x1        OBJ        1.500000000000e+01
    x2        c_1        1.000000000000e+00
    x2        c_2        2.000000000000e+00
    x2        OBJ        1.800000000000e+01
    x3        c_1        1.000000000000e+00
    x3        c_2        1.000000000000e+00
    x3        OBJ        3.000000000000e+01
RHS
    RHS       c_1        6.000000000000e+01
    RHS       c_2        6.000000000000e+01
BOUNDS
 UP BND       x1         1.000000000000e+100
 UP BND       x2         1.000000000000e+100
 UP BND       x3         3.000000000000e+01
ENDATA
print(model)
lo1:
MAXIMIZE
15*x1 + 18*x2 + 30*x3 + 0
SUBJECT TO
c_1: 2 x1 + x2 + x3 <= 60

c_2: x1 + 2 x2 + x3 <= 60

VARIABLES
x1 <= 1e+100 Continuous
x2 <= 1e+100 Continuous
x3 <= 30 Continuous

例題 (双対問題)


資源(豚肉,鶏肉,牛肉)100グラムの価値を推定

主問題

$$ \begin{array}{l c c c c c} maximize & 15 x_1 & + 18 x_2 & +30 x_3 & & \\ s.t. & 2x_1 & + x_2 & + x_3 & \leq & 60 \\ & x_1 & + 2 x_2 & + x_3 &\leq & 60 \\ & & & x_3 &\leq & 30 \\ & & &x_1,x_2,x_3 & \geq & 0 \end{array} $$
from gurobipy import GRB, Model, quicksum, multidict, tuplelist

model = Model()

x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
x3 = model.addVar(ub=30.0, name="x3")

model.update()

c1 = model.addConstr(2 * x1 + x2 + x3 <= 60)
c2 = model.addConstr(x1 + 2 * x2 + x3 <= 60)
model.setObjective(15 * x1 + 18 * x2 + 30 * x3, GRB.MAXIMIZE)

model.optimize()

print("Opt. Value=", model.ObjVal)
for v in model.getVars():
    print(v.VarName, v.X)
Academic license - for non-commercial use only - expires 2022-04-03
Using license file /Users/mikiokubo/gurobi.lic
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 2 rows, 3 columns and 6 nonzeros
Model fingerprint: 0x1da73bc0
Coefficient statistics:
  Matrix range     [1e+00, 2e+00]
  Objective range  [2e+01, 3e+01]
  Bounds range     [3e+01, 3e+01]
  RHS range        [6e+01, 6e+01]
Presolve time: 0.01s
Presolved: 2 rows, 3 columns, 6 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    3.3000000e+31   3.000000e+30   3.300000e+01      0s
       3    1.2300000e+03   0.000000e+00   0.000000e+00      0s

Solved in 3 iterations and 0.01 seconds
Optimal objective  1.230000000e+03
Opt. Value= 1230.0
x1 10.0
x2 10.0
x3 30.0

双対問題

$$ \begin{array}{l c c c c c} minimize & 60 \pi_1 & + 60 \pi_2& +30 \pi_3 & & \\ s.t. & 2\pi_1 & + \pi_2 & & \geq & 15 \\ & \pi_1 & + 2\pi_2 & &\geq & 18 \\ & \pi_1 & +\pi_2 & +\pi_3 &\geq & 30 \\ & & & \pi_1,\pi_2,\pi_3 & \geq & 0 \end{array} $$

最適双対変数は 4, 7, 19

豚肉は百グラム$400$ 円,鶏肉は百グラム $700$ 円, 牛肉は百グラム $1900$ 円の価値をもつ

for c in model.getConstrs():
    print(c.ConstrName, c.Pi)
R0 4.0
R1 7.0
print(c1, "dual=", c1.Pi)
<gurobi.Constr R0> dual= 4.0

例題 (輸送問題)


quicksummultidictを用いた一般的な記述法

顧客 $i$ 1 2 3 4 5
需要量 $d_i$ 80 270 250 160 180
工場 $j$ 輸送費用 $c_{ij}$ 容量 $M_j$
1 4 5 6 8 10 500
2 6 4 3 5 8 500
3 9 7 4 3 4 500

$x_{ij}= \mbox{工場 $j$ から顧客 $i$ に輸送される量} $

定式化

$$ \begin{array}{l l l} minimize & \displaystyle\sum_{i \in I} \displaystyle\sum_{j \in J} c_{ij} x_{ij} & \\ s.t. & \displaystyle\sum_{j \in J} x_{ij} =d_i & \forall i \in I \\ & \displaystyle\sum_{i \in I} x_{ij} \leq M_j & \forall j \in J \\ & x_{ij} \geq 0 & \forall i \in I; j \in J \end{array} $$
name, height, weight = multidict(
    {"Taro": [145, 30], "Hanako": [138, 34], "Simon": [150, 45]}
)
print(name)
print(height)
print(weight)
['Taro', 'Hanako', 'Simon']
{'Taro': 145, 'Hanako': 138, 'Simon': 150}
{'Taro': 30, 'Hanako': 34, 'Simon': 45}
model = Model()
a = [5, 4, 2]
x = [model.addVar(name=f"x({i})") for i in range(3)]
L = quicksum(a[i] * x[i] for i in range(3))
model.update()
print(L)
<gurobi.LinExpr: 5.0 x(0) + 4.0 x(1) + 2.0 x(2)>
I, d = multidict({1: 80, 2: 270, 3: 250, 4: 160, 5: 180})  # demand
J, M = multidict({1: 500, 2: 500, 3: 500})  # capacity
c = {
    (1, 1): 4,
    (1, 2): 6,
    (1, 3): 9,  # cost
    (2, 1): 5,
    (2, 2): 4,
    (2, 3): 7,
    (3, 1): 6,
    (3, 2): 3,
    (3, 3): 4,
    (4, 1): 8,
    (4, 2): 5,
    (4, 3): 3,
    (5, 1): 10,
    (5, 2): 8,
    (5, 3): 4,
}

model = Model("transportation")
x = {}
for i in I:
    for j in J:
        x[i, j] = model.addVar(vtype="C", name=f"x({i},{j})")
model.update()

for i in I:
    model.addConstr(
        quicksum(x[i, j] for j in J if (i, j) in x) == d[i], name=f"Demand({i})"
    )
for j in J:
    model.addConstr(
        quicksum(x[i, j] for i in I if (i, j) in x) <= M[j], name=f"Capacity({j})"
    )
model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), GRB.MINIMIZE)

model.optimize()
print("Optimal value:", model.ObjVal)

EPS = 1.0e-6
for (i, j) in x:
    if x[i, j].X > EPS:
        # print('{0:>5} from factory {1:>2} to customer {2:>2}'.format(x[i,j].X,j,i) )
        print(f"{x[i,j].X:>5} from factory {j:>2} to customer {i:>2}")

print(f'{"Constr.Name":>15}: {"Slack":>8} , {"Dual":>4}')
for c in model.getConstrs():
    print(f"{c.ConstrName:>15}: {c.Slack:>8} , {c.Pi:>4}")
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 8 rows, 15 columns and 30 nonzeros
Model fingerprint: 0xa225ad1d
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [3e+00, 1e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [8e+01, 5e+02]
Presolve time: 0.00s
Presolved: 8 rows, 15 columns, 30 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    3.3500000e+03   2.000000e+01   0.000000e+00      0s
       1    3.3700000e+03   0.000000e+00   0.000000e+00      0s

Solved in 1 iterations and 0.01 seconds
Optimal objective  3.370000000e+03
Optimal value: 3370.0
 80.0 from factory  1 to customer  1
 20.0 from factory  1 to customer  2
250.0 from factory  2 to customer  2
250.0 from factory  2 to customer  3
160.0 from factory  3 to customer  4
180.0 from factory  3 to customer  5
    Constr.Name:    Slack , Dual
      Demand(1):      0.0 ,  4.0
      Demand(2):      0.0 ,  5.0
      Demand(3):      0.0 ,  4.0
      Demand(4):      0.0 ,  3.0
      Demand(5):      0.0 ,  4.0
    Capacity(1):    400.0 ,  0.0
    Capacity(2):      0.0 , -1.0
    Capacity(3):    160.0 ,  0.0

例題 (多品種輸送問題)


疎な問題の扱い方と tuplelist

$ x_{ijk}= 工場 j から顧客 i に製品 k が輸送される量 $

工場1では製品2,4を,工場2では製品1,2,3を,工場3では製品2,3,4を製造可能

produce = {1:[2,4], 2:[1,2,3], 3:[2,3,4]}

定式化

$$ \begin{array}{l l l} minimize & \displaystyle\sum_{i \in I} \displaystyle\sum_{j \in J} \displaystyle\sum_{k \in K} c_{ijk} x_{ijk} & \\ s.t. & \displaystyle\sum_{j \in J} x_{ijk} =d_{ik} & \forall i \in I; k \in K \\ & \displaystyle\sum_{i \in I} \displaystyle\sum_{k \in K} x_{ijk} \leq M_j & \forall j \in J \\ & x_{ijk} \geq 0 & \forall i \in I; j \in J; k \in K \end{array} $$
T = tuplelist(
    [("Sara", "Apple"), ("Taro", "Pear"), ("Jiro", "Orange"), ("Simon", "Apple")]
)
print(T.select("*", "Apple"))
<gurobi.tuplelist (2 tuples, 2 values each):
 ( Sara  , Apple )
 ( Simon , Apple )
>
d = {
    (1, 1): 80,
    (1, 2): 85,
    (1, 3): 300,
    (1, 4): 6,
    (2, 1): 270,
    (2, 2): 160,
    (2, 3): 400,
    (2, 4): 7,
    (3, 1): 250,
    (3, 2): 130,
    (3, 3): 350,
    (3, 4): 4,
    (4, 1): 160,
    (4, 2): 60,
    (4, 3): 200,
    (4, 4): 3,
    (5, 1): 180,
    (5, 2): 40,
    (5, 3): 150,
    (5, 4): 5,
}

I = set([i for (i, k) in d])
K = set([k for (i, k) in d])
J, M = multidict({1: 3000, 2: 3000, 3: 3000})  # capacity
produce = {
    1: [2, 4],
    2: [1, 2, 3],
    3: [2, 3, 4],
}  # products that can be produced in each facility
weight = {1: 5, 2: 2, 3: 3, 4: 4}  # {commodity: weight<float>}

cost = {
    (1, 1): 4,
    (1, 2): 6,
    (1, 3): 9,  # {(customer,factory): cost<float>}
    (2, 1): 5,
    (2, 2): 4,
    (2, 3): 7,
    (3, 1): 6,
    (3, 2): 3,
    (3, 3): 4,
    (4, 1): 8,
    (4, 2): 5,
    (4, 3): 3,
    (5, 1): 10,
    (5, 2): 8,
    (5, 3): 4,
}

c = {}
for i in I:
    for j in J:
        for k in produce[j]:
            c[i, j, k] = cost[i, j] * weight[k]
model = Model("multi-commodity transportation")

x = {}
for (i, j, k) in c:
    x[i, j, k] = model.addVar(vtype="C", name=f"x({i},{j},{k})")
model.update()
arcs = tuplelist([(i, j, k) for (i, j, k) in x])
for i in I:
    for k in K:
        model.addConstr(
            quicksum(x[i, j, k] for (i, j, k) in arcs.select(i, "*", k)) == d[i, k],
            f"Demand({i},{k})",
        )
for j in J:
    model.addConstr(
        quicksum(x[i, j, k] for (i, j, k) in arcs.select("*", j, "*")) <= M[j],
        f"Capacity({j})",
    )

model.setObjective(quicksum(c[i, j, k] * x[i, j, k] for (i, j, k) in x), GRB.MINIMIZE)
model.optimize()
print("Optimal value:", model.ObjVal)

EPS = 1.0e-6
for (i, j, k) in x:
    if x[i, j, k].X > EPS:
        print(f"{x[i,j,k].X} units of commodity {k} from {j} to {i}")
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 23 rows, 40 columns and 80 nonzeros
Model fingerprint: 0x77b0a6c4
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [6e+00, 4e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [3e+00, 3e+03]
Presolve removed 23 rows and 40 columns
Presolve time: 0.00s
Presolve: All rows and columns removed
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    4.3536000e+04   0.000000e+00   0.000000e+00      0s

Solved in 0 iterations and 0.01 seconds
Optimal objective  4.353600000e+04
Optimal value: 43536.0
85.0 units of commodity 2 from 1 to 1
6.0 units of commodity 4 from 1 to 1
80.0 units of commodity 1 from 2 to 1
300.0 units of commodity 3 from 2 to 1
7.0 units of commodity 4 from 1 to 2
270.0 units of commodity 1 from 2 to 2
160.0 units of commodity 2 from 2 to 2
400.0 units of commodity 3 from 2 to 2
250.0 units of commodity 1 from 2 to 3
130.0 units of commodity 2 from 2 to 3
350.0 units of commodity 3 from 2 to 3
4.0 units of commodity 4 from 3 to 3
160.0 units of commodity 1 from 2 to 4
60.0 units of commodity 2 from 3 to 4
200.0 units of commodity 3 from 3 to 4
3.0 units of commodity 4 from 3 to 4
180.0 units of commodity 1 from 2 to 5
40.0 units of commodity 2 from 3 to 5
150.0 units of commodity 3 from 3 to 5
5.0 units of commodity 4 from 3 to 5
arcs.select(1, 1, "*")
<gurobi.tuplelist (2 tuples, 3 values each):
 ( 1 , 1 , 2 )
 ( 1 , 1 , 4 )
>
d = {
    (1, 1): 80,
    (1, 2): 85,
    (1, 3): 300,
    (1, 4): 6,
    (2, 1): 270,
    (2, 2): 160,
    (2, 3): 400,
    (2, 4): 7,
    (3, 1): 250,
    (3, 2): 130,
    (3, 3): 350,
    (3, 4): 4,
    (4, 1): 160,
    (4, 2): 60,
    (4, 3): 200,
    (4, 4): 3,
    (5, 1): 180,
    (5, 2): 40,
    (5, 3): 150,
    (5, 4): 5,
}

I = set([i for (i, k) in d])
K = set([k for (i, k) in d])
J, M = multidict({1: 3000, 2: 3000, 3: 3000})  # capacity
produce = {
    1: [2, 4],
    2: [1, 2, 3],
    3: [2, 3, 4],
}  # products that can be produced in each facility
weight = {1: 5, 2: 2, 3: 3, 4: 4}  # {commodity: weight<float>}

cost = {
    (1, 1): 4,
    (1, 2): 6,
    (1, 3): 9,  # {(customer,factory): cost<float>}
    (2, 1): 5,
    (2, 2): 4,
    (2, 3): 7,
    (3, 1): 6,
    (3, 2): 3,
    (3, 3): 4,
    (4, 1): 8,
    (4, 2): 5,
    (4, 3): 3,
    (5, 1): 10,
    (5, 2): 8,
    (5, 3): 4,
}

c = {}
for i in I:
    for j in J:
        for k in produce[j]:
            c[i, j, k] = cost[i, j] * weight[k]
model = Model("multi-commodity transportation")

x = {}
for (i, j, k) in c:
    x[i, j, k] = model.addVar(vtype="C", name=f"x({i},{j},{k})")
model.update()

for i in I:
    for k in K:
        model.addConstr(
            quicksum(x[i, j, k] for j in J if (i, j, k) in x) == d[i, k],
            f"Demand({i},{k})",
        )
for j in J:
    model.addConstr(
        quicksum(x[i, j, k] for i in I for k in K if (i, j, k) in x) <= M[j],
        f"Capacity({j})",
    )

model.setObjective(quicksum(c[i, j, k] * x[i, j, k] for (i, j, k) in x), GRB.MINIMIZE)
model.optimize()
print("Optimal value:", model.ObjVal)

EPS = 1.0e-6
for (i, j, k) in x:
    if x[i, j, k].X > EPS:
        print(f"{x[i,j,k].X} units of commodity {k} from {j} to {i}")
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 23 rows, 40 columns and 80 nonzeros
Model fingerprint: 0x77b0a6c4
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [6e+00, 4e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [3e+00, 3e+03]
Presolve removed 23 rows and 40 columns
Presolve time: 0.02s
Presolve: All rows and columns removed
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    4.3536000e+04   0.000000e+00   0.000000e+00      0s

Solved in 0 iterations and 0.02 seconds
Optimal objective  4.353600000e+04
Optimal value: 43536.0
85.0 units of commodity 2 from 1 to 1
6.0 units of commodity 4 from 1 to 1
80.0 units of commodity 1 from 2 to 1
300.0 units of commodity 3 from 2 to 1
7.0 units of commodity 4 from 1 to 2
270.0 units of commodity 1 from 2 to 2
160.0 units of commodity 2 from 2 to 2
400.0 units of commodity 3 from 2 to 2
250.0 units of commodity 1 from 2 to 3
130.0 units of commodity 2 from 2 to 3
350.0 units of commodity 3 from 2 to 3
4.0 units of commodity 4 from 3 to 3
160.0 units of commodity 1 from 2 to 4
60.0 units of commodity 2 from 3 to 4
200.0 units of commodity 3 from 3 to 4
3.0 units of commodity 4 from 3 to 4
180.0 units of commodity 1 from 2 to 5
40.0 units of commodity 2 from 3 to 5
150.0 units of commodity 3 from 3 to 5
5.0 units of commodity 4 from 3 to 5

例題 (栄養問題: 実行不可能性の取り扱い)


あなたは,某ハンバーガーショップの調査を命じられた健康オタクの諜報員だ.

あなたは任務のため,毎日ハンバーガーショップだけで食事をしなければならないが, 健康を守るため,なるべく政府の決めた栄養素の推奨値を遵守しようと考えている.

考慮する栄養素は,カロリー(Cal),炭水化物(Carbo),タンパク質(Protein), ビタミンA(VitA),ビタミンC(VitC),カルシウム(Calc),鉄分(Iron)であり, 1 日に必要な量の上下限は,以下の表の通りとする.

現在,ハンバーガーショップで販売されている商品は, CQPounder, Big M, FFilet, Chicken, Fries, Milk, VegJuice の6 種類だけであり, それぞれの価格と栄養素の含有量は,以下の表のようになっている.

さらに,調査費は限られているので,なるべく安い商品を購入するように命じられている. さて,どの商品を購入して食べれば,健康を維持できるだろうか?

栄養素 $N$ Cal Carbo Protein VitA VitC Calc Iron 価格 
商品名 $F$ $n_{ij}$ $c_j$ 
CQPounder 556 39 30 147 10 221 2.4 360
Big M 556 46 26 97 9 142 2.4 320
FFilet 356 42 14 28 1 76 0.7 270
Chicken 431 45 20 9 2 37 0.9 290
Fries 249 30 3 0 5 7 0.6 190
Milk 138 10 7 80 2 227 0 170
VegJuice 69 17 1 750 2 18 0 100
上限 $a_i$ 3000 375 60 750 100 900 7.5
下限 $b_i$ 2000 300 50 500 85 660 6.0

定式化

$x_j$ は商品 $j$ の購入量(実数)

$$ \begin{array}{l l l} minimize & \displaystyle\sum_{j \in F} c_j x_j & \\ s.t. & a_i \leq \displaystyle\sum_{j \in F} n_{ij} x_j \leq b_i & i \in N \\ & x_j \geq 0 & j \in F \end{array} $$

$d_i$: 不足変数

$s_i$: 超過変数

改良した制約

$$a_i - d_i \leq \displaystyle\sum_{j \in F} n_{ij} x_j \leq b_i +s_i \ \ i \in N$$

変更した目的関数(Mは大きな数)

$$ \begin{array}{l l l} minimize & \displaystyle\sum_{j \in F} c_j x_j + \displaystyle\sum_{i \in N} M (d_i+s_i) & \end{array} $$
F, c, n = multidict(
    {
        "CQPounder": [
            360,
            {
                "Cal": 556,
                "Carbo": 39,
                "Protein": 30,
                "VitA": 147,
                "VitC": 10,
                "Calc": 221,
                "Iron": 2.4,
            },
        ],
        "Big M": [
            320,
            {
                "Cal": 556,
                "Carbo": 46,
                "Protein": 26,
                "VitA": 97,
                "VitC": 9,
                "Calc": 142,
                "Iron": 2.4,
            },
        ],
        "FFilet": [
            270,
            {
                "Cal": 356,
                "Carbo": 42,
                "Protein": 14,
                "VitA": 28,
                "VitC": 1,
                "Calc": 76,
                "Iron": 0.7,
            },
        ],
        "Chicken": [
            290,
            {
                "Cal": 431,
                "Carbo": 45,
                "Protein": 20,
                "VitA": 9,
                "VitC": 2,
                "Calc": 37,
                "Iron": 0.9,
            },
        ],
        "Fries": [
            190,
            {
                "Cal": 249,
                "Carbo": 30,
                "Protein": 3,
                "VitA": 0,
                "VitC": 5,
                "Calc": 7,
                "Iron": 0.6,
            },
        ],
        "Milk": [
            170,
            {
                "Cal": 138,
                "Carbo": 10,
                "Protein": 7,
                "VitA": 80,
                "VitC": 2,
                "Calc": 227,
                "Iron": 0,
            },
        ],
        "VegJuice": [
            100,
            {
                "Cal": 69,
                "Carbo": 17,
                "Protein": 1,
                "VitA": 750,
                "VitC": 2,
                "Calc": 18,
                "Iron": 0,
            },
        ],
    }
)
N, a, b = multidict(
    {
        "Cal": [2000, 3000],
        "Carbo": [300, 375],
        "Protein": [50, 60],
        "VitA": [500, 750],
        "VitC": [85, 100],
        "Calc": [660, 900],
        "Iron": [6.0, 7.5],
    }
)
model = Model("modern_diet")
x, s, d = {}, {}, {}
for j in F:
    x[j] = model.addVar(vtype="C", name=f"x({j})")
for i in N:
    s[i] = model.addVar(vtype="C", name=f"surplus({i})")
    d[i] = model.addVar(vtype="C", name=f"deficit({i})")
model.update()
for i in N:
    model.addConstr(quicksum(n[j][i] * x[j] for j in F) >= a[i] - d[i], f"NutrLB({i})")
    model.addConstr(quicksum(n[j][i] * x[j] for j in F) <= b[i] + s[i], f"NutrUB({i})")
model.setObjective(
    quicksum(c[j] * x[j] for j in F) + quicksum(9999 * d[i] + 9999 * s[i] for i in N),
    GRB.MINIMIZE,
)
model.optimize()
status = model.Status
if status == GRB.Status.OPTIMAL:
    print("Optimal value:", model.ObjVal)
    for j in F:
        if x[j].X > 0:
            print(j, x[j].X)
    for i in N:
        if d[i].X > 0:
            print(f"deficit of {i} ={d[i].X}")
        if s[i].X > 0:
            print(f"surplus of {i} ={s[i].X}")
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 14 rows, 21 columns and 106 nonzeros
Model fingerprint: 0x3629e088
Coefficient statistics:
  Matrix range     [6e-01, 8e+02]
  Objective range  [1e+02, 1e+04]
  Bounds range     [0e+00, 0e+00]
  RHS range        [6e+00, 3e+03]
Presolve time: 0.00s
Presolved: 14 rows, 21 columns, 106 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    0.0000000e+00   2.553750e+02   0.000000e+00      0s
      11    2.6511919e+05   0.000000e+00   0.000000e+00      0s

Solved in 11 iterations and 0.01 seconds
Optimal objective  2.651191876e+05
Optimal value: 265119.1875926751
CQPounder 0.013155307054175128
Fries 10.422664500064354
Milk 2.5154631133990755
VegJuice 0.7291054943881469
deficit of VitC =26.265987213562028

例題 (混合整数最適化問題)

宇宙ステーションに、地球人(頭ひとつで足2本)、犬星人(頭ひとつで足4本)、タコ星人(頭1つで足8本)が何人かいる。

頭の数の合計は32、足の数の合計は80である。人間以外の星人の数が最小の場合の、各星人の人数を求めよ。

model = Model("puzzle")
x = model.addVar(vtype="I", name="x")
y = model.addVar(vtype="I", name="y")
z = model.addVar(vtype="I", name="z")
model.update()

model.addConstr(x + y + z == 32, "Heads")
model.addConstr(2 * x + 4 * y + 8 * z == 80, "Legs")

model.setObjective(y + z, GRB.MINIMIZE)

model.optimize()

print("Opt. Val.=", model.ObjVal)
print("(x,y,z)=", (x.X, y.X, z.X))
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 2 rows, 3 columns and 6 nonzeros
Model fingerprint: 0x9688a20c
Variable types: 0 continuous, 3 integer (0 binary)
Coefficient statistics:
  Matrix range     [1e+00, 8e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [3e+01, 8e+01]
Presolve removed 2 rows and 3 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.01 seconds
Thread count was 1 (of 16 available processors)

Solution count 1: 4 

Optimal solution found (tolerance 1.00e-04)
Best objective 4.000000000000e+00, best bound 4.000000000000e+00, gap 0.0000%
Opt. Val.= 4.0
(x,y,z)= (28.0, 2.0, 2.0)

例題 (多制約ナップサック問題)

あなたは,ぬいぐるみ専門の泥棒だ. ある晩,あなたは高級ぬいぐるみ店にこっそり忍び込んで,盗む物を選んでいる.

狙いはもちろん,マニアの間で高額で取り引きされているクマさん人形だ. クマさん人形は,現在4体販売されていて, それらの値段と重さと容積は, 以下のようになっている.

  • 極小クマ: 16, 2, 3000
  • 小クマ: 19, 3, 3500
  • 中クマ: 23, 5100
  • 大クマ: 28, 5, 7200

あなたは,転売価格の合計が最大になるようにクマさん人形を選んで逃げようと思っているが, あなたが逃走用に愛用しているナップサックはとても古く, $7$ kgより重い荷物を入れると,底がぬけてしまうし, 10000 $cm^3$ を超えた荷物を入れると破けてしまう.

さて,どのクマさん人形をもって逃げれば良いだろうか?


  • $n$ 個のアイテム,$m$本の制約
  • 各々のアイテム $j=1,2,\ldots,n$ の価値 $v_j (\geq 0)$,
  • アイテム $j$ の制約 $i=1,2,\ldots,m$ に対する重み $a_{ij} (\geq 0)$
  • 制約 $i$ に対する制約の上限値 $b_i (\geq 0)$
  • アイテムの番号の集合 $J =\{1,2,\ldots,n \}$
  • 制約の番号の集合 $I =\{1,2,\ldots,m\}$
  • アイテム $j$ をナップサックに詰めるとき $1$, それ以外のとき $0$ になる $0$-$1$ 変数 $x_j$

定式化

$$ \begin{array}{l l l} maximize & \displaystyle\sum_{j \in J} v_j x_j & \\ s.t. & \displaystyle\sum_{j \in J} a_{ij} x_j \leq b_i & \forall i \in I \\ & x_j \in \{0,1\} & \forall j \in J \end{array} $$
J, v = multidict({1: 16, 2: 19, 3: 23, 4: 28})
a = {
    (1, 1): 2,
    (1, 2): 3,
    (1, 3): 4,
    (1, 4): 5,
    (2, 1): 3000,
    (2, 2): 3500,
    (2, 3): 5100,
    (2, 4): 7200,
}
I, b = multidict({1: 7, 2: 10000})

model = Model("mkp")
x = {}
for j in J:
    x[j] = model.addVar(vtype="B", name=f"x({j})")
model.update()

for i in I:
    model.addConstr(quicksum(a[i, j] * x[j] for j in J) <= b[i], f"Capacity({i})")

model.setObjective(quicksum(v[j] * x[j] for j in J), GRB.MAXIMIZE)

model.optimize()
print("Optimal value=", model.ObjVal)

EPS = 1.0e-6
for v in model.getVars():
    if v.X > EPS:
        print(v.VarName, v.X)
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 2 rows, 4 columns and 8 nonzeros
Model fingerprint: 0xc3121b0b
Variable types: 0 continuous, 4 integer (4 binary)
Coefficient statistics:
  Matrix range     [2e+00, 7e+03]
  Objective range  [2e+01, 3e+01]
  Bounds range     [1e+00, 1e+00]
  RHS range        [7e+00, 1e+04]
Found heuristic solution: objective 35.0000000
Presolve removed 2 rows and 4 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.01 seconds
Thread count was 1 (of 16 available processors)

Solution count 2: 42 35 

Optimal solution found (tolerance 1.00e-04)
Best objective 4.200000000000e+01, best bound 4.200000000000e+01, gap 0.0000%
Optimal value= 42.0
x(2) 1.0
x(3) 1.0

例題 論理条件

上の多制約 $0$-$1$ ナップサック問題の例題に対して, 以下の付加条件をつけた問題を定式化して解を求めよ.

  1. 小クマと中クマは仲が悪いので,同時に持って逃げてはいけない.
  2. 小クマは極小クマが好きなので,小クマを持って逃げるときには必ず極小クマももって逃げなければならない.
  3. 極小クマ,小クマ,大クマのうち,少なくとも2つは持って逃げなければならない.
J, v = multidict({1: 16, 2: 19, 3: 23, 4: 28})
a = {
    (1, 1): 2,
    (1, 2): 3,
    (1, 3): 4,
    (1, 4): 5,
    (2, 1): 3000,
    (2, 2): 3500,
    (2, 3): 5100,
    (2, 4): 7200,
}
I, b = multidict({1: 7, 2: 10000})

model = Model("mkp")
x = {}
for j in J:
    x[j] = model.addVar(vtype="B", name=f"x({j})")
model.update()

for i in I:
    model.addConstr(quicksum(a[i, j] * x[j] for j in J) <= b[i], f"Capacity({i})")

model.setObjective(quicksum(v[j] * x[j] for j in J), GRB.MAXIMIZE)

# いずれかのコメントを外す
# 1.
# model.addConstr( x[2] + x[3] <=1 )
# 2.
# model.addConstr( x[2] <= x[1])
# 3.
model.addConstr(x[1] + x[2] + x[4] >= 2)

model.optimize()
print("Optimal value=", model.ObjVal)

EPS = 1.0e-6
for v in model.getVars():
    if v.X > EPS:
        print(v.VarName, v.X)
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 3 rows, 4 columns and 11 nonzeros
Model fingerprint: 0xc27c0c02
Variable types: 0 continuous, 4 integer (4 binary)
Coefficient statistics:
  Matrix range     [1e+00, 7e+03]
  Objective range  [2e+01, 3e+01]
  Bounds range     [1e+00, 1e+00]
  RHS range        [2e+00, 1e+04]
Found heuristic solution: objective 35.0000000
Presolve removed 3 rows and 4 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.00 seconds
Thread count was 1 (of 16 available processors)

Solution count 1: 35 

Optimal solution found (tolerance 1.00e-04)
Best objective 3.500000000000e+01, best bound 3.500000000000e+01, gap 0.0000%
Optimal value= 35.0
x(1) 1.0
x(2) 1.0

問題 線形最適化(1)

あなたは業務用ジュースの販売会社の社長だ. いま,原料の ぶどう原液 $200$ リットルとりんご原液 $100$ リットルを使って 2 種類のミックスジュース(商品名はA,B)を作ろうと思っている.

ジュースAを作るにはぶどう $3$ リットルとりんご $1$ リットルが必要で, ジュースBを作るにはぶどう $2$ リットルとりんご $2$ リットルが必要である. (なんとジュースは $4$ リットル入りの特大サイズなのだ!)

ジュースAは $3$ 千円,ジュースBは$4$ 千円の値段をつけた.

さて,ジュースAとジュースBをそれぞれ何本作れば,利益が最大になるだろうか?

問題 線形最適化(2)

AさんとBさんとCさんがいくらかずつお金を持っている. AさんがBさんに $100$ 円をわたすとすると,AさんBさんの所持金が等しくなる. また,BさんがCさんに $300$ 円わたすと,BさんCさんの所持金が等しくなる. 上の条件を満たす中で,3 人の所持金の和が最小になるものを求めよ.

問題 線形最適化(3)

マラソン大会に出場した裕一郎君の証言をもとに,彼がどれくらい休憩していたかを推定せよ. ただし時間はすべて整数でなく実数で測定するものとする.

  • 証言1: フルマラソンの $42.195$kmはきつかったです.でもタイムは $6$時間 $40$分と自己ベスト更新です.
  • 証言2: できるだけ走ろうと頑張りましたが,ときどき歩いたり,休憩をとっていました.
  • 証言3: 歩いている時間は走ってる時間のちょうど2倍でした.
  • 証言4: 僕の歩く速度は分速$70$メートルで,走る速度は分速 $180$メートルです.

問題 線形最適化(4)

あなたは丼チェーンの店長だ.店の主力製品は, トンコケ丼,コケトン丼,ミックス丼,ビーフ丼の4 種類で, トンコケ丼を1 杯作るには,$200$ グラムの豚肉と $100$ グラムの鶏肉, コケトン丼を1 杯作るには,$100$ グラムの豚肉と $200$ グラムの鶏肉, ミックス丼を1 杯作るには,豚肉,鶏肉,牛肉を $100$ グラムずつ, 最後のビーフ丼は,牛肉だけを$300$ グラム使う. ただし,ビーフ丼は限定商品のため1 日 $10$ 杯しか作れない.

原料として使用できる豚,鶏,牛の肉は,最大1 日あたり $9$ キログラム,$9$ キログラム,$6$ キログラムで, 販売価格は,トンコケ丼1 杯 $1500$ 円,コケトン丼1 杯 $1800$ 円,ミックス丼1 杯 $2000$ 円, そしてビーフ丼は $5000$ 円だ.

さて,お店の利益を最大にするためには,あなたは丼を何杯ずつ作るように指示を出せばよいのだろうか?

また、上の問題において, 豚肉,鶏肉,牛肉の $100$ グラムあたりの価値はいくらになるか計算せよ.

問題 鶴亀蛸キメラ算

鶴と亀と蛸とキメラ(伝説に出てくる空想生物)が何匹ずつかいる. 頭の数を足すと $32$,足の数を足すと $99$ になる. キメラの頭の数は $2$,足の数は $3$ としたときに, 鶴と亀と蛸の数の和を一番小さくするような匹数を求めよ.

問題 輸送問題と双対性

あなたは,スポーツ用品販売チェインのオーナーだ. あなたは,店舗展開をしている5つの顧客 (需要地点)に対して, 3 つの自社工場で生産した1種類の製品を運ぶ必要がある. 工場の生産可能量(容量)と顧客の需要量は表のようになっている.

$$ \begin{array}{c c c | c c c c c } 工場生産可能容量 & & & 顧需要客 & & & & \\ \hline 1 & 2 & 3 & 1 & 2 & 3 & 4 & 5 \\ \hline 500 & 500 & 500 & 80 & 270 & 250 & 160 & 180 \\ \hline \end{array} $$

ただし,工場から顧客へ製品を輸送する際は必ず2箇所の倉庫のいずれかを経由しなければならない. 工場と倉庫間,倉庫と顧客間の輸送費用は, 以下の表のようになっているとしたとき, どのような輸送経路を選択すれば,総費用が最小になるであろうか?

$$ \begin{array}{c | c c c | c c c c c } & 工場 & & & 顧客 & & & & \\ \hline 倉庫 & 1 & 2 & 3 & 1 & 2 & 3 & 4 & 5 \\ \hline 1 & 1 & 2 & 3 & 4 & 5 & 6 & 8 & 10 \\ 2 & 3 & 1 & 2 & 6 & 4 & 3 & 5 & 8 \\ \hline \end{array} $$

上の輸送問題において,各顧客の追加注文は $1$ 単位あたりいくらの費用の増加をもたらすのか,双対性の概念を用いて計算せよ.

費用削減のためには,どの工場を拡張すれば良いかを,やはり双対性を用いて考えよ.

問題 実行不可能性

上の輸送問題において,各顧客の需要がすべて2 倍になった場合を考えよ. この問題は実行不可能になるので,それを回避する定式化を示せ.

ヒント: 工場容量の逸脱を許すか,顧客需要の不足を許す定式化を考えれば良い.

問題 ゼロ和ゲーム (1)

サッカーのPK(ペナルティキック)の最適戦略を考える.

いま,ゴールキーパーは(キッカーから見て)左に飛ぶか右に飛ぶかの2 つの戦略を持っており, キッカーは左に蹴るか右に蹴るかの2 つの戦略を持っているものとする. 得点が入る確率は,両選手の得意・不得意から以下のような確率になっているものとする(行がキーパーで,列がキッカーの戦略である).

$$ \begin{array}{ c | c c} & 左 & 右 \\ \hline 左 & 0.9 & 0.5 \\ 右 & 0.6 & 0.8 \\ \hline \end{array} $$

キーパーは得点が入る確率を最小化したいし,キッカーは得点が入る確率を最大化したい. さて,両選手はどのような行動をとれば良いだろうか?

ヒント: これは,ゼロ和ゲームの混合戦略を求める問題となり,確率的な行動をとることが最適戦略となる. (つまりどちらに飛ぶかはサイコロを振って決めるのだ.) キーパーが左に行く確率を変数 $L$,右に行く確率を変数 $R$ として線形最適化問題として定式化を行う. キッカーの戦略を(たとえば左に)固定した場合には,ゲームの値(お互いが最適戦略をとったときにゴールが決まる確率)$V$ は, $$ V \geq 0.9 L + 0.6 R $$ を満たす.キーパーはなるべく $V$ が小さくなるように変数を決め,キッカーは $V$ の値を最大化する.

問題 ゼロ和ゲーム (2)

上の問題において,キッカーの戦略を変数として定式化を行い求解せよ. 元の(キーパーの戦略を変数とした)問題の最適双対変数が, この問題の最適解になっていることを確認せよ.

例題 最大値の最小化

実数変数 $x_1, x_2$ に対する2つの線形関数 $3x_1+4x_2$ と $2x_1+7x_2$ の大きい方を小さくしたい場合を考える. これは最大値の最小化問題であり,新しい実数変数 $z$ を導入し, $$ 3x_1 + 4x_2 \leq z, \ \ 2x_1 + 7x_2 \leq z $$ の制約を加えた後で, $z$ を最小化することによって,通常の線形最適化に帰着できる. これは,制約が何本あっても同じである.また,最大値の最小化も同様に定式化できる.

整数変数に対する線形式に対しては,「最大値の最小化」をしたいときに上のアプローチをとることは推奨されない.(小規模な問題例は別である.) これは,混合整数最適化ソルバーで用いている分枝限定法が,このタイプの制約に対して弱く,限界値の改善が難しいためである.

例題: 以下の制約の下で, $3x_1+4x_2$ と $2x_1+7x_2$ の大きい方を最小化せよ.

$$ x_1 + 2 x_2 \geq 12 \\ 2 x_2 + x_2 \geq 15 $$
model = Model("min-max")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
z = model.addVar(name="z")
model.update()
model.addConstr(x1 + 2 * x2 >= 12)
model.addConstr(2 * x1 + x2 >= 15)
model.addConstr(3 * x1 + 4 * x2 <= z)
model.addConstr(2 * x1 + 7 * x2 <= z)
model.setObjective(z, GRB.MINIMIZE)
model.optimize()
print(model.ObjVal)
print(x1.X, x2.X)
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 4 rows, 3 columns and 10 nonzeros
Model fingerprint: 0x3c2903f1
Coefficient statistics:
  Matrix range     [1e+00, 7e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [1e+01, 2e+01]
Presolve time: 0.00s
Presolved: 4 rows, 3 columns, 10 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    0.0000000e+00   1.350000e+01   0.000000e+00      0s
       4    3.1200000e+01   0.000000e+00   0.000000e+00      0s

Solved in 4 iterations and 0.01 seconds
Optimal objective  3.120000000e+01
31.200000000000003
7.2 2.4

最小値の最小化

最小値の最小化はさらにやっかいである. 最大値の最小化の場合と同様に,新しい実数変数 $z$ を導入し,今度は最小値と一致するように以下の制約を加える.

$$ 3x_1 + 4x_2 \geq z, \ \ 2x_1 + 7x_2 \geq z $$

これらの制約だけで $z$ を最小化すると目的関数値は $-\infty$ になる. これを避けるためには,上の制約において,左辺の小さい方と $z$ が等しくなるという条件が必要になる. 大きな数 $M$ と $0$-$1$ 変数 $y$ を用いると,この条件は以下の制約で記述することができる.

$$ 3x_1 + 4x_2 \leq z+ My, \ \ 2x_1 + 7x_2 \leq z + M(1-y) $$

$y=0$ のときには,左側の制約だけが意味をもち,$3x_1 + 4x_2 \leq z$ となる. これを元の制約 $3x_1 + 4x_2 \geq z$ と合わせることによって,$3x_1 + 4x_2 = z$ を得る. $y=1$ のときには,右側の制約だけが意味をもち, $2x_1 + 7x_2 \leq z$ となる. これを元の制約 $2x_1 + 7x_2 \geq z$ と合わせることによって $2x_1 + 7x_2 =z$ を得る.

例題: 以下の制約の下で, $3x_1+4x_2$ と $2x_1+7x_2$ の小さい方を最小化せよ.

$$ x_1 + 2 x_2 \geq 12 \\ 2 x_2 + x_2 \geq 15 $$
model = Model("min-min")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
z = model.addVar(name="z")
y = model.addVar(name="y", vtype="B")
M = 99999
model.update()
model.addConstr(x1 + 2 * x2 >= 12)
model.addConstr(2 * x1 + x2 >= 15)
model.addConstr(3 * x1 + 4 * x2 >= z)
model.addConstr(2 * x1 + 7 * x2 >= z)
model.addConstr(3 * x1 + 4 * x2 <= z + M * y)
model.addConstr(2 * x1 + 7 * x2 <= z + M * (1 - y))
model.setObjective(z, GRB.MINIMIZE)
model.optimize()
print(model.ObjVal)
print(x1.X, x2.X)
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 6 rows, 4 columns and 18 nonzeros
Model fingerprint: 0x80228eea
Variable types: 3 continuous, 1 integer (1 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+05]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+01, 1e+05]
Presolve time: 0.00s
Presolved: 6 rows, 4 columns, 18 nonzeros
Variable types: 3 continuous, 1 integer (1 binary)

Root relaxation: objective 0.000000e+00, 3 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0    0.00000    0    1          -    0.00000      -     -    0s
H    0     0                      30.0000000    0.00000   100%     -    0s
H    0     0                      24.0000000    0.00000   100%     -    0s
     0     0    0.00000    0    1   24.00000    0.00000   100%     -    0s
     0     2    0.00000    0    1   24.00000    0.00000   100%     -    0s

Explored 3 nodes (5 simplex iterations) in 0.01 seconds
Thread count was 16 (of 16 available processors)

Solution count 2: 24 30 

Optimal solution found (tolerance 1.00e-04)
Best objective 2.400000000000e+01, best bound 2.400000000000e+01, gap 0.0000%
24.0
12.0 0.0

例題 絶対値

実数変数 $x$ の絶対値 $|x|$ を「最小化」したいという要求は実務でしばしば現れる. 「最小化」なら,変数を1つ追加するだけで処理できる. まず,$|x|$ を表す新しい変数 $z$ を追加し, $z \geq x$ と $z \geq -x$ の2つの制約を追加する. $z$ を最小化すると,$x$ が非負のときには $z \geq x$ の制約が効いて,$x$ が負のときには $z \geq -x$ の制約が効いてくるので,$z$ は $|x|$ と一致する.

別の方法もある.この方法は,後述するように絶対値の「最大化」にも拡張できる. 2つの新しい非負の実数変数 $y$ と $z$ を導入する. まず,変数 $x$ を $x=y-z$ と2つの非負条件をもつ実数変数の差として記述する. $x$ の絶対値 $|x|$ を $y+z$ と記述すると, $y \geq 0, z \geq 0$ の制約が付加されているので,$x$ が負のときには $z$ が正(このとき $y$ は $0$ )となり, $x$ が正のときには $y$ が正(このとき $z$ は $0$)になる. つまり,定式化内の $x$ をすべて $y-z$ で置き換え,最小化したい目的関数内の $|x|$ をすべて $y+z$ で置き換えれば良い.

例題: 以下の線形最適化問題を,上で説明した2通りの定式化で解け.

$$ \begin{array}{ll} minimize & 2 x_1 + 3 x_2 + |x_1-x_2| \\ s.t. &x_1 + 2 x_2 \geq 12 \\ & 2 x_2 + x_2 \geq 15 \\ & x_1, x_2 \geq 0 \end{array} $$
model = Model("absolute(1)")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
A = model.addVar(name="A", lb=-1000000000.0)  # 差を表す実数変数
z = model.addVar(name="z")  # |A|を表す補助変数
model.update()
model.addConstr(x1 + 2 * x2 >= 12)
model.addConstr(2 * x1 + x2 >= 15)
model.addConstr(x1 - x2 == A)
model.addConstr(z >= A)
model.addConstr(z >= -A)
model.setObjective(2 * x1 + 3 * x2 + z, GRB.MINIMIZE)
model.optimize()
print(model.ObjVal)
print(x1.X, x2.X, A.X)
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 5 rows, 4 columns and 11 nonzeros
Model fingerprint: 0x6de7dfc6
Coefficient statistics:
  Matrix range     [1e+00, 2e+00]
  Objective range  [1e+00, 3e+00]
  Bounds range     [1e+09, 1e+09]
  RHS range        [1e+01, 2e+01]
Presolve time: 0.00s
Presolved: 5 rows, 4 columns, 11 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    0.0000000e+00   1.350000e+01   0.000000e+00      0s
       3    2.4000000e+01   0.000000e+00   0.000000e+00      0s

Solved in 3 iterations and 0.01 seconds
Optimal objective  2.400000000e+01
24.0
6.0 3.0 3.0
model = Model("absolute(2)")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
y = model.addVar(name="y")  # x_1-x_2 = y-z
z = model.addVar(name="z")
model.update()
model.addConstr(x1 + 2 * x2 >= 12)
model.addConstr(2 * x1 + x2 >= 15)
model.addConstr(x1 - x2 == y - z)
model.setObjective(2 * x1 + 3 * x2 + y + z, GRB.MINIMIZE)
model.optimize()
print(model.ObjVal)
print(x1.X, x2.X, y.X, z.X)
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 3 rows, 4 columns and 8 nonzeros
Model fingerprint: 0xd0287f74
Coefficient statistics:
  Matrix range     [1e+00, 2e+00]
  Objective range  [1e+00, 3e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [1e+01, 2e+01]
Presolve time: 0.00s
Presolved: 3 rows, 4 columns, 8 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    0.0000000e+00   1.350000e+01   0.000000e+00      0s
       3    2.4000000e+01   0.000000e+00   0.000000e+00      0s

Solved in 3 iterations and 0.00 seconds
Optimal objective  2.400000000e+01
24.0
6.0 3.0 3.0 0.0

絶対値の最大化

絶対値を「最大化」したい場合には,$y$ と $z$ のいずれか一方だけを正にすることができるという制約が必要になる. そのためには,大きな数 $M$ と $0$-$1$ 変数 $\xi$ を使い,以下の制約を付加する. $$ y \leq M \xi $$ $$ z \leq M (1-\xi) $$ 変数 $\xi$ が $1$ のときには,$y$ が正になることができ,$0$ のときには $z$ が正になることができる.

例題: 以下の線形最適化問題を解け.

$$ \begin{array}{ll} maximize & -2 x_1 - 3 x_2 + |x_1-x_2| \\ s.t. &x_1 + 2 x_2 \geq 12 \\ & 2 x_2 + x_2 \geq 15 \\ & x_1, x_2 \geq 0 \end{array} $$
model = Model("absolute(3)")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
y = model.addVar(name="y")  # |x_1-x_2| = y-z
z = model.addVar(name="z")
xi = model.addVar(name="xi", vtype="B")
M = 999999.0
model.update()
model.addConstr(x1 + 2 * x2 >= 12)
model.addConstr(2 * x1 + x2 >= 15)
model.addConstr(x1 - x2 == y - z)
model.addConstr(y <= M * xi)
model.addConstr(z <= M * (1 - xi))
model.setObjective(-2 * x1 - 3 * x2 + y + z, GRB.MAXIMIZE)
model.optimize()
print(model.ObjVal)
print(x1.X, x2.X, y.X, z.X)
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 5 rows, 5 columns and 12 nonzeros
Model fingerprint: 0x6f41ab3b
Variable types: 4 continuous, 1 integer (1 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+06]
  Objective range  [1e+00, 3e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+01, 1e+06]
Presolve time: 0.00s
Presolved: 5 rows, 5 columns, 12 nonzeros
Variable types: 4 continuous, 1 integer (1 binary)

Root relaxation: objective 9.999780e+05, 5 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 999978.000    0    1          - 999978.000      -     -    0s
H    0     0                     -25.0000000 999978.000      -     -    0s
H    0     0                     -12.0000000 999978.000      -     -    0s
     0     0     cutoff    0       -12.00000  -12.00000  0.00%     -    0s

Explored 1 nodes (6 simplex iterations) in 0.01 seconds
Thread count was 16 (of 16 available processors)

Solution count 2: -12 -25 
No other solutions better than -12

Optimal solution found (tolerance 1.00e-04)
Best objective -1.200000000000e+01, best bound -1.200000000000e+01, gap 0.0000%
-12.0
12.0 0.0 12.0 0.0

問題 スーパーの位置(絶対値)

2 次元の格子上にある4 つの点 $(0,1), (2,0), (3,3), (1,3)$ に住んでいる人たちが, 最もアクセスが良くなる地点にスーパーを配置しようと考えている. 格子上での移動時間が,$x$ 座標の差と $y$ 座標の差の和で与えられるとしたとき,最適なスーパーの位置を求めよ.

ヒント: スーパーの位置を $X,Y$ としたとき,(たとえば)点 $(1,3)$ からの距離は, $|X-1|+|Y-3|$ と絶対値を用いて計算できる.

問題 消防署の位置(最大値の最小化)

上と同じ地点に住む人たちが,今度は消防署を作ろうと考えている. 最も消防署から遠い地点に住む人への移動距離を最小にするには,どこに消防署を配置すれば良いだろうか?

例題 離接制約

either or の条件 : 何れか1つが満たされていなければならないという条件

  • $2x_1 + x_2 \leq 30, x_1 + 2x_2 \leq 30$ の2 本の制約の何れかが成立するという条件

    • 大きな数 $M$ と $0$-$1$変数 $y$ を用いて,

$$2x_1 + x_2 \leq 30 +My, \ \ \ x_1 + 2x_2 \leq 30+M(1-y)$$

  • $2x_1 + x_2 \leq 30, \ \ \ x_1 + 2x_2 \leq 30, \ \ \ 5x_1+x_2 \leq 50$ の3 本の制約の何れかが成立するという条件

    • 3 つの$0$-$1$変数 $y_1, y_2, y_3$ を用いて,

$$2x_1 + x_2 \leq 30 +M(1-y_1), \ \ \ x_1 + 2x_2 \leq 30+M(1-y_2), \ \ \ 5x_1+x_2 \leq 50+M(1-y_3)$$

$$y_1 +y_2 + y_3 \geq 1$$

例題: 以下の線形最適化問題を解け.

$$ \begin{array}{ll} maximize & x_1 + x_2 \\ s.t. & 2x_1 + x_2 \leq 30  もしくは  x_1 + 2x_2 \leq 40 \\ & x_1, x_2 \geq 0 \end{array} $$
model = Model("disjunctive")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
y = model.addVar(name="y", vtype="B")
M = 99999
model.update()
model.addConstr(2 * x1 + x2 <= 30 + M * y)
model.addConstr(x1 + 2 * x2 <= 40 + M * (1 - y))
model.setObjective(x1 + x2, GRB.MAXIMIZE)
model.optimize()
print(model.ObjVal)
print(x1.X, x2.X)
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 2 rows, 3 columns and 6 nonzeros
Model fingerprint: 0x4e541ef5
Variable types: 2 continuous, 1 integer (1 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+05]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [3e+01, 1e+05]
Found heuristic solution: objective 15.0000000
Found heuristic solution: objective 40.0000000
Presolve time: 0.00s
Presolved: 2 rows, 3 columns, 6 nonzeros
Variable types: 2 continuous, 1 integer (1 binary)

Root relaxation: objective 4.400000e+01, 2 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0     cutoff    0        40.00000   40.00000  0.00%     -    0s

Explored 0 nodes (2 simplex iterations) in 0.01 seconds
Thread count was 16 (of 16 available processors)

Solution count 2: 40 15 

Optimal solution found (tolerance 1.00e-04)
Best objective 4.000000000000e+01, best bound 4.000000000000e+01, gap 0.0000%
40.0
40.0 0.0

例題 if A then B 条件

「制約Aが成立している場合には,制約Bも成立しなければならない」という条件

実際の問題を考える際には、論理制約がしばしば出てくる。例えば、「倉庫を建設しなければそこに保管することができない」とか「この作業とこの作業は同時刻に処理することはできない」などが代表的だ。ここでは、その取り扱いについて解説する。

事象が真か嘘かは0-1変数で表現できる。0-1変数x,yに対して、「xが真ならばyも真である」という論理制約は、x<=y で表現できる。何れかが真であるという論理制約は x+y=1 となる。 もう少し難しいものとして、「「この条件が成立しているときには、この制約が必要だ」(if ... then ... )という論理制約がある。 「制約Aが成立している場合には,制約Bも成立しなければならない」という条件は, 「Aが成立しない or Bが成立する」と同値である。 0−1変数x,yのいずれかが1になるという制約は x+y=1 であるので、「NOT A or B」は0−1変数を用いて記述することができる。

「if A then B」は「NOT A or B」と同値

A B if A then B NOT A NOT A or B

例として,整数変数 $x_1,x_2$ に対する以下の論理制約を考えよう。

  • if.. then の条件:

$$if \ \ \ x_1 +x_2 \leq 10 \ \ \ then \ \ \ x_1 \leq 5 $$

この条件は,

$$x_1+x_2 \geq 11 \ \ \ or \ \ \ x_1 \leq 5 $$

と同値であるので, 0-1変数 $y$ と大きな数 $M$ を用いて

  • 離接制約で表現:
$$x_1 + x_2 > 10 - My, \ \ \ x_1 \leq 5 + M(1-y)$$

と書くことができる. $x_1, x_2$ が整数なので, 第1式は以下と同値

$$x_1 + x_2 \geq 11 - My$$

  • $x_1,x_2$ が整数でなく実数変数の場合の問題点:

数理最適化ソルバーでは,より大きい($>$)と以上($\geq$)の制約を区別しない

  • $x_1 + x_2 > 10$ は $x_1 + x_2 \geq 10$ と同値
  • 微少な値 $\epsilon >0$ を用いて,
$$x_1 + x_2 \geq 10 +\epsilon - My$$

例題: 以下の論理条件を含む整数最適化問題を解け.

$$ \begin{array}{ll} minimize & x_1 + x_2 \\ s.t. & x_2 \geq 2 \\ & if \ \ \ x_1 \leq 2 \ \ \ then \ \ \ x_2 \geq 6 \\ & x_1, x_2 は非負の整数 \end{array} $$
model = Model("if_then")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
y = model.addVar(name="y", vtype="B")
M = 99999
model.update()
model.addConstr(x2 >= 2)
model.addConstr(x1 >= 3 - M * y)
model.addConstr(x2 >= 6 - M * (1 - y))
model.setObjective(x1 + x2, GRB.MINIMIZE)
model.optimize()
print(model.ObjVal)
print(x1.X, x2.X)
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 3 rows, 3 columns and 5 nonzeros
Model fingerprint: 0xb7b497a0
Variable types: 2 continuous, 1 integer (1 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+05]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [2e+00, 1e+05]
Found heuristic solution: objective 5.0000000
Presolve removed 3 rows and 3 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.00 seconds
Thread count was 1 (of 16 available processors)

Solution count 1: 5 

Optimal solution found (tolerance 1.00e-04)
Best objective 5.000000000000e+00, best bound 5.000000000000e+00, gap 0.0000%
5.0
3.0 2.0

例題 嘘つき島パズル

ある島には正直族と嘘つき族と呼ばれる2 種類の人たちが仲良く住んでいる. 正直族は必ず本当のことを言い,嘘つき族は必ず嘘をつく.

あなたは,この島の人たちが正直族か嘘つき族なのかの調査を依頼された.

  1. 最初の家の旦那に聞いたところ「夫婦は両方とも嘘つき族だよ」という答えだった.

  2. 次の家に言って旦那に「ご夫婦は両方とも嘘つき族ですか?」と聞いたところ「少なくとも1人はね」という答えだった.

  3. その次の家に行ったところ旦那が答えた「もし私が正直族なら,私の家内も正直族です」

さて,上の情報から,調査結果を報告せよ.

# 1. x=1 <=> x+y=0
model = Model()
x = model.addVar(name="x", vtype="B")
y = model.addVar(name="y", vtype="B")
z = model.addVar(name="z", vtype="B")
# if x=1, then x+y=0
# x=0 or x+y =0
model.addConstr(x + (x + y) <= 1)
# if x+y=0, then x=1
# x+y>=1 or x>=1
model.addConstr(x + y >= 1 - 9999 * z)
model.addConstr(x >= 1 - 9999 * (1 - z))
model.setObjective(x, GRB.MAXIMIZE)
model.optimize()
print(x.X, y.X)  # 旦那は嘘つき,奥さんは正直
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 3 rows, 3 columns and 7 nonzeros
Model fingerprint: 0xd3674ad5
Variable types: 0 continuous, 3 integer (3 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+04]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+04]
Presolve removed 3 rows and 3 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.00 seconds
Thread count was 1 (of 16 available processors)

Solution count 1: -0 
No other solutions better than -0

Optimal solution found (tolerance 1.00e-04)
Best objective -0.000000000000e+00, best bound -0.000000000000e+00, gap 0.0000%
0.0 1.0
# 2. x=1 <=> x+y<=1
model = Model()
x = model.addVar(name="x", vtype="B")
y = model.addVar(name="y", vtype="B")
z1 = model.addVar(name="z1", vtype="B")
z2 = model.addVar(name="z2", vtype="B")
# if x=1, then x+y<=1
# x<=0 or x+y <=1
model.addConstr(x <= 9999 * z1)
model.addConstr(x + y <= 1 + 9999 * (1 - z1))
# if x+y<=1, then x=1
# x+y>=2 or x>=1
model.addConstr(x + y >= 2 - 9999 * z2)
model.addConstr(x >= 1 - 9999 * (1 - z2))
model.setObjective(x, GRB.MAXIMIZE)
model.optimize()
print(x.X, y.X)  # 旦那は正直,奥さんは嘘つき
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 4 rows, 4 columns and 10 nonzeros
Model fingerprint: 0x47462c39
Variable types: 0 continuous, 4 integer (4 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+04]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [2e+00, 1e+04]
Presolve removed 4 rows and 4 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.00 seconds
Thread count was 1 (of 16 available processors)

Solution count 1: 1 

Optimal solution found (tolerance 1.00e-04)
Best objective 1.000000000000e+00, best bound 1.000000000000e+00, gap 0.0000%
1.0 0.0
# 3. x=1 <=> (if x=1, then y=1) => x<=y
model = Model()
x = model.addVar(name="x", vtype="B")
y = model.addVar(name="y", vtype="B")
z = model.addVar(name="z", vtype="B")
# if x=1, then x<=y
model.addConstr(x <= y)
# if x<=y, then x=1
# x>=y+1 or x>=1
model.addConstr(x >= y + 1 - 9999 * z)
model.addConstr(x >= 1 - 9999 * (1 - z))
model.setObjective(x, GRB.MAXIMIZE)
model.optimize()
print(x.X, y.X)  # 両方とも正直
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 3 rows, 3 columns and 7 nonzeros
Model fingerprint: 0x5cd4f006
Variable types: 0 continuous, 3 integer (3 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+04]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+04]
Found heuristic solution: objective 1.0000000

Explored 0 nodes (0 simplex iterations) in 0.00 seconds
Thread count was 1 (of 16 available processors)

Solution count 1: 1 

Optimal solution found (tolerance 1.00e-04)
Best objective 1.000000000000e+00, best bound 1.000000000000e+00, gap 0.0000%
1.0 1.0

問題 論理パズル (1)

ある島には正直族と嘘つき族と呼ばれる2 種類の人たちが仲良く住んでいる. 正直族は必ず本当のことを言い,嘘つき族は必ず嘘をつく.

またこの島には,夜になると狼に変身して人を襲う狼男が紛れ込んでいる. 狼男もこの島の住民なので,正直族か嘘つき族の何れかに属する. あなたは,この島の人たちが狼男なのかの調査を依頼された.

3人のうち1 人が狼男であることが分かっている. A,B,C の3人組への証言は以下の通りである.

  • A: 「わたしは狼男です.」
  • B: 「わたしも狼男です.」
  • C: 「わたしたちの中の高々1人が正直族です.」

さて,狼男は誰か?また誰が正直族で誰が嘘つき族か?

問題 論理パズル (2)

同じ島でまた別の3 人組 A,B,C の証言を得た.

  • A: 「わたしたちの中の少なくとも1 人が嘘つき族です.」
  • B: 「Cさんは正直族です.」

この3 人組も彼らのうちの1人が狼男で,彼は正直族であることが分かっているとき,狼男は誰か?また誰が正直族で誰が嘘つき族か?

問題 絶対値

7 組の家族がみんなで食事会をしようと考えている. 4 人がけのテーブルが7 卓であり,親睦を深めるために,同じ家族の人たちは同じテーブルに座らないようにしたい.

家族の構成は以下の通りとしたとき,各テーブルの男女比をなるべく均等にする(女性と平均人数との差の絶対値の和を最小化する)座り方を求めよ. ただし,女性には (F) の記号を付してある.

  • 磯野家: 波平,フネ (F),カツオ,ワカメ(F)
  • バカボン家: バカボンパパ,バカボンママ (F),バカボン,ハジメ
  • 野原家: ひろし,みさえ(F),しんのすけ,ひまわり(F)
  • のび家: のび助,玉子(F),のび太
  • 星家: 一徹,明子(F),飛雄馬
  • レイ家: テム,カマリア(F),アムロ
  • ザビ家: デギン,ナルス(F),ギレン,キシリア(F),サスロ,ドズル,ガルマ

    また、男女比をなるべく均等にしない(女性と平均人数との差の絶対値の和を最大化する)座り方を求めよ。

例題 多期間生産計画

あなたの会社では,2つの製品(P1, P2)を2台の機械(M1, M2)で製造している. 今後3ヶ月の需要が決まっており,製品P1に対しては 4000,8000,3000, 製品P2に対しては 1000,5000,5000である.

各月の機械を使用できる時間(h)も決まっており,機械M1に対しては700,300,1000,機械M2に対しては1500,400,300である.

製品は機械によって製造速度が異なり, 製品を1単位製造するのに必要な時間(h)は,以下の表のように与えられているものとする.

M1    M2
P1  0.15   0.16
P2  0.12    0.14

製造の変動費用は,製品や機械によらず1時間あたり5ドル,製品1単位を1ヶ月持ち越したときの生じる在庫費用は0.1ドルとする. なお,初期在庫は0,最終在庫も0と仮定する.

  1. 最適な生産計画を立案せよ.

  2. 最終期に機械M1の使用可能な時間を100時間増やせる設備を100ドルで購入可能である. 購入すべきか?

  3. 機械のメンテナンスをするなら,どの月にどの機械に対してやるべきか?

from mypulp import Model, GRB, quicksum

T = 3
products = ["P1", "P2"]
machines = ["M1", "M2"]
demand = {"P1": [4000, 8000, 3000], "P2": [1000, 5000, 5000]}
capacity = {"M1": [700, 300, 1000], "M2": [1500, 400, 300]}
rate = {("P1", "M1"): 0.15, ("P1", "M2"): 0.16, ("P2", "M1"): 0.12, ("P2", "M2"): 0.14}
variable_cost = 5.0
inventory_cost = 0.1

model = Model()
x, I = {}, {}
for p in products:
    for m in machines:
        for t in range(T):
            x[p, m, t] = model.addVar(name=f"x({p},{m},{t})")
for p in products:
    I[p, -1] = I[p, T - 1] = 0.0
    for t in range(T - 1):
        I[p, t] = model.addVar(name=f"I({p},{t})")
model.update()

capacity_constraints = {}
for m in machines:
    for t in range(T):
        capacity_constraints[m, t] = model.addConstr(
            quicksum(rate[p, m] * x[p, m, t] for p in products) <= capacity[m][t],
            name=f"Capacity Constraint({m},{t})",
        )

demand_constraints = {}
for p in products:
    for t in range(T):
        demand_constraints[p, t] = model.addConstr(
            I[p, t - 1] + quicksum(x[p, m, t] for m in machines)
            == demand[p][t] + I[p, t]
        )

model.setObjective(
    quicksum(
        variable_cost * rate[p, m] * x[p, m, t]
        for p in products
        for m in machines
        for t in range(T)
    )
    + quicksum(inventory_cost * I[p, t] for p in products for t in range(T)),
    GRB.MINIMIZE,
)
model.optimize()
print("Obj. Val.=", model.ObjVal)
for t in range(T):
    for p in products:
        for m in machines:
            print(p, m, t, x[p, m, t].X)
Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - /Users/mikiokubo/Library/Caches/pypoetry/virtualenvs/analytics-v-sH3Dza-py3.8/lib/python3.8/site-packages/pulp/apis/../solverdir/cbc/osx/64/cbc /var/folders/69/5y96sdc94jxf6khgc8mlmxrr0000gn/T/a30e33bff56c4f9a87928084e1af31a1-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/69/5y96sdc94jxf6khgc8mlmxrr0000gn/T/a30e33bff56c4f9a87928084e1af31a1-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 17 COLUMNS
At line 66 RHS
At line 79 BOUNDS
At line 96 ENDATA
Problem MODEL has 12 rows, 16 columns and 32 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Presolve 12 (0) rows, 16 (0) columns and 32 (0) elements
0  Obj 0 Primal inf 26000 (6)
12  Obj 19173.333
Optimal - objective value 19173.333
Optimal objective 19173.33333 - 12 iterations time 0.002
Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

Obj. Val.= 19173.333354000002
P1 M1 0 1866.6667
P1 M2 0 7633.3333
P2 M1 0 3500.0
P2 M2 0 0.0
P1 M1 1 0.0
P1 M2 1 2500.0
P2 M1 1 2500.0
P2 M2 1 0.0
P1 M1 2 2666.6667
P1 M2 2 333.33333
P2 M1 2 5000.0
P2 M2 2 0.0

余裕変数(slack)や双対変数(Pi)をみると,機械M1の1時間あたりの価値は0.33ドルなので,購入すべきでない.

余裕変数をみるとメインテナンスは,最初の月か最後の月にすべきである.

for c in capacity_constraints:
    print(c, capacity_constraints[c].slack, capacity_constraints[c].Pi)
('M1', 0) -0.0 -0.33333333
('M1', 1) -0.0 -1.1666667
('M1', 2) -0.0 -0.33333333
('M2', 0) 278.6667 0.0
('M2', 1) -0.0 -0.625
('M2', 2) 246.666667 0.0

問題 2段階製造問題

ある会社が,3つの製品A,B,Cを製造している。 製造工程には、「切断」と「プレス」という2つの工程がある. すべての製品はこの2つの瞬間を通過しなければならない。

1日8時間使用可能な切削部門の生産能力は,製品Aは1時間に2000個、製品Bは1時間に1600個,製品Cは1時間に1100個である.

1日8時間使用可能なプレス部門の生産能力は,製品Aは1時間に1000個、製品Bは1時間に1500個、製品Cは1時間に2400個である.

製品切り替えの段取り時間(費用)はかからないものとする.

製品A,B,Cの1個あたりの利益をそれぞれ12,9,8(万円)としたとき,利益を最大化する製造計画を立案せよ.

問題 多期間ロットサイズ決定問題

1つの製品の生産をしている工場を考える. 在庫費用は1日あたり,1トンあたり1万円とする. いま7日先までの需要が分かっていて,7日分の生産量と在庫量を決定したい. 各日の需要は,以下の表のようになっている. 工場の1日の稼働時間は $8$時間,製品1トンあたりの製造時間は$1$時間としたとき, 稼働時間上限を満たした最小費用の生産・在庫量を決定せよ.

日  1   2   3  4  5   6   7 
需要量  5   7 8 2 9 1 3

また,生産をした日には,段取り費用 $10$ 万円がかかると仮定したときの生産・在庫量を求めよ.

問題 ナプキンのクリーニング(難)

あなたはホテルの宴会係だ.あなたは1週間に使用するナプキンを手配する必要がある. 各日の綺麗なナプキンの需要量は平日は $100$枚,土曜日と日曜日は $125$枚だ. 新しいナプキンを購入するには $100$円かかる. 使用したナプキンはクリーニング店で洗濯して綺麗なナプキンにすることができるが, 早いクリーニング店だと1日で1枚あたり $30$円かかり, 遅いクリーニング店だと2日で1枚あたり $10$円かかる. 月曜の朝のナプキンの在庫が $0$ としたとき,需要を満たす最適なナプキンの購入ならびにクリーニング計画をたてよ.

さらに,日曜末の在庫を月曜の朝に使うことができると仮定したときの最適なナプキンのクリーニング計画をたてよ.