GLPKで線形計画法

このページ(GLPKで楽しく最適化しよう!)を参考に、修論の研究でしているネットワークのスループット最適化をしてみます。また後半では、GLPKを自分で書いたCのプログラムに組み込んで用いる方法を示します。

問題

以下のようなスイッチ・ノードのネットワークにおいて、リンクのバンド幅l0-l6が与えられているものとします。ここで、二つの転送t0とt1が存在するとします。t0は二つのノードに跨るパイプライン転送で、t1は1ノードから1ノードへの単純な転送です。ここで、各ノードに到着するスループットを最大化するような最適化を考えてみます。t0とt1のバンド幅はl3とl4でリンクを共有しているため和で制約されているほか、経路中の各リンクの転送量には最大値(1Gbpsなど)があります。ここで、それぞれの転送にどれだけの転送速度を割り振ることでスループットが最大化されるかを考えてみます。
l0-l6までのバンド幅は適当に下のように設定してみます。
l0 l1 l2 l3 l4 l5 l6
20Mbps 30Mbps 25Mbps 32Mbps 27Mbps 19Mbps 23Mbps

モデル化

まず数式で表さないといけません。 GLPKを初めとするソルバというソフトでは、モデルとバラメータに分けて考えます。 t0とt1が使うリンクを表すために、リンクと転送の関係の行列を書いてみます。 その転送があるリンクを使うときには1, 使わない時には0を書きます。
- l0 l1 l2 l3 l4 l5 l6
t0 1 1 0 1 1 1 0
t1 0 0 1 1 1 0 1
制約式は これをモデルファイル
# trans.mod
set Transfers; # 独立要素 (商品:後で詳しく定義しないといけない)
set Links;

param bandwidth{Links}; # Prodに依存するパラメータ (定数)
param uselink {Links,Transfers}; #
param weights {Transfers};
var throughput{Transfers} >= 0; # 変数集合 (単独の変数ならvar x;)

maximize THROUGHPUT: sum {t in Transfers} throughput[t] * weights[t];
# for p in Prod: X[p] * value[p]

s.t. BANDWIDTH{l in Links}: sum{t in Transfers} throughput[t] * uselink[l,t] <= bandwidth[l];

# 全てのl(リンク) について、全てのt (転送)の合計がbandwidth(l)を超えない
# - Linksの要素lは値ではないから、l = 100みたいな式は現れないし、l*10みたいな式も無意味
# - bandwidth[l]やuselink[l,t]は具体的な値だが、予め与えられる (定数)
データです。
# trans.dat
set Links := l0 l1 l2 l3 l4 l5 l6;
set Transfers := t0 t1;

# バンド幅
param bandwidth :=
 l0 20
 l1 30
 l2 25
 l3 32
 l4 27
 l5 19
 l6 23;


param uselink: t0 t1 :=
  l0 1 0
  l1 1 0
  l2 0 1
  l3 1 1
  l4 1 1
  l5 1 0
  l6 0 1;
# 行列を一致させないといけない

param weights:=
   t1 2
   t2 1;

実行方法

$ glpsol -m trans.mod -d trans.dat -o trans.out
これでtrans.outに結果が出ます。

実行結果

こんな感じで結果が出ます。
Problem:    trans
Rows:       8
Columns:    2
Non-zeros:  11
Status:     OPTIMAL
Objective:  THROUGHPUT = 46 (MAXimum)

   No.   Row name   St   Activity     Lower bound   Upper bound    Marginal
------ ------------ -- ------------- ------------- ------------- -------------
     1 THROUGHPUT   B             46
     2 BANDWIDTH[l0]
                    B             19                          20
     3 BANDWIDTH[l1]
                    B             19                          30
     4 BANDWIDTH[l2]
                    B              8                          25
     5 BANDWIDTH[l3]
                    B             27                          32
     6 BANDWIDTH[l4]
                    NU            27                          27             1
     7 BANDWIDTH[l5]
                    NU            19                          19             1
     8 BANDWIDTH[l6]
                    B              8                          23

   No. Column name  St   Activity     Lower bound   Upper bound    Marginal
------ ------------ -- ------------- ------------- ------------- -------------
     1 throughput[t0]
                    B             19             0
     2 throughput[t1]
                    B              8             0

t0には19Mbps, t1には8Mbps割り振ったらいいようです。

自分のCのプログラムに組み込んでみる

こうして線形計画法を解けるようになりましたが、自分で書いたプログラムに組み込むことで、もっと柔軟な問題を解けるようになります。そのためにGLPKには色々な関数とヘッダファイルが用意されていて、これをリンクすることで簡単に線形計画法を用いるプログラムを書くことができます。
例として、GLPK上に上げた問題を解くCのプログラムを示します。 このプログラムではデータをstruct transfer_problemという構造体に格納するものとして、プログラム中央付近のmake_sample_data()で問題を作成しています。そしてこの(勝手に定義した)データをGLPKに渡して解いているのがsolve()です。 lpx_simplex(lp)までがデータのロード、その後は解いたデータの取り出しです。
このプログラムをコンパイルするには、以下のようなオプションをつけてください。
$ gcc -I (glpkのインストール先)/include (glpkのインストール先)/lib/libglpk.a -lm -lglpk  glpk_sample.c
例えば以下のようなMakefileを用いるとよいです。
CFLAGS=-I /home/kay/local/glpk/include
LINK= /home/kay/local/glpk/lib/libglpk.a -lm -lglpk

trans: trans.o
	gcc $^ $(LINK) -o $@ 

.c.o:
	gcc $(CFLAGS) -c $^
/* glpk_sample.c */
#include <stdio.h>
#include <stdlib.h>
#include "glpk.h"


struct transfer_problem {
  int n_transfers;      // 
  int n_links;          // Number of links
  double *link_bws;     // Bandwidth for each link
  int *weights;         // Number of destination nodes in each transfer
  double *transfer_bws; // Bandwidth assigned to each transfer (output)
  
  // Connction matrix (Sparse expression, the index starts from 1)
  int *idx_links;
  int *idx_transfers;
  double *cm; 
  int n_cm_elems;
};


struct transfer_problem *new_transfer_problem(int n_transfers, int n_links, int n_cm_elems){
  struct transfer_problem *tp =\
    (struct transfer_problem *)malloc(sizeof(struct transfer_problem));
  
  tp->n_transfers = n_transfers;
  tp->n_links = n_links;
  tp->weights = (int *)malloc(sizeof(int) * n_transfers);
  tp->transfer_bws = (double *)malloc(sizeof(double) * n_transfers);
  tp->link_bws = (double *)malloc(sizeof(double) * n_links);

  // Connction matrix
  tp->n_cm_elems = n_cm_elems;
  tp->idx_links = (int *)malloc(sizeof(int) * (n_cm_elems+1));
  tp->idx_transfers = (int *)malloc(sizeof(int) * (n_cm_elems+1));
  tp->cm = (double *)malloc(sizeof(double) * (n_cm_elems+1));
  return tp;
}


struct transfer_problem *make_sample_data(void){
  // (link, transfer, value)
  //1 (1, 1, 1) 
  //2 (2, 1, 1)
  //3 (3, 2, 1)
  //4 (4, 1, 1)
  //5 (4, 2, 1)
  //6 (5, 1, 1)
  //7 (5, 2, 1)
  //8 (6, 1, 1)
  //9 (7, 2, 1)

  static int WEIGHTS[2]        = {2,1};
  static double LINK_BWS[7]    = {20, 30, 25, 32, 27, 19, 23};
  static int IDX_LINKS[9]     = {1, 2, 3, 4, 4, 5, 5, 6, 7};
  static int IDX_TRANSFERS[9] = {1, 1, 2, 1, 2, 1, 2, 1, 2};
  static int CM[9]            = {1, 1, 1, 1, 1, 1, 1, 1, 1};
  struct transfer_problem *tp;
  int i;

  tp = new_transfer_problem(2, 7, 9);

  for(i = 0; i < tp->n_transfers; i++)
    tp->weights[i] = WEIGHTS[i];
      
  for(i = 0; i < tp->n_links; i++)
    tp->link_bws[i] = LINK_BWS[i];

  for(i = 0; i < tp->n_cm_elems; i++){
    tp->idx_links[i+1] = IDX_LINKS[i];
    tp->idx_transfers[i+1] = IDX_TRANSFERS[i];
    tp->cm[i+1] = CM[i];
  }
  return tp;
}


double solve(struct transfer_problem *tp){
  int i;
  char label[128];
  LPX *lp;
  double Z;
  int cnt;

  lp = lpx_create_prob();
  lpx_set_prob_name(lp, "trans");
  lpx_set_obj_dir(lp, LPX_MAX);

  lpx_add_rows(lp, tp->n_links);
  for(i = 0; i < tp->n_links; i++){
    sprintf(label, "link_%d", i);
    lpx_set_row_name(lp, i+1, label);
    lpx_set_row_bnds(lp, i+1, LPX_UP, 0.0, tp->link_bws[i]);
  }

  lpx_add_cols(lp, tp->n_transfers);
  for(i = 0; i < tp->n_transfers; i++){
    sprintf(label, "transfer_%d", i);
    lpx_set_col_name(lp, i+1, label);
    lpx_set_col_bnds(lp, i+1, LPX_LO, 0.0, 0.0);
    lpx_set_obj_coef(lp, i+1, tp->weights[i]);
  }
  
  lpx_load_matrix(lp, tp->n_cm_elems, tp->idx_links, tp->idx_transfers, tp->cm);
  lpx_simplex(lp);
  Z = lpx_get_obj_val(lp);
  
  for(i = 0; i < tp->n_transfers; i++){
    tp->transfer_bws[i] = lpx_get_col_prim(lp, i+1);
  }

  lpx_delete_prob(lp);
  return Z;
}


int main(void){
  int i;
  double Z;
  struct transfer_problem *tp;
  
  tp = make_sample_data();
  Z = solve(tp);
  printf("Z = %5.3d\n", Z);
  for(i = 0; i < tp->n_links; i++){
    printf("link %d:%5.3f  ", i, tp->link_bws[i]);
  }
  printf("\n");
}