单纯形法的C语言实现_c实现单纯形法

单纯形法原理

线性规划问题中，若存在最优解那么一定在某个基本可行解处取得，而我们可以将所有基本可行解找出来，代入目标函数后找到最优解，但是随着问题规模的增大，寻找所有基本可行解就变得异常困难。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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195