![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
Rで数理計画
(例題など、どんどん付け足してください。)
ここでは、Rの関数を利用した数理計画問題の解法を説明します。 Rでも既存の関数を利用して、制約条件付き最適化問題を解くことができます。
LP。目的関数が1次(linear)で、制約条件が1次不等式で表現される数理計画問題。
パッケージ boot の関数 simplex() を利用します。
シンプレックス法を用いてLPを解く。linprogパッケージには，MPS形式のファイルを読みこむ，readMPS関数なども含まれている。 LPの解(最適解，最大(小)値)とともに，双対変数(シャドウバリュー)や感度分析の結果を得ることもできる。 linprog(線形計画/最適化)パッケージ中のオブジェクト一覧
ソルバーとして，lp_solveを利用している。整数計画法(および０‐１整数計画問題)にも対応*1。 lpSolve パッケージには割り当て問題や輸送問題の関数もある。lpSolve(線形・整数問題解法用 Lp_solve インターフェース)パッケージ中のオブジェクト一覧
機能は上と同じ？
glpk(Gnu Linear Programming Kit)へのインターフェイスを提供。
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 線形計画キットインターフェース)パッケージ中のオブジェクト一覧
subplex 最適化アルゴリズム
QP。目的関数が2次で制約条件が1次の不等式で表現される数理計画問題。
解法アルゴリズムには，Goldfarb-Idnani法が採用されているので狭義凸2次計画問題(いわゆる行列Qが正定値対称行列の場合)を解くために利用できます。
Weighted and lexicographical goal programming and optimization.
制約条件が線形不等式の場合には、パッケージ base の関数 constrOptim() も利用することができます(線形不等式制約付きの最適化関数 constrOptim 参照)。
orloca(オペレーションリサーチの立地分析モデルを扱う)パッケージ中のオブジェクト一覧
Optimization and Mathematical Programming
x1 + 9x2 + 3 x3
上の式の最大化
x1 + 2 x2 + 3 x3 <= 9 3 x1 + 2 x2 + 2 x3 <= 15
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
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
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)
#０です
tmp <- matrix(0,nrow=n,ncol=n)
#ショートなし制約のための行列
diag(tmp) <- 1
Amat <- cbind(rep(1,n),apply(return,2,mean),tmp)
#↑制約条件行列（ウエイトの合計が１／各銘柄の期待収益率／ショートなし）
bvec <- c(1.0, 0.030, rep(0,n))
#↑制約条件行列（ウエイトの合計が１／ポートフォリオの期待収益率３％／ショートなし）
sol <- solve.QP(Dmat,dvec,Amat,bvec=bvec,meq=2)
#二次計画法実行
sum(sol$sol)
#ウエイトの合計が１になっているかどうかの確認
t(matrix(sol$sol))%*%Dmat%*%matrix(sol$sol)
#リスク
t(matrix(sol$sol))%*%matrix(apply(return,2,mean))
#リターン（ちゃんと３％になってます）
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)}))
実行例（発送地が８, 目的地が５の場合）
>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
flontier.jpeg