这是我从网上找的一个求一般线性规划方程单纯形法的程序
可以解决一般线性规划的程序,可以运行,但结果有问题,有没有大神了解这个算法的帮帮我改正确。可以谈价格,拜托了。
#include <stdio.h>
#include <math.h>
#include
using namespace std;
float matrix[100][100],x[100]; /* 记录总方程的数组,解的数组 /
int a[100]; / 记录基础,非基础的解的情况,0:非基础,1:基础 /
int m,n,s,type; / 方程变量,约束数,求最大最小值的类型,0:最小 1:最大 /
int indexe,indexl,indexg; / 剩余变量,松弛变量,人工变量 */
void Jckxj() /*basic feasible solution */
{
int i,j;
for(i=0;i <n;i++)
for(j=0;j <s;j++)
if(matrix[i][j]==1&&a[j]==1)
{
x[j]=matrix[i][s];
j=s;
}
for(i=0;i <s;i++)
if(a[i]==0)
x[i]=0;
}
int Rj()
{
int i;
for(i=0;i <s;i++)
if(fabs(matrix[n][i])>=0.000001)
if(matrix[n][i] <0)
return 0;
return 1;
}
int Min()
{
int i,temp=0;
float min=matrix[n][0];
for(i=1;i <s;i++)
if(min> matrix[n][i])
{
min=matrix[n][i];
temp=i;
}
return temp;
}
void JustArtificial()
{
int i;
for(i=m+indexe+indexl;i <s;i++)
if(fabs(x[i])>=0.000001)
{
printf( "No Answer\n ");
return;
}
}
int Check(int in)
{
int i;
float max1=-1;
for(i=0;i <n;i++)
if(fabs(matrix[i][in])>=0.000001&&max1 <matrix[i][s]/matrix[i][in])
max1=matrix[i][s]/matrix[i][in];
if(max1 <0)
return 1;
return 0;
}
int SearchOut(int *temp,int in)
{
int i;
float min=10000;
for(i=0;i <n;i++)
if(fabs(matrix[i][in])>=0.000001&&(matrix[i][s]/matrix[i][in]>=0)&&min> matrix[i][s]/matrix[i][in])
{
min=matrix[i][s]/matrix[i][in];
*temp=i;
}
for(i=0;i <s;i++)
if(a[i]==1&&matrix[*temp][i]==1)
return i;
}
void Mto(int in,int temp)
{
int i;
for(i=0;i <=s;i++)
if(i!=in)
matrix[temp][i]=matrix[temp][i]/matrix[temp][in];
matrix[temp][in]=1;
}
void Be(int temp,int in)
{
int i,j;
float c;
for(i=0;i <=n;i++)
{
c=matrix[i][in]/matrix[temp][in];
if(i!=temp)
for(j=0;j <=s;j++)
matrix[i][j]=matrix[i][j]-matrix[temp][j]*c;
}
}
void Achange(int in,int out)
{
int temp=a[in];
a[in]=a[out];
a[out]=temp;
}
void Print()
{
int i,j,k,temp=0;
for(i=0;i <n;i++)
{
for(k=temp;k <s;k++)
if(a[k]==1)
{
printf( "X%d ",k);
temp=k+1;
k=s;
}
for(j=0;j <=s;j++)
printf( "%8.2f ",matrix[i][j]);
printf( "\n ");
}
printf( "Rj ");
for(j=0;j <=s;j++)
printf( "%8.2f ",matrix[n][j]);
printf( "\n ");
}
void InitPrint()
{
int i;
printf( "X ");
for(i=0;i <s;i++)
printf( " a%d ",i);
printf( " b\n ");
Print();
printf( "\n ");
}
void Result()
{
int i;
printf( " ( ");
for(i=0;i <s;i++)
printf( "%8.2f ",x[i]);
printf( " ) ");
if(type==1)
printf( " Zmax=%f\n\n ",matrix[n][s]);
else
printf( " Zmin=%f\n\n ",matrix[n][s]);
}
void PrintResult()
{
if(type==0)
printf( "The Minimal :%f\n ",-matrix[n][s]);
else
printf( "The Maximum :%f\n ",matrix[n][s]);
}
void Merge(float nget[][100],float nlet[][100],float net[][100],float b[])
{
int i,j;
for(i=0;i <n;i++)
{
for(j=m;j <m+indexe;j++)
if(nget[i][j-m]!=-1)
matrix[i][j]=0;
else
matrix[i][j]=-1;
for(j=m+indexe;j <m+indexe+indexl;j++)
if(nlet[i][j-m-indexe]!=1)
matrix[i][j]=0;
else
matrix[i][j]=1;
for(j=m+indexe+indexl;j <s;j++)
if(net[i][j-m-indexe-indexl]!=1)
matrix[i][j]=0;
else
matrix[i][j]=1;
matrix[i][s]=b[i];
}
for(i=m;i <m+indexe+indexl;i++)
matrix[n][i]=0;
for(i=m+indexe+indexl;i <s;i++)
matrix[n][i]=100;
matrix[n][s]=0;
}
void ProcessA()
{
int i;
for(i=0;i <m+indexe;i++)
a[i]=0;
for(i=m+indexe;i <s;i++)
a[i]=1;
}
void Input(float b[],int code[])
{
int i=0,j=0;
printf( "The equator Variable and Restrictor\n "); /* 输入方程变量和约束数 /
cin>> m>> n;
for(i=0;i <n;i++)
{
printf( "Input b[] and Restrictor code 0: <= 1:= 2:> =\n "); / 输入方程右边的值,code的值 /
cin>> b[i]>> code[i];
printf( "The XiShu\n ");
for(j=0;j <m;j++)
cin>> matrix[i][j]; / 输入方程 /
}
printf( "The Type 0:Min 1:Max \n "); / 输入求最大值还是最小值 /
do
{
cin>> type;
if(type!=0&&type!=1)
printf( "Error,ReInput\n ");
}
while(type!=0&&type!=1);
printf( "The Z\n "); / 输入z */
for(i=0;i <m;i++)
cin>> matrix[n][i];
if(type==1)
for(i=0;i <m;i++)
matrix[n][i]=-matrix[n][i];
}
void Xartificial()
{
int i,j,k;
if(indexg!=0)
{
for(i=m+indexe+indexl;i <s;i++)
{
for(j=0;j <n;j++)
if(matrix[j][i]==1)
{
for(k=0;k <=s;k++)
matrix[n][k]=matrix[n][k]-matrix[j][k]*100;
j=n;
}
}
}
}
void Process(float c[][100],int row,int vol)
{
int i;
for(i=0;i <n;i++)
if(i!=row)
c[i][vol]=0;
}
void Sstart(float b[],int code[])
{
int i;
float nget[100][100],nlet[100][100],net[100][100]; /* 剩余变量数组,松弛变量数组,人工变量数组 /
indexe=indexl=indexg=0;
for(i=0;i <n;i++)
{
if(code[i]==0)
{
nlet[i][indexl++]=1;
Process(nlet,i,indexl-1);
}
if(code[i]==1)
{
net[i][indexg++]=1;
Process(net,i,indexg-1);
}
if(code[i]==2)
{
net[i][indexg++]=1;
nget[i][indexe++]=-1;
Process(net,i,indexg-1);
Process(nget,i,indexe-1);
}
}
s=indexe+indexl+indexg+m;
Merge(nget,nlet,net,b); / 合并 /
ProcessA(); / 初始化a[] /
InitPrint(); / 初始化打印 /
Xartificial(); / 消去人工变量 */
}
void Simplix() /* 单纯型算法 /
{
int in,out,temp=0;
while(1)
{
Jckxj(); / 基础可行解 /
Print(); / 打印 /
Result(); / 打印结果 /
if(!Rj())
in=Min(); / 求换入基 /
else
{
if(indexg!=0)
JustArtificial(); / 判断人工变量 /
PrintResult(); / 打印最后结果 /
return;
}
if(Check(in))
{ / 判断无界情况 /
printf( "No Delimition\n ");
return;
}
out=SearchOut(&temp,in); / 求换出基 /
Mto(in,temp); / 主元化1 /
Be(temp,in); / 初等变换 /
Achange(in,out); / 改变a[]的值 */
}
}
int main()
{
int code[100]; /* 输入符号标记 /
float b[100]; / 方程右值 /
Input(b,code); / 初始化 /
Sstart(b,code); / 化标准型 /
Simplix(); / 单纯型算法 */
}