带组的R二进制整数优化

问题描述

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

我有三个等长的向量（52），每个向量都将与目标函数中的二进制参数向量相乘。

library(Rsolnp)
a <- c(251,179,215,251,63,45,54,47,34,40,141,101,121,94,67,81,157,108,133,126,85,106,110,74,92,52,63)
b <- c(179,0)
c <- c(179,118,30,22,44,71,56,49,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),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求解器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)是2d集，具有组和项之间的映射。

对于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代码。但这可能会给您一些思考。

,

我设置了一个R笔记本作为混合整数线性程序来解决（或尝试解决）该问题，使用CPLEX作为MIP求解器，并使用Rcplex包作为其接口。结果并不明显。经过五分钟的磨合，CPLEX的解决方案略逊于Erwin的解决方案（-886.8748 v。他的-889.346），其差距超过146％（考虑到Erwin的结果，该差距大多只是非常缓慢地收敛的上限）。我很高兴与大家分享显示线性化的笔记本，但是要使用它，您需要安装CPLEX。

更新：我还有第二本笔记本，使用GA遗传算法软件包，可以在不到五秒钟的时间内始终接近Erwin的解决方案（并偶尔点击它）。结果是随机的，因此重新运行可能会更好（或更坏），并且没有最优的证据。

,

在CPLEX中，您可以尝试保罗编写的数学编程，但也可以使用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,47,34,40,141,101,121,94,67,81,157,108,133,126,85,106,110,74,92,52,63];
int b[r]=[179,0];
int c[r]=[179,118,30,22,44,71,56,49,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秒内

NB：您可以从R调用OPL CPLEX或依靠任何其他CPLEX API

在python中，您可以编写相同的文字

from docplex.cp.model import CpoModel

n=52

r=range(0,n)

a =[251,63]
b =[179,0]
c =[179,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