AMPLによる最適化問題のモデリング

線形最適化と混合整数最適化と非線形最適化のモデリングのコツについて述べる．

数理最適化とモデリング言語 AMPL

ここでは，AMPLを用いて数理最適化のモデリングの基本について学ぶ．

  • 線形最適化
    • 線形最適化
    • 簡単な例題と双対問題
    • 栄養問題
  • 混合整数最適化問題
    • 簡単な例題（パズル）
    • 多制約ナップサック問題
  • モデリングのコツ
    • 最大値の最小化
    • 絶対値
    • 離接制約
    • if A then B 条件
from amplpy import AMPL, Environment, tools
ampl = AMPL(Environment("/Users/mikiokubo/Documents/ampl/"))

Google Colabの場合には，以下のセルを実行して amplpy とAMPLをインストールする （途中でUUIDを入れるように促されるので，コピーする）．

#!pip install amplpy
# from amplpy import ampl_notebook
# ampl = ampl_notebook(
#     modules=["highs","gurobi","cbc","scip","coin","gecode"],  #coin includes ipopt, couenne, bonmin
#     license_uuid=None)

線形最適化問題

線形最適化問題の一般形

線形最適化問題の一般形は以下のように書ける．

\[ \begin{array}{ll} minimize & c^T x \\ s.t. & A x \leq b \\ & x \in \mathbf{R}^n \end{array} \]

ここで，\(x\) は変数を表す \(n\)次元実数ベクトル， \(A\) は制約の左辺係数を表す \(m \times n\) 行列，\(c\) は目的関数の係数を表す \(n\)次元ベクトル， \(b\) は制約式の右辺定数を表す \(m\)次元ベクトルである．

実際問題を解く際の注意

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

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

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

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

簡単な例題

簡単な例題で線形最適化問題を説明する．

トンコケ，コケトン，ミックスの丼を販売している． 各丼の販売数を変数とし， トンコケ丼 \(x_1\)，コケトン丼 \(x_2\) ，ミックス丼 \(x_3\) とする． 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} \]

これをAMPLを用いて解く．

マジックコマンド%%ampl_evalを入れると、そのままAMPLのモデルとコマンドが記述できる。 （記述したamplファイルは、Pythonのamplインスタンスに保管される。）

amplインスタンスはグローバルに保管されるので、最初にAMPLのresetコマンドで初期化をしておく（最初の場合は省略可）。

reset;

var x1 >=0;
var x2 >=0;
var x3 >=0, <=30;

maximize obj: 15 * x1 + 18 * x2 + 30 * x3;
con1: 2 * x1 + x2 + x3 <= 60;
con2: x1 + 2 * x2 + x3 <= 60;

option solver gurobi; #ソルバーを指定する。 scip, highsでも可
solve; 

display x1, x2, x3; #displayコマンドで変数の値を得ることができる
Gurobi 12.0.1: optimal solution; objective 1230
3 simplex iterations
x1 = 10
x2 = 10
x3 = 30

モデルの確認

  • expandコマンドで式を展開できる。式（目的関数）につけた名前をつけると、一部だけ展開できる。
  • showコマンドでモデルに含まれる要素を確認できる。
expand con1;
subject to con1:
    2*x1 + x2 + x3 <= 60;

show;

variables:   x1   x2   x3

constraints:   con1   con2

objective:   obj

例題 双対問題

上の線形最適化問題において，資源（豚肉，鶏肉，牛肉）の価値を推定したい．

これは，先程の定式化を主問題としたときの双対問題を解けば良い．

主問題 \[ \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} \]

双対問題

\[ \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} \]

実は，主問題を解けば双対問題の最適解も同時に得ることができる．

AMPLでは、制約名（もしくは制約名.dual）を表示すると最適双対変数が得られる。 変数 \(x_3\) に対しては上限制約 \(x_3 \leq 30\) を明示的に入れていないので、 変数 \(x_3\) の被約費用 (reduced cost) を（変数名.rcで）みることによって制約の双対変数を得ることができる。

display con1, con2.dual, x3.rc;
con1 = 4
con2.dual = 7
x3.rc = 19

一般的な記述法

実際の定式化では，多くの変数と制約を記述する必要がある．

たとえば，以下のような制約を定義したいとしよう．

\(a = [5,4,2]\) に対して \[ \sum_{i=1}^3 a_i x_ i \]

これはsumを用いて以下のように記述できる．

最初にmodelでモデルを定義する。インデックス（添字）は集合idxとして定義しておく。 パラメータ \(a\) と変数 \(x\) は、どちらもインデックスidxを添え字としてもつ。 パラメータ \(a\) の実際の値は、データ記述部（dataの後）で定義される。

reset;

model;
set idx := {1,2,3};
param a{idx};
var x{idx} >= 0;
con1: sum{i in idx} a[i] * x[i] =10;

data;
param a:= 1 5 
          2 4
          3 2;

expand con1;
subject to con1:
    5*x[1] + 4*x[2] + 2*x[3] = 10;

最初の例題を一般的に記述してみよう。

amplではモデルとデータを分離することが推奨されている。

reset;

model;

set Bowls;
set Meats;
param price{Bowls};
param ub{Meats};
param usage{Meats, Bowls};
var x{Bowls} >=0;

maximize obj: sum{j in Bowls}  price[j]*x[j];
subject to resource{i in Meats}: sum{j in Bowls} usage[i,j]*x[j] <= ub[i];

data;

set Bowls := Tonkoke Koketon Mix;
set Meats := Pork  Chiken Beef;
param price :=  Tonkoke 15 Koketon 18  Mix 30; 
param ub := Pork 60 Chiken 60 Beef 30;
param usage:   Tonkoke Koketon Mix := 
       Pork       2      1      1
       Chiken     1      2      1  
       Beef       0      0      1    
     ;

option solver gurobi; 
solve; 

display x;
Gurobi 12.0.1: optimal solution; objective 1230
3 simplex iterations
x [*] :=
Koketon  10
    Mix  30
Tonkoke  10
;

例題：栄養問題

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

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

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

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

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

栄養素 \(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 35 660 6.0

商品の集合を \(F\) （Foodの略），栄養素の集合を \(N\)（Nutrientの略）とする． 栄養素 \(i\) の1日の摂取量の下限を \(a_i\)， 上限を \(b_i\) とし， 商品 \(j\) の価格を \(c_j\)，含んでいる栄養素 \(i\) の量を \(n_{ij}\) とする． 商品 \(j\) を購入する個数を非負の実数変数 \(x_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} \]

モデルファイルをファイルに書き出す

マジックコマンド %%writefile を用いて，栄養問題のモデルファイル diet.mod を保存する．

set NUTRITIONS;  # 栄養素 (Cal, Carbo, Protein, VitA, VitC, Calc, Iron)
set FOODS;       # 食品 (CQPounder, BigM, FFilet, Chicken, Fries, Milk, VegJuice)

param lower{NUTRITIONS};  # 栄養素の下限 (b_i)
param upper{NUTRITIONS};  # 栄養素の上限 (a_i)
param nutrition_content{FOODS, NUTRITIONS};  # 栄養素の含有量 (n_{ij})
param cost{FOODS};        # 食品の価格 (c_j)

var x{FOODS} >= 0;  # 各食品の購入量（非負連続変数）

minimize Total_Cost: sum{j in FOODS} cost[j] * x[j];

subject to Nutrition_Bounds{i in NUTRITIONS}:
    lower[i] <= sum{j in FOODS} nutrition_content[j,i] * x[j]  <= upper[i];
Overwriting diet.mod

データファイルを書き出す

マジックコマンド %%writefile diet.dat でデータファイルを書き出す。

data;

# 栄養素と食品の定義
set NUTRITIONS := Cal Carbo Protein VitA VitC Calc Iron;
set FOODS := CQPounder BigM FFilet Chicken Fries Milk VegJuice;

# 栄養素の下限・上限
param:  lower upper :=
    Cal     2000      3000
    Carbo   300       375
    Protein 50        60
    VitA    500       750
    VitC    35        100
    Calc    660       900
    Iron    6.0       7.5;

param nutrition_content : Cal Carbo Protein VitA VitC Calc Iron :=
    CQPounder  556  39    30      147  10   221  2.4
    BigM       556  46    26      97   9    142  2.4
    FFilet     356  42    14      28   1    76   0.7
    Chicken    431  45    20      9    2    37   0.9
    Fries      249  30    3       0    5    7    0.6
    Milk       138  10    7       80   2    227  0
    VegJuice   69   17    1       750  2    18   0;

param cost :=
    CQPounder  360
    BigM       320
    FFilet     270
    Chicken    290
    Fries      190
    Milk       170
    VegJuice   100;
Overwriting diet.dat

コマンドで最適化をする

マジックコマンド %%ampl_eval でコマンドを実行する． ソルバーはオープンソースのcbcを指定し，solveコマンドで求解する．その後，displayコマンドで変数Buyを表示する．

reset; 
model diet.mod;
data diet.dat;
option solver cbc;
solve;
display x;
cbc 2.10.12: optimal solution; objective 2133.04522
0 simplex iterations
x [*] :=
     BigM  0.506981
CQPounder  0
  Chicken  0
   FFilet  0.735637
    Fries  7.11383
     Milk  2.07028
 VegJuice  0.686137
;

問題 線形最適化(1）

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

ジュースAを作るにはぶどう \(3\) \(\ell\)とりんご \(1\) \(\ell\)が必要で， ジュースBを作るにはぶどう \(2\) \(\ell\)とりんご \(2\) \(\ell\)が必要である． （なんとジュースは \(4\) \(\ell\)入りの特大サイズなのだ！）

ジュース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\)mで，走る速度は分速 \(180\)mです．」

問題 線形最適化(4）

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

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

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

また、上の問題において， 豚肉，鶏肉，牛肉の \(100\)gあたりの価値はいくらになるか計算せよ．

問題 輸送問題と双対性

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

\[ \begin{array}{c c c c c c c c } \hline 工場生産可能容量 & & & 顧需要客 & & & & \\ \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 & 工場 & & & 顧客 & & & & \\ \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 & 左 & 右 \\ \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)

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

混合整数最適化問題

（混合）整数最適化問題の一般形

一般の整数線形最適化問題は以下のように書ける．

\[ \begin{array}{ll} minimize & c^T x \\ s.t. & A x \leq b \\ & x \in \mathbf{Z}^n \end{array} \]

ここで，\(x\) は変数を表す \(n\)次元ベクトル，\(A\) は制約の左辺係数を表す \(m \times n\) 行列， \(c\) は目的関数の係数を表す \(n\)次元ベクトル， \(b\) は制約式の右辺定数を表す \(m\)次元ベクトルである．

変数の一部が整数でなく実数であることを許した問題を混合整数最適化問題とよぶ．

例題： 鶴亀算の拡張

簡単な例題で整数線形最適化問題を説明する．

宇宙ステーションに、地球人（頭1つで足２本）、犬星人（頭1つで足４本）、タコ星人（頭１つで足８本）が何人かいる。

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

整数変数はintegerと宣言する。

reset;
var x integer >=0;  
var y integer >=0;
var z integer >=0;

minimize obj: y + z;
head: x + y+ z = 32;
leg: 2*x + 4*y + 8*z = 80;

option solver gurobi;
solve;
display x,y,z;
Gurobi 12.0.1: optimal solution; objective 4
0 simplex iterations
x = 28
y = 2
z = 2

#  hide
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, 5, 100
  • 大クマ: 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} \]

\(0\)-\(1\)整数変数（2値変数）はbinaryと宣言する。

reset;
model;
set I;
set J;
var x{J} binary;
param v{J};
param a{I,J};
param b{I};

maximize obj: sum{j in J} v[j]*x[j];

s.t. constraint{i in I}:
    sum{j in J} a[i,j]*x[j] <= b[i];

data;

set I := 1  2;
set J := 1  2  3  4;
param v := 1  16  2  19  3  23  4  28;

param a   : 1      2    3    4 :=
      1     2      3    4    5
      2    3000 3500 5100 7200;

param b := 1  7  2  10000;

option solver gurobi;
solve;
display x;
Gurobi 12.0.1: optimal solution; objective 42
0 simplex iterations
x [*] :=
1  0
2  1
3  1
4  0
;

例題：論理条件

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

  1. 小クマと中クマは仲が悪いので，同時にもって逃げてはいけない．
  2. 小クマは極小クマが好きなので，小クマを持って逃げるときには必ず極小クマももって逃げなければならない．
  3. 極小クマ，小クマ，大クマのうち，少なくとも2つは持って逃げなければならない．
# いずれかのコメントを外す
# 1.
add1:  x[2] + x[3] <=1;
# 2.
#add2:  x[2] <= x[1];
# 3.
#add3: x[1] + x[2] + x[4] >= 2;

option solver gurobi;
solve;

display x;
Gurobi 12.0.1: optimal solution; objective 39
0 simplex iterations
x [*] :=
1  1
2  0
3  1
4  0
;

問題 鶴亀蛸キメラ算

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

最大値の最小化

実数変数 \(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_1 + x_2 \geq 15 \]

reset;
set idx := {1,2};
var x{idx} >=0;
minimize obj: max( 3*x[1] + 4*x[2], 2*x[1]+7*x[2]);
con1: x[1]+2*x[2] >= 12;
con2: 2*x[1]+x[2] >= 15;

option solver gurobi;
solve;
display x;
Gurobi 12.0.1: optimal solution; objective 31.2
2 simplex iterations
1 branching node
x [*] :=
1  7.2
2  2.4
;

#  hide
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_1 + x_2 \geq 15 \]

reset;
set idx := {1,2};
var x{idx} >=0;
minimize obj: min( 3*x[1] + 4*x[2], 2*x[1]+7*x[2]);
con1: x[1]+2*x[2] >= 12;
con2: 2*x[1]+x[2] >= 15;

option solver gurobi;
solve;
display x;
Gurobi 12.0.1: optimal solution; objective 24
4 simplex iterations
3 branching nodes
x [*] :=
1  12
2   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\) で置き換えれば良い．

例題

以下の絶対値が入った線形最適化問題を解け．AMPLを使う場合はabs関数を使えば簡単だが、上の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_1 + x_2 \geq 15 \\ & x_1, x_2 \geq 0 \end{array} \]

reset;
set idx := {1,2};
var x{idx} >=0;
minimize obj: 2*x[1] + 3*x[2] + abs(x[1]-x[2]);
con1: x[1]+2*x[2] >= 12;
con2: 2*x[1]+x[2] >= 15;

option solver gurobi;
solve;
display x;
Gurobi 12.0.1: optimal solution; objective 24
3 simplex iterations
1 branching node
x [*] :=
1  6
2  3
;

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_1 + x_2 \geq 15 \\ & x_1, x_2 \geq 0 \end{array} \]

reset;
set idx := {1,2};
var x{idx} >=0;
maximize obj: -2*x[1] - 3*x[2] + abs(x[1]-x[2]);
con1: x[1]+2*x[2] >= 12;
con2: 2*x[1]+x[2] >= 15;

option solver gurobi;
solve;
display x;
Gurobi 12.0.1: optimal solution; objective -12
6 simplex iterations
3 branching nodes
x [*] :=
1  12
2   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} \]

reset;
set idx := {1,2};
var x{idx} >=0;
maximize obj: x[1] + x[2];
con: 2*x[1]+x[2] <= 30 or x[1]+2*x[2] <= 40;

option solver gurobi;
solve;
display x;
Gurobi 12.0.1: optimal solution; objective 40
5 simplex iterations
3 branching nodes
x [*] :=
1  40
2   0
;

if A then B 条件

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

事象が真か嘘かは0-1変数で表現できる。0-1変数 \(x,y\) に対して、「\(x\)が真ならば\(y\)も真である」という論理制約は、\(x \leq y\) で表現できる。 いずれかが真であるという論理制約は \(x+y=1\) となる。 もう少し難しいものとして、「「この条件が成立しているときには、この制約が必要だ」（if … then … ）という論理制約がある。 「制約Aが成立している場合には，制約Bも成立しなければならない」という条件は， 「Aが成立しない or Bが成立する」と同値である。 ０−１変数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\]

AMPLだと ==> とすると自動的に変換してくれる。

例題

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

\[ \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} \]

reset;
set idx := {1,2};
var x{idx} integer >=0;
minimize obj: x[1] + x[2];
con1: x[2]>=2;
con2: x[1] <= 2 ==> x[2] >=6;

option solver gurobi;
solve;
display x;
Gurobi 12.0.1: optimal solution; objective 5
0 simplex iterations
x [*] :=
1  3
2  2
;

例題 嘘つき島パズル

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

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

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

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

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

# x=1: 旦那, y=1: 奥さん が正直族
# 1. x=1 <=> x+y=0
reset;
var x binary;
var y binary;
constraint: x=1 <==> x+y=0;

option solver gurobi;
solve;
display x,y;
Gurobi 12.0.1: optimal solution
0 simplex iterations
Objective = find a feasible point.
x = 0
y = 1

# x=1: 旦那, y=1: 奥さん が正直族
# 2. x=1 <=> x+y<=1
reset;
var x binary;
var y binary;
constraint: x=1 <==> x+y<=1;

option solver gurobi;
solve;
display x,y;
Gurobi 12.0.1: optimal solution
0 simplex iterations
Objective = find a feasible point.
x = 1
y = 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)，サスロ，ドズル，ガルマ

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

非線形最適化問題

線形でない一般的な関数を扱う最適化問題を非線形最適化問題と呼ぶ。

Rosenbrock関数

まず，無制約の問題の例として，代表的なベンチマーク問題であるRosenbrock関数(Rosenbrock function) を最小化してみよう．

Rosenbrock関数 \[ f(x_1, x_2, \ldots, x_n) = \sum_{i=1}^{n-1} (100(x_i^2 - x_{i+1})^2 + (1-x_i)^2) \] は, \[ -2.048 \leq x_i \leq 2.048 \ \ \ i=1,2,\ldots,n \] で定義され，最適値は \(0\)，最適解は \(x_i=1 (i=1,2,\ldots,n)\) であることが知られている．

reset;
param n = 5;
var x{1..n};

minimize obj: sum{i in 1..n-1} (100*(x[i]^2-x[i+1])^2 + (1-x[i])^2);

option solver ipopt;  #or knitro (for small instances pnly) 
solve;
display x;
Ipopt 3.12.13: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.13, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        9

Total number of variables............................:        5
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  4.0000000e+00 0.00e+00 2.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  3.9091833e+00 0.00e+00 1.40e+01  -1.0 1.00e+00    -  1.00e+00 2.50e-01f  3
   2  3.1769017e+00 0.00e+00 2.74e+00  -1.0 1.16e-01    -  1.00e+00 1.00e+00f  1
   3  2.8732384e+00 0.00e+00 9.24e+00  -1.0 8.59e-01    -  1.00e+00 2.50e-01f  3
   4  2.1171829e+00 0.00e+00 4.21e+00  -1.0 1.59e-01    -  1.00e+00 1.00e+00f  1
   5  1.6905818e+00 0.00e+00 4.71e+00  -1.0 2.82e-01    -  1.00e+00 5.00e-01f  2
   6  1.1447415e+00 0.00e+00 5.84e+00  -1.0 1.90e-01    -  1.00e+00 1.00e+00f  1
   7  7.0591600e-01 0.00e+00 3.59e+00  -1.0 1.34e-01    -  1.00e+00 1.00e+00f  1
   8  4.5380777e-01 0.00e+00 5.23e+00  -1.0 1.62e-01    -  1.00e+00 1.00e+00f  1
   9  2.4021869e-01 0.00e+00 1.68e+00  -1.0 1.19e-01    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.7145893e-01 0.00e+00 5.90e+00  -1.0 2.10e-01    -  1.00e+00 1.00e+00f  1
  11  5.0746431e-02 0.00e+00 3.60e-01  -1.0 9.36e-02    -  1.00e+00 1.00e+00f  1
  12  2.3320006e-02 0.00e+00 1.54e+00  -1.7 2.37e-01    -  1.00e+00 5.00e-01f  2
  13  6.1134648e-03 0.00e+00 9.87e-01  -1.7 1.11e-01    -  1.00e+00 1.00e+00f  1
  14  8.5477417e-04 0.00e+00 4.42e-01  -1.7 7.52e-02    -  1.00e+00 1.00e+00f  1
  15  3.3161505e-05 0.00e+00 9.82e-02  -1.7 3.59e-02    -  1.00e+00 1.00e+00f  1
  16  7.8871311e-08 0.00e+00 5.03e-03  -2.5 8.20e-03    -  1.00e+00 1.00e+00f  1
  17  5.0688121e-13 0.00e+00 1.27e-05  -3.8 4.14e-04    -  1.00e+00 1.00e+00f  1
  18  2.1081641e-23 0.00e+00 8.25e-11  -8.6 1.05e-06    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 18

                                   (scaled)                 (unscaled)
Objective...............:   2.1081641190848290e-23    2.1081641190848290e-23
Dual infeasibility......:   8.2542639389890805e-11    8.2542639389890805e-11
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   8.2542639389890805e-11    8.2542639389890805e-11


Number of objective function evaluations             = 41
Number of objective gradient evaluations             = 19
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 18
Total CPU secs in IPOPT (w/o function evaluations)   =      0.013
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
 
Ipopt 3.12.13: Optimal Solution Found

suffix ipopt_zU_out OUT;
suffix ipopt_zL_out OUT;
x [*] :=
1  1
2  1
3  1
4  1
5  1
;

Weber問題

制約なしの非線形最適化のもう1つの例として、次の問題を考えよう．

砂漠に7 組の家がある．以下の表に各家の位置と 人数が書かれている． この地域には井戸がなく， 水は何キロも離れたところへ，毎日取りにいかねばならないのが 重労働であった． そこでみんなで出資して新しく井戸を掘ることになった． なるべく公平に井戸を位置を決めたいので， 「各家が水を運ぶ距離 × 各家の水の消費量」の総和が最小になる場所に 井戸を掘ることにした．どこを掘っても水が出るものとしたとき，どのようにして掘る場所を決めれば良いだろうか． ただし, 各家の水の消費量は人数に比例するものとする．

家 x座標 y座標 人数
1 24 54 2
2 60 63 1
3 1 84 2
4 23 100 3
5 84 48 4
6 15 64 5
7 52 74 4

この問題は Weber問題 (Weber problem) とよばれ，一般的には以下のように書ける．

Weber問題

家の集合を \(H\) とする． 各家の位置を \((x_i, y_i)\) \((i\in H)\) とする． 井戸の位置を \((X, Y)\) とすれば，家 \(i\) から井戸までの距離は\[ \sqrt{(x_i - X)^2 + (y_i - Y)^2} \] である．家 \(i\) が \(1\) 日に必要とする水の量を \(w_i\) としたとき， \[ \sum_{i\in H} w_i \sqrt{(x_i-X)^2 + (y_i-Y)^2} \] を最小にする \((X, Y)\) を求めよ．

Weber問題は次のように定式化することができる． \[ \begin{array}{ll} \text{minimize} & \sum_{i\in H} w_i z_i \\ \text{subject to} & \sqrt{(x_i-X)^2 + (y_i-Y)^2} \leq z_i \ \ \forall i\in H \\ \end{array} \] この問題の変数は \((X, Y)\) および \(z_i\ (i\in H)\) である．

reset;

model;
param n = 7;
param w{1..n};
param x{1..n};
param y{1..n};

var X;
var Y;

minimize obj: sum{i in 1..n} w[i]* sqrt( (x[i]-X)^2 +(y[i]-Y)^2 );

data;
param:   x  y  w :=
  1   24   54   2  
  2   60   63   1  
  3   1   84   2  
  4   23   100   3  
  5   84   48   4  
  6   15   64   5  
  7   52   74   4  ;

option solver ipopt;
solve;

display X,Y;
Ipopt 3.12.13: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.13, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        3

Total number of variables............................:        2
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.7584643e+03 0.00e+00 1.76e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  7.6545977e+02 0.00e+00 1.13e+01  -1.0 7.22e+02    -  1.00e+00 1.25e-01f  4
   2  6.8813995e+02 0.00e+00 7.99e+00  -1.0 3.39e+01    -  1.00e+00 1.00e+00f  1
   3  6.2797629e+02 0.00e+00 2.13e+00  -1.0 1.57e+01    -  1.00e+00 1.00e+00f  1
   4  6.2366446e+02 0.00e+00 1.13e-01  -1.0 3.47e+00    -  1.00e+00 1.00e+00f  1
   5  6.2364929e+02 0.00e+00 5.03e-05  -2.5 2.16e-01    -  1.00e+00 1.00e+00f  1
   6  6.2364929e+02 0.00e+00 2.71e-10  -5.7 2.16e-04    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 6

                                   (scaled)                 (unscaled)
Objective...............:   6.2364929164124828e+02    6.2364929164124828e+02
Dual infeasibility......:   2.7066260344099646e-10    2.7066260344099646e-10
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.7066260344099646e-10    2.7066260344099646e-10


Number of objective function evaluations             = 14
Number of objective gradient evaluations             = 7
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 6
Total CPU secs in IPOPT (w/o function evaluations)   =      0.006
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
 
Ipopt 3.12.13: Optimal Solution Found

suffix ipopt_zU_out OUT;
suffix ipopt_zL_out OUT;
X = 32.7843
Y = 68.8547

ポートフォリオ最適化問題

制約のある非線形最適化の例として、次の問題を考えよう．

\(100\) 万円のお金を 5 つの株に分散投資したいと考えている． 株の価格は，現在はすべて \(1\) 株あたり \(1\) 円だが， 証券アナリストの報告によると，それらの株の \(1\) 年後の価格と分散は、それぞれ以下の表のように確率的に変動すると予測されている． 目的は\(1\)年後の資産価値を最大化することである． しかしながら，よく知られているように，\(1\) つの株式に集中投資するのは 危険であり，大損をすることがある．「うまく」分散投資するにはどうすれば良いだろうか．

株式 1 2 3 4 5
期待値 (\(r_i\)) 1.01 1.05 1.08 1.10 1.20
標準偏差 (\(\sigma_i\)) 0.07 0.09 0.1 0.2 0.3

この問題では，「うまく」分散投資をする，とあいまいに書かれているので， このままでは最適化問題として定式化することはできない． そこで以下では，この問題に対してMarkowitzモデルと呼ばれる古典的なモデルを考える。

Markowitzモデル

資産 \(M\) 円を持つとき， \(1\) 年後の期待資産価値が \(\alpha M\) 円以上（ただし \(\alpha >1\)）という制約のもとで, 「リスク」を最小化することを考える． ただし，株を $i=1,2,,n $ で表し， 株 \(i\) の \(1\) 年後の価値は期待値 \(r_i\), 分散 \(\sigma^2_i\) の確率分布に したがうと仮定する．

注:上の問題では，それぞれの株の1年後の価格が独立, すなわち，株 \(i\) と株 \(j\) とはお互いに関連がない，という仮定を置いている． A 社と B 社が同じ分野であったりすると，一般には A 社の株の価格と B 社の株の価格には（どちらかが上がるともう一方が下がるというような） 関連がある． このような場合には，各株式 \(i\) の分散を考えるだけでは情報が足りず， 株式 \(i\) と \(j\) の間の共分散 \(\sigma_{ij}\) を考える必要があるが、 基本的な考え方は同じである。

株を \(1, 2, \ldots, n\) とし， 各株 \(i\) の \(1\) 年後の価格を，確率変数 \(R_i\) で表す． ただし，\(R_i\) の期待値は \(r_i\), 分散は \(\sigma_i^2\) とする． 期待値をとる操作を \(\mathbf{E}[\cdot]\) で書けば，次の関係が成り立っている． \[ \mathbf{E}[R_i] = r_i, \ \mathbf{E}[(R_i- r_i)^2] = \sigma_i^2 \ \ \forall i=1,2,\ldots,n \] 株 \(i\) に投資する割合を \(x_i\) とする．すると \(x_1, x_2, \ldots, x_n\) は 非負で和が \(1\) になっているはずである． \[ \sum_{i=1}^n x_i =1, \quad x_i\geq 0, \ \forall i=1,2,\ldots, n \] この投資比率で投資したとき， \(1\) 年後の財産価格は確率変数 \(R_i\) を用いて $ M _{i=1}^n R_i x_i $ と書けるので，期待値が $ M$以上という制約は次のようになる．

\[ M \sum_{i=1}^n r_i x_i \geq \alpha M \]

\(M\) は両辺に現れるので消去でき，結局以下のようになる．

\[ \sum_{i=1}^n r_i x_i \geq \alpha \]

最後に「リスク」とは何かを考えなければならない． Markowitz は「リスク」とは期待値からのずれと解釈し， 確率でいう「分散」をリスクと定義した． 期待値からのずれの \(2\) 乗の期待値は以下のように計算できる．

\[ \begin{aligned} \mathbf{E} \left[\left(M\sum_{i=1}^n R_i x_i - M\sum_{i=1}^n r_i x_i\right)^2 \right] &= M^2 \mathbf{E}\left[\left(\sum_{i=1}^n (R_i - r_i) x_i\right)^2\right] \\ &= M^2 \sum_{i=1}^n \mathbf{E}[(R_i - r_i)^2] x_i^2 \\ &= M^2 \sum_{i=1}^n \sigma_i^2 x_i^2 \end{aligned} \]

Markowitz のモデルではこれが最小化すべき関数であるが， やはり \(M\) は定数なので省略して良く， 結局解くべき最適化問題は以下となる． \[ \begin{array}{lrcl} \text{minimize} &\displaystyle\sum_{i=1}^n \sigma_i^2 x_i^2 \\ \text{subject to} & \displaystyle\sum_{i=1}^n r_i x_i &\geq& \alpha \\ & \displaystyle\sum_{i=1}^n x_i &=& 1 \\ & x_i &\geq& 0 \ \ \ \forall i=1, 2, \ldots, n \end{array} \]

reset;
model;
param n=5;
param sigma{1..n}>=0;
param r{1..n}>=0;
param alpha;
var x{1..n}>=0;

minimize total_risk: sum{i in 1..n} sigma[i]^2*x[i]^2;

subject to expected_return: sum{i in 1..n} r[i]*x[i] >= alpha;
subject to sum_to_one: sum{i in 1..n} x[i] = 1;

data;
param alpha = 1.05;

param : r  sigma  := 
 1   1.01   0.07
 2   1.05   0.09
 3   1.08   0.1
 4   1.10   0.2
 5   1.20   0.3;

option solver ipopt;
solve;
display x;
Ipopt 3.12.13: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.13, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:        5
Number of nonzeros in inequality constraint Jacobian.:        5
Number of nonzeros in Lagrangian Hessian.............:        5

Total number of variables............................:        5
                     variables with only lower bounds:        5
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        1
Total number of inequality constraints...............:        1
        inequality constraints with only lower bounds:        1
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.5299969e-05 9.96e-01 1.09e-01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  6.3818176e-03 1.11e-16 9.74e+00  -1.0 1.97e-01    -  9.22e-02 1.00e+00f  1
   2  1.1843933e-02 0.00e+00 5.10e+00  -1.0 2.44e-01    -  1.00e+00 5.00e-01f  2
   3  9.7480091e-03 0.00e+00 1.00e-06  -1.0 4.12e-02    -  1.00e+00 1.00e+00f  1
   4  9.0505017e-03 0.00e+00 7.95e-03  -2.5 1.53e-02    -  9.95e-01 1.00e+00f  1
   5  4.2737368e-03 0.00e+00 9.12e-03  -3.8 1.23e-01    -  8.86e-01 1.00e+00f  1
   6  2.8077110e-03 0.00e+00 1.50e-09  -3.8 9.49e-02    -  1.00e+00 1.00e+00f  1
   7  2.3378605e-03 0.00e+00 2.75e-03  -5.7 6.77e-02    -  8.12e-01 1.00e+00f  1
   8  2.2249894e-03 1.11e-16 1.84e-11  -5.7 3.27e-02    -  1.00e+00 1.00e+00f  1
   9  2.1956434e-03 0.00e+00 1.84e-11  -5.7 1.18e-02    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  2.1897797e-03 2.22e-16 4.57e-06  -8.6 3.14e-03    -  9.93e-01 1.00e+00f  1
  11  2.1894981e-03 1.11e-16 2.51e-14  -8.6 1.66e-04    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 11

                                   (scaled)                 (unscaled)
Objective...............:   2.1894981151316462e-03    2.1894981151316462e-03
Dual infeasibility......:   2.5059815333960955e-14    2.5059815333960955e-14
Constraint violation....:   1.1102230246251565e-16    1.1102230246251565e-16
Complementarity.........:   3.4371223264188027e-09    3.4371223264188027e-09
Overall NLP error.......:   3.4371223264188027e-09    3.4371223264188027e-09


Number of objective function evaluations             = 14
Number of objective gradient evaluations             = 12
Number of equality constraint evaluations            = 14
Number of inequality constraint evaluations          = 14
Number of equality constraint Jacobian evaluations   = 12
Number of inequality constraint Jacobian evaluations = 12
Number of Lagrangian Hessian evaluations             = 11
Total CPU secs in IPOPT (w/o function evaluations)   =      0.009
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
 
Ipopt 3.12.13: Optimal Solution Found

suffix ipopt_zU_out OUT;
suffix ipopt_zL_out OUT;
x [*] :=
1  0.391754
2  0.270308
3  0.239192
4  0.063172
5  0.0355738
;

\(\alpha = 1.05\), つまり 5% 以上 の利益を期待してこれを解くと，株式を 1 から順に \(0.392\), \(0.270\), \(0.239\), \(0.063\), \(0.036\) の割合で購入するのが最適となることが分かる．

問題：シャッターチャンス

鳥が放物線 \(y=x^2 + 10\) を描いて飛んでいる．いま \((10,0)\) の位置にいるカメラマンが，鳥との距離が最小の地点で シャッターを押そうとしている．鳥がどの座標に来たときにシャッターを押せば良いだろうか？

問題：彗星接近

たくさんの彗星が地球に接近している．地球を原点 \((0,0,0)\) としたとき彗星の描く軌道は曲面 \(2x^2 + y^2 + z^2 = 1000\) 上にあることが 予測されている．地球に最も接近するときの距離を求めよ．

問題：体積最大（最小）の直方体

断面の面積の和が \(100 \mbox{cm}^2\) の直方体で，体積最大のものは何か？また体積最小のものは何か？

問題: 交通割り当て

\(A\) 町から \(B\) 町までは4 本の道が通っている．これらの道は交わることなく2 つの町を繋いでいるのだが， それぞれ移動時間と交通容量が異なっている． 各道の移動時間は， 基本移動時間（台数 \(0\) のときの移動時間），交通容量と通過する車の台数 \(x\) に対する非線形関数になっており，以下の式で得られるものとする． \[ \text{移動時間} = \text{基本移動時間} \times \left(1+ \left(\frac{\text{台数}}{\text{交通容量}}\right)^4 \right) \] 調査の結果，以下の表のようなデータが得られたとき， \(5000\) 台の車が移動する総時間を最小化するには，どのように車を道に割り振ったら良いだろうか？ またそのとき，道ごとの移動時間はどのようになっているだろうか？ また，通過する車の台数が増えたときにはどのようになるだろうか？

道 基本移動時間 交通容量
1 15 1000
2 20 2000
3 30 3000
4 35 4000

問題: 学校配置

あなたの町に新しく2つの学校を作ろうと考えている．現在の学区は5つに分かれており， その中心位置と学生数は以下の表のようになっている． 学生たちの歩く距離の合計を最小にするように新しい学校の位置と学生の学校への割当を決めたい． 学区の中心から学校までの距離はEuclid距離で測定し，それに学生数を乗じた値が歩く距離だとすると， 平面上のどこに学校を作れば良いだろうか？また各学区内の学生をどの学校に割り当てれば良いだろうか？

学区 x座標 y座標 人数
1 0 0 40
2 0 100 40
3 100 0 40
4 100 100 40
5 50 50 40

問題: 経済発注量

以下の仮定に基づく経済発注量問題(economic lot sizing problem)を考える．

  1. 品目（商品，製品）は一定のスピードで消費されており， その使用量（これを需要量とよぶ）は \(1\) 日あたり \(d (>0))\) 単位である．
  2. 品目の品切れは許さない．
  3. 品目は発注を行うと同時に調達される． 言い換えれば発注リード時間（注文してから品目が到着するまでの時間）は \(0\) である．
  4. 発注の際には，発注量によらない固定的な費用（これを発注費用とよぶ） \(F (>0)\) 円が課せられる．
  5. 在庫保管費用は保管されている在庫量に比例してかかり， 品目 \(1\) 個あたりの保管費用は \(1\) 日で \(h (>0)\) 円とする．
  6. 考慮する期間は無限期間とする．
  7. 初期在庫は \(0\) とする．

最適方策は周期的に発注を繰り返すことを示すことができる． 発注間隔を \(T\) としたとき，1 日あたりの平均費用は以下のようになる． \[ \frac{F}{T} + \frac{hdT}{2} \]

発注固定費用 \(F=300\)，在庫費用 \(h=1\)，需要量 \(d=100\) としたとき，費用を最小にする発注間隔を求めよ．