nigoblog

技術系会社のCEOブログ~私的編~

JSで挑戦!線形計画法をプログラミングで解くためのアルゴリズム実装

アルゴリズム JavaScript Github

今回は線形計画法をプログラミングで解いてみました。

その時のプログラムとアルゴリズムを紹介します。

  1. 線形計画法とは?
  2. 線形計画法のアルゴリズム
  3. 線形計画法数式モデル
  4. ピボット列、ピボット行の選択
  5. 掃き出し法の実装
  6. 最終結果
  7. まとめ

このような流れで書いていきます。

線形計画法とは?

簡単に説明すると、
「ある量を最大化させるための変数を、制限の中で決定すること。」
これが一番一言でしっくりくる気がします。
もう少し具体的に例を出すと、
あるwebアプリを2つ作る際、クオリティによって報酬が変わってくるとする。
基本単価が1つを400、もう1つを300とする。
その時、クオリティを上げるための工数がx1, x2とし、報酬をsalesとすると

400*x1 + 300*x2 = sales

となります。ここでwebサイトはそれぞれruby、PHPで書かれているとする。
作成する人は3人、全員ruby、PHPはできるが、そのレベルは違い、単位工数当たりでクオリティを上げる能力がバラバラである。また、それぞれのリソース(全工数)もバラバラであり、次の表のようになります。

人 単位クオリティ毎にかかる工数Ruby PHP リソース
A 60 40 3800
B 20 30 2100
C 20 10 1200

式で表すと

60*x1 + 40*x2 <= 3800
20*x1 + 30*x2 <= 2100
20*x1 + 10*x2 <= 1200

このような状態の時に
Aさん、Bさん、Cさんの工数をどう分配していけば報酬を最大化できるのか。
これを計算によって解くのが線形計画法です。

実際のことを考えると、Aさんは工数はたくさん使えるが、非効率。Cさんは効率的だけど使える時間が少ない。

以上簡単な説明でしたが、線形計画法でした。
この例題を元に実際にそれぞれのwebサイトにどのように工数を分配していけばよいかをチェックするため、線形計画法を解いて行きましょう。

ちなみに数字とアルゴリズムですが、こちらのサイトを参考にしました。わかり易かったです。
見てる途中全然気づきませんでしたが、自分の卒業した学校、学科の資料でした。笑
どこかで見たこと在ると思ったら…
システム工学

線形計画法のアルゴリズム

アルゴリズムを次に示します。言葉は後に説明します。

  1. ピボット列の選択
  2. ピボット行の選択
  3. 掃き出し法によりピボット行、ピボット列の変数を求める
  4. 売上にかかる係数が全て0より大きいか調べる。
    1. 0より大きい数があれば1に
    2. 全て0以下であれば終了

以上アルゴリズムはこんな感じです。
ただこれだけだと伝わらないことも多いので、次に線形計画法の数式モデルを示します。

線形計画法数式モデル

それでは最初の例題を数式に変換します。
すると次のようになります。

60x1 + 40x2 <= 3800
20x1 + 30x2 <= 2100
20x1 + 10x2 <= 1200
400x1 + 300x2 = sales

ここで不等式を等式にするためにスラック変数というものを導入します。
すると次のようになります。

60x1 + 40x2 + y1 = 3800
20x1 + 30x2 + y2 = 2100
20x1 + 10x2 + y3 = 1200
400x1 + 300x2 = 0

新たに追加されたy1, 2, 3がスラック変数です。
これの物理的な意味はAさんBさんCさんの工数余りです。
つまりこれ以上クオリティを上げても報酬がこれ以上あがらないという時にそれぞれが余力を残すための数値です。
よくwebサイトとかアプリとか見てると無駄な機能とかありますよね。
線形計画法で余力まで計算してあげればそういうこともなくなってくるのかなと思います。

これらをプログラミングで表すために配列で表します。
次のように定義します。

var x    = new Array();
var y    = new Array();
var profit = new Array();
var stock  = new Array();

配列xは各メンバーの係数、yはスラック変数、profitはx1, x2の単価、stockはメンバーのリソースと報酬
こちらを数値に従い、初期化すると、

function xArrayInit(){
  x[0] = [60, 40];
  x[1] = [20, 30];
  x[2] = [20, 10];
}
function yArrayInit(){
  for(var i=0; i<4;i++){
    for(var j=0; j<3;j++){
      y[i] = [0, 0, 0];
    }
  }
  for(var i=0; i<4;i++){
    for(var j=0; j<3;j++){
      if(i == j){
        y[i][j] = 1;
      }
    }
  }
}
function profitArrayInit(){
  profit[0] = 400;
  profit[1] = 300;
}
function stockArrayInit(){
  stock[0] = 3800;
  stock[1] = 2100;
  stock[2] = 1200;
  stock[3] = 0;
}

以上のようになります。
yは

60x1 + 40x2 + y1 + 0y2 + 0y3 = 3800
20x1 + 30x2 + 0y1 + y2 + 0y3 = 2100
20x1 + 10x2 + 0y1 + 0y2 + y3 = 1200
400x1 + 300x2 + 0y1 + 0y2 + 0y3= 0

のようなイメージだと思って下さい。
ではモデル化をしましたので、いよいよ計算に入っていきます。
(※売上は「最大化」させると説明しましたが、計算上では「最小化」になります。)

ピボット列、ピボット行の選択

ピボット列は
「最初に上げていくべき数値」
という意味で捉えて下さい。
どういうことかといいますと、

400x1 + 300x2 = 0

となっております。単価がそれぞれ400, 300なのでまず400から上げていくべきだと単純に考えられます。なのでとりあえず400にかかる変数をどれだけ上げられるかをチェックします。
そしてこの列こそがピボット列といいます。
次にピボット行の選択です。
まずはピボット列に注目。
メンバーごとにリソースをピボット列で割っていきます。
するとAさんは

3800 / 60 = 63.3333...

Bさん

2100 / 20 = 105

Cさん

1200 / 20 = 60

というようにCさんが先に力尽きてしまいます。なのでCさんに合わせるためここからCさん基準で計算をすすめていきます。つまりピボット行はCさんの行ということになります。

まとめると
ピボット列は
「単価と売上の式より、最大の係数のものから選択」
ピボット行は
「最大限使えるリソースに対するピボット列にかかる係数の商で最小のもの」
となります。

実装には次の関数を作成しました

function pivotCol(line, profitArray, yArray){
  var mProfit = profitArray[0];
  var pivo = 0;
  var count = profitArray.length;
  for(var i=1; i < count;i++){
    if(mProfit < profitArray[i]){
      mProfit = profitArray[i];
      pivo = i;
    }
  }
  var yCount = yArray[0].length;
  for(var j=0; j < yCount;j++){
    if(mProfit < yArray[line][j]){
      mProfit = yArray[line][i];
      pivo = count + j;
    }
  }
  return pivo;
}
function max(array){
  var max = 0;
  var count = array.length;
  for(var i=0; i<count;i++){
    if(max < array[i]){
      max = array[i];
    }
  }
  return max;
}
function pivotLine(pivoCol, xArray, yArray, stockArray){
  var pivoLine = 0;
  if(pivoCol < xArray[0].length){
    var empty = max(stockArray);
    var count = xArray.length;
    for(var i=1; i<count; i++){
      if(empty > stockArray[i]/xArray[i][pivoCol] && 0 < stockArray[i]/xArray[i][pivoCol]){
        empty = stockArray[i]/xArray[i][pivoCol];
        pivoLine = i;
      }
    }
  }else{
    pivoCol = pivoCol - xArray[0].length;
    var empty = max(stockArray);
    var yCount = yArray[0].length;
    for(var j=0; j<yCount; j++){
      if(empty > stockArray[j]/yArray[j][pivoCol] && 0 < stockArray[j]/yArray[j][pivoCol]){
        empty = stockArray[j]/yArray[j][pivoCol];
        pivoLine = j;
      }
    }
  }
  return pivoLine;
}

pivotCol()がピボット列を求める関数。
pivotLine()がピボット行を求める関数。
max()はpivotLine()の中で仕様しております。
pivotCol()、pivotLine()共に、スラック変数がピボットになった際にも対応出来ます。

掃き出し法の実装

掃き出し法は細かく説明しませんが、
ピボット列、ピボット行を基準に計算し、
ピボット列の中でピボット行以外の係数を0にしていく計算方法です。
実装結果は次のようになります。

function normalize(pivoLine, pivoCol, xArray, yArray, stockArray, profitArray)
{
  if(pivoCol < xArray[0].length){
    var pivoCoef = xArray[pivoLine][pivoCol];
  }else{
    pivoCol = pivoCol - xArray[0].length;
    var pivoCoef = yArray[pivoLine][pivoCol];
  }
  var xCount = xArray[0].length;
  for(var i=0; i<xCount;i++){
    xArray[pivoLine][i] = xArray[pivoLine][i]/pivoCoef;
  }
  var yCount = yArray[0].length;
  for(var i=0; i<yCount;i++){
    yArray[pivoLine][i] = yArray[pivoLine][i]/pivoCoef;
  }
  stock[pivoLine] = stock[pivoLine]/pivoCoef;
}
function swapMatrix(pivoLine, xArray, yArray, stockArray){
  var xCount = xArray[0].length;
  for(var i=0; i<xCount;i++){
    swapDoubleArray(0, i, pivoLine, i, xArray);
  }
  var yCount = yArray[0].length;
  for(var i=0; i<yCount;i++){
    swapDoubleArray(0, i, pivoLine, i, yArray);  
  }
  swapArray(0, pivoLine, stockArray);
}
function swapDoubleArray(line, col, line2, col2, targetArray){
  var temp;
  temp = targetArray[line][col];
  targetArray[line][col] = targetArray[line2][col2];
  targetArray[line2][col2] = temp;
}
function swapArray(line, line2, targetArray){
  var temp;
  temp = targetArray[line];
  targetArray[line] = targetArray[line2];
  targetArray[line2] = temp;
}
function transMatrix(line, pivoLine, pivoCol, xArray, yArray, stockArray){
  if(pivoCol < 3){
    var pivo = xArray[line][pivoCol];
  }else{
    pivoCol = pivoCol - xArray[0].length;
    var pivo = yArray[line][pivoCol];
  }
  var xCount = xArray[0].length;
  for(var i=0;i<xCount;i++){
    xArray[line][i] = xArray[line][i] - xArray[pivoLine][i]*pivo;
  }
  var yCount = yArray[0].length;
  for(var i=0;i<yCount;i++){
    yArray[line][i] = yArray[line][i] - yArray[pivoLine][i]*pivo;
  }
  stockArray[line] = stockArray[line] - stockArray[pivoLine]*pivo;
}
function transProfitMatrix(line, pivoLine, pivoCol, xArray, yArray, stockArray, profitArray){
  if(pivoCol < 3){
    var pivo = profitArray[pivoCol];
  }else{
    pivoCol = pivoCol - xArray[0].length;
    var pivo = yArray[line][pivoCol];
  }
  var pCount = profitArray.length;
  for(var i=0; i<pCount;i++){
    profitArray[i] = profitArray[i] - xArray[pivoLine][i]*pivo;
  }
  var yCount = yArray[0].length;
  for(var i=0; i<yCount;i++){
    yArray[line][i] = yArray[line][i] - yArray[pivoLine][i]*pivo;
  }
  stockArray[line] = stockArray[line] - stockArray[pivoLine]*pivo;
}

まず一つ一つ説明します。
normalize()はピボット行のピボット列の係数を1にする関数。
行列なのでそれに伴い、他の列の係数も変化します。
つまりこの関数の後

(A) 60x1 + 40x2 + 1y1 + 0y2 + 0y3 = 3800
(B) 20x1 + 30x2 + 0y1 + 1y2 + 0y3 = 2100
(C) 1x1 + 0.5x2 + 0y1 + 0y2 + 0.05y3 = 60
400x1 + 300x2 + 0y1 + 0y2 + 0y3 = 0

となります。
次にswapMatrix()ですが、ピボット行を一番上に持ってくる関数です。

(C) 1x1 + 0.5x2 + 0y1 + 0y2 + 0.05y3 = 60
(A) 60x1 + 40x2 + 1y1 + 0y2 + 0y3 = 3800
(B) 20x1 + 30x2 + 0y1 + 1y2 + 0y3 = 2100
400x1 + 300x2 + 0y1 + 0y2 + 0y3 = 0

となります。
最後にtransMatrix()とtransProfitMatrix()ですが、ピボット行以外のピボット列の係数を0にする関数です。

(C) 1x1 + 0.5x2 + 0y1 + 0y2 + 0.05y3 = 60
(A) 0x1 + 20x2 + 0y1 + 1y2 - 1y3 = 900
(B) 0x1 + 10x2 + 1y1 + 0y2 - 3y3 = 200
0x1 + 100x2 + 0y1 + 0y2 - 20y3 = -24000

以上掃き出し法の部分を終了条件を満たすまで繰り返します。

最終結果

最終結果の前に終了条件の実装を示します。
それは次のようになります。

function finishCheck(profitArray, yArray, line){
  var check = false;
  var count = profitArray.length;
  for(var i=0;i<count;i++){
    if(profitArray[i] > 0){
      check = true;
      return check;
    }
  }
  var yCount = yArray[0].length;
  for(var j=0;j<yCount;j++){
    if(yArray[line][j] > 0){
      check = true;
      return check;
    }
  }
  return check;
}

つまり、元々次の式

400x1 + 300x2 + 0y1 + 0y2 + 0y3 = 0

だったものが掃き出し法の変形により
係数全てが0以下になるとこれ以上変化させると最小値(最大値)ではなくなってしまうということになります。
この条件になるまで計算した結果が

0x1 + 0x2 - 0.4y1 + 0.2y2 + 1y3 = 100
0x1 + 1x2 - 0.02y1 + 0.06y2 + 0y3 = 50
1x1 + 0x2 + 0.03y1 - 0.04y2 + 0y3 = 30
0x1 + 0x2 - 6y1 - 2y2 + 0y3 = -27000

となりました。この結果より

sales = 27000
x1 = 30
x2 = 50
y3 = 100

つまり

項目 最終結果
売上 27000
rubyのアプリに費やす全リソース 30
phpのアプリに費やす全リソース 50
Cさんの余力 100

という結果になりました。
このように工数を使えば売上を最大化させるということを計算で求めることができました。

まとめ

よく会社では優先度とかそういう言葉が使われますが、
うまく売上単価や工数をモデリングさえ出来れば線形計画法を使い無駄を減らした開発ができるかと考えます。
つまり費用対効果や投資対効果の最適化です。
まぁ現実問題こんな簡単には行かないと思いますが、
論理的な根拠から優先度を算出できるのは説得力を上げられるのでプレゼントかには使えそうですよね!
それ以外にも汎用性の高いものなので是非やっておきましょう!!

また最適化に利用できるということは機械学習などにも応用できる可能性を秘めています。
(機械学習用にはもっとよい最適化手法があるけど)

というわけで線形計画法のプログラミングでした。

関連記事

jQueryプラグインを作成するための3つのポイント - nigoblog
これから先もプログラマーを続けるために重要なこと - nigoblog
プログラマーは今こそアルゴリズムを書くべき!!~モンテカルロ法でπを計算してみよう(コードあり)~ - nigoblog
プログラマーは今こそアルゴリズムを書くべき!!2~再帰アルゴリズムでハノイの塔を解く~ - nigoblog
TDD(テスト駆動開発)でハノイの塔の実装をしてみる~TDD超入門~ - nigoblog
ルービックキューブにハマる~そしてアルゴリズムへ~ - nigoblog
ルービックキューブの実装をPHPでしてみた~ルービックキューブクラス~ - nigoblog

関連書籍

データ分析の基本と業務 (仕組みが見えるゼロからわかる)

データ分析の基本と業務 (仕組みが見えるゼロからわかる)


様々な分析法が書かれており、そのうちの一つに線形計画法も!