2016-12-01 22 views
5

我有一個線性規劃問題,我試圖從大量的二進制資源中選擇優化值,基本上是一個揹包問題。我遇到的問題是不同的資源具有共同的特徵,我想確保我的最終解決方案具有0或2個具有特定特徵的資源。有什麼辦法可以做到這一點?儘管大量搜索,我還是無法想象或找到一個。在我的數據中,決策變量是資源,約束是這些資源的特徵。請看下面的代碼:有條件限制的線性規劃在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。我不知道如何編寫代碼這個或者是否有可能。任何幫助,將不勝感激。我不是一個線性編程專家,所以請原諒關於這種方法的任何誤解。

+2

+1只是因爲我不明白你的問題,但我希望有人做和答案,所以我可以滿足我的好奇心 – natario

回答

2

我理解的是每一個的最後7個約束爲大於2或等於零,即不1.

1)只有7這樣的約束條件所以有2^7 = 128個可能性足夠小,以至於我們可以運用所提供的公式來運行每個程序,而無需過多運行時間,然後找到最大程度的這些可能性。

dec2bin以10爲底數(即十進制數)並將其轉換爲0和1的二進制向量。在0到127之間的每個數字上運行它會給出二進制數,例如1表示約束> = 2(其餘等於0)。

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個約束之一。也就是說,如果C是原始最後7個約束矩陣,則這些14個約束替換那些原來7個約束:

Cx + D1 y <= 0 
Cx + D2 y >= 0 

其中,D1是一個對角矩陣,其對角元素是任何足夠大的負數和D2是一個對角線元素都是-2的對角矩陣。在這裏,我們正在優化xy二元變量的向量。 x變量與問題中一樣,並且有7個新的y二元變量,例如y [i]爲0,以將最後7個原始約束的第i個約束限制爲0或1,以將其約束爲2或更多。 (1)中的y變量被稱爲bin。目標變量的係數全部爲零。

在lpSolve R代碼裏面術語:

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) 

給予的62相同的值(1)中。 y變量(最後7)都是0,這也對應於(1)。這也提供了一個雙重檢查,因爲現在兩種方法已經給出一致的答案。

> 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連續可變(或整數,如果你只想要整數值),並添加約束:

2*y <= x <= U*y 

其中U是一個上限x

+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 = 0x6 + 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