您的位置:首页 > 编程语言 > C语言/C++

一般线性规划求最小值+c语言

2015-09-16 15:18 645 查看
求最小值的线性规划,仍然使用旋转算法。

对于如下的线性函数,求最优解:



初始表数据:






算法步骤:





















一个实例:


























c语言程序运行结果:






c语言代码如下:

#include "stdio.h"

#include "stdlib.h"

#include "math.h"

#include "vector"

using namespace std;

#define INF 9999

#define varNum 4 //变量数

#define equNum 0 //等式约束数

#define MAXITER 1000 //最大迭代次数

typedef vector<double> dim1Vector;

typedef vector<dim1Vector> dim2Vector;

dim2Vector getFileINFO(char *File1, char *File2); //获取文本信息

void Cal_LP(dim2Vector info); //解线性方程组

dim2Vector Cal_Spin(dim1Vector C, dim2Vector A, dim1Vector DEV, int laber_i, int laber_j); //旋转计算

void main()

{

char *File1 = "info.txt";

char *File2 = "s.t.txt";

dim2Vector info;

//获取文件信息

info = getFileINFO(File1, File2);

//解线性方程组

Cal_LP(info);

}

//解线性方程组

void Cal_LP(dim2Vector info)

{

int i, j, k;

double sum;

dim1Vector dim1Temp;

dim1Vector DEV; //偏差

dim1Vector C; //目标函数的系数向量

dim1Vector X; //初始基本解

dim2Vector A; //约束条件的系数向量

dim1Vector Col;

dim1Vector Row;

//将info中的数据进行拆分

for(i=0; i<info[0].size(); i++)

{

C.push_back(info[0][i]);

Col.push_back(i);

}

for(i=1; i<info.size()-2; i++)

{

dim1Temp.clear();

for(j=0; j<info[i].size(); j++)

dim1Temp.push_back(info[i][j]);

A.push_back(dim1Temp);

Row.push_back(INF);

}

for(j=0; j<info[i].size(); j++)

X.push_back(info[i][j]);

for(j=0; j<info[info.size()-1].size(); j++)

DEV.push_back(info[info.size()-1][j]);

//初始化偏差

for(i=0; i<A.size(); i++)

{

sum = 0;

for(j=0; j<A[i].size(); j++)

sum += A[i][j]*X[j];

DEV[i] = sum-DEV[i];

}

printf("初始表:\n");

for(i=0; i<C.size(); i++)

printf("%lf ", C[i]);

printf("\n");

for(i=0; i<A.size(); i++)

{

for(j=0; j<A[i].size(); j++)

printf("%lf ", A[i][j]);

printf("%lf\n", DEV[i]);

}

printf("初始解向量:\n");

for(i=0; i<X.size(); i++)

printf("X%d=%lf ", i+1, X[i]);

printf("\n");

printf("\n");

/**********************************************************************/

//开始进入计算

/*********************************************************************/

int count;

int laber_i, laber_j; //记录选择的系数的下标

double maxNum, minNum, tempNum;

dim2Vector Spin_Result;

int flag;

//预处理等式约束

for(i=0; i<equNum; i++)

{

flag = 0;

laber_i = i;

count = 0;

maxNum = -INF;

minNum = INF;

//偏差大于0

if(DEV[i]>0)

{

for(j=0; j<A[i].size(); j++)

if(A[i][j]<0 && A[i][j]/C[j]>maxNum)

{

maxNum = A[i][j]/C[j];

laber_j = j;

count++;

}

if(count==0)

{

printf("线性方程组无解!!!\n");

exit(0);

}

}

//偏差小于0

if(DEV[i]<0)

{

for(j=0; j<A[i].size(); j++)

if(A[i][j]>0 && A[i][j]/C[j]<minNum)

{

minNum = A[i][j]/C[j];

laber_j = j;

count++;

}

if(count==0)

{

printf("线性方程组无解!!!\n");

exit(0);

}

}

//偏差等于0

if(DEV[i]==0)

{

if(flag == 0)

{

for(j=0; j<A[i].size(); j++)

if(A[i][j]>0)

{

flag = 1;

count++;

break;

}

}

if(flag == 0)

{

for(j=0; j<A[i].size(); j++)

if(A[i][j]<0)

{

flag = 2;

count++;

break;

}

}

//情况1

if(flag==1)

{

for(j=0; j<A[i].size(); j++)

if(A[i][j]>0 && A[i][j]/C[j]<minNum)

{

minNum = A[i][j]/C[j];

laber_j = j;

}

}

//情况2

if(flag==2)

{

for(j=0; j<A[i].size(); j++)

if(A[i][j]<0 && A[i][j]/C[j]>maxNum)

{

maxNum = A[i][j]/C[j];

laber_j = j;

}

}

//情况3

if(count==0)

{

printf("线性方程组无解!!!\n");

exit(0);

}

}

tempNum = Row[laber_i];

Row[laber_i] = Col[laber_j];

Col[laber_j] = tempNum;

//旋转计算

Spin_Result = Cal_Spin(C, A, DEV, laber_i, laber_j);

for(k=0; k<Spin_Result.size()-2; k++)

for(j=0; j<Spin_Result[k].size(); j++)

A[k][j] = Spin_Result[k][j];

for(j=0; j<Spin_Result[k].size(); j++)

C[j] = Spin_Result[k][j];

for(j=0; j<Spin_Result[k+1].size(); j++)

DEV[j] = Spin_Result[k+1][j];

}

/***********************************************/

//迭代运算

/***********************************************/

int iter=0;

while(iter<MAXITER)

{

count = 0;

minNum = INF;

for(i=0; i<DEV.size(); i++)

if(DEV[i]<0)

{

sum = 0;

for(j=0; j<A[i].size(); j++)

sum += A[i][j]*A[i][j];

sum = DEV[i]/sqrtf(sum);

if(sum<minNum)

laber_i = i;

count++;

}

//所以偏差均为正数,结束迭代

if(count==0)

break;

count = 0;

minNum = INF;

for(j=0; j<A[laber_i].size(); j++)

{

if(A[laber_i][j]<0)

continue;

if(C[j]/A[laber_i][j]<minNum)

{

minNum = C[j]/A[laber_i][j];

laber_j = j;

}

count++;

}

if(count==0)

{

printf("线性方程组无解!!!\n");

exit(0);

}

tempNum = Row[laber_i];

Row[laber_i] = Col[laber_j];

Col[laber_j] = tempNum;

//旋转计算

Spin_Result = Cal_Spin(C, A, DEV, laber_i, laber_j);

for(k=0; k<Spin_Result.size()-2; k++)

for(j=0; j<Spin_Result[k].size(); j++)

A[k][j] = Spin_Result[k][j];

for(j=0; j<Spin_Result[k].size(); j++)

C[j] = Spin_Result[k][j];

for(j=0; j<Spin_Result[k+1].size(); j++)

DEV[j] = Spin_Result[k+1][j];

//计步器

iter++;

}

for(i=0; i<Row.size(); i++)

if(Row[i]!=INF)

{

X[Row[i]] = X[Row[i]]+DEV[i];

}

printf("\n迭代次数:%d\n", iter);

for(i=0; i<C.size(); i++)

printf("%lf ", C[i]);

printf("\n");

for(i=0; i<A.size(); i++)

{

for(j=0; j<A[i].size(); j++)

printf("%lf ", A[i][j]);

printf("%lf\n", DEV[i]);

}

printf("最优解为:\n");

for(i=0; i<X.size(); i++)

printf("X%d=%lf ", i+1, X[i]);

printf("\n");

sum = 0;

for(i=0; i<X.size(); i++)

sum += info[0][i]*X[i];

printf("最优目标函数值:%lf\n", sum);

}

//旋转计算

dim2Vector Cal_Spin(dim1Vector C, dim2Vector A, dim1Vector DEV, int laber_i, int laber_j)

{

int i, j;

dim1Vector temp_C;

dim2Vector temp_A;

dim1Vector temp_DEV;

dim2Vector dst;

temp_A = A;

temp_C = C;

temp_DEV = DEV;

//计算旋转向量C

C[laber_j] = temp_C[laber_j]/temp_A[laber_i][laber_j];

for(j=0; j<C.size(); j++)

{

if(j==laber_j)

continue;

C[j] = temp_C[j]-(temp_C[laber_j]/temp_A[laber_i][laber_j])*temp_A[laber_i][j];

}

//计算旋转向量DEV

DEV[laber_i] = -temp_DEV[laber_i]/temp_A[laber_i][laber_j];

for(i=0; i<DEV.size(); i++)

{

if(i==laber_i)

continue;

DEV[i] = temp_DEV[i]-(temp_A[i][laber_j]/temp_A[laber_i][laber_j])*temp_DEV[laber_i];

}

//计算旋转矩阵A

A[laber_i][laber_j] = 1/temp_A[laber_i][laber_j];

for(j=0; j<A[laber_i].size(); j++)

{

if(j==laber_j)

continue;

A[laber_i][j] = -temp_A[laber_i][j]/temp_A[laber_i][laber_j];

}

for(i=0; i<A.size(); i++)

{

if(i==laber_i)

continue;

A[i][laber_j] = temp_A[i][laber_j]/temp_A[laber_i][laber_j];

}

for(i=0; i<A.size(); i++)

for(j=0; j<A[i].size(); j++)

{

if(i==laber_i || j==laber_j)

continue;

A[i][j] = temp_A[i][j]-(temp_A[i][laber_j]/temp_A[laber_i][laber_j])*temp_A[laber_i][j];

}

dst = A;

dst.push_back(C);

dst.push_back(DEV);

return dst;

}

//获取文本信息

dim2Vector getFileINFO(char *File1, char *File2)

{

int i=1;

double num;

dim1Vector temp;

dim2Vector dst;

FILE *fp1 = fopen(File1, "r");

FILE *fp2 = fopen(File2, "r");

if(fp1==NULL || fp2==NULL)

{

printf("Open File Error!!\n");

exit(0);

}

while(fscanf(fp1, "%lf", &num)!=EOF)

{

temp.push_back(num);

if(i%varNum==0)

{

dst.push_back(temp);

temp.clear();

}

i++;

}

temp.clear();

while(fscanf(fp2, "%lf", &num)!=EOF)

{

temp.push_back(num);

}

dst.push_back(temp);

return dst;

}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: