个人数学建模算法库之整数线性规划模型_混合整数线性规划csdn

模型描述

整数线性规划模型是线性规划模型的子集，它要求决策变量部分或全部为整数。
$$\begin{gathered} \min z=c^{T}x \\ s.t.\{\begin{array}{l}Ax\le b\\ Aeq\cdot x=beq\\ lb\le x\le ub\end{array} \\ c、x、b、lb、ub为列向量 \\ A为矩阵 \\ 其中,x部分或全部分量为整数 \end{gathered}$$
整数线性规划可分为

• 纯整数规划:决策变量全为整数
• 全整数规划:所有参数中的元素全为整数
• 混合整数规划:决策变量既有整数，也有非整数
• 0-1整数规划:决策变量取值为0或1

算法及代码

以下面的问题为例:
max ⁡ z = x 1 + x 2 { 14 x 1 + 9 x 2 ≤ 51 − 6 x 1 + 11 x 2 ≤ 1 x 1 , x 2 ≥ 0 且为整数 \max z=x_1+x_2\\

$$\begin{cases} 14x_1+9x_2\le 51\\ -6x_1+11x_2\le 1\\ x_1,x_2\ge 0 且为整数 \end{cases}$$ maxz=x1​+x2​⎩ ⎨ ⎧​14x1​+9x2​≤51−6x1​+11x2​≤1x1​,x2​≥0且为整数​
首先化问题为标准形式
$$\begin{gathered} c=-\left[\begin{array}{c}1\\ 1\end{array}\right] \\ A=\left[\begin{array}{cc}14& 9\\ -6& 11\end{array}\right],b=\left[\begin{array}{c}51\\ 1\end{array}\right] \\ lb=\left[\begin{array}{c}0\\ 0\end{array}\right] \end{gathered}$$

分支定界算法

• 不考虑整数限制，求出对应松弛线性规划问题的最优解。
• 任选一个不满足整数条件的变量$x_i$,在上下界约束中更新 $[x_i]\le x_i$ 或 $x_i\le [x_i]+1$ ，形成两个新的线性规划子问题
• 在缩小的可行域中，分别求出两个子问题的最优解
• 重复上述步骤
• 直至子问题无解或有整数最优解(称为子问题被查清)

若所有子问题无解，考虑其他算法；若子问题有解，则在可行整数解中寻找最优解，即为原整数线性规划的最优解

Matlab版

% int

function [x,fval,status] = intprog(f,A,b,intcon,Aeq,beq,lb,ub,e)
% 整数线性规划分支定界算法
% f为目标函数向量
% intcon为整数约束
% A,b为不等式约束
% Aeq,beq为整式约束
% lb,ub为上下界约束
% e为终止误差

% 控制输入参数
if nargin < 9
    e=0.00001;
    if nargin<8
        ub=[];
        if nargin < 7
            lb=[];
            if nargin < 6
                beq=[];
                if nargin < 5
                    Aeq=[];
                    if nargin < 4
                        intcon= [1:length(f)];
                    end
                end
            end
        end
    end
end

% 求解整数规划对应的松散线性规划,并判断是否有解
options=optimset('display','off');
[x0,fval0,status0]=linprog(f,A,b,Aeq,beq,lb,ub,options);

if status0<0
    disp('没有合适的整数解');
    x=x0;
    fval=fval0;
    status=status0;
    return;
else
    % 采用分支定界法
    bound=inf;
    [x,fval,status]=branchbound(f,A,b,intcon,x0,fval0,bound,Aeq,beq,lb,ub,e);
end

end
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

% branchbound.m
function [newx,newfval,status,newbound] = branchbound(f,A,b,intcon,x,fval,bound,Aeq,beq,lb,ub,e)
% 分支定界法
% f,A,b,Aeq,Beq,lb,ub与线性规划模型相同
% bound为当前满足整数约束的最优值的上界
% intcon为整数限制
% x为初始解,fval为初始值

options=optimset('display','off');
[x0,fval0,status0]=linprog(f,A,b,Aeq,beq,lb,ub,[],options);

% 递归退出条件
% 无解或比上界大返回原解
if status0<=0 || fval0>=bound
    newx=x;
    newfval=fval;
    newbound=bound;
    status=status0;
    return;
end

% 是否为整数解,如果是则直接返回
intindex=find( abs(x0(intcon)-round(x0(intcon)) )>e, 1);
if isempty(intindex)
    newx=x0;
    newx(intcon)=round(x0(intcon));
    newfval=fval0;
    newbound=fval0;
    status=1;
    return
end

% 找到第一个不满足整数约束的变量
n=intcon(intindex(1));
addA=zeros(1,length(f));
addA(n)=1;

% 添加约束，构造第一个分支并求分支的解
A=[A;addA];
b=[b,floor(x(n))];
[x1,fval1,status1,bound1]=branchbound(f,A,b,intcon,x0,fval0,bound,Aeq,beq,lb,ub,e);
A(end,:)=[];
b(:,end)=[];

% 若第一个分支的解为更优解则替换，否则保持原状
status=status1;
if status1 >0 && bound1<bound
    newx=x1;
    newfval=fval1;
    bound=fval1;
    newbound=bound1;
else
    newx=x0;
    newfval=fval0;
    newbound=bound;
end

% 构造第二分支
A=[A;-addA];
b=[b,-ceil(x(n))];
[x2,fval2,status2,bound2]=branchbound(f,A,b,intcon,x0,fval0,bound,Aeq,beq,lb,ub,e);

% 若第二个分支的解为更优解则替换，否则保持原状
if status2 >0 && bound2<bound
    status=status2;
    newx=x2;
    newfval=fval2;
    newbound=bound2;
end

end
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

测试结果如下:

% IntprogTest.m

% 准备参数
f=-[1 1];
A=[14 9;-6 11];
b=[51 1];
lb=[0 0];

% 整数线性规划求解
[x,fval,status]= intprog(f,A,b,[1,2],[],[],lb);
display(x);
zmax=-fval;
display(zmax);
1
2
3
4
5
6
7
8
9
10
11
12
13

x =

     3
     1


zmax =

     4
1
2
3
4
5
6
7
8
9

Python版

# 导库
from scipy import optimize
import numpy as np


class Intprog:
    """整数线性规划类"""
    def __init__(self,
                 f,
                 A,
                 b,
                 intcon=None,
                 Aeq=None,
                 beq=None,
                 lb=None,
                 ub=None,
                 e=0.00001):
        # 初始化参数
        
        self.rows, self.cols = np.shape(A)
        if intcon.all() == None:
            intcon = np.array(list(range(1, cols + 1)))
        self.intcon = intcon
        if lb == None:
            lb = np.ones((self.cols, 1)) * -np.inf
        if ub == None:
            ub = np.ones((self.cols, 1)) * np.inf
        self.f, self.A, self.b, self.Aeq, self.beq, self.e = f, A, b, Aeq, beq, e
        self.bounds = np.array(list(zip(lb, ub)), dtype=object)
        self.x, self.fval, self.status = None, None, False

    def solve(self):
        # 启动求解过程

        # 求解整数规划对应的松散线性规划
        result = optimize.linprog(self.f, self.A, self.b, self.Aeq, self.beq,
                                  self.bounds)

        # 无解退出
        if result.success == False:
            return

        # 有解则采用分支定界法
        x = result.x
        fval = result.fun
        bound = np.inf
        x, fval, status, bound = self.__branchbound(self.A, self.b, x, fval,
                                                    bound)
        if np.abs(fval - np.round(fval)) < self.e:
            fval = np.round(fval)
        self.x, self.fval, self.status = x, fval, status
    
    def __branchbound(self, A, b, x, fval, bound):
        """分支定界法"""
        result = optimize.linprog(self.f, A, b, self.Aeq, self.beq,
                                  self.bounds)
        x0, fval0, status0 = result.x, result.fun, result.success  # 初始解

        # 递归退出条件
        # 无解或比上界大返回初始解
        if status0 == False or fval0 > bound:
            return x, fval, status0, bound

        # 判断是否为整数解
        for i in range(len(x0)):
            all_int = True
            if np.abs(x0[i] -
                      np.round(x0[i])) > self.e and i + 1 in self.intcon:
                all_int = False
                first_not_int = i
                break

        # 为整数解,直接返回
        if all_int:
            x0[self.intcon - 1] = np.round(x0[self.intcon - 1])
            return x0, fval0, status0, fval0

        newx, newfval, status, newbound = x0, fval0, status0, bound

        # 不为整数解,添加约束,构造第一个分支并求解
        addA = np.zeros((1, len(self.f)))
        addA[0, first_not_int] = 1
        A1 = np.concatenate((A, addA))
        b1 = np.concatenate((b, [np.floor(x[first_not_int])]))
        x1, fval1, status1, bound1 = self.__branchbound(
            A1, b1, x0, fval0, bound)

        # 如果第一个分支的解为更优解,则替换;否则保持
        if status1 == True and bound1 < bound:
            status = status1
            newx = x1
            newfval = fval1
            bound = fval1  # 更新新上界
            newbound = bound1

        # 构造第二分支并求解
        A2 = np.concatenate((A, -addA))
        b2 = np.concatenate((b, [-np.ceil(x[first_not_int])]))
        x2, fval2, status2, bound2 = self.__branchbound(
            A2, b2, x0, fval0, bound)

        # 如果第二个分支的解为更优解,则替换;否则保持
        if status2 == True and bound2 < bound:
            status = status2
            newx = x2
            newfval = fval2
            bound = fval2
            newbound = bound2
        return newx, newfval, status, newbound
    
    def show(self, for_min=True):
        """展示解"""
        if self.status == False:
            print("无满足条件的解")
			return
			
        if for_min:
            print(f"当x取 {self.x} 时,得最小值{self.fval}")
        else:
            print(f"当x取 {self.x} 时,得最大值{-self.fval}")
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

测试结果如下:

# 准备参数
f=-np.array([1 ,1])
A=np.array([[14 ,9],[-6 ,11]])
b=np.array([51 ,1])
lb=[0 ,0]

# 整数线性规划求解
pro=Intprog(f,A,b,np.array([1,2]),lb=[0,0])
pro.solve()
pro.show(for_min=False)
1
2
3
4
5
6
7
8
9
10

当x取 [3. 1.] 时,得最大值4.0
1

库函数

Matlab版

基于求解器

[x,fval]=intlinprog(f,intcon,A,b,Aeq,beq,lb,ub)
1

• intcon参数设置整数约束，指定整数约束变量的位置
• 其他参数同linprog

测试结果如下:

% 准备参数
f=-[1 ;1];
A=[14 9;-6 11];
b=[51 ;1];
lb=[0 ;0];

% 整数线性规划求解
[x,fval]=intlinprog(f,[1,2],A,b,[],[],[0,0],[]);
display(x);
zmax=-fval;
display(zmax);
1
2
3
4
5
6
7
8
9
10
11

LP:                Optimal objective value is -4.000000.                                            


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value,
options.AbsoluteGapTolerance = 0 (the default value). The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05 (the default value).


x =

     3
     1


zmax =

     4

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

基于问题

clc,clear

% 建立问题
prob = optimproblem('ObjectiveSense','max');

% 创建决策变量
x=optimvar('x',2,'Type','integer','LowerBound',0);

% 目标函数
prob.Objective = x(1)+x(2);

% 约束条件
prob.Constraints.con1 = 14*x(1)+9*x(2)<=51;
prob.Constraints.con2 = -6*x(1)+11*x(2)<=1;

% 求解
[result,zmax]=solve(prob);

% 展示结果
display(result.x);
display(zmax);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

Solving problem using intlinprog.
LP:                Optimal objective value is -4.000000.                                            


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value,
options.AbsoluteGapTolerance = 0 (the default value). The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05 (the default value).

     3
     1


zmax =

     4
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

Python版

# 导库
import pulp
import numpy as np

# 初始化线性规划问题模型
pro = pulp.LpProblem()

# 创建整数型决策变量
lb = [0, 0]
ub = [None, None]
bounds = np.array(list(zip(lb, ub)))
x = np.array([
    pulp.LpVariable(f'x{i}',
                    lowBound=bounds[i - 1][0],
                    upBound=bounds[i - 1][1],
                    cat="Integer") for i in [1, 2]
])

#建立目标函数
c = -np.array([1, 1])
pro += pulp.lpDot(c, x)

# 不等式约束
A = np.array([[14, 9], [-6, 11]])
b = np.array([51, 1])
for i in range(len(A)):
    pro += (pulp.lpDot(A[i], x) <= b[i])  # 添加限制条件

#求解
pro.solve()

#输出结果
print(pro)
print(f'最优值：{-pulp.value(pro.objective)}')
print(f'最优解：{[pulp.value(var) for var in 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

NoName:
MINIMIZE
-1*x1 + -1*x2 + 0
SUBJECT TO
_C1: 14 x1 + 9 x2 <= 51

_C2: - 6 x1 + 11 x2 <= 1

VARIABLES
0 <= x1 Integer
0 <= x2 Integer

最优值：4.0
最优解：[3.0, 1.0]
1
2
3
4
5
6
7
8
9
10
11
12
13
14

• 重点在于创建决策变量时，加入cat='Integer'参数，限制决策变量