我在 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