单纯形法原理
线性规划问题中,若存在最优解那么一定在某个基本可行解处取得,而我们可以将所有基本可行解找出来,代入目标函数后找到最优解,但是随着问题规模的增大,寻找所有基本可行解就变得异常困难。1947年G.B.Bantzig提出了单纯形法,成为了求解线性规划问题的基本方法。即先找到一个基本可行解,判断是否为最优解,若不是,从该可行解出发,寻找另一个基本可行解,再进行最优性判断,经此有限次步骤后便能够找到LP的最优解。
参考文献
[1]于松南.长春工业大学运筹学课程设计[Z]. https://wenku.baidu.com/view/4380168552ea551811a6872f.html,2015.
[2]尚有林.运筹学(第一版)[M]. 北京:高等教育出版社,2001.
Code
#include <stdio.h>
#include <math.h>
#pragma warning(disable:4996)
#define m 3 //定义约束方程个数
#define n 5 //定义决策变量个数
double M = 1e6;
double A[m][n]; //存储A约束矩阵
double C[n]; //存储目标函数价值系数C矩阵
double b[m]; //存储约束条件资源系数b矩阵
double CB[m]; //存储基变量系数矩阵
double seta[m]; //存放出基和入基的变化情况
double sigma[n]; //存储检验数矩阵σ
double x[n]; //存储决策变量
int num[m]; //存储出基和入基变量的情况
double Z = 0; //记录目标函数值
void input();
void print();
int FindSwapInVar(); //确定入基变量
int FindSwapOutVar(int a); //确定出基变量
void Iterate(int a, int b); //换基迭代
//确定入基变量,对于所有检验数均大于等于0时,当前解为最优解
int FindSwapInVar()
{
int i, k = 0;//用k带出入基变量
int flag = 1;
double min = 0;
for (i = 0; i < n; i++)
{
if (sigma[i] < 0)
{
flag = 0; break;
}
}
if (flag == 1)//找到了咯
return -1;
for (i = 0; i < n; i++)//找出最大的检验数σ
{
if (sigma[i] < min)//求min问题按照sigma最小的为入基变量
{
min = sigma[i];
k = i;
}
}
return k;
}
//确定出基变量
int FindSwapOutVar(int a)
{
int i, j;
int flag = 1;
int k = a;//a入基变量
for (i = 0; i < m; i++)//如果某个数小于0的检验数,对应的列向量中所有元素≤0 该问题有无界解
{
if (A[i][k] > 0)
{
flag = 0; break;
}
}
if (flag == 1)
{
printf("该线性规划问题有无界解、无最优解!\n");
return -1;
}
for (i = 0; i < m; i++)
{
if (A[i][k] > 0)
seta[i] = b[i] / A[i][k];
else seta[i] = M;//当A[i][k]≤0的时候是不需要考虑的
//然而是根据比值最小原则整的 所以给对应的卡一个很大的M
}
double min = M;
for (i = 0; i < m; i++)//得到换出变量
{
if (min >= seta[i])
{
min = seta[i];
j = i;
}
}
num[j] = k + 1;
CB[j] = C[k];
return j;
}
//迭代运算,计算新的单纯形表
void Iterate(int p, int q)
{
int i, j, r, c;//row行 column列(r,l)就是转轴元
r = p;//行号 p是出基变量
c = q;//列号 q是入基变量
double temp1 = A[r][c];
double temp2, temp3;
//标准化该行
b[r] /= temp1;
for (j = 0; j < n; j++)
A[r][j] /= temp1;
for (i = 0; i < m; i++)//标准化其他行
{
if (i != r)
if (A[i][c] != 0)
{
temp2 = A[i][c];
b[i] -= temp2 * b[r];//b[r]转轴元对应b
for (j = 0; j < n; j++)
A[i][j] -= temp2 * A[r][j];//A[r][j]是转轴元对应行
}
}
//σ的迭代计算
temp3 = sigma[c];
for (i = 0; i < n; i++)
{
sigma[i] -= A[r][i] * temp3;
}
}
//输入函数、约束条件的技术系数矩阵A、初始基变量的数字代码
//输入约束条件的资源系数b、目标函数价值系数C
void input()
{
int i, j;
printf("请按顺序输入方程组的系数矩阵A(共%d行%d列):\n", m, n);
for (i = 0; i < m; i++)
for (j = 0; j < n; j++)
scanf("%lf", &A[i][j]);
printf("请按顺序输入初始基变量所在列数的数字代码num矩阵:\n");
for (i = 0; i < m; i++)
scanf("%d", &num[i]);
printf("请按顺序输入约束条件的资源系数b矩阵:\n");
for (i = 0; i < m; i++)
scanf("%lf", &b[i]);
printf("请按顺序输入目标函数的价值系数C矩阵:\n");
for (i = 0; i < n; i++)
scanf("%lf", &C[i]);
for (i = 0; i < n; i++)
sigma[i] = C[i];
for (i = 0; i < m; i++)
CB[i] = C[num[i] - 1];//表示基变量的系数,因为将人类输入的转化为数组序号
}
//输出
void print()
{
int i, j;
printf("\n--------------------------------------------------------------------------\n");
for (i = 0; i < m; i++)
{
printf("%8.2f\tX(%d) %8.2f ", CB[i], num[i], b[i]);
for (j = 0; j < n; j++)
printf("%8.2f ", A[i][j]);
if (i != m - 1)
printf("\n");
}
printf("\n--------------------------------------------------------------------------\n");
printf("\t\t σ ");
for (i = 0; i < n; i++)
printf(" %8.2f", sigma[i]);
printf("\n--------------------------------------------------------------------------\n");
}
int main()
{
int i, j, k = 1;
int p, q;//q是换入变量
input();
printf("\n--------------------------------------------------------------------------\n");
printf("\tCB\tXB\tb\t");
for (i = 0; i < n; i++)
printf(" X(%d)\t", i + 1);
for (i = 0; i < n; i++)
x[i] = 0;
while (true)
{
q = FindSwapInVar();
if (q == -1)
{
print();
printf("\n恭喜您得到最优解!\n");
printf("最有解为:");
for (j = 0; j < m; j++)
x[num[j] - 1] = b[j];//因为x[n]对应第n+1个决策变量
for (i = 0; i < n; i++)
{
printf("x%d=%.2f ", i + 1, x[i]);
Z = Z + x[i] * C[i];
}
printf("\n最优值为:Z* = %.2f", Z);
break;
}
print();
p = FindSwapOutVar(q);
printf("\n进行第%d次迭代,迭代元位置为(%d,%d)\n", k++, p + 1, q + 1);
if (q == -1) break;
Iterate(p, q);
}
}
// momentum