一般化線形モデル

一般化線形モデル(generalized linear model; GLM)はPythonの統計解析パッケージstatsmodelsの中核を成す重要なモデルである.しかし,市販の日本語のテキストだと間違いだらけだったり,定義が不足していたり,動機付けなしで結果だけ羅列してあったりしていまいちなので,簡単に説明しておく.

GLMは名前の通り,指数関数族で表現できる様々なモデルを一般的に扱うためのものであるが,背景にあるのは線形モデルだ.

基本となる線形回帰だと,独立変数 x^{(i)} を用いて従属変数 y^{(i)} を推定する.上付き添え字の(i) はトレーニングデータのインデックスを表す.評価関数は最小自乗誤差であり,それを最小にするような重みベクトル w を求める.

通常の線形回帰(最小自乗モデル)は,一般化線形モデル的に見直すと以下のように解釈できる.

  1. 従属変数 y^{(i)} は平均 \mu^{(i)},標準偏差 \sigma の正規分布 N(\mu^{(i)},\sigma^2) にしたがう.
  2. 線形予測子 z^{(i)} を独立変数 x^{(i)} を用いて z^{(i)} = w x^{(i)} と定義する.ここで w は最適化するパラメータ(重み)である.
  3.   リンク関数 g を用いて  \mu^{(i)} と  z^{(i)} を繋ぐが,線形モデルでは g(\mu) =\mu である.

ロジスティック回帰は,タイタニック号で生き残ったか否かを予測するために用いられる.つまり,従属変数が 01 かになるので,一般化線形モデルでの仮定は以下のようになる.

  1. 従属変数 y^{{i}} は平均 \mu^{(i)} (表が出る確率)のコイン投げの分布にしたがう.
  2. 線形予測子 z^{(i)} を独立変数 x^{(i)} を用いて z^{(i)} = w x^{(i)} と定義する.(この部分は全部共通)
  3.   リンク関数 g を用いて  \mu^{(i)} と  z^{(i)} を繋ぐが,\mu は確率なので [0,1] の範囲しかとらない,一方, z は線形予測子なので [-\infty,+\infty] の定義域をもつ.これを繋ぐために以下のリンク関数 g を用いる.

        \[  z = g(\mu) = \log \left( \frac{\mu}{1-\mu} \right) \]

    歴史的な都合で g\mu から z への写像となっているが,逆写像として考えた方がわかりやすい.すなわち,線形予測子 z から分布の平均 \mu を逆写像 g^{-1} で写すのである.この関数は

        \[ \mu = \frac{ \exp (z) }{ 1+\exp (z)} \]

    となり,いわゆるロジスティック関数である.

Poisson回帰は救急車の出動回数などを予測する際に用いられる.このの場合には,従属変数が 0 以上の値になるので,一般化線形モデルでの仮定は以下のようになる.

  1. 従属変数 y^{{i}} は平均 \mu^{(i)} のPoisson分布にしたがう.
  2. 線形予測子 z^{(i)} を独立変数 x^{(i)} を用いて z^{(i)} = w x^{(i)} と定義する.(この部分は全部共通)
  3.   リンク関数 g を用いて  \mu^{(i)} と  z^{(i)} を繋ぐが,\mu0 以上で z[-\infty,+\infty] の定義域をもつ.これを繋ぐために以下のリンク関数 g を用いる.

        \[  z = g(\mu) = \log (\mu) \]

    g の逆写像は指数関数

        \[ \mu = \exp (z) \]

    である.

 

新技術がSCMに与える影響

新しい技術がSCMに与える影響について考察する.

1. モバイルコンピューティング:インターネットに接続した携帯端末の普及は急激に進んでいる.これを利用してSCMを効率化する方法はいくらでもあるのだが,現状では巧く使っている企業は少ない.これは今すぐにでも使える技術であり,その影響は多大である.たとえば,宅配便の配達における再配達が問題になっているが,訪問予定時刻を携帯端末に通信し,不在の場合には訪問希望時間を入力できるような簡易アプリを配布しておくだけで,再配達は大幅に減少するだろう.訪問時間の変更のためにはポイントもしくは別料金が必要で,変更なしで在宅した場合にはポイントが増え,配達の料金に使えるような工夫をすれば,多くのユーザがアプリをインストールするだろう.このようなシステムを稼働させるためには,きちんと配送計画をたてておく必要があり,さらには情報が変化した場合の再ルーティングなどが必要になる.現在の配送計画の技術を使えば,これはそれほど難しいことではない.

2.自動運転車:自動運転車はすでに実用の寸前まできている.後は法律をどうクリアするかだが,高速道路から規制が緩和され,長距離便が自動運転トラックで運用できるようになると考えられる.この技術はおそらく欧米では10年後くらいに実用になるが,日本だけ法律の立ち後れが危惧される.これにより,長距離便の費用が低減されるので,サプライ・チェイン全体の設計をし直す必要が出てくる.

3.ドローン:一時期流行したドローンだが,サプライ・チェインに与える影響は小さい.配達をドローンで行うというと聞こえは良いが,積載可能な荷物の大きさに制限があることや,バッテリーがもたないため長短距離輸送しかできないことなどがネックになり,実用化は遠いと考えられる.超短距離の超高速輸送で,かつ軽量のものなら可能であるが,やはりサプライ・チェイン全体に影響は小さいと言わざるを得ない.アマゾンなどは,充電機能付きの宅配ポストをいたるところに配置するとか,飛行船をとばしてそこからドローンを発進させるなどの特許を出しているが,実用になるとは考えにくい.

4.3Dプリンタ:これはサプライ・チェインというより製造に革命をもたらす「可能性」がある.ただ今すぐに革命が起きるという訳ではなく,しばらくは楽しい玩具として普及し,やがて一部の製造に利用する時代がくると考えられる.

5.クラウド:これは新しい技術という訳ではなく,時代の流れで当然のことなのだが,自社でサーバーをたてて運用する時代は過去のものになりつつある.クラウドは情報の漏洩が心配だというのは逆であり,自社サーバーの方がはるかに危険だ.古い技術に固執した情報担当がいる会社は,新しい技術が安全で安価なことを勉強して移行していかないといけない.

 

Service network design problem

サービスネットワーク設計問題とは,宅配便などの長期レベルのネットワークを設計するためのモデルである.このようなモデルは諸外国では成功しているが,我が国ではデータ未整備や意思決定におけるストラテジックなビジョンが不足しているため,未だ使用されていない.

Consider a service network:

Center (i) => Base (k) => Base (\ell ) =>Center (j)

Number of centers is more than 1000.

Number of bases is around 200.

Given:

* c_{ik} : transportation cost of a vehicle between points i and k
* Q: capacity of vehicles
* D_{ij}: demand from origin i to destination j

Variable:

* x_{ik}: =1 if center i is assigned to base k, =0 otherwise
* y_{k \ell}: number of vehicles between k and \ell

Since the number of variables is too large, we need to restrict the search domain in the following way:

For center i, the candidate bases to be assigned are restricted to the “near” bases.

For center i, let C^*_i is the cost (or distance) to the nearest base.

Given a parameter \alpha > 1, we restrict the set of bases to be assigned to center i:

    \[\{ k \in Base | C_{ik} \leq \alpha C^*_i \}\]

If x_{ik}=1 and x_{j \ell}=1, the flow volume between bases j and \ell is increaded by D_{ij}.

We get:

    \[\sum_{i,j} D_{ij} x_{ik} x_{j \ell} \leq Q y_{k \ell} \ \ \ \forall k,\ell.\]

We also get the number of vehicles between center i and base k:

    \[\sum_{j} D_{ij} x_{ik} \leq Q y_{ik} \ \ \ \forall i,k\]

and

    \[\sum_{j} D_{ji} x_{ik} \leq Q y_{ki} \ \ \ \forall i,k\]

Capacity of base k is less than or equal to CAP_k:

    \[\sum_{i} \sum_{j} (D_{ij} + D_{ji}) x_{ik} \leq CAP_k \ \ \ \forall k\]

The objective function

    \[\sum_{i,j} c_{ij} y_{ij}\]

should be minimized under the constraints:

    \[\sum_{k} x_{ik} =1 \ \ \forall i\]

For the MIP solver, we need to rewrite

    \[\sum_{i,j} D_{ij} x_{ik} x_{j \ell} \leq Q y_{k \ell} \ \ \ \forall k,\ell.\]

into linear constraints:

By intruducing a new 0-1 variable z_{ikj\ell} that is equal to 1 if and only if x_{ik}=x_{j\ell}=1,
we get:

    \[x_{ik} + x_{j \ell} -1 \leq z_{ikj\ell}\]

and

    \[\sum_{i,j} D_{ij} z_{ikj\ell} \leq Q y_{k \ell} \ \ \ \forall k,\ell.\]

バケツリレー方式

人道支援ロジスティクス(災害時物流)における輸送について考える.

直列多段階モデルの仮定は,以下の通り.

地点は被災地に向かい 1,2,\ldots,n と番号付けられており,被災地を表す点を n とする.
i= 1,2,\ldots,n-1 に対して,点i から点 i+1 の間だけ輸送を行う.(以下の節ではこの条件を緩和する.)
点~1 には支援物資が無限にあるものと仮定し,運搬車はピストン輸送で下流(地点 n 方向)に物資を輸送する.
n では一定の速度で需要が発生し,支援された物資を消費する.
輸送は地点間ごとに容量が決められた運搬車で行い,運搬車は満載で輸送を行い,
到着地点ですべての荷を降ろし,出発地点に再び戻る.
運搬車の作業速度は与えられており,これには地点間の移動と荷を積載したり降ろしたりする時間も含まれるものとする.

地点 i から地点 j (>i+1) への輸送を許した場合を考える.実際には,最適な輸送経路は分からないので,地点 i から地点 i+1 への輸送を行った際に,ある条件が成立していたら,地点 i+1 を跳ばして,地点 i+2 に向かうものとする.これがバケツリレー方式である.

バケツリレー方式を評価するために,最適な運搬車数を求めるモデルを構築する.

記号の定義は以下の通り.

[d:] 地点 n(被災地)における単位期間(たとえば日)における需要量
[K:] 運搬車のタイプの集合(たとえば5トン車と10トン車など)

[x_{ij}:] 地点 i から j への輸送を行う運搬車の数

[C_{ij}:] 地点 i から j への輸送を行う運搬車の積載容量上限
[t_{ij}:] 地点 i から j への輸送を行う運搬車の作業時間.これは単位期間(日)に対する割合として表されているものとする.
たとえば,運搬車が1~日に3~回転可能な場合には,t_i1/3 となる.
[\lambda_{ij}:] 地点 i から j への,運搬車1~台あたりの単位期間あたりの輸送量 =C_{ij}/t_{ij}
[U_{k}:] 運搬車タイプ k の台数上限
[A_k:] 運搬車タイプ k で輸送可能な地点対の集合
[M_i:] 地点 i にある荷さばき施設における積み替え可能な運搬車の台数の上限.
これは入庫運搬車数と出庫運搬車数の和の上限であると仮定する.

    \[ \begin{array}{ l l l } \mbox{max.} & d & \\ \mbox{s.t.} & \sum_j \lambda_{ji} x_{ji} - \sum_j \lambda_{ij} x_{ij} = \left\{ \begin{array}{l l } d & i=n \\ 0 & \forall i=2,3,\ldots,n-1 \\ -d & i=1 \end{array} \right. \\ & \sum_{ (i,j) \in A_k} x_{ij} \leq U_k & \forall k \in K \\ & \sum_{j} x_{ji} + \sum_{j} x_{ij} \leq M_i & \forall i=1,2,\ldots,n \\ & x_{ij} \geq 0 & \forall i,j, i<j \\ \end{array} \]

実際問題の定式化(1)

実際問題を定式化するにはちょっと頭を使う必要があります.
ここでは,そのような実際問題を定式化(モデル化)するための,ちょっとしたパズルをご紹介します.

こんな問題を考えます.

大型,中型,小型のトラックを,3種類のバースで荷捌きしています.
バースは大型用,中型用,小型用の3種類があって,大型トラックは大型用のみ,中型トラックは大型用と中型用で,小型トラックはどのバースでも荷捌きができます.

バースの数が大型は5台分,中型は2台分,小型は4台分しかないとき,時間帯別に到着するトラックの台数の上限を表す制約をどのように表現するかを考えます.

トラックをどのバースに割り当てるのかを表す0-1変数を用いるのが普通の考え方ですが,大規模な問題の一部としてそのような表現で定式化するのは,あまり効率が良くありません.

トラックのサイズによるバースへの割り当てが入れ子構造になっていることを着目すると,以下のように定式化すれば良いことが分かります.

大型トラックの台数 ≦ 5
大型トラックの台数+中型トラックの台数 ≦ 5+2=7
大型トラックの台数+中型トラックの台数 +小型トラックの台数 ≦ 5+2+4=11

大型のバース(5台分)には大型車のみなので5以下ですが,もし大型トラックが1台しかこないときには,中型は大型バースに4台入れることができるので,合計5台入れることになります.

このように,入れ子構造を利用して無駄な0-1変数を入れずに定式化することで,格段に大きな問題を解くことができるようになります.

Knapsack Problem

複数制約ナップサック問題は比較的MIPソルバーで解きやすい問題である.比較的古い論文でChu-Beasleyが遺伝的アルゴリズムを, Glover-Kochenbergerが禁断探索法のテストをするために作成したベンチマーク問題例は,Gurobiである程度までは解くことができる.

計算時間のグラフを以下に示す.
問題のサイズと計算時間

これを見ると,変数の数が2500を超えるようになると厳密解を出すのに多大な時間を要することが分かる.しかし,どの問題例も短時間で上界と下界のギャップは比較的小さくなっており,途中で計算を打ち切って得られる近似解は,十分に実用に耐えうると考えられる.

最大の問題例であるmk_gk11_100x2500(この問題例は250000個の変数を含む!)と呼ばれる問題の計算ログの一部を以下に示す.

Found heuristic solution: objective 93813.000000
Root relaxation: objective 9.527749e+04, 1007 iterations, 0.37 seconds
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time

0 0 95277.4913 0 100 93813.0000 95277.4913 1.56% – 0s
H 0 0 95155.000000 95277.4913 0.13% – 1s
H 0 0 95218.000000 95277.4913 0.06% – 5s
0 2 95277.4912 0 100 95218.0000 95277.4912 0.06% – 12s

これによると分枝の前のヒューリスティクスで93813という解を見つけ,分枝開始から12秒後に誤差0.06%の解を得ている.

他のメタヒューリスティクス(詳細については,「メタヒューリスティクスの数理」(共立出版)を参照)との比較はまだ行っていないが,近似解法としてのGurobiに対抗するのは難しいと考えられる.

SCOP Excel Interfaces

従来では,SCOP (Solver for Constraint Programming)は,テキストファイル,Python,C++やVBからの呼び出しなどのインターフェイスを準備してあったが,Excelインターフェイスも準備された.これによって,Excelのソルバーアドインと同じように,制約最適化ソルバーを呼び出すことができるようになる.

以下に示すのは,看護婦スケジューリングのExcelインターフェイスの例である.

看護婦スケジューリングの例

サンプルやお試し版は準備中であるが,ベータ版のテストなどを希望される方は,LOGOPTの問い合わせフォームからご連絡ください.

Experimental Analysis 12

あたらしい数理最適化 -GurobiとPython言語で解く- (近代科学社) 2012年 久保幹雄,ジョア(ジョアン)・ペドロ・ペドロソ,村松正和,アブドル・レイス 著では,実験の結果は示していない.今回は,実験的解析の12回目であり,非線形の設置費用を持つ施設配置問題に対する実験的解析である.

非線形関数は規模の経済や安全在庫費用を近似した平方根の関数であるが,これを区分的線形で近似して扱っている.個々での目的は,様々な区分的線形近似の手法の中で,どれが最も高速化を現実的な問題で調べることである.

以下に,区分数を100とした場合の結果を示す.

非線形施設配置問題に対する実験(区分数100)

図中の記号は,以下の通りである.

mselect 多重選択定式化(multiple selection model)
cc_dis 非集約型凸結合定式化(disaggregated convex combination model)
cc_dis_log 対数個の変数を用いた非集約型凸結合定式化(disaggregated convex combination model with a logarithmic number of variables)
cc_agg 集約型凸結合定式化(aggregated convex combination model)
cc_agg_log  対数個の変数を用いた集約型凸結合定式化(aggregated convex combination model with a logarithmic number of variables)
sos タイプ2の特殊順序集合を用いた定式化 (model using sos constraints of type 2)

結果としては,単純なタイプ2の特殊順序集合を用いたものが最も良いので,これを推奨する.以前の研究では,対数個の変数を用いた定式化が良いという結論を出しているものもあるが,Gurobiのようなモダンなソルバーでは,単純にSOSと宣言する方が,良い結果を出すと考えられる.

次いで,対数個の変数を用いた集約型凸結合定式化,対数個の変数を用いた非集約型凸結合定式化(であるが,これらは複雑な制約を付加する必要がある.3変数以上の(分離不能な)非線形関数の近似の際には,これらの定式化も推奨される(SOSは使えないからだ).ほぼ同じ性能で,集約型凸結合定式化も有効である.この定式化は下界が弱いので使うべきではないと,従来の研究では結論づけられていたが,変数の少なさや簡潔さがそれを補うらしい.

良い下界を算出し,理論的に優れていると言われていた非集約型凸結合定式化は,点数200を超えると解けなくなる.多重選択定式化も同じである.他にも線形最適化の生みの親であるDantzigらが提案した定式化もあるが,やはり論外であり,実務的には使うべきではない.

Experimental Analysis 11

Weber問題と呼ばれる施設配置問題は,二次錘最適化問題として定式化できる.

平面上に一様ランダムに発生させた点に対するWeber問題に対する実験を以下に示す.

Wer問題に対する計算機実験

点の数は10から10万まで増やしているが,計算時間はほとんど線形に増加していることが見て取れる.Gurobiの二次錘最適化は高速であり,ポートフォリオ問題も同様に,ヒューリスティクスなどに頼るべきではないと考えられる.

Experimental Analysis 10

あたらしい数理最適化 -GurobiとPython言語で解く- (近代科学社) 2012年 久保幹雄,ジョア(ジョアン)・ペドロ・ペドロソ,村松正和,アブドル・レイス 著では,実験の結果は示していない.今回は,実験的解析の10回目であり,スケジューリング問題に対する実験的解析である.

実験に用いたのはランダムに生成した1機械リース時刻付きの納期遅れ最小化問題である.すべてのデータは[1,99]の一様乱数にしたがって生成した.以下に示すのは Intel(R) Xeon(R) CPU E5-2687W 0 @ 3.10GHz RAM: 64 GB でのCPU時間である.

スケジューリング問題に対する計算時間の変化

図中で,lopは線形順序付け定式化,tiは時刻添え字定式化,disjは施設定式化を表す.結果から,最も良いのは線形順序付け定式化,次いで時刻添え字定式化,最も悪いのは大きな数(M)を用いた離接定式化である.

ただし,時刻添え字定式化は,時刻を表す数字が大きくない場合には,より高速になることが予想される.

いずれにせよ,1機械の簡単な問題に対してもGurobiでは比較的長い計算時間がかかることから,スケジューリング問題に対しては,制約プログラミングや専用のメタ解法に基づくソルバー(たとえばSCOPやOptSeq)が実務的に推奨されることが分かる.