This is not the end.

0%

线性规划

Posted on

一些定义

线性函数、线性约束、线性规划问题

已知一组实数$a_1,a_2\cdots a_n$和一组变量$x_1,x_2,\cdots x_n$，在这些变量上的一个线性函数定义为$f(x) = \sum a_ix_i$。

$f(x_1,x_2,\cdots x_n)=b$、$f(x_1,x_2\cdots x_n)\le b$、$f(x_1,x_2\cdots x_n)\ge b$统称为线性约束。

线性规划问题要求最小化或者最大化一个受限于一组有限的线性约束的线性函数。

解

称满足所有限制条件的解$x_1,x_2,\cdots x_n$为可行解，称使目标函数达到最优的解为最优解，称所有可行解构成的区域为解空间。

考虑我们的解空间，它是若干个限制的解空间的交。由于每一个限制的解空间都是一个凸形区域，所以我们的解空间也一定是一个凸形区域。因此，局部最优解只有一个，我们可以用类似爬山的贪心法来求解最优解。

表示方法

标准型

要求最大化$f(x_1,x_2\cdots x_n) = \sum_{i=1}^n c_i x_i$
约束为$\forall i\in [1,m], \sum_{j=1}^n a_{i,j} x_j \le b_i $
且要求$\forall i\in [1,n], x_i \ge 0$

所有的线性规划都可以用标准型描述。

  • 如果取值范围是$(-\infty,+\infty)$，可以把$x$拆成两个取值范围为$[0,+\infty)$的变量$a,b$，令$x=a-b$。
  • 如果限制中是大于等于，或者目标函数要求最小化，可以把所有的系数全部取反。
  • 如果限制中有等式，可以拆成大于等于和小于等于。

标准型也可以用矩阵描述：

最大化$c^Tx$
约束为$A^Tx\le b$且$x\ge 0$

其中向量$x\le y$当且仅当对于每一维都有$x_i \le y_i$

松弛型

主要是为了避免系数为负的时候，不等式变号的问题。

要求最大化$f(x_1,x_2\cdots x_n)= \sum_{i=1}^n c_i x_i$
约束为$\forall i\in [1,m], x_{i+n} = b_i - \sum_{j=1}^n a_{i,j}x_j$
且要求$\forall i\in [1,n+m],x_i\ge 0$

相当于我们强行设了一个大于等于$0$的变量，令它等于$b_i - \sum_{j=1}^n a_{i,j}x_j$。

单纯形法

基本步骤

1）找到一个基本解（如果找不到则无解）。
2）进行转轴操作，使其向解空间内使目标函数增大的方向移动。重复这个步骤直到最优解。

一些定义

松弛型左边的变量称为基变量，右边的变量称为非基变量。初始的时候，$x_1,x_2\cdots x_n$为非基变量，$x_{n+1},x_{n+2},\cdots x_{n+m}$为基变量。

一组基本解的定义是，我们通过一些操作使得所有限制中等式右边的常量都大于等于$0$，这样所有的非基变量取$0$，所有的基变量取右边的常量，就是一组基本解。

转轴操作

转轴操作的本质是，我们增大这个函数中某一个非基变量的值，让它取到它可以取的值的上界。这就要求我们必须把这个变量变成基变量。

首先我们选择这个非基变量，选择的方法是从$c_i>0$的所有的$i$中选择下标最小的$x_i$（如果不强制选择的这个位置的下标最小似乎会死循环？）。如果此时已经没有大于$0$的$c_i$了，这就说明所有的非基变量取$0$是最优的，此时直接返回目标函数的常数项就可以了。

然后我们考虑它最大能够取到多大。我们考虑所有$a_{j,i}>0$的限制，这个上界就应该是$\min \{ {b_j\over a_{j,i}} \}$。如果所有的限制的$a_{j,i}$都小于等于$0$，那么答案是无限大。

设这个式子在$j$处取到最小值。

原先我们有

移项、除系数：

然后我们再把这个式子带入其它的等式。显然由于我们选择了$b_j\over a_{j,i}$最小的$j$，所以带入过后所有的等式的常数项仍然是非负的。

整个操作相当于交换$x_i$和$x_{j+n}$这两个变量，使得$x_{j+n}$变成了非基变量，而$x_i$变成了基变量。

初始化

如果当前还有常数项小于$0$，我们从中随意选择一个等式，然后随机选择这个等式中$a_{i,j}$小于$0$的一个变量，然后对它进行转轴操作。重复这个操作直到所有限制中的常量都非负。

如果我们选择的这个等式中不存在任何一个$a_{i,j}$小于$0$的变量，则无解。

时间复杂度

进行一次转轴操作的复杂度是$O(nm)$的。

而理论上进行转轴操作的上限是${n+m\choose n}$，可以理解为我们决定哪些变量作为基变量、哪些变量作为非基变量的方案数是${n+m\choose n}$。但是通常这个复杂度是远远达不到的。ZJK说期望复杂度大约是在$O(n^2m)$或者$O(nm^2)$左右（？）。

模板 uoj179

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdlib>
#include <ctime>
#define ll long long
#define db long double
using namespace std;
template <class T>
inline void rd(T &x)
{
	x=0; char c=getchar(); int f=1;
	while(!isdigit(c)){if(c=='-')f=-1; c=getchar();}
	while(isdigit(c)) x=x*10-'0'+c,c=getchar(); x*=f;
}
const db eps=1e-10;
const int N=55;
db fabs(db x) { return x>0?x:-x; }
int dcmp(db x) { if(fabs(x)<eps) return 0; return x>0?1:-1; }
db a[N][N];
int n,m,id[N<<1],pos[N<<1],is_inf;
void pivot(int r,int c) {
	swap(id[r+n],id[c]);
	db fr=-a[r][c]; a[r][c]=-1;
	for(int i=1;i<=n+1;++i) a[r][i]/=fr;
	for(int i=1;i<=m+1;++i) if(i!=r) {
		db xi=a[i][c]; if(fabs(xi)<eps) continue; 
		a[i][c]=0;
		for(int j=1;j<=n+1;++j) a[i][j]+=a[r][j]*xi;
	}
}	
db LP() {
	while(1) {
		int p=0;
		for(int i=1;i<=n;++i)
			if(dcmp(a[m+1][i])>0) { p=i; break; }
		if(p==0) return a[m+1][n+1];
		db mi=1e24; int q=0;
		for(int i=1;i<=m;++i) {
			if(dcmp(a[i][p])>=0) continue;
			if(a[i][n+1]/(-a[i][p])<mi) mi=a[i][n+1]/(-a[i][p]),q=i;
		}
		if(q==0) { is_inf=1; return 0; }
		pivot(q,p);
	}
}
bool init() {
	while(1) {
		int r=0,c=0;
		for(int i=1;i<=m;++i) if(dcmp(a[i][n+1])<0&&(!r||rand()&1)) r=i;
		if(!r) return 1;
		for(int i=1;i<=n;++i) if(dcmp(a[r][i])>0&&(!c||rand()&1)) c=i;
		if(!c) return 0;
		pivot(r,c);
	}
}
int main() {
	srand(time(0));
	int ty; rd(n),rd(m),rd(ty);
	for(int i=1;i<=n;++i) scanf("%Lf",&a[m+1][i]);
	for(int i=1;i<=m;++i)
		for(int j=1;j<=n+1;++j) {
			scanf("%Lf",&a[i][j]);
			if(j!=n+1) a[i][j]=-a[i][j];
		}
	for(int i=1;i<=n+m;++i) id[i]=i;
	if(!init()) printf("Infeasible");
	else {
		db Ans=LP();
		if(is_inf) printf("Unbounded");
		else {
			printf("%.8Lf\n",Ans);
			if(ty==0) return 0;
			for(int i=1;i<=n+m;++i) pos[id[i]]=i;
			for(int i=1;i<=n;++i) {
				if(pos[i]<=n) printf("0.0000000 ");
				else printf("%.8Lf ",a[pos[i]-n][n+1]);
			}
		}
	}
	return 0;
}

对偶问题

要求最大化$f(x_1,x_2\cdots x_n) = \sum_{i=1}^n c_i x_i$
约束为$\forall i\in [1,m], \sum_{j=1}^n a_{i,j} x_j \le b_i $
且要求$\forall i\in [1,n] x_i \ge 0$

它的对偶问题是：

要求最小化$f(y_1,y_2\cdots y_n) = \sum_{i=1}^m b_i y_i$
约束为$\forall i\in [1,n], \sum_{j=1}^m a_{j,i} y_j \ge c_i $
且要求$\forall i\in [1,m], y_i \ge 0$

定理1

原问题的最优解等于对偶问题的最优解。

定理2

若$X,Y$分别是原问题和对偶问题的可行解，那么它们都是最优解当且仅当

$\forall i\in [1,n]$，满足$x_i = 0 $或者$\sum a_{j,i}y_j = c_i$。
$\forall i\in [1,m]$，满足$y_i = 0$或者$\sum a_{i,j}x_j = b_i$。

全幺模矩阵

定义：定义一个矩阵为全幺模矩阵，当且仅当它的任意一个子方阵的行列式都是$\pm 1$或者$0$。

如果$A$是全幺模矩阵且$b$中的元素都是整数，那么这个线性规划的最优解，必定可以通过一组整数解达到。

例题

codechef FLYDIST
poj3689 Equations
BZOJ3118 Orz the MST
codechef CHEFBOOK