單純形法學習筆記
(懶得維護以前的學習筆記的格式了)
用來解決鬆弛型線性規劃,也就是這個形式的線性規劃:
\[\begin{matrix} \text{maximize}&\sum_{j=1}^n c_jx_j&\\ \text{s.t.}&x_{i+n}=b_i-\sum_{j=1}^n a_{i,j}x_j,&i=1,2,\cdots,m\\ &x_j\ge 0,&j=1,2,\cdots,n+m \end{matrix} \]其中 \(x_{i+n}\) 是用來處理小於等於號的,叫做基變數,其他的叫做非基變數。
當 \(b_i\) 非負時可以直接令非基變數為 0 來得到一組可行解,否則後面再說。
得到可行解之後就有一個神奇操作叫轉軸 (pivot) 操作,其實就是把一個非基變數用一個基變數加一堆非基變數再加常數表示出來,然後這個非基變數就變成了基變數,那個基變數就變成非基變數。
比如選擇 \(x_B=b_i-\sum_{j=1}^n a_{i,j}x_j\) 把 \(x_N\) 幹掉,就獲得了 \(x_N={1\over a_{i,N}}(b_i-\sum_{j\ne N}a_{i,j}x_j-x_B)\) ,再把其他等式中的 \(x_N\) 全部換掉,把答案函式中的 \(c_Nx_N\) 也替換掉。
這裡直接給出最優化過程:
- 選擇一個 \(c_e>0\) 的非基變數 \(x_e\)
- 找到滿足 \(a_{l,e}>0\) 且 \(b_l/a_{l,e}\) 最小的基變數 \(x_{n+l}\) ,如果找不到則答案為 \(+\infty\) 。
- 用 \(x_{n+l}\) 替換 \(x_e\) (pivot) ,然後返回第一步。
考慮這個過程在幹些什麼。
- 先找到一個 \(x_e\) ,使得從它這裡可以再增大目標函式。如果所有 \(c_e\) 都非正,那麼此時取所有非基變數為 0 顯然就是最優解。
- 找到 \(\min \{b_l/a_{l,e}\mid a_{l,e}>0\}\) 和對應的 \(l\) 。如果所有 \(a_{l,e}\) 都小於 0 那麼 \(x_e\)
- 替換的時候,由於找到的是最小的 \(b_l/a_{l,e}\) ,所以替換之後仍然滿足 \(b_i\ge 0\) ,取所有非基變數為 0 仍然是一組合法解。替換之後會給目標函式帶來 \(b_l/a_{l,e}\times c_e\) 的收益。
由於某些神奇原因(可能是因為這個演算法在可行解組成的凸包上有實際意義),這個過程一定會終止。複雜度是指數級的,但為了把它卡到指數級需要指數級的權值,所以還是很實用的。 1s 之內衝幾十上百的資料不成問題。
但是好像漏了一件事: \(b_i\) 中存在負數怎麼辦?
構造一個輔助線性規劃
\[\begin{matrix} \text{maximize}&-x_0&\\ \text{s.t.}&x_{i+n}=b_i-\sum_{j=1}^n a_{i,j}x_j+x_0,&i=1,2,\cdots,m\\ &x_j\ge 0,&j=0,1,2,\cdots,n+m \end{matrix} \]這個線性規劃是比較容易搞到 \(b_i\ge 0\) 的。找到最小的 \(b_i\) ,用 \(x_{i+n}\) 替換 \(x_0\) 即可。此時所有 \(b_i\) 非負,可以用上面的方法做一遍。
如果做出來的答案有 \(\max\{-x_0\}<0\) 那麼顯然原線性規劃無解。否則,如果我們強制 \(x_0=0\) ,並且 \(x_0\) 是非基變數,那麼直接把 \(x_0\) 刪去之後得到的這組限制其實與原限制等價,只需要把目標函式換回來即可。
如果 \(x_0\) 是基變數,那麼它對應的 \(b\) 一定是 0 ,所以只需要執行一次轉軸操作把 \(x_0\) 變成非基變數,然後再把它刪掉。
但是這個做法還是有點麻煩。有一種神奇 init 演算法,直接在原來的線性規劃上面操作,每次找到一個 \(b_i<0\) ,此時如果右邊的 \(a_i\) 都非負那麼顯然無解,否則任意找一個負的 \(a_{i,j}\) 進行轉軸操作。不知道為什麼一定會停下。
為了寫起來方便,程式碼裡面的 \(a_{i,j}\) 的意義變為 \(x_{i+n}=b_i+\sum_{j=1}^n a_{i,j}x_j\) ,即取相反數。
#include<bits/stdc++.h>
clock_t __t=clock();
namespace my_std{
using namespace std;
#define pii pair<int,int>
#define fir first
#define sec second
#define MP make_pair
#define rep(i,x,y) for (int i=(x);i<=(y);i++)
#define drep(i,x,y) for (int i=(x);i>=(y);i--)
#define go(x) for (int i=head[x];i;i=edge[i].nxt)
#define templ template<typename T>
#define sz 233
typedef long long ll;
typedef double db;
mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
templ inline T rnd(T l,T r) {return uniform_int_distribution<T>(l,r)(rng);}
templ inline bool chkmax(T &x,T y){return x<y?x=y,1:0;}
templ inline bool chkmin(T &x,T y){return x>y?x=y,1:0;}
templ inline void read(T& t)
{
t=0;char f=0,ch=getchar();double d=0.1;
while(ch>'9'||ch<'0') f|=(ch=='-'),ch=getchar();
while(ch<='9'&&ch>='0') t=t*10+ch-48,ch=getchar();
if(ch=='.'){ch=getchar();while(ch<='9'&&ch>='0') t+=d*(ch^48),d*=0.1,ch=getchar();}
t=(f?-t:t);
}
template<typename T,typename... Args>inline void read(T& t,Args&... args){read(t); read(args...);}
char __sr[1<<21],__z[20];int __C=-1,__zz=0;
inline void Ot(){fwrite(__sr,1,__C+1,stdout),__C=-1;}
inline void print(int x)
{
if(__C>1<<20)Ot();if(x<0)__sr[++__C]='-',x=-x;
while(__z[++__zz]=x%10+48,x/=10);
while(__sr[++__C]=__z[__zz],--__zz);__sr[++__C]='\n';
}
void file()
{
#ifdef NTFOrz
freopen("a.in","r",stdin);
#endif
}
inline void chktime()
{
#ifdef NTFOrz
cout<<(clock()-__t)/1000.0<<'\n';
#endif
}
#ifdef mod
ll ksm(ll x,int y){ll ret=1;for (;y;y>>=1,x=x*x%mod) if (y&1) ret=ret*x%mod;return ret;}
ll inv(ll x){return ksm(x,mod-2);}
#else
ll ksm(ll x,int y){ll ret=1;for (;y;y>>=1,x=x*x) if (y&1) ret=ret*x;return ret;}
#endif
// ll mul(ll x,ll y){ull s=1.0*x/mod*y;ll res=1ull*x*y-s*mod;return (res%mod+mod)%mod;}
}
using namespace my_std;
const db eps=1e-8;
int n,m,tp;
db a[sz][sz]; int id[sz],pos[sz];
void pivot(int B,int N)
{
swap(id[B+n],id[N]);
db w=-1/a[B][N]; a[B][N]=-1; rep(i,0,n) a[B][i]*=w;
rep(i,0,m) if (i!=B) { w=a[i][N],a[i][N]=0; rep(j,0,n) a[i][j]+=w*a[B][j]; }
}
void work()
{
rep(i,1,n) id[i]=i;
while (233)
{
int pos1=0,pos2=0; db w=-eps;
rep(i,1,m) if (chkmin(w,a[i][0])) pos1=i; if (!pos1) break;
rep(i,1,n) if (a[pos1][i]>eps) pos2=i; if (!pos2) return puts("Infeasible"),void();
pivot(pos1,pos2);
}
while (233)
{
int e=0,l=0; db w;
w=eps; rep(i,1,n) if (chkmax(w,a[0][i])) e=i; if (!e) break;
w=1e18; rep(i,1,m) if (a[i][e]<-eps&&chkmin(w,-a[i][0]/a[i][e])) l=i; if (!l) return puts("Unbounded"),void();
pivot(l,e);
}
printf("%.9lf\n",a[0][0]);
rep(i,n+1,n+m) pos[id[i]]=i-n;
if (tp) rep(i,1,n) printf("%.9lf ",pos[i]?a[pos[i]][0]:0);
}
int main()
{
file();
read(n,m,tp);
rep(i,1,n) read(a[0][i]);
rep(i,1,m) { rep(j,1,n) read(a[i][j]),a[i][j]*=-1; read(a[i][0]); }
work();
return 0;
}