本专栏内容为：数学建模原理 记录学习数学建模

💓博主csdn个人主页：小小unicorn
⏩专栏分类：数学建模
🚚代码仓库：小小unicorn的代码仓库🚚
🌹🌹🌹关注我带你学习编程知识

从线性代数->线性规划

线性规划（回忆）

在高中阶段我们其实就学习过线性规划的知识，但是当时我们可能没有觉得它有多重要所以可能忘掉了。我们通过一个例子回忆一下：

例3.1 若𝑥，𝑦满足条件

则𝑧=3𝑥+𝑦的最小值为（）
A. 18 B. 10 C. 6 D. 4

我想这个问题大家不会感到陌生了。答案很简单，选择C。血脉觉醒的高中生们会将不等式组中的不等式两两配对成三个二元一次方程组求解，得到三个点，最后把三个点分别代入z=3x+y最后得到最小值。这是比较快的方法，在90%的情况下这种策略是很奏效的。而学霸们也会采取画图的策略，平移直线来求解。这个问题的可行域如下图所示：

这个题目的解其实很明显，就是把不等式组里面每个不等式在平面直角坐标系中表示出来，然后根据不等号的方向确定可行域。将目标函数进行移项以后转化为𝑦=−3𝑥+𝑧，通过平移直线的方式找到可行域内使目标函数截距最大的点，就是正确答案。大家回忆起来了吗？

注意：平移直线法是一种通用方法，解方程组得到的最优解很多情况下确实有效但不排除有些特殊情况下可行域是开放的，这个时候不一定存在最小值或最大值。

我为什么把线性规划作为第一个知识点呢？因为线性规划是真的非常重要也非常实用，也最贴近一个入门级小白认识数学建模所具备的数学基础。线性规划的实际应用有很多，比如说：我们可能见过这样一种问题，去运输一批货物的时候大车能运五箱，小车只能运三箱，但我们大车和小车数量都有限，怎么安排运输方案能够在车辆够用的情况下运费还能最小？

如果我们把大车数量记为𝑥，小车数量记为𝑦，那么除了𝑥和𝑦的范围，5𝑥+3𝑦也有自己的一个范围，算上运费作为优化目标，这不就构成了一个线性规划吗？背景熟悉吧，甚至于有一种小学应用题的恐惧感。

再打个比方，生产原料问题，生产产品A、B需要原材料甲、乙、丙；生产一吨A需要多少多少甲、多少多少乙、多少多少丙，这样就有了对AB的三个原料约束。再来一个利润最大作为目标，这也是线性规划。我猜很多读者看到这些例子可能就会暗想：“这就是数学建模？我怎么感觉梦回小学应用题？难不成我被骗了？”，是的，数学建模其实没有那么恐怖，小时候我们做的是应用题，到大了我们只不过是需要用更多知识和编程方法解背景更产业化更学术化的应用题。因为这是一门应用数学学科。

我们大概可以总结出，中学的线性规划通常就两个变量𝑥和𝑦，约束条件三个不等式，最后一个线性的形如
𝑧=𝑎𝑥+𝑏𝑦（这里𝑎和𝑏都是常数）的目标函数。这样的式子我们解方程可以解，画图也可以解，总能在两分钟之内算出正确结果。但在实际情况中，问题真的有这么容易吗？其实不然。同样是拿生产问题做文章，如果我们这里生产的原料不止甲乙丙三种呢（通常在有机化合物合成的时候原料可能有十几种甚至上百种），产出的产品也不止AB而是能够产出数十种化合物，还能简单地用高中的方法写吗？

所以我们说，中学的规划存在这样一些局限性：

• 决策变量（如果不好理解暂且称之为自变量吧）往往不止两三个；
• 当变量个数超过三个的时候还能在直角坐标系里面画图吗？不能了；
• 约束条件往往不止三个不等式，不等式可能比变量更多一些；
• 当变量较多的时候，还可能出现方程形式下的约束；
• 中学阶段我们只讨论了线性规划，但如果不等式或者目标函数非线性呢。

这么多情况不知读者朋友会不会被吓到。如果十几个甚至几十个不等式方程组成约束条件，那我草稿纸甚至不知道写不写的下，况且中学阶段没有接触过高维问题。于是，为了以更简单的形式描述更一般的线性规划，我们需要借助一样数学工具——线性代数。

线性规划基本形式

前面我们已经看到，中学线性规划的局限性在于难以描述多约束、多变量，但无论是目标函数还是方程还是不等式，我们都可以看成是一个系数向量与变量向量在做乘法（例如，2𝑎+3𝑏−𝑐 实际上可以看成向量 [2,3,−1]与向量 [𝑎,𝑏,𝑐] 做内积）。多个约束无非就是把多个向量拼接在一块做成了一个矩阵而已。

我们把所有的方程约束中系数做成系数矩阵𝐴_eq ，等号右边的常数作为列向量𝑏_eq；不等式约束中的系数矩阵𝐴和不等号右边的常数 𝑏 ，为了方便起见通常将不等式统一为小于等于；变量𝑥在向量 $l$_b 到 𝑢_𝑏 之间取值；目标函数的系数向量为 𝑐 ，那么线性规划的标准形式就如下所示：

使用Numpy进行矩阵运算

线性代数基本知识

线性代数是一门基础数学科目，基本上所有理工科学生大一的时候都得学线性代数。如果是数学系可能学的就是高等代数了。这门课主要是研究矩阵与向量的数学理论，也会探究线性方程组的解等相关问题。

相信大多数同学高中毕业是记得向量这个概念的，但中学阶段我们也仅仅是接触到了三维向量。事实上向量的维数可以是很多维，从代数的意义上你可以认为向量是一个集合，从几何的意义上你又可以认为向量是一个𝑛维 Euclid 空间中的一个点：


和二维、三维空间中的向量一样，高维空间中的向量同样可以进行加减运算、数量乘运算和数乘运算。但毕竟这是一门应用数学课程，我们不打算把太多精力放在任何一本线代课本里面都能找到的公式上，使matlab举例子恐怕会更加直观。从程序设计的角度来看，如果读者接触过C语言应该会了解数组的概念，而在C++语言中STL里面已经包含了vector类型。

数通过集合形成了向量，那向量集合以后又会变成什么呢？如果向量只是沿着同一个方向进行拼接，那么得到的只不过是一个更长一些的向量。但如果是在纵向上做拼接，那么我们或许可以把一个向量排成表格：

这个数表要求每个矩阵的维度相同，排成的这个表格就可以称作一个矩阵。那么矩阵作为向量的集合，自然也保留了向量的一些特性，包括行列相同的矩阵的加减法、数量乘。比较有趣的是矩阵的乘法，它把两个矩阵分别按行、列规约：



需要修改，改成“平行多面体的体积” 那我们能否类比向量的模，提出“矩阵的模”这样一个概念呢？在线性代数中确实存在这样一个类似的概念，这个概念叫行列式：

行列式虽然同样是排成了一个表，但注意，矩阵是一个表格，行列式是一个数，它的值是可以算出来的！有关行列式的计算方法有很多，但最经典最通用的方法是代数余子式展开法。代数余子式的本质就是递归式求解，将行列式A中下标为（i,j）的元素所在行和所在列全部去除以后求新的行列式，再乘上对应的符号。而对于行列式A，计算定义为：

注意：行列式除了按某一行展开也可以按某一列展开，这一展开行或展开列的选取是任意的，方便计算即可。另外，矩阵可以不要求行列数相等但行列式必须行列数相等！

而递归到最后，我们逃不开最低的二阶行列式求解。二阶行列式的定义为：
将𝑛阶行列式降低到𝑛−1阶，𝑛−1阶再降低为n-2阶，逐层展开到最后二阶，整个行列式求解就可以完成了。不过高阶行列式如果不是特殊行列式计算会有些复杂，我们可以将这个过程交给计算机程序来完成。

这样我们就可以在命令行输出矩阵A对应的行列式的值。有了行列式的概念，我们可以用它定义矩阵的逆矩阵计算方法。一个行列数均为n的矩阵的逆矩阵满足这样的定义：

Numpy的基本使用

使用numpy之前首先需要导入模块：

import numpy as np

接下来我们需要了解一下numpy的几种属性:

• ndim：维度
• shape：行数和列数
• size：元素个数

array=np.array([[1,2,3],[4,5,6]])
print('number of dim:',array.ndim)  # 维度
# number of dim: 2
print('shape :',array.shape)    # 行数和列数
# shape : (2, 3)
print('size:',array.size)   # 元素个数
# size: 6

numpy的创建array

• array：创建数组
• dtype：指定数据类型
• zeros：创建数据全为0
• ones：创建数据全为1
• empty：创建数据接近0
• arrange：按指定范围创建数据
• linspace：创建线段

创建数组

a = np.array([2,23,4])  # list 1d
print(a)
# [2 23 4]

指定数据类型

a = np.array([2,23,4],dtype=np.int64)
print(a.dtype)
# int 64

创建特定数据

a = np.array([[2,23,4],[2,32,4]])  # 2d 矩阵 2行3列
print(a)
"""
[[ 2 23  4][ 2 32  4]]
"""

numpy基本运算

import numpy as np
a=np.array([10,20,30,40])   # array([10, 20, 30, 40])
b=np.arange(4)              # array([0, 1, 2, 3])
c=a-b  # array([10, 19, 28, 37])
c=a+b   # array([10, 21, 32, 43])
c=a*b   # array([  0,  20,  60, 120])
c=b**2  # array([0, 1, 4, 9])
c=10*np.sin(a)  
# array([-5.44021111,  9.12945251, -9.88031624,  7.4511316 ])
print(b<3)  
# array([ True,  True,  True, False], dtype=bool)
a=np.array([[1,1],[0,1]])
b=np.arange(4).reshape((2,2))print(a)
# array([[1, 1],
#       [0, 1]])print(b)
# array([[0, 1],
#       [2, 3]])
c_dot = np.dot(a,b)
# array([[2, 4],
#       [2, 3]])
c_dot_2 = a.dot(b)
# array([[2, 4],
#       [2, 3]])a=np.random.random((2,4))
print(a)
# array([[ 0.94692159,  0.20821798,  0.35339414,  0.2805278 ],
#       [ 0.04836775,  0.04023552,  0.44091941,  0.21665268]])
np.sum(a)   # 4.4043622002745959
np.min(a)   # 0.23651223533671784
np.max(a)   # 0.90438450240606416
print("a =",a)
# a = [[ 0.23651224  0.41900661  0.84869417  0.46456022]
# [ 0.60771087  0.9043845   0.36603285  0.55746074]]print("sum =",np.sum(a,axis=1))
# sum = [ 1.96877324  2.43558896]print("min =",np.min(a,axis=0))
# min = [ 0.23651224  0.41900661  0.36603285  0.46456022]print("max =",np.max(a,axis=1))
# max = [ 0.84869417  0.9043845 ]A = np.arange(2,14).reshape((3,4)) 
# array([[ 2, 3, 4, 5]
#        [ 6, 7, 8, 9]
#        [10,11,12,13]])print(np.argmin(A))    # 0
print(np.argmax(A))    # 11print(np.mean(A))        # 7.5
print(np.average(A))     # 7.5print(np.cumsum(A)) 
# [2 5 9 14 20 27 35 44 54 65 77 90]print(np.diff(A))    
# [[1 1 1]
#  [1 1 1]
#  [1 1 1]]A = np.arange(14,2, -1).reshape((3,4)) 
# array([[14, 13, 12, 11],
#       [10,  9,  8,  7],
#       [ 6,  5,  4,  3]])print(np.sort(A))
# array([[11,12,13,14]
#        [ 7, 8, 9,10]
#        [ 3, 4, 5, 6]])print(np.transpose(A))    
print(A.T)
# array([[14,10, 6]
#        [13, 9, 5]
#        [12, 8, 4]
#        [11, 7, 3]])
# array([[14,10, 6]
#        [13, 9, 5]
#        [12, 8, 4]
#        [11, 7, 3]])

numpy的合并

import numpy as np
A = np.array([1,1,1])
B = np.array([2,2,2])print(np.vstack((A,B)))    # vertical stack
"""
[[1,1,1][2,2,2]]
"""D = np.hstack((A,B))       # horizontal stack
print(D)
# [1,1,1,2,2,2]C = np.concatenate((A,B,B,A),axis=0)
print(C)
"""
array([[1],[1],[1],[2],[2],[2],[2],[2],[2],[1],[1],[1]])
"""
D = np.concatenate((A,B,B,A),axis=1)
print(D)
"""
array([[1, 2, 2, 1],[1, 2, 2, 1],[1, 2, 2, 1]])
"""

numpy的分割

import numpy as np
A = np.arange(12).reshape((3, 4))
print(A)
"""
array([[ 0,  1,  2,  3],[ 4,  5,  6,  7],[ 8,  9, 10, 11]])
"""print(np.split(A, 2, axis=1))
"""
[array([[0, 1],[4, 5],[8, 9]]), array([[ 2,  3],[ 6,  7],[10, 11]])]
"""print(np.split(A, 3, axis=0))
# [array([[0, 1, 2, 3]]), array([[4, 5, 6, 7]]), array([[ 8,  9, 10, 11]])]print(np.array_split(A, 3, axis=1)) #不等量分割
"""
[array([[0, 1],[4, 5],[8, 9]]), array([[ 2],[ 6],[10]]), array([[ 3],[ 7],[11]])]
"""

线性规划算法原理

单纯形法

在之前我们提出了线性规划的标准形式，但注意这种标准形式是针对程序设计工具而言的。如果读者有凸优化理论的背景可能会感到狐疑，说：为什么我看到的标准形式和你这里写的不太一样？是的，经典的凸优化教材会把模型写成另外一种形式：

为了和2.1当中提出的标准形式区分，我们暂且把这种形式称作规范形式。规范形式求的是函数的极大值，并且把不等关系和等式关系统一为等式关系方便求解。读者朋友可能会有些疑惑，说：不等式怎么可以充当为方程呢？这就是一种数学思想。可能读者朋友可以理解方程是不等式的特例，但不一定理解不等式也可以视作方程的特例，我举个例子。比如对于不等式 2𝑎+3𝑏+𝑐<10 ，左边比右边小，但是小多少呢？我们把这个差额记作d，左边如果补上这个差额就可以写作 2𝑎+3𝑏+𝑐+𝑑=10 ，这样就转化成了等式。这里的d被称为松弛变量。包括决策变量的上下界 𝑙_𝑏 和 𝑢_𝑏 也会被转化为不等关系引入松弛量。

在单纯形法中，我们解决问题通常从理论上都会把问题转换为规范形式来求解，对每一个不等式都引入一个松弛变量去增广我们的原问题。但这些松弛变量不会出现在目标函数当中。

注意：在程序设计中我们输入的是它的标准形式，而在matlab底层以规范形式进行运算，只是从标准形式到规范形式这个操作我们看不见。不等式条件中增广了n个松弛变量的同时等式条件也会增广，只不过在增广后的等式条件中松弛变量的系数都是0。

例题1

此时，将𝑥4替换为𝑥1，将𝑥1系数5填入表格中，原来的表格都减去第二个约束条件除以𝑥1在第二行约束中的系数2得到新的系数与常数项。此时，𝑧等于𝑥3的系数0乘以第一行新约束加x_1的系数5乘以第二行新约束，得到新的z结果。𝑢=𝑓−𝑧得到𝑢的结果，所有的数都是非正数，此时迭代就结束了。常数项对应的1,4其实就是解，计算可以得到极值。 单纯形法的实现代码如下：

import numpy as np
def pivot(d,bn):l = list(d[0][:-2])jnum = l.index(max(l)) #转入编号m = []for i in range(bn):if d[i][jnum] == 0:m.append(0.)else:m.append(d[i][-1]/d[i][jnum])inum = m.index(min([x for x in m[1:] if x!=0]))  #转出下标s[inum-1] = jnumr = d[inum][jnum]d[inum] /= rfor i in [x for x in range(bn) if x !=inum]:r = d[i][jnum]d[i] -= r * d[inum]  
#定义基变量函数
def solve(d,bn):flag = Truewhile flag:if max(list(d[0][:-1])) <= 0: #直至所有系数小于等于0flag = Falseelse:pivot(d,bn)            
def printSol(d,cn):for i in range(cn - 1):if i in s:print("x"+str(i)+"=%.2f" % d[s.index(i)+1][-1])else:print("x"+str(i)+"=0.00")print("objective is %.2f"%(-d[0][-1]))
d = np.array([[5,1,0,0,0],[1,1,1,0,5],[2,1/2,0,1,8]])
(bn,cn) = d.shape
s = list(range(cn-bn,cn-1)) #基变量列表
solve(d,bn)
printSol(d,cn)

此时我们可以看到:𝑥1=4，𝑥2=0，𝑥3=1，𝑥4=0，极值为20。

内点法

单纯形法之所以需要遍历所有顶点才能获得最优解，归根结底还是在于单纯形算法的搜索过程是从一个顶点出发，然后到下一个顶点，就这样一个顶点一个顶点的去搜寻最优解。单纯形算法的搜索路径始终是沿着多面体的边界的。显然，当初始点离最优点的距离很远的时候单纯形算法的搜索效率就会大大降低。 能否直接从多边形内部打进来呢？这就需要用到内点法，如图所示：


使用内点法实现上述用例的代码如下：

import numpy as np
def Interior_Point(c,A,b):# 当输入的c,A,b有缺失值时，输出错误原因，函数执行结束if c.shape[0] != A.shape[1]:print("A和C形状不匹配")return 0if b.shape[0] != A.shape[0]:print("A和b形状不匹配")return 0# 初值的设置x=np.ones((A.shape[1],1)) # 将x的初值设为1v=np.ones((b.shape[0],1))  # 将v的初值设为1lam=np.ones((x.shape[0],1))  # 将lam的初值设为1one=np.ones((x.shape[0],1))mu=1 # 将mu的初值设为1n=A.shape[1]x_=np.diag(x.flatten()) # 将x转换为对角矩阵lam_=np.diag(lam.flatten()) # 将lam转换为对角矩阵# 初始的F，r=Fr1=np.matmul(A,x)-br2=np.matmul(np.matmul(x_,lam_),one)-mu*oner3=np.matmul(A.T,v)+c-lamr=np.vstack((r1,r2,r3))F=r# 求得r1、r2、r3的初始范数n1=np.linalg.norm(r1)n2=np.linalg.norm(r2)n3=np.linalg.norm(r3)# nablaF中零矩阵和单位阵的设置zero11=np.zeros((A.shape[0],x.shape[0]))zero12=np.zeros((A.shape[0],A.shape[0]))zero22=np.zeros((x.shape[0],A.shape[0]))zero33=np.zeros((A.shape[1],A.shape[1]))one31=np.eye(A.shape[1])tol=1e-8 # 设置最优条件的容忍度t=1alpha = 0.5while max(n1,n2,n3)>tol:print("-----------------step",t,"-----------------")# F的Jacobian矩阵nablaF = np.vstack((np.hstack((zero11, zero12, A)), np.hstack((x_, zero22, lam_)), np.hstack((-one31, A.T, zero33))))# F+nablaF@delta=0,解方程nablaF@delta=-rdelta = np.linalg.solve(nablaF, -r)  # 解方程，求出delta的值delta_lam = delta[0:lam.shape[0], :]delta_v = delta[lam.shape[0]:lam.shape[0] + v.shape[0], :]delta_x = delta[lam.shape[0] + v.shape[0]:, :]# 更新lam、v、x、mualpha=Alpha(c,A,b,lam,v,x,alpha,delta_lam,delta_v,delta_x)lam=lam+alpha*delta_lamv=v+alpha*delta_vx=x+alpha*delta_xx_ = np.diag(x.flatten())  # 将x转换为对角矩阵lam_ = np.diag(lam.flatten())  # 将lam转换为对角矩阵mu=(0.1/n)*np.dot(lam.flatten(),x.flatten()) #更新mu的值# 计算更新后的Fr1 = np.matmul(A, x) - br2 = np.matmul(np.matmul(x_, lam_), one) - mu * oner3 = np.matmul(A.T, v) + c - lamr = np.vstack((r1, r2, r3))F = r# 计算更新后F的范数n1 = np.linalg.norm(r1)n2 = np.linalg.norm(r2)n3 = np.linalg.norm(r3)t=t+1print("x的取值",x.flatten())print("v的取值",v.flatten())print("lam的取值",lam.flatten())print("mu的取值",mu)print("alpha的取值",alpha)z=(c.T @ x).flatten()[0]print("值为",z)print("##########################找到最优点##########################")print("x的取值",x.flatten())print('最优值为',z)# 寻找alphadef Alpha(c,A,b,lam,v,x,alpha,delta_lam,delta_v,delta_x):alpha_x=[]alpha_lam=[]for i in range(x.shape[0]):if delta_x.flatten()[i]<0:alpha_x.append(x.flatten()[i]/-delta_x.flatten()[i])if delta_lam.flatten()[i]<0:alpha_lam.append(lam.flatten()[i]/-delta_lam.flatten()[i])if len(alpha_x)==0 and len(alpha_lam)==0:return alphaelse:alpha_x.append(np.inf)alpha_lam.append(np.inf)alpha_x = np.array(alpha_x)alpha_lam= np.array(alpha_lam)alpha_max = min(np.min(alpha_x), np.min(alpha_lam))alpha_k = min(1,0.99*alpha_max)return alpha_kc = np.array([-5, -1, 0,0]).reshape(-1, 1)
A = np.array([[1, 1, 1, 0], [2, 0.5, 0, 1]])
b = np.array([5, 8]).reshape(-1, 1)
Interior_Point(c,A,b)