带组的 R 二进制整数优化
R Binary Integer Optimization with Groups
我试图让 Rsolnp 将我的参数限制为二进制整数或几乎相同的小数（例如，.999 足够接近 1）。
我有三个等长的向量 (52)，每个向量都将乘以我的 objective 函数中的二进制参数向量。
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 函数。
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 returns 我的解决方案向量：
[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。
如果这个包无法做到这一点，有人可以告诉我如何将我的 objective 函数转换为与其他优化包一起工作吗？据我所知，它们需要将 objective 函数转换为系数向量。任何帮助表示赞赏。谢谢！
我尝试了几个求解器。使用 MINLP 求解器 Couenne 和 Baron，我们可以直接解决这个问题。使用 Gurobi，我们需要将 objective 分解为两个二次部分。所有这些求解器给出：
---- 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 代码。但它可能会给你一些思考。
我设置了一个 R 笔记本来解决（或尝试解决）作为混合整数线性程序的问题，使用 CPLEX 作为 MIP 求解器，并使用 Rcplex 包作为它的接口。结果并不引人注目。经过五分钟的磨合，CPLEX 的解决方案比 Erwin 得到的解决方案稍逊一筹（-886.8748 对他的 -889.346），差距超过 146%（考虑到 Erwin 的结果，这主要是上限收敛非常缓慢）。我很高兴分享笔记本，它显示了线性化，但要使用它，您需要安装 CPLEX。
更新：我有第二个笔记本，使用 GA 遗传算法包，它在五秒内始终接近 Erwin 的解决方案（偶尔会命中）。结果是随机的，因此重新运行可能会做得更好（或更差），并且没有最优性证明。
在 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
我试图让 Rsolnp 将我的参数限制为二进制整数或几乎相同的小数（例如，.999 足够接近 1）。
我有三个等长的向量 (52)，每个向量都将乘以我的 objective 函数中的二进制参数向量。
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 函数。
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 returns 我的解决方案向量：
[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。
如果这个包无法做到这一点，有人可以告诉我如何将我的 objective 函数转换为与其他优化包一起工作吗？据我所知，它们需要将 objective 函数转换为系数向量。任何帮助表示赞赏。谢谢！
我尝试了几个求解器。使用 MINLP 求解器 Couenne 和 Baron，我们可以直接解决这个问题。使用 Gurobi，我们需要将 objective 分解为两个二次部分。所有这些求解器给出：
---- 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 代码。但它可能会给你一些思考。
我设置了一个 R 笔记本来解决（或尝试解决）作为混合整数线性程序的问题，使用 CPLEX 作为 MIP 求解器，并使用 Rcplex 包作为它的接口。结果并不引人注目。经过五分钟的磨合，CPLEX 的解决方案比 Erwin 得到的解决方案稍逊一筹（-886.8748 对他的 -889.346），差距超过 146%（考虑到 Erwin 的结果，这主要是上限收敛非常缓慢）。我很高兴分享笔记本，它显示了线性化，但要使用它，您需要安装 CPLEX。
更新：我有第二个笔记本，使用 GA 遗传算法包，它在五秒内始终接近 Erwin 的解决方案（偶尔会命中）。结果是随机的，因此重新运行可能会做得更好（或更差），并且没有最优性证明。
在 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