数学建模专题1.1线性规划模型（LP问题）

一:前言

数学规划是运筹学的一个重要分支，而线性规划又是数学规划中的一部分主要内容，很多实际问题都可以归结为"线性规划（Linear Programming，LP）"问题。

二:解决步骤

规划模型由3个要素组成：

1. 决策变量:问题中要确定的未知量，用于表明规划问题中的方案，措施等
2. 目标函数:是决策变量的函数，优化目标通常是求该函数的最大值或最小值
3. 约束条件:是决策变量的约束和限制条件，通常由等式和不等式组成

解决步骤:

1. 第一步:分析问题，找出决策变量
2. 第二步:找出约束条件，即决策变量必须满足的一组线性等式或不等式约束
3. 第三步:根据问题的目标，构造一个关于决策变量的线性函数，即目标函数

三:线性规划模型的一般形式

max(或min)

       

四:例题求解及其代码

主要利用cvxpy库求解,命令行输入以下代码即可下载cvxpy库（若想下载其他库，修改cvxpy即可）

pip install -i https://pypi.tuna.tsinghua.edu.cn/simple cvxpy

例一

某企业利用两种原材料A和B生产三种产品P₁、P₂和P₃。已知每生产 1 公斤的产品所消耗的原材料A、B的数量 (单位：公斤) 和花费的加工时间C (单位：小时)，每公斤产品销售后所带来的利润 (单位：元) 以及每天可用的资源的数量如表 4.1 所示，则该企业应该如何制订每天的生产计划，才能使所获利润达到最大?

模型建立

 代码实现

import cvxpy as cp
from numpy import array
A = array([[2,4,3],[3,1,5],[7,3,5]]) #约束矩阵
b = array([150,160,200])             #约束条件的右边向量
c = array([70,50,60])                #目标向量
x = cp.Variable(3,pos=True) #定义决策变量范围为非负数
obj = cp.Maximize(c@x) #构造目标函数
cons = [A@x<=b] #构造约束条件
prob =cp.Problem(obj,cons)
prob.solve(solver='GLPK_MI') #求解问题
print("最优解为:",x.value)
print("最优值为:",prob.value)

例二

某部门在今后五年内考虑给下列项目投资，已知：
项目A，从第一年到第四年每年年初需要投资，并于次年末回收本利115%;
　　项目B，从第三年初需要投资，到第五年末能回收本利 125%，但规定最大投资额不超过4万元；
　　项目C，第二年初需要投资，到第五年末能回收本利 140%，但规定最大投资额不超过3万元；
项目D，五年内每年初可购买公债，于当年末归还， 并加利息 6%
　　该部门现有资金 10 万元， 问它应如何确定给这些项目每年的投资额，使到第五年末拥有的资金的本利总额为最大?

模型建立

 

 

代码实现

安装cvxop库可以使用GLPK_MI求解器

import cvxpy as cp
from numpy import array

c = array([1.15,1.40,1.25,1.06])
x = cp.Variable((5,4),pos=True) #Xij(5,4),非负数
obj = cp.Maximize(1.15*x[3,0]+1.4*x[1,2]+1.25*x[2,1]+1.06*x[4,3])
cons = [x[0,0]+x[0,3]==100000,#第一年
        x[1,0]+x[1,2]+x[1,3]==1.06*x[0,3],# 第二年
        x[2,0]+x[2,1]+x[2,3]==1.15*x[0,0]+1.06*x[1,3],#第三年
        x[3,0]+x[3,3]==1.15*x[1,0]+1.06*x[2,3],#第四年
        x[4,3]==1.15*x[2,0]+1.06*x[3,3],#第五年
        x[2,1]<=40000,
        x[1,2]<=30000
]#约束条件
prob = cp.Problem(obj,cons)
prob.solve(solver='CPLEX')#对应的求解器也可以是GLPK_MI
print("最优解:\n",x.value)
print("最优值:",prob.value)

例三

 捷运公司在下一年度的1~4月的4个月内拟租用仓库堆放物资。已知各月份所需仓库面积列于错误！未找到引用源。。仓库租借费用随合同期而定，期限越长，折扣越大，具体数字见表4.2。租借仓库的合同每月初都可办理，每份合同具体规定租用面积和期限。因此该公司可根据需要，在任何一个月初办理租借合同。每次办理时可签一份合同也可签若干份租用面积和租借期限不同的合同，试确定该公司签订租借合同的最优决策，目的是使所付租借费用最小

 模型建立

 代码实现

import cvxpy as cp
from numpy import array
x =cp.Variable((4,4),pos=True)
obj =cp.Minimize(2800*(x[0,0]+x[1,0]+x[2,0]+x[3,0])+4500*(x[0,1]+x[1,1]+x[2,1])+6000*(x[0,2]+x[1,2])+7300*x[0,3])
cons =[x[0,0]+x[0,1]+x[0,2]+x[0,3]>=15,
       x[0,1]+x[0,2]+x[0,3]+x[1,0]+x[1,1]+x[1,2]>=10,
       x[0,2]+x[0,3]+x[1,1]+x[1,2]+x[2,0]+x[2,1]>=20,
       x[0,3]+x[1,2]+x[2,1]+x[3,0]>=12]
'''上述代码可以简写为
obj = cp.Minimize(2800*sum(x[:,0])+4500*sum(x[:3,1])+6000*sum(x[:2,2])+7300*x[0,3])
cons =[sum(x[0,:])>=15,
       sum(x[0,1:])+sum(x[1,:3])>=10,
       sum(x[0,2:])+sum(x[1,1:3])+sum(x[2,:2])>=20,
       x[0,3]+x[1,2]+x[2,1]+x[3,0]>=12
]
'''
prob =cp.Problem(obj,cons)
prob.solve(solver='GLPK_MI')#对应的求解器也可以是CPLEX
print("最优解:\n",x.value)
print("最优值:",prob.value)

例四

代码实现

可以通过txt或者xlsx读入

import cvxpy as cp
import numpy as np
import pandas as pd
c = np.genfromtxt("data0.txt",dtype=float,max_rows=6,usecols=range(8))#读取前6行前8列数据
e = np.genfromtxt("data0.txt",dtype=float,max_rows=6,usecols=8)#产量,读取最后一列数据
d = np.genfromtxt("data0.txt",dtype=float,skip_header=6)#销量，读取最后一行数据
# 上述代码也可写为
# data0 =pd.read_excel("data1.xlsx",header=None)#读取没有表头的表格
# data = data0.values
# c = data[:-1,:-1]
# e = data[:-1,-1]
# d = data[-1,:-1]
x = cp.Variable((6,8),pos=True)#6行8列数据，非负数
obj =cp.Minimize(cp.sum(cp.multiply(c,x)))#目标的最小值函数
#此处只能用cp.multipy(),不能用*
con = [
    cp.sum(x,axis=0)==d,#对应列的和
    cp.sum(x,axis=1)<=e,#对应行的和
]
prob =cp.Problem(obj,con)
prob.solve(solver='GLPK_MI')
print("最优解:\n",x.value)
print("最优值:",prob.value)
#将数据写到Excel
res = pd.DataFrame(x.value)
res.to_excel("result.xlsx")

五：细节与反思

python中矩阵乘法

1. 使用NumPy库进行矩阵乘法：

   import numpy as np
   
   # 定义两个矩阵
   A = np.array([[1, 2], [3, 4]])
   B = np.array([[5, 6], [7, 8]])
   
   # 使用np.dot进行矩阵乘法
   result = np.dot(A, B)
   print(result)

2. 使用矩阵乘法运算符@：

   import numpy as np
   
   # 定义两个矩阵
   A = np.array([[1, 2], [3, 4]])
   B = np.array([[5, 6], [7, 8]])
   
   # 使用矩阵乘法运算符@
   result = A @ B
   print(result)

两个规格相同的矩阵对应元素相乘 

1. 使用*运算符进行对应元素相乘：

   import numpy as np
   
   # 定义两个矩阵
   A = np.array([[1, 2], [3, 4]])
   B = np.array([[5, 6], [7, 8]])
   
   # 对应元素相乘
   result = A * B
   print(result)

2. 使用numpy.multiply()函数进行对应元素相乘：

   import numpy as np
   
   # 定义两个矩阵
   A = np.array([[1, 2], [3, 4]])
   B = np.array([[5, 6], [7, 8]])
   
   # 对应元素相乘
   result = np.multiply(A, B)
   print(result)

还应该关注线性规划问题的灵敏度，可行域之类问题