运筹教学|快速掌握单纯形法(附java代码)
什么是线性规划？
线性规划（Linear programming, 简称LP）是运筹学中研究较早、发展较快、应用广泛、方法较成熟的一个重要分支，它辅助人们进行科学管理、寻找线性约束条件下线性目标函数的极值。它广泛应用于军事作战、经济分析、经营管理和工程技术等领域，为合理地利用有限的人力、物力、财力等资源做出最优决策，提供科学依据。线性规划一般由决策变量、约束条件和目标函数三个部分组成。
线性规划的标准形式
由于线性规划的三个部分都有不同的形式，例如约束条件可以是等式也可以是不等式，目标函数可以是最大化或者最小化某个关于决策变量的线性函数的函数值。为了方便单纯形法的讲述，约定线性规划的标准形式为包含以下三个特征的的线性规划：
1. 目标函数统一求极大值(或者极小值)
2. 所有的约束条件都必须转化为等式，并且约束的右端项的值必须为非负
3. 所有变量的取值必须为非负值
下面给出一个线性规划的标准形式的例子
其中
代表决策变量，
代表
的价值系数，
代表第
种资源的数量，通常称为右端项，
代表技术系数或者工艺系数，即
每增加一个单位所消耗的资源
的数量。
如果找到一个解
满足所有的约束条件，则X为该线性规划问题的一个可行解，令c（X）为其目标函数值，如果在上述例子中存在一个可行解X*对于其余的所有可行解X都有
，则X*为该线性规划问题的最优解。
标准形式的转换
对线性规划的各个不满足标准形式的部分，都可以通过以下步骤进行转化：
1. 目标函数值的转换
统一求极大值，对于求极小值的目标函数
可以转换为
.
2. 变量的转换
(1) 如果变量满足非负约束，则不用做转换。
(2) 对于变量
，定义一个新的变量
替换变量x并且要求
。
(3) 对于变量
也即取值无约束，定义两个新的非负变量
通过令
对变量x进行替换，同样要求
。
(4) 在 (2) 和 (3) 中增加的新变量，其目标函数的系数与被替换的变量的系数相同。
3. 约束条件的右端项
对于右端项
，通过将约束两端的式子同乘 -1 (初等行变换) 可以完成转换。
4. 约束条件的转换
(1) 对于
的约束，在左端加入松弛变量
从而有
。
(2) 对于
的约束，在左端加入松弛变量
从而有
。
(3) 在 (1) 和 (2) 中定义的松弛变量的目标函数系数设为0，在实际问题中代表未被充分利用的资源或者是缺少的资源，因此松弛变量的目标函数的系数为0。
此外，给定线性模型的标准形式，为了构造出初始基变量，约束条件还可能需要加上人工变量。人工变量最终必须等于0才能保持原问题性质不变。为保证人工变量为0，在目标函数中令其系数为M。M为无限大的正数，这是一个惩罚项，倘若人工变量不为零，则目标函数就永远达不到最优，所以必须将人工变量逐步从基变量中替换出去。如若到最终表中人工变量仍没有置换出去，那么这个问题就没有可行解，当然亦无最优解。【Tips: 若原约束条件中已有线性无关的基向量，可以不需要再加入人工变量】
此处给出一个转换的例子供大家尝试
给定原问题：
则该问题的标准形式转换为：
为了方便运算的表示，我们用矩阵的形式来表示上述标准形式，即令
,
。同时，约束条件中的变量的系数矩阵表示为
，并且变量
的约束系数对应的列向量的集合表示为
。假设A的秩为m， 若A中存在一个m x m阶的非奇异子矩阵
，则称该子矩阵为基，用N表示非基部分的向量的组合。B中每一列系数称为一个基向量，每个基向量对应一个基变量。对于每一个基，如果我们令X中的非基变量取零，则原问题对应的方程组总存在一个唯一解(克莱姆法则)，我们称这个解为基解，但是这个基解不一定可行，有可能违反变量的非负约束，因此可行的基解称为基可行解。线性规划理论有以下定理：
定理1. 如果线性规划问题有最优解，则一定存在一个基可行解是它的最优解。
由于篇幅有限，这里不做证明。基于这个定理，如果我们把所有的基可行解找出来，一个一个试一下，如果存在最优解的话，我们就能够找出来。那么基解有多少个呢? 答案是
。当问题的规模上去以后，这样遍历显然耗费的时间非常长，而单纯形法，实际上就是基于一个基可行解进行搜索，并且迭代到搜索过程中得到的更优的基可行解，再不断重复这个搜索迭代，直到不存在更优的基可行解为止。
单纯形法步骤
1. 确定初始可行基和初始基可行解，建立初始单纯形表
2. 进行最优性检验，如果当前解为最优解，则算法停止，否则转入下一步
3. 如果存在非基变量其系数矩阵中对应的列向量中的系数均小于等于零，则此问题存在无界解，算法停止，否则转入下一步
4. 选取非负检验数中，数值最大的检验数(这里针对目标函数为求最大值的情况) 对应的非基变量作为入基变量，这里标记为
。
5. 根据
规则进行计算，令
为出基变量，则
应当满足
。 注意，这里的
和
指的是当前迭代中的单纯形表中的右端项和第i个约束中变量
的系数。
6. 通过高斯消元法，将
在A中的系数列
转换为单位列向量，如下图所示，然后转入第二步。
单纯形法举例
下面我们给出一个具体的例子来走一遍这个算法流程，具体的代码也会在文末给出
对于线性规划问题：
加入松弛变量，转化为标准形式得：
于是我们可以构造单纯形表，其中最后一行有星号的列为基变量。初始基可行解为(x_4, x_5, x_6, x_7)。
在单纯形表中，我们发现非基变量x的系数大于零，因此可以通过增加这些x的值，来使目标函数增加。
上表中c_2最大，因此我们选择x_2作为新的基变量。按照θ规则，x_7出基。通过高斯变换得到的新的单纯形表为：
继续计算，我们得到：
此时我们发现，所有非基变量的系数全部非正，即增大任何基变量的值并不能使得目标函数增大。于是我们可以断定该问题的最优解是z = 32, X = (0, 1, 3, 0, 2, 0, 0)。
恭~ 喜~ 大~ 家~
到这里单纯形法的原理就
搞！定！啦！
代码示例
代码分为两个文件，一个是读取和存储算例数据的代码，另一个则是根据算例数据进行计算的代码。
首先给出数据读取和存储部分的代码：
import java.io.BufferedReader;
import java.io.FileInputStream;
import java.io.InputStreamReader;
import java.util.ArrayList;
import java.util.Stack;
public class Data {
/**
* 系数矩阵的行数
*/
public int m;
/**
* 系数矩阵的列数
*/
public int n;
/**
* 系数矩阵
*/
public double[][] matrix;
/**
* 目标函数系数向量
*/
public int[] c;
/**
* 右端常数项向量
*/
public double[] b;
/**
* 基变量的索引
*/
public int[] xb;
/**
* 非基变量的索引
*/
public int[] xn;
/**
* 检验数
*/
public double[] z;
public double M = - 999999.9;
/**
* 读取数据
* @param path 数据文件路径
* @return 算例数据
*/
public static Data read(String path){
Data data = new Data();
try {
BufferedReader reader = new BufferedReader(new InputStreamReader(new FileInputStream(path)));
String str = reader.readLine();
String[] strs = str.split(" ");
data.m = Integer.parseInt(strs[0]);
data.n = Integer.parseInt(strs[1]);
data.matrix = new double[data.m][data.n];
data.c = new int[data.n];
data.b = new double[data.m];
data.xb = new int[data.m];
data.xn = new int[data.n - data.m];
data.z = new double[data.n];
strs = reader.readLine().split(" ");
for (int i = 0; i < data.n; i++) {
data.c[i] = Integer.parseInt(strs[i]);
}
for (int i = 0; i < data.m; i++) {
strs = reader.readLine().split(" ");
for (int j = 0; j < data.n; j++) {
data.matrix[i][j] = Integer.parseInt(strs[j]);
}
data.b[i] = Integer.parseInt(strs[strs.length - 1]);
}
}catch (Exception e){
e.printStackTrace();
}
return data;
}
}
单纯形法部分代码：
public class Simplex {
Data data;
/**
* 输入算例文件，进行求解并输出结果
* @param path 算例文件路径
*/
public void solve(String path){
data = Data.read(path);
print_info();
get_init();//获取初始解
while (true){
get_z();//计算检验数
if (check() != 0){//进行最优性检验
break;
}
if (update()){//对单纯形表进行更新
System.out.println("存在无界解");
break;
}
print_matrix();//输出当前单纯形表
print_sol();//输出当前解
}
}
private void print_info(){
System.out.println("初始信息");
StringBuilder s = new StringBuilder();
for (int i = 0; i < data.n; i++) {
s.append(data.c[i]).append("; ");
}
System.out.println("系数项 C: " + s);
for (int i = 0; i < data.m; i++) {
s.append(data.b[i]).append("; ");
}
System.out.println("右端项 b: " + s);
System.out.println("系数矩阵 A Matrix: ");
for (int i = 0; i < data.matrix.length; i++) {
s = new StringBuilder();
for (int j = 0; j < data.matrix[i].length; j++) {
s.append(data.matrix[i][j]).append(" ");
}
System.out.println(s);
}
s = new StringBuilder();
for (int i = 0; i < data.n; i++) {
s.append(data.z[i]).append("; ");
}
System.out.println("检验数 Z: " + s);
System.out.println("---------------------");
}
private void print_sol(){
System.out.println();
double val = 0;
for (int i = 0; i < data.m; i++) {
val += data.b[i] * data.c[data.xb[i]];
}
System.out.println("当前函数值: " + val + " ,当前解如下：");
StringBuilder s = new StringBuilder();
for (int i = 0; i < data.xb.length; i++) {
s.append("X_{").append(data.xb[i] + 1).append("}:").append(data.b[i]).append("; ");
}
s.append(", 其余变量取值为0");
System.out.println(s);
System.out.println("---------------------");
}
private void print_matrix(){
System.out.println("---------------------");
StringBuilder s ;
s = new StringBuilder();
for (int i = 0; i < data.m; i++) {
s.append(data.b[i]).append("; ");
}
System.out.println("当前右端项 b': " + s);
System.out.println("当前系数矩阵A: ");
for (int i = 0; i < data.matrix.length; i++) {
s = new StringBuilder();
for (int j = 0; j < data.matrix[i].length; j++) {
s.append(data.matrix[i][j]).append(" ");
}
System.out.println(s);
}
s = new StringBuilder();
for (int i = 0; i < data.n; i++) {
s.append(data.z[i]).append("; ");
}
System.out.println("当前检验数 Z: " + s);
s = new StringBuilder();
for (int i = 0; i < data.m; i++) {
s.append(data.xb[i] + 1).append("; ");
}
System.out.println("当前基变量X_{b}: " + s);
s = new StringBuilder();
for (int i = 0; i < data.n - data.m; i++) {
s.append(data.xn[i] + 1).append("; ");
}
System.out.println("当前非基变量X_{n}: " + s);
}
/**
* 获取初始解，初始基变量对应的系数均为单位向量
*/
private void get_init(){
int count = 0;
while (count < data.m){
out: for (int i = 0; i < data.n; i++) {
if (data.matrix[count][i] == 1){
for (int j = 0; j < data.m; j++){
if ( j == count){
continue ;
}
if (data.matrix[j][i] != 0){
continue out;
}
}
data.xb[count++] = i;
break ;
}
}
}
count = 0;
while (count < data.n - data.m){
out : for (int i = 0; i < data.n; i++) {
for (int j = 0; j < data.xb.length ; j++ ){
if (data.xb[j] == i){
continue out;
}
}
data.xn[count++] = i;
}
}
}
/**
* 计算检验数
*/
private void get_z(){
for (int i = 0; i < data.n; i++) {
double val = data.c[i];
for (int j = 0; j < data.m; j++) {
val -= data.c[data.xb[j]] * data.matrix[j][i];
}
data.z[i] = val;
}
}
/**
* 对检验数进行判断
* @return 0 : 存在可以进基的检验数; 1:有无穷多最优解; 2: 无解; 3：已经找到最优解
*/
private int check(){
boolean minus = false;//判断是否存在非基变量，其检验数为正数，如果存在，则赋值为 true
boolean zero = false;//判断是否存在非基变量，其检验数为零，如果存在，则赋值为 true
boolean nb = false;//判断基变量中是否存在非零人工变量，如果存在，则赋值为 true
for (int i = 0; i < data.xn.length; i++) {
if (data.z[data.xn[i]] > 0){
minus = true;
}
if (data.z[data.xn[i]] == 0){
zero = true;
}
}
for (int i = 0; i < data.xb.length; i++) {
//此处是判断基变量中是否有非零人工变量，此处通过目标函数系数进行判断，读者可以根据实际情况作调整
if (data.c[i] == data.M && data.b[i] > 0 ){
nb = true;
break;
}
}
if (zero){
System.out.println("有无穷多最优解");
return 1;
}else if (nb){
System.out.println("无解");
return 2;
}else if (!minus){
System.out.println("已经求得最优解");
return 3;
}else {
return 0;
}
}
/**
* 确定进基和出基变量
*/
private boolean update(){
//选择进基变量
int in = -1;
double max = -1;
for (int i = 0; i < data.xn.length; i++) {
if (data.z[data.xn[i]] > max){
in = i;
max = data.z[data.xn[i]];
}
}
//检查基变量的判别数
boolean um = true;
int out = -1;
double min = 9999999.9;
for (int i = 0; i < data.m; i++) {
double val = data.b[i] / data.matrix[i][data.xn[in]];
if (val >= 0){//用于无界解的判定
um = false;
}
if (val < min){
min = val;
out = i;
}
}
int val = data.xn[in];
data.xn[in] = data.xb[out];
data.xb[out] = val;
gaussian(data.xb[out], out);
return um;
}
/**
* 通过高斯消元法更新单纯形表
* @param col 检查的列
* @param row 正则化行
*/
private void gaussian(int col, int row){
double norm = data.matrix[row][col];
for (int i = 0; i < data.n; i++) {//对行进行处理
data.matrix[row][i] /= norm;
}
data.b[row] /= norm;
for (int count = 0; count < data.m; count++) {
if (count == row){
continue;
}
if (data.matrix[count][col] > 0){
double val = data.matrix[count][col];
data.b[count] -= (data.b[row] * val) ;
for (int j = 0; j < data.n; j++) {
data.matrix[count][j] -= (data.matrix[row][j] * val) ;
}
}
}
}
}
-The End-
文案：庄浩城(华中科技大学管理学院博士生一年级）
排版：程欣悦
指导老师：秦虎老师 （华中科技大学管理学院）
如对文中内容有疑问，欢迎交流。PS：部分资料来自网络。