注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

fanhq666的博客

Fan-Fun

 
 
 

日志

 
 

Day2《读心》题目的线性规划解法  

2011-02-17 12:52:12|  分类: 程序 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
先回忆一下什么是线性规划问题:

找出一组非负实数X[0..M-1],使
最大化:
         Z[0]*X[0]+Z[1]*X[1]+...+Z[M-1]*X[M-1]
满足约束
        A[0][0]*X[0]+A[0][1]*X[1]+...+A[0][M-1]*X[M-1] <=R[0]
        A[1][0]*X[0]+A[1][1]*X[1]+...+A[1][M-1]*X[M-1] <=R[1]
        ...
        A[N-1][0]*X[0]+A[N-1][1]*X[1]+...+A[N-1][M-1]*X[M-1] <=R[N-1]
其中,最大化也可改为最小化(将Z[i]取相反数即可),小于等于也可改成大于等于(等式两边同乘-1即可),X[i]也可以没有非负约束(用X[i]=X'[i]-X''[i]代换即可)

容易证明,在《读心》这个问题中,是没有可能ExpScore>0的,于是,我们的目标就是让对手所有的决策下我们的收益非负。
其实就是要解决这样一个问题:
寻找一组非负实数P[i],满足
           sigma(P[i])=1
           for i=0..N-1,   sigma(A[j][i]*P[j])>=0
这个可以通过一些技巧转化为线性规划问题:
          最大化  sigma(P[i])
          满足:for i=0..N-1,sigma(-A[j][i]*P[j])<=0
                    sigma(P[i])<=1
易知,计算出的最大值一定是1,求出的P[0..M-1]全部乘以100000即为输出。

在实践中,为了线性规划的速度和稳定性,一般把第一行的约束条件的右边写成<=0.00001.

线性规划有一个“古典”的方法:单纯形法。详见算法导论。
它的核心思想:不断调整当前的可行解。
例(来自算法导论):
最大化3x1+x2+2x3,满足
x1+x2+3x3<=30
2x1+2x2+5x3<=24
4x1+x2+2x3<=36
(x1,x2,x3均非负)
为每个不等式增加一个变量,使其变为等式;
最大化
z=3x1+x2+2x3,
满足
x4=30-x1-x2-3x3
x5=24-2x1-2x2-5x3
x6=36-4x1-x2-2x3
(x1~x6非负)
在等式左边的x4~x6为“基本变量”,右边的为“非基本变量”。
把所有非基本变量设为0显然是一组解(对于不是的情况,我们会讨论处理的方法),我们称之为“基本解”。
目前的基本解得到的z的值是0.我们想办法让一个基本变量和非基本变量互换地位,使得它产生的基本解的z上升。
例如,改写第一个式子,让x1来到左边,
x1=9-x2/4-x3/2-x6/4
代入其他式子,得到:
最大化
z=27+x2/4+x3/2-3x6/4
满足
x1=9-x2/4-x3/2-x6/4
x4=21-3x2/4-5x3/2+x6/4
x5=6-3x2/2-4x3+x6/2
现在的基本变量是x1,x4和x5,非基本变量是x2,x3,x6,非基本变量都是零的情况下z值为27,可见,目标值上升了。
改写第3个式子,让x4来到左边,并代入其他式子
最大化
z=111/4+x2/16-x5/8-11x6/16
满足
x1=33/4-x2/16+x5/8-5x6/16
x3=3/2-3x2/8-x5/4+x6/8
x4=69/4+3x2/16+5x5/8-x6/16
继续改写第2个式子,让x2来到左边,并代入其他式子
最大化
z=28-x3/6-x5/6-2x6/3
满足
x1=8+x3/6+x5/6-x6/3
x2=4-8x3/3-2x5/3+x6/3
x4=18-x3/2+x5/2
因为x3>=0,x5>=0,x6>=0,所以z=28-x3/6-x5/6-2x6/3<=28。于是,我们得到了最优解:
z=28。令非基本变量均为0,得到x1=8,x2=4,x3=0。这就是这个例子的答案。

在每一步中,如何选择需要改写的式子以及让哪个变量来到左边是有一套精细的法则的。具体的细节可以参看代码。

对于把x1..xM均设为0而不满足非负性要求的情况,
例如
z=2x1-x2,
2x1-x2<=2
x1-5x2<=-4
我们通过增加一个“人工变量”x0,建立一个“辅助线性规划”
最大化
Z=-x0
满足
x3=2-2x1+x2+x0
x4=-4-x1+5x2+x0
将第二个式子改写,x0换到左边,有
Z=-4-x1+5x2-x4
x0=4+x1-5x2+x4
x3=6-x1-4x2+x4
现在基本解是合法的了。像往常那样进行变形,可得:
Z=-x0
x2=4/5-x0/5+x1/5+x4/5
x3=14/5+4x0/5-9x1/5+x4/5
如果问题有解,那么Z的最优值一定是0.此时x0也是0.删除所有x0,把Z换成z,得到
z=2x1-x2
x2=4/5-x0/5+x1/5+x4/5
x3=14/5+4x0/5-9x1/5+x4/5
把基本变量的表达式代入z,得到正规的形式:
z=4/5+9x1/5-x4/5
x2=4/5+x1/5+x4/5
x3=14/5-9x1/5+x4/5
现在,终于可以按照普通的线性规划的方法继续求解了。

代码的实现是需要技巧的。我写了一个版本,和以往写的不是很一样。
http://cid-354ed8646264d3c4.office.live.com/self.aspx/.Public/LinearProgramming.cpp


#include <cstdio>
#include <vector>
using namespace std;
double Abs(double a){return a>=0?a:-a;}
namespace LP{
vector<vector<double> >A;//系数矩阵A
vector<double> R,Z,X;//R表示等式中的常数项,Z表示z的系数,X用于存储最终答案中x的值
double Now;//Now记录当前z的表达式的常数项。当运行结束时,Now的值就是答案
vector<int> I,P;//I用于标记第i个式子的左边的变量是谁,P用于标记第x[i]在那个等式的左边(如果它是基本变量的话)
int N,M,nn,mm;
void subtract(vector<double> &a,double &r,vector<double> &b,double s,int c){//一个式子减另一个式子
double o=a[c];
for (int i=0;i<mm;i++)a[i]-=o*b[i];
a[c]=0;
r-=o*s;
}
void exchange(int a,int b){//改写第a行的式子,把x[b]换到左边,并代入
P[I[a]]=-1;I[a]=b;P[b]=a;
double o=A[a][b];
for (int i=0;i<mm;i++)A[a][i]/=o;
A[a][b]=1;
R[a]/=o;
for (int i=0;i<nn;i++)if (a!=i)subtract(A[i],R[i],A[a],R[a],b);
subtract(Z,Now,A[a],-R[a],b);
}
int PIVOT(){//挑选用于改写的式子和变量
int u=-1;
                //挑选x[u]作为等式左边的变量,要求是它在z的表达式中的系数为正。一般来说,挑选z中系数最大的那个
int c=rand()%2;//一个小技巧,防止某种特殊条件下的死循环
for (int i=0;i<mm;i++)if (P[i]==-1 && Z[i]>1e-9)
if (u==-1 || (c?Z[i]>Z[u]:0))u=i;
if (u==-1)return -1;
int k=-1;
                //挑选第k行为等待改写的式子,要求A[k][u]>0,且改写它之后所有的R[i]都是非负的。这个通过挑选最小的R[i]/A[i][u]来实现。
double bound;
for (int i=0;i<nn;i++)if (A[i][u]>1e-9){
double h=R[i]/A[i][u];
if (k==-1 || bound>h)k=i,bound=h;
}
if (k==-1)return -2;
exchange(k,u);
return 0;
}
int LP(int n,int m,vector<double> z,vector<vector<double> > a,vector<double> r){//LP的主过程
M=m;N=n;
                mm=N+M+1;//M个原有的变量,N个“松弛变量”,和一个人工变量,共计mm个
                nn=N;
X.assign(M,0);
A.assign(nn,vector<double>(mm,0));
for (int i=0;i<N;i++)for (int j=0;j<M;j++)A[i][j]=a[i][j];
R=r;
I.assign(nn,0);
P.assign(mm,-1);
                //构造各个等式
for (int i=0;i<nn;i++){
A[i][M+i]=1;A[i][mm-1]=-1;
I[i]=M+i;P[M+i]=i;
}
                //准备进行辅助线性规划
Z.assign(mm,0);
Z[mm-1]=-1;
Now=0;
int u=0;
for (int i=1;i<nn;i++)if (R[i]<R[u])u=i;
if (R[u]<0){//真的必须进行辅助线性规划了
exchange(u,mm-1);//通过改写令基本解合法
while (PIVOT()==0);//增广使x0变成0
if (Now<-1e-9)return -1;//No Solution
if (P[mm-1]!=-1){
u=-1;
for (int i=0;i<mm;i++)if (P[i]==-1)
if (u==-1 || Abs(A[P[mm-1]][i])>Abs(A[P[mm-1]][u]))u=i;
exchange(P[mm-1],u);
}
}
mm--;//删除人工变量
                //构造z的表达式
for (int i=0;i<M;i++)Z[i]=z[i];
for (int i=M;i<mm;i++)Z[i]=0;
Now=0;
for (int i=0;i<M;i++)if (P[i]!=-1)
subtract(Z,Now,A[P[i]],-R[P[i]],i);


                //这个才是LP的主过程
while ((u=PIVOT())==0);


if (u==-2)return -2;//Unbounded
for (int i=0;i<nn;i++)if (I[i]<M)X[I[i]]=R[i];
return 0;
}
}
int main()
{
int n,m;
scanf("%d%d",&n,&m);
vector<vector<double> > a(n,vector<double>(m,0));
vector<double> r(n),z(m);
for (int i=0;i<m;i++)scanf("%lf",&z[i]);
for (int i=0;i<n;i++){
for (int j=0;j<m;j++)scanf("%lf",&a[i][j]);
scanf("%lf",&r[i]);
}
int t=LP::LP(n,m,z,a,r);
if (t==-1)puts("No Solution");
else if (t==-2)puts("Un Bounded");
else{
printf("%f\n",LP::Now);
for (int i=0;i<m;i++)printf("%f ",LP::X[i]);
puts("");
}
return 0;
}

  评论这张
 
阅读(1130)| 评论(2)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017