在有限资源和竞争约束（人、财、物、时间）下取得最佳经济效益--线性规划（下）

第四部分 线性规划建模工具

$scipy.optimize.linprog$ 参数说明

$scipy.optimize.linprog$ $( c ,\ A_{ub}=None ,\ b_{ub}=None ,\ A_{eq}=None ,\ b_{eq}=None ,\ bounds=None ,\ method='simplex' ,\ callback=None ,\ options=None )$

  • $c$ ：求最大值的函数的系数数组
  • $A_{ub}$ ：不等式未知量的系数矩阵
  • $B_{ub}$ ：不等式的右边
  • $A_{eq}$ ：其中等式的未知量系数矩阵
  • $B_{eq}$ ：等式右边
  • $bounds$ ：每个未知量的范围

  这个模块只能找到线性规划的最小值，因此，如果约束条件的不等式（变量除外）中含有“大于等于”的需要修改为“小于等于”。

$PuLp$ 参数说明

该函数库只能用于线性模型

安装PuLp
pip install pulp

导入PuLp
from pulp import *

定义线性规划问题
PB = LpProblem(problem name, sense)

  • 构造函数，用来构造一个LP问题实例，其中:
    • name 指定问题名（输出信息用）
    • sense 值是 LpMinimize 或 LpMaximize 中的一个，用来指定目标函数是求最大值还是最小值

定义决策变量
DV = LPVariable(decision variable name,lowbound,upbound,category)

  • decision variable name 指定变量名
  • lowBound和upBound是下界和上界， 默认分别是负无穷到正无穷
  • cat用来指定变量是离散(LpInteger,LpBinary)还是连续(LpContinuous)

添加目标函数
PB += linear objective in equantion from objective name

添加约束条件
PB += linear objective in equantion from constraint name

写入LP文件
PB.writeLP(filename)

模型求解
PB.slove())

结果显示
check status : LpStatus[status]

线性规划原理之图解法

问题

  $M1$和$M2$两种原料用于生产内外墙涂料，$M1$日最大可用量24吨，$M2$日最大可用量为6吨，外墙涂料每吨需要6吨$M1$，1吨$M2$，内墙涂料每吨需要4吨$M1$，2吨$M2$，外墙涂料每吨利润5个单位，内墙涂料每吨利润4个单位。且市场需求调查数据得出，内墙日需求量不超过外墙的日需求量+1吨，内墙最大日需求量为2吨。

建模

设外墙生产$x_1$吨,内墙生产$x_2$吨,设利润为$z$，要得到$z$的最大化，也就是最优解

$$ max\ z=5x_1+4x_2 \\ $$

$$ s.t.\ \ \left\{ \begin{aligned} 6x_1+4x_2 \leq 24 \\ x_1+2x_2 \leq \ \ 6 \\ -x_1+x_2 \leq \ \ 1 \\ x_2 \leq \ \ 2 \\ x_1,x_2 \geq \ \ 0 \\ \end{aligned} \right. $$

几何图解

  $Z$ 函数表示一条截距为$\dfrac{z}{4}$斜率为$-\dfrac{5}{4}$的线性直线，要求$Z$最大化的最优解，就是在所有的可行区域内找到可以满足$Z$曲线截距最大的点。实际上也就是在多面形的各个顶点逐一求解。 $$ x_2=-\dfrac{5}{4}x_1+\dfrac{z}{4} $$   从下图看出顶点为 $6x_1+4x_2=24$ 和 $x_1+2x_2=6$ 的交点，计算可得为[3, 1.5]

In [301]:
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.font_manager import *
myfont = FontProperties(fname='/jupyterfile/simsun.ttf')
figure, ax=plt.subplots(figsize=(10, 6))
ax.set_xlim(0,6)
ax.set_ylim(0,6)
#设x1为自变量x，x2为因变量y。
x=np.linspace(0,6,700)
y1=6-x*6/4
y2=3-x*1/2
y3=1+x
y4=[2*i/i for i in range(1,701)]
y0=[0*i/i for i in range(1,701)]
y6=[6*i/i for i in range(1,701)]
yz1=-x*5/4+3.25
yz2=-x*5/4+4.5
yz3=-x*5/4+5.25
plt.plot(x, y1 , "y", lw=1 )
plt.plot(x, y2 , "k", lw=1 )   
plt.plot(x, y3 , "g", lw=1 )   
plt.plot(x, y4 , "b", lw=1 )   
plt.plot(x, y0 , "r", lw=1 )
plt.plot(x, y6 , "r", lw=1 )
plt.plot(x, yz1 , "r--", lw=2 )  
plt.plot(x, yz2 , "r--", lw=2 )  
plt.plot(x, yz3 , "r--", lw=2 )
ax.fill_between(x,y1,y0,where=y1>y0,color='g',alpha=0.1)
ax.fill_between(x,y2,y0,where=y1>y0,color='b',alpha=0.1)
ax.fill_between(x,y3,y0,where=y1>y0,color='k',alpha=0.1)
ax.fill_between(x,y4,y0,where=y1>y0,color='y',alpha=0.1)
ax.fill_between(x,y1,y6,where=y1<y6,color='w',alpha=1)
ax.fill_between(x,y2,y6,where=y1<y6,color='w',alpha=1)
ax.fill_between(x,y3,y6,where=y1<y6,color='w',alpha=1)
ax.fill_between(x,y4,y6,where=y1<y6,color='w',alpha=1)
plt.text(1, 0.75, '所有的可行解', fontproperties=myfont, fontsize=24, color='b') 
plt.text(-1, 5.2, '截距(max z)', fontproperties=myfont, fontsize=16, color='b') 
plt.text(3, 3, 'www.jasper.wang', fontproperties=myfont, fontsize=24, color='k') 
plt.show()

工具验证

$SciPy.optimize.linprog$

In [302]:
from scipy import optimize as op
import numpy as np
c=np.array([5,4])
A_ub=np.array([[6,4], [1,2], [-1, 1], [0, 1]])
B_ub=np.array([24, 6, 1, 2])
x1=(0,3) # x2最大为2，x1最大不会超过x21个单位，所以2+1=3，x1最大不超过3
x2=(0,2)
res=op.linprog(-c, A_ub, B_ub, bounds=(x1, x2))
print('最大利润为：',-res.fun,'个单位')
print('生产M1',res.x[0],'吨', '\n生产M2',res.x[1],'吨')

最大利润为： 21.0 个单位
生产M1 3.0 吨 
生产M2 1.5 吨

$pulp$

In [303]:
from pulp import *
prob = LpProblem("problem1",LpMaximize)
x1 = LpVariable("x1",0,None,LpContinuous)
x2 = LpVariable("x2",0,None,LpContinuous)
prob += 5 * x1 + 4 * x2
prob += 6 * x1 + 4 * x2 <= 24
prob += x1 + 2 * x2 <= 6
prob += -1 * x1 + x2 <= 1
prob += 0 * x1 + x2 <= 2
prob.writeLP("problem1.lp")
prob.solve()
print("最大利润为",value(prob.objective),"个单位")
for v in prob.variables():
    print("生产",v.name,":",v.varValue,"吨")

最大利润为 21.0 个单位
生产 x1 : 3.0 吨
生产 x2 : 1.5 吨

第五部分 单纯形算法推导

主要思想

  可以理解成”不等式的高斯消元”
  单纯形通过不断改变基本变量与非基本变量的集合，不断重写松弛形，直到变成最优解显而易见的形式
  基变量 是运筹学的一个术语。在线性规划问题约束条件方程组中，系数矩阵中的基向量对应的变量称为基变量
  非基变量 是运筹学的一个术语。在线性规划中除基变量以外的变量称为非基变量
  松弛变量 为了把一般线性规划模型变为标准型，需要把不等式约束条件变为等式约束条件，于是引入松弛变量和剩余变量

化标准型

线性规划表示方法标准型

$$ 最大化： \sum_{i=1}^n c_ix_i $$

$$ 满足约束：\sum_{j=1}^n a_{ij} * x_j \leq b_i \ \ \ i=1, 2, ..., m \ \ \ \ 参考矩阵与向量相乘 $$

$$ x_i \geq 0 \ \ \ i=1, 2, ..., m $$

化标准型方法

  • 若目标是要最小化，取负号 $"-"$ 化为最大化
  • 若有变量没有非负约束，拆成 $x_i=x_i^{'} − x_i^{''}\ \ 且\ \ \ x_i^{'},\ \ x_i^{''} \geq 0$
  • 若有等式，就变成 $\geq$ 和 $\leq$ 两个约束
  • 若有 $\geq$ 约束，就取反

化单纯形法标准型

  • 所有的变量都是非负数（$\geq 0$）
  • 所有的约束都是等式( $=$ )(非负限制除外)，且具有非负的右端项
  • 根据需要而加入的变量称为松弛变量(非负)

接前一示例的单纯形法求解过程

$$ max \ \ z=5x_1+4x_2+0s_1+0s_2+0s_3+0s_4 $$ $$ s.t.\ \ \left\{ \begin{aligned} 6x_1+4x_2+s_1=24 \\ x_1+2x_2+s_2=\ \ 6 \\ -x_1+x_2+s_3=\ \ 1 \\ x_2+s_4=\ \ 2 \\ \end{aligned} \right. $$
$$s_1, s_2, s_3, s_4\ 表示松弛变量(非负)，根据 \geq 、\leq 来决定他们的正负号$$

对于标准化形式

  • 设列举的约束方程个数为 $m$
  • 设参数个数为 $n$
  • 当 $m = n$ 时，方程组有唯一解
  • 当 $m < n$ 时，方程组有无穷个解（即解是一个区域）
  • 当 $m > n$ 时，方程组无解

确定基变量、非基变量

  • 令 $(n-m)$ 个变量为 $0$ 的这些变量为非基变量，令其余的 $m$ 个变量为基变量（$n=6, \ \ m=4$）
  • 从原点，也就是从非松弛变量全为0的时候开始
    $$ \left\{ \begin{array}{aligned} s_1=24-6x_1-4x_2 \\ s_2=6-x_1-2x_2 \\ s_3=1+x_1-x_2 \\ s_4=2-x_2 \\ \end{array} \right. $$
  • 右边的变量为非基变量，左边为基变量
  • 因为非基变量为0，所以左边的基变量的解值就是右边的常系数，也就是截距
基变量 z系数 x1系数 x2系数 s1系数 s2系数 s3系数 s4系数 解值
$min\ z$ 1 -5 -4 0 0 0 0 0
$s_1$ 0 6 4 1 0 0 0 24
$s_2$ 0 1 2 0 1 0 0 6
$s_3$ 0 -1 1 0 0 1 0 1
$s_4$ 0 0 1 0 0 0 1 2

换基，选择进基于离基

  • 选择 $x_1$ 为进基变量，因为其在目标方程中的求max时的系数最大，为 5。（或求min时的系数最小）
  • 选择 $s_1$ 为离基变量，因$S_1$行的解（24）与$x_1$系数（-6）之比（$\dfrac{24}{6}$）为最小非负数
  • 进行矩阵运算，使得矩阵的第一行中，代表$x_1,s_2,s_3,s_4$的系数为0，$s_1$不为0。
  • 对于迭代后的新枢纽行的计算规则为:
    • 进基变量所在列为枢纽列，离基变量所在行为枢纽行，行列交叉位置为枢纽元素(上表中的6)
    • 在基列中，以进基变量替换离基变量，$新的枢纽行=\dfrac{当前枢纽行}{枢纽元素}$
    • ($\dfrac{0}{6},\dfrac{6}{6},\dfrac{4}{6},\dfrac{1}{6},\dfrac{0}{6},\dfrac{0}{6},\dfrac{0}{6},\dfrac{24}{6}$) $\Longrightarrow$ ($0,\ 1,\ \dfrac{2}{3},\ \dfrac{1}{6},\ 0,\ 0,\ 0,\ 4$)$\ \ $备注：当$\dfrac{2}{3},\ \dfrac{1}{6}$移动到基列的等号右侧时变为负号
    • 其他所有行，包括Z行的计算规则为，$新的行 = 当前行 - 当前行枢纽列的系数 \times 新的枢纽行$
      • 对于 $max\ z$ 行，$当前行枢纽列的系数 \times 新的枢纽行$ 为 $-5 \times (0,\ 1,\ \dfrac{2}{3},\ \dfrac{1}{6},\ 0,\ 0,\ 0,\ 4)$ $\Longrightarrow$ $(0, -5, \dfrac{10}{3}, -\dfrac{5}{6}, 0, 0, 0, -20)$
      • 对于 $max\ z$ 行，新行值为$(1, -5, -4, 0, 0, 0, 0, 0) - (0, -5, \dfrac{10}{3}, -\dfrac{5}{6}, 0, 0, 0, -20)=(1, 0, -\dfrac{2}{3}, \dfrac{5}{6}, 0, 0, 0, 20)$
基变量 z系数 x1系数 x2系数 s1系数 s2系数 s3系数 s4系数 解值
$min\ z$ 1 0 $-\dfrac{2}{3}$ $\dfrac{5}{6}$ 0 0 0 20
$x_1$ 0 1 $-\dfrac{2}{3}$ $-\dfrac{1}{6}$ 0 0 0 4
$s_2$ 0 0 $\dfrac{4}{3}$ $-\dfrac{1}{6}$ 1 0 0 2
$s_3$ 0 0 $\dfrac{5}{3}$ $\dfrac{1}{6}$ 0 1 0 5
$s_4$ 0 0 1 0 0 0 1 2
  • 上表结果，与方程组移项推导结果一致：

$$ \left\{ \begin{array}{aligned} x_1=4-\dfrac{1}{6}s_1-\dfrac{2}{3}x_2 \\ s_2=2+\dfrac{1}{6}s_1-\dfrac{4}{3}x_2 \\ s_3=5-\dfrac{1}{6}s_1-\dfrac{5}{3}x_2 \\ s_4=2-x_2 \\ \end{array} \right. $$

  • 上述过程一直循环迭代 直到没有最优条件时，此时 $z$ 就是最优解，基变量的值就是 $z$ 最优时的变量的解。

爱与彼岸财经分析团队编制

 爱与彼岸财经分析