本文章为上篇建模学习打卡第二天的续
文章目录
一、本次问题
二、本题理解
三、问题求解
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）先抛除整数约束条件对问题求解
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代码实现如下：详解请看上方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 pulppython代码实现：
1.导入库
import pulp as pulp2.定义解决问题的函数
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）加入整数约束条件求解实现 分枝界定法求解
何为分枝界定法，请看详解
python代码实现：
1.导入库
from scipy.optimize import linprog
import numpy as np
import math
import sys
from queue import Queue2.定义整数线性规划类
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；符合题意
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
