//佐藤(2003-11-5)
COLOR(red){SIZE(30){Rで数理計画}}
(例題など、どんどん付け足してください。)
ここでは、Rの関数を利用した数理計画問題の解法を説明します。
Rでも既存の関数を利用して、制約条件付き最適化問題を解くことができます。
// "*", "**", "***" による段落に応じて、ハイパーリンク付きの目次一覧をこの位置に置く
// 行頭が "//" で始まる行はコメントで表示されません。
#contents
// SIZE(num){テキスト} はテキスト文字の大きさを指定します
// COLOR{color}{テキスト} はテキスト文字の色を指定します
// 行頭が "**" で始まると、直前の "*" で始まる段落の副段落として整形されます
// 二重括弧 "((...))" で囲ったテキストは注釈になります。注釈はページ末尾に置かれ、ハイパーリンクが張られます
// 行頭に半角空白を置くとこの通り枠で囲って見易く整形されます
// 行頭に半角空白がある引き続く行は一まとめに整形されます
// コードを書くのに便利です
// 行頭が "***" で始まると、直前の "**" で始まる副段落の副段落として整形されます
// 行頭に "----" を置くとその場所に水平線が引かれます
----
// 行頭に "*" を置くと、段落見出しとして整形されます
*線形計画法/整数計画法 [#k79bb2f1]
LP。目的関数が1次(linear)で、制約条件が1次不等式で表現される数理計画問題。
**simplex 関数(bootパッケージ) [#bc9dc2f0]
パッケージ boot の関数 simplex() を利用します。
**solveLP 関数(linprogパッケージ) [#h57b4eda]
シンプレックス法を用いてLPを解く。linprogパッケージには，MPS形式のファイルを読みこむ，readMPS関数なども含まれている。
LPの解(最適解，最大(小)値)とともに，双対変数(シャドウバリュー)や感度分析の結果を得ることもできる。
[[linprog(線形計画/最適化)パッケージ中のオブジェクト一覧]]
**lp関数(lpSolve パッケージ) [#bb6b975e]
ソルバーとして，[[lp_solve:ftp://ftp.es.ele.tue.nl/pub/lp_solve/]]を利用している。整数計画法(および０‐１整数計画問題)にも対応((lp 関数の引数 int.vec で設定))。
lpSolve パッケージには割り当て問題や輸送問題の関数もある。[[lpSolve(線形・整数問題解法用 Lp_solve インターフェース)パッケージ中のオブジェクト一覧]]
**lpSolveAPI(Interface to lp_solve v. 5.5) [#ac36e411]
機能は上と同じ？
**glpk パッケージ [#ua7d5a74]
glpk([[Gnu Linear Programming Kit:http://www.gnu.org/software/glpk/glpk.html]])へのインターフェイスを提供。
**Rglpk パッケージ [#vb3ccfce]
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 パッケージ [#m37949f4]
subplex 最適化アルゴリズム
**Rcplex パッケージ [#ke76f81e]
**DEA パッケージ 包絡分析 [#f9e9486c]
**[[FEAR: Frontier Efficiency Analysis with R:http://www.clemson.edu/economics/faculty/wilson/Software/FEAR/fear.html]] [#uc2e1e93]
**limSolve: Solving Linear Inverse Models [#u4a31a53]
*2次計画法 [#y928225b]
QP。目的関数が2次で制約条件が1次の不等式で表現される数理計画問題。
**パッケージ quadprog の関数 solve.QP() を利用します。 [#a3ff23e0]
解法アルゴリズムには，Goldfarb-Idnani法が採用されているので狭義凸2次計画問題(いわゆる行列Qが正定値対称行列の場合)を解くために利用できます。
**パッケージ kernlab の ipop 関数 2次計画法(Quadratic Programming) ソルバー [#u197f44e]
-Karatzoglou, Alexandros; Smola, Alex; Hornik, Kurt; Zeileis, Achim(2004): [[kernlab - An S4 package for kernel methods in R. Research Report Series / Department of Statistics and Mathematics, Nr. 9, Institut fur Statistik und Mathematik :http://epub.wu-wien.ac.at/dyn/virlib/wp/showentry?ID=epub-wu-01_766]] 解説あり
**Rcplex パッケージ [#s24bd10c]
**LowRankQP: Low Rank Quadratic Programming [#y67ab392]
*非線形計画法 [#c9b75dc6]
**BB: Solving and Optimizing Large-Scale Nonlinear Systems [#g0ace6d0]
*目標計画法 [#bad8517f]
**goalprog パッケージ [#cabc2788]
Weighted and lexicographical goal programming and optimization.
*線形および2次混合線形計画 [#hc5b9c84]
**Rcplex パッケージ [#p4e41dba]
*線形不等式制約条件付き最適化問題 [#o850bad6]
制約条件が線形不等式の場合には、パッケージ base の関数 constrOptim()
も利用することができます([[線形不等式制約付きの最適化関数 constrOptim]] 参照)。
*割り当て問題 [#veb7971d]
**lp.assign.R 関数(lpSolveパッケージ) [#h1f30602]
*輸送問題 [#jbb54a76]
**lp.transport.R 関数(lpSolveパッケージ) [#i02a497a]
*立地配分 [#d4f99855]
**orloca: オペレーションリサーチの立地分析モデルを扱う [#z6b3ff04]
[[orloca(オペレーションリサーチの立地分析モデルを扱う)パッケージ中のオブジェクト一覧]]
*滑らかな非線形関数の最適化 [#rcc64f5b]
-[[Rdonlp2>http://arumat.net/Rdonlp2/]]
-制約条件が非線形関数の場合に使えます．目的関数はもちろん非線形ＯＫ
-まだ、正式なパッケージではありません．
-[[「donlp2」を使用した最適化ルーティン]]
*その他 [#k251f292]
-MSVAR: Markov Switching VAR
C言語やFORTRANのルーチンを利用。
optim()なども少し関連。
-nleqslv: 非線形連立方程式を解く
*CRAN TASK VIEWS [#w73a431d]
[[Optimization and Mathematical Programming:http://cran.md.tsukuba.ac.jp/web/views/Optimization.html]]
----
*例題1(LP/IP/0-1IP) [#b529ec14]
**線形計画問題 [#rd841374]
- 目的関数
x1 + 9x2 + 3 x3
上の式の最大化
- 制約式
x1 + 2 x2 + 3 x3 <= 9
3 x1 + 2 x2 + 2 x3 <= 15
*** lpSolve [#j9808d5e]
- 目的関数
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)
**整数計画問題 [#vf8561c4]
- 目的関数
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 [#cc414bd9]
- 目的関数
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整数計画問題 [#f9a1e157]
- 目的関数
- 制約式
*** lpSolve [#h4efd094]
*例題2(QP) [#n067e010]
**Rで有効フロンティア [#jd637cf2]
-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)
#０です
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)}))
-ご注意
--結果は毎回違います
--nが1001より大きいと上記のバージョンでは動きません。
--Dmatがランク落ちしていたのが原因。
--乱数の種が1000以上に対応していないのかも。。。そのため共分散行列（Dmat）がランク落ち？
--実際のリターン系列を使用して、やってみたいですね。
-結果の有効フロンティアの例
#Ref(flontier.jpeg)
--補足ですが、tseriesパッケージにsolve.QPを利用したportfolio.optimという
関数があります(by Y.)
--補足ありがとうございます (。^_^。)／追加補足です。「portfolio.optim()」のウエイトの制約条件は、?ショートするかしないか?ウエイトの上限下限、を指定できるようになっています（もちろんウエイトの合計が１になるという条件はコーディングされているはず、コード見てないのでなんともいえませんが...）。一方で、たとえば４資産のポートフォリオで、あるＡ資産とＢ資産の合計が３５％を超えないようにするという制約条件は指定できません。そんな場合は「solve.QP()」を使用してAmatとbvecを指定する必要があります。（違いますでしょうか。。。）(byどうむ)
--で、個人的な問題意識から[[こんなページ>「donlp2」を使用した最適化ルーティン]]を作成してみましたのでみなさん（Y.さんもよろしく^^）ぜひ情報をお願いいたします（←これが言いたかった(^_^;)）。(byどうむ)
*例題3(AP) [#f0691db0]
**割り当て問題 [#a839cbaa]
*例題4(TP) [#w67361a9]
**輸送問題 [#xe0e78ec]
***lpSolve [#ta940948]
実行例（発送地が８, 目的地が５の場合）
-輸送費用
>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
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
