如何更改 lpSolve 中的重量限制？答案

【问题标题】：How to change weight restrictions in lpSolve?如何更改 lpSolve 中的重量限制？
【发布时间】：2015-09-26 13:50:33
【问题描述】：

我正在应用 CCR 数据包络分析模型对股票数据进行基准测试。为此，我正在运行发布在here 的 DEA 论文中的 R 代码。本文档使用 lpSolve 解决线性问题。 lpSolve 的文档是here。

本文档附带有关如何在 R 中实现以下模型的分步说明。

背景信息：

数据包络分析 (=DEA) 是一种面向数据的方法，用于评估称为决策单元的一组实体的性能，这些实体将多个数据输入转换为多个数据输出。DEA 用于评估和比较决策单元的性能(DMU) 通过考虑生产单元与效率边界的接近程度来确定生产单元的相对效率。

我们假设有 n 个 DMU 需要评估。每个 DMU 消耗不同数量的 ????不同的投入生产？？？不同的输出。具体来说？？？？？？？？？？？？？？？？消耗量????????????输入？？？并产生数量????????????输出？？？？？？我们进一步假设????????????，???????????? ≥ 0 且每个 DMU 至少有一个正输入和一个正输出值。在 CCR 模型中，目标是确定权重 ????,????使用线性规划来最大化每个决策单元的复合效率。这是由比率定义的：

然后变成这样的线性问题：

应用：

我正在将此效率模型应用于股票数据 CAPM 风险因素：

这是一个 5 个输入 - 1 个输出的情况

这里的目标是为分析的每个资产找到“最佳”权重集（????1、????1…????5）。此处使用的术语“最佳”表示当将这些权重分配给所有 DMU / Stocks (1...10) 的输入和输出时，上述比率相对于所有其他比率最大化。

这是经过标准化的原始样本 z 得分数据，所有值都介于 0 和 1 之间。

> dput(testdfst)
structure(list(Name = structure(1:10, .Label = c("Stock1", "Stock2", 
"Stock3", "Stock4", "Stock5", "Stock6", "Stock7", "Stock8", "Stock9", 
"Stock10"), class = "factor"), Date = structure(c(14917, 14917, 
14917, 14917, 14917, 14917, 14917, 14917, 14917, 14917), class = "Date"), 
    `(Intercept)` = c(0.454991569278089, 1, 0, 0.459437188169979, 
    0.520523252955415, 0.827294243132907, 0.642696631099892, 
    0.166219881886161, 0.086341470900152, 0.882092217743293), 
    rmrf = c(0.373075150411683, 0.0349067218712968, 0.895550280607866, 
    1, 0.180151549474574, 0.28669170468735, 0.0939821798173586, 
    0, 0.269645291515763, 0.0900619760898984), smb = c(0.764987877309785, 
    0.509094491489323, 0.933653313048327, 0.355340700554647, 
    0.654000372286503, 1, 0, 0.221454091364611, 0.660571586102851, 
    0.545086931342479), hml = c(0.100608151187926, 0.155064872867367, 
    1, 0.464298576152336, 0.110803875258027, 0.0720803195598597, 
    0, 0.132407005239869, 0.059742053684015, 0.0661623383303703
    ), rmw = c(0.544512524466665, 0.0761995312858816, 1, 0, 0.507699534880555, 
    0.590607506295898, 0.460148690870041, 0.451871218073951, 
    0.801698199214685, 0.429094840372901), cma = c(0.671162426988512, 
    0.658898571758625, 0, 0.695830176886926, 0.567814542084284, 
    0.942862571603074, 1, 0.37571611336359, 0.72565234813082, 
    0.636762557753099), Returns = c(0.601347600017365, 0.806071701848376, 
    0.187500487065719, 0.602971876359073, 0.470386289298666, 
    0.655773224143057, 0.414258177255333, 0, 0.266112191477882, 
    1)), .Names = c("Name", "Date", "(Intercept)", "rmrf", "smb", 
"hml", "rmw", "cma", "Returns"), row.names = c("Stock1.2010-11-04", 
"Stock2.2010-11-04", "Stock3.2010-11-04", "Stock4.2010-11-04", 
"Stock5.2010-11-04", "Stock6.2010-11-04", "Stock7.2010-11-04", 
"Stock8.2010-11-04", "Stock9.2010-11-04", "Stock10.2010-11-04"
), class = "data.frame")

我现在通过从开头提到的论文中提取的 R 代码运行上述数据。请注意，testdfst 数据框中的截距变量不包括在计算中。

require(lpSolve)

dea_results <- list()

namesDMU <- testdfst[1]
inputs <- testdfst[c(4,5,6,7,8)]
outputs <- testdfst[9]

N <- dim(testdfst)[1] # number of DMU
s <- dim(inputs)[2] # number of inputs
m <- dim(outputs)[2] # number of outputs

f.rhs <- c(rep(0,N),1) # RHS constraints
f.dir <- c(rep("<=",N),"=") # directions of the constraints
aux <- cbind(-1*inputs,outputs) # matrix of constraint coefficients in (6)
for (i in 1:N) {
  f.obj <- c(0*rep(1,s),outputs[i,]) # objective function coefficients
  f.con <- rbind(aux, c(as.matrix(inputs[i,]), rep(0, m))) # add LHS of bTz=1
  results <-lp("max",f.obj,f.con,f.dir,f.rhs,scale=1,compute.sens=TRUE) # solve LPP
  multipliers <- results$solution # input and output weights
  efficiency <- results$objval # efficiency score
  duals <- results$duals # shadow prices
  if (i==1) {
    weights <- multipliers
    effcrs <- efficiency
    lambdas <- duals [seq(1,N)]
  } else {
    weights <- rbind(weights,multipliers)
    effcrs <- rbind(effcrs , efficiency)
    lambdas <- rbind(lambdas,duals[seq(1,N)])
  }
}

report <- as.data.frame(cbind(effcrs,weights))
colnames(report) <- c('efficiency',names(inputs),
                      names(outputs)) # header
rownames(report) <- namesDMU[,1]

代码说明：

可能与我的问题有关： 注意 z >= 0 部分？如果我错了，请纠正我，但 lpSolve 必须默认采用这个，因为这个限制不直接包含在代码中？如果我能正确识别这部分，我可以对其进行调整以解决我的问题。

问题：

这是最终结果：

 > report

         efficiency rmrf      smb        hml       rmw       cma   Returns
Stock1   0.5674100    0 0.000000  0.1769187 0.0000000 1.4634319 0.9435640
Stock2   1.0000000    0 0.000000  0.0000000 0.7713486 1.4284803 1.2405844
Stock3   1.0000000    0 1.071061  0.0000000 0.0000000 7.4588210 5.3333195
Stock4   1.0000000    0 1.930968  0.0000000 0.7427269 0.4510419 1.6584521
Stock5   0.5218197    0 0.000000  0.2080023 0.0000000 1.7205486 1.1093429
Stock6   0.5498426    0 0.000000  9.3299443 0.0000000 0.3473408 0.8384645
Stock7   1.0000000    0 3.260381  0.0000000 0.0000000 1.0000000 2.4139536
Stock8   0.0000000    0 0.000000  0.0000000 2.2130199 0.0000000 0.0000000
Stock9   0.2756548    0 0.000000 11.5264372 0.0000000 0.4291132 1.0358592
Stock10  1.0000000    0 0.000000  0.1875005 0.0000000 1.5509620 1.0000000

在所有情况下，RMRF 输入的权重均为 0。为了使这个效率度量与我应用它的上下文相关，每次计算都必须考虑所有 5 个输入。

我要解决的问题是进一步限制权重。我希望将所有 5 个输入合并到问题中，所以我想像这样限制“v”权重：

0.2 <= V(n) / V(n+1) <=5

与此相同

1/5 <= V(n) / V(n+1) <= 5/1

通过这种方式，所有权重都应限制为最大为任何其他权重的 5 倍，并且最小为 5 倍。

我一直在检查模型中的每个变量，并尝试了几个不同的数据集。我还交换了数据帧中 rmrf 和 smb 的位置，试图查看是代码有问题还是数据有问题。

        efficiency       smb      rmrf       hml      rmw        cma  Returns
Stock1   1.0000000 0.2540156 0.0000000 1.1812905 1.043781 0.03466956 1.000000
Stock2   1.0000000 1.2150960 0.0000000 0.5443910 2.854127 0.00000000 2.422061
Stock3   1.0000000 0.0000000 0.0000000 1.0000000 0.000000 1.48815164 1.104670
Stock4   1.0000000 1.2348830 0.0000000 0.0000000 4.088438 0.89216307 3.674087
Stock5   0.8151261 0.0000000 0.0000000 0.5921236 1.173539 0.61261329 1.231220
Stock6   0.1491688 0.0000000 0.0000000 2.8715984 1.262761 0.00000000 1.148347
Stock7   1.0000000 2.0599619 0.0000000 0.0000000 0.000000 1.03785902 1.520484
Stock8   1.0000000 0.8346388 0.1206426 0.0000000 1.862840 0.00000000 1.535768
Stock9   0.0000000 0.0000000 0.0000000 0.0000000 1.167498 0.00000000 0.000000
Stock10  0.8065724 0.0000000 0.6973118 2.9529776 1.318778 0.00000000 1.357774

在尝试了不同的数据集之后，似乎存在包含 RMRF 的例外情况（如上面的那个）。我开始相信我的问题的答案可能非常简单。该程序发现最好不将我的 RMRF 输入包含在计算中。解决这个问题的方法是对权重进行额外限制。

这是一个具有非标准化输入的解决方案示例（数据中的数字分布完全不同）

           efficiency smb      rmrf hml rmw cma  Returns
Stock1           1   0 1.0775390   0   0   0 346.1943
Stock2           0   0 1.3576882   0   0   0   0.0000
Stock3           0   0 0.7443477   0   0   0   0.0000
Stock4           0   0 0.6879011   0   0   0   0.0000
Stock5           0   0 1.2115507   0   0   0   0.0000
Stock6           0   0 1.1305046   0   0   0   0.0000
Stock7           0   0 1.2472639   0   0   0   0.0000
Stock8           0   0 1.4552312   0   0   0   0.0000
Stock9           0   0 1.1937843   0   0   0   0.0000
Stock10          0   0 1.3569824   0   0   0   0.0000

【问题讨论】：

标签： r linear-programming restriction lpsolve

【解决方案1】：

首先解决您的第二个问题，通过乘以V(n+1)，约束0.2 <= V(n) / V(n+1) 可以转换为：

0.2 * V(n+1) <= V(n)

由于这是一个线性约束，可以直接添加到模型中。类似地，不等式的另一边可以建模为线性约束

V(n) <= 5 * V(n+1)

关于您的第一个问题，因为除了链接到另一篇描述变量的论文的论文之外，您还没有真正描述算法，因此您至少没有提供足够的上下文让我能够回答。如果您希望 rmrf 系数取正值，您可以添加一个约束，要求 rmrf 系数的总和超过您选择的某个正值；这将强制至少一个系数为正。

【讨论】：

嗨 josilber，我编辑了问题并添加了有关问题实际试图解决的更多信息。当我第一次写这个问题时，我没有意识到我错过了主要想法。谢谢你提出来。至于我想添加的额外限制，将其分成两部分肯定会有所帮助。
嗨 josil，通过反复试验，我发现了一些我的 RMRF 输入包含在计算中的情况。所以我开始相信问题可能不是代码，但对于计算来说根本不包括这个变量是最佳的。你认为那可能是真的吗？我再次更新了问题以包含此信息。谢谢！