📄 linesearch_sm.cpp
字号:
/*单纯形法求解线性规划
bug:1 约束条件右端不能出现负值
2 不能用分数，只能用小数
3 默认自变量非负，如果有负值需要平移为正值，加一个正数，最优值不变，但是最优解需要减掉此正数平移量
4 目标函数中常量不被计算在内，如果含有常量，需要计算后自己加上
*/ #include<stdio.h> #include<math.h> #include<iostream.h>
#include <stdlib.h> #include"define.h"
static double matrix[100][100],x[100]; /* 记录总方程的数组,解的数组 */ static int a[100]; /* 记录基础,非基础的解的情况,0:非基础,1:基础 */ static int m,n,s,type; /* 方程变量,约束数,求最大最小值的类型,0:最小 1:最大 */ static int indexe,indexl,indexg; /* 剩余变量,松弛变量,人工变量 */
double sm_min;
double sm_x[NUM_variable];
/*问题1的导函数*/
/*
double* fun1_grads(int dim)
{
double *grads = malloc(sizeof(double)*dim);
for(int i=0;i<dim;i++)//求梯度(dim维的向量)的每个分量值
{
for(int )//对
*grad =
}
return 4*x[0]-2*x[1]-4+2*x[1]*x[1]-6*x[1];
}
*/
void Jckxj() { 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; double 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; double 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; double 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; printf("\n*****************\n [SearchOut]error!!!!!\n*****************\n\n"); return 0; } 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; double 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(double nget[][100],double nlet[][100],double net[][100],double 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(double b[],int code[]) { int i=0,j=0;
#if 0 printf("The number of equator Variables and Restrictors: m n \n"); /* 输入方程变量和约束数 */ cin>>m>>n;
#else
m = 2;
n = 2;
#endif
for(i=0;i<n;i++){
#if 0 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]; /* 输入方程 */
#else
b[i]= 0.0;
code[i]= 0;
#endif }
#if 0
#else
matrix[0][0]=-1.0;
matrix[0][1]=0.0;
matrix[1][0]=0.0;
matrix[1][1]=-1.0;
#endif
#if 0 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);
#else
type= 0;
#endif
#if 0 printf("The Z\n"); /* 输入z */ for(i=0;i<m;i++) cin>>matrix[n][i];
#else//对于可行方向法，线性规划子过程的目标函数为原目标函数在XK点的梯度
double x0[2] = {0.0,0.0};
fun_grads(x0,matrix[n]);
#endif
if(type==1) for(i=0;i<m;i++) matrix[n][i]=-matrix[n][i]; }
void sm_init(double b[],int code[],int M,int N,double B[],double MATRIX[][NUM_variable],double XK[])
{
int i=0,j=0;
//The number of equator Variables and Restrictors: m n
m = M;
n = N;
for(i=0;i<n-1;i++){
b[i]= B[i];//0.0;
code[i]= 0;
}
b[n-1] = 1.0;//最后一个约束条件
code[n-1] = 0;//需要转化为 <= 1
#if 0
for(i=0;i<n;i++){
for(j=0;j<m;j++)
{
matrix[i][j]=MATRIX[i*m+j];
}
}
#else
for(i=0;i<n-1;i++){
for(j=0;j<m;j++)
{
matrix[i][j]=MATRIX[i][j];
}
}
fun_grads(XK,matrix[n-1]);//最后一个约束条件
for(j=0;j<m;j++)//需要转化为 <= 1
{
matrix[n-1][j] = -matrix[n-1][j];
}
#endif
type= 0;
//对于可行方向法，线性规划子过程的目标函数为原目标函数在XK点的梯度
fun_grads(XK,matrix[n]);
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(double c[][100],int row,int vol) { int i; for(i=0;i<n;i++) if(i!=row) c[i][vol]=0; } void Sstart(double b[],int code[]) { int i; double 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(); /* 消去人工变量 */ } bool Simplix() /* 单纯型算法 */ { int in,out,temp=0;
#if 1
int nstep = 0;//
while(nstep <100){
nstep ++;
if(nstep ==100)
{
//PrintResult(); /* 打印最后结果 */
/*存储单纯形法求解线性规划得到的结果：最优解和最优值*/
for(int i=0;i<m;i++)
{
sm_x[i]=x[i];
}
if(type==0)
sm_min=-matrix[n][s];
else
sm_min=matrix[n][s];
//打印单纯形法求dk结果，返回
printf("单纯形法求dk结果：\n dk = (");
for(i = 0 ; i<NUM_variable ; i++)
{
printf(" %f ",sm_x[i]);
}
printf(")");
printf("\n min z = %f",sm_min);
return true;
}
#else while(1){
#endif Jckxj(); /* 基础可行解 */ //Print(); /* 打印 */ //Result(); /* 打印结果 */ if(!Rj()) in=Min(); /* 求换入基 */ else { if(indexg!=0) JustArtificial(); /* 判断人工变量 */ //PrintResult(); /* 打印最后结果 */
/*存储单纯形法求解线性规划得到的结果：最优解和最优值*/
for(int i=0;i<m;i++)
{
sm_x[i]=x[i];
}
if(type==0)
sm_min=-matrix[n][s];
else
sm_min=matrix[n][s];
return true; } if(Check(in)){ /* 判断无界情况 */ printf("No Delimition\n"); return false; } out=SearchOut(&temp,in); /* 求换出基 */ Mto(in,temp); /* 主元化1 */ Be(temp,in); /* 初等变换 */ Achange(in,out); /* 改变a[]的值 */ } }
#if 0 void main() { int code[100]; /* 输入符号标记 */ double b[100]; /* 方程右值 */ Input(b,code); /* 初始化 */ Sstart(b,code); /* 化标准型 */ Simplix(); /* 单纯型算法 */ }
#endif
//The number of equator Variables and Restrictors: m n
void sm(int M,int N,double B[],double MATRIX[][NUM_variable],double XK[])
{
int code[100]; /* 输入符号标记 */
double b[100]; /* 方程右值 */
//Input(b,code); /* 初始化 */
sm_init(b,code, M, N, B, MATRIX,XK);//The number of equator Variables and Restrictors: m n
Sstart(b,code); /* 化标准型 */
Simplix(); /* 单纯型算法 */
}
/*
void main()
{
double B[] = {1.0,2.2};
int num_sm_restrictor =2;
// int NUM_variable = 2;
double A1[2][2] = {1.0,2.2};
double *MATRIX = (double*)malloc(num_sm_restrictor*NUM_variable*sizeof(double));
for(int i=0;i<num_sm_restrictor;i++)
{
for(int j=0;j<NUM_variable;j++)
{
MATRIX[i*NUM_variable+j] = A1[i][j];
}
}
//double MATRIX[2][2] = {1.0,2.2,1.0,2.2};
double XK[2] = {1.0,1.0};
sm(2,2, B, MATRIX,XK);
}
*/
#if 0
for example:
min z = -7x/3 - 13y/3
s.t. x+5y <=0,
-1<=x<=1,
-1<=y<=1
The number of equator Variables and Restrictors: m n
2 3
Input b[] and Restrictor code 0:<= 1:= 2:>=
6 0
The XiShu
1 5
Input b[] and Restrictor code 0:<= 1:= 2:>=
2 0
The XiShu
1 0
Input b[] and Restrictor code 0:<= 1:= 2:>=
2 0
The XiShu
0 1
The Type 0:Min 1:Max
0
The Z
-2.3333333 -4.33333333 6.666666667
X a0 a1 a2 a3 a4 b
X2 1.00 5.00 1.00 0.00 0.00 6.00
X3 1.00 0.00 0.00 1.00 0.00 2.00
X4 0.00 1.00 0.00 0.00 1.00 2.00
Rj -2.33 -4.33 0.00 0.00 0.00 0.00
X2 1.00 5.00 1.00 0.00 0.00 6.00
X3 1.00 0.00 0.00 1.00 0.00 2.00
X4 0.00 1.00 0.00 0.00 1.00 2.00
Rj -2.33 -4.33 0.00 0.00 0.00 0.00
( 0.00 0.00 6.00 2.00 2.00 ) Zmin=0.000000
X1 0.20 1.00 0.20 0.00 0.00 1.20
X3 1.00 0.00 0.00 1.00 0.00 2.00
X4 -0.20 0.00 -0.20 0.00 1.00 0.80
Rj -1.47 0.00 0.87 0.00 0.00 5.20
( 0.00 1.20 0.00 2.00 0.80 ) Zmin=5.200000
X0 0.00 1.00 0.20 -0.20 0.00 0.80
X1 1.00 0.00 0.00 1.00 0.00 2.00
X4 0.00 0.00 -0.20 0.20 1.00 1.20
Rj 0.00 0.00 0.87 1.47 0.00 8.13
( 2.00 0.80 0.00 0.00 1.20 ) Zmin=8.133333
The Minimal :-8.133333
Press any key to continue
#endif⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -