整数线性规划实现（lingo，python分枝界定法）

本文章为上篇建模学习打卡第二天的续

文章目录

一、本次问题

 二、本题理解

三、问题求解

1.lingo实现

（1）先抛除整数约束条件对问题求解

 （2）加入整数约束条件求解

2.python实现求解

（1）先抛除整数约束条件对问题求解

 （2）加入整数约束条件求解实现   通过 pulp 库求解

 （3）加入整数约束条件求解实现  分枝界定法求解

---

一、本次问题

 二、本题理解

目标函数：

max = 40x1+90x2

一级约束条件：

9x1+7x2<=56
7x1+20x2<=70
x1,x2 >= 0

二级约束条件：

x1，x2全为整数

三、问题求解

1.lingo实现

lingo编写代码时，每行代码结束后必须以 ‘ ; ’ 结束，否则无法运行。

（1）先抛除整数约束条件对问题求解

基础线性规划实现（matlab，lingo）_菜菜笨小孩的博客-CSDN博客

        lingo代码实现：（l无其他条件下，ingo中默认变量大于等于0）

max = 40*x1+90*x2;
9*x1+7*x2<=56;
7*x1+20*x2<=70;

        结果：最优解 x1=4.80916 , x2 = 1.816794 ; 最优值为355.8779；显然不符合题意

 （2）加入整数约束条件求解

首先，需要引出lingo的变量界定函数 @gin(x) --- 将x限制为整数条件

       ingo代码实现:通过变量界定函数将x1，x2限制为整数约束

max = 40*x1+90*x2;
9*x1+7*x2<=56;
7*x1+20*x2<=70;
@gin(x1);
@gin(x2);

        结果：最优解 x1=4 , x2 = 2 ; 最优值为340；符合题意

 lingo实现求解到此结束。

2.python实现求解

（1）先抛除整数约束条件对问题求解

基础线性规划实现---python_菜菜笨小孩的博客-CSDN博客

        python代码实现如下：详解请看上方python基础线性规划的文章

#导入包
from scipy import optimize as opt
import numpy as np

#确定c,A,b,Aeq,beq
c = np.array([40,90]) #目标函数变量系数
A = np.array([[9,7],[7,20]]) #不等式变量系数
b = np.array([56,70]) #不等式变量值
Aeq = np.array([[0,0]]) #等式变量系数
beq = np.array([0]) #等式变量值
#限制
lim1=(0,None)
lim2=(0,None)
#求解
res = opt.linprog(-c,A,b,Aeq,beq,bounds=(lim1,lim2))
#输出结果
print(res)

        结果：最优解 x1=4.80916 , x2 = 1.816794 ; 最优值为355.8779；显然不符合题意

 （2）加入整数约束条件求解实现   通过 pulp 库求解

安装库：我用python3.8安装成功了

pip install pulp

python代码实现：

1.导入库

import pulp as pulp

2.定义解决问题的函数

def solve_ilp(objective , constraints) :
    print(objective)
    print(constraints)
    prob = pulp.LpProblem('LP1' , pulp.LpMaximize)
    prob += objective
    for cons in constraints :
        prob += cons
    print(prob)
    status = prob.solve()
    if status != 1 :
        return None
    else :
        return [v.varValue.real for v in prob.variables()]

3.设置目标函数和约束条件等

V_NUM = 2 #本题变量个数
# 变量，直接设置下限
variables = [pulp.LpVariable('X%d' % i, lowBound=0, cat=pulp.LpInteger) for i in range(0, V_NUM)]
# 目标函数
c = [40, 90]
objective = sum([c[i] * variables[i] for i in range(0, V_NUM)])
# 约束条件
constraints = []
a1 = [9, 7]
constraints.append(sum([a1[i] * variables[i] for i in range(0, V_NUM)]) <= 56)
a2 = [7, 20]
constraints.append(sum([a2[i] * variables[i] for i in range(0, V_NUM)]) <= 70)

4.求解：

res = solve_ilp(objective, constraints)
print(res) #输出结果

完整代码如下：

# -*- coding: utf-8 -*-
import pulp as pulp

def solve_ilp(objective, constraints):
    print(objective)
    print(constraints)
    prob = pulp.LpProblem('LP1', pulp.LpMaximize)
    prob += objective
    for cons in constraints:
        prob += cons
    print(prob)
    status = prob.solve()
    if status != 1:
        # print 'status'
        # print status
        return None
    else:
        # return [v.varValue.real for v in prob.variables()]
        return [v.varValue.real for v in prob.variables()]

V_NUM = 2
# 变量，直接设置下限
variables = [pulp.LpVariable('X%d' % i, lowBound=0, cat=pulp.LpInteger) for i in range(0, V_NUM)]
# 目标函数
c = [40, 90]
objective = sum([c[i] * variables[i] for i in range(0, V_NUM)])
# 约束条件
constraints = []
a1 = [9, 7]
constraints.append(sum([a1[i] * variables[i] for i in range(0, V_NUM)]) <= 56)
a2 = [7, 20]
constraints.append(sum([a2[i] * variables[i] for i in range(0, V_NUM)]) <= 70)
# print (constraints)

res = solve_ilp(objective, constraints)
print(res) #输出解

结果：最优解 x1=4 , x2 = 2 ; 最优值为340；符合题意

 （3）加入整数约束条件求解实现  分枝界定法求解

何为分枝界定法，请看详解https://blog.csdn.net/qq_25990967/article/details/121211474

python代码实现：

1.导入库

from scipy.optimize import linprog
import numpy as np
import math
import sys
from queue import Queue

2.定义整数线性规划类

class ILP()

3.定义分枝界定法函数

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))

4.定义求解问题中的变量级约束条件

def test():
    """ 此测试的真实最优解为 [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 's result:", solver.opt_val, solver.opt_x)
    print("Test 's true optimal x: [4, 2]\n")

5.求解并输出

if __name__ == '__main__':
    test()

完整代码如下：

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 test():
    """ 此测试的真实最优解为 [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 's result:", solver.opt_val, solver.opt_x)
    print("Test 's true optimal x: [4, 2]\n")



if __name__ == '__main__':
    test()

结果：最优解 x1=4 , x2 = 2 ; 最优值为340；符合题意