数学建模算法(1)—线性规划模型

这里写自定义目录标题

什么是规划问题

在人们的生产实践中,经常会遇到如何利用现有资源来安排生产,以取得最大经济效益的问题。此类问题构成了运筹学的一个重要分支—数学规划,而线性规划(Linear Programming 简记LP)则是数学规划的一个重要分支。

例1: 某机床厂生产甲、乙两种机床,每台销售后的利润分别为4000元与3000 元。生产甲机床需用 A、B机器加工,加工时间分别为每台 2 小时和 1 小时;生产乙机床需用 A、B、C三种机器加工,加工时间为每台各一小时。若每天可用于加工的机器时数分别为A 机器10 小时、B 机器8 小时和C 机器7 小时,问该厂应生产甲、乙机床各几台,才能使总利润最大?

上述问题的数学模型:设该厂生产x1x_1x1​台甲机床和x2x_2x2​台乙机床时总利润最大,则x1,x2x_1,x_2x1​,x2​应满足
maxz=4x1+3x2s.t.{2x1+x210x1+x28x27x1,x20 \max z = 4x_1 + 3x_2 \\ \textbf{s.t.} \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. maxz=4x1​+3x2​s.t.⎩⎪⎪⎨⎪⎪⎧​2x1​+x2​≤10x1​+x2​≤8x2​≤7x1​,x2​≥0​

这里变量x1,x2x_1,x_2x1​,x2​称之为决策变量,maxz=4x1+3x2\max z = 4x_1 + 3x_2maxz=4x1​+3x2​ 被称为问题的目标函数,下面的几个不等式是问题的约束条件,记为s.t.(即subject to)。

简单来说,线性规划问题是在一组线性约束条件的限制下,求一线性目标函数最大或最小的问题,无论是约束条件还是目标函数出现非线性项,那么就问题就变成了非线性规划。

规划问题的分类

类别 定义
线性规划 在一组线性约束条件的限制下,求一线性目标函数最大或最小的问题
非线性规划 无论是约束条件还是目标函数出现非线性项1,那么规划问题就变成了非线性规划
整数规划 当约束条件加强,要求所有的自变量必须是整数时,成为整数规划(特别地,自变量只能为0或1时称为0-1规划)
多目标规划 在一组约束条件的限制下,求多个目标函数最大或最小的问题

本文仅关心线性规划的数学形式及其求解方法。

线性规划及其Python解法

线性规划是所有规划问题的基础,因此我们首先讨论线性规划问题及其Python实现。对于一个规划类问题来说,我们通常用如下的步骤建立线性规划:

1. 选择适当的决策变量
在解决实际问题时,把问题归结成一个线性规划数学模型是很重要的一步,但往往也是困难的一步,模型建立得是否恰当,直接影响到求解。而选适当的决策变量,是我们建立有效模型的关键之一。

2. 将求解目标简化为求一个目标函数的最大/最小值
能把要求解的问题简化为一个最值问题是能否使用线性规划模型的关键,如果这一点不能达到,之后的工作都有没有意义的。

3. 根据实际要求写出约束条件(正负性,资源约束等)
线性规划的约束条件针对不同的问题有不同的形式,总结来说有以下三种:

  • 等式约束
  • 不等式约束
  • 正负性约束

单纯形法是求解线性规划问题的最常用、最有效的算法之一。这里我们就不介绍单纯形法的数学原理,有兴趣的同学可以参看其它线性规划书籍。下面我们直接介绍线性规划的Python求解方法。以下面的规划问题为例:
maxz=2x1+3x25x3s.t.{x1+x2+x3=72x15x2+x310x1+3x2+x312x1,x2,x30\max z=2 x_{1}+3 x_{2}-5 x_{3} \\ \textbf{s.t.} \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.maxz=2x1​+3x2​−5x3​s.t.⎩⎪⎪⎨⎪⎪⎧​x1​+x2​+x3​=72x1​−5x2​+x3​≥10x1​+3x2​+x3​≤12x1​,x2​,x3​≥0​

实际的线性规划问题千变万化,约束条件大于、小于、等于的情况均有发生,为了简化程序的输入,我们首先定义线性规划的标准型
minxcTx s.t. {Axb Aeq x=beqlbxub \begin{array}{l} {\displaystyle \min_x c^{T} x} \\ {\textbf { s.t. }\left\{\begin{array}{l}{A x \leq b} \\ {\text { Aeq } \cdot x=b e q} \\ {l b \leq x \leq u b}\end{array}\right.}\end{array} xmin​cTx s.t. ⎩⎨⎧​Ax≤b Aeq ⋅x=beqlb≤x≤ub​​

其中,标准型要求目标函数都要化成求最小值的形式,对于需要求最大值的问题,直接取负即可。此外,标准型还要求不等式约束都用"\leq≤"给出。

本例对应的标准型中各个参数为
c=[2,3,5]T,A=[251131],b=[10,12]T,Aeq=[1,1,1]beq=[7] c= [-2, -3 ,5]^T, \quad A =\begin{bmatrix} -2 & 5 & -1\\1 & 3 & 1 \end{bmatrix} , \quad b= [-10,12]^T, \quad Aeq = [1,1,1] \quad beq = [7] c=[−2,−3,5]T,A=[−21​53​−11​],b=[−10,12]T,Aeq=[1,1,1]beq=[7]

这里我们使用scipy中的op库中的linprog方法进行线性规划的计算,代码如下

from scipy import optimize as op
import numpy as np
c=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=(0,7)								  #定义变量x_1的范围
x2=(0,7)								  #定义变量x_2的范围
x3=(0,7)								  #定义变量x_3的范围
res=op.linprog(-c,A_ub,B_ub,A_eq,B_eq,bounds=(x1,x2,x3))   #调用函数进行求解

非线性规划及其Python解法

对于如下的非线性规划问题:

minf(x)=x12+x22+x32+8 s.t. {x12x2+x320x1+x22+x3320x1x22+2=0x2+2x32=3x1,x2,x30{\min f(x)=x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+8} \\ \textbf { s.t. } \left\{ \begin{array}{c} \displaystyle { x_{1}^{2}-x_{2}+x_{3}^{2} \geq 0} \\ {x_{1}+x_{2}^{2}+x_{3}^{3} \leq 20} \\ {-x_{1}-x_{2}^{2}+2=0} \\ {x_{2}+2 x_{3}^{2}=3} \\ {x_{1}, x_{2}, x_{3} \geq 0}\end{array} \right.minf(x)=x12​+x22​+x32​+8 s.t. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​x12​−x2​+x32​≥0x1​+x22​+x33​≤20−x1​−x22​+2=0x2​+2x32​=3x1​,x2​,x3​≥0​

我们可以使用scipy.optimize中的minmize方法进行求解,注意这里的标准型与上面的线性规划略有不同,对于非线性不等式约束采用的是左侧大于等于右侧。

minf(x){Aeqx=BeqC(x)0Ceq(x)=0 \min f(x)\\ \left\{ \begin{array}{l} {\text {Aeq} \cdot x=\operatorname{Beq}} \\ {C(x) \geq 0} \\ {\operatorname{Ceq}(x)=0} \end{array} \right. minf(x)⎩⎨⎧​Aeq⋅x=BeqC(x)≥0Ceq(x)=0​

对应于以上的问题,求解代码如下:

from scipy import  optimize as opt
import numpy as np
from scipy.optimize import minimize
# 目标函数
def objective(x):
    return x[0] ** 2 + x[1]**2 + x[2]**2 +8

# 约束条件
def constraint1(x):
    return x[0] ** 2 - x[1] + x[2]**2  # 不等约束

def constraint2(x):
    return -(x[0] + x[1]**2 + x[2]**2-20)  # 不等约束,注意minimize方法的标准型为左边大于等于右边

def constraint3(x):
    return -x[0] - x[1]**2 + 2        # 等式约束

def constraint4(x):
    return x[1] + 2*x[2]**2 -3           # 等式约束

# 边界约束
b = (0.0, None)
bnds = (b, b ,b) 

con1 = {'type': 'ineq', 'fun': constraint1}
con2 = {'type': 'ineq', 'fun': constraint2}
con3 = {'type': 'eq', 'fun': constraint3}
con4 = {'type': 'eq', 'fun': constraint4}
cons = ([con1, con2, con3,con4])  # 4个约束条件

# 计算
x0=np.array([0, 0, 0])
solution = minimize(objective, x0, method='SLSQP', \
                    bounds=bnds, constraints=cons)
x = solution.x

print('目标值: ' + str(objective(x)))
print('答案为')
print('x1 = ' + str(x[0]))
print('x2 = ' + str(x[1]))
print('x3 = ' + str(x[2]))

  1. 非线性项即为不满足叠加原理的项,例如:x2,1x,log2(x)x^2,\dfrac{1}{x},\log_2(x)x2,x1​,log2​(x),他们均为非线性项。 ↩︎

上一篇:一种在程序中求两直线交点的简单数学方法


下一篇:带默认形参值的函数 C++练习