Rで数理計画

(例題など、どんどん付け足してください。)

ここでは、Rの関数を利用した数理計画問題の解法を説明します。 Rでも既存の関数を利用して、制約条件付き最適化問題を解くことができます。


線形計画法/整数計画法

LP。目的関数が1次(linear)で、制約条件が1次不等式で表現される数理計画問題。

simplex 関数(bootパッケージ)

パッケージ boot の関数 simplex() を利用します。

solveLP 関数(linprogパッケージ)

シンプレックス法を用いてLPを解く。linprogパッケージには,MPS形式のファイルを読みこむ,readMPS関数なども含まれている。 LPの解(最適解,最大(小)値)とともに,双対変数(シャドウバリュー)や感度分析の結果を得ることもできる。 linprog(線形計画/最適化)パッケージ中のオブジェクト一覧?

lp関数(lpSolve パッケージ)

ソルバーとして,lp_solveを利用している。整数計画法(および0‐1整数計画問題)にも対応*1。 lpSolve パッケージには割り当て問題や輸送問題の関数もある。lpSolve(線形・整数問題解法用 Lp_solve インターフェース)パッケージ中のオブジェクト一覧?

lpSolveAPI(Interface to lp_solve v. 5.5)

機能は上と同じ?

glpk パッケージ

glpk(Gnu Linear Programming Kit)へのインターフェイスを提供。

Rglpk パッケージ

R interface to the GNU Linear Programing Kit (GLPK version 4.25). GLPK is open source software for solving large-scale linear programming (LP), mixed integer linear programming (MILP) and other related problems. Rglpk(R/GNU 線形計画キットインターフェース)パッケージ中のオブジェクト一覧?

sublex パッケージ

subplex 最適化アルゴリズム

Rcplex パッケージ

DEA パッケージ 包絡分析

FEAR: Frontier Efficiency Analysis with R

limSolve: Solving Linear Inverse Models

2次計画法

QP。目的関数が2次で制約条件が1次の不等式で表現される数理計画問題。

パッケージ quadprog の関数 solve.QP() を利用します。

解法アルゴリズムには,Goldfarb-Idnani法が採用されているので狭義凸2次計画問題(いわゆる行列Qが正定値対称行列の場合)を解くために利用できます。

パッケージ kernlab の ipop 関数 2次計画法(Quadratic Programming) ソルバー

Rcplex パッケージ

LowRankQP: Low Rank Quadratic Programming

非線形計画法

BB: Solving and Optimizing Large-Scale Nonlinear Systems

目標計画法

goalprog パッケージ

Weighted and lexicographical goal programming and optimization.

線形および2次混合線形計画

Rcplex パッケージ

線形不等式制約条件付き最適化問題

制約条件が線形不等式の場合には、パッケージ base の関数 constrOptim() も利用することができます(線形不等式制約付きの最適化関数 constrOptim 参照)。

割り当て問題

lp.assign.R 関数(lpSolveパッケージ)

輸送問題

lp.transport.R 関数(lpSolveパッケージ)

立地配分

orloca: オペレーションリサーチの立地分析モデルを扱う

orloca(オペレーションリサーチの立地分析モデルを扱う)パッケージ中のオブジェクト一覧?

滑らかな非線形関数の最適化

その他

  • MSVAR: Markov Switching VAR C言語やFORTRANのルーチンを利用。 optim()なども少し関連。
  • nleqslv: 非線形連立方程式を解く

CRAN TASK VIEWS

Optimization and Mathematical Programming


例題1(LP/IP/0-1IP)

線形計画問題

  • 目的関数
x1 + 9x2 + 3 x3

 上の式の最大化

  • 制約式
x1 + 2 x2 + 3 x3  <= 9
3 x1 + 2 x2 + 2 x3 <= 15

lpSolve

  • 目的関数
f.obj <- c(1, 9, 3) 
  • 制約式
f.con <- matrix (c(1, 2, 3, 3, 2, 2), nrow=2, byrow=TRUE)
f.dir <- c("<=", "<=")
f.rhs <- c(9, 15)

整数計画問題

  • 目的関数
x1 + 9 x2 + 3 x3

 上の式の最大化

  • 制約式
x1 + 2 x2 + 3 x3  <= 9
3 x1 + 2 x2 + 2 x3 <= 15
  • 計算実行
lp ("max", f.obj, f.con, f.dir, f.rhs)
  • 解の表示
lp ("max", f.obj, f.con, f.dir, f.rhs)$solution

lpSolve

  • 目的関数
f.obj <- c(1, 9, 3) 
  • 制約式
f.con <- matrix (c(1, 2, 3, 3, 2, 2), nrow=2, byrow=TRUE)
f.dir <- c("<=", "<=")
f.rhs <- c(9, 15)
  • 計算実行
lp ("max", f.obj, f.con, f.dir, f.rhs, int.vec=1:3)
  • 解の表示
lp ("max", f.obj, f.con, f.dir, f.rhs)$solution

0-1整数計画問題

  • 目的関数
  • 制約式

lpSolve

例題2(QP)

Rで有効フロンティア

  • Rで金融の世界における有効フロンティアを描画してみました(byどうむ)
  • quadprogを使用しています。
  • ↓ソースコード(WINDOWS版、Xwindow版両方での動作を確認)
    library(quadprog)
        #二次計画法ルーティンのロード
    n <- 100
        #銘柄数(自由に変更してください)
    return <- matrix(rnorm(1000*n,0.03,0.30), ncol=n)
    return[,n] <- return[,n]+0.05
        #共分散行列を計算するためにリターンの系列を乱数で発生させてます
    
    ##################################################
    # n変数の有効フロンティア(リターン条件あり/ショートなし)
    # min(-d^T b + 1/2 b^T D b)
    # with the constraints A^T b >= b_0. 
    # リターン一定で最小のリスクを算出しています(←ちょっと説明わかりにくい)
    ##################################################
    Dmat <- cov(return)
        #共分散行列
    dvec <- rep(0,n)
        #0です
    tmp <- matrix(0,nrow=n,ncol=n)
        #ショートなし制約のための行列
    diag(tmp) <- 1
    Amat <- cbind(rep(1,n),apply(return,2,mean),tmp)
        #↑制約条件行列(ウエイトの合計が1/各銘柄の期待収益率/ショートなし)
    bvec <- c(1.0, 0.030, rep(0,n))
        #↑制約条件行列(ウエイトの合計が1/ポートフォリオの期待収益率3%/ショートなし)
    sol <- solve.QP(Dmat,dvec,Amat,bvec=bvec,meq=2)
        #二次計画法実行
    sum(sol$sol)
        #ウエイトの合計が1になっているかどうかの確認
    t(matrix(sol$sol))%*%Dmat%*%matrix(sol$sol)
        #リスク
    t(matrix(sol$sol))%*%matrix(apply(return,2,mean))
        #リターン(ちゃんと3%になってます)
    sol <- list(sol)
        #上の結果をリスト型として保存
    
    for(iret in seq(0.03,0.05,by=0.0005)){
    	#リターンを3%から5%までうごかす
    	bvec[2] <- iret
    	sol <- c(sol,list(solve.QP(Dmat,dvec,Amat,bvec=bvec,meq=2)))
    	    #計算結果を変数「sol」に保存
    }
    RetSigma<- 
    sapply(sol,
    	function(x){
    		c(
    			sigma=sqrt(t(matrix(x$sol))%*%Dmat%*%matrix(x$sol)),
    			ret=t(matrix(x$sol))%*%matrix(apply(return,2,mean))
    		)
    	}
    )
        #↑結果を元にリターンとリスクを計算
    RetSigma <- t(RetSigma)
    plot(RetSigma)
        #↑表示
    #system.time(apply(matrix(1:100),1,function(x){sol <-
                solve.QP(Dmat,dvec,Amat,bvec=bvec,meq=2)}))
  • ご注意
    • 結果は毎回違います
    • nが1001より大きいと上記のバージョンでは動きません。
    • Dmatがランク落ちしていたのが原因。
    • 乱数の種が1000以上に対応していないのかも。。。そのため共分散行列(Dmat)がランク落ち?
    • 実際のリターン系列を使用して、やってみたいですね。
  • 結果の有効フロンティアの例
    flontier.jpeg
    • 補足ですが、tseriesパッケージにsolve.QPを利用したportfolio.optimという 関数があります(by Y.)
    • 補足ありがとうございます (。^_^。)/追加補足です。「portfolio.optim()」のウエイトの制約条件は、?ショートするかしないか?ウエイトの上限下限、を指定できるようになっています(もちろんウエイトの合計が1になるという条件はコーディングされているはず、コード見てないのでなんともいえませんが...)。一方で、たとえば4資産のポートフォリオで、あるA資産とB資産の合計が35%を超えないようにするという制約条件は指定できません。そんな場合は「solve.QP()」を使用してAmatとbvecを指定する必要があります。(違いますでしょうか。。。)(byどうむ)
    • で、個人的な問題意識からこんなページを作成してみましたのでみなさん(Y.さんもよろしく^^)ぜひ情報をお願いいたします(←これが言いたかった(^_^;))。(byどうむ)

例題3(AP)

割り当て問題

例題4(TP)

輸送問題

lpSolve

 実行例(発送地が8, 目的地が5の場合)

  • 輸送費用
>costs <- matrix (100000, 8, 5);#最大輸送費用 10,000
>costs[4,1] <- costs[-4,5] <- 0
>costs[1,2] <- costs[2,3] <- costs[3,4] <- 7; costs[1,3] <- costs[2,4] <- 7.7
>costs[5,1] <- costs[7,3] <- 8; costs[1,4] <- 8.4; costs[6,2] <- 9
>costs[8,4] <- 10; costs[4,2:4] <- c(.7, 1.4, 2.1)
  • 発送地
>row.signs <- rep ("<", 8)
>row.rhs <- c(200, 300, 350, 200, 100, 50, 100, 150)
  • 目的地
>col.signs <- rep (">", 5)
>col.rhs <- c(250, 100, 400, 500, 200)
  • 実行
>lp.transport (costs, row.signs, row.rhs, col.signs, col.rhs)
  • 計算結果(最小の輸送費用)
>lp.transport (costs, row.signs, row.rhs, col.signs, col.rhs)$solution
         [,1] [,2] [,3] [,4] [,5]
    [1,]    0  100    0  100    0
    [2,]    0    0  300    0    0
    [3,]    0    0    0  350    0
    [4,]  200    0    0    0    0
    [5,]   50    0    0    0   50
    [6,]    0    0    0    0   50
    [7,]    0    0  100    0    0
    [8,]    0    0    0   50  100

*1 lp 関数の引数 int.vec で設定

添付ファイル: fileflontier.jpeg 3061件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Google
WWW を検索 OKADAJP.ORG を検索
Last-modified: 2015-03-01 (日) 01:15:59 (1514d)