Python之建模规划篇--整数规划_整数规划 python 手搓

基本介绍

规划中的变量（部分或全部）限制为整数时，称为整数规划。若在线性规划模型中，变量限制为整数，则称为整数线性规划。目前所流行的求解整数规划的方法，往往只适用于整数线性规划。目前还没有一种方法能有效地求解一切整数规划。

整数规划的分类

如不加特殊说明，一般指整数线性规划。对于整数线性规划模型大致可分为两类：

• 变量全限制为整数时，称纯（完全）整数规划。
• 变量部分限制为整数的，称混合整数规划。

整数规划的特点

• 原线性规划有最优解，当自变量限制为整数后，其整数规划解出现下述情况：
  ①原线性规划最优解全是整数，则整数规划最优解与线性规划最优解一致。
  ②整数规划无可行解
  
• 整数规划最优解不能按照实数最优解简单取整而获得。

求解方法分类

• 分枝定界法—可求纯或混合整数线性规划。
• 割平面法—可求纯或混合整数线性规划。
• 隐枚举法—求解“0-1”整数规划：
  ①过滤隐枚举法；
  ②分枝隐枚举法。
• 匈牙利法—解决指派问题（“0-1”规划特殊情形）。
• 蒙特卡洛法—求解各种类型规划。

0 - 1 型整数规划

0 −1型整数规划是整数规划中的特殊情形，它的变量 x_j 仅取值0或1。这时x_j 称为0−1变量，或称二进制变量。x_j 仅取值0 或1 这个条件可由下述约束条件：

所代替，是和一般整数规划的约束条件形式一致的。在实际问题中，如果引入 0 −1变量，就可以把有各种情况需要分别讨论的线性规划问题统一在一个问题中讨论了。我们先介绍引入0 −1变量的实际问题，再研究解法。
比如有一些相互排斥的约束条件,就是一种0-1问题,如运输方式只能选择一种,用车或者用船等类似的
除此之外,还有关于固定费用的问题,在讨论线性规划时，有些问题是要求使成本为最小。那时总设固定成本为常数，并在线性规划的模型中不必明显列出。但有些固定费用（固定成本）的问题不能用一般线性规划来描述，但可改变为混合整数规划来解决

蒙特卡洛法 (随机取样法)

蒙特卡洛方法也称为计算机随机模拟方法,它源于世界著名的赌城一摩纳哥的Monte Carlo(蒙特卡洛)。它是基于对大量事件的统计结果来实现–些确定性问题的计算。使用蒙特卡洛方法必须使用计算机生成相关分布的随机数,Matlab和python等各种编程语言都给出了生成各种随机数的命令。
如,给个例子

前面介绍的常用的整数规划求解方法，主要是针对线性整数规划而言，而对于非线性整数规划目前尚未有一种成熟而准确的求解方法，因为非线性规划本身的通用有效解法尚未找到，更何况是非线性整数规划。
然而，尽管整数规划由于限制变量为整数而增加了难度；然而又由于整数解是有限个，于是为枚举法提供了方便。当然，当自变量维数很大和取值范围很宽情况下，企图用显枚举法（即穷举法）计算出最优值是不现实的，但是应用概率理论可以证明，在一定的计算量的情况下，完全可以得出一个满意解。
如

如果用显枚举法试探，共需计算(100)⁵ = 10¹⁰个点，其计算量非常之大。然而应用蒙特卡洛去随机计算10⁶个点，便可找到满意解，那么这种方法的可信度究竟怎样
呢？
下面就分析随机取样采集10⁶个点计算时，应用概率理论来估计一下可信度。
不失一般性，假定一个整数规划的最优点不是孤立的奇点。
假设目标函数落在高值区的概率分别为 0.01，0.00001，则当计算10⁶个点后，有
任一个点能落在高值区的概率分别为

首先编写M 文件mente.m 定义目标函数f 和约束向量函数g，程序如下：

function [f,g]=mengte(x);
f=x(1)^2+x(2)^2+3*x(3)^2+4*x(4)^2+2*x(5)-8*x(1)-2*x(2)-3*x(3)-...
x(4)-2*x(5);
g=[sum(x)-400
x(1)+2*x(2)+2*x(3)+x(4)+6*x(5)-800
2*x(1)+x(2)+6*x(3)-200
x(3)+x(4)+5*x(5)-200]
1
2
3
4
5
6
7

编写M文件mainint.m如下求问题的解

rand('state',sum(clock));  %初始化随机数发生器
p0=0;
tic    %计时开始
for i=1:10^6
   x=randi([0,99],1,5); %产生一行五列的区间[0,99]上的随机整数
   [f,g]=mengte(x);
   if all(g<=0)
       if p0<f
           x0=x; p0=f; %记录下当前较好的解
       end
   end
end
x0,p0
toc    %计时结束
1
2
3
4
5
6
7
8
9
10
11
12
13
14

由于是随机模拟,所以最后的答案都是不同的

整数线性规划的计算机求解

整数规划问题的求解使用Lingo等专用软件比较方便。对于整数线性规划问题,也可以使用Matlab的intlinprog函数求解,但使用Matlab软件求解数学规划问题有–个缺陷，即必须把所有的决策变量化成一-维决策向量,实际上对于多维变量的数学规划问题,用Matlab软件求解，需要做–个变量替换，把多维决策变量化成–维决策向量，变量替换后,约束条件很难写出;而使用Lingo软件求解数学规划问题是不需要做变换的,使用起来相对比较容易。


这里用一个指派问题作为例子进行解决

若用Matlab解决,代码如下

clc, clear
c=[3 8 2 10 3;8 7 2 9 7;6 4 2 7 5
   8 4 2 3 5;9 10 6 9 10];
c=c(:); a=zeros(10,25); intcon=1:25;
for i=1:5
   a(i,(i-1)*5+1:5*i)=1;
   a(5+i,i:5:25)=1;
end
b=ones(10,1); lb=zeros(25,1); ub=ones(25,1);
x=intlinprog(c,intcon,[],[],a,b,lb,ub);
x=reshape(x,[5,5])
1
2
3
4
5
6
7
8
9
10
11

这样可以得到最优的指派结果

分枝定界法

对有约束条件的最优化问题（其可行解为有限数）的所有可行解空间恰当地进行系统搜索，这就是分枝与定界内容。通常，把全部可行解空间反复地分割为越来越小的子集，称为分枝；并且对每个子集内的解集计算一个目标下界（对于最小值问题），这称为定界。在每次分枝后，凡是界限超出已知可行解集目标值的那些子集不再进一步分枝，这样，许多子集可不予考虑，这称剪枝。这就是分枝定界法的主要思路。

分枝定界法可用于解纯整数或混合的整数规划问题。在本世纪六十年代初由 Land Doig 和Dakin 等人提出的。由于这方法灵活且便于用计算机求解，所以现在它已是解整数规划的重要方法。目前已成功地应用于求解生产进度问题、旅行推销员问题、工厂选址问题、背包问题及分配问题等。
设有最大化的整数规划问题 A ，与它相应的线性规划为问题B ，从解问题B 开始，若其最优解不符合 A的整数条件，那么B的最优目标函数必是 A的最优目标函数z的上界，记作z1 ；而 A的任意可行解的目标函数值将是z的一个下界z2 。分枝定界法就是将B的可行域分成子区域的方法。逐步减小z 1和增大z2 ，最终求到z*。现用下例来说明：



从以上解题过程可得用分枝定界法求解整数规划（最大化）问题的步骤为：
开始，将要求解的整数规划问题称为问题 A ，将与它相应的线性规划问题称为问题B 。
（i）解问题B 可能得到以下情况之一：
（a） B 没有可行解，这时A 也没有可行解，则停止．
（b） B 有最优解，并符合问题A 的整数条件， B 的最优解即为A 的最优解，则停止。
（c） B 有最优解，但不符合问题A 的整数条件，记它的目标函数值为z1 。
（ii）用观察法找问题 A的一个整数可行解，一般可取x _j , j = 0,1.2,…,n,试探，求得其目标函数值，并记作z2。以z*表示问题 A的最优目标函数值；这时有

z2 ≤ z* ≤ z1 进行迭代。 第一步：分枝，在 B 的最优解中任选一个不符合整数条件的变量x _j ，其值为b _j，以b[ j ] 表示小于b _j的最大整数。构造两个约束条件 x _j ≤ [b _j] 和 x _j ≥ [b _j] + 1 将这两个约束条件，分别加入问题B ，求两个后继规划问题B ₁ 和B ₂。不考虑整数条件求解这两个后继问题。
定界，以每个后继问题为一分枝标明求解的结果，与其它问题的解的结果中，找出 最优目标函数值最大者作为新的上界z 1。从已符合整数条件的各分支中，找出目标函数 值为最大者作为新的下界z2，若无作用z 不变。
第二步：比较与剪枝，各分枝的最优目标函数中若有小于z2 者，则剪掉这枝，即 以后不再考虑了。若大于z2 ，且不符合整数条件，则重复第一步骤。一直到最后得到 z* = z2 为止。得最优整数解x _j* , j = 1,2,.....,n

Python 实现 (分支定界代码)

• 整数规划的模型与线性规划基本相同，只是额外增加了部分变量为整数的约束
• 整数规划求解的基本框架是分支定界法，首先去除整数约束得到“松弛模型”，使用线性规划的方法求解。
• 若有某个变量不是整数，在松弛模型.上分别添加约束: x≤floor(A)和x≥ceil(A),然后再分别求解，这个过程叫做分支。当节点求解结果中所有变量都是整数时，停止分支。这样不断迭代，形成了一颗树。
• 所谓定界，指的是叶子节点产生后，相当于给问题定了一个下界。之后在求解过程中一旦某个节点的目标函数值小于这个下界，那就直接pass，不再进行分支了;每次新产生叶子节点，则更新下界。

from scipy.optimize import linprog
import numpy as np
import math
import sys
from queue import Queue
 
 
class ILP():
    def __init__(self, c, A_ub, b_ub, A_eq, b_eq, bounds):
        # 全局参数
        self.LOWER_BOUND = -sys.maxsize
        self.UPPER_BOUND = sys.maxsize
        self.opt_val = None
        self.opt_x = None
        self.Q = Queue()
 
        # 这些参数在每轮计算中都不会改变
        self.c = -c
        self.A_eq = A_eq
        self.b_eq = b_eq
        self.bounds = bounds
 
        # 首先计算一下初始问题
        r = linprog(-c, A_ub, b_ub, A_eq, b_eq, bounds)
 
        # 若最初问题线性不可解
        if not r.success:
            raise ValueError('Not a feasible problem!')
 
        # 将解和约束参数放入队列
        self.Q.put((r, A_ub, b_ub))
 
    def solve(self):
        while not self.Q.empty():
            # 取出当前问题
            res, A_ub, b_ub = self.Q.get(block=False)
 
            # 当前最优值小于总下界，则排除此区域
            if -res.fun < self.LOWER_BOUND:
                continue
 
            # 若结果 x 中全为整数，则尝试更新全局下界、全局最优值和最优解
            if all(list(map(lambda f: f.is_integer(), res.x))):
                if self.LOWER_BOUND < -res.fun:
                    self.LOWER_BOUND = -res.fun
 
                if self.opt_val is None or self.opt_val < -res.fun:
                    self.opt_val = -res.fun
                    self.opt_x = res.x
 
                continue
 
            # 进行分枝
            else:
                # 寻找 x 中第一个不是整数的，取其下标 idx
                idx = 0
                for i, x in enumerate(res.x):
                    if not x.is_integer():
                        break
                    idx += 1
 
                # 构建新的约束条件（分割
                new_con1 = np.zeros(A_ub.shape[1])
                new_con1[idx] = -1
                new_con2 = np.zeros(A_ub.shape[1])
                new_con2[idx] = 1
                new_A_ub1 = np.insert(A_ub, A_ub.shape[0], new_con1, axis=0)
                new_A_ub2 = np.insert(A_ub, A_ub.shape[0], new_con2, axis=0)
                new_b_ub1 = np.insert(
                    b_ub, b_ub.shape[0], -math.ceil(res.x[idx]), axis=0)
                new_b_ub2 = np.insert(
                    b_ub, b_ub.shape[0], math.floor(res.x[idx]), axis=0)
 
                # 将新约束条件加入队列，先加最优值大的那一支
                r1 = linprog(self.c, new_A_ub1, new_b_ub1, self.A_eq,
                             self.b_eq, self.bounds)
                r2 = linprog(self.c, new_A_ub2, new_b_ub2, self.A_eq,
                             self.b_eq, self.bounds)
                if not r1.success and r2.success:
                    self.Q.put((r2, new_A_ub2, new_b_ub2))
                elif not r2.success and r1.success:
                    self.Q.put((r1, new_A_ub1, new_b_ub1))
                elif r1.success and r2.success:
                    if -r1.fun > -r2.fun:
                        self.Q.put((r1, new_A_ub1, new_b_ub1))
                        self.Q.put((r2, new_A_ub2, new_b_ub2))
                    else:
                        self.Q.put((r2, new_A_ub2, new_b_ub2))
                        self.Q.put((r1, new_A_ub1, new_b_ub1))
 
 
def test1():
    """ 此测试的真实最优解为 [4, 2] """
    c = np.array([40, 90])
    A = np.array([[9, 7], [7, 20]])
    b = np.array([56, 70])
    Aeq = None
    beq = None
    bounds = [(0, None), (0, None)]
 
    solver = ILP(c, A, b, Aeq, beq, bounds)
    solver.solve()
 
    print("Test 1's result:", solver.opt_val, solver.opt_x)
    print("Test 1's true optimal x: [4, 2]\n")
 
 
def test2():
    """ 此测试的真实最优解为 [2, 4] """
    c = np.array([3, 13])
    A = np.array([[2, 9], [11, -8]])
    b = np.array([40, 82])
    Aeq = None
    beq = None
    bounds = [(0, None), (0, None)]
 
    solver = ILP(c, A, b, Aeq, beq, bounds)
    solver.solve()
 
    print("Test 2's result:", solver.opt_val, solver.opt_x)
    print("Test 2's true optimal x: [2, 4]\n")
 
if __name__ == '__main__':
    test1()
    test2()
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

可以得到结果，就可以求出整数规划了（个人感觉整数规划没有matlab那么好，matlab有直接的函数）

每日一句
Never underestimate your power to change yourself!(永远不要低估你改变自我的能力! )