本文章為上篇建模學(xué)習(xí)打卡第二天的續(xù)

文章目錄

一、本次問題

?二、本題理解

三、問題求解

1.lingo實(shí)現(xiàn)

（1）先拋除整數(shù)約束條件對(duì)問題求解

?（2）加入整數(shù)約束條件求解

2.python實(shí)現(xiàn)求解

（1）先拋除整數(shù)約束條件對(duì)問題求解

?（2）加入整數(shù)約束條件求解實(shí)現(xiàn)? ?通過 pulp 庫(kù)求解

?（3）加入整數(shù)約束條件求解實(shí)現(xiàn)? 分枝界定法求解

---

一、本次問題

?二、本題理解

目標(biāo)函數(shù)：

max = 40x1+90x2

一級(jí)約束條件：

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

二級(jí)約束條件：

x1，x2全為整數(shù)

三、問題求解

1.lingo實(shí)現(xiàn)

lingo編寫代碼時(shí)，每行代碼結(jié)束后必須以 ‘ ; ’ 結(jié)束，否則無法運(yùn)行。

（1）先拋除整數(shù)約束條件對(duì)問題求解

基礎(chǔ)線性規(guī)劃實(shí)現(xiàn)（matlab，lingo）_菜菜笨小孩的博客-CSDN博客

? ? ? ? lingo代碼實(shí)現(xiàn)：（l無其他條件下，ingo中默認(rèn)變量大于等于0）

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

????????結(jié)果：最優(yōu)解 x1=4.80916 , x2 = 1.816794 ; 最優(yōu)值為355.8779；顯然不符合題意

?（2）加入整數(shù)約束條件求解

首先，需要引出lingo的變量界定函數(shù) @gin(x) --- 將x限制為整數(shù)條件

? ? ? ?ingo代碼實(shí)現(xiàn):通過變量界定函數(shù)將x1，x2限制為整數(shù)約束

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

????????結(jié)果：最優(yōu)解 x1=4 , x2 = 2 ; 最優(yōu)值為340；符合題意

?lingo實(shí)現(xiàn)求解到此結(jié)束。

2.python實(shí)現(xiàn)求解

（1）先拋除整數(shù)約束條件對(duì)問題求解

基礎(chǔ)線性規(guī)劃實(shí)現(xiàn)---python_菜菜笨小孩的博客-CSDN博客

????????python代碼實(shí)現(xiàn)如下：詳解請(qǐng)看上方python基礎(chǔ)線性規(guī)劃的文章

#導(dǎo)入包 from scipy import optimize as opt import numpy as np#確定c,A,b,Aeq,beq c = np.array([40,90]) #目標(biāo)函數(shù)變量系數(shù) A = np.array([[9,7],[7,20]]) #不等式變量系數(shù) b = np.array([56,70]) #不等式變量值 Aeq = np.array([[0,0]]) #等式變量系數(shù) beq = np.array([0]) #等式變量值 #限制 lim1=(0,None) lim2=(0,None) #求解 res = opt.linprog(-c,A,b,Aeq,beq,bounds=(lim1,lim2)) #輸出結(jié)果 print(res)

????????結(jié)果：最優(yōu)解 x1=4.80916 , x2 = 1.816794 ; 最優(yōu)值為355.8779；顯然不符合題意

?（2）加入整數(shù)約束條件求解實(shí)現(xiàn) ? 通過 pulp 庫(kù)求解

安裝庫(kù)：我用python3.8安裝成功了

pip install pulp

python代碼實(shí)現(xiàn)：

1.導(dǎo)入庫(kù)

import pulp as pulp

2.定義解決問題的函數(shù)

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

3.設(shè)置目標(biāo)函數(shù)和約束條件等

V_NUM = 2 #本題變量個(gè)數(shù) # 變量，直接設(shè)置下限 variables = [pulp.LpVariable('X%d' % i, lowBound=0, cat=pulp.LpInteger) for i in range(0, V_NUM)] # 目標(biāo)函數(shù) 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) #輸出結(jié)果

完整代碼如下：

# -*- coding: utf-8 -*- import pulp as pulpdef solve_ilp(objective, constraints):print(objective)print(constraints)prob = pulp.LpProblem('LP1', pulp.LpMaximize)prob += objectivefor cons in constraints:prob += consprint(prob)status = prob.solve()if status != 1:# print 'status'# print statusreturn Noneelse:# return [v.varValue.real for v in prob.variables()]return [v.varValue.real for v in prob.variables()]V_NUM = 2 # 變量，直接設(shè)置下限 variables = [pulp.LpVariable('X%d' % i, lowBound=0, cat=pulp.LpInteger) for i in range(0, V_NUM)] # 目標(biāo)函數(shù) 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) #輸出解

結(jié)果：最優(yōu)解 x1=4 , x2 = 2 ; 最優(yōu)值為340；符合題意

?（3）加入整數(shù)約束條件求解實(shí)現(xiàn)? 分枝界定法求解

何為分枝界定法，請(qǐng)看詳解https://blog.csdn.net/qq_25990967/article/details/121211474

python代碼實(shí)現(xiàn)：

1.導(dǎo)入庫(kù)

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

2.定義整數(shù)線性規(guī)劃類

class ILP()

3.定義分枝界定法函數(shù)

def __init__(self, c, A_ub, b_ub, A_eq, b_eq, bounds):# 全局參數(shù)self.LOWER_BOUND = -sys.maxsizeself.UPPER_BOUND = sys.maxsizeself.opt_val = Noneself.opt_x = Noneself.Q = Queue()# 這些參數(shù)在每輪計(jì)算中都不會(huì)改變self.c = -cself.A_eq = A_eqself.b_eq = b_eqself.bounds = bounds# 首先計(jì)算一下初始問題r = linprog(-c, A_ub, b_ub, A_eq, b_eq, bounds)# 若最初問題線性不可解if not r.success:raise ValueError('Not a feasible problem!')# 將解和約束參數(shù)放入隊(duì)列self.Q.put((r, A_ub, b_ub))def solve(self):while not self.Q.empty():# 取出當(dāng)前問題res, A_ub, b_ub = self.Q.get(block=False)# 當(dāng)前最優(yōu)值小于總下界，則排除此區(qū)域if -res.fun < self.LOWER_BOUND:continue# 若結(jié)果 x 中全為整數(shù)，則嘗試更新全局下界、全局最優(yōu)值和最優(yōu)解if all(list(map(lambda f: f.is_integer(), res.x))):if self.LOWER_BOUND < -res.fun:self.LOWER_BOUND = -res.funif self.opt_val is None or self.opt_val < -res.fun:self.opt_val = -res.funself.opt_x = res.xcontinue# 進(jìn)行分枝else:# 尋找 x 中第一個(gè)不是整數(shù)的，取其下標(biāo) idxidx = 0for i, x in enumerate(res.x):if not x.is_integer():breakidx += 1# 構(gòu)建新的約束條件（分割new_con1 = np.zeros(A_ub.shape[1])new_con1[idx] = -1new_con2 = np.zeros(A_ub.shape[1])new_con2[idx] = 1new_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)# 將新約束條件加入隊(duì)列，先加最優(yōu)值大的那一支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.定義求解問題中的變量級(jí)約束條件

def test():""" 此測(cè)試的真實(shí)最優(yōu)解為 [4, 2] """c = np.array([40, 90])A = np.array([[9, 7], [7, 20]])b = np.array([56, 70])Aeq = Nonebeq = Nonebounds = [(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 Queueclass ILP():def __init__(self, c, A_ub, b_ub, A_eq, b_eq, bounds):# 全局參數(shù)self.LOWER_BOUND = -sys.maxsizeself.UPPER_BOUND = sys.maxsizeself.opt_val = Noneself.opt_x = Noneself.Q = Queue()# 這些參數(shù)在每輪計(jì)算中都不會(huì)改變self.c = -cself.A_eq = A_eqself.b_eq = b_eqself.bounds = bounds# 首先計(jì)算一下初始問題r = linprog(-c, A_ub, b_ub, A_eq, b_eq, bounds)# 若最初問題線性不可解if not r.success:raise ValueError('Not a feasible problem!')# 將解和約束參數(shù)放入隊(duì)列self.Q.put((r, A_ub, b_ub))def solve(self):while not self.Q.empty():# 取出當(dāng)前問題res, A_ub, b_ub = self.Q.get(block=False)# 當(dāng)前最優(yōu)值小于總下界，則排除此區(qū)域if -res.fun < self.LOWER_BOUND:continue# 若結(jié)果 x 中全為整數(shù)，則嘗試更新全局下界、全局最優(yōu)值和最優(yōu)解if all(list(map(lambda f: f.is_integer(), res.x))):if self.LOWER_BOUND < -res.fun:self.LOWER_BOUND = -res.funif self.opt_val is None or self.opt_val < -res.fun:self.opt_val = -res.funself.opt_x = res.xcontinue# 進(jìn)行分枝else:# 尋找 x 中第一個(gè)不是整數(shù)的，取其下標(biāo) idxidx = 0for i, x in enumerate(res.x):if not x.is_integer():breakidx += 1# 構(gòu)建新的約束條件（分割new_con1 = np.zeros(A_ub.shape[1])new_con1[idx] = -1new_con2 = np.zeros(A_ub.shape[1])new_con2[idx] = 1new_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)# 將新約束條件加入隊(duì)列，先加最優(yōu)值大的那一支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():""" 此測(cè)試的真實(shí)最優(yōu)解為 [4, 2] """c = np.array([40, 90])A = np.array([[9, 7], [7, 20]])b = np.array([56, 70])Aeq = Nonebeq = Nonebounds = [(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()

結(jié)果：最優(yōu)解 x1=4 , x2 = 2 ; 最優(yōu)值為340；符合題意

總結(jié)

以上是生活随笔為你收集整理的整数线性规划实现（lingo，python分枝界定法）的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。