1. 全栈程序员必看首页

python线性规划模型求解_python线性规划算法

全栈程序员-用户IM • • 未分类

python线性规划模型求解_python线性规划算法各位童鞋们大家好，我是小小明，前几天我给大家分享了一个SMT求解器z3，链接地址见：https://xxmdmst.blog.csdn.net/article/details/120279521虽然SMT求解器很强大，能够解逻辑题、解数独、解方程、甚至解决逆向问题，但是有个缺点就是只能找出一个可行解，如果我想要找出可行解的最大值或最小值就不行，无法完成类似Excel的规划求解的功能。前文中已经提到了scipy这个库可以进行线性规划求解，可惜我在这周的实际测试中发现，不支持整数约束，只能求解出实数。差.

大家好，又见面了，我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全系列IDE使用 1年只要46元 售后保障 童叟无欺

python线性规划模型求解_python线性规划算法

各位童鞋们大家好，我是小小明，前几天我给大家分享了一个SMT求解器z3，链接地址见：

https://xxmdmst.blog.csdn.net/article/details/120279521

虽然SMT求解器很强大，能够解逻辑题、解数独、解方程、甚至解决逆向问题，但是有个缺点就是只能找出一个可行解，如果我想要找出可行解的最大值或最小值就不行，无法完成类似Excel的规划求解的功能。

前文中已经提到了scipy这个库可以进行线性规划求解，可惜我在这周的实际测试中发现，不支持整数约束，只能求解出实数。差点放弃写这篇文章，不过后面我又发现了PuLP这个库，简直发现了新大陆，原来有这么个专门进行规划求解的库，于是又花了几天时间专门来研究它。

今天呢，我将首先演示scipy这个库的用法，然后再演示PuLP库如何进行线性规划求解。

当然，PuLP库也不是万能的，虽然可以解决线性规划问题，但是相对Excel的规划求解还有很多难以实现的功能，例如非线性的规划求解。如果真的必须要进行非线性的规划求解时，可以考虑研究一下CVXPY这个库，官网地址是：https://www.cvxpy.org/

当然对于规划求解，95%以上的场景都是线性规划求解，PuLP就足够应对我们需要应对的场景。

认真看看本文噢，这都是近50小时研究，分享出来的精华噢，基本上你花一个小时就能掌握50小时的心血，而且本文文末还有技术彩蛋，带你掌握更多简单高效的技巧。

文章目录

使用scipy库进行规划求解

scipy与numpy离线文档的下载地址：https://docs.scipy.org/doc/

在线文档地址：https://docs.scipy.org/doc/scipy/reference/

scipy库支持线性规划的函数是scipy.optimize.linprog

函数申明为：

cipy.optimize.linprog(c, A_ub = None, b_ub = None, A_eq = None, b_eq = None, bounds = None, method = ‘interior-point’, callback = None, options = None, x0 = None)

该函数用于在如下约束条件下，寻找目标函数( c T x c^{T} x cTx)的最小值：
A u b x ≤ b u b A e q x = b e q l ≤ x ≤ u \begin{array}{r} A_{u b} x \leq b_{u b} \\ A_{e q} x=b_{e q} \\ l \leq x \leq u \end{array} Aub​x≤bub​Aeq​x=beq​l≤x≤u​
主要参数说明：

  • c：目标函数的未知数系数 {一维数组}
  • A_up：不等式约束的未知数系数 {二维数组}
  • b_up：不等式约束值的上限 {一维数组}
  • A_eq：等式约束的未知数系数 {二维数组}
  • B_eq：等式约束值的上限 {一维数组}
  • bounds：决策变量的最小值和最大值

示例1：

假设要求目标函数 z = 4 x 1 + 3 x 2 z=4x_1+3x_2 z=4x1​+3x2​的最大值。

约束条件为：
{ 2 x 1 + x 2 ≤ 10 x 1 + x 2 ≤ 8 x 2 ≤ 7 x 1 , x 2 ≥ 0 \left\{\begin{array}{l} 2 x_{1}+x_{2} \leq 10 \\ x_{1}+x_{2} \leq 8 \\ x_{2} \leq 7 \\ x_{1}, x_{2} \geq 0 \end{array}\right. ⎩⎪⎪⎨⎪⎪⎧​2x1​+x2​≤10x1​+x2​≤8x2​≤7x1​,x2​≥0​

可以将目标函数理解为求-z的最小值，然后转换为linprog需要的参数传入：

import numpy as np
from scipy import optimize as op

z = np.array([4, 3])

A_ub = np.array([[2, 1], [1, 1]])
b_ub = np.array([10, 8])
# 元组分布表示x1和x2的边界范围
bounds = [0, None], [0, 7]
# 求目标函数取反后的最小值
res = op.linprog(-z, A_ub, b_ub, bounds=bounds)
print(f"目标函数的最大值z={ 
     -res.fun:.2f}，此时目标函数的决策变量为{ 
     res.x.round(2)}")

目标函数的最大值z=26.00，此时目标函数的决策变量为[2. 6.]

示例2：

目标函数： max ⁡ z = 2 x 1 + 3 x 2 − 5 x 3 \max z=2 x_{1}+3 x_{2}-5 x_{3} maxz=2x1​+3x2​−5x3​

约束条件：
{ x 1 + x 2 + x 3 = 7 2 x 1 − 5 x 2 + x 3 ≥ 10 x 1 + 3 x 2 + x 3 ≤ 12 x 1 , x 2 , x 3 ≥ 0 \left\{\begin{array}{l} x_{1}+x_{2}+x_{3}=7 \\ 2 x_{1}-5 x_{2}+x_{3} \geq 10 \\ x_{1}+3 x_{2}+x_{3} \leq 12 \\ x_{1}, x_{2}, x_{3} \geq 0 \end{array}\right. ⎩⎪⎪⎨⎪⎪⎧​x1​+x2​+x3​=72x1​−5x2​+x3​≥10x1​+3x2​+x3​≤12x1​,x2​,x3​≥0​

首先转换一下约束条件为：

{ − 2 x 1 + 5 x 2 − x 3 ≤ − 10 x 1 + 3 x 2 + x 3 ≤ 12 x 1 + x 2 + x 3 = 7 x 1 , x 2 , x 3 ≥ 0 \left\{\begin{array}{l} -2 x_{1}+5 x_{2}-x_{3} \leq -10 \\ x_{1}+3 x_{2}+x_{3} \leq 12 \\ x_{1}+x_{2}+x_{3}=7 \\ x_{1}, x_{2}, x_{3} \geq 0 \end{array}\right. ⎩⎪⎪⎨⎪⎪⎧​−2x1​+5x2​−x3​≤−10x1​+3x2​+x3​≤12x1​+x2​+x3​=7x1​,x2​,x3​≥0​

然后编写出如下代码：

z = np.array([2, 3, -5])

A_ub = np.array([[-2, 5, -1], [1, 3, 1]])
b_ub = np.array([-10, 12])
A_eq = np.array([[1, 1, 1]])
b_eq = np.array([7])
# 元组分别表示x1、x2和x3的边界范围
bounds = [0, None], [0, None], [0, None]
# 求目标函数取反后的最小值即目标函数的最大值
res = op.linprog(-z, A_ub, b_ub, A_eq, b_eq, bounds=bounds)
print(f"目标函数的最大值z={ 
     -res.fun:.2f}，此时目标函数的决策变量为{ 
     res.x.round(2)}")

结果：

目标函数的最大值z=14.57，此时目标函数的决策变量为[6.43 0.57 0.  ]

此时假如需要将变量约束为整数，scipy库就无能为力。

咱们再介绍专门用来进行线性规划的库PuLP的使用：

使用PuLP进行规划求解

PULP是一个线性规划(Linear Programming, LP)问题求解库。它将优化问题描述为数学模型，生成MPS或者LP文件，然后调用LP求解器，如CBC、GLPK、CPLEX、Gurobi等来进行求解。

安装完pulp库默认就拥有了CBC求解器，其他求解器需要额外安装才能使用。

安装pulp库：

pip install pulp

详见：https://pypi.org/project/PuLP/

官方文档：https://coin-or.github.io/pulp/

使用步骤：

  1. 通过pulp.LpProblem建立优化问题：
prob = pulp.LpProblem(name='NoName', sense=LpMinimize)

参数：

  • name：在lp文件中写入的问题名称；
  • sense：最大或者最小，可为LpMinimize\LpMaximize二者之一。
  1. 通过LpVariable函数创建优化变量：
pulp.LpVariable(name, lowBound=None, upBound=None, cat='Continuous', e=None)

参数：

  • name：变量名，会输出到.lp文件；
  • lowBound：变量下界；
  • upBound：变量上界；
  • cat：变量类型，可以为LpInteger\LpBinary\LpContinuous三者之一（整型、二进制或连续）；
  • e：指明变量是否在目标函数和约束中存在，主要用来实现列生成算法。
  1. 通过+=向优化问题添加目标函数和约束，例如：
prob += x1 + x2 == 100, "目标函数"
prob += 0.100*x1 + 0.200*x2 >= 8.0, "约束1"
  1. 使用求解器求解，默认为CBC：
prob.solve()

假如本机已经安装其他的LP求解器，就可以指定求解器，例如GLPK求解器：

prob.solve(GLPK(msg = 0))

GLPK求解器的下载地址：https://sourceforge.net/projects/winglpk/files/latest/download

假如需要使用默认的CBC求解MIP问题，可以指定use_mps=False：

prob.solve(use_mps=False)

下面我们看看如何使用PuLP解决上面的要求整数约束的一个线性规划问题：

目标函数:
max ⁡ z = 2 x 1 + 3 x 2 − 5 x 3 \max z=2 x_{1}+3 x_{2}-5 x_{3} maxz=2x1​+3x2​−5x3​

约束条件：

{ x 1 + x 2 + x 3 = 7 2 x 1 − 5 x 2 + x 3 ≥ 10 x 1 + 3 x 2 + x 3 ≤ 12 x 1 , x 2 , x 3 ≥ 0 \left\{\begin{array}{l} x_{1}+x_{2}+x_{3}=7 \\ 2 x_{1}-5 x_{2}+x_{3} \geq 10 \\ x_{1}+3 x_{2}+x_{3} \leq 12 \\ x_{1}, x_{2}, x_{3} \geq 0 \end{array}\right. ⎩⎪⎪⎨⎪⎪⎧​x1​+x2​+x3​=72x1​−5x2​+x3​≥10x1​+3x2​+x3​≤12x1​,x2​,x3​≥0​

代码如下：

from pulp import *

prob = LpProblem('max_z', sense=LpMaximize)
x1 = LpVariable('x1', 0, None, LpInteger)
x2 = LpVariable('x2', 0, None, LpInteger)
x3 = LpVariable('x3', 0, None, LpInteger)

# 设置目标函数
prob += 2*x1+3*x2-5*x3
# 约束条件
prob += x1+x2+x3 == 7
prob += 2*x1-5*x2+x3 >= 10
prob += x1+3*x2+x3 <= 12

status = prob.solve()
# print("求解状态:", LpStatus[prob.status])
print(f"目标函数的最大值z={ 
     value(prob.objective)}，此时目标函数的决策变量为:",
      { 
   v.name: v.varValue for v in prob.variables()})

结果：

目标函数的最大值z=14.0，此时目标函数的决策变量为: {'x1': 7.0, 'x2': 0.0, 'x3': 0.0}

再举个例子：

生产两种风机（风机A和风机B）
生产风机A，需要工时3小时，用电4千瓦，钢材9吨；
生产风机B，需要工时7小时，用电5千瓦，钢材5吨。
公司可提供的工时为300小时，可提供的用电量为250千瓦，可提供的钢材为420吨。
假设，两种产品的单位利润分别为200万元和210万元。怎样安排两种产品的生产量，所获得的利润最大？

解：设风机A产量为x，风机B产量为y，最大利润为Pmax,则

x,y>=0
3x+7y<=300
4x+5y<=250
9x+5y<=420
Pmax=200x+210y

由于生产量必须为整数，所以必须有整数约束，编码如下：

from pulp import *

prob = LpProblem('myProblem', sense=LpMaximize)

x = LpVariable('x', 0, 100, LpInteger)
y = LpVariable('y', 0, 100, LpInteger)

# 设置目标函数
prob += 200*x+210*y
# 约束条件
prob += 3*x+7*y <= 300
prob += 4*x+5*y <= 250
prob += 9*x+5*y <= 420

status = prob.solve()
# print("求解状态:", LpStatus[prob.status])
print(f"目标函数的最大值z={ 
     value(prob.objective)}，此时目标函数的决策变量为:",
      { 
   v.name: v.varValue for v in prob.variables()})

目标函数的最大值z=11460.0，此时目标函数的决策变量为: {'x': 30.0, 'y': 26.0}

可以看到pulp已经在整数约束下求解出了该问题。

关于LpStatus涉及的常量：

LpStatus = { 
    
LpStatusNotSolved:"Not Solved",
LpStatusOptimal:"Optimal",
LpStatusInfeasible:"Infeasible",
LpStatusUnbounded:"Unbounded",
LpStatusUndefined:"Undefined",
}
  • Not Solved：还没有调研solve()函数前的状态。
  • Optimal：找到了最优解。
  • Infeasible：问题没有可行解(比如定义了constraints x <= 1并且x >=2这样的约束)。
  • Unbounded：约束条件是无界的(not bounded)，最大化会导致无穷大(比如只有一个x >= 3这样的约束)。
  • Undefined：最优解可能存在但是没有求解出来。

下面我们多学习编码几个案例巩固所学的知识：

PuLP求解线性规划案例

组合药物配料使总成本最低

某药厂生产A、B、C三种药物，可供选择的原料有甲、乙、丙、丁四种，成本分别是毎公斤5元、6元、7元、8元。每公斤不同原料所能提供的各种药物用量如下。若该药厂每天要生产A药品恰好100克，B药品至少530克，C药品至多160克。怎样选择各种配料才能满足生产又能使总成本最低？

image-20210916110833532

首先，我们计算出目标函数和约束条件。

目标函数：
min ⁡ z = 5 x 1 + 6 x 2 + 7 x 3 + 8 x 4 \min z=5x_{1}+6x_{2}+7x_{3}+8×4 minz=5x1​+6x2​+7x3​+8x4

约束条件：

{ x 1 + x 2 + x 3 + x 4 = 100 5 x 1 + 4 x 2 + 5 x 3 + 6 x 4 ≥ 530 2 x 1 + x 2 + x 3 + 2 x 4 ≤ 160 x 1 , x 2 , x 3 , x 4 ≥ 0 \left\{\begin{array}{l} x_{1}+x_{2}+x_{3}+x_4=100 \\ 5x_1+4x_2+5x_3+6x_4 \geq 530 \\ 2x_1+x_2+x_3+2x_4 \leq 160 \\ x_1, x_2, x_3, x_4 \geq 0 \end{array}\right. ⎩⎪⎪⎨⎪⎪⎧​x1​+x2​+x3​+x4​=1005x1​+4x2​+5x3​+6x4​≥5302x1​+x2​+x3​+2x4​≤160x1​,x2​,x3​,x4​≥0​

下面我们使用python进行求解：

from pulp import *

prob = LpProblem('目标函数和约束', sense=LpMinimize)

x1 = LpVariable('x1', 0, None, LpInteger)
x2 = LpVariable('x2', 0, None, LpInteger)
x3 = LpVariable('x3', 0, None, LpInteger)
x4 = LpVariable('x4', 0, None, LpInteger)

# 设置目标函数
prob += 5*x1+6*x2+7*x3+8*x4
# 约束条件
prob += x1+x2+x3+x4 == 100
prob += 5*x1+4*x2+5*x3+6*x4 >= 530
prob += 2*x1+x2+x3+2*x4 <= 160

print(prob)
status = prob.solve(use_mps=True)
print(f"目标函数的最小值z={ 
     value(prob.objective)}，此时目标函数的决策变量为:",
      { 
   v.name: v.varValue for v in prob.variables()})

结果：

目标函数和约束:
MINIMIZE
5*x1 + 6*x2 + 7*x3 + 8*x4 + 0
SUBJECT TO
_C1: x1 + x2 + x3 + x4 = 100

_C2: 5 x1 + 4 x2 + 5 x3 + 6 x4 >= 530

_C3: 2 x1 + x2 + x3 + 2 x4 <= 160

VARIABLES
0 <= x1 Integer
0 <= x2 Integer
0 <= x3 Integer
0 <= x4 Integer

目标函数的最小值z=670.0，此时目标函数的决策变量为: {'x1': 30.0, 'x2': 0.0, 'x3': 40.0, 'x4': 30.0}

医院每天护士人数分配

某县新建一家医院，根据各个科室要求需要配备护士，周一到周日分别最小需要34、25、36、30、28、31、32人，按照规定，一个护士一周要连续上班5天。这家医院至少需要多少个护士？

解答：

一个护士一周要连续上班5天，我们假设有从周一开始连续上班的护士人数为x1，周二x2，…，以此类推。对于周一这天，除了从周二和周三开始上班的护士外，其他护士都会在这天上班，其他日期也同理。那么我们就可以用线性规划解决很简单的这个问题：

from pulp import *

# 设周1到周日开始上班的护士人数
x = [LpVariable(f"x{ 
     i}", lowBound=0, upBound=18, cat=LpInteger)
     for i in range(1, 8)]
min_nums = [34, 25, 36, 30, 28, 31, 32]
prob = LpProblem('目标函数和约束', sense=LpMinimize)
prob += lpSum(x)
for i, num in enumerate(x):
    prob += lpSum(x)-x[(i+5) % 7]-x[(i+6) % 7] >= min_nums[i]
status = prob.solve()
print("最少护士人数 z=", value(prob.objective))

print("周1到周日开始上班的护士人数分别为：", [value(x[i]) for i in range(7)])
print("周一到周日上班人数分别为：", [
      value(lpSum(x)-x[(i+5) % 7]-x[(i+6) % 7]) for i in range(7)])

最少护士人数 z= 44.0
周1到周日开始上班的护士人数分别为： [8.0, 0.0, 14.0, 1.0, 12.0, 0.0, 9.0]
周一到周日上班人数分别为： [35.0, 27.0, 36.0, 30.0, 29.0, 31.0, 32.0]

结论上最少需要44个护士，当然护士的分配这并不是唯一解，还有很多种分配方式，上面给出了其中一种可行的分配方式。

数据包络分析(DEA)

数据包络分析方法(Data Envelopment Analysis,DEA)是运筹学、管理科学与数理经济学交叉研究的一个新领域。它是根据多项投入指标和多项产出指标，利用线性规划的方法，对具有可比性的同类型单位进行相对有效性评价的一种数量分析方法。

例如，某银行的某区域有六个分理处，它们的投入产出情况如下：

image-20210916145836997

现在需要分析分理处1是否可以通过其他分理处的组合达到投入更少产出更多的效果。

对于这个问题，其实还比较复杂，我们可以先看看如何用Excel解决这个问题。

首先补充Excel表为如下结构：

image-20210916182250664

（投入）B11单元格的公式为：=VLOOKUP($A11,$A$3:$F$8,COLUMN(),0)*$H$9

（产出）D11单元格的公式为：=VLOOKUP($A11,$A$3:$F$8,COLUMN(),0)

image-20210916182347424

（组合后的投入产出）B12单元格的公式为：=SUMPRODUCT(B$3:B$8,$H$3:$H$8)

image-20210916182621750

x之和设置为六个分理处的权重之和，H2单元格的公式为：=SUM(H3:H8)

然后设置规划求解参数：

image-20210916183822005

最终结果：

image-20210916184017986

那么我们如何使用Python来解决这个问题：

首先读取数据：

import pandas as pd
import numpy as np

df = pd.read_excel(r"C:\Users\ASUS\Desktop\规划求解.xlsx",
                   index_col=0, usecols="A:F", sheet_name=1,
                   skiprows=1, nrows=6)
data = df.values
data

array([[  20,  149, 1300,  636, 1570],
       [  18,  152, 1500,  737, 1730],
       [  23,  140, 1500,  659, 1320],
       [  22,  142, 1500,  635, 1420],
       [  22,  129, 1200,  626, 1660],
       [  25,  142, 1600,  775, 1590]], dtype=int64)

然后开始求解：

from pulp import *

prob = LpProblem('目标函数和约束', sense=LpMinimize)
# 定义6个变量
x = [LpVariable(f"x{ 
     i}", lowBound=0, upBound=1) for i in range(1, 7)]
# 定义期望E
e = LpVariable("e", lowBound=0, upBound=1)

# 办事处1的数据
office1 = data[0]
# 各办事处的加权平均值
office_wavg = np.sum(data*np.array(x)[:, None], axis=0)

# 定义目标变量，期望存储在e变量中
prob += e
# 权重之和为1
prob += lpSum(x) == 1
# 投入更少
for i in range(2):
    prob += office_wavg[i] <= office1[i]*e
# 产出更多
for i in range(2, data.shape[1]):
    prob += office_wavg[i] >= office1[i]

print(prob)
status = prob.solve()
print("求解状态:", LpStatus[prob.status])
print(f"目标函数的最小值z={ 
     value(prob.objective)}，此时目标函数的决策变量为:",
      { 
   v.name: v.varValue for v in prob.variables()})
print("组合后的投入和产出：",[f"{ 
     value(office_wavg[i]):.2f}" for i in range(data.shape[1])])

结果：

目标函数和约束:
MINIMIZE
1*e + 0
SUBJECT TO
_C1: x1 + x2 + x3 + x4 + x5 + x6 = 1

_C2: - 20 e + 20 x1 + 18 x2 + 23 x3 + 22 x4 + 22 x5 + 25 x6 <= 0

_C3: - 149 e + 149 x1 + 152 x2 + 140 x3 + 142 x4 + 129 x5 + 142 x6 <= 0

_C4: 1300 x1 + 1500 x2 + 1500 x3 + 1500 x4 + 1200 x5 + 1600 x6 >= 1300

_C5: 636 x1 + 737 x2 + 659 x3 + 635 x4 + 626 x5 + 775 x6 >= 636

_C6: 1570 x1 + 1730 x2 + 1320 x3 + 1420 x4 + 1660 x5 + 1590 x6 >= 1570

VARIABLES
e <= 1 Continuous
x1 <= 1 Continuous
x2 <= 1 Continuous
x3 <= 1 Continuous
x4 <= 1 Continuous
x5 <= 1 Continuous
x6 <= 1 Continuous

求解状态: Optimal
目标函数的最小值z=0.96780303，此时目标函数的决策变量为: {'e': 0.96780303, 'x1': 0.0, 'x2': 0.66098485, 'x3': 0.0, 'x4': 0.0, 'x5': 0.33901515, 'x6': 0.0}
组合后的投入和产出： ['19.36', '144.20', '1398.30', '699.37', '1706.27']

可以看到求解出来的结果与Excel一致。

当然，针对任何一个分理处进行计算都可以得到一致的结果：

from pulp import *

# 定义6个变量
x = [LpVariable(f"x{ 
     i}", lowBound=0, upBound=1) for i in range(1, 7)]
# 定义期望E
e = LpVariable("e", lowBound=0, upBound=1)

# 办事处1的数据
for num, office in enumerate(data, 1):
    prob = LpProblem('目标函数和约束', sense=LpMinimize)
    # 各办事处的加权平均值
    office_wavg = np.sum(data*np.array(x)[:, None], axis=0)

    # 定义目标变量，期望存储在e变量中
    prob += e
    # 权重之和为1
    prob += lpSum(x) == 1
    # 投入更少
    for i in range(2):
        prob += office_wavg[i] <= office1[i]*e
    # 产出更多
    for i in range(2, data.shape[1]):
        prob += office_wavg[i] >= office1[i]

    status = prob.solve()
    print(f"---------办事处{ 
     num}----------")
    print("求解状态:", LpStatus[prob.status])
    print(f"目标函数的最小值z={ 
     value(prob.objective)}，此时目标函数的决策变量为:",
          { 
   v.name: v.varValue for v in prob.variables()})
    print("组合后的投入和产出：", [
          f"{ 
     value(office_wavg[i]):.2f}" for i in range(data.shape[1])])

---------办事处1----------
目标函数的最小值z=0.96780303，此时目标函数的决策变量为: {'e': 0.96780303, 'x1': 0.0, 'x2': 0.66098485, 'x3': 0.0, 'x4': 0.0, 'x5': 0.33901515, 'x6': 0.0}
组合后的投入和产出： ['19.36', '144.20', '1398.30', '699.37', '1706.27']
---------办事处2----------
目标函数的最小值z=0.96780303，此时目标函数的决策变量为: {'e': 0.96780303, 'x1': 0.0, 'x2': 0.66098485, 'x3': 0.0, 'x4': 0.0, 'x5': 0.33901515, 'x6': 0.0}
组合后的投入和产出： ['19.36', '144.20', '1398.30', '699.37', '1706.27']
---------办事处3----------
目标函数的最小值z=0.96780303，此时目标函数的决策变量为: {'e': 0.96780303, 'x1': 0.0, 'x2': 0.66098485, 'x3': 0.0, 'x4': 0.0, 'x5': 0.33901515, 'x6': 0.0}
组合后的投入和产出： ['19.36', '144.20', '1398.30', '699.37', '1706.27']
---------办事处4----------
目标函数的最小值z=0.96780303，此时目标函数的决策变量为: {'e': 0.96780303, 'x1': 0.0, 'x2': 0.66098485, 'x3': 0.0, 'x4': 0.0, 'x5': 0.33901515, 'x6': 0.0}
组合后的投入和产出： ['19.36', '144.20', '1398.30', '699.37', '1706.27']
---------办事处5----------
目标函数的最小值z=0.96780303，此时目标函数的决策变量为: {'e': 0.96780303, 'x1': 0.0, 'x2': 0.66098485, 'x3': 0.0, 'x4': 0.0, 'x5': 0.33901515, 'x6': 0.0}
组合后的投入和产出： ['19.36', '144.20', '1398.30', '699.37', '1706.27']
---------办事处6----------
目标函数的最小值z=0.96780303，此时目标函数的决策变量为: {'e': 0.96780303, 'x1': 0.0, 'x2': 0.66098485, 'x3': 0.0, 'x4': 0.0, 'x5': 0.33901515, 'x6': 0.0}
组合后的投入和产出： ['19.36', '144.20', '1398.30', '699.37', '1706.27']

从结果看，所有结果都一致。

旅行商问题

某邮递员每天从邮局（位于A小区）出发需要到7个不同位置的小区送信件，然后从这7个地方的邮筒收集信件回到邮局。通过长时间的观察记录，邮递员对7个区域之间骑车所需的平均时间进行整理，如下表所示。如何规划邮递员一天的路线才能使得路上花费的时间最少？

image-20210917172804723

这题用Python可以通过暴力的方式解决，思维也非常直接简单。

不过，今天我们主要的目标是深入讲解Python进行线性规划求解，所以通过这题演示一下线性规划中的二值变量的使用。

首先，我们使用基本的约束尝试一下：

  1. 不能是对角线
  2. 出发地唯一性
  3. 目标唯一性
from pulp import *
import numpy as np
import pandas as pd

distances = np.array([
    [0, 19, 20, 16, 25, 20, 19],
    [19, 0, 20, 17, 11, 12, 8],
    [20, 20, 0, 25, 17, 18, 20],
    [16, 17, 25, 0, 11, 8, 9],
    [25, 11, 17, 11, 0, 13, 14],
    [20, 12, 18, 8, 13, 0, 15],
    [19, 8, 20, 9, 14, 15, 0]
])
# 创建二值变量，i表示列（起点），j表示行（目标），xij表示当前路径是否被选中
x = np.array([
    [LpVariable(f"{ 
     i}{ 
     j}", cat=LpBinary)
     for i in "ABCDEFG"]
    for j in "ABCDEFG"
])
prob = LpProblem('目标函数和约束', sense=LpMinimize)
# 所需时间求最小
z = (distances*x).sum()
prob += z
# 不能是对角线，即自己到自己
for i in range(x.shape[0]):
    prob += x[i, i] == 0
# 出发地唯一性
for i in x.sum(axis=0):
    prob += i == 1
# 目标唯一性
for i in x.sum(axis=1):
    prob += i == 1

status = prob.solve()
print("z=", value(prob.objective))
print("路径选择：")
result = [[int(value(cell)) for cell in row]for row in x]
result = pd.DataFrame(result, index=list("ABCDEFG"), columns=list("ABCDEFG"))
print(result)

结果：

z= 88.0
路径选择：
   A  B  C  D  E  F  G
A  0  0  1  0  0  0  0
B  0  0  0  0  0  0  1
C  1  0  0  0  0  0  0
D  0  0  0  0  1  0  0
E  0  0  0  0  0  1  0
F  0  0  0  1  0  0  0
G  0  1  0  0  0  0  0

可以看到求解的结果路径中出现了，A<->C，B<->G，D->F->E->D这种循环路径的问题。

如何增加约束条件避免循环路径的问题呢？

假设 U i U_i Ui​表示第i个小区是第i个到达的，有理论研究（暂不懂其数学原理）表明当约束为：
{ u i − u j + n x i j ≤ n − 1 1 ≤ i ≠ j ≤ n − 1 \left\{\begin{array}{l} u_i-u_j+nx_{ij} \leq n-1 \\ 1 \leq i \neq j \leq n-1 \end{array}\right. {
ui​−uj​+nxij​≤n−11≤i​=j≤n−1​
其中n表示目的地的个数，即7个， x i j x_{ij} xij​表示指定路径是否被选择，有0或1两种可能。

可以避免上述循环路径的问题，那么我们继续在我们的程序中增加上述约束测试一下：

paths = [LpVariable(f"x{ 
     _}", lowBound=0, upBound=6,
                    cat=LpInteger) for _ in range(7)]
prob += paths[0] == 1
for j in range(1, 7):
    for i in range(1, 7):
        if i == j:
            continue
        prob += paths[i]-paths[j]+7*x[j, i] <= 6
status = prob.solve()
print("z=", value(prob.objective))

print("到达顺序：", [int(value(path)) for path in paths])

print("路径选择：")
result = [[int(value(cell)) for cell in row]for row in x]
result = pd.DataFrame(result, index=list("ABCDEFG"), columns=list("ABCDEFG"))
print(result)

结果：

z= 93.0
到达顺序： [1, 2, 0, 5, 1, 6, 3]
路径选择：
   A  B  C  D  E  F  G
A  0  0  0  0  0  1  0
B  0  0  0  0  1  0  0
C  1  0  0  0  0  0  0
D  0  0  0  0  0  0  1
E  0  0  1  0  0  0  0
F  0  0  0  1  0  0  0
G  0  1  0  0  0  0  0

从结果可以看到，最少花费时间为93，路径为A->C->E->B->G->D->F->A。

完整代码为：

from pulp import *
import numpy as np
import pandas as pd

distances = np.array([
    [0, 19, 20, 16, 25, 20, 19],
    [19, 0, 20, 17, 11, 12, 8],
    [20, 20, 0, 25, 17, 18, 20],
    [16, 17, 25, 0, 11, 8, 9],
    [25, 11, 17, 11, 0, 13, 14],
    [20, 12, 18, 8, 13, 0, 15],
    [19, 8, 20, 9, 14, 15, 0]
])
# i表示列（起点），j表示行（目标）
x = np.array([
    [LpVariable(f"{ 
     i}{ 
     j}", cat=LpBinary)
     for i in "ABCDEFG"]
    for j in "ABCDEFG"
])
prob = LpProblem('目标函数和约束', sense=LpMinimize)
# 所需时间求最小
z = (distances*x).sum()
prob += z
# 不能是对角线，即自己到自己
for i in range(x.shape[0]):
    prob += x[i, i] == 0
# 出发地唯一性
for i in x.sum(axis=0):
    prob += i == 1
# 目标唯一性
for i in x.sum(axis=1):
    prob += i == 1

paths = [LpVariable(f"x{ 
     _}", lowBound=0, upBound=6,
                    cat=LpInteger) for _ in range(7)]
prob += paths[0] == 1
for j in range(1, 7):
    for i in range(1, 7):
        if i == j:
            continue
        prob += paths[i]-paths[j]+7*x[j, i] <= 6
status = prob.solve()
print("z=", value(prob.objective))

print("到达顺序：", [int(value(path)) for path in paths])

print("路径选择：")
result = [[int(value(cell)) for cell in row]for row in x]
result = pd.DataFrame(result, index=list("ABCDEFG"), columns=list("ABCDEFG"))
print(result)

那么线性规划求解出来的结果是否正确呢？我们可以通过通过程序暴力遍历，获取最准确的最小值对比一下：

import itertools
import pandas as pd
import numpy as np

distances = np.array([
    [0, 19, 20, 16, 25, 20, 19],
    [19, 0, 20, 17, 11, 12, 8],
    [20, 20, 0, 25, 17, 18, 20],
    [16, 17, 25, 0, 11, 8, 9],
    [25, 11, 17, 11, 0, 13, 14],
    [20, 12, 18, 8, 13, 0, 15],
    [19, 8, 20, 9, 14, 15, 0]
])

distances = pd.DataFrame(distances, index=list(
    "ABCDEFG"), columns=list("ABCDEFG"))

min_time = distances.sum().sum()
min_paths = []
for path in itertools.permutations("BCDEFG", 6):
    path = ("A",)+path+("A",)
    distance_cur = 0
    for src, dest in zip(path[:-1], path[1:]):
        distance_cur += distances.at[dest, src]
    if min_time > distance_cur:
        min_time = distance_cur
        min_paths.clear()
    if min_time >= distance_cur:
        min_paths.append(path)
print("最短路径有：")
for min_path in min_paths:
    print("->".join(min_path))
print("最短耗时为：", min_time)

最短路径有：
A->C->E->B->G->D->F->A
A->F->D->G->B->E->C->A
最短耗时为： 93

从结果看到线性规则正确的求解出了其中一个正确的答案。

彩蛋：特殊技巧使Z3求解器可以求出多解

感谢还在认真看文的你，赠送一个技巧，让Z3约束求解器也可以获取所有的可行解。

前文在介绍z3求解器（详见文首地址）时，提到z3的一个缺点是只能找到一个可行解，无法找出所有的可行解。但其实只要我们采用一点小技巧，即可实现让z3找出所有的可行解。思路：每找出一个解就将其添加到约束中，不允许结果正好满足这个解，直到再也找不出可行解，就找到了所有的可行解。

以前面的八皇后问题为例进行演示，首先我们先看看普通的回溯算法求解的结果：

result = [0] * 8  # 角标代表皇后所在的行数，值存储皇后所在的列


def print_eight_queen(result):
    for column in result:
        for i in range(8):
            if i == column:
                print(end="Q ")
            else:
                print(end="* ")
        print()


n = 1  # 第n个满足条件的情况


def cal8queen(row: int = 0):
    global n
    if row == 8:
        print(f"----------{ 
     n}----------")
        print_eight_queen(result)
        n += 1
    for column in range(8):
        # 只有在第row行column列放置不与之前已经放置的冲突时才进行放置
        if is_ok(row, column):
            result[row] = column
            cal8queen(row + 1)


def is_ok(row, column):
    """ 校验在第row行将皇后放置在第column列是否满足条件 通过校验的条件是第row行前面所有行，都没有与当前列相同的列数，也不在当前列的对角线上 """
    leftup, rightup = column - 1, column + 1
    for i in range(row - 1, -1, -1):
        if column == result[i] or leftup == result[i] or rightup == result[i]:
            return False
        leftup -= 1
        rightup += 1
    return True


cal8queen(0)

截取首尾的结果打印如下：

----------1----------
Q  *  *  *  *  *  *  *  
*  *  *  *  Q  *  *  *  
*  *  *  *  *  *  *  Q  
*  *  *  *  *  Q  *  *  
*  *  Q  *  *  *  *  *  
*  *  *  *  *  *  Q  *  
*  Q  *  *  *  *  *  *  
*  *  *  Q  *  *  *  *  
...
----------92----------
*  *  *  *  *  *  *  Q  
*  *  *  Q  *  *  *  *  
Q  *  *  *  *  *  *  *  
*  *  Q  *  *  *  *  *  
*  *  *  *  *  Q  *  *  
*  Q  *  *  *  *  *  *  
*  *  *  *  *  *  Q  *  
*  *  *  *  Q  *  *  *

然后我们让z3约束求解器帮我们找出这92个解：

from z3 import *

# 每个皇后必须在不同的行中，记录每行对应的皇后对应的列位置
Q = [Int(f'Q_{ 
     i}') for i in range(8)]

# 每个皇后在列 0,1,2,...,7
val_c = [And(0 <= Q[i], Q[i] <= 7) for i in range(8)]

# 每列最多一个皇后
col_c = [Distinct(Q)]

# 对角线约束
diag_c = [If(i == j,
             True,
             And(Q[i] - Q[j] != i - j, Q[i] - Q[j] != j - i))
          for i in range(8) for j in range(i)]


def print_eight_queen(result):
    for column in result:
        for i in range(8):
            if i == column:
                print(end="Q ")
            else:
                print(end="* ")
        print()


s = Solver()
s.add(val_c + col_c + diag_c)
n = 1  # 第n个满足条件的情况
while s.check() == sat:
    result = s.model()
    result = [result[Q[i]].as_long() for i in range(8)]
    print(f"----------{ 
     n}----------")
    print_eight_queen(result)
    s.add(Not(And([qx == p for qx, p in zip(Q, result)])))
    n += 1

----------1----------
*  *  *  Q  *  *  *  *  
*  Q  *  *  *  *  *  *  
*  *  *  *  *  *  *  Q  
*  *  *  *  *  Q  *  *  
Q  *  *  *  *  *  *  *  
*  *  Q  *  *  *  *  *  
*  *  *  *  Q  *  *  *  
*  *  *  *  *  *  Q  *  
...
----------92----------
*  *  *  *  *  *  Q  *  
*  Q  *  *  *  *  *  *  
*  *  *  *  *  Q  *  *  
*  *  Q  *  *  *  *  *  
Q  *  *  *  *  *  *  *  
*  *  *  Q  *  *  *  *  
*  *  *  *  *  *  *  Q  
*  *  *  *  Q  *  *  *

可以看到也是92的解，虽然顺序不一样，但这92个解都找出来了。

版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 举报，一经查实，本站将立刻删除。

发布者：全栈程序员-用户IM，转载请注明出处：https://javaforall.cn/194021.html原文链接：https://javaforall.cn

【正版授权，激活自己账号】： Jetbrains全家桶Ide使用，1年售后保障，每天仅需1毛

【官方授权 正版激活】： 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

赞 (0)
全栈程序员-用户IM全栈程序员-用户IM
0 0


相关推荐

  • 硬核！SpringBoot连接MySQL数据库，十分钟啃透「建议收藏」

    硬核！SpringBoot连接MySQL数据库，十分钟啃透「建议收藏」

    硬核！SpringBoot连接MySQL数据库，十分钟啃透「建议收藏」目录1.pom.xml添加依赖2.application.properties添加数据库配置3.添加实体类4.添加Dao5.添加Controller6.新建数据库7.测试上篇只需两步！Eclipse+Maven快速构建第一个SpringBoot项目已经构建了一个SpringBoot项目，本文在此基础上进行连接MySQL数据库的操作。1.pom.xml添加依赖<dependency><groupId&.

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年6月25日
  • 图像检索评价指标：PR曲线的计算与绘制

    图像检索评价指标：PR曲线的计算与绘制

    图像检索评价指标：PR曲线的计算与绘制#@filename:test2.py#@brief:如何绘制PR曲线#@author:liupc#@date:2021/8/2importnumpyasnpfromtqdmimporttqdmimportmatplotlib.pyplotasplt#计算汉明距离。有几位不同，距离就为几。defCalcHammingDist(B1,B2):q=B2.shape[1]distH=.

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年6月23日
  • 更新pip版本(anaconda更新pip)

    更新pip版本(anaconda更新pip)

    @更新pip版本#window：正常windows下更新pip版本的指令无法更新时，可以尝试如下指令python-mpipinstall–upgradepip#Ubuntu更新pip版本指令：sudo-Hpython-mpipinstall–upgradepip…

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年4月11日
  • 在树莓派4B上使用YOLO v3 Tiny进行实时目标检测

    在树莓派4B上使用YOLO v3 Tiny进行实时目标检测

    在树莓派4B上使用YOLO v3 Tiny进行实时目标检测首先尝试使用yolo官网yolo给的教程，在树莓派上测试，但是在运行时出现段错误，尝试很多方法无法解决。在国外的网站找到darknet-nnpack这个东西，可以完美的在树莓派上运行。参考的文章darknet-nnpackhttps://github.com/zxzhaixiang/darknet-nnpackFunofDIYhttp://funofdiy.blogspot.c…

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年6月9日
  • 一名正在学习Python的新手。「建议收藏」

    一名正在学习Python的新手。「建议收藏」

    一名正在学习Python的新手。「建议收藏」1.写代码，有如下列表，按照要求实现每一个功能计算列表的长度并输出li=["alex","WuSir","ritian",&q

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年7月5日
  • 关于网络制式（例如：TD-SCDMA /GSM）

    关于网络制式（例如：TD-SCDMA /GSM）

    关于网络制式（例如：TD-SCDMA /GSM）在中国，TD-SCDMA指的是中国移动的3G，WCDMA指的是中国联通的3G，CDMA2000指的是中国电信的3G，GSM是2G网络（中国联通和中国移动都有），简言之，CDMA2000是CDMA演进的3G！CDMA(CodeDivisionMultipleAccess)又称码分多址，是在无线通讯上使用的技术。在国内，是中国电信的手机网络在用的方式。三大运营商网络制式…

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年10月3日

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号