单纯形法学习笔记
mzgwty
·
·
个人记录
首先，在高中的学习中，比较简单的会有如下问题：
\left\{
\begin{array}{l}
x+y\le 5 \\
x-y\le 5 \\
z=2x+y
\end{array}
\right.
然后叫你求z的最大值/最小值
为了便于理解，这里使用三维线性规划
\left\{
\begin{array}{l}
2x+y\le 10\\
x+z\le 5\\
x\le y\\
c=3x+2y+4z
\end{array}
\right.
线性规划有一个标准形式，即我们需要满足的全部都是等式而不是不等式，并且我们要求的是最大值，且每一个未知数都必须是非负数
如果要求最小值的话，直接先对z各项系数取负求出最大值后在再输出其相反数
然后为了方便，我们规定所有不等式之间线性无关
我们可以通过增加新的变量将一切不等式化成等式。
如果是x+y\le c，那么写成x+y+nx=c
如果是x+y\ge c，那么写成x+y-nx=c
比如上面的例子
\left\{
\begin{array}{l}
2x_1+x_2+x_4=10\\
x_1+x_3+x_5=5\\
x_1-x_2+x_6=0\\
c=3x_1+2x_2+4x_3\\
x_1,x_2,x_3,x_4,x_5,x_6\le0
\end{array}
\right.
我们还可以将它表示成这个样子
\max c=3x_1+2x_2+4x_3+0x_4+0x_5+0x_6
\begin{pmatrix}x_4\\x_5\\x_6\end{pmatrix}+
\begin{pmatrix}2\\1\\1\end{pmatrix}x_1+
\begin{pmatrix}1\\0\\-1\end{pmatrix}x_2+
\begin{pmatrix}0\\1\\0\end{pmatrix}x_3=
\begin{pmatrix}10\\5\\0\end{pmatrix}
在这个式子中我们可以用x_1,x_2,x_3表示出x_4,x_5,x_6，所以我们称x_4,x_5,x_6为基变量，x_1,x_2,x_3为非基变量，然后始终用非基变量来表示c，由于保证x_i非负，即下界为0，我们只需要让最后表示c的非基变量前面的系数都为负数就行了，最大值就是那个常数
我们发现在c中此时x_1前的系数是正的，所以我们考虑把x_1加到基里面去
首先要将x_1在那一个位置的系数变为1，其他位置上的系数为0，这一步通过行消元来解决
那么我们就必须让基里面的元素出来一个，那这个元素要怎么找呢？
其实只需要找对应等式右侧常数除以对应x_1系数最小的就行了，证明就是行消元之后其他等式的常数项不能为负数，所以要找这个值最小的。当然，前提是x_1系数非负，如果是负数的话式子右侧的那个常数就会变成负数。
那么我们就找到了x_6
对原式进行行变换得到
\begin{pmatrix}x_4-2x_6\\x_5-x_6\\x_6\end{pmatrix}+
\begin{pmatrix}0\\0\\1\end{pmatrix}x_1+
\begin{pmatrix}3\\1\\-1\end{pmatrix}x_2+
\begin{pmatrix}0\\1\\0\end{pmatrix}x_3=
\begin{pmatrix}10\\5\\0\end{pmatrix}
交换x_1和x_6的位置，得到
\begin{pmatrix}x_4\\x_5\\x_1\end{pmatrix}+
\begin{pmatrix}-2\\-1\\1\end{pmatrix}x_6+
\begin{pmatrix}3\\1\\-1\end{pmatrix}x_2+
\begin{pmatrix}0\\1\\0\end{pmatrix}x_3=
\begin{pmatrix}10\\5\\0\end{pmatrix}
c=3x_1+2x_2+4x_3
=3(x_2-x_6)+2x_2+4x_3
=5x_2+4x_3-3x_6
然后c中x_2的系数仍然为正，继续上述操作，先找到常数项除以x_2系数最小的那一个，即x_4
先将x_2该项前系数变为1，得到
\begin{pmatrix}\frac13x_4\\x_5\\x_1\end{pmatrix}+
\begin{pmatrix}-\frac23\\-1\\1\end{pmatrix}x_6+
\begin{pmatrix}1\\1\\-1\end{pmatrix}x_2+
\begin{pmatrix}0\\1\\0\end{pmatrix}x_3=
\begin{pmatrix}\frac{10}3\\5\\0\end{pmatrix}
进行行变换，得到
\begin{pmatrix}\frac13x_4\\x_5-\frac13x_4\\x_1+\frac13x_4\end{pmatrix}+
\begin{pmatrix}-\frac23\\-\frac13\\\frac13\end{pmatrix}x_6+
\begin{pmatrix}1\\0\\0\end{pmatrix}x_2+
\begin{pmatrix}0\\1\\0\end{pmatrix}x_3=
\begin{pmatrix}\frac{10}3\\\frac53\\\frac{10}3\end{pmatrix}
交换x_2和x_4的位置有
\begin{pmatrix}x_2\\x_5\\x_1\end{pmatrix}+
\begin{pmatrix}-\frac23\\-\frac13\\\frac13\end{pmatrix}x_6+
\begin{pmatrix}\frac13\\-\frac13\\\frac13\end{pmatrix}x_4+
\begin{pmatrix}0\\1\\0\end{pmatrix}x_3=
\begin{pmatrix}\frac{10}3\\\frac53\\\frac{10}3\end{pmatrix}
c=5x_2+4x_3-3x_6
=5(\frac{10-x_4+2x_6}3)+4x_3-3x_6
=\frac{50}3+4x_3-\frac53x_4+\frac13x_6
那么继续让x_3进基，将x_5出基，那么得到
\begin{pmatrix}x_2\\x_3\\x_1\end{pmatrix}+
\begin{pmatrix}-\frac23\\-\frac13\\\frac13\end{pmatrix}x_6+
\begin{pmatrix}\frac13\\-\frac13\\\frac13\end{pmatrix}x_4+
\begin{pmatrix}0\\1\\0\end{pmatrix}x_5=
\begin{pmatrix}\frac{10}3\\\frac53\\\frac{10}3\end{pmatrix}
c=\frac{50}3+4x_3-\frac53x_4+\frac13x_6
=\frac{50}3+4(\frac{5-3x_5+x_4+x_6}3)-\frac53x_4+\frac13x_6
=\frac{70}3-\frac13x_4-4x_5+\frac53x_6
继续让x_6进基，x_1出基，得到
\begin{pmatrix}x_2\\x_3\\x_6\end{pmatrix}+
\begin{pmatrix}2\\1\\3\end{pmatrix}x_1+
\begin{pmatrix}1\\0\\1\end{pmatrix}x_4+
\begin{pmatrix}0\\1\\0\end{pmatrix}x_5=
\begin{pmatrix}10\\5\\10\end{pmatrix}
则
c=\frac{70}3-\frac13x_4-4x_5+\frac53x_6
=\frac{70}3-\frac13x_4-4x_5+\frac{5(10-x_4-3x_1)}3
=40-5x_1-2x_4-4x_5
此时所有系数都小于0，得到最大值c=40
此时X=\begin{pmatrix}0,10,5,0,0,10\end{pmatrix}
以上是对于单纯形法过程的描述，下面给出代码
#include<stdio.h>
#include<string.h>
#include<math.h>
#include<algorithm>
#include<iostream>
#include<queue>
using namespace std;
const int N=3505;
const double INF=1e12,eps=1e-8;
int n,m;
double a[N][N];
inline void pivot(int l,int e) {
double t=a[l][e];a[l][e]=1;
for(int i=0 ; i<=m ; ++i) a[i][e]/=t;
for(int i=0 ; i<=m ; ++i)
if(i!=l&&fabs(a[i][e])>eps) {
t=a[i][e],a[i][e]=0;
for(int j=0 ; j<=n ; ++j) a[i][j]-=t*a[l][j];
}
}
inline double simplex() {
int l,e;
double t;
while(1) {
e=n+1;
for(int i=1 ; i<=n ; ++i) if(a[0][i]>eps) {e=i;break ;}
if(e==n+1) break ;
t=INF;
for(int i=1 ; i<=m ; ++i)
if(a[i][e]>eps&&t>a[i][0]/a[i][e])
l=i,t=a[i][0]/a[i][e];
if(t==INF) return INF;
pivot(l,e);
}
return a[0][0];
}
int main() {
scanf("%d%d",&n,&m);
for(int i=1 ; i<=m ; ++i) {
for(int j=1 ; j<=n ; ++j) scanf("%lf",&a[i][j]);
scanf("%lf",&a[i][0]);
}
for(int i=1 ; i<=n ; ++i) scanf("%lf",&a[0][i]);
printf("%.6lf\n",-simplex());
return 0;
}