2016-12-01 33 views
5

私は、値を最適化するためにいくつかのバイナリリソースから基本的にはナップザックの問題を選択しようとしていますが、線形プログラミングの問題があります。私が抱えている問題は、さまざまなリソースが共通の特性を持ち、最終的なソリューションが特定の特性を持つリソースが0または2かどうかを確認したいということです。これを達成するための方法はありますか?私は広範囲の調査にもかかわらず、1つを考えたり、1つを見つけることができませんでした。私のデータでは、決定変数はリソースであり、制約はそれらのリソースの特性です。次のコードを考えてみましょう:条件付き制約付きリニアプログラミングR

library(lpSolve) 
const_mat_so<-matrix(c(
    c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,1,0,0,1,0,1) 
    ,c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0,1,1,0,0,1,1) 
    ,c(0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1,0,1,0,1,0,0) 
    ,c(1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0,0,0,0,0,0,0) 
    ,c(8800, 8500, 7600, 8600, 8400, 7500, 7000, 8500, 8800, 7700, 6700,5500,1200,6700,9500,8700,6500) 
    ,c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,0,0,0,0,0,0) 
    ,c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,0,0,0,0,0,0) 
    ,c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0) 
    ,c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1,0,0,1,0,1,0) 
    ,c(0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0,1,1,0,0,0,0) 
    ,c(0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0,0,0,0,0,0,0) 
    ,c(0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0,1,1,1,0,1,0) 
    ,c(0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,0,0,0,0,1,0) 
    ,c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,0,0,0,1,0,0) 
    ,c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1,0,0,0,0,0,0) 
    ),nrow=15,byrow = TRUE) 

const_dir_so<-c("=","=","=","=","<=","<=","<=","<=","<=","<=","<=","<=","<=","<=","<=") 

max_cost_so = 25000 

objective_so = c(21.0, 19.3, 19.2, 18.8, 18.5, 16.6, 16.4, 16.4, 16.0, 16.0, 14.9, 14.6, 14.0, 13.9,12.0,5.5,24.6) 

const_rhs_so<-c(1,1,1,1,25000,3,3,3,2,2,2,2,2,2,2) 

x = lp ("max", objective_so, const_mat_so, const_dir_so, const_rhs_so,  all.bin=TRUE, all.int=TRUE 
    ) 

> x 
Success: the objective function is 68.1 

> x$solution 
[1] 1 0 1 0 0 0 0 0 0 0 0 0 1 1 0 0 0 

上記のソリューションを生成するが、私は実際に> = 2または0に最後の7つの制約をしたいので、それは私が望むソリューションではありません、私はこれをコーディングする方法を見当もつかないまたはそれが可能かどうか。どんな助けもありがとう。私は線形プログラミングの人ではないので、このアプローチに関する誤解を許してください。

+2

+1私はあなたの質問を理解していないという理由だけではなく、私は私の好奇心 – natario

答えて

2

私の理解は、最後の7つの制約のそれぞれは、1

1)あるようのみ7ような制約があり、2^7 =、ゼロ2以上であるとされていない、すなわちということです128の可能性は十分に小さいので、過剰な実行時間なしに疑問を抱かせた方程式を使用して1つ1つを実行し、その最大値を見つけることができます。

dec2binは、10進数の10進数を取り、それを0と1のバイナリベクトルに変換します。 0と127の間の各数値でそれを実行すると、1が> = 2(残りは0に等しい)の制約に対応するように2進数が与えられます。

dec2bin <- function(dec, digits = 7) { 
    # see http://stackoverflow.com/questions/6614283/converting-decimal-to-binary-in-r 
    tail(rev(as.integer(intToBits(dec))), digits) 
} 

runLP <- function(i) { 
    bin <- dec2bin(i) 
    n <- length(const_rhs_so) # 15 
    ix <- seq(to = n, length = length(bin)) # indexes of last 7 constraints, i.e. 9:15 
    const_dir_so[ix] <- ifelse(bin, ">=", "=") 
    const_rhs_so[ix] <- 2*bin 
    lp("max", objective_so, const_mat_so, const_dir_so, const_rhs_so, all.bin = TRUE) 
} 

lpout <- lapply(0:127, runLP) 
ixmax <- which.max(sapply(lpout, "[[", "objval")) 
ans <- lpout[[ixmax]] 
ans 
ans$solution 
tail(c(const_mat_so %*% ans$solution), 7) 

寄付:

> ans 
Success: the objective function is 62 
> ans$solution 
[1] 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 
> tail(c(const_mat_so %*% ans$solution), 7) # last 7 constraint values 
[1] 0 0 0 0 0 0 0 

2) @Erwin Kalvelagenの第二の代替では、制約の変数を参照するが、私はどのようなものだったことは彼の答えでxはのLHSの値であることだったと思います最後の7つの制約の1つ。すなわち、Cは、元の最後の7つの制約の行列である場合、これらの14個の制約を有するもの元7つの制約を置き換える、ある:

D1は対角要素の任意の十分に大きな負の数である対角行列であり、D2は
Cx + D1 y <= 0 
Cx + D2 y >= 0 

対角要素が全て-2である対角行列。ここでは、xyバイナリ変数のベクトルを最適化しています。 x変数は疑問と同じであり、最後の7個の元の制約のi番目の制約を0または1に制限して2以上に制約するために、y [i]が0であるような7個の新しいバイナリ変数があります。 yの変数は、(1)のbinと呼ばれます。目的変数内の変数yの係数はすべてゼロです。 lpSolve Rコードの観点から

:のように62の同一の値を与える

objective_so2 <- c(objective_so, numeric(7)) 
const_mat_so2 <- cbind(rbind(const_mat_so, const_mat_so[9:15, ]), 
         rbind(matrix(0, 8, 7), diag(-100, 7), diag(-2, 7))) 
const_dir_so2 <- c(const_dir_so, rep(">=", 7)) 
const_rhs_so2 <- c(const_rhs_so[1:8], numeric(14)) 
x2 = lp ("max", objective_so2, const_mat_so2, const_dir_so2, const_rhs_so2, all.bin = TRUE) 

(1)。 y変数(last7)はすべて0で、これも(1)に対応します。これはまた、2つのメソッドが一貫した回答を与えているので、二重チェックを提供します。

> x2 
Success: the objective function is 62 
> x2$solution 
[1] 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 
+0

本当に興味深いソリューションを満たすことができるので、私は誰が行う願うと回答。これらのタイプの制約の数が多くなるにつれ、パフォーマンスの問題はあるものの、機能します。どうもありがとう。これは決して私には起こりませんでした。 –

+0

downvote?もしそうなら、意図せず。いくぶん新しいスタックオーバーフロー –

+0

上記のすべては完全に機能します。この質問に答えた皆さん、ありがとう。大いに感謝します –

1

私はLpSolveが半連続変数をサポートしていると信じています。下限Lと上限Uを持つ半連続変数は、0またはLとUの間の値をとることができます。RパッケージlpSolveがこの変数型をサポートしているかどうかはわかりません。

ただし、余分なバイナリ変数yと余分な制約を使ってこれをシミュレートできます。だから、あなたは(あなたが整数値のみを希望する場合または整数)あなたのx変数連続作ると制約を追加する必要があります。

Uxの上限である
2*y <= x <= U*y 

+0

私はあなたを誤解しているかもしれませんが、制約は変数ではなく制約です。変数は自由に1の値をとることができますが、制約は1に等しくないようにする必要があります。したがって、これは半連続変数ではなく半連続制約です。私が誤解した場合は私を修正してください –

+0

それほど難しいことではありません。 'lhs = 0またはlhs> = 2'の代わりに、lhsが制約の左側にある場合は、半連続的な変数' x'を追加し、 'lhs-x = 0'として制約を書いてください。 LpSolveは半連続変数をサポートしているため、理論的にはこれが最善の実装となります。 –

1

lpSolveAPIパッケージは、 "lp_solve"のより高度なインターフェイスを提供します。 @Erwin Kalvelagenが述べたように、 "lp_solve"とlpSolveAPIは半連続変数をサポートしています(半連続の決定変数は、上限と下限の間で許容される値とゼロを取ります)。また、制約行列を使用すると、9〜15番目の制約式の出力を18〜24番目の変数に転送できます。例えば、(約9番目の制約)、x6 + x11 + x14 + x16 - x18 = 0,x6 + x11 + x14 + x16 = x18の場合。ですから、x18という半連続変数でx6 + x11 + x14 + x16を制御できると思います。

library(lpSolveAPI) 

    ## add 18-24th cols to define the 18-24th variables 
const_mat_so2 <- cbind(const_mat_so, rbind(matrix(0, nrow = 8, ncol = 7), diag(-1, 7))) 

    ## [EDITED] make a model and set a constraint matrix and objective coefs 
model <- make.lp(nrow(const_mat_so2), 0) 

for(i in 1:ncol(const_mat_so2)) add.column(model, const_mat_so2[,i]) 
set.constr.type(model, c(const_dir_so[-c(9:15)], rep("=", 7))) 
set.rhs(model, c(const_rhs_so[-c(9:15)], rep(0, 7))) # each original output - 18-24th = 0 

set.objfn(model, c(objective_so, rep(0, 7)))   # 18-24th are 0 

    ## define semi-continuous and bounds. 
set.semicont(model, col = 18:24) 
set.bounds(model, lower = rep(1.9, 7), col = 18:24) # default upper is Inf. 

    ## define other things 
set.type(model, col = 1:17, type = "binary")  # original variable 
set.type(model, col = 18:24, type = "integer") # outputs of original constraint formulas 
lp.control(model, sense = "max")     # do maximize 

# write.lp(model, "filename.lp", "lp") # if you want to watch the whole model 
solve(model) 
get.variables(model) 
# [1] 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 [18] 0 0 0 0 0 0 0 
get.objective(model) 
# [1] 62 
t(const_mat_so %*% res[1:17]) 
#  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] 
# [1,] 1 1 1 1 22300 1 0 0 0  0  0  0  0  0  0