3

我在 R 中使用 lpSolve。我的模型(数据包络分析)在我的 MAC 上运行良好,但是当我尝试在 UNIX 集群上运行它时,发现许多模型是退化的。两个系统上的 lp.control 选项相同。通过使用 presolve 和 anti.degen 选项,我能够减少但不能消除退化的数量。

请注意,当我使用预构建的 R 包(Benchmarking、nonparaeff)来解决相同的线性规划模型时,会出现完全相同的问题。

有谁知道为什么 UNIX 集群上的退化错误?

干杯,

彼得

如果有人感兴趣,代码如下。基本上,它为三百个代理中的每一个创建了一个线性规划模型并解决了这个问题。在我的 MAC 上,所有问题都解决得很好,但在集群上,90% 的问题被发现是退化的:

library(lpSolveAPI)
set.seed(198302)

##############Create data
x=matrix(rnorm(1200,5,3),300,4)
y1=x%*%c(.4,.2,.7,.8)+rnorm(300,4,.5)
y2=x%*%c(.5,.8,.2,.3)+rnorm(300,4,.5)
y=cbind(y1,y2)

##############Write DEA function
xref=x
yref=y

##Define dimensions
mx<-ncol(xref)
my<-ncol(yref)
nref<-nrow(xref)
nobs<-nrow(x)

##Define empty matrices for efficiency scores, lambdas and slacks
eff<-rep(0,nobs) 
lambda<-matrix(0,nobs,nref)
slacks<-matrix(0,nobs,mx)

##Start the model, noting that there will be as many constraints as there are inputs, outputs and one additional constraint on lambda.  And as many decision variables as there are producers (lambdas) + one (efficiency score)
deamodel<-make.lp(nrow=mx+my+1,ncol=nref+1)
    ## For each input and output set the row as a constraint
    for (h in 1:mx) set.row(deamodel, h, c(0, -xref[, h]))
    for (h in 1:my) set.row(deamodel, mx + h, c(0, yref[, h]))
    ##Set a single objective function
    set.objfn(deamodel, 1, 1)
    ##Set the type of constraints.  Inputs and outputs are all greater than, lambdas is equal to
    set.constr.type(deamodel, c(rep(">=", mx + my),"="))
    ##Set another row as a constraint on lambdas
    set.row(deamodel, (mx+my+1), c(0,rep(1,nref)))
    set.rhs(deamodel, 1, mx+my+1)

##In a loop create a lp model for each producer
for (k in 1:nobs){
    ##Set the right hand side equal to output of the individual producer 
    set.rhs(deamodel, c(rep(0,mx), y[k, ]), 1:(mx + my))
    ##Set first column equal to input of producer
    set.column(deamodel, 1, c(1,x[k,]), 0:mx)
    ##Set some presolve options
    lp.control(deamodel, presolve=c("rows", "columns","rowdominate","coldominate","bounds"))
    ##Solve the model
    a=solve(deamodel)
    if (a==0){
        eff[k]<-get.objective(deamodel)
        lambda[k,]<-get.variables(deamodel)[-1]
        slacks[k,]<-get.constraints(deamodel)[1:mx]}
    if (a!=0){
        eff[k]<-NA
        lambda[k,]<-NA
        slacks[k,]<-NA
    }}

eff
4

1 回答 1

1

看起来您正在解决 300 个本质上高度退化的小问题(301 个变量,7 个约束)。Presolve 和 anti.degen 只能带您到此为止。

lpSolve 常见问题解答

这些代码的主要问题与缩放、使用显式逆和缺乏逆反以及退化处理有关。即使是病态或退化的小问题也会使这些画面代码中的大部分都屈服。

您的 Unix 集群上的 lpSolve 实现(识别退化的解决方案)似乎与您的 Mac(从 R 调用)正在调用的不同。

第一次测试:写出 300 个 MIP ( write.lp) 并查看它们在 R 和您的 unix 集群中是否相同。(您正在使用rnorm,即使是非常小的舍入也可以摆脱其中一些高度退化的问题。)

如果您的目标只是摆脱退化,请尝试扰动 rhs,尤其是您的目标函数。

如果您真的想了解这两个系统为何不同的根源,我建议您自己编译 lpSolve *.c 文件(参见此处)并从 R 以及您的 Unix 集群中调用该版本以查看是否还有结果的变化。

希望能帮助你前进。

于 2013-06-25T19:22:37.307 回答