Logic-based Benders分解学习笔记(一)

发表于 更新于 阅读次数： 本文字数： 15k 阅读时长 ≈ 27 分钟
Logic-based Benders分解学习过程记录

[TOC]

学习目标

1.掌握Logic-based Benders分解的原理、学会如何应用

2.理解Logic-based Benders分解与Benders分解的区别

3.将Logic-based Benders分解应用到毕业论文的第五章

理论介绍

核心思想： 从错误中学习 (Learn from mistake)。

Logic-based Benders分解与Benders分解的区别: 逻辑 Benders分解（Logic-based Benders decomposition）是一种基于经典 Benders 分解的扩展方法，通过逻辑推理代替传统的线性对偶来生成 Benders 切割，从而适用于更广泛的问题类型。这种方法的关键在于定义一种“推理对偶(inference duality)”而非线性对偶，通过逻辑推理在约束下生成最强的下界。这种方法不仅适用于线性优化，还可处理可满足性问题（satisfiability）、0-1 编程和机器调度问题等。

推理对偶

推理对偶是逻辑benders分解的核心。为了学习推理对偶，需要引入语义蕴含(semantic implication)这个概念。

语义蕴含

这里还没太明白，下次再写

logic-based benders求解

针对参考文献中这篇文章的问题。

具体的benders cut可以设置为：$C_{\max} \geq C_{\max}^{hi’^*} - \sum_{j \in N_{i’}^h} (1 - x_{i’j}) \theta_{hi’j}$

为什么这么设置？因为其具有两大属性。

首先，这个cut必须从主问题的可行空间中移除当前的解；其次，这个cut不会将全局最优解从解空间中移除。

并行机调度问题

2013年的state-of-the-art model

image-20241106161912566

为了适合LBBD的算法框架，我们建立的数学模型如下：

image-20241106162000709

先直接用求解器求解，不使用逻辑Benders分解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
"""
Description: 求解并行机调度问题
Author: TUUG
Date: 2024/11/04 16:48
Version: V5.0
"""

import gurobipy as gp
from gurobipy import GRB
import numpy as np
import pandas as pd

processing_time = pd.read_csv('工件加工时间.csv',index_col=0).T
relation = pd.read_csv('工件优先关系.csv',index_col=0)  # 按列截取再转化就是前序加工集合，按行截取再转化就是后续加工集合
set_up = pd.read_csv('工件切换时间.csv',index_col=0)
jobs = ['job1','job2','job3','job4','job5']
serus = ['seru1','seru2','seru3']
num_serus = len(serus)
num_jobs = len(jobs)
batch_info = pd.read_csv('批量大小.csv',index_col=0)
batch_info['size'] = np.random.normal(loc=batch_info['mean'],scale=batch_info['std']).astype(int)
batch_size = batch_info['size'].to_list()

relation = [('job1','job2')]  # 这里存储一些不能够先后加工的约束。
dummy_job = ['job0'] # 表示两个虚拟工件

def get_pre(job):
    """返回该工件最大的切换时间,benders cut需要用到"""
    return max(set_up[job])

max_pre_dict = dict(zip(jobs+dummy_job,[get_pre(job) for job in jobs+dummy_job]))

def get_theta(ind,job):
    """返回theta值,benders cut需要用到"""
    return max_pre_dict[job]+processing_time.loc[ind,job]

master_model = gp.Model('master_problem')
x = master_model.addVars(serus,jobs+dummy_job,vtype=GRB.BINARY,name='x')
y = master_model.addVars(serus,jobs+dummy_job,jobs+dummy_job,vtype=GRB.BINARY,name='y')
C_max = master_model.addVar(vtype=GRB.CONTINUOUS, name="C_max")
C = master_model.addVars(jobs+dummy_job,vtype=GRB.CONTINUOUS,name='C')
xi = master_model.addVars(serus,vtype=GRB.CONTINUOUS,name='xi')
master_model.setObjective(C_max, GRB.MINIMIZE)
M = 10086

# 约束22
for i in serus:
    master_model.addConstr(gp.quicksum(x[i,j]*processing_time.loc[i,j] for j in jobs)+xi[i] <= C_max)

# 约束23
for j in jobs:
    master_model.addConstr(gp.quicksum(x[i,j] for i in serus) == 1)

# 约束24
for i in serus:
    master_model.addConstr(x[i,'job0'] == 1)

# 约束25
for i in serus:
    master_model.addConstr(xi[i] == gp.quicksum(y[i,j,k]*set_up.loc[j,k] for j in jobs+dummy_job for k in jobs+dummy_job))

# 约束26
for i in serus:
    for k in jobs+dummy_job:
        master_model.addConstr(x[i,k] == gp.quicksum(y[i,j,k] for j in jobs+dummy_job))

# 约束27
for i in serus:
    for j in jobs+dummy_job:
        master_model.addConstr(x[i,j] == gp.quicksum(y[i,j,k] for k in jobs+dummy_job))

# 约束28
for i in serus:
    for j in jobs+dummy_job:
        for k in jobs:
            master_model.addConstr(C[k] -C[j] +M*(1-y[i,j,k]) >= set_up.loc[j,k]+processing_time.loc[i,k])

# 约束29
master_model.addConstr(C['job0'] == 0)

# 先后关系约束
for i in serus:
    for pair in relation:
        master_model.addConstr(y[i,pair[0],pair[1]] == 0)

master_model.optimize()
master_model.write('right.lp')
if master_model.status == GRB.OPTIMAL:
    assignment = np.array([[1 if x[m, j].X > 0.1 else 0 for j in jobs+dummy_job] for m in serus])
    print('----------','x的取值情况','--------------------')
    for row in serus:
        for col in jobs:
            if x[row, col].X > 0:
                print(x[row, col])
    
    print('----------','y的取值情况','----------')
    for row in serus:
        for col in jobs:
            for k in jobs:
                if y[row, col, k].X > 0:
                    print(y[row,col,k])
    
    print('----------','C的取值情况','--------------------')
    for col in jobs:
        if C[col].X > 0:
            print(C[col])
    
    print('最大完工时间为：',str(C_max.X))

求解结果如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (win64)

CPU model: AMD Ryzen 9 7945HX with Radeon Graphics, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 16 physical cores, 32 logical processors, using up to 32 threads

Optimize a model with 144 rows, 136 columns and 613 nonzeros
Model fingerprint: 0x5891c35e
Variable types: 10 continuous, 126 integer (126 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+04]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+04]
Presolve removed 32 rows and 30 columns
Presolve time: 0.00s
Presolved: 112 rows, 106 columns, 503 nonzeros
Variable types: 5 continuous, 101 integer (100 binary)
Found heuristic solution: objective 41.0000000
Found heuristic solution: objective 35.0000000
Found heuristic solution: objective 33.0000000

Root relaxation: objective 1.378378e+01, 59 iterations, 0.00 seconds (0.00 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0   13.78378    0   25   33.00000   13.78378  58.2%     -    0s
H    0     0                      25.0000000   13.78378  44.9%     -    0s
H    0     0                      18.0000000   13.78378  23.4%     -    0s
H    0     0                      17.0000000   14.20076  16.5%     -    0s
     0     0   15.19231    0   30   17.00000   15.19231  10.6%     -    0s
     0     0   15.19231    0   31   17.00000   15.19231  10.6%     -    0s
     0     0   15.78689    0   25   17.00000   15.78689  7.14%     -    0s

Cutting planes:
  Learned: 2
  Gomory: 14
  Cover: 1
  Implied bound: 3
  Clique: 6
  MIR: 2
  Zero half: 2
  RLT: 4

Explored 1 nodes (124 simplex iterations) in 0.03 seconds (0.01 work units)
Thread count was 32 (of 32 available processors)

Solution count 6: 17 18 25 ... 41

Optimal solution found (tolerance 1.00e-04)
Best objective 1.700000000000e+01, best bound 1.700000000000e+01, gap 0.0000%
---------- x的取值情况 --------------------
<gurobi.Var x[seru1,job3] (value 1.0)>
<gurobi.Var x[seru2,job1] (value 1.0)>
<gurobi.Var x[seru2,job5] (value 1.0)>
<gurobi.Var x[seru3,job2] (value 1.0)>
<gurobi.Var x[seru3,job4] (value 1.0)>
---------- y的取值情况 ----------
<gurobi.Var y[seru2,job5,job1] (value 1.0)>
<gurobi.Var y[seru3,job4,job2] (value 1.0)>
---------- C的取值情况 --------------------
<gurobi.Var C[job1] (value 14.0)>
<gurobi.Var C[job2] (value 17.0)>
<gurobi.Var C[job3] (value 16.0)>
<gurobi.Var C[job4] (value 8.0)>
<gurobi.Var C[job5] (value 10.0)>
最大完工时间为： 17.0

采用逻辑benders分解后求解结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
# %%
"""
Description: 复现Informs journal on computing上的论文。
Author: TUUG
Date: 2024/11/06 10:48
Version: V6.0
"""

import gurobipy as gp
from gurobipy import GRB
import numpy as np
import pandas as pd

processing_time = pd.read_csv('工件加工时间.csv',index_col=0).T
relation = pd.read_csv('工件优先关系.csv',index_col=0)  # 按列截取再转化就是前序加工集合，按行截取再转化就是后续加工集合
set_up = pd.read_csv('工件切换时间.csv',index_col=0)
jobs = ['job1','job2','job3','job4','job5']
serus = ['seru1','seru2','seru3']
num_serus = len(serus)
num_jobs = len(jobs)
batch_info = pd.read_csv('批量大小.csv',index_col=0)
batch_info['size'] = np.random.normal(loc=batch_info['mean'],scale=batch_info['std']).astype(int)
batch_size = batch_info['size'].to_list()

relation = [('job1','job2')]  # 这里存储一些不能够先后加工的约束。
dummy_job = ['job0'] # 表示两个虚拟工件

def get_pre(h,i,job,assign_jobs_list_record):
    ans_list = [set_up.loc[i,job] for i in assign_jobs_list_record[h][i]]
    return max(ans_list)

def get_theta(h,i,job,assign_jobs_list_record):
    return get_pre(h,i,job,assign_jobs_list_record)+processing_time.loc[serus[i],job]


# %%
def create_master_model(gen=None,C_max_current_record=None, assign_jobs_list_record=None):
    master_model = gp.Model('master problem')
    x = master_model.addVars(serus,jobs+dummy_job,vtype=GRB.BINARY,name='x')
    C_max = master_model.addVar(vtype=GRB.CONTINUOUS, name="C_max")
    C = master_model.addVars(jobs+dummy_job,vtype=GRB.CONTINUOUS,name='C')
    xi = master_model.addVars(serus,vtype=GRB.CONTINUOUS,name='xi')
    master_model.setObjective(C_max, GRB.MINIMIZE)
    master_model.setParam('OutputFlag', 0)

    # 约束22
    for i in serus:
        master_model.addConstr(gp.quicksum(x[i,j]*processing_time.loc[i,j] for j in jobs)+xi[i] <= C_max)

    # 约束23
    for j in jobs:
        master_model.addConstr(gp.quicksum(x[i,j] for i in serus) == 1)

    # 约束24
    for i in serus:
        master_model.addConstr(x[i,'job0'] == 1)
    
    # 约束29
    master_model.addConstr(C['job0'] == 0)
    
    if gen:
        print('********','增加了benders cut','*********')
                
        for h in range(gen):
            for i in range(len(serus)):
                assign_jobs = assign_jobs_list_record[h][i]
                if 'job0' in assign_jobs:
                    assign_jobs.remove('job0')
                master_model.addConstr(C_max >= C_max_current_record[h][i] - gp.quicksum((1-x[serus[i],j])*get_theta(h,i,j,assign_jobs_list_record) for j in assign_jobs),name='cut')

    master_model.optimize()
    master_model.write("master_model.lp")
    
    if master_model.status == GRB.OPTIMAL:
        print(master_model.objVal,'===主模型的解为=====')
        assignment = np.array([[1 if x[m, j].X > 0.1 else 0 for j in jobs+dummy_job] for m in serus])
        assignment = pd.DataFrame(assignment, columns=jobs+dummy_job, index=serus)
        print('----------','x的取值情况','--------------------')
        for row in serus:
            for col in jobs:
                if x[row, col].X > 0:
                    print(x[row, col])
    
        print('----------','C的取值情况','--------------------')
        for col in jobs:
            if C[col].X > 0:
                print(C[col])
        return assignment, master_model.objVal
    else:
        return None, None
    
def get_sub_job(assignment,ind):
    row = assignment.loc[ind]
    columns_with_1 = row[row == 1].index.tolist()
    return columns_with_1

# %%
def solve_subproblem(ind,assignment):
    """求解第ind个机器上的完工时间"""
    i = ind
    x = assignment
    
    print('--------------求解第' + str(i) + '个机器上的完工时间--------------------')
    # assign_jobs = [i for i,j in enumerate(list(assignment.loc[ind])) if j==1]  # 分配到这个seru的工件集合
    assign_jobs = get_sub_job(assignment,ind)
    sub_model = gp.Model('sub_problem')
    
    C = sub_model.addVars(jobs+dummy_job,vtype=GRB.CONTINUOUS,name='C')
    y = sub_model.addVars(serus,jobs+dummy_job,jobs+dummy_job,vtype=GRB.BINARY,name='y')
    xi = sub_model.addVars(serus,vtype=GRB.CONTINUOUS,name='xi')  # 该机器上的总切换时间
    C_machine = sub_model.addVars(serus,vtype=GRB.CONTINUOUS,name='C_machine')  # 该机器上的最大完工时间
    sub_model.setObjective(C_machine[i],sense=GRB.MINIMIZE)
    M = 10086
    
    # 自己的约束，子问题里添加的约束
    for j in jobs:
        sub_model.addConstr(C_machine[i] >= C[j]*x.loc[i,j])
    
    # 约束25
    sub_model.addConstr(xi[i] == gp.quicksum(y[i,j,k]*set_up.loc[j,k] for j in jobs+dummy_job for k in jobs+dummy_job))

    # 约束26
    for k in jobs+dummy_job:
        sub_model.addConstr(gp.quicksum(y[i,j,k] for j in jobs+dummy_job) == x.loc[i,k])

    # 约束27
    for j in jobs+dummy_job:
        sub_model.addConstr(gp.quicksum(y[i,j,k] for k in jobs+dummy_job) == x.loc[i,j])

    # 约束28
    for j in jobs+dummy_job:
        for k in jobs:
            sub_model.addConstr(C[k] -C[j] +M*(1-y[i,j,k]) >= set_up.loc[j,k]+processing_time.loc[i,k])

    # 自己的约束：不能违反的先后关系
    for pair in relation:
        sub_model.addConstr(y[i,pair[0],pair[1]] == 0)
    
    sub_model.setParam('OutputFlag',0)
            
    sub_model.optimize()
    if sub_model.status == GRB.OPTIMAL:
        print('----------','y的取值情况','----------')
        for row in serus:
            for col in jobs:
                for k in jobs:
                    if y[row, col, k].X > 0:
                        print(y[row,col,k])
        result = np.array([[1 if y[i,m,j].X > 0 else 0 for j in jobs] for m in jobs])
        
        return result,sub_model.objVal,assign_jobs
    else:
        print('未求得最优解')
        return None, None

# %%
assignment,obj = create_master_model()
gen = 1
result_obj_list = [int(obj*3)]*num_serus  # 确保第一次循环能够进入
C_max_list = [0]*num_serus
result_obj_list_record = []
assign_jobs_list_record = []
theta_record = []
LB_record = []
UB_record = []
UB_min = result_obj_list[0]

while max(result_obj_list)-obj >= 1e-6:
    
    LB_record.append(obj)
    UB_record.append(max(result_obj_list))

    print('===============第'+str(gen)+'次循环=====================')
    result_list = [solve_subproblem(i,assignment=assignment) for i in serus]
    result_assign_list = [i[0] for i in result_list]
    result_obj_list = [i[1] for i in result_list]
    assign_jobs_list = [i[2] for i in result_list]
    result_obj_list_record.append(result_obj_list)
    assign_jobs_list_record.append(assign_jobs_list)

    assignment,obj = create_master_model(gen=gen, C_max_current_record=result_obj_list_record,assign_jobs_list_record=assign_jobs_list_record)
    print('主问题的目标值',str(obj))
    gen += 1
    if gen == 100:
        print('迭代次数超过100,退出循环')
        break

print('Logic-based benders求解结束')
print('目标值为:',obj)

# %%
from matplotlib import pyplot as plt
import matplotlib.cm as cm

plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False

UB_record.append(obj)
LB_record.append(obj)

def update_to_historical_min(UB_record):
    # 初始化一个变量用于存储历史最小值
    min_value = UB_record[0]
    # 遍历数组并更新每个位置的值为历史最小值
    for i in range(len(UB_record)):
        # 更新当前位置的最小值
        min_value = min(min_value, UB_record[i])
        # 将当前位置的值改为历史最小值
        UB_record[i] = min_value
    return UB_record
# UB_record = update_to_historical_min(UB_record)

# 创建图表
plt.figure(figsize=(10, 6))
colors = cm.viridis(np.linspace(0, 1, 7))
plt.plot(UB_record, label="子问题", color='purple',alpha=0.5, marker='s',linestyle='--', linewidth=2)
plt.plot(LB_record, label="主问题", color='cornflowerblue',alpha=0.8, marker='x',linestyle='-.', linewidth=2)

# 添加标题和标签
# plt.title("逻辑benders分解迭代过程中主问题与子问题解的变化情况")
plt.xlabel("迭代次数")
plt.ylabel("最大完工时间")

# 显示图例和网格
plt.legend()
plt.grid(True,alpha=0.5,linestyle='--',color='gray')
# 显示图表
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
===============第17次循环=====================
--------------求解第seru1个机器上的完工时间--------------------
---------- y的取值情况 ----------
--------------求解第seru2个机器上的完工时间--------------------
---------- y的取值情况 ----------
<gurobi.Var y[seru2,job5,job1] (value 1.0)>
--------------求解第seru3个机器上的完工时间--------------------
---------- y的取值情况 ----------
<gurobi.Var y[seru3,job4,job2] (value 1.0)>
******** 增加了benders cut *********
17.0 ===主模型的解为=====
---------- x的取值情况 --------------------
<gurobi.Var x[seru1,job3] (value 1.0)>
<gurobi.Var x[seru2,job1] (value 1.0)>
<gurobi.Var x[seru2,job5] (value 1.0)>
<gurobi.Var x[seru3,job2] (value 1.0)>
<gurobi.Var x[seru3,job4] (value 1.0)>
---------- C的取值情况 --------------------
主问题的目标值 17.0
Logic-based benders求解结束
目标值为: 17.0

求解过程图如下

image-20241106152925359

参考文献

Tran T T, Araujo A, Beck J C. Decomposition methods for the parallel machine scheduling problem with setups[J]. INFORMS Journal on Computing, 2016, 28(1): 83-95.