R：等滲迴歸最小化

Question

我想盡量減少以下等式：

F=SUM{u 1:20}sum{w 1:10} Quw(ruw-yuw) 

有以下限制：

yuw >= yu,w+1 
yuw >= yu-1,w 
y20,0 >= 100 
y0,10 >= 0 

我有一個20×10 RUW和20 * 10 QUW矩陣，我現在需要生成一個遵守約束的yuw矩陣。我在R編碼，熟悉lpsolve和optimx軟件包，但不知道如何將它們用於這個特定的問題。

Answer 1

因爲Quw和ruw都是數據，所以在yuw決策變量中，所有約束以及目標都是線性的。因此，這是一個線性編程問題，可以通過lpSolve包解決。

爲了抽象出一點，我們假設R=20和C=10描述了輸入矩陣的維數。然後有R*C決策變量，我們可以將它們分配的順序爲y11, y21, ... yR1, y12, y22, ... yR2, ..., y1C, y2C, ..., yRC，讀取變量矩陣的列。

目標中每個yuw變量的係數爲-Quw;請注意，求和中的Quw*ruw項是常數（不受我們爲決策變量選擇的值的影響），因此不會輸入到線性規劃求解器。有趣的是，這意味着ruw實際上對優化模型解決方案沒有影響。

第一個R*(C-1)約束對應於yuw >= yu,w+1約束，並且下一個(R-1)*C約束對應於約束yuw >= yu-1,w。最後兩個約束對應於y20,1 >= 100和y1,10 >= 0約束。

我們可以用下面的R命令將這個模型輸入到lpsolve包中，輸入一個簡單的Q矩陣，其中每個條目都是-1（最終的解決方案應該將所有決策變量設置爲0，除了左下角之外，這應該是100）：

# Sample data 
Quw <- matrix(-1, nrow=20, ncol=10) 
R <- nrow(Quw) 
C <- ncol(Quw) 

# Build constraint matrix 
part1 <- matrix(0, nrow=R*(C-1), ncol=R*C) 
part1[cbind(1:(R*C-R), 1:(R*C-R))] <- 1 
part1[cbind(1:(R*C-R), (R+1):(R*C))] <- -1 
part2 <- matrix(0, nrow=(R-1)*C, ncol=R*C) 
pos2 <- as.vector(sapply(2:R, function(r) r+seq(0, R*(C-1), by=R))) 
part2[cbind(1:nrow(part2), pos2)] <- 1 
part2[cbind(1:nrow(part2), pos2-1)] <- -1 
part3 <- rep(0, R*C) 
part3[R] <- 1 
part4 <- rep(0, R*C) 
part4[(C-1)*R + 1] <- 1 
const.mat <- rbind(part1, part2, part3, part4) 

library(lpSolve) 
mod <- lp(direction = "min", 
      objective.in = as.vector(-Quw), 
      const.mat = const.mat, 
      const.dir = rep(">=", nrow(const.mat)), 
      const.rhs = c(rep(0, nrow(const.mat)-2), 100, 0)) 

我們現在可以訪問模型解決方案：

mod 
# Success: the objective function is 100 
matrix(mod$solution, nrow=R) 
#  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] 
# [1,] 0 0 0 0 0 0 0 0 0  0 
# [2,] 0 0 0 0 0 0 0 0 0  0 
# [3,] 0 0 0 0 0 0 0 0 0  0 
# [4,] 0 0 0 0 0 0 0 0 0  0 
# [5,] 0 0 0 0 0 0 0 0 0  0 
# [6,] 0 0 0 0 0 0 0 0 0  0 
# [7,] 0 0 0 0 0 0 0 0 0  0 
# [8,] 0 0 0 0 0 0 0 0 0  0 
# [9,] 0 0 0 0 0 0 0 0 0  0 
# [10,] 0 0 0 0 0 0 0 0 0  0 
# [11,] 0 0 0 0 0 0 0 0 0  0 
# [12,] 0 0 0 0 0 0 0 0 0  0 
# [13,] 0 0 0 0 0 0 0 0 0  0 
# [14,] 0 0 0 0 0 0 0 0 0  0 
# [15,] 0 0 0 0 0 0 0 0 0  0 
# [16,] 0 0 0 0 0 0 0 0 0  0 
# [17,] 0 0 0 0 0 0 0 0 0  0 
# [18,] 0 0 0 0 0 0 0 0 0  0 
# [19,] 0 0 0 0 0 0 0 0 0  0 
# [20,] 100 0 0 0 0 0 0 0 0  0 

請注意，你的模型很容易變得不可行，如果改變Quw（例如，如果我們用1我填這而不是-1）。在這些情況下，模型將以狀態3退出（您可以通過運行mod或mod$status來查看）。