带组的 R 二进制整数优化答案

【问题标题】：R Binary Integer Optimization with Groups带组的 R 二进制整数优化
【发布时间】：2020-12-14 11:46:54
【问题描述】：

我试图让 Rsolnp 将我的参数限制为二进制整数或几乎相同的小数（例如，.999 足够接近 1）。

我有三个等长 (52) 的向量，每个向量都将乘以我的目标函数中的二进制参数向量。

library(Rsolnp)
a <- c(251, 179, 215, 251, 63, 45, 54, 63, 47, 34, 40, 47, 141, 101, 121, 141, 47, 34, 40, 47, 94, 67, 81, 94, 47, 34, 40, 47, 157, 108, 133, 157, 126, 85, 106, 126, 126, 85, 106, 126, 110, 74, 92, 110, 110, 74, 92, 110, 63, 40, 52, 63)
b <- c(179, 251, 215, 0, 45, 63, 54, 0, 34, 47, 40, 0, 101, 141, 121, 0, 34, 47, 40, 0, 67, 94, 81, 0, 34, 47, 40, 0, 108, 157, 133, 0, 85, 126, 106, 0, 85, 126, 106, 0, 74, 110, 92, 0, 74, 110, 92, 0, 40, 63, 52, 0)
c <- c(179, 179, 118, 179, 45, 45, 30, 45, 34, 34, 22, 34, 101, 101, 67, 101, 34, 34, 22, 34, 67, 67, 44, 67, 34, 34, 22, 34, 108, 108, 71, 108, 85, 85, 56, 85, 85, 85, 56, 85, 74, 74, 49, 74, 74, 74, 49, 74, 40, 40, 27, 40)

x 是我的参数向量，下面是我的目标函数。

objective_function = function(x){
  -(1166 *  sum(x[1:52] * a) / 2000) *
    (((sum(x[1:52] * b)) / 2100) + .05) *
    (((sum(x[1:52] * c))/1500) + 1.5)
}

我基本上希望每组 4 中的 1 个参数等于 1，其余为 0，我不确定如何为此创建正确的约束，但我相信我需要将这些总和约束与另一种类型的约束也是如此。这是我的限制：

eqn1=function(x){
  z1=sum(x[1:4])
  z2=sum(x[5:8])
  z3=sum(x[9:12])
  z4=sum(x[13:16])
  z5=sum(x[17:20])
  z6=sum(x[21:24])
  z7=sum(x[25:28])
  z8=sum(x[29:32])
  z9=sum(x[33:36])
  z10=sum(x[37:40])
  z11=sum(x[41:44])
  z12=sum(x[45:48])
  z13=sum(x[49:52])
  return(c(z1,z2,z3,z4,z5,z6,z7,z8,z9,z10,z11,z12,z13))
}

最后，这是我的函数调用：

opti<-solnp(pars=rep(1,52), fun = objective_function, eqfun = eqn1, eqB = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), LB=rep(0,52))

调用 opti$pars 会返回我的解向量：

 [1] 7.199319e-01 2.800680e-01 6.015388e-08 4.886578e-10 5.540961e-01 4.459036e-01 2.906853e-07 4.635970e-08 5.389325e-01
[10] 4.610672e-01 2.979195e-07 3.651954e-08 6.228346e-01 3.771652e-01 1.980380e-07 3.348488e-09 5.389318e-01 4.610679e-01
[19] 2.979195e-07 3.651954e-08 5.820231e-01 4.179766e-01 2.099869e-07 2.624076e-08 5.389317e-01 4.610680e-01 2.979195e-07
[28] 3.651954e-08 6.499878e-01 3.500120e-01 1.959133e-07 1.059012e-08 6.249098e-01 3.750900e-01 2.588037e-07 1.752927e-08
[37] 6.249106e-01 3.750892e-01 2.588037e-07 1.752927e-08 6.095743e-01 3.904254e-01 2.741968e-07 2.233806e-08 6.095743e-01
[46] 3.904254e-01 2.741968e-07 2.233806e-08 5.679608e-01 4.320385e-01 6.821224e-07 3.997882e-08

正如我们所看到的，权重在每组 4 个中的多个变量之间分配，而不是被强制为 1，其余为 0。

如果这个包无法做到这一点，有人可以告诉我如何将我的目标函数转换为与其他优化包一起使用吗？据我所见，它们需要将目标函数转换为系数向量。任何帮助表示赞赏。谢谢！

【问题讨论】：

看起来像一个非线性目标函数、线性约束和二元变量的问题。这称为 MINLP 问题（混合整数非线性规划）。
也许你想看看 gurobi。
请提供最低限度！可重现的例子。
@Christoph 您应该能够运行提供的代码。理想情况下，我的解决方案向量看起来像 0,1,0,0,1,0,0,0...，每四个参数中就有一个等于 1。

标签： r optimization binary linear-programming mixed-integer-programming

【解决方案1】：

我尝试了几个求解器。使用 MINLP 求解器 Couenne 和 Baron，我们可以直接解决这个问题。使用 Gurobi，我们需要将目标分解为两个二次部分。所有这些求解器都给出：

----    119 VARIABLE x.L  

i1  1.000,    i5  1.000,    i9  1.000,    i14 1.000,    i17 1.000,    i21 1.000,    i25 1.000,    i29 1.000
i34 1.000,    i38 1.000,    i41 1.000,    i46 1.000,    i49 1.000


----    119 VARIABLE z.L                   =     -889.346  obj

此处不打印零。

我使用了 GAMS（商业），但如果你想使用免费工具，你可以使用 Pyomo(Python) + Couenne。我不确定 R 的 MINLP 求解器，但可以从 R 中使用 Gurobi。

请注意，组约束很简单：

groups(g)..  sum(group(g,i),x(i)) =e= 1;

其中g 是组，group(g,i) 是具有组和项目之间映射的二维集。

对于 Gurobi，您需要执行以下操作（在伪代码中）：

z1 = 1166 *  sum(i,x(i)*a(i)) / 2000        (linear)
z2 = ((sum(i, x(i)*b(i))) / 2100) + .05     (linear) 
z3 = ((sum(i, x(i)*c(i)))/1500) + 1.5       (linear) 
z23 = z2*z3                                 (non-convex quadratic)
obj = -z1*z23                               (non-convex quadratic)

并告诉 Gurobi 使用非凸 MIQCP 求解器。

抱歉，这里没有 R 代码。但它可能会让你有所思考。

【讨论】：

【解决方案2】：

在 CPLEX 中，您可以尝试 Paul 所写的数学编程，但您也可以使用 Constraint Programming。

在 OPL（CPLEX 建模语言）中

using CP;

execute
{
  cp.param.timelimit=5; // time limit 5 seconds
  
}

int n=52;
range r=1..n;

int a[r]=[251, 179, 215, 251, 63, 45, 54, 63, 47, 34, 40, 47, 141, 101, 121, 141, 47, 34, 40, 47, 94, 67, 
81, 94, 47, 34, 40, 47, 157, 108, 133, 157, 126, 85, 106, 126, 126, 
85, 106, 126, 110, 74, 92, 110, 110, 74, 92, 110, 63, 40, 52, 63];
int b[r]=[179, 251, 215, 0, 45, 63, 54, 0, 34, 47, 40, 0, 101, 141, 121, 0, 
34, 47, 40, 0, 67, 94, 81, 0, 34, 47, 40, 0, 108, 157, 133, 0, 85, 126, 106, 0, 85, 
126, 106, 0, 74, 110, 92, 0, 74, 110, 92, 0, 40, 63, 52, 0];
int c[r]=[179, 179, 118, 179, 45, 45, 30, 45, 34, 34, 22, 34, 101, 101, 67, 101, 
34, 34, 22, 34, 67, 67, 44, 67, 34, 34, 22, 34, 108, 108, 71, 108, 85, 85, 56, 85,
 85, 85, 56, 85, 74, 74, 49, 74, 74, 74, 49, 74, 40, 40, 27, 40];

// decision variable 
 dvar boolean x[r];
 
 // objective
 dexpr float obj=
 
  -(1166 *  sum(i in r) (x[i]*a[i]) / 2000) *
    (((sum(i in r) (x[i]* b[i])) / 2100) + .05) *
    (((sum(i in r) (x[i]*c[i]))/1500) + 1.5);
    
minimize obj;

subject to
{
  // one and only one out of 4 is true
  forall(i in 1..n div 4) count(all(j in 1+(i-1)*4..4+(i-1)*4)x[j],1)==1;
  
}

给予

// solution with objective -889.3463
x = [1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1
         0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0
         0 0];

5 秒内

注意：您可以从 R 调用 OPL CPLEX 或依赖任何其他 CPLEX API

在python中你也可以这样写

from docplex.cp.model import CpoModel

n=52

r=range(0,n)

a =[251, 179, 215, 251, 63, 45, 54, 63, 47, 34, 40, 47, 141, 101, 121, 141, 47, 34, 40, 47, 94, 67, 81, 94, 47, 34, 40, 47, 157, 108, 133, 157, 126, 85, 106, 126, 126, 85, 106, 126, 110, 74, 92, 110, 110, 74, 92, 110, 63, 40, 52, 63]
b =[179, 251, 215, 0, 45, 63, 54, 0, 34, 47, 40, 0, 101, 141, 121, 0, 34, 47, 40, 0, 67, 94, 81, 0, 34, 47, 40, 0, 108, 157, 133, 0, 85, 126, 106, 0, 85, 126, 106, 0, 74, 110, 92, 0, 74, 110, 92, 0, 40, 63, 52, 0]
c =[179, 179, 118, 179, 45, 45, 30, 45, 34, 34, 22, 34, 101, 101, 67, 101, 34, 34, 22, 34, 67, 67, 44, 67, 34, 34, 22, 34, 108, 108, 71, 108, 85, 85, 56, 85, 85, 85, 56, 85, 74, 74, 49, 74, 74, 74, 49, 74, 40, 40, 27, 40]

mdl = CpoModel(name='x')

#decision variables
mdl.x = {i: mdl.integer_var(0,n,name="x"+str(i+1)) for i in r}

mdl.minimize(-1166 *  sum(mdl.x[i]*a[i] / 2000 for i in r) \
             *((sum(mdl.x[i]* b[i] / 2100  for i in r) +0.05)) \
             *((sum(mdl.x[i]*c[i]/1500  for i in r) +1.5)) )

for i in range(0,n // 4):
    mdl.add(1==sum( mdl.x[j] for j in range(i*4+0,i*4+4)))

msol=mdl.solve(TimeLimit=5)

# Dislay solution
for i in r:
    if (msol[mdl.x[i]]==1):
        print(i+1," ")

这给了

! Best objective         : -889.3464 
 
1  
5  
9  
13  
17  
22  
25  
30  
34  
38  
41  
45  
49

【讨论】：

知道了！谢谢亚历克斯
@Alex Fleischer：很好的答案！你能看看我的问题吗？ stackoverflow.com/questions/69319495/… 谢谢！

【解决方案3】：

我设置了一个 R 笔记本来解决（或尝试解决）作为混合整数线性程序的问题，使用 CPLEX 作为 MIP 求解器，并使用 Rcplex 包作为它的接口。结果并不引人注目。经过 5 分钟的打磨后，CPLEX 得到的解比 Erwin 得到的解（-886.8748 与他的 -889.346）略逊一筹，差距超过 146%（考虑到 Erwin 的结果，这主要是收敛非常缓慢的上限）。我很高兴分享显示线性化的笔记本，但要使用它，您需要安装 CPLEX。

更新：我有第二个笔记本，使用 GA 遗传算法包，它在五秒内始终接近 Erwin 的解决方案（并且偶尔会达到它）。结果是随机的，因此重新运行可能会做得更好（或更差），并且没有最优性证明。

【讨论】：

我目前正在用 python 做这件事，但如果你能和我分享这两个笔记本，那就太棒了@prubin
我用我使用的 cplex python 代码更新了我的答案
第一个笔记本（使用 CPLEX）：rubin.msu.domains/blog/groupselect.nb.html.
第二个笔记本（使用 GA）：rubin.msu.domains/blog/groupselect2.nb.html.