建模代码python —— 规划类问题

叮叮车不是车

已于 2022-09-02 10:05:40 修改

阅读量1.4k

点赞数 6

分类专栏：数学建模文章标签： python

于 2022-08-02 20:08:35 首次发布

本文链接：https://blog.csdn.net/sinat_56238820/article/details/126125333

版权

数学建模专栏收录该内容

10 篇文章 8 订阅

订阅专栏

一、规划类问题

规划是运筹学的一个重要分支，主要研究数值最优化问题。三个主要构成要素为决策变量、目标函数以及约束条件。

1.线性规划

线性规划变量的取值范围是实数范围，也就是说最优解所对应的变量取值可以是整数也可以是小数。

即目标函数（max,min）和约束条件（s.t.）都是线性的规划

求解前化成标准形式：

$s.t.\begin{Bmatrix} min c^{T}x\\ Ax\leq b\\ Aeq*X=beq\\ lb\leq x\leq ub\end{Bmatrix}$

（ $lb$ 和 $ub$ 分别为x的上界和下界）

①scipy库求解

from scipy import optimize
import numpy as np
#求解函数
res=optimize.linprog(c,A,b,Aeq,beq,LB,UB,X0,OPTIONS)
#目标函数最小值
print(res.fun)
#最优解
print(res.x)

例子：

from scipy import optimize
import numpy as np

# 约束式（大于等于和小于等于的式子）都写成二维格式
c=np.array([2,3,-5])
A=np.array([[-2,5,-1],[1,3,1]])
B=np.ar ray([-10,12])
Aeq=np.array([[1,1,1]]) #与A匹配矩阵相乘所以为二维的
Beq=np.array([7])

# 求解
res=optimize.linprog(-c,A,B,Aeq,Beq)
print(res)

输出：

#目标函数最小值 fun
  fun: -14.571428571428571

       message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'
           nit: 3
         slack: array([0.        , 3.85714286])
        status: 0
       success: True

#最优解 x
    x: array([6.42857143, 0.57142857, 0.        ])

②pulp库求解（不用变成标准形式）

例子：

from scipy import optimize
import pulp

#目标函数的系数
z=[2,3,1]
a=[[1,4,2],[3,2,0]]
b=[8,6]

#确定最大最小化问题，最大化只要把Min改成Max即可
m=pulp.LpProblem(sense=pulp.LpMinimize)

#定义三个变量放到列表中
x=[pulp.LpVariable(f'x{i}',lowBound=0) for i in [1,2,3]]

#定义目标函数，lpDot可以将两个列表的对应位相乘再加和
#相当于z[0]*x[0]+z[1]*x[1]+z[2]*x[2]
m+=pulp.lpDot(z,x)

#设置约束条件
for i in range(len(a)):
    m+=(pulp.lpDot(a[i],x)>=b[i])

    #求解
    m.solve()

    #输出结果
    print(f'优化结果：{pulp.value(m.objective)}')
    print(f'参数取值：{[pulp.value(var) for var in x]}')

③运输问题

某商品有m个产地、n个销地，各产地的产量分别为a1,....，am，各销地的需求量分别为b1，...，bn。若该商品由i产地运到j销地的单位运价为cij，问应该如何调用才能使运费最省

引入变量xij，其取值为由i产地运往j销地的该产品数量，模型为：

题目：

import pulp
import numpy as np
from pprint import pprint
# 写函数
def transportation_problem(costs, x_max, y_max):  # cost=cij,x_max和y_max是对i、j的大小约束
    row = len(costs)
    col = len(costs[0])
    prob = pulp.LpProblem('Transportation Problem', sense=pulp.LpMaximize)
    var = [[pulp.LpVariable(f'x{i}{j}', lowBound=0, cat=pulp.LpInteger) for j in range(col)] for i in
           range(row)]  # 先i后j lowBound是下界
    # 重点 将多维转变为一维，如果函数的类型是列表，就进行中括号中的操作，若不是则先将x转化为列表形式
    flatten = lambda x: [y for l in x for y in flatten(l)] if type(x) is list else [x]
    prob += pulp.lpDot(flatten(var), costs.flatten())
    for i in range(row):
        prob += (pulp.lpSum(var[i]) <= x_max[i])
    for j in range(col):
        prob += (pulp.lpSum([var[i][j] for i in range(row)]) <= y_max[j])
        prob.solve()
        return {'objective': pulp.value(prob.objective), # 最大值
                'var': [[pulp.value(var[i][j]) for j in range(col)] for i in range(row)]}  # 目标解


# 调用函数
if __name__ == '__main__':
    costs = np.array([[500, 550, 630, 1000, 800, 700],
                      [800, 700, 600, 950, 900, 930],
                      [1000, 960, 840, 650, 600, 700],
                      [1200, 1040, 980, 860, 880, 780]])
    max_plant = [76, 88, 96, 40]
    max_cultivation = [42, 56, 44, 39, 60, 59]
    res = transportation_problem(costs, max_plant, max_cultivation)
    print(f'最大值为{res["objective"]}')
    print('各变量的取值为：')
    pprint(res['var'])

2.整数规划

在线性规划的基础上增加了部分变量为整数的约束，即整数规划变量的取值范围是整数范围
求解的基本框架使分支定界法
去除整数约束得到“松弛模型”，使用线性规划的方法求解
若有某个变量不是整数：在松弛模型上分别添加约束：x≤floor(A)和x≥ceil(A)，然后再分别求解，这个过程叫做分支。当节点求解过程中所有的变量都是整数时，停止分支。这样的迭代，形成了一棵树。
叶子节点产生后，相当于给问题定了一个下界。
每次新产生叶子节点就更新下界
求解过程中一旦某个节点的目标函数小于这个下界，就直接舍弃

分支定界代码

import math
from scipy.optimize import linprog
import sys  # 导入系统库

# 函数
def integerPro(c, A, b, Aeq, beq, t=1.0E-12):  # t无限趋近于0
    res = linprog(c, A_ub=A, b_ub=b, A_eq=Aeq, b_eq=beq)
    if (type(res.x) is float):  # 如果最优解中出现float型
        bestX = [sys.maxsize] * len(c)
    else:
        bestX = res.x
    bestVal = sum([x * y for x, y in zip(c, bestX)])  # c和bestX按位配比

    # x-小于x的最大整数 看是不是0 ；大于x的最大整数-x看是不是0 ---- 判断是否为整数
    # 因为有可能-完=1，所以用or
    if all(((x - math.floor(x)) < t or (math.ceil(x) - x < t) for x in bestX)):
        return (bestVal, bestX)
    else:
        # 获取不是整数的位置
        ind = [i for i, x in enumerate(bestX) if (x - math.floor(x)) > t and (math.ceil(x) - x) > t][0]
        newCon1 = [0] * len(A[0])
        newCon2 = [0] * len(A[0])
        newCon1[ind] = -1
        newCon2[ind] = 1
        newA1 = A.copy()
        newA2 = A.copy()
        newA1.append(newCon1)
        newA2.append(newCon2)
        newB1 = b.copy()
        newB2 = b.copy()
        newB1.append(-math.ceil(bestX[ind]))
        newB2.append(math.floor(bestX[ind]))
        r1 = integerPro(c, newA1, newB1, Aeq, beq)
        r2 = integerPro(c, newA2, newB2, Aeq, beq)
        if r1[0] < r2[0]:
            return r1
        else:
            return r2

# 值
c = [3, 4, 1]
A = [[-1, -6, -2], [-2, 0, 0]]
b = [-5, -3]
Aeq = [[0, 0, 0]]
beq = [0]
print(integerPro(c, A, b, Aeq, beq))

3.非线性规划

目标为凸函数or非凸函数
凸函数：常用库
非凸函数：①求导求极值 ②神经网络、深度学习 ③scipy.optimize.minimize
scipy.optimize.minimize(fun,x0,args=(),method=None,
                        jac=None,hess=None,hessp=None,bounds=None,constraints=(),
                        tol=None,callbacks=None,options=None)
fun：求最小值的目标函数
args：常数值
method：求极值方法，一般默认
constraints：约束条件
x0：变量的初始猜测值，minimize是局部最优

① 计算1/x+x的最小值

②计算（2+x1）/（1+x2）- 3 * x1 + 4 * x3

二、数值逼近类问题

1. 一维插值

插值函数经过样本点，拟合函数一般基于最小二乘法尽量靠近所有样本点穿过。常见插值方法：

拉格朗日插值法：当节点数n较大时，拉格朗日插值多项式次数较高，可能收敛不一致，且计算复杂。（高插值带来误差的震动现象称为龙格现象
分段插值法：虽然收敛，但光滑性较差
样条插值法：由于样条插值可以使用低阶多项式样条实现较小的插值误差，从而避免龙格现象