单纯形法学习笔记

（懒得维护以前的学习笔记的格式了）

用来解决松弛型线性规划，也就是这个形式的线性规划：

\[\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\) 就可以取 \(+\infty\) ，所以答案没有上界。
• 替换的时候，由于找到的是最小的 \(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;
}

标签：非基,ch,变量,单纯形法,int,ll,笔记,学习,rep
来源： https://www.cnblogs.com/p-b-p-b/p/14170326.html