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
blank

相关推荐

  • Django（26）HttpResponse对象和JsonResponse对象「建议收藏」

    Django（26）HttpResponse对象和JsonResponse对象「建议收藏」

    Django（26）HttpResponse对象和JsonResponse对象「建议收藏」HttpResponse对象Django服务器接收到客户端发送过来的请求后，会将提交上来的这些数据封装成一个HttpRequest对象传给视图函数。那么视图函数在处理完相关的逻辑后，也需要返回一个响

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年7月30日
  • Kettle Spoon入门教程「建议收藏」

    Kettle Spoon入门教程「建议收藏」

    Kettle Spoon入门教程「建议收藏」Kettle是一款国外开源的ETL工具，纯java编写，可以在Window、Linux、Unix上运行，数据抽取高效稳定。其中，Spoon是Kettle中的一个组件，其他组件有PAN，CHEF，Encr和KITCHEN等。Spoon通过图形化的页面，方便直观的让你完成数据转换的操作。1、安装配置点我下载最新版本是7.1，小编这里使用的是6.1。下载完毕后，解压即可。因为是Java编写的，需要jdk环

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年5月23日
  • Javadeque_deque接口

    Javadeque_deque接口

    Javadeque_deque接口Queue也是Java集合框架中定义的一种接口，直接继承自Collection接口。除了基本的Collection接口规定测操作外，Queue接口还定义一组针对队列的特殊操作。通常来说，Queue是按照先进先出(FIFO)的方式来管理其中的元素的，但是优先队列是一个例外。Deque接口继承自Queue接口，但Deque支持同时从两端添加或移除元素，因此又被成为双端队列。鉴…

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年9月20日
  • Java必备常见单词

    Java必备常见单词

    Java必备常见单词资源共享学习交流群号:769674658（快满）qq交流二群（296389054）（一）Java基础 public公有的 private私有的 protected保护的 …

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年10月9日
  • textview可复制_android长按点击

    textview可复制_android长按点击

    textview可复制_android长按点击有这么一个需求,用户在浏览文本信息时希望长按信息就能弹出复制的选项方便保存或者在别的页面使用这些信息.类似的,就像长按WebView或者EditText的内容就自动弹出复制选项.这里面主要是2个特点:1,用户只能浏览文本信息而不能编辑这些文本信息;2,用户对着文本信息长时间点按可以弹出”复制”选项实现复制;网上有好多种方法可实现,也比较零散,此处做个小结,希望有所帮助.1,通过继承EditT…

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年9月29日
  • MBUS系列产品特点（科慧铭远）[通俗易懂]

    MBUS系列产品特点（科慧铭远）[通俗易懂]

    MBUS系列产品特点（科慧铭远）[通俗易懂]     北京科慧铭远自控技术有限公司联合国际标准化组织、计量中心、热力集团、清华大学检测与电子技术研究所，成立国内首家M-BUS通信设备检测中心，对于M-BUS主站、从站的通信设备全方位的检测其是否符合国际和国内标准，并予以认证。北京科慧铭远自控技术有限公司有着在M-BUS领域最全面的的研发和生产能力，获得国际标准化组织的认可，产品在欧洲、亚洲、中国获得全面应用。其生产的MBUS设备的主要…

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

发表回复

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

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