线性规划单纯形法

一.解释

本文主要介绍单纯形法的基本原理，及程序的实现原理。本人在网上查阅资料时发现网上关于单纯形法的文章对不了解运筹学的人不太友好，因此我想在这写一个小白也能看懂的文章（默认学过线代）。本方法与传统的单纯形表法有所不同，进行了一些改进以便于程序运算。

二.线性规划的一般形式

一些朋友应该知道，求线性规划问题一般都是使用matlab或者是lingo，这些语言中有着对应的函数，python在导入scipy库之后也能找到函数linprog，但是这些函数都需要实现规定标准型以便于调用，这里拿matlab的函数linprog的一般形式来举例。

 其中：f,x,b,beq,lb,ub为列向量；A，Aeq为矩阵。一般，将f称为价值向量；x称为决策向量；b称为资源向量。

我们可以将类似的形式称作一般形式，一般形式中的约束方程可以由不等式或者等式确立。而决策向量x的上下限lb和ub可以使用一些手段来变成x‘≥0，其中x‘是新的决策向量，新旧之间存在着一定的转化关系以互相转化。

我们在此定义我们所用的一般形式：

 正常的线性规划问题，不是一般形式的都可以想些办法变成一般形式，求最小值可以通过取负变为求最大值，决策变量可以通过平移变为≥0。

三.线性规划标准形式

我们规定一般形式是为了方便代入，更容易编成函数后调用，而规定标准形式则是为了计算的便捷。

我们首先要明确的是，不等式可以通过引入新的受约束变量（一般令其≥0）来将不等式方程变为等式方程。例如：

 我们可以设一个x3，使x3满足：x3=3-x1+x2，同时约束x3≥0。故可以得到：

 因此，我们可以将上述一般形式，通过类似变换变为标准形式：

 显然，我们的变量相比一般形式要更多了，增加的个数等于约束方程中不等式的个数，这使得决策变量的个数要多于约束方程的个数，故m×n的矩阵A，应满足m＜n。所以由线性代数的知识可知，Ax=b这个式子的解应该有无穷多个（我们认为A是行满秩矩阵，即约束方程中的系数组合间线性无关）。

四. 基本原理

1.单纯形

我们简要介绍一下单纯形的概念。单纯形被称为最简单的拓扑模型，其与二维空间的三角形类似，在欧氏空间中，单纯形是一个凸集。我们都了解过凸函数的概念，即二阶导数大于等于0。而凸集也类似，首先我们先定义顶点，所谓的顶点，就是一个集合的边缘。而凸集，就是集合上任意两个顶点间所连线段上的任何一点也属于这个集合，类似的图形如圆、正方形。

2.基与基解

由上节可知m＜n，故对m×n阶矩阵A而言，必定能找到一些m×m阶方阵，使得这些方阵满秩，我们称这样的方阵为这个方程的一个基，记为B。B由一组属于A的列向量构成，因为A中的每个列向量都能对应一个决策变量x，故我们可以将基B所对应的变量组xB提出，组合成为新的方程：

 因为B是一个满秩矩阵，因此我们可以获得一组唯一解：

 这些确定的变量我们便称为基变量，而其他的变量构成的组合设为xN，这些变量我们称其为非基变量，易知这些非基变量都是不确定的变量，他们可以取任何值，如果我们让这些非基变量全部取零，那便构成了方程组的一个基解，我们知道在m×n阶矩阵A中可以找到m和n的组合数个m×m阶方阵，因此一般我们的基解都不止一个，而满足条件x≥0（非负约束）的基解，我们称其为基可行解，基解对应处在整个凸集的顶点，而我们所要找到使目标函数值最大的可行解，一定在凸集的某一顶点上，因此我们将使目标函数达到最大的基可行解称作

最优解，如若不存在基可行解，则称该方程无解，如果最优解对应无穷多个解，那称该方程有无穷多个解，在下面我们会给出两种情况的判断。

3.线性规划矩阵扩增形式

这个形式是我独创的形式，是为了方便使用Gauss消元来进行程序计算。线性代数的方程中提到了增广矩阵的形式，方程Ax=b的增广矩阵为[A,b],现如今我们将前面所述的标准形式在增广的基础上将目标函数的系数向量也增广在其中，即：

 易知扩增后的矩阵P为m+1×n+1阶矩阵，即相对A，在右侧及下方分别添加一排数据，在程序中不妨令所求值z为0，在经过初等变换的过程便可以得到最后的最优值。为表示方便，现在，我们用C来表示fT，C为行向量。因为我们需要引入新的变量来将不等式变为等式，因此目标函数的fT也需要扩充，我们如今用CB来表示基变量对应目标函数的系数组合，CN表示非基变量对应系数组合。下面将对程序原理，及最优检验进行说明。

我们可以看到，初等行变换并不会改变这个矩阵的目标值，因此经过初等变换所得到的矩阵与原矩阵在求解问题上是等价的，而处理过程中最好不使用初等列变换，因为这样会改变决策变量的顺序，然而在程序执行中，却可以引入一维整型数组，该数组用于储存列向量（决策变量）的顺序，然后用该数组调用矩阵，便可以起到保留变量顺序的作用，并且还不改变矩阵中数据的存储顺序。因而我们可以在程序中使用初等列变换中的列互换来进行求解，我们便可以放心的将基变量对应列向量放在前m列，这样可以化成标准的单位矩阵，更方便后面的运算。

4.基可行解的选取

想进行特定程序求解，就需要先选取一个初始的基可行解。对于（1）的矩阵P，我们首先对其进行排序处理，将目标向量系数C从小至大进行排列，这样处理的目的是更容易接近最优解，从而减少后面的转轴次数。因为基是m阶的方阵，我们首先选取排序处理后的矩阵前m-1个列向量作为基的一部分，随后将这一部分通过Gauss消元法化为单位矩阵，并将该单位矩阵的下方元素均通过Gauss消元法变为0。即变为：

 如上形式，单位矩阵为m-1阶。为了寻找到所需的基可行解，首先要将αm变为正（一般常系数均为正，如果是0则在m-1项中寻找一个正的进行替换），然后我们需要利用初等变换将前m-1个常系数α变为0，即用第m行依次向上减。

 最后应该变成(4)的形式，这种形式下我们便可以寻找基可行解。

定理1：利用初等变换化基为单位向量后，如若常系数均大于零，则该基解为一基可行解。

这个定理的证明很简单，由(1)得，B=E，因此解xB=b，故可知如若b有负项，则xB存在小于0的项，不满足非负约束，因此不是基可行解。

由这一定理我们可知，如果将(4)中m行乘以某些倍数分别加至前m-1行之后，可以产生一个单位矩阵的同时常系数要保持非负，这时便找到了所求基可行解，而若满足这一条件，应该对(4)有：

任意的i∈N，m≤i≤n-1，存在βmi＞0，使得任意j∈N，1≤j≤m-1，都有βji≤0   (5)

如若找不到满足这一条件的i，则在由β1m至β(m-1)n组成的矩阵中任意选择一个正数，交换前面单位矩阵中对应的列向量，重新将前面化为单位矩阵，然后不断重复这个过程，直至找到满足条件(5)的一组基。

如果αm＞0，且max(βmm，βm(m+1),……,βmn)≤0，则该线性规划方程无解。

5.检验数

检验数是程序的终止条件。上面所述找到一组可行解之后，要将其组合成单位矩阵，同时下方对应系数都变为0。这样做的目的是为了查看非基变量目标函数的系数是否均小于等于0，这是为了查看检验数在σj=zj-cj是否均小于等于0。为了方便叙述，我们先设σ为一行向量，以向量的形式来表述。我们由前面所述可知，目标函数z=C*x，C为一行向量，如果我们将其分成基变量组合xB以及非基变量组合xN，则有z=CBxB+CNxN。因为xB是可以求出来的，因此我们将CBxB用z0表示，而可以将CB变为0，与CN扩充，变成和C同阶的行向量，扩充之后的行向量用σ表示，xN也同样扩充变成一基可行解。因此我们的目标函数可表示为：

 因为基可行解x中的非基变量xN是可以取任意非负值的，如若σ<0成立，则非基变量均取0时变是可行域内的一个顶点，即基可行解，这一解也是我们所求的最优解。

在单纯形表中的检验数即为如上的原理，由Ax=b，以及BxB=b可知：

 再通过z，可以将检验数σ表示出来：

 

 以上式子提出x（只研究x的系数）后便是传统单纯形表法的检验数规则，这种方法相对比较麻烦。而我所设计的程序中使用的检验原理相同，虽然较复杂，但是更适合程序运行，可以在一定程度上降低时间复杂度。

定理2：对(2)扩增矩阵进行倍加行变换，或者对增广矩阵[A,b]进行初等行变换，不改变检测数及最优解。

这个定理的证明很简单，有一些线代基础的就可以进行验证。大致思路是先验证没有最下面一行的增广矩阵，这个比较简单，然后是算上最后一行的扩增矩阵，设某行的一个行向量的某一倍数，加到最下面一行，然后进行验证即可。

根据如上原理，我们只需用消元法对扩增矩阵进行处理，将目标函数中基变量的系数全部变为0，保留非基变量的系数，直接对非基变量的系数进行检验即可，如上(6)，，使用这种方法便可以略去传统过程中的求逆和矩阵乘积运算，直接使用Gauss消元法来替代。当目标函数中非基变量的系数均小于等于0时，我们便找到了最优解，最优解即为此时的z0。因为x有非负约束，故此时σx所能取到的最大值也只有0。若非基变量中存在σj=0，则该方程有无穷多个解，这些解均对应一个最优值。

6.转轴

一般而言，选取的第一个基可行解能直接满足检验条件成为最优解的概率极低，因此，我们需要对基变量中的一个变量进行替换，换走非基变量中的一个变量，组成一个新的基可行解，这一过程称为转轴，每次转轴之后都进行一次检验数检验，直到转到满足检验数要求的情况跳出循环。但是转轴又不能随意的转，程序又不知道，搞不好就会转入一个死循环，这时就需要一个规则来约束。

首先是新的基变量的选取，我们选择在目标函数当前系数中最大的一个，这样可以使替换之后的基变量尽可能大，按照前面所述规则，这个系数应该是正数。这是因为目标函数的值随着该系数所对变量的增大而增加的最快，因此选择他作为新的基变量。

然后是被替换的基变量的选取，这个时候就要看比值了，为了保证替换后的最右端系数列的所有值依旧为正，即保证始终为基可行解，如图(9)所示。

 假设Cs为所有目标函数系数中的最大值，根据上面新基变量选取规则，应该选取该系数所对应非基变量作为下一轮的新基变量。为了使用Gauss消元来将该列变为单位向量，需要选取一个主元，最后这一列只有该主元为1，其他均会变为0。因此，我们应保证转变后的常系数α始终大于0，因此我们应选择所有的

 中最小的一个，这样在消元的过程中就只会产生基可行解，从而保证了运算的稳定进行。

五.C语言算法数据结构

1.梳子

梳子这个东西在我家里话叫long梳，这种结构是一种树，我们调用矩阵的开端是梳子的柄，然后梳子的每一个刺都是一行数据。这种结构是依靠指针实现的二维数据，为什么用指针而不用c自带的数组呢，因为数组大小的定义是固定的，我们只能定义一个固定大小的数组，如果用来计算小的线性规划问题倒是可以，但是如果计算的是超大型的呢，程序要如何执行？c++为了克服这个问题而引入了vector，这是一种很好的东西，用起来很灵活，而我使用c直接编写，就会遇到很大的困难。但是没有c语言办不到的事，c语言的精华是什么呢，当然是指针了，利用指针，再使用stdlib.h库中的内存分配函数calloc，便可以克服上述问题，实现不固定大小的动态数据使用。

我们首先要定义一个二级指针，然后给这个二级指针分配m个连续的浮点数指针大小的内存，接下来利用for循环循环m次，每次分配n个连续的浮点数大小的内存，具体操作如下：

 而调用这个矩阵的i行j列的元素的指令为：*(*(A+i)+j)。

2.映射

我引用了函数的映射关系，这里是集合的映射，是为了对列互换的方便使用，这种方法没有改变数据的储存结构，只是在数据的调用上动手脚。

我们首先定义了一个矩阵，我们的操作相当于对这个类进行的，对矩阵的调用有特殊的格式，如上一节。我们的映射就是再定义一个整型数组（我喜欢指针，用的指针分配内存）row，这个row是用来存储列的序号的，没进行过列互换的矩阵记为A，在运算过程中进行了列互换的矩阵记作A‘，显然，存储的矩阵始终都是A，从未发生过变化，那么这个数组row的作用就是建立起A与A‘之间的桥梁，这样可以在保留原来列向量顺序的前提下进行一系列的运算。

输入现在矩阵的列序号j，通过row可以得到一个序号*(row+j)，这个序号对应着实际储存的矩阵的列序号，这个关系很像集合间的映射，因此我称其为映射。因此我们在运算过程中使用的矩阵i行j列元素调用命令为：*(*(A+i)+(*(row+j)))。

这个是c特有的复杂结构，其他语言都没有，这也是c独特魅力的一部分。