r - 使用 lpSolve 时的简并性

Question

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

请注意，当我使用预构建的 R 包(基准测试、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

Answer 1

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

来自lpSolve FAQ :

The main problems with these codes have to do with scaling, use of explicit inverses and lack of reinversion, and handling of degeneracy. Even small problems that are ill-conditioned or degenerate can bring most of these tableau codes to their knees.

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

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

如果您的目标只是摆脱简并性，请尝试扰动 rhs，尤其是您的目标函数。

如果您确实想找到两个系统不同的根源，我建议您自己编译 lpSolve *.c 文件(请参阅 here )并从 R 调用该版本以及从您的 Unix 集群查看结果是否仍然存在变化。

希望能帮助您前进。