数学建模（线性规划）python版

数学建模——规划模型

数学规划是运筹学的一个重要分支、而线性规划又是数学规划的一部分主要内容，所有实际问题都可以归总为“线性规划”问题。线性规划（linear programming，LP)有比较完善的理论基础和有效的解决方法。在实际问题中有极其广泛的应用。

一、线性规划模型

1、建立线性规划模型的步骤

规划问题的数学模型由三个要素组成
①决策变量：是问题中要确定的未知量，用于表明规划问题中用数量表示的方案、措施等，可由决策者决定和控制
②目标函数：是决策变量的函数，优化目标通常是求该函数的最大值或最小值。
③约束条件：是决策变量的取值所受到的约束和限制条件，通常用含有决策变量的等式或不等式表示。

建立线性规划模型通常需要以下三个步骤

（1）分析问题，找到决策变量

（2）根据问题所给的条件，找到决策变量必须满足的一组线性等式或者不等式约束，即确定约束条件

（3）根据问题的目标、构造关于决策变量的一个线性函数，即目标函数

有了决策变量、约束条件和目标函数这三个要素之后，一个线性规划模型就建立起来了。

2、线性规划模型的形式

线性规划模型的一般形式为
max(或min) z=$c1x1+c2x2+....+c_nx_n$
$s.t.\left\{\begin{aligned}a_{11}x_{1}+a_{12}x_{2}+...+a_{1n}x_{n}\leq(或=、\geq)b_{1}\\a_{11}x_{1}+a_{12}x_{2}+...+a_{1n}x_{n}\leq(或=、\geq)b_{2}\\a_{m1}x_{1}+a_{m2}x_{2}+...+a_{mn}x_{n}\leq(或=、\geq)b_{m}\\x{1},x{2},...x{n}\leq0\end{aligned}\right.$

或者简写为

max（或min) z=$\displaystyle \sum^{n}_{y=1}{C_{j}x_{j}}$

$s.t.\left\{\begin{aligned}\\\sum^{n}_{j=1}a_{ij}\leq(或=,\geq)b_{i} ,i=1,2,3,...m,\\x_{j}\geq0,j=1,2,3,...n\end{aligned}\right.$

其向量表示形式为

max(或min) z=C^Tx

$s.t.\left\{\begin{aligned}\\\sum^{n}_{j=1}P_{ij}\leq(或=,\geq)b\\x\geq0\end{aligned}\right.$

其矩阵表示形式为

max(或min) z=C^Tx

$s.t.\left\{\begin{aligned}\\Ax\leq(或=,\geq)b\\x\geq0\end{aligned}\right.$X

其中c=$(c_1,c2,...c_n)^T$为目标函数的系数向量；又称为价值向量；x=$(x_1,x_2,...x_n)^T$称为决策向量；A=$(a_{ij})_{mxn}$称为约束方程组的系数矩阵；而P$j$=$(a_{1j},a_{2j},a_{3j},....a_{mj})^T$,j=1,2,3,4...n,n为A的列向量，又称为约束方程组的系数向量；b=$(b_{1j},b_{2j}....b_{mj})^T$称为约束方程组的常数向量。

3、用python求解线性规划模型

（1）cvxpy库，用于解决凸优化问题，即用于求解凸规划模型,具体用法，可参考cvxpy 常用功能汇总 - 哔哩哔哩 (bilibili.com)

（2）scipy库中的optimize linprog模块，可参考scipy.optimize.linprog函数参数最全详解_佐佑思维的博客-CSDN博客_scipy.optimize.linprog

例如：max z=$70_{x1}+50_{x2}+60_{x3}$

$s.t.\left\{\begin{aligned}\\2x_{1}+4x_{2}+3x_{3}\leq150,\\3x_{1}+x_{2}+53_{3}\leq160,\\7x_{1}+3x_{2}+5x_{3}\leq200,\\x_{i}\geq0,i=1,2,3\end{aligned}\right.$
cvxpy的具体链接：cvxpy 常用功能汇总 - 哔哩哔哩 (bilibili.com)

# 使用cvxpy来实现线性规划问题 
import cvxpy as cp
import numpy as np
c=np.array([70,50,60])  # 定义目的向量
a=np.array([[2,4,3],[3,1,5],[7,3,5]]) #定义约束矩阵
b=np.array([150,160,200])           # 定义约束条件的右边向量
x=cp.Variable(3,pos=True)           # 定义3个决策变量(x1,x2,x3)
obj=cp.Maximize(c@x)                # 构造约束条件最大值(70x1+50x2+60x3)
cons=[a@x<=b]                       # 构造约束条件
pro=cp.Problem(obj,cons)            # 求解问题
pro.solve(solver='GLPK_MI',verbose=True)			
#用不同的求解器求出来的几个是不一样的，用GLPKMI解决问题，verbose=True表明显示具体求解信息
print('最优解是',x.value)
print('最优值是',pro.value)
# Minimize 表示最小值，Maximize表示最大值，出现等于的约束条件可以将约束条件卸载cons里面

# 输出的结果：
# 最优的解 [15.90909091 29.54545455 -0.        ]
# 最优的值 2590.909090909091

cvxpy的求解器(solver)

LP:指线性规划 QP:指二次规划（二次函数） SOCP:指二次锥规划

SDP:半正定规划 EXP:指数规划 POW:幂规划 MIP:混合整数规划

scipy的具体链接：scipy.optimize.linprog函数参数最全详解_佐佑思维的博客-CSDN博客_scipy.optimize.linprog

# 使用scipy库实现线性规划
from scipy import optimize as op
import numpy as np
c=np.array([-70,50,60])      # 目标函数
A_ub=np.array([[2,4,3],[3,1,5],[7,3,5]])    # 约束条件的左边
B_ub=np.array([[150,160,200]])              # 约束条件的右边
x1=(0,None)                                 # 三个变量及其范围(0到正无穷）
x2=(0,None)
x3=(0,None)
res=op.linprog(-c,A_ub,B_ub,bounds=(x1,x2,x3),method='revised simplex')      # 使用method='revised simplex'高精度运算
print(res)

# 输出的结果：
#    con: array([], dtype=float64)
#     fun: -2590.909090909091
# message: 'Optimization terminated successfully.'
#     nit: 2
#   slack: array([ 0.        , 82.72727273,  0.        ])    
#  status: 0
# success: True
#       x: array([15.90909091, 29.54545455,  0.        ])    
#  求解的值 2590.909090909091
#  求解的x [15.90909091 29.54545455  0.        ]

二、整数线性规划

1、整数规划的类型

在一些问题（例如汽车企业在制定生产计划是，要求所生产的不同类型的汽车数量必须是整数）等在问题上，我们必须要求一部分或者全部决策变量必须是整数值的规划问题称为整数规划（interger programming IP)。

从决策变量的取值范围来看，整数规划通常分为以下几种类型

（1）纯整数规划：全部决策变量都必须整数值的整数规划模型

（2）混合整数规划：决策变量中优一部分必须去整数值，另一部分可以不取整数的整数规划模型

（3）0-1整数规划：决策变量只能取0或1的整数规划

2、整数线性规划模型一般形式为

max（或min) z=$\displaystyle \sum^{n}_{y=1}{C_{j}x_{j}}$

$s.t.\left\{\begin{aligned}\sum^n_{j=1}a_{ij}x{j}\leq(或=,\geq)b_{i},i=1,2,....m\\x_{j}\geq0,j=1,2,...n\\x_1,x_2,....x_n 中部分或者全部取整数\end{aligned}\right.$

3、python求解模型

整数处理方式：

1、纯整数规划问题及混合整数规划问题

x=cp.Variable(3,pos=True,interger=True)  # interger 化为整数形式

2、0-1整数规划

# 在约束条件中加入 0<=x,x<=1
x=cp.Variable(3,pos=True,interger=True)
cons=[a@x<=b,0<=x,x<=1]   

例题：

min z=$\displaystyle\sum^4_{i=1}\sum^5_{j=1}c_{ij}x{ij}$

$s.t.\left\{\begin{aligned}\sum^4_{i=1}a_{ij}\leq=1,i=1,2,....5\\\sum^5_{i=1}x_{ij}\leq0,j=1,2,3,4\\x_{ij}=0或1,i=1,2,3,4,j=1,2,...5\end{aligned}\right.$

cvxpy的解法

import cvxpy as cp
import numpy as np
c=np.array([[15,13.8,12.5,11,14.3],[14.5,14,13.2,10.5,15],[13.8,13,12.8,11.3,14.6],[14.,13.6,13,11.6,14]])
x=cp.Variable((4,5),integer=True)      # 定义目标变量，全为整数
obj=cp.Minimize(cp.sum(cp.multiply(c,x)))       # 求解最小值 构造目标函数
con=[cp.sum(x,axis=1)<=2,cp.sum(x,axis=0)==1,0<=x,x<=1]     # 构造约束条件
# cp.sum(x,axis=1)表示在每一列中不能大于2，0<=x,x<=1,表示0-1整数问题处理方法
pro=cp.Problem(obj,con)
pro.solve(solver='GLPK_MI')
print('最优解：',pro.value)
print('最优值',x.value)

三、非线性规划模型

我们遇到的规划模型中的目标函数和约束条件均为决策变量的线性函数，我们称其为线性规划模型。如果在目标函数或者约束条件中含有决策变量的非线性函数，则该模型为非线性模型（nonlinear progamming,NIP)。

1、非线性规划模型形式

非线性规划模型的一般形式为

min f(x)

$s.t.\left\{\begin{aligned}g_i(x)\leq0,i=1,2,...,p\\h_j(x)\geq0,i=1,2,...q\end{aligned}\right.$

其中x=$(x_1,x_2,...x_n)^T$是n维欧氏空间$R^n$中的一个n位向量，函数f（x),$g_i(x)$,i=1,2,...,p和$hi(x)$,j=1,2...q是定义在$R^n$上的实值函数，且其中至少一个是x的非线性函数。

特别，若在模型中，p=q=0，我们称其为无约束非线性规划或无约束优化，即求多元函数的极值问题；若目标函数f(x)是一个二次函数，而$g_i(x)$,i=1,2,...p和$h_j(x)$,j=1,2,...q都是线性函数，我们称其为二次归规划。若目标函数f(x)为凸函数,则均为线性函数，我们称其为凸规划。

2、用python求解非线性规划模型

对于一般的非线性问题，由于不是凸规划，就不能用cvxpy库求解，求解一般的非线性规划问题使用scipy.optimize模块的minimize函数求解。[minimize的使用](python 非线性规划方式(scipy.optimize.minimize)_python_脚本中心 - 编程客栈 (cppcns.com))

min T=$\displaystyle\sum^3_{i_=1}t_i$

$s.t.\left\{\begin{aligned}(t1+t3)1.5cosa_1=400\\t_21.5a_2=760\\(1.47-1.5sina_1)(t_1+t_3)+(2.11-1.5sina_1)t2=1000\\t1=t3\\t_i\geq0,i=1,2,3\\o<a_i<Π/2,i=1,2\end{aligned}\right.$

import numpy as np
from scipy.optimize import minimize
obj=lambda x:2*x[0]+x[1]        # 定义目标的函数
con=[{'type':'eq','fun':lambda x:2*x[0]*1.5*np.cos(x[2]-400)},
{'type':'eq','fun':lambda x:x[1]*1.5*np.cos(x[3])-760},
{'type':'eq','fun':lambda x:(1.46-1.5*np.sin(x[2]))*2*x[0]+(2.11-1.5*np.sin(x[3]))*x[1]-1000}]      # 定义约束条件
b1=[0]*4
b2=[np.inf,np.inf,np.pi/2,np.pi/2]
bd=list(zip(b1,b2))                         # 将两个绑定在一起，以元组的形式
x=minimize(obj,np.random.rand(4),constraints=con,bounds=bd)  
# 根据minimize的格式，np.random.rand(4)表示初始化四个数，constraints约束条件，method表示方法
print(x)
print('角度',x.x[2:]*180/np.pi)