第29章：线性规划

这一章讲解了如何用单纯行算法解决线性规划问题。下面给出用c++实现的代码，输入为标准型线性规划。

先给出单纯行算法的驱动部分，代码如下：

//求解一个输入为标准型线性规划的问题。
double simplex(vector< vector<double> >& A,vector<double>& b,vector<double>& c,vector<double>& solution)
{
int nonBasicNum=A[0].size();
int basicNum=A.size();

vector<int> nonBasic; //nonBasic数组存储的是非基本变量下标的集合；
vector<int> basic; //basic数组存储的是基本变量下标的集合；
double v;

initializeSimplex(nonBasic,basic,A,b,c,v);

int l,e;
basicSolutions(nonBasic,basic,A,b,c,v,basicNum,l,e);

//solution是最终的解
solution.resize(nonBasicNum+1);
for(int i=1;i<=nonBasicNum;++i)
if(isBasic(i,basic))
solution[i]=b[i];
else
solution[i]=0;

return v;//返回的是标准型线性规划的最大目标值。
}

在上述代码中，有两个函数，一个是函数initializeSimplex()，其主要功能是判断输入标准型线性规划问题是否可行，如果不可行，则返回一个消息说明不可行，然后终止，否则，返回一个初始基本解可行的松弛型。还有一个函数basicSolutions()，当线性规划可行时，该函数负责找出此线性规划一系列基本可行解并且进行相应的转动，直至找到最优解。下面给出这两个函数的代码。

函数initializeSimplex:

//输入是标准型线性规划，判断该线性规划是否有可行解。并且如果有可行解，则返回初始解可行的松弛型。
void initializeSimplex(vector<int>& nonBasic,vector<int>& basic,vector< vector<double> >& A,vector<double>& b,vector<double>& c,double& v)
{
int basicNum=A.size();
int nonBasicNum=A[0].size();

nonBasic.resize(nonBasicNum);
for(int i=0;i!=nonBasic.size();++i)
nonBasic[i]=i+1;

basic.resize(basicNum);
for(int i=0;i!=basic.size();++i)
basic[i]=i+1+nonBasicNum;

v=0;

//将数组A扩充为维数为(basicNum+nonBasicNum+1)*(basicNum+nonBasicNum+1)的数组
A.resize(basicNum+nonBasicNum+1);
for(int i=0;i!=A.size();++i)
{
A[i].resize(basicNum+nonBasicNum+1,0);
for(int j=nonBasicNum-1;j>=0;--j)
A[i][j+1]=A[i][j];
A[i][0]=0;
}

for(int i=nonBasicNum+1;i!=A.size();++i)
for(int j=0;j!=A[i].size();++j)
{
A[i][j]=A[i-nonBasicNum-1][j];
A[i-nonBasicNum-1][j]=0;
}

//将数组b扩充为维数为(basicNum+nonBasicNum+1)的数组;
b.resize(basicNum+nonBasicNum+1);
for(int i=b.size()-1;i>=nonBasicNum+1;--i)
{
b[i]=b[i-nonBasicNum-1];
b[i-nonBasicNum-1]=0;
}
//将数组c扩充为维数为(basicNum+nonBasicNum+1)的数组；
c.resize(basicNum+nonBasicNum+1);
for(int i=nonBasicNum-1;i>=0;--i)
c[i+1]=c[i];
c[0]=0;

//找出b中最小的元素及其对应的指标；
double minValue=DBL_MAX;
int minValueIndex;
for(int i=nonBasicNum+1;i!=b.size();++i)
if(b[i]<minValue){
minValue=b[i];
minValueIndex=i;
}

//线性规划基本解是可行的
if(minValue>=0)
return;

//线性规划基本解不可行，构造辅助线性规划；

//将指标0加入到非基本变量指标集合中并将其排序
nonBasic.push_back(0);
sort(nonBasic.begin(),nonBasic.end());

//将A[i][0](i=nonBasicNum+1,...nonBasicNum+basicNum)的值加入到数组A中
for(int i=nonBasicNum+1;i!=A.size();++i)
A[i][0]=-1;

//构造目标函数的系数
vector<double> c1;
c1.resize(nonBasicNum+basicNum+1,0);
c1[0]=-1;

int l=minValueIndex;
int e=0;
pivot(nonBasic,basic,A,b,c1,v,l,e);

basicSolutions(nonBasic,basic,A,b,c1,v,basicNum,l,e);

if(b[0]==0){   //辅助线性规划最优解为0
if(basic[0]==0){     //判断X0是否是基本变量，因为basic是排序后的数组，所有可以用basic[0]=0;
//将X0变成非基本变量
for(int i=1;i!=A[0].size();++i)
if(A[0][i]>0){
e=i;

e496
break;
}
pivot(nonBasic,basic,A,b,c1,v,0,e);
}

for(int i=1;i<=nonBasicNum;++i) //替换原始目标函数，使其只包含非基本变量
if(isBasic(i,basic)){
l=i;
v=v+c[l]*b[l];
for(int j=0;j!=nonBasic.size();++j)
c[nonBasic[j]]=c[nonBasic[j]]-c[l]*A[l][nonBasic[j]];
}

//将非基本变量指标集合中的0删除掉；
for(int i=0;i!=nonBasic.size()-1;++i)
nonBasic[i]=nonBasic[i+1];
nonBasic.pop_back();
}
else
throw runtime_error("the linear-programming problem can't be solved");
}

函数basicSolutions()

//寻找一系列基本可行解并且进行相应的转动；
void basicSolutions(vector<int>& nonBasic,vector<int>& basic,vector< vector<double> >& A,vector<double>& b,vector<double>& c,double& v, int basicNum,int l,int e)
{
vector<double> increment;
increment.resize(basicNum);

while(isIncreaseValue(nonBasic,c,e)){
for(int i=0;i!=basic.size();++i)
if(A[basic[i]][e]>0)
increment[i]=b[basic[i]]/A[basic[i]][e];
else
increment[i]=DBL_MAX;

double minIncreaseValue=DBL_MAX;
int minBasicIndex;
for(int i=0;i!=increment.size();++i)
{
if(increment[i]<minIncreaseValue){
minIncreaseValue=increment[i];
minBasicIndex=i;
}
}

if(minIncreaseValue==DBL_MAX)
throw runtime_error("the linear-programming problem is unbounded");

l=basic[minBasicIndex];
pivot(nonBasic,basic,A,b,c,v,l,e);
}
}

在上面两个函数中出现了pivot()，isIncreaseValue()和isBasic()函数。pivot()函数功能是当我们找出了替入变量和替出变量后执行一次转动，isIncreaseValue()函数是为了在目标函数中找出能够使得目标值增加的具有最小下标的非基本变量，isBasic函数是判断某一个变量是否属于基本变量。这三个函数代码如下：

pivot函数

//将替入变量Xe与替出变量Xl进行转动，输入为松弛型；
void pivot(vector<int>& nonBasic,vector<int>& basic,vector< vector<double> >& A, vector<double>& b,vector<double>& c, double& v,int l,int e)
{
//创建一个跟A维数一样的矩阵；
vector< vector<double> > newA;
newA.resize(A.size());
for(int i=0;i!=newA.size();++i)
newA[i].resize(A[i].size(),0);

//创建一个跟b维数一样的矩阵；
vector<double> newb;
newb.resize(b.size(),0);

//创建一个跟c维数一样的矩阵；
vector<double> newc;
newc.resize(c.size(),0);

//对替出和替入变量所在的那一个束缚条件的b与A做如下相应操作
newb[e]=b[l]/A[l][e];
for(int j=0;j!=nonBasic.size();++j)
if(nonBasic[j]==e)
continue;
else
newA[e][nonBasic[j]]=A[l][nonBasic[j]]/A[l][e];
newA[e][l]=1/A[l][e];

//对其它束缚条件的b和A做如下相应操作；
for(int i=0;i!=basic.size();++i)
if(basic[i]==l)
continue;
else{
newb[basic[i]]=b[basic[i]]-A[basic[i]][e]*newb[e];

for(int j=0;j!=nonBasic.size();++j)
if(nonBasic[j]==e)
continue;
else
newA[basic[i]][nonBasic[j]]=A[basic[i]][nonBasic[j]]-A[basic[i]][e]*newA[e][nonBasic[j]];
newA[basic[i]][l]=-A[basic[i]][e]*newA[e][l];
}

//对目标函数的c和v做如下相应操作;
v=v+c[e]*newb[e];

for(int j=0;j!=nonBasic.size();++j)
if(nonBasic[j]==e)
continue;
else
newc[nonBasic[j]]=c[nonBasic[j]]-c[e]*newA[e][nonBasic[j]];
newc[l]=-c[e]*newA[e][l];

//对非基本变量和基本变量指标集合进行更改；
for(int j=0;j!=nonBasic.size();++j)
if(nonBasic[j]==e){
nonBasic[j]=l;
break;
}

for(int i=0;i!=basic.size();++i)
if(basic[i]==l){
basic[i]=e;
break;
}

//对非基本变量和基本变量指标集合进行排序方便在simplex函数中选择具有最小下标的变量；
sort(nonBasic.begin(),nonBasic.end());//
sort(basic.begin(),basic.end());

//将newA,newb,newc数组中的值复制给A,b,c数组；
for(int i=0;i!=newA.size();++i)
for(int j=0;j!=newA[i].size();++j)
A[i][j]=newA[i][j];

for(int i=0;i!=newb.size();++i)
b[i]=newb[i];

for(int i=0;i!=newc.size();++i)
c[i]=newc[i];
}

isIncreaseValue()函数

bool isIncreaseValue(const vector<int>& nonBasic,const vector<double>& c,int& minIndex) // 在目标函数中找出能够使得目标值增加的具有最小下标的非基本变量
{
for(int j=0;j!=nonBasic.size();++j)
if(c[nonBasic[j]]>0){
minIndex=nonBasic[j];
return true;
}

return false;
}

isBasic()函数

bool isBasic(int index,const vector<int>& basic)//判断某一个变量是否属于基本变量
{
for(int i=0;i!=basic.size();++i)
if(basic[i]==index)
return true;

return false;
}