單純形法與線性規劃入門
前言
說實話其實這個演算法真的不難，只要你有耐心把它啃完。。。
線性規劃
線性規劃的標準型長這個樣子：
\[\begin{align}max\ \ &\sum_{j=1}^nc_jx_j&\\s.t.\ \ &\sum_{j=1}^na_{i,j}x_j\le b_i&i=1,2,...,m\\&x_j\ge0&j=1,2,...,n\end{align} \]
它的鬆弛型則是長這個樣子：
\[\begin{align}max\ \ &\sum_{j=1}^nc_jx_j&\\s.t.\ \ &x_{n+i}=b_i-\sum_{j=1}^na_{i,j}x_j&i=1,2,...,m\\&x_j\ge0&j=1,2,...,n+m\end{align} \]
考慮和原本的標準型相比，鬆弛型有什麼優勢。
容易發現，原本限制條件中的不等號，被我們轉化成了等號。
這一步轉化是之後所有操作的基礎。
\(Pivot\)（轉軸）
考慮到我們需要最大化的式子\(\sum_{j=1}^nc_jx_j\)。
如果我們能夠通過某些手段把\(x\)的係數全部都變成負數，那麼因為\(x_j\ge0\)，只要令\(x_{1\sim n}=0\)即可得到最大值。
但我們該怎樣把\(x\)的係數都變成負數呢？
其實，只要通過等量代換就可以了。
考慮\(x_j\ge 0\)並不是僅針對\(x_{1\sim n}\)，而同時要求\(x_{n+1\sim n+m}\)也滿足這個條件。
而它們之間存在一定的等量關係，就可以試著通過等量代換來調整所需最大化的式子。
我們規定把\(x_e\)用\(x_{n+l}\)代換的過程叫做\(Pivot(l,e)\)。
已知：
\[x_{n+l}=b_l-\sum_{j=1}^na_{l,j}x_j \]
我們移出右式中的\(x_e\)，將它和\(x_{n+l}\)調換位置得到：
\[x_e=\frac{b_l}{a_{l,e}}-\frac1{a_{l,e}}x_{n+l}-\sum_{j=1}^n[j\not=e]\frac{a_{l,j}}{a_{l,e}}x_{j} \]
然後我們再把其餘式子中的\(x_e\)全部用右式代換，就徹底抹消了\(x_e\)的存在，完成了一次轉軸。
\(Init\)（初始化）
考慮如果存在某個\(b_i<0\)
仔細觀察轉軸操作，可發現這一操作過程中並不會改變\(b_i\)的正負性，因此只要我們在一開始把所有\(b_i\)都強制修改為正數即可。
為了達成這一目標，每次我們找出一個\(b_l<0\)，並找到一個\(a_{l,e}<0\)，然後只要執行\(Pivot(l,e)\)即可。
如果找不到，顯然無解。
\(Simplex\)（單純形法）
單純形法主過程的核心問題，就是每次應該對哪一對\(l,e\)執行轉軸操作：
-
首先，需要滿足\(c_e>0\)，且轉軸後可以讓新的\(c_e<0\)（這是我們執行轉軸操作的初衷）。
-
其次，在所有可選的\(l\)中，我們要使\(\frac{b_l}{a_{l,e}}\)最小（要選擇最緊的限制）。
只要不斷進行\(Pivot(l,e)\)，直到無法繼續就可以了。
程式碼
【UOJ179】線性規劃（據說因為出題人卡常卡精度，只要拿到\(97\)分就可以了？）
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 20
#define DB double
#define eps 1e-8
using namespace std;
int n,m,op;DB a[N+5][N+5];
namespace SimplexMethod//單純形法
{
int id[2*N+5];DB s[N+5];
I void Pivot(CI l,CI e)//轉軸
{
RI i,j;DB t=a[l][e];swap(id[n+l],id[e]);//交換編號
for(a[l][e]=1,i=0;i<=n;++i) a[l][i]/=t;//對這個式子變形
for(i=0;i<=m;++i) if(i^l&&fabs(a[i][e])>eps)//列舉其他式子
for(t=a[i][e],a[i][e]=j=0;j<=n;++j) a[i][j]-=t*a[l][j];//等量代換
}
I bool Init()//初始化，使所有b[i]為正
{
RI i,l,e;W(1)//不斷操作
{
for(l=e=0,i=1;i<=m;++i) a[i][0]<-eps&&(!l||rand()&1)&&(l=i);if(!l) break;//找到b[l]<0
for(i=1;i<=n;++i) a[l][i]<-eps&&(!e||rand()&1)&&(e=i);//找到a[l][e]<0
if(!e) return puts("Infeasible"),0;Pivot(l,e);//找不到則無解，否則轉軸
}return 1;
}
I bool Simplex()//單純形法主過程
{
RI i,l,e;DB Mn;W(1)//不斷操作
{
for(l=e=0,Mn=1e9,i=1;i<=n;++i) if(a[0][i]>eps) {e=i;break;}if(!e) break;//找到c[e]>0
for(i=1;i<=m;++i) a[i][e]>eps&&a[i][0]/a[i][e]<Mn&&(Mn=a[i][0]/a[i][e],l=i);//找到限制最緊的l
if(!l) return puts("Unbounded"),0;Pivot(l,e);//如果所有a[l][e]<0無窮大，否則轉軸
}return 1;
}
I void Solve()//單純形法
{
RI i;for(srand(20050521),i=1;i<=n;++i) id[i]=i;if(Init()&&Simplex())
{
if(printf("%.10lf\n",-a[0][0]),!op) return;//輸出答案
for(i=1;i<=m;++i) s[id[n+i]]=a[i][0];for(i=1;i<=n;++i) printf("%.10lf ",s[i]);//輸出方案
}
}
}
int main()
{
RI i,j;for(scanf("%d%d%d",&n,&m,&op),i=1;i<=n;++i) scanf("%lf",a[0]+i);//讀入
for(i=1;i<=m;++i) {for(j=1;j<=n;++j) scanf("%lf",a[i]+j);scanf("%lf",a[i]);}//讀入
return SimplexMethod::Solve(),0;//單純形法
}