混合整数非线性规划_(FOCUS)对于混合整数规划问题的 Capacity scaling method

参考文献

benders​www.arvinzyy.cn

这篇分享用的模型是上面这个链接中的。

---

Capacity scaling method

这里和大家分享的 Capacity scaling method 是根据带有固定费用的优化问题的特殊构造而产生的算法。

这里我们来考虑一个一般的带有固定费用的混合整数规划 (MIP) 问题。

在

中

是变动费用，

是固定费用。

(Capacity scaling method)

1. 设定平滑化参数
   ，把变化的容量参数
   的初期值设定为
   。
2. 反复执行以下操作，直到满足适当的收束条件后停止。(比如，线性规划松弛问题的解
   距离 0 或者1 十分接近)

• 把0-1变量
  松弛为
  ，约束条件变为
  ，解线性规划松弛问题。
• 得到的线性松弛规划问题的解
  来更新容量

Capacity scaling method 求解速度很快，并且得到是局部最优解。如果想跳出局部最优解，需要另外结合其他的方法。

Exercise

我们给出一个经典的带固定费用的运输问题，模型如下

是一个很大的实数，可看作为

的上界，这里我们可设

。

在

中，我们找出带有容量限制的约束条件，即为

首先，我们初始化

(也可尝试其他数值)，

，然后替换上边的约束条件为

其次，我们把变量

松弛为

。

这样我们得到新的模型，如下

下面我们调用 Gurobi 来实现这一算法

import numpy as np
from gurobipy import *

example = 1

if example == 1:
    # example 1, optimal 350
    num_s = 4
    num_d = 3
    supply = [10, 30, 40, 20]
    demand = [20, 50, 30]
    c = np.array([
        [2, 3, 4],
        [3, 2, 1],
        [1, 4, 3],
        [4, 5, 2]
    ])
    f = np.array(
        [[10, 30, 20] for _ in range(4)]
    )

elif example == 2:
    # example 2, optimal 4541
    num_s = 8
    num_d = 7
    supply = [20, 20, 20, 18, 18, 17, 17, 10]
    demand = [20, 19, 19, 18, 17, 16, 16]
    c = np.array([
        [31, 27, 28, 10, 7, 26, 30],
        [15, 19, 17, 7, 22, 17, 16],
        [21, 17, 19, 29, 27, 22, 13],
        [9, 19, 7, 15, 20, 17, 22],
        [19, 7, 18, 10, 12, 27, 23],
        [8, 16, 10, 10, 11, 13, 15],
        [14, 32, 22, 10, 22, 15, 19],
        [30, 27, 24, 26, 25, 15, 19]
    ]).reshape(num_s, num_d)
    f = np.array([
        [649, 685, 538, 791, 613, 205, 467],
        [798, 211, 701, 506, 431, 907, 945],
        [687, 261, 444, 264, 443, 946, 372],
        [335, 385, 967, 263, 423, 592, 939],
        [819, 340, 233, 889, 211, 854, 823],
        [307, 620, 845, 919, 223, 854, 823],
        [560, 959, 782, 417, 358, 589, 383],
        [375, 791, 720, 416, 251, 887, 235]
    ])

else:
    print("Example not found")

M = np.ones((num_s, num_d))
for i in range(num_s):
    for j in range(num_d):
        M[i, j] = min(supply[i], demand[j])

params = dict()
params["num_s"] = num_s
params["num_d"] = num_d
params["supply"] = supply
params["demand"] = demand
params["c"] = c
params["f"] = f
params["M"] = M

mip = Model()
x = mip.addVars([i for i in range(num_s)], [j for j in range(num_d)], vtype=GRB.CONTINUOUS, name="x")
y = mip.addVars([i for i in range(num_s)], [j for j in range(num_d)], vtype=GRB.CONTINUOUS,lb=0.0,ub=1.0, name="y")
mip.setObjective(
    quicksum(
            x[i,j] * c[i,j] + y[i,j] * f[i,j]
            for j in range(num_d)
        for i in range(num_s)
    ), GRB.MINIMIZE
)
for i in range(num_s) :
    mip.addConstr((
        quicksum(
            x[i,j] for j in range(num_d)
        ) <= supply[i]
        
    ))
cons=[]
for j in range(num_d):
     cons.append(mip.addConstr((
         quicksum(
             x[i,j] for i in range(num_s)
         ) >= demand[j]
   
     )))

obj_original=quicksum(x[i,j] * c[i,j] + y[i,j] * f[i,j] for j in range(num_d)for i in range(num_s))

mip.setObjective(obj_original,GRB.MINIMIZE)
mip.update()
N_0,N_1={},{}
iterlimit=100
lambda_value=0.10
UB=100000
cons1=[]
M_bar= np.ones((num_s, num_d))
M_bar=M.copy()
for iter in range(iterlimit):
    for i in range(num_s):
        crelax1=[]
        for j in range(num_d):
            crelax1.append(mip.addConstr((
                x[i,j] <= M_bar[i,j] * y[i,j]  
            )))
        cons1.append(crelax1)
    print("================LP model==================")
    mip.update()
    mip.optimize()
    for i in range(num_s):
        for j in range(num_d):
            M_bar[i, j] = M_bar[i, j]*y[i,j].x*lambda_value+(1-lambda_value)*M_bar[i, j]
            y[i,j].ub=M[i, j]/M_bar[i, j]
            print(y[i,j].x)

    for i in range(num_s):
        for j in range(num_d):
            mip.remove(cons1[i][j])
    cons1=[]
for i in range(num_s):
    for j in range(num_d):
        if y[i,j].x>=0.99:
            N_1[i,j]=1
        else:
            N_0[i,j]=0
for i in range(num_s):
    crelax=[]
    for j in range(num_d):
        if (i,j) in N_0:
            mip.addConstr((
                y[i,j]==0  
            ))
            flag=1
        if (i,j) in N_1:
            mip.addConstr((
                y[i,j]==1  
            ))
for i in range(num_s):
    for j in range(num_d):
        mip.addConstr((
            x[i,j] <= M[i,j] * y[i,j]  
        ))
print("================integer model==================")
mip.update()
mip.optimize()
if mip.status==GRB.INFEASIBLE:
    print("****model is infeasible****")
if mip.status==GRB.OPTIMAL:
    print("****model is optimal****")
    if mip.Objval<=UB:
        UB=mip.Objval
print("The upper bound is:")
print(UB)

这里我们来看example 1，线性规划松弛问题的目标函数值以及变量

的值如下

然后我们经过100次迭代后，线性规划松弛问题的目标函数值以及变量

值如下

我们可以看出

的值充分接近0 或者 1，然后我们将其固定为0 或者 1，解

，结果如下

我们得到了一个近似值为360.0，而真正的最优值为350.0。由此可知这个算例得到的近似值不算差。大家也可尝试验证以下其他的算例。

大家如果手中有自己的模型之类的一定要尝试一下上述的算法！

如果有什么疑问或者建议，可以给我发邮件。➡ zhaoyou728 at outlook.com