trivival风格的一篇线规博客
单纯形算法
标准型
线性规划是对有限制的最优化线性问题的求解。
人话来说，给出若干个限制方程，其中的所有变量次数都为一次。要求在满足所有限制的情况下，安排各变量的取值，使得一个由这些线性变量组合得到的目标函数取得最值。
更加形式化的给出上述描述
这个形式也被称为线规的标准式，可以写成矩阵的形式
如果不等号的方向与标准型不同，或者要求最小化木匾函数，系数取负就可以了。
如果决策变量$x_i$没有非负限制，将$x_i$拆为两个变量$x_{i1} \leq 0,x_{i2} \leq 0 ,x_i=x_{i1}-x_{i2}$
由于不等式不如等式方便处理，将$\Sigma_j c_{i,j} x_{j} \leqslant b_i$改成$\Sigma_j c_{i,j} x_{j} + f_i = b_i$
$x_i$作主元，参与最优化的决策变量。$f_i$是自由变量。
如果初始有n个决策变量，m个限制方程，那么新扩展的线性方程组就有$n+m$个变量。
算法
如果$n > m$会存在自由的决策变量，答案是无界的。下面只考虑$n \leq m$的情况。
对于上述限制，考虑会有确定的最值的情况。
如果$z=\Sigma_i a_i x_i,a_i\leqslant0$，那么在各个方向上都存在单调性，$x_i$尽可能去满足限制的较小值，最优情况下取到$0$。
由于有辅助变量的存在，我们可以在上述方程组中进行换元，即用一个辅助变量代替主元，成为决策变量。从实现上来说，就是在决策变量的方程中把主元$x_i$用$\Sigma_{j\neq i}v_j x_j + v_t f_t$代替，并且将$f_t$在原$x_i$的方程中系数定为1，调整参数矩阵。这个操作命名为pivot，意为转轴，将主元和自由元转换。
很容易想到一个逼近最优解的办法，尝试不断的将$a_i > 0 $的主元替换，直到所有$a_i \leq 0$。那么该如何选取自由变量$f_t$呢？只有$x_i = v_{i,t} f_t ,v_{i,t} \leq 0$的情况下，换元后新的$a_t$会非正。考虑最优化每次选取$f_t$的收益，我们找$v_{i,t}$最大的$f_t$。
综上可以得到一个算法，每次选一个$a_i>0$的主元，用最大$v_{i,t}\geq 0$的自由元代替。如果不存在$a_i>0$的主元，代表此时已经找到$z$的最优解。如果找不到$v_{i,t}\geq 0$的自由元，那么$x_i$没有上界，答案无界。
至于可能存在的无解情况，考虑在$\Sigma_j c_{i,j} x_{j} \leqslant b_i$中,如果$c_{i,j}>0,b_i<0$，那么这个限制就不可能被满足。
代码
辅助变量并不需要显式的表示在方程中表示，本质就是不等式的常数部分。
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
|
ld ar[M][N];int n,m;
void pivot(int fre,int bas){
ld t=-ar[fre][bas];ar[fre][bas]=-1;
ffor(j,0,n)ar[fre][j]/=t;
ffor(i,0,m)if(fabs(ar[i][bas])>EPS&&i!=fre){
t=ar[i][bas],ar[i][bas]=0;
ffor(j,0,n)ar[i][j]+=t*ar[fre][j];
}
}
void simplex(){
while(1){
int bas=0,fre=0;ld mn=-EPS;
ffor(i,1,m)if(ar[i][0]<mn)mn=ar[fre=i][0];if(!fre)break;
ffor(j,1,n)if(ar[fre][j]>EPS){bas=j;break;}
if(!bas){cout<<"Infeasible";return ;}
pivot(fre,bas);
}
while(1){
int bas=0,fre=0;ld mn=EPS,mx=INF;
ffor(j,1,n)if(ar[0][j]>mn)mn=ar[0][bas=j];if(!bas)break;
ffor(i,1,m)if(ar[i][bas]<-EPS&&(-ar[i][0]/ar[i][bas])<mx)
mx=-ar[i][0]/ar[i][bas],fre=i;
if(!fre){cout<<"Unbounded";return ;}
pivot(fre,bas);
}
}
|
几何意义
这个线性不等式的形式含义是什么呢？观察不等式$\Sigma_j c_{i,j} x_{j} \leqslant b_i$，很自然的想到这是一个半平面的表示。即一个n维超平面将空间分割为两部分，并去下半部分。
那么不等式$CX \leqslant B$就是若干个这样超平面的交，对于向量组$X$，满足不等式的范围也是$X$的合法取值。
再看$z=\Sigma_i a_i x_i$，即$\Sigma_i a_i x_i - z = 0$也是一个超平面，在各维的导数是$a_i$，截距是$a_i/z$。
那么原问题可以理解为，在一个满足限制的半平面交中，平移一个平面，使得该平面始终与半平面交相交，并最大化截距的绝对值。
例题
uoj 线规模板
题目链接:https://uoj.ac/problem/179
题意:给出一组标准型线规，求解或判断线规解的类型
代码:
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
| #include<bits/stdc++.h>
#define ffor(i,a,b) for(int i=a;i<=b;++i)
using namespace std;
using ld=long double;
const int N=2e2+5;
const ld EPS=1e-9;
const ld INF=1e15;
ld ar[N][N];
int id[N];
int n,m;
void pivot(int fre,int bas){
swap(id[fre+n],id[bas]);
ld t=-ar[fre][bas];ar[fre][bas]=-1;
ffor(j,0,n)ar[fre][j]/=t;
ffor(i,0,m)if(fabs(ar[i][bas])>EPS&&i!=fre){
t=ar[i][bas],ar[i][bas]=0;
ffor(j,0,n)ar[i][j]+=t*ar[fre][j];
}
}
int ty;
void simplex(){
ffor(i,1,n) id[i]=i;
while(1){
int bas=0,fre=0;ld mn=-EPS;
ffor(i,1,m)if(ar[i][0]<mn)mn=ar[fre=i][0];if(!fre)break;
ffor(j,1,n)if(ar[fre][j]>EPS){bas=j;break;}
if(!bas){ty=-1;cout<<"Infeasible";return ;}
pivot(fre,bas);
}
while(1){
int bas=0,fre=0;ld mn=EPS,mx=INF;
ffor(j,1,n)if(ar[0][j]>mn)mn=ar[0][bas=j];if(!bas)break;
ffor(i,1,m)if(ar[i][bas]<-EPS&&(-ar[i][0]/ar[i][bas])<mx)
mx=-ar[i][0]/ar[i][bas],fre=i;
if(!fre){ty=-1;cout<<"Unbounded";return ;}
pivot(fre,bas);
}
}
int rec[N];
void solve(){
simplex();
if(ty==-1)return ;
cout<<fixed<<setprecision(10);
cout<<ar[0][0]<<"\n";
if(ty==0)return ;
ffor(i,n+1,n+m)rec[id[i]]=i-n;
ffor(i,1,n){
cout<<(rec[i]?ar[rec[i]][0]:0)<<" ";
}
}
void input(){
cin>>n>>m>>ty;
ffor(j,1,n)cin>>ar[0][j];
ffor(i,1,m){
ffor(j,1,n)cin>>ar[i][j],ar[i][j]*=-1;
cin>>ar[i][0];
}
}
signed main(){
ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
input();
solve();
return 0;
}
|
NOI2008 志愿者招募
题目链接:https://www.luogu.com.cn/problem/P3980
题面:项目需要$n$天才能完成，其中第$i$天至少需要$a_i$个人。布布通过了解得知，一共有$m$类志愿者可以招募,每类志愿者的数量都是无限多的。其中第$i$类可以从第$s_i$天工作到第$t_i$天，招募费用是每人$c_i$元。要求在完成项目的前提下，最小化花费。
题意:$n$天即有$n$个约束方程，$m$种人即$m$个决策变量，记为$x_i$。记每类志愿者选$k_i$人，最小化目标函数$z=\Sigma k_i c_i x_i$。
对于工作时间的限制,可以转化为此人在当前约束方程的系数。问题就被化归到一个标准型线规
代码:
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
| #include<bits/stdc++.h>
#define ffor(i,a,b) for(int i=a;i<=b;++i)
using namespace std;
using ld=double;
const int N=1e3+5,M=1e4+5;
const ld EPS=1e-15;
const ld INF=1e15;
ld ar[M][N],b[M],c[N],ans;
int n,m;
void pivot(int fre,int bas){
ans+=c[bas]*b[fre];
ld t=ar[fre][bas];b[fre]/=t;
ffor(j,1,n)ar[fre][j]/=t;
ffor(j,1,n)if(j!=bas)c[j]-=c[bas]*ar[fre][j];
c[bas]*=-ar[fre][bas];
ffor(i,1,m)if(i!=fre&&fabs(ar[i][bas])>0){
t=ar[i][bas];b[i]-=t*b[fre];
ffor(j,1,n)if(j!=bas)ar[i][j]-=t*ar[fre][j];
ar[i][bas]*=-ar[fre][bas];
}
}
double simplex(){
for(;;){
int bas=0,fre=0;
ffor(i,1,n)if(c[i]>EPS){bas=i;break;}
if(!bas)return ans;
ld mn=INF;
ffor(i,1,m){
if(ar[i][bas]>EPS&&mn>b[i]/ar[i][bas])
mn=b[i]/ar[i][bas],fre=i;
}
if(!fre)return INF;
pivot(fre,bas);
}
}
void solve(){
cout<<int(simplex()+0.5);
}
void input(){
cin>>n>>m;
ffor(i,1,n)cin>>c[i];
ffor(i,1,m){
int l,r;cin>>l>>r>>b[i];
ffor(j,l,r)ar[i][j]=1;
}
}
signed main(){
ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
input();
solve();
return 0;
}
|