数学建模part（1）：线性规划

线性规划的目标：寻求决策变量最优的组合方案

模型建立步骤：

1. 模型假设：需要假设不同决策变量之间相互独立，以及约束条件与价值向量的值是固定
2. 求解线性规划
   • 非线性尽可能转化为线性
   • 多目标规划：
     • 控制变量：逐个目标考虑（先固定一个目标）
     • 同时考虑两个目标，设置两个目标之间的权重分配
3. 灵敏性分析：对系数的波动影响可视化结果

符号说明

• 🍁 不要求掌握
• ✅要求掌握

解决工具：

• matlab
• python
• lingo

线性规划（linear programming） 三要素：

• 决策变量 ：控制量
• 目标函数：决策变量的线性函数
• 约束条件：限制范围

线性规划模型的一般形式：
$$\begin{gathered} max(min)z=\sum _{j=1}^n{c}_j{x}_j \\ s.t.\{\begin{array}{rlr}\sum _j^n{a}_{ij}{x}_j& \le (\ge ,=)b_i,& i=1,2,\dots ,m,\\ x_j\ge 0,& & i=1,2,\dots ,m,\end{array} \end{gathered}$$

• 目标函数是$x_j$ j个变量的线性组合
• 约束条件的形式：
  • 不等式
  • 等式
  • 整数（整数规划）

矩阵形式：
$$\begin{gathered} max(min)z=c^{T}x \\ s.t.\{\begin{array}{rl}Ax\le & (\ge ,=)b,\\ x\ge 0,& \end{array} \end{gathered}$$

• x 为$[x_1,x_2...]$ ,矩阵A的列数等于决策变量个数，每行为线性组合的系数:$Ax={\sum }_{j=1}^n{p}_j{x}_j$
• matlab中以矩阵为单位求解线性规划问题(不等号统一用$\le$)
• $c^T$为价值向量，$b$ 为资源向量（约束资源），$x$ 为决策向量，A 为决策向量的系数矩阵

线性规划模型标准型
$$\begin{gathered} maxz=c^{T}x \\ s.t.\{\begin{array}{r}Ax=b,\\ x\ge 0.\end{array} \end{gathered}$$

---

线性规划定理：

1. 满足全部约束条件的决策向量 $x\in R^n$ 称为可行解
2. 全部可行解的集合称为可行域
3. 使目标函数达到最大值的可行解为最优解

• 存在最优解时：
  • 无穷多个最优解（对应的目标函数值都相等）
  • 一个最优解
  • 有限个最优解
• 无最优解
  • 可行域为空集
  • 可行域非空，但目标函数无界

4. 求解线性规划问题：单纯型法

线性规划问题转换

1. 绝对值变换：
   $$\begin{gathered} min|x_1|+|x_2|+\dots |x_n| \\ s.t.Ax\le b \end{gathered}$$
   变量替换：$x_i=u_i-v_i,|x_i|=u_i+v_i$

   ​ $u_i=\frac{x_i+|x_i|}{2},v_i=\frac{|x_i|-x_i}{2}$

   表达式变换为：
   $$\begin{gathered} min\sum _{i=1}^n(u_i+v_i) \\ s.t.[A,-A]\left[\begin{array}{r}{u}_i\\ v_i\end{array}\right]\le b \\ {u}_i\ge 0,v_i\ge 0 \end{gathered}$$
   变量$u_i,v_i$都大于零则$|x|\ge 0$恒成立，将绝对值问题化为两个决策变量之和

2. 最大（小）值变换：
   m i n   { m a x   c i x i } s . t . A x ≤ b min~\{max~c_{i}x_i\}\\ s.t. Ax\le b min {max ci​xi​}s.t.Ax≤b
   令$x_{n+1}=max{c}_i{x}_i$,得到:
   $$\begin{gathered} min{x}_{n+1} \\ Ax\le b \\ {c}_i{x}_i-x_{n+1}\le 0,i=1,2,..n \end{gathered}$$

3. 引入松弛变量：将不等式转化为等式

---

求解线性规划问题

线性规划标准型：
$$\begin{gathered} minz=c^{T}x \\ s.t.\{\begin{array}{rl}Ax\le b,& \\ Aeqx=& beq,\\ Lb\le x\le & Ub.\end{array} \end{gathered}$$

• matlab

  • [x,fval]=linprog(f,A,b)
  • [x,fval]=linprog(f,A,b,Aeq,beq)
  • [x,fval]=linprog(f,A,b,Aeq,beq,lb,ub)
  • lb 和 ub 分别对应决策变量的下界，上界列向量
  • 如果没有等式约束:Aeq=[],beq=[] 设置为空矩阵

---

• python

  • scipy.optimize求解器（求解简单的线性规划和非线性规划）

  • 函数：linprog (c,A,b,Aeq,beq,bounds,method) ✅

    • c：目标向量
    • A：约束矩阵
    • b:资源向量
    • Aeq,beq ：等式矩阵
    • bounds：x的取值，默认为$0\le x\le +\infty$
    • method:算法（单纯型法‘simplex’，内点法…）

  • 代码：

    from scipy.optimize import linprog
    res=linprog(c,A,b,Aeq,beq,bounds,method)
    res.fun#目标函数最小值
    res.x#最有解
    

  • 对于max 用负号转化为最小值标准型

  范例：
  $$\begin{gathered} minz=-x_1+4x_2 \\ s.t.\{\begin{array}{r}-3x_1+x_2\le 6\\ x_1+2x_2\le 4\\ x_2\ge 3\end{array} \end{gathered}$$

  from scipy.optimize import linprog
  c = [-1, 4]; A = [[-3, 1], [1, 2]]
  b = [[6], [4]]
  bound=((None,None),(-3,None))# bound ((x1min,x1max),(x2min,x2max))
  res=linprog(c,A,b,None,None,bound)
  print("目标函数的最小值：",res.fun)
  print("最优解为：",res.x)
  
  #result:
  目标函数的最小值： -21.99999984082492
  最优解为： [ 9.99999989 -2.99999999]
  

  • 对于无约束的边界设置为None

  • 向量也可以传入数组代替列表（numpy）

  • 注意向量的维度

  • 构造bounds也可以像这样：

    Lb=[None,None]
    Ub=[-3,None]
    bound=tuple(zip(Lb,Ub))
    

  • cvxopt.solvers 模块求解线性规划：🍁
    $$\begin{gathered} minz=c^{T}x \\ s.t.\{\begin{array}{rl}Ax\le b,& \\ Aeqx=& beq,\end{array} \end{gathered}$$

    • 必须加载numpy 库
    • 构造矩阵 cvxpot.matrix()
    • 求解线性规划：cvxpot.solvers.lp(c,A,b) （x默认大于0）
    • 最优值cvxpot.solvers.lp(c,A,b)[‘x’]
    • 最优解cvxpot.solvers.lp(c,A,b) [‘primal objective’]

  • 凸优化库cvxpy 求解线性规划✅(最专业，可以求解线性，非线性，整数规划)

    • 构造决策变量矩阵
    • cvxpy.Variable((m,n))
    • cvxpy.Variable(n)
    • 构造目标函数 cvxpy.Minimize()
    • 构造约束条件（列表）
    • 构造模型 p=cvxpy.Problem()
    • 求解模型 p.solve(solver,verbose): 求解器(str)，是否展示过程（bool）

    求解参照：

    import cvxpy as cp
    import numpy as np
    import pandas as pd
    d1=pd.read_excel("Pdata5_6.xlsx",header=None)
    d2=d1.values; #获取DataFrame 的值
    c=d2[:-1,:-1]# 构造价值矩阵
    d=d2[-1,:-1].reshape(1,-1); e=d2[:-1,-1].reshape(-1,1)
    x=cp.Variable((6,8))#决策变量6*8
    obj=cp.Minimize(cp.sum(cp.multiply(c,x)))  #构造目标函数
    con=[cp.sum(x,axis=1,keepdims=True)<=e,
    cp.sum(x,axis=0,keepdims=True)==d,x>=0]  #构造约束条件
    prob=cp.Problem(obj,con)  #构造模型
    prob.solve(solver='GLPK_MI',verbose =True)#求解模型
    
    print("最优值为：",prob.value)
    print("最优解为：\n",x.value)
    

    灵敏度分析

  灵敏度分析：对系统因周围变化显示出来的敏感程度的分析

  分析如下：

  1. 如果系数 $a_{ij},bi(资源约束),c_j(价值系数)$ 发生改变，当前的最优方案会发生什么改变
  2. 将这些系数限制在什么范围内仍是最优解

  • 可以使用计算机迭代的方法（控制变量）寻求变换规律（$f(x)$）

    • 控制步长，控制变量循环求最优解
    • 绘制散点图，可视化结果