本文提供线性规划的单纯形法的简易代码，并给出一个退化现象的死循环算例，做简要分析。

单纯形法的简易代码

有关单纯形法的原理和迭代过程不做赘述，提供简易代码如下：

from typing import Tuple, Listimport numpy as np
import pandas as pdclass Simplex(object):"""Simplex method, iteration mode"""def __init__(self, a: List[List[float]], b: List[float], c: List[float], if_min: bool = True):"""Simplex method, initialise:param a:  A matrix:param b:  constraint params:param c:  cost params:param if_min:  if objective direction is minimisation,  default True"""self.a = aself.b = bself.c = cself.if_min = if_minself.num_var = len(c)  # number of variablesself.num_cons = len(b)  # number of constraintsdef run(self):"""Simplex method, main function:return:  Nothing"""# initialise simplex tablec = [-i for i in self.c] + [0]a = [self.a[i] + [self.b[i]] for i in range(len(self.a))]basic = list(range(self.num_var))[-self.num_cons:]num_iter = 0  # current number of iterations# print the initial simplex tabledic = {"basic": ['z'] + ["x_" + str(i) for i in basic]}for j in range(self.num_var):dic["x_" + str(j)] = [c[j]] + [a[i][j] for i in range(self.num_cons)]dic["solution"] = [c[-1]] + [a[i][-1] for i in range(self.num_cons)]df_iter = pd.DataFrame({k: dic[k]for k in ["basic"] + ["x_" + str(j) for j in range(self.num_var)] + ["solution"]})print("iteration {}:".format(num_iter))print(df_iter, '\n')num_iter += 1# simplex iterationsis_correct = Truewhile is_correct and min(c[: -1]) < 0:print("iteration {}:".format(num_iter))is_correct, c, a, basic = self.iteration(c=c, a=a, basic=basic)num_iter += 1# solutionsdict_solution = {}for i in range(len(basic)):dict_solution["x_" + str(basic[i])] = a[i][-1]for i in sorted(list(dict_solution.keys())):print("{0}:  {1}".format(i, dict_solution[i]))print()opt_objective = c[-1]print("optimal objective:  {}".format(opt_objective))def iteration(self, c: List[float], a: List[List[float]],basic: List[int]) -> Tuple[bool, List[float], List[List[float]], List[int]]:"""Simplex method, a single iteration:param c:  extended cost row, before current iteration:param a:  extended A matrix, before current iteration:param basic:  current basic variable indices, before current iteration:return: c:  extended cost row, after current iteration:return: a:  extended A matrix, after current iteration:return: basic:  current basic variables, after current iteration"""# entering variablevar_in = c[: -1].index(min(c[: -1]))pivot_row_index, ratio = -1, np.Inf  # todo: initialise pivot rowfor i in range(len(a)):if a[i][var_in] > 0:if a[i][-1] / a[i][var_in] < ratio:pivot_row_index, ratio = i, a[i][-1] / a[i][var_in]# error case: unboundedif pivot_row_index == -1:print("no leaving variable, model unbounded!")return False, c, a, basic# update A matrix and cost rowa[pivot_row_index] = [a[pivot_row_index][j] * (1 / a[pivot_row_index][var_in]) if j != var_in else 1for j in range(len(a[pivot_row_index]))]for i in range(len(a)):if i != pivot_row_index:if a[i][var_in]:a[i] = [a[i][j] - a[pivot_row_index][j] * (a[i][var_in] / a[pivot_row_index][var_in])for j in range(len(a[i]))]c = [c[j] - a[pivot_row_index][j] * (c[var_in] / a[pivot_row_index][var_in]) for j in range(len(c))]# update basic variable indicesbasic = [basic[i] if i != pivot_row_index else var_in for i in range(len(basic))]# print the simplex tabledic = {"basic": ['z'] + ["x_" + str(i) for i in basic]}for j in range(self.num_var):dic["x_" + str(j)] = [c[j]] + [a[i][j] for i in range(self.num_cons)]dic["solution"] = [c[-1]] + [a[i][-1] for i in range(self.num_cons)]df_iter = pd.DataFrame({k: dic[k]for k in ["basic"] + ["x_" + str(j) for j in range(self.num_var)] + ["solution"]})print(df_iter, '\n')return True, c, a, basic

注：该代码仅可计算 ≤\leq≤ 约束的线性规划模型，即不含处理人工变量的两阶段法，且输入须为标准格式。

退化现象的死循环算例

下面给出一个退化算例，算法在迭代过程中会出现死循环：

数学模型

数学模型如下：


max⁡\maxmax	34x1−20x2+12x3−6x4\frac{3}{4} x_1 - 20 x_2 + \frac{1}{2} x_3 - 6 x_443x1−20x2+21x3−6x4
s.t.s.t.s.t.	14x1−8x2−x3+9x4≤0\frac{1}{4} x_1 - 8 x_2 - x_3 + 9 x_ 4 \leq 041x1−8x2−x3+9x4≤0
	12x1−12x2−12x3+3x4≤0\frac{1}{2} x_1 - 12 x_2 - \frac{1}{2} x_3 + 3 x_ 4 \leq 021x1−12x2−21x3+3x4≤0
	x3≤1x_3 \leq 1x3≤1
	x1,x2,x3,x4≥0x_1, x_2, x_3, x_4 \geq 0x1,x2,x3,x4≥0

算法调用代码如下：

import Simplexa = [[1/4, -8, -1, 9, 1, 0, 0],[1/2, -12, -1/2, 3, 0, 1, 0],[0, 0, 1, 0, 0, 0, 1]]
b = [0, 0, 1]
c = [3/4, -20, 1/2, -6, 0, 0, 0]
if_min = Falsemodel = Simplex(a=a, b=b, c=c, if_min=if_min)
model.run()

迭代过程的单纯形表

iteration 0:

basic	x1x_1x1	x2x_2x2	x3x_3x3	x4x_4x4	x5x_5x5	x6x_6x6	x7x_7x7	solution
zzz	−34-\frac{3}{4}−43	202020	−12-\frac{1}{2}−21	666	000	000	000	000
x5x_5x5	14\frac{1}{4}41	−8-8−8	−1-1−1	999	111	000	000	000
x6x_6x6	12\frac{1}{2}21	−12-12−12	−12-\frac{1}{2}−21	333	000	111	000	000
x7x_7x7	000	000	111	000	000	000	111	111

iteration 1:

basic	x1x_1x1	x2x_2x2	x3x_3x3	x4x_4x4	x5x_5x5	x6x_6x6	x7x_7x7	solution
zzz	000	−4-4−4	−72-\frac{7}{2}−27	333333	333	000	000	000
x1x_1x1	111	−32-32−32	−4-4−4	363636	444	000	000	000
x6x_6x6	000	444	32\frac{3}{2}23	−15-15−15	−2-2−2	111	000	000
x7x_7x7	000	000	111	000	000	000	111	111

iteration 2:

basic	x1x_1x1	x2x_2x2	x3x_3x3	x4x_4x4	x5x_5x5	x6x_6x6	x7x_7x7	solution
zzz	000	000	−2-2−2	181818	111	111	000	000
x1x_1x1	111	000	888	−84-84−84	−12-12−12	888	000	000
x2x_2x2	000	111	38\frac{3}{8}83	−154-\frac{15}{4}−415	−12-\frac{1}{2}−21	14\frac{1}{4}41	000	000
x7x_7x7	000	000	111	000	000	000	111	111

iteration 3:

basic	x1x_1x1	x2x_2x2	x3x_3x3	x4x_4x4	x5x_5x5	x6x_6x6	x7x_7x7	solution
zzz	14\frac{1}{4}41	000	000	−3-3−3	−2-2−2	333	000	000
x3x_3x3	18\frac{1}{8}81	000	111	−212-\frac{21}{2}−221	−32-\frac{3}{2}−23	111	000	000
x2x_2x2	−364-\frac{3}{64}−643	111	000	316\frac{3}{16}163	116\frac{1}{16}161	−18-\frac{1}{8}−81	000	000
x7x_7x7	−18-\frac{1}{8}−81	000	000	212\frac{21}{2}221	32\frac{3}{2}23	−1-1−1	111	111

iteration 4:

basic	x1x_1x1	x2x_2x2	x3x_3x3	x4x_4x4	x5x_5x5	x6x_6x6	x7x_7x7	solution
zzz	−12-\frac{1}{2}−21	161616	000	000	−1-1−1	111	000	000
x3x_3x3	−52-\frac{5}{2}−25	565656	111	000	222	−6-6−6	000	000
x4x_4x4	−14-\frac{1}{4}−41	163\frac{16}{3}316	000	111	13\frac{1}{3}31	−23-\frac{2}{3}−32	000	000
x7x_7x7	52\frac{5}{2}25	−56-56−56	000	000	−2-2−2	666	111	111

iteration 5:

basic	x1x_1x1	x2x_2x2	x3x_3x3	x4x_4x4	x5x_5x5	x6x_6x6	x7x_7x7	solution
zzz	−74-\frac{7}{4}−47	444444	12\frac{1}{2}21	000	000	−2-2−2	000	000
x5x_5x5	−54-\frac{5}{4}−45	282828	12\frac{1}{2}21	000	111	−3-3−3	000	000
x4x_4x4	16\frac{1}{6}61	−4-4−4	−16-\frac{1}{6}−61	111	000	13\frac{1}{3}31	000	000
x7x_7x7	000	000	111	000	000	000	111	111

iteration 6:

basic	x1x_1x1	x2x_2x2	x3x_3x3	x4x_4x4	x5x_5x5	x6x_6x6	x7x_7x7	solution
zzz	−34-\frac{3}{4}−43	202020	−12-\frac{1}{2}−21	666	000	000	000	000
x5x_5x5	14\frac{1}{4}41	−8-8−8	−1-1−1	999	111	000	000	000
x6x_6x6	12\frac{1}{2}21	−12-12−12	−12-\frac{1}{2}−21	333	000	111	000	000
x7x_7x7	000	000	111	000	000	000	111	111

第 6 次迭代后的单纯形表与初始表相同，循环长度为 6

产生死循环的原因分析

每次迭代至少有两个基变量取 0，且其中至少有一个为候选的出基变量。

解决办法

摄动法：给约束右侧的 bbb 列叠加一个足够小的数，使其不为 0，可以从根本上消除退化现象，且不影响求解结果。

例如，可将上述算例中 bbb 列的第一个参数 0 修改为 10−810^{-8}10−8，即：


max⁡\maxmax	34x1−20x2+12x3−6x4\frac{3}{4} x_1 - 20 x_2 + \frac{1}{2} x_3 - 6 x_443x1−20x2+21x3−6x4
s.t.s.t.s.t.	14x1−8x2−x3+9x4≤10−8\frac{1}{4} x_1 - 8 x_2 - x_3 + 9 x_ 4 \leq 10^{-8}41x1−8x2−x3+9x4≤10−8
	12x1−12x2−12x3+3x4≤0\frac{1}{2} x_1 - 12 x_2 - \frac{1}{2} x_3 + 3 x_ 4 \leq 021x1−12x2−21x3+3x4≤0
	x3≤1x_3 \leq 1x3≤1
	x1,x2,x3,x4≥0x_1, x_2, x_3, x_4 \geq 0x1,x2,x3,x4≥0

则循环消失，迭代过程的单纯形表如下：

iteration 0:

basic	x1x_1x1	x2x_2x2	x3x_3x3	x4x_4x4	x5x_5x5	x6x_6x6	x7x_7x7	solution
zzz	−34-\frac{3}{4}−43	202020	−12-\frac{1}{2}−21	666	000	000	000	000
x5x_5x5	14\frac{1}{4}41	−8-8−8	−1-1−1	999	111	000	000	10−810^{-8}10−8
x6x_6x6	12\frac{1}{2}21	−12-12−12	−12-\frac{1}{2}−21	333	000	111	000	000
x7x_7x7	000	000	111	000	000	000	111	111

iteration 1:

basic	x1x_1x1	x2x_2x2	x3x_3x3	x4x_4x4	x5x_5x5	x6x_6x6	x7x_7x7	solution
zzz	000	222	−54-\frac{5}{4}−45	212\frac{21}{2}221	000	32\frac{3}{2}23	000	000
x5x_5x5	000	−2-2−2	−34-\frac{3}{4}−43	152\frac{15}{2}215	111	−12-\frac{1}{2}−21	000	10−810^{-8}10−8
x1x_1x1	111	−24-24−24	−1-1−1	666	000	222	000	000
x7x_7x7	000	000	111	000	000	000	111	111

iteration 2:

basic	x1x_1x1	x2x_2x2	x3x_3x3	x4x_4x4	x5x_5x5	x6x_6x6	x7x_7x7	solution
zzz	000	222	000	212\frac{21}{2}221	000	32\frac{3}{2}23	54\frac{5}{4}45	54\frac{5}{4}45
x5x_5x5	000	−2-2−2	000	152\frac{15}{2}215	111	−12-\frac{1}{2}−21	34\frac{3}{4}43	34\frac{3}{4}43
x1x_1x1	111	−24-24−24	000	666	000	222	111	000
x3x_3x3	000	000	111	000	000	000	111	111

最优解：
zmax⁡=54x1=1x3=1(x5=34+10−8)z_{\max} = \frac{5}{4} \quad x_1 = 1 \quad x_3 = 1 \quad (x_5 = \frac{3}{4} + 10^{-8}) zmax=45x1=1x3=1(x5=43+10−8)

单纯形法的代码实现与退化算例相关推荐

配电网重构|基于新颖的启发式算法SOE的随机(SDNR)配电网重构（Matlab代码实现）【算例33节点、84节点、119节点、136节点、417节点】
详解 Benders 分解与一个算例的 python 代码
听说过 benders 分解几年了,在生产管理.路径规划与选址问题中经常应用,一直没有细看,最近论文里面也见到,还是有必要了解一下它的基本思想与用法的. 目录 1. 基本思想 2. 线性规划模型与对偶 ...
中文文本纠错算例实现（有算例完整代码）
概述文本纠错又称为拼写错误或者拼写检查,由于纯文本往往来源于手打或者OCR识别,很可能存在一些错误,因此此技术也是一大关键的文本预处理过程,一般存在两大纠错类型. 1拼写错误第一种是Non-wor ...
动态规划原理介绍（附7个算例，有代码讲解）
动态规划思想动态规划(Dynamicprogramming)是一种通过把原问题分解为相对简单的子问题的方式求解复杂问题的方法. 动态规划常常适用于有重叠子问题和最优子结构性质的问题,动态规划方法所耗 ...
微电网日前优化调度。算例有代码（1）
个人电气博文传送门学好电气全靠它,个人电气博文目录(持续更新中-) 符号说明问题1 求解经济性评估方案: 若微网中蓄电池不作用,且微网与电网交换功率无约束,在无可再生能源和可再生能源全额利用两 ...
HMM算例 python 有代码
原理原理文字来源于 https://www.cnblogs.com/lcj1105/p/4936103.html 隐马尔可夫(HMM). 还是用最经典的例子,掷骰子.假设我手里有三个不同的骰子.第一 ...
矩形障碍算例（附Fortran计算代码及MATLAB后处理代码）
! 我又来啦,做了个小算例,在100*100个格子的空腔中加了几个矩形障碍 !上代码!西边界为水流进入区 parameter (n=100,m=100) real f(0:8,0:n,0:m) rea ...
万字字符长文带你了解遗传算法（有几个算例源码）
一.遗传算法的基本概念简单而言,遗传算法使用群体搜索技术,将种群代表一组问题解, 通过对当前种群施加选择.交叉和变异等一系列遗传操作来产生新-一代的种群,并逐步使种群进化到包含近似最优解的状态.由于 ...
qr分解求线性方程组_梯度下降求解线性方程组算例设计
凸二次优化问题 Theory. 设是实对称正定矩阵, ,则求解凸二次优化问题等价于求解线性方程组 . Proof. 二次型二阶可导,极小值点处梯度为零,现对优化的目标函数求梯度.二次型本质上具有: ...
DTW（动态时间归整）算法与DTW，fastDTW的python算例
文章目录 DTW算法思路 DTW算例原生代码计算 dtaidistance库计算DTW fastDTW算例参考资料 DTW算法思路论文地址:https://irit.fr/~Julien.Pin ...

单纯形法的代码实现与退化算例