问题
标准型问题：
最大化
∑
j
=
1
n
c
j
x
j
\sum\limits_{j=1}^{n}c_jx_j
j=1∑ncjxj
满足约束
∑
j
=
1
n
a
i
j
⋅
x
j
≤
b
i
(
i
=
1
,
2
,
.
.
.
,
m
)
\sum\limits_{j=1}^{n}a_{ij}\cdot x_j\le b_i\ (i=1,2,...,m)
j=1∑naij⋅xj≤bi (i=1,2,...,m)
且
x
j
≥
0
(
j
=
1
,
2
,
3
,
.
.
.
m
)
x_j\geq 0\ (j=1,2,3,...m)
xj≥0 (j=1,2,3,...m)
可以转换为松弛型问题：
最大化
∑
j
=
1
n
c
j
x
j
\sum\limits_{j=1}^{n}c_jx_j
j=1∑ncjxj
满足约束
x
i
+
n
=
b
i
−
∑
j
=
1
n
a
i
j
⋅
x
j
(
i
=
1
,
2
,
.
.
.
,
m
)
x_{i+n}=b_i-\sum\limits_{j=1}^{n}a_{ij}\cdot x_j\ (i=1,2,...,m)
xi+n=bi−j=1∑naij⋅xj (i=1,2,...,m)
且
x
j
≥
0
(
j
=
1
,
2
,
3
,
.
.
.
n
+
m
)
x_j\geq 0\ (j=1,2,3,...n+m)
xj≥0 (j=1,2,3,...n+m)
把约束中 x i + n x_{i+n} xi+n称为非基变量， x j x_j xj称为基变量
算法
基本思路
初始找出一个满足所有约束的初始解，通过这个解不断的更新，找到更优的解，直到找不到位置（类似爬山算法）。
可以发现这样一定能找到最优解，因为如果把求解的x当做n维空间的一个点，约束条件一定在n维空间中框出一个n维凸“多面体”，一个解为凸“多面体”的一个顶点。凸面体从哪一个发现看都是单峰的，用爬山的算法是肯定能找到最优解的。
例子：
最大化：
S
=
3
x
1
+
x
2
+
2
x
3
S=3x_1+x_2+2x_3
S=3x1+x2+2x3
约束：
x
4
=
30
−
x
1
−
x
2
−
3
x
3
x
5
=
24
−
2
x
1
−
2
x
2
−
5
x
3
x
6
=
36
−
4
x
1
−
x
2
−
2
x
3
x
j
≥
0
(
j
=
1
,
2
,
3
)
\begin{aligned} x_4 &= 30 − x_1 − x_2 − 3x_3\\ x_5 &= 24 − 2x_1 − 2x_2 − 5x_3\\ x_6 &= 36 − 4x_1 − x_2 − 2x_3\\ x_j&\geq 0\ (j=1,2,3) \end{aligned}
x4x5x6xj=30−x1−x2−3x3=24−2x1−2x2−5x3=36−4x1−x2−2x3≥0 (j=1,2,3)
显然可以看出初始解(0,0,0,30,24,36)
把第三个式子化为：
x
1
=
9
−
1
4
x
6
−
1
4
x
2
−
1
2
x
3
x_1=9-\frac 1 4 x_6-\frac 1 4 x_2-\frac 1 2 x_3
x1=9−41x6−41x2−21x3
再将
x
1
x_1
x1代入其它每个式子中去，可以发现
S
=
27
−
3
4
x
6
+
3
4
x
2
+
3
2
x
3
S=27-\frac 3 4 x_6+\frac 3 4 x_2+\frac 3 2x_3
S=27−43x6+43x2+23x3
如此转换，如果我们把S式中所有的x系数转为负数，就可以把常数项当做最优解。
这种把一个非基变量和基变量交换的操作叫做转轴（pivot）。
转轴
转轴操作可以把交换一个非基变量和一个基变量
设交换非基变量
x
N
x_N
xN和基变量
x
B
x_B
xB
x
B
=
b
i
−
∑
j
=
1
n
a
i
j
⋅
x
j
x_B=b_i-\sum\limits_{j=1}^na_{ij}\cdot x_j
xB=bi−j=1∑naij⋅xj
则转为
x
N
=
(
b
i
−
x
B
−
∑
j
=
1
,
j
≠
N
n
a
i
j
⋅
x
j
)
/
a
i
N
x_N=(b_i-x_B-\sum\limits_{j=1,j\ne N}^na_{ij}\cdot x_j)/a_{iN}
xN=(bi−xB−j=1,j̸=N∑naij⋅xj)/aiN
上面的例子转轴
x
1
,
x
6
x_1,x_6
x1,x6之后为
最大化：
S
=
27
−
3
4
x
6
+
3
4
x
2
+
3
2
x
3
S=27-\frac 3 4 x_6+\frac 3 4 x_2+\frac 3 2x_3
S=27−43x6+43x2+23x3
约束：
x
4
=
21
−
3
4
x
6
+
1
4
x
2
+
1
2
x
3
x
5
=
6
+
1
2
x
6
−
3
2
x
2
−
4
x
3
x
1
=
9
−
1
4
x
6
−
1
4
x
2
−
1
2
x
3
x
j
≥
0
(
j
=
1
,
2
,
3
)
\begin{aligned} x_4 &= 21 − \frac 3 4 x_6 + \frac 1 4 x_2 + \frac 1 2 x_3\\ x_5 &= 6 +\frac 1 2 x_6− \frac 3 2 x_2 − 4x_3\\ x_1&=9-\frac 1 4 x_6-\frac 1 4 x_2-\frac 1 2 x_3\\ x_j&\geq 0\ (j=1,2,3) \end{aligned}
x4x5x1xj=21−43x6+41x2+21x3=6+21x6−23x2−4x3=9−41x6−41x2−21x3≥0 (j=1,2,3)
可以发现S至少为27了
simplex操作
即不停的转轴，使得S式中所有的参数都为负。
每次寻找S中参数最大的 x j x_j xj（贪心的想，用最大的参数可以使答案扩得更大），然后从约束条件中，找到对 x j x_j xj限制最大的约束 i i i，即 b i a i j \frac {b_i} {a_{ij}} aijbi最小，转轴 x j x_j xj与 x i + n x_{i+n} xi+n（选择限制最大的约束，使得转轴后，其余约束条件中的 b i b_i bi非负）
寻找初始解
当满足所有 b i ≥ 0 b_i\geq 0 bi≥0，显然可以把 x j = 0 , x i + n = b i x_j=0,x_{i+n}=b_i xj=0,xi+n=bi当做一个可行的初始解。
如果存在 b i < 0 b_i<0 bi<0，有两种方法
法一
论文方法：
可以通过引入辅助变量
x
0
x_0
x0 ，经过转轴操作将
b
i
b_i
bi转换为
≥
0
\geq 0
≥0。
最大化
−
x
0
-x_0
−x0
满足约束
x
i
+
n
=
b
i
−
∑
j
=
1
n
a
i
j
⋅
x
j
+
x
0
(
i
=
1
,
2
,
.
.
.
,
m
)
x_{i+n}=b_i-\sum\limits_{j=1}^{n}a_{ij}\cdot x_j+x_0\ (i=1,2,...,m)
xi+n=bi−j=1∑naij⋅xj+x0 (i=1,2,...,m)
如果得到最优解
x
0
=
0
x_0=0
x0=0，则这个
x
0
x_0
x0对原来问题的答案没有影响，否则原问题无解。
选择最小的
b
i
b_i
bi，设为第k个约束，把那个约束中
x
0
x_0
x0进行转轴操作换出来，则
x
0
=
−
b
k
+
x
k
+
n
+
∑
j
=
1
n
a
k
j
⋅
x
j
x_0=-b_k+x_{k+n}+\sum\limits_{j=1}^{n}a_{kj}\cdot x_j
x0=−bk+xk+n+j=1∑nakj⋅xj
其它约束变为
x
i
+
n
=
b
i
−
b
k
+
x
k
+
n
+
∑
j
=
1
n
(
a
k
j
−
a
i
j
)
⋅
x
j
x_{i+n}=b_i-b_k+x_{k+n}+\sum\limits_{j=1}^{n}(a_{kj}-a_{ij})\cdot x_j
xi+n=bi−bk+xk+n+j=1∑n(akj−aij)⋅xj
因为
b
k
b_k
bk为最小的b，所以这个式子中
b
i
−
b
k
b_i-b_k
bi−bk一定为非负数，从而把所有的约束中常数项变为非负。
此时再从要求最大化的式子中，选一个参数为正的 x p x_p xp，将其与 x 0 x_0 x0转轴。
重复这样的操作，直到所有 b i b_i bi全为正。
实现简单的法二
每次随机找到一个<0的
b
i
b_i
bi，在从
b
i
b_i
bi的约束条件中找一个系数大于0的
x
j
x_j
xj，转轴
x
i
+
n
,
x
j
x_{i+n},x_j
xi+n,xj，就可以把当前
b
i
b_i
bi变为正的（虽然有可能把其它条件的
b
i
b_i
bi变为负数，所以要随机。。。）
多次执行随机操作，基本上可以把所有
b
i
b_i
bi变为非负数。（UOJ拿97分）
法三
把法二的 b i b_i bi，每次选择最小的负数 b i b_i bi来进行操作也能过，并不知道为什么。。。
代码
UOJ179
#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
#include<algorithm>
using namespace std;
const int MAXN=55;
const double EPS=1e-8;
int dcmp(double a,double b)
{
if(fabs(a-b)<=EPS)
return 0;
return a<b?-1:1;
}
int n,m,type;
double A[MAXN][MAXN];
int id[MAXN],tp[MAXN];
void pivot(int r,int c)
{
swap(id[r+n],id[c]);
double tmp=-A[r][c];
A[r][c]=-1;
for(int j=0;j<=n;j++)
A[r][j]/=tmp;
for(int i=0;i<=m;i++)
if(i!=r&&dcmp(A[i][c],0)!=0)
{
tmp=A[i][c];
A[i][c]=0;
for(int j=0;j<=n;j++)
A[i][j]+=tmp*A[r][j];
}
}
void solve()
{
for(int i=1;i<=n;i++)
id[i]=i;
for(;;)
{
int i=0,j=0;
for(int k=1;k<=m;k++)
if(dcmp(A[k][0],0)==-1&&(!i||rand()&1))
i=k;
if(i==0)
break;
for(int k=1;k<=n;k++)
if(dcmp(A[i][k],0)==1&&(!j||rand()&1))
j=k;
if(j==0)
{
puts("Infeasible");
return;
}
pivot(i,j);
}
for(;;)
{
int i=0,j=0;
double mx=0,mn=1e9;
for(int k=1;k<=n;k++)
if(dcmp(A[0][k],mx)==1)
mx=A[0][k],j=k;
if(j==0)
break;
for(int k=1;k<=m;k++)
if(dcmp(A[k][j],0)==-1&&dcmp(-A[k][0]/A[k][j],mn)==-1)
mn=-A[k][0]/A[k][j],i=k;
if(i==0)
{
puts("Unbounded");
return;
}
pivot(i,j);
}
printf("%.9f\n",A[0][0]);
for(int i=n+1;i<=n+m;i++)
tp[id[i]]=i-n;
if(type)
for(int i=1;i<=n;i++)
printf("%.9f%c",tp[i]?A[tp[i]][0]:0," \n"[i==n]);
}
int main()
{
scanf("%d%d%d",&n,&m,&type);
for(int i=1;i<=n;i++)
scanf("%lf",&A[0][i]);
for(int i=1;i<=m;i++)
{
for(int j=1;j<=n;j++)
{
scanf("%lf",&A[i][j]);
A[i][j]*=-1;
}
scanf("%lf",&A[i][0]);
}
solve();
return 0;
}
6489
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
