• Python数学建模—线性规划

            我是猿童学，本文是根据司守奎老师《数学建模算法与程序》的书本内容编写，使用其书中案例，书中的编程语言是MATLAB、Lingo，我将使用Python来解决问题。接下来的一个月我将学习使用Python来解决数学建模问题。将记录下学习中的笔记和部分案例。

     

    1、线性规划

    1.1 线性规划的实例与定义

    例   1 某机床厂生产甲、乙两种机床每台销售后的利润分别为 4000 元与3000元。 生产甲机床需用 A、B 机器加工，加工时间分别为每台 2 小时和1小时生产乙机床需用 A、B、C 三种机器加工，加工时间为每台各一小时。若每天可用于加工的机器时数分别为 A 机器 10 小时、B 机器 8 小时和C 机器 7 小时，问该厂应生产甲、乙机床各几台，才能使总利润最大？

    上述问题的数学模型：设该厂生产 x_{1}台甲机床和x_{2}乙机床时总利润最大，则x_{1}, x_{2}应满足：

    这里变量 x_{1}, x_{2} 称之为决策变量，（1）式被称为问题的目标函数，（2）中的几个不等式 是问题的约束条件，记为 s.t.(即 subject to)。由于上面的目标函数及约束条件均为线性 函数，故被称为线性规划问题。总之，线性规划问题是在一组线性约束条件的限制下，求一线性目标函数最大或最 小的问题。

    1.2 使用Python中的scipy库求解列 1 线性规划问题：

    线性规划的标准型为

    \begin{gathered} \min c^{T} x \\ \text { s.t. }\left\{\begin{array}{l} A x \leq b \\ A e q \cdot x=b e q \\ l b \leq x \leq u b \end{array}\right. \end{gathered}                                                  （3）

    其中，c是价值向量；A和b对应线性不等式约束；Aeq和beq对应线性等式约束；bounds对应公式中的lb和ub，决策向量的下界和上界。

    scipy.optimize.linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=None, method='simplex', callback=None, options=None)

    输入参数解析：

    • c：线性目标函数的系数
    • A_ub（可选参数）：不等式约束矩阵，\(A_{ub}\)的每一行指定x上的线性不等式约束的系数
    • b_ub（可选参数）：不等式约束向量，每个元素代表\(A_{ub}\)x的上限
    • A_eq（可选参数）：等式约束矩阵，\(A_{eq}\)的每一行指定x xx上的线性等式约束的系数
    • b_eq（可选参数）：等式约束向量，\(A_{eq}\)的每个元素必须等于\(b_{eq }\)的对应元素
    • bounds（可选参数）：定义决策变量x的最小值和最大值
      • 数据类型： (min, max)序列对
      • None：使用None表示没有界限，默认情况下，界限为（0，None）（所有决策变量均为非负数）
      • 如果提供一个元组(min, max)，则最小值和最大值将用作所有决策变量的界限。

    输出参数解析：

    • x：在满足约束的情况下将目标函数最小化的决策变量的值
    • fun：目标函数的最优值（c^{^{T}}x）
    • slack：不等式约束的松弛值（名义上为正值）b_{ub}-A_{ub}x
    • con：等式约束的残差（名义上为零）\(b_{eq}-A_{eq}x\)
    • success：当算法成功找到最佳解决方案时为 True
    • status：表示算法退出状态的整数
      • 0 : 优化成功终止
      • 1 : 达到了迭代限制
      • 2 : 问题似乎不可行
      • 3 : 问题似乎是不收敛
      • 4 : 遇到数值困难
    • nit：在所有阶段中执行的迭代总数
    • message：算法退出状态的字符串描述符

    使用scipy库来解决例1：

    1. from scipy.optimize import linprog
    2. c = [-4, -3]   #目标函数的系数
    3. A= [[2, 1], [1, 1],[0,1]]  #<=不等式左侧未知数系数矩阵
    4. b = [10,8,7]    # <=不等式右侧常数矩阵
    5. #A_eq = [[]]   #等式左侧未知数系数矩阵
    6. #B_eq = []     #等式右侧常数矩阵
    7. x1 = [0, None]  #未知数取值范围
    8. x2 = [0, None]  #未知数取值范围
    9. res =linprog(c, A, b, bounds=(x1, x2))#默认求解最小值，求解最大值使用-c并取结果相反数
    10. print(res)
    11.

    输出结果：

    1.      con: array([], dtype=float64)
    2.      fun: -25.999999999841208
    3.  message: 'Optimization terminated successfully.'
    4.      nit: 5
    5.    slack: array([8.02629074e-11, 3.92663679e-11, 1.00000000e+00])
    6.   status: 0
    7.  success: True
    8.        x: array([2., 6.])

    fun为目标函数的最优值，slack为松弛变量，status表示优化结果状态，x为最优解。

    此模型的求解结果为：当x1=2，x2=6时，函数取得最优值25.999999999841208。


    例 2 求解下列线性规划问题：

    \(\begin{gathered} \max z=2x_{1}+3x_{2}-5x_{3} \\ \text { s.t. }\left\{\begin{array}{l}x _{1}+x _{2}+x _{3}=7 \\2x _{1}-5x _{2}+x _{3}\geq 10 \\x _{1}+3x _{2}+x _{3}\leq 12 \\ x _{1},x _{2},x _{3}\geq 0 \end{array}\right. \end{gathered}\)

    注意要先标准化！

    这里使用另一种比较规范的写法：

    1. from scipy import optimize
    2. import numpy as np
    3. c = np.array([2,3,-5])
    4. A = np.array([[-2,5,-1],[1,3,1]])
    5. B = np.array([-10,12])
    6. Aeq = np.array([[1,1,1]])
    7. Beq = np.array([7])
    8. x1 = [0, None]  #未知数取值范围
    9. x2 = [0, None]  #未知数取值范围
    10. x3 = [0, None]  #未知数取值范围
    11. res = optimize.linprog(-c,A,B,Aeq,Beq,bounds=(x1, x2, x3))
    12. print(res)

    输出结果： 

    1.      con: array([1.80713489e-09])
    2.      fun: -14.571428565645054
    3.  message: 'Optimization terminated successfully.'
    4.      nit: 5
    5.    slack: array([-2.24579466e-10,  3.85714286e+00])
    6.   status: 0
    7.  success: True
    8.        x: array([6.42857143e+00, 5.71428571e-01, 2.35900788e-10])

    此模型的求解结果为：当x1=6.4，x2=5.7，x3=2.3时，函数取得最优值14.571428565645054

    例3 求解线性规划问题：

    \(\begin{gathered} \min z=2x_{1}+3x_{2} + x_{3} \\ \text { s.t. }\left\{\begin{array}{l}x _{1}+4x _{2}+2x _{3} \geq 8 \\ 3x _{1} + 2x _{2} \geq 6 \\ x _{1},x _{2},x _{3}\geq 0 \end{array}\right. \end{gathered}\)

    使用Python编写：

    1. from scipy import optimize
    2. import numpy as np
    3. c = np.array([2,3,1])
    4. A = np.array([[1,4,2],[3,2,0]])
    5. B = np.array([8,6])
    6. x1 = [0, None]  #未知数取值范围
    7. x2 = [0, None]  #未知数取值范围
    8. x3 = [0, None]  #未知数取值范围
    9. res = optimize.linprog(c,-A,-B,bounds=(x1, x2, x3))
    10. print(res)

    输出结果：

    1.      con: array([], dtype=float64)
    2.      fun: 6.999999994872993
    3.  message: 'Optimization terminated successfully.'
    4.      nit: 3
    5.    slack: array([ 3.85261245e-09, -1.41066261e-08])
    6.   status: 0
    7.  success: True
    8.        x: array([1.17949641, 1.23075538, 0.94874104])

    此模型的求解结果为：当x1=1.1，x2=1.2，x3=0.9时，函数取得最优值6.999999994872993

    2、运输问题

    2.1产销平衡的数学建模过程：

            例 6 某商品有\(m\)个产地、\(n\)个销地，各产地的产量分别为\(a_{1},... ,a_{m}\)，各销地的需求量分别为\(b_{1},... ,b_{n}\) 。若该商品由 \(i\) 产地运到 \(j\) 销地的单位运价为\(c_{ij}\)，问应该如何调运才能使总运费最省？

            解：引入变量 \(x_{ij}\) ，其取值为由 \(i\) 产地运往 \(j\) 销地的该商品数量，数学模型为

    \(\begin{array}{l} \min \sum_{i=1}^{m} \sum_{j=1}^{n} c_{i j} x_{i j} \\ \text { s.t. }\left\{\begin{array}{l} \sum_{j=1}^{n} x_{i j}=a_{i}, \quad i=1, \cdots, m \\ \sum_{i=1}^{m} x_{i j}=b_{j}, \quad j=1,2, \cdots, n \\ x_{i j} \geq 0 \end{array}\right. \end{array}\)

    显然是一个线性规划问题，当然可以用单纯形法求解。

    对产销平衡的运输问题，由于有以下关系式存在：

    \(\sum_{j=1}^{n} b_{j}=\sum_{i=1}^{m}\left(\sum_{j=1}^{n} x_{i j}\right)=\sum_{j=1}^{n}\left(\sum_{i=1}^{m} x_{i j}\right)=\sum_{i=1}^{m} a_{i}\)

    其约束条件的系数矩阵相当特殊，可用比较简单的计算方法，习惯上称为表上作业法（由 康托洛维奇和希奇柯克两人独立地提出，简称康—希表上作业法），下面使用例子来讲解。

    2.2 产销平衡的运输问题

    1.问题描述

    某公司下属有甲、乙、丙三个工厂，分别向A、B、C、D四个销地提供产品，产量、需求量及工厂到销售地的运价（单位：元/吨）如下表所示，求使费用最少的最佳运输方案。

    工厂 \ 销地ABCD产量（吨）
    甲8610918
    乙91213118
    丙141216519
    销量（吨）161571755

    2.问题解析

    总产量=18+18+19=55

    总销量=16+15+7+17=55

    产量等于销量，即这是产销平衡的运输问题。直接采用表上作业法进行求解。

    3.使用python求解：

    1. # 运输问题求解：使用Vogel逼近法寻找初始基本可行解
    2. import numpy as np
    3. import copy
    4. import pandas as pd
    5.  
    6.  
    7. def main():
    8.     # mat = pd.read_csv('表上作业法求解运输问题.csv', header=None).values
    9.     # mat = pd.read_excel('表上作业法求解运输问题.xlsx', header=None).values
    10.     a = np.array([18, 18, 19])  # 产量
    11.     b = np.array([16, 15, 7, 17])  # 销量
    12.     c = np.array([[8, 6, 10, 9], [9, 12, 13, 7], [14, 12, 16, 5]])
    13.     # [c, x] = TP_vogel(mat)
    14.     [c,x]=TP_vogel([c,a,b])
    15.  
    16.  
    17. def TP_split_matrix(mat):  # 运输分割矩阵
    18.     c = mat[:-1, :-1]
    19.     a = mat[:-1, -1]
    20.     b = mat[-1, :-1]
    21.     return (c, a, b)
    22.  
    23.  
    24. def TP_vogel(var):  # Vogel法代码，变量var可以是以numpy.ndarray保存的运输表，或以tuple或list保存的(成本矩阵,供给向量,需求向量)
    25.     import numpy
    26.     typevar1 = type(var) == numpy.ndarray
    27.     typevar2 = type(var) == tuple
    28.     typevar3 = type(var) == list
    29.     if typevar1 == False and typevar2 == False and typevar3 == False:
    30.         print('>>>非法变量<<<')
    31.         (cost, x) = (None, None)
    32.     else:
    33.         if typevar1 == True:
    34.             [c, a, b] = TP_split_matrix(var)
    35.         elif typevar2 == True or typevar3 == True:
    36.             [c, a, b] = var
    37.         cost = copy.deepcopy(c)
    38.         x = np.zeros(c.shape)
    39.         M = pow(10, 9)
    40.         for factor in c.reshape(1, -1)[0]:
    41.             p = 1
    42.             while int(factor) != M:
    43.                 if np.all(c == M):
    44.                     break
    45.                 else:
    46.                     print('第1次迭代！')
    47.                     print('c:\n', c)
    48.                     # 获取行/列最小值数组
    49.                     row_mini1 = []
    50.                     row_mini2 = []
    51.                     for row in range(c.shape[0]):
    52.                         Row = list(c[row, :])
    53.                         row_min = min(Row)
    54.                         row_mini1.append(row_min)
    55.                         Row.remove(row_min)
    56.                         row_2nd_min = min(Row)
    57.                         row_mini2.append(row_2nd_min)
    58.                     # print(row_mini1,'\n',row_mini2)
    59.                     r_pun = [row_mini2[i] - row_mini1[i] for i in range(len(row_mini1))]
    60.                     print('行罚数：', r_pun)
    61.                     # 计算列罚数
    62.                     col_mini1 = []
    63.                     col_mini2 = []
    64.                     for col in range(c.shape[1]):
    65.                         Col = list(c[:, col])
    66.                         col_min = min(Col)
    67.                         col_mini1.append(col_min)
    68.                         Col.remove(col_min)
    69.                         col_2nd_min = min(Col)
    70.                         col_mini2.append(col_2nd_min)
    71.                     c_pun = [col_mini2[i] - col_mini1[i] for i in range(len(col_mini1))]
    72.                     print('列罚数：', c_pun)
    73.                     pun = copy.deepcopy(r_pun)
    74.                     pun.extend(c_pun)
    75.                     print('罚数向量：', pun)
    76.                     max_pun = max(pun)
    77.                     max_pun_index = pun.index(max(pun))
    78.                     max_pun_num = max_pun_index + 1
    79.                     print('最大罚数：', max_pun, '元素序号：', max_pun_num)
    80.                     if max_pun_num <= len(r_pun):
    81.                         row_num = max_pun_num
    82.                         print('对第', row_num, '行进行操作：')
    83.                         row_index = row_num - 1
    84.                         catch_row = c[row_index, :]
    85.                         print(catch_row)
    86.                         min_cost_colindex = int(np.argwhere(catch_row == min(catch_row)))
    87.                         print('最小成本所在列索引：', min_cost_colindex)
    88.                         if a[row_index] <= b[min_cost_colindex]:
    89.                             x[row_index, min_cost_colindex] = a[row_index]
    90.                             c1 = copy.deepcopy(c)
    91.                             c1[row_index, :] = [M] * c1.shape[1]
    92.                             b[min_cost_colindex] -= a[row_index]
    93.                             a[row_index] -= a[row_index]
    94.                         else:
    95.                             x[row_index, min_cost_colindex] = b[min_cost_colindex]
    96.                             c1 = copy.deepcopy(c)
    97.                             c1[:, min_cost_colindex] = [M] * c1.shape[0]
    98.                             a[row_index] -= b[min_cost_colindex]
    99.                             b[min_cost_colindex] -= b[min_cost_colindex]
    100.                     else:
    101.                         col_num = max_pun_num - len(r_pun)
    102.                         col_index = col_num - 1
    103.                         print('对第', col_num, '列进行操作：')
    104.                         catch_col = c[:, col_index]
    105.                         print(catch_col)
    106.                         # 寻找最大罚数所在行/列的最小成本系数
    107.                         min_cost_rowindex = int(np.argwhere(catch_col == min(catch_col)))
    108.                         print('最小成本所在行索引：', min_cost_rowindex)
    109.                         # 计算将该位置应填入x矩阵的数值（a,b中较小值）
    110.                         if a[min_cost_rowindex] <= b[col_index]:
    111.                             x[min_cost_rowindex, col_index] = a[min_cost_rowindex]
    112.                             c1 = copy.deepcopy(c)
    113.                             c1[min_cost_rowindex, :] = [M] * c1.shape[1]
    114.                             b[col_index] -= a[min_cost_rowindex]
    115.                             a[min_cost_rowindex] -= a[min_cost_rowindex]
    116.                         else:
    117.                             x[min_cost_rowindex, col_index] = b[col_index]
    118.                             # 填入后删除已满足/耗尽资源系数的行/列，得到剩余的成本矩阵，并改写资源系数
    119.                             c1 = copy.deepcopy(c)
    120.                             c1[:, col_index] = [M] * c1.shape[0]
    121.                             a[min_cost_rowindex] -= b[col_index]
    122.                             b[col_index] -= b[col_index]
    123.                     c = c1
    124.                     print('本次迭代后的x矩阵：\n', x)
    125.                     print('a:', a)
    126.                     print('b:', b)
    127.                     print('c:\n', c)
    128.                 if np.all(c == M):
    129.                     print('【迭代完成】')
    130.                     print('-' * 60)
    131.  
    132.                 else:
    133.                     print('【迭代未完成】')
    134.                     print('-' * 60)
    135.  
    136.                     p += 1
    137.                     print(f"第{p}次迭代！")
    138.         total_cost = np.sum(np.multiply(x, cost))
    139.         if np.all(a == 0):
    140.             if np.all(b == 0):
    141.                 print('>>>供求平衡<<<')
    142.             else:
    143.                 print('>>>供不应求，需求方有余量<<<')
    144.         elif np.all(b == 0):
    145.             print('>>>供大于求，供给方有余量<<<')
    146.         else:
    147.             print('>>>无法找到初始基可行解<<<')
    148.         print('>>>初始基本可行解x*：\n', x)
    149.         print('>>>当前总成本：', total_cost)
    150.         [m, n] = x.shape
    151.         varnum = np.array(np.nonzero(x)).shape[1]
    152.         if varnum != m + n - 1:
    153.             print('【注意：问题含有退化解】')
    154.     return (cost, x)
    155.  
    156.  
    157. if __name__ == '__main__':
    158.     main()

     输出结果：

    1. D:\Anaconda\python.exe D:/PyCharm/数学建模/线性规划/www.py
    2. 第1次迭代！
    3. c:
    4.  [[ 8  6 10  9]
    5.  [ 9 12 13  7]
    6.  [14 12 16  5]]
    7. 行罚数： [2, 2, 7]
    8. 列罚数： [1, 6, 3, 2]
    9. 罚数向量： [2, 2, 7, 1, 6, 3, 2]
    10. 最大罚数： 7 元素序号： 3
    11. 对第 3 行进行操作：
    12. [14 12 16  5]
    13. 最小成本所在列索引： 3
    14. 本次迭代后的x矩阵：
    15.  [[ 0.  0.  0.  0.]
    16.  [ 0.  0.  0.  0.]
    17.  [ 0.  0.  0. 17.]]
    18. a: [18 18  2]
    19. b: [16 15  7  0]
    20. c:
    21.  [[         8          6         10 1000000000]
    22.  [         9         12         13 1000000000]
    23.  [        14         12         16 1000000000]]
    24. 【迭代未完成】
    25. ------------------------------------------------------------
    26. 第2次迭代！
    27. 第1次迭代！
    28. c:
    29.  [[         8          6         10 1000000000]
    30.  [         9         12         13 1000000000]
    31.  [        14         12         16 1000000000]]
    32. 行罚数： [2, 3, 2]
    33. 列罚数： [1, 6, 3, 0]
    34. 罚数向量： [2, 3, 2, 1, 6, 3, 0]
    35. 最大罚数： 6 元素序号： 5
    36. 对第 2 列进行操作：
    37. [ 6 12 12]
    38. 最小成本所在行索引： 0
    39. 本次迭代后的x矩阵：
    40.  [[ 0. 15.  0.  0.]
    41.  [ 0.  0.  0.  0.]
    42.  [ 0.  0.  0. 17.]]
    43. a: [ 3 18  2]
    44. b: [16  0  7  0]
    45. c:
    46.  [[         8 1000000000         10 1000000000]
    47.  [         9 1000000000         13 1000000000]
    48.  [        14 1000000000         16 1000000000]]
    49. 【迭代未完成】
    50. ------------------------------------------------------------
    51. 第3次迭代！
    52. 第1次迭代！
    53. c:
    54.  [[         8 1000000000         10 1000000000]
    55.  [         9 1000000000         13 1000000000]
    56.  [        14 1000000000         16 1000000000]]
    57. 行罚数： [2, 4, 2]
    58. 列罚数： [1, 0, 3, 0]
    59. 罚数向量： [2, 4, 2, 1, 0, 3, 0]
    60. 最大罚数： 4 元素序号： 2
    61. 对第 2 行进行操作：
    62. [         9 1000000000         13 1000000000]
    63. 最小成本所在列索引： 0
    64. 本次迭代后的x矩阵：
    65.  [[ 0. 15.  0.  0.]
    66.  [16.  0.  0.  0.]
    67.  [ 0.  0.  0. 17.]]
    68. a: [3 2 2]
    69. b: [0 0 7 0]
    70. c:
    71.  [[1000000000 1000000000         10 1000000000]
    72.  [1000000000 1000000000         13 1000000000]
    73.  [1000000000 1000000000         16 1000000000]]
    74. 【迭代未完成】
    75. ------------------------------------------------------------
    76. 第4次迭代！
    77. 第1次迭代！
    78. c:
    79.  [[1000000000 1000000000         10 1000000000]
    80.  [1000000000 1000000000         13 1000000000]
    81.  [1000000000 1000000000         16 1000000000]]
    82. 行罚数： [999999990, 999999987, 999999984]
    83. 列罚数： [0, 0, 3, 0]
    84. 罚数向量： [999999990, 999999987, 999999984, 0, 0, 3, 0]
    85. 最大罚数： 999999990 元素序号： 1
    86. 对第 1 行进行操作：
    87. [1000000000 1000000000         10 1000000000]
    88. 最小成本所在列索引： 2
    89. 本次迭代后的x矩阵：
    90.  [[ 0. 15.  3.  0.]
    91.  [16.  0.  0.  0.]
    92.  [ 0.  0.  0. 17.]]
    93. a: [0 2 2]
    94. b: [0 0 4 0]
    95. c:
    96.  [[1000000000 1000000000 1000000000 1000000000]
    97.  [1000000000 1000000000         13 1000000000]
    98.  [1000000000 1000000000         16 1000000000]]
    99. 【迭代未完成】
    100. ------------------------------------------------------------
    101. 第5次迭代！
    102. 第1次迭代！
    103. c:
    104.  [[1000000000 1000000000 1000000000 1000000000]
    105.  [1000000000 1000000000         13 1000000000]
    106.  [1000000000 1000000000         16 1000000000]]
    107. 行罚数： [0, 999999987, 999999984]
    108. 列罚数： [0, 0, 3, 0]
    109. 罚数向量： [0, 999999987, 999999984, 0, 0, 3, 0]
    110. 最大罚数： 999999987 元素序号： 2
    111. 对第 2 行进行操作：
    112. [1000000000 1000000000         13 1000000000]
    113. 最小成本所在列索引： 2
    114. 本次迭代后的x矩阵：
    115.  [[ 0. 15.  3.  0.]
    116.  [16.  0.  2.  0.]
    117.  [ 0.  0.  0. 17.]]
    118. a: [0 0 2]
    119. b: [0 0 2 0]
    120. c:
    121.  [[1000000000 1000000000 1000000000 1000000000]
    122.  [1000000000 1000000000 1000000000 1000000000]
    123.  [1000000000 1000000000         16 1000000000]]
    124. 【迭代未完成】
    125. ------------------------------------------------------------
    126. 第6次迭代！
    127. 第1次迭代！
    128. c:
    129.  [[1000000000 1000000000 1000000000 1000000000]
    130.  [1000000000 1000000000 1000000000 1000000000]
    131.  [1000000000 1000000000         16 1000000000]]
    132. 行罚数： [0, 0, 999999984]
    133. 列罚数： [0, 0, 999999984, 0]
    134. 罚数向量： [0, 0, 999999984, 0, 0, 999999984, 0]
    135. 最大罚数： 999999984 元素序号： 3
    136. 对第 3 行进行操作：
    137. [1000000000 1000000000         16 1000000000]
    138. 最小成本所在列索引： 2
    139. 本次迭代后的x矩阵：
    140.  [[ 0. 15.  3.  0.]
    141.  [16.  0.  2.  0.]
    142.  [ 0.  0.  2. 17.]]
    143. a: [0 0 0]
    144. b: [0 0 0 0]
    145. c:
    146.  [[1000000000 1000000000 1000000000 1000000000]
    147.  [1000000000 1000000000 1000000000 1000000000]
    148.  [1000000000 1000000000 1000000000 1000000000]]
    149. 【迭代完成】
    150. ------------------------------------------------------------
    151. >>>供求平衡<<<
    152. >>>初始基本可行解x*：
    153.  [[ 0. 15.  3.  0.]
    154.  [16.  0.  2.  0.]
    155.  [ 0.  0.  2. 17.]]
    156. >>>当前总成本： 407.0
    157.  
    158. 进程已结束，退出代码为 0

    上面使用的方法是表上作业法，属下所示：

     最终最优的总成本： 407.0

    2.3 产大于销的运输问题

    1.问题描述

    有三个牧业基地向4个城市提供鲜奶，4个城市每日的鲜奶需求量、3个基地的每日鲜奶供应量以及运送每千升鲜奶的费用如下表所示，试确定最经济的鲜奶运输方案。

    城市\基地ABcD供应量
    甲324530
    乙235340
    丙142450
    需求量16302430

    2.问题解析

    总供应量=30+40+50=120

    总需求量=16+30+24+30=100

    供过于求，即产量大于销量，这是一个产大于销的运输问题。要将问题转化为产销平衡的运输问题，需进行以下几个方面的调整。

    （1）增加一个虚拟销地E，使其总需求量为120-100=20。

    （2）由于销地是虚拟的，实际上是产量过剩的物资在产地就地储存，所以不会产生实际的运输，即不会产生运费。

    因此新的供需量及单位运价表如下所示。

    城市\基地ABcDE供应量
    甲3245030
    乙2353040
    丙1424050
    需求量1630243020120

    3.问题求解

    使用python来实现： 

    这里也可以使用（最小元素法）来解决:

    1. '----------------------------------------最小元素法------------------------------------' \
    2. '最小元素法是指有限满足单位运价最小的供销服务'
    3. import numpy as np
    4. import copy
    5.  
    6. supplyNum = 3  # 供应商数量
    7. demandNum = 5  # 需求商数量
    8. A = np.array([30,40,50])  # 产量
    9. B = np.array([16,30,24,30,20])  # 销量
    10. C = np.array([[3,2,4,5,0], [2,3,5,3,0], [1,4,2,4,0]])  # 成本
    11. X = np.zeros((supplyNum, demandNum))  # 初始解
    12. c = copy.deepcopy(C)
    13. maxNumber = 10000  # 极大数
    14.  
    15.  
    16. def pivot(A, B, c, X):
    17.     index = np.where(c == np.min(c))  # 寻找最小值的索引
    18.     minIndex = (index[0][0], index[1][0])  # 确定最小值的索引
    19.     # 确定应该优先分配的量，并且改变价格
    20.     if A[index[0][0]] > B[index[1][0]]:
    21.         X[minIndex] = B[index[1][0]]
    22.         c[:, index[1][0]] = maxNumber  # 改变价格
    23.         A[index[0][0]] = A[index[0][0]] - B[index[1][0]]
    24.         B[index[1][0]] = 0  # 改变销量
    25.     else:
    26.         X[minIndex] = A[index[0][0]]
    27.         c[index[0][0], :] = maxNumber
    28.         B[index[1][0]] = B[index[1][0]] - A[index[0][0]]
    29.         A[index[0][0]] = 0
    30.     return A, B, c, X
    31.  
    32.  
    33. # 最小元素法
    34. def minimumElementMethod(A, B, c, X):
    35.     while (c < maxNumber).any():
    36.         A, B, c, X = pivot(A, B, c, X)
    37.     return X
    38.  
    39.  
    40. solutionMin = minimumElementMethod(A, B, c, X)
    41. print('最小元素法得到的初始解为：')
    42. print(solutionMin)
    43. print('最小元素法得到的总运费为:', ((solutionMin * C).sum()))

    输出结果：

    1. 最小元素法得到的初始解为：
    2. [[ 0. 10.  0.  0. 20.]
    3.  [ 0. 20.  0. 20.  0.]
    4.  [16.  0. 24. 10.  0.]]
    5. 最小元素法得到的总运费为: 244.0

    2.4 产小于销的运输问题

    1.问题描述

    某三个煤炭厂供应4个地区，假定等量的煤炭在这些地区使用效果相同，已知各煤炭厂年产量，各地区的需要量及从各煤炭厂到各地区的单位运价表如下所示，试决定最优的调运方案。

    销地\产地B1B2B3B4产量
    A15310490
    A2269640
    A314105770
    销量305010040

    2.问题解析

    总供应量=30+40+50=120

    总需求量=16+30+24+30=100

    供过于求，即产量大于销量，这是一个产大于销的运输问题。要将问题转化为产销平衡的运输问题，需进行以下几个方面的调整。

    （1）增加一个虚拟销地E，使其总需求量为120-100=20。

    （2）由于销地是虚拟的，实际上是产量过剩的物资在产地就地储存，所以不会产生实际的运输，即不会产生运费。

    因此新的供需量及单位运价表如下所示。

    销地\产地BBB3B4产量
    A,5310490
    A2269640
    A314105770
    A4000020
    销量305010040220

    3.问题求解

    使用python求解：

    这里可以使用（Vogel逼近法）寻找初始基本可行解

    1. '-----------------------------------Vogel法---------------------------------'
    2. supplyNum = 4  # 供应商数量
    3. demandNum = 4  # 需求商数量
    4. A = np.array([90,40,70,20])  # 产量
    5. B = np.array([30,50,100,40])  # 销量
    6. C = np.array([[5,3,10,4], [2,6,9,6], [14,10,5,7],[0,0,0,0]])  # 成本
    7. X = np.zeros((supplyNum, demandNum))  # 初始解
    8. c = copy.deepcopy(C)
    9. numPunish = [0] * (supplyNum + demandNum)  # 整体罚数
    10.  
    11.  
    12. def pivot(X, A, B, C):
    13.     # 求解罚数，其中列罚数排在行罚数后面
    14.     for i in range(supplyNum):
    15.         sortList = np.sort(C[i, :])
    16.         numPunish[i] = sortList[1] - sortList[0]
    17.     for i in range(demandNum):
    18.         sortList = np.sort(C[:, i])
    19.         numPunish[i + supplyNum] = sortList[1] - sortList[0]
    20.  
    21.     maxIndex = np.argmax(numPunish)  # 寻找最大罚数的索引
    22.  
    23.     # 若最大的罚数属于行罚数
    24.     if maxIndex < supplyNum:
    25.         minIndex = np.argmin(C[maxIndex, :])
    26.         index = (maxIndex, minIndex)  # 得到最大罚数的一行中最小的运价
    27.         # 若产量大于销量
    28.         if A[maxIndex] > B[minIndex]:
    29.             X[index] = B[minIndex]  # 更新解
    30.             C[:, minIndex] = maxNumber  # 将已经满足的一列中的运价增大替代删除操作
    31.             A[maxIndex] -= B[minIndex]  # 更新剩余产量
    32.             B[minIndex] = 0  # 更新已经满足的销量
    33.         # 若销量大于产量
    34.         else:
    35.             X[index] = A[maxIndex]
    36.             C[maxIndex, :] = maxNumber
    37.             B[minIndex] -= A[maxIndex]  # 更新销量
    38.             A[maxIndex] = 0
    39.  
    40.     # 若最大的罚数为列罚数
    41.     else:
    42.         minIndex = np.argmin(C[:, maxIndex - supplyNum])  # 这时maxIndex-supplyNum是罚数最大的列
    43.         index = (minIndex, maxIndex - supplyNum)
    44.         if A[minIndex] > B[maxIndex - supplyNum]:  # 这里是产量大于销量，因此有限满足销量
    45.             X[index] = B[maxIndex - supplyNum]
    46.             C[:, maxIndex - supplyNum] = maxNumber
    47.             A[minIndex] -= B[maxIndex - supplyNum]
    48.             B[maxIndex - supplyNum] = 0
    49.         else:
    50.             X[index] = A[minIndex]
    51.             C[minIndex, :] = maxNumber
    52.             B[maxIndex - supplyNum] -= A[minIndex]
    53.             A[minIndex] = 0
    54.     return X, A, B, C
    55.  
    56.  
    57. # 沃格尔法得到初始解
    58. def Vogel(A, B, C):
    59.     X = np.zeros((supplyNum, demandNum))  # 初始解
    60.     numPunish = [0] * (supplyNum + demandNum)  # 整体罚数
    61.     # 迭代，直到所有的产量和销量全部满足
    62.     while (A != 0).any() and (B != 0).any():
    63.         X, A, B, C = pivot(X, A, B, C)
    64.     return X
    65.  
    66.  
    67. # 上面对于条件的判断，还需要对销量和需求量进行改变
    68. solutionVogel = Vogel(A, B, c)
    69. print('Vogel法得到的初始解为：')
    70. print(solutionVogel)
    71. print('Vogel法得到的总运费为:', (C * solutionVogel).sum())

     输出结果：

    1. Vogel法得到的初始解为：
    2. [[ 0. 50.  0. 40.]
    3.  [30.  0. 10.  0.]
    4.  [ 0.  0. 70.  0.]
    5.  [ 0.  0. 20.  0.]]
    6. Vogel法得到的总运费为: 810.0

    上面三种产销问题是常见的情况，除了产销问题还有弹性需求的运输问题、中间转运的运输问题。

    这里就不罗嗦了，自己查资料噢！

    3、指派问题

    3.1 指派问题的数学模型

            例 7 拟分配 \(n\) 人去干 \(n\)项工作，每人干且仅干一项工作，若分配第 \(i\) 人去干第 \(j\) 项工作，需花费 \(c_{ij}\) 单位时间，问应如何分配工作才能使工人花费的总时间最少？

            容易看出，要给出一个指派问题的实例，只需给出矩阵C=(c_{ij}) ， C被称为指派 问题的系数矩阵。 引入变量 x_{ij} ，若分配 i 干 j 工作，则取 x_{ij} = 1，否则取 x_{ij}= 0 。上述指派问题的数学模型为：

    \begin{aligned} &\min \sum_{i=1}^{n} \sum_{j=1}^{n} c_{i j} x_{i j} \\ &\text { s.t. } \sum_{j=1}^{n} x_{i j}=1 \\ &\sum_{i=1}^{n} x_{i j}=1 \\ &x_{i j}=0 \text { or } 1 \end{aligned}

            述指派问题的可行解可以用一个矩阵表示，其每行每列均有且只有一个元素为 1，其余元素均为 0；可以用\(1,...,n\)中的一个置换表示。

             问题中的变量只能取 0 或 1，从而是一个 0-1 规划问题。一般的 0-1 规划问题求解 极为困难。但指派问题并不难解，其约束方程组的系数矩阵十分特殊（被称为全单位模矩阵，其各阶非零子式均为 ±1），其非负可行解的分量只能取 0 或 1，故约束 x_{ij} = 0或1 可改写为 x_{ij} ≥ 0而不改变其解。此时，指派问题被转化为一个特殊的运输问题，其中 \(m=n\) ，\(a_{i}=b_{j}=1\)。

    3.2 求解指派问题的匈牙利算法

    算法主要依据以下事实：如果系数矩阵 C=(c_{ij}) 一行（或一列）中每 一元素都加上或减去同一个数，得到一个新矩阵 \(B=(b_{ij})\)，则以C 或 B 为系数矩阵的 指派问题具有相同的最优指派。

    例 8 求解指派问题，其系数矩阵为

    \(C=\left[\begin{array}{llll} 16 & 15 & 19 & 22 \\ 17 & 21 & 19 & 18 \\ 24 & 22 & 18 & 17 \\ 17 & 19 & 22 & 16 \end{array}\right]\)

    作业法解：

     将第一行元素减去此行中的最小元素 15，同样，第二行元素减去 17，第三行 元素减去 17，最后一行的元素减去 16，得：

    \(B_{1}=\left[\begin{array}{llll} 1 & 0 & 4 & 7 \\ 0 & 4 & 2 & 1 \\ 7 & 5 & 1 & 0 \\ 1 & 3 & 6 & 0 \end{array}\right]\)

    再将第 3 列元素各减去 1，得

    \(B_{2}=\left[\begin{array}{cccc} 1 & 0^{*} & 3 & 7 \\ 0^{*} & 4 & 1 & 1 \\ 7 & 5 & 0^{*} & 0 \\ 1 & 3 & 5 & 0^{*} \end{array}\right]\)

    以\(B_{2}\) 为系数矩阵的指派问题有最优指派

    \(\left(\begin{array}{llll} 1 & 2 & 3 & 4 \\ 2 & 1 & 3 & 4 \end{array}\right)\)

    由等价性，它也是C 的最优指派。

    python解：

    1. from scipy.optimize import linear_sum_assignment
    2. import numpy as np
    3. cost = np.array([[16,15,19,22], [17,21,19,18], [24,22,18,17],[17,19,22,16]])
    4. row_ind, col_ind = linear_sum_assignment(cost)
    5. print('row_ind:', row_ind)  # 开销矩阵对应的行索引
    6. print('col_ind:', col_ind)  # 对应行索引的最优指派的列索引
    7. print('cost:', cost[row_ind, col_ind])  # 提取每个行索引的最优指派列索引所在的元素，形成数组
    8. print('cost_sum:', cost[row_ind, col_ind].sum())  # 最小开销

     输出结果：

    1. row_ind: [0 1 2 3]
    2. col_ind: [1 0 2 3]
    3. cost: [15 17 18 16]
    4. cost_sum: 66

    例 9 求解系数矩阵C 的指派问题

    C=\left[\begin{array}{lllll} 12 & 7 & 9 & 7 & 9 \\ 8 & 9 & 6 & 6 & 6 \\ 7 & 17 & 12 & 14 & 12 \\ 15 & 14 & 6 & 6 & 10 \\ 4 & 10 & 7 & 10 & 6 \end{array}\right]

    python实现：

    1. from scipy.optimize import linear_sum_assignment
    2. import numpy as np
    3. cost = np.array([[12, 7, 9 ,7, 9], 
    4.                  [8, 9, 6, 6, 6], 
    5.                  [7, 17, 12, 14, 12],
    6.                  [15,14, 6 ,6, 10],
    7.                  [4,10 ,7 ,10 ,6]])
    8. row_ind, col_ind = linear_sum_assignment(cost)
    9. print('row_ind:', row_ind)  # 开销矩阵对应的行索引
    10. print('col_ind:', col_ind)  # 对应行索引的最优指派的列索引
    11. print('cost:', cost[row_ind, col_ind])  # 提取每个行索引的最优指派列索引所在的元素，形成数组
    12. print('cost_sum:', cost[row_ind, col_ind].sum())  # 最小开销

    输出结果：

    1. row_ind: [0 1 2 3 4]
    2. col_ind: [1 2 0 3 4]
    3. cost: [7 6 7 6 6]
    4. cost_sum: 32


    4、对偶理论与灵敏度分析

    4.1 原始问题和对偶问题

    考虑下列一对线性规划模型：

    \(\max c^{T} x \quad \text { s.t. } \quad A x \leq b, x \geq 0\)                                                                                         （P）

    和

    \(\min b^{T} y \quad \text { s.t. } \quad A^{T} y \geq c, y \geq 0\)                                                                                        （D）

    称（P）为原始问题，（D）为它的对偶问题。

    不太严谨地说，对偶问题可被看作是原始问题的“行列转置”：

    （1） 原始问题中的第 j 列系数与其对偶问题中的第 j 行的系数相同；

    （2） 原始目标函数的各个系数行与其对偶问题右侧的各常数列相同；

    （3） 原始问题右侧的各常数列与其对偶目标函数的各个系数行相同；

    （4） 在这一对问题中，不等式方向和优化方向相反。

    考虑线性规划：

    \(\max c^{T} x \quad \text { s.t. } \quad A x \leq b, x \geq 0\)

     把其中的等式约束变成不等式约束，可得

    \(\min c^{T} x \quad \text { s.t. }\left[\begin{array}{c} A \\ -A \end{array}\right] x \geq\left[\begin{array}{c} b \\ -b \end{array}\right], x \geq 0\)

    它的对偶问题是

    \(\max \left[\begin{array}{ll} b^{T} & -b^{T} \end{array}\right]\left[\begin{array}{l} y_{1} \\ y_{2} \end{array}\right] \quad \text { s.t. }\left[\begin{array}{ll} A^{T} & -A^{T} \end{array}\right]\left[\begin{array}{l} y_{1} \\ y_{2} \end{array}\right] \leq c\)

    其中 \(y_{1}\) 和\(y_{2}\)分别表示对应于约束 \(Ax\geq b\) 和 \(-Ax\geq -b\) 的对偶变量组。令\(y=y_{1}-y_{2}\)，则上式又可写成

    \(\min b^{T} y \quad \text { s.t. } \quad A^{T} y \geq c, y \geq 0\)

     原问题和对偶的对偶约束之间的关系：

    口诀：

    不等号方向相同：Max变量的不等号与min约束条件的不等号方向相同

    不等号方向相反：Max约束条件的不等号与min变量的不等号方向相反等号对应无约束:

    约束条件的“=”对应变量“无约束”
     

    例题1：

    试求下列线性规划的对偶问题。

    \(\begin{aligned} &\min Z=2 x_{1}+3 x_{2}-5 x_{3}+x_{4} \\ &\left\{\begin{array}{l} x_{1}+x_{2}-3 x_{3}+x_{4} \geq 5 \\ 2 x_{1}+2 x_{3}-x_{4} \leq 4 \\ x_{2}+x_{3}+x_{4}=6 \\ x_{1} \leq 0 ; x_{2}, x_{3} \geq 0 ; x_{4} \epsilon R\end{array}\right. \end{aligned}\)

    解:设对应于三个约束的对偶变量分别为y1、y2、y3，则原问题的对偶问题为:
    \(\begin{aligned} &\max Z^{\prime}=5 y_{1}+4 y_{2}+6 y_{3} \\ &\left\{\begin{array}{l} y_{1}+2 y_{2} \geq 2 \\ y_{1}+y_{3} \leq 3 \\ -3 y_{1}+2 y_{2}+y_{3} \leq-5 \\ y_{1}-y_{2}+y_{3}=1 \\ y_{1} \geq 0 ; y_{2} \leq 0 ; y_{3} \epsilon R \end{array}\right. \end{aligned}\)

    4.2 对偶问题的基本性质

    原问题和对偶问题：

    \(\begin{array}{ll} \max Z=C X & \min Z^{\prime}=Y b \\ \begin{cases}A X \leq b \\ X \geq 0\end{cases} & \left\{\begin{array}{l} Y A \geq C \\ Y \geq 0 \end{array}\right. \end{array}\)

    (1）对称性
    对偶问题的对偶是原问题
    (2）弱对偶性
    若\(\overline{X}\)是原问题(max)的可行解，\(\overline{Y}\)是对偶问题(min)的可行解，则\(C\overline{X}\leq \overline{Y}b\)，证：

    max问题：\(A\overline{X}\leq b\)    左乘\(\overline{Y}\)     \(\overline{Y}A\overline{X}\leq \overline{Y}b\)

    min问题：\(\overline{Y}A\geq C\)     右乘\(\overline{X}\)     \(\overline{Y}A\overline{X}\geq C\overline{X}\)

    (3)无界性
    若原问题（对偶问题）为无界解，则其对偶问题(原问题）无可行解
     

     (4）最优性
    设\(\widehat{X}\)是原问题的可行解，\(\widehat{Y}\)是对偶问题的可行解，当\(C\widehat{X}=\widehat{Y}b\)=府时，\(\widehat{X},\widehat{Y}\)是最优解

    (5）对偶定理
    若原问题有最优解，则对偶问题也有最优解，且日标函数值相等

    (6）互补松弛性
    设\(\widehat{X},\widehat{Y}\)分别是原问题和对偶问题的可行解，\(X_{s},Y_{s}\)分别是其松弛变量的可行解，则\(\widehat{X},\widehat{Y}\)是最优解当且仅当\(Y_{s}\widehat{X}=0\)和\(X_{s}\widehat{Y}=0\)

     例 10 已知线性规划问题:

    已知其对偶问题的最优解为\(y_{1}^{*}=\frac{4}{5}, y_{2}^{*}=\frac{3}{5}, z=5\)。试用对偶理论找出原问题的最优

    解 先写出它的对偶问题

    \(\begin{aligned} &\max z=4 y_{1}+3 y_{2} \\ &\left\{\begin{array}{l} y_{1}+2 y_{2} \leq 2 \\ y_{1}-y_{2} \leq 3 \\ 2 y_{1}+3 y_{2} \leq 5 \\ y_{1}+y_{2} \leq 2 \\ 3 y_{1}+y_{2} \leq 3 \\ y_{1}, y_{2} \geq 0 \end{array}\right. \end{aligned}\)     添加松弛变量    \(\left\{\begin{array}{l} y_{1}+2 y_{2}+y_{s_{1}}=2 \\ y_{1}-y_{2}+y_{s_{2}}=3 \\ 2 y_{1}+3 y_{2}+y_{s_{3}}=5 \\ y_{1}+y_{2}+y_{s_{4}}=2 \\ 3 y_{1}+y_{2}+y_{s_{5}}=3 \\ y_{1}, y_{2}, y_{s_{1}}, y_{s_{2}}, y_{s_{3}}, y_{s_{4}}, y_{s_{5}} \geq 0 \end{array}\right.\)    代入最优解 

     \(\left\{\begin{array}{l} y_{s_{1}}=0 \\ y_{s_{2}} \neq 0 \\ y_{s_{3}} \neq 0 \\ y_{s_{4}} \neq 0 \\ y_{s_{5}}=0 \end{array}\right.\)    互补松弛性       \(\left\{\begin{array}{l} y_{s_{2}} x_{2}^{*}=0 \\ y_{s_{3}} x_{3}^{*}=0 \\ y_{s_{4}} x_{4}^{*}=0 \end{array}\right.\)         求解      \(\left\{\begin{array}{l} x_{2}^{*}=0 \\ x_{3}^{*}=0 \\ x_{4}^{*}=0 \end{array}\right.\)

    再分析原问题:
    \(\left\{\begin{array}{l} x_{1}+x_{2}+2 x_{3}+x_{4}+3 x_{5}-x_{s_{1}}=4 \\ 2 x_{1}-x_{2}+3 x_{3}+x_{4}+x_{5}-x_{s_{2}}=3 \\ x_{1}, x_{2}, x_{3}, x_{4}, x_{5}, x_{s_{1}}, x_{s_{2}} \geq 0 \end{array}\right.\)     由于 \(y_{1}^{*} \neq 0, y_{2}^{*} \neq 0\)      \(\left\{\begin{array}{l} x_{s_{1}}=0 \\ x_{s_{2}}=0 \end{array}\right.\)     将 

    \(x_{s_{1}}, x_{s_{2}}, x_{2}^{*}, x_{3}^{*}, x_{4}^{*}\)代入得\(\left\{\begin{array}{l} x_{1}^{*}+3 x_{5}^{*}=4 \\ 2 x_{1}^{*}+x_{5}^{*}=3 \end{array}\right.\)   解   \(\left\{\begin{array}{l} x_{1}^{*}=1 \\ x_{5}^{*}=1 \end{array}\right.\)  最优解 ：  \(\begin{aligned} &X^{*}=(1,0,0,0,1)^{T} \\ &\omega^{*}=5 \end{aligned}\) 

    4.3 灵敏度分析

            提出这样两个问题：当这些系数有一个或几个发生变化时，已求得的线性规划问题的最优解会有什么变化；或者这些系数在什么范围内变化时，线性规划问题的最优解或最优基不变。这里我们就不讨论了。

    4.4 参数线性规划

    参数线性规划是研究\(a_{ij},b_{i},c_{j}\) 这些参数中某一参数连续变化时，使最优解发生变化 的各临界点的值。即把某一参数作为参变量，而目标函数在某区间内是这参变量的线性 函数，含这参变量的约束条件是线性等式或不等式。因此仍可用单纯形法和对偶单纯形法进行分析参数线性规划问题。

    4.4.1 单纯形法：

    例题：用单纯形法求解线性规划问题
    \(\begin{aligned} &\max z=2 x_{1}+x_{3} \\ &\text { s.t. }\left\{\begin{array}{c} 5 x_{2} \leq 15 \\ 6 x_{1}+2 x_{2} \leq 24 \\ x_{1}+x_{2} \leq 5 \\ x_{1}, x_{2} \geq 0 \end{array}\right. \end{aligned}\)

    （1）标准化

    \(\begin{aligned} &\max z=2 x_{1}+x_{2}+0 x_{3}+0 x_{4}+0 x_{5} \\ &\text { s.t. }\left\{\begin{array}{c} 5 x_{2}+x_{3}=15 \\ 6 x_{1}+2 x_{2}+x_{4}=24 \\ x_{1}+x_{2}+x_{5}=5 \\ x_{1}, x_{2}, x_{3}, x_{4}, x_{5} \geq 0 \end{array}\right. \end{aligned}\)

    （2）使用python编程：

    1. import numpy as np
    2. class Simplex:
    3.     def __init__(self) -> None:
    4.         self.solutionStatus = 0
    5.  
    6.     def set_param(self, A, b, c, baseInd):
    7.  
    8.         self.A = A
    9.         self.m, self.n = self.A.shape
    10.         self.b = b
    11.         self.c = c
    12.         self.baseInd = baseInd
    13.         self.nonbaseInd = [i for i in range(self.n) if i not in self.baseInd]
    14.         self.B = self.A[:, self.baseInd]
    15.         self.N = self.A[:, self.nonbaseInd]
    16.         self.B_inv = self.B.I
    17.  
    18.     def main(self):
    19.         # 计算非基变量的检验数，基变量检验数为0
    20.         reducedCost = [0] * self.n
    21.         X_B = None
    22.         simplexTable = SimplexTable()
    23.         while (self.solutionStatus == 0):
    24.             for i in range(self.n):
    25.                 reducedCost[i] = float(self.c[:, i] - np.dot(self.c[:, self.baseInd], self.A[:, i]))
    26.                 # p_j已经全部左乘B_inv，因此这里c_B可以直接与更新后的p_j相乘
    27.             self.solutionStatus = self.check_solutionStatus(reducedCost)
    28.             simplexTable.show(self.A, self.b, self.c, reducedCost, self.baseInd, self.solutionStatus)
    29.             if self.solutionStatus == 0:
    30.                 inVarInd = reducedCost.index(max(reducedCost))
    31.                 outVarInd = self.get_outVar(inVarInd)
    32.                 self.baseInd[self.baseInd.index(outVarInd)] = inVarInd
    33.                 self.nonbaseInd[self.nonbaseInd.index(inVarInd)] = outVarInd
    34.                 # 得到新基
    35.                 self.B = self.A[:, self.baseInd]
    36.                 self.B_inv = self.B.I
    37.                 self.b = np.dot(self.B_inv, self.b)
    38.                 X_B = self.b.T.tolist()[0]
    39.                 self.A[:, inVarInd] = self.A[:, outVarInd]
    40.                 self.A[:, self.nonbaseInd] = np.dot(self.B_inv, self.A[:, self.nonbaseInd])
    41.             else:
    42.                 break
    43.         if self.solutionStatus == 1:
    44.             solution = {}
    45.             for i in range(self.n):
    46.                 if i in self.baseInd:
    47.                     solution[i] = X_B[self.baseInd.index(i)]
    48.                 if i in self.nonbaseInd:
    49.                     solution[i] = 0
    50.             objctive = float(np.dot(np.dot(self.c[:, self.baseInd], self.B_inv), self.b))
    51.  
    52.             return solution, objctive, self.solutionStatus
    53.         else:
    54.             solution = None
    55.             objctive = None
    56.             return solution, objctive, self.solutionStatus
    57.  
    58.     def check_solutionStatus(self, reducedCost):
    59.         if all(x < 0 or x == 0 for x in reducedCost):
    60.             # 所有检验数<=0
    61.             if any(reducedCost[i] == 0 for i in self.nonbaseInd):
    62.                 # 存在非基变量的检验数为0
    63.                 return 2  # 无穷多最优解
    64.             else:
    65.                 return 1  # 唯一最优解
    66.         else:
    67.             if all(all(x < 0 for x in self.A[:, i]) for i in self.nonbaseInd if reducedCost[i] > 0):
    68.                 return 3  # 无界解
    69.             else:
    70.                 # 选择最大检验数对应的非基变量入基
    71.                 return 0
    72.  
    73.     def get_outVar(self, inVarInd):
    74.         inVarCoef = self.A[:, inVarInd]
    75.         ratio = [np.inf] * self.m
    76.         for i in range(len(inVarCoef)):
    77.             if float(inVarCoef[i, :]) > 0:
    78.                 ratio[i] = float(self.b[i, :]) / float(inVarCoef[i, :])
    79.         outVarInd = self.baseInd[ratio.index(min(ratio))]
    80.         return outVarInd
    81.  
    82. class SimplexTable:
    83.     def __init__(self):
    84.         # 左端表头格式设置
    85.         self.setting_left = '{:^8}'+'{:^8}'+'{:^8}|'
    86.         # 右端变量格式设置
    87.         self.setting_var = '{:^8}|'
    88.         # 求解状态格式设置
    89.         self.setting_status = '{:^24}|'+'{:^44}|'
    90.     def show(self,A,b,c,reducedCost,baseInd,solutionStatus):
    91.         var_num = A.shape[1]
    92.         # 打印系数行
    93.         c_j_list = c.flatten().tolist()[0]
    94.         print('='*80)
    95.         print(self.setting_left.format('','c_j',''),end = '')
    96.         for i in c_j_list:
    97.             print(self.setting_var.format(i),end = '')
    98.         print()
    99.         # 打印变量行
    100.         print(self.setting_left.format('c_B','x_B','b'),end = '')
    101.         for i in range(var_num):
    102.             print(self.setting_var.format('x'+str(i)),end = '')
    103.         print()
    104.         # 打印变量与矩阵
    105.         for i in range(len(baseInd)):
    106.             varInd = baseInd[i]
    107.             varCeof = round(float(c[:,varInd]),2)
    108.             varValue = round(float(b[i,:]),2)
    109.             row = A[i,:].flatten().tolist()[0]
    110.             print(self.setting_left.format(str(varCeof),'x'+str(varInd),str(varValue)),end='')
    111.             for i in range(var_num):
    112.                 print(self.setting_var.format(round(row[i],2)),end = '')
    113.             print()
    114.         # 打印检验数
    115.         print(self.setting_left.format('','c_j-z_j',''),end='')
    116.         for i in reducedCost:
    117.             print(self.setting_var.format(round(i,2)),end = '')
    118.         print()
    119.         # 显示求解状态
    120.         if solutionStatus == 0:
    121.             print(self.setting_status.format('status','...Iteration continues...'))
    122.         if solutionStatus == 1:
    123.             print(self.setting_status.format('status','Unique Optimal Solution'))
    124.             B = A[:,baseInd]
    125.             B_inv =  B.I
    126.             objctive = float(np.dot(np.dot(c[:, baseInd], B_inv), b))
    127.             print(self.setting_status.format('objctive', round(objctive,2)))
    128.             solution = [0]*var_num
    129.             X_B = b.T.tolist()[0]
    130.             for i in range(var_num):
    131.                 if i in baseInd:
    132.                     solution[i] = X_B[baseInd.index(i)]
    133.             print(self.setting_left.format('', 'solution', ''), end='')
    134.             for i in solution:
    135.                 print(self.setting_var.format(round(i, 2)), end='')
    136.             print()
    137.  
    138.         if solutionStatus == 2:
    139.             print(self.setting_status.format('status','Multiple Optimal Solutions'))
    140.             B = A[:, baseInd]
    141.             B_inv = B.I
    142.             objctive = float(np.dot(np.dot(c[:, baseInd], B_inv), b))
    143.             print(self.setting_status.format('objctive', round(objctive, 2)))
    144.         if solutionStatus == 3:
    145.             print(self.setting_status.format('status','Unboundedness'))
    146.  
    147. '''=============================================================================='''
    148. A = np.matrix([[0,5,1,0,0],
    149.                [6,2,0,1,0],
    150.                [1,1,0,0,1]], dtype=np.float64)
    151.                     # 必须设置精度，否则在替换矩阵列时会被自动圆整
    152. b = np.matrix([[15], [24], [5]], dtype=np.float64)
    153. c = np.matrix([2,1,0,0,0], dtype=np.float64)
    154. baseInd = [2,3,4]
    155.  
    156. s=Simplex()
    157. s.set_param( A, b, c, baseInd)
    158. s.main()

    以此形式进行列表求解,满足单纯形法的基本条件，具体如下：

    1. ================================================================================
    2.           c_j           |  2.0   |  1.0   |  0.0   |  0.0   |  0.0   |
    3.   c_B     x_B      b    |   x0   |   x1   |   x2   |   x3   |   x4   |
    4.   0.0      x2     15.0  |  0.0   |  5.0   |  1.0   |  0.0   |  0.0   |
    5.   0.0      x3     24.0  |  6.0   |  2.0   |  0.0   |  1.0   |  0.0   |
    6.   0.0      x4     5.0   |  1.0   |  1.0   |  0.0   |  0.0   |  1.0   |
    7.         c_j-z_j         |  2.0   |  1.0   |  0.0   |  0.0   |  0.0   |
    8.          status         |         ...Iteration continues...          |
    9. ================================================================================
    10.           c_j           |  2.0   |  1.0   |  0.0   |  0.0   |  0.0   |
    11.   c_B     x_B      b    |   x0   |   x1   |   x2   |   x3   |   x4   |
    12.   0.0      x2     15.0  |  0.0   |  5.0   |  1.0   |  0.0   |  0.0   |
    13.   2.0      x0     4.0   |  1.0   |  0.33  |  0.0   |  0.17  |  0.0   |
    14.   0.0      x4     1.0   |  0.0   |  0.67  |  0.0   | -0.17  |  1.0   |
    15.         c_j-z_j         |  0.0   |  0.33  |  0.0   | -0.33  |  0.0   |
    16.          status         |         ...Iteration continues...          |
    17. ================================================================================
    18.           c_j           |  2.0   |  1.0   |  0.0   |  0.0   |  0.0   |
    19.   c_B     x_B      b    |   x0   |   x1   |   x2   |   x3   |   x4   |
    20.   0.0      x2     7.5   |  0.0   |  0.0   |  1.0   |  1.25  |  -7.5  |
    21.   2.0      x0     3.5   |  1.0   |  0.0   |  0.0   |  0.25  |  -0.5  |
    22.   1.0      x1     1.5   |  0.0   |  1.0   |  0.0   | -0.25  |  1.5   |
    23.         c_j-z_j         |  0.0   |  0.0   |  0.0   | -0.25  |  -0.5  |
    24.          status         |          Unique Optimal Solution           |
    25.         objctive        |                    8.5                     |
    26.         solution        |  3.5   |  1.5   |  7.5   |   0    |   0    |

    最优解：\({X}^{*}=(3.5,1.5,7.5,0,0)^{T}\)

    最优值：\(z^{*}=8.5\)

    单纯形法的进一步讨论：

    1.人工变量法（大M法）

    一些线性规划问题在化为标准形后约束条件系数矩阵不存在单位矩阵，就需要在约束条件中添加人工变量构造一个新的线性规划问题。将人工变量在目标函数中的系数取任意大的负数，在求解过程中可以防止人工变量出现在最优解中。若在最终单纯形表中所有检验数都小于等于零，但基变量中仍存在不为零的人工变量，则问题无解。下面，给出人工变量法的计算步骤。

    2.两阶段法

    用大M法处理人工变量，在用电子计算机求解时，对M只能在计算机内输入一个机器最大字长的数字。如果线性规划问题中的参数值与这个代表M的数比较接近，或远远小于这个数字，由于计算机计算时取值上的误差，有可能使计算结果发生错误。

    为了克服这个困难，可以对添加人工变量后的线性规划问题分两个阶段来计算，称两阶段法。

    进一步了解：运筹说 第16期 | 线性规划硬核知识点梳理—单纯形法 - 知乎

    4.4.2 对偶单纯形法：

    例题：用对偶单纯形法求解LP：

    \(\begin{aligned} &\operatorname{Min} W=2 x_{1}+3 x_{2}+4 x_{3} \\ &\text { s.t. }\left\{\begin{array}{c} x_{1}+2 x_{2}+x_{3} \geq 3 \\ 2 x_{1}-x_{2}+3 x_{3} \geq 4 \\ x_{1}, x_{2}, x_{3} \geq 0 \end{array}\right. \end{aligned}\)

    （1）标准化：

    \(\begin{aligned} &\operatorname{Max} z=-2 x_{1}-3 x_{2}-4 x_{3} \\ &\text { s.t. }\left\{\begin{array}{c} x_{1}+2 x_{2}+x_{3}-x_{4}=3 \\ 2 x_{1}-x_{2}+3 x_{3}-x_{5}=4 \\ x_{1}, x_{2}, x_{3}, x_{4}, x_{5} \geq 0 \end{array}\right. \end{aligned}\)

    将两个等式约束两边分别乘以-1，得

    \(\begin{aligned} &\operatorname{Max}z=-2 x_{1}-3 x_{2}-4 x_{3} \\ &\text { s.t. }\left\{\begin{array}{c} -x_{1}-2 x_{2}-x_{3}+x_{4}=-3 \\ -2 x_{1}+x_{2}-3 x_{3}+x_{5}=-4 \\ x_{1}, x_{2}, x_{3}, x_{4}, x_{5} \geq 0 \end{array}\right. \end{aligned}\)

    （2）使用python编程：

    1. import re
    2. import pandas as pd
    3. from fractions import Fraction
    4.  
    5. # data = pd.read_csv("data2.csv", header=0, index_col=0)
    6. # data=pd.DataFrame({"x1":[-2,-2,-1,-9],"x2":[-2,-3,-1,-12],"x3":[-1,-1,-5,-15],"x4":[1,0,0,0],"x5":[0,1,0,0],"x6":[0,0,1,0],"b":[-10,-12,-14,0]},index = ['x4', 'x5', 'x6',"sigma"])
    7. data=pd.DataFrame({"x1":[-1,-2,-2],"x2":[-2,1,-3],"x3":[-1,-3,-4],"x4":[1,0,0],"x5":[0,1,0],"b":[-3,-4,0]},index = ['x4', 'x5' ,"sigma"])
    8. print("=======原始表=======")
    9.  
    10. print(data)
    11.  
    12.  
    13. def do_it(data):
    14.     in_base_index = data.loc["sigma", :][data.loc['sigma', :] < 0].index
    15.     rec = {}
    16.     res_val = []
    17.     for i in in_base_index:
    18.         out_base_index = data.loc[:, i][data.loc[:, i] < 0].index.drop("sigma")
    19.         for j in out_base_index:
    20.             res = Fraction(data.loc["sigma", i], data.loc[j, i])
    21.             res_val.append(res)
    22.             rec[str(j) + "," + str(i)] = res
    23.  
    24.     def get_key(dic, value):
    25.         return [k for k, v in dic.items() if v == value]
    26.  
    27.     s = get_key(rec, min(res_val))
    28.     s = re.split(r"\,", s[0])
    29.  
    30.     # 这里是将本身变成1
    31.     param1 = Fraction(1, data.loc[s[0], s[1]])
    32.     data.loc[s[0], :] = data.loc[s[0], :] * param1
    33.     # 将其他变为0
    34.     for row in data.index.drop(s[0]):
    35.         target = data.loc[row, s[1]]
    36.         param2 = Fraction(-target, 1)
    37.         data.loc[row, :] = data.loc[s[0], :] * param2 + data.loc[row, :]
    38.  
    39.     data = data.rename(index={s[0]: s[1]})
    40.     print("================================")
    41.     print(data)
    42.     if (data["b"].drop("sigma") < 0).any():
    43.         print("需要继续迭代！")
    44.         do_it(data)
    45.     else:
    46.         print("迭代结束！")
    47.     return data
    48.  
    49.  
    50. # 如何b中的任何一个数小于零，都需要进一步操作
    51. if (data["b"].drop("sigma") < 0).any():
    52.     print("需要继续迭代！")
    53.     do_it(data)
    54. else:
    55.     print("迭代结束！")

    以此形式进行列表求解,满足对偶单纯形法的基本条件，具体如下：

    1. =======原始表=======
    2.        x1  x2  x3  x4  x5  b
    3. x4     -1  -2  -1   1   0 -3
    4. x5     -2   1  -3   0   1 -4
    5. sigma  -2  -3  -4   0   0  0
    6. 需要继续迭代！
    7. ================================
    8.       x1    x2   x3 x4    x5   b
    9. x4     0  -5/2  1/2  1  -1/2  -1
    10. x1     1  -1/2  3/2  0  -1/2   2
    11. sigma  0    -4   -1  0    -1   4
    12. 需要继续迭代！
    13. ================================
    14.       x1 x2    x3    x4    x5     b
    15. x2     0  1  -1/5  -2/5   1/5   2/5
    16. x1     1  0   7/5  -1/5  -2/5  11/5
    17. sigma  0  0  -9/5  -8/5  -1/5  28/5
    18. 迭代结束！

    最优解：\({X}^{*}=(11 / 5,2 / 5,0,0,0) ^{T}\)

    最优值：\(\begin{aligned} &\operatorname{minW}=-\max { }^{*}= -[11 / 5 \times(-2)+2 / 5 \times(-3)]=28 / 5 \end{aligned}\)
     

  • 相关阅读:
    k8s常用命令
    代码随想录训练营第36天|LeetCode 435. 无重叠区间、763.划分字母区间、 56. 合并区间
    当三年前端开发掌握了工程化，真就无敌了？
    信息系统项目管理师---第七章项目成本管理历年考题
    牛客刷题之图论-最小生成树
    C语言17---计算机的存储规则
    小波神经网络的基本原理,小波神经网络设计方案
    JAVA毕业设计服装连锁店后台管理系统计算机源码+lw文档+系统+调试部署+数据库
    华为欧拉 openEuler 23.09 一键安装 Oracle 12CR2 单机
    初级永磁直线电机双动子电流镜像容错控制
  • 原文地址：https://blog.csdn.net/qq_21402983/article/details/126313619