转载自:https://mp.weixin.qq.com/s/7xv9runzEGblvDi40ismBQ
线性规划(Linear programming, 简称LP)是运筹学中研究较早、发展较快、应用广泛、方法较成熟的一个重要分支,它辅助人们进行科学管理、寻找线性约束条件下线性目标函数的极值。它广泛应用于军事作战、经济分析、经营管理和工程技术等领域,为合理地利用有限的人力、物力、财力等资源做出最优决策,提供科学依据。
标准式VS矩阵式
标准式
线性规划问题必须有一个标准形式,主要包括以下三个特征:
1)目标函数统一为求极大值(或极小值);
2)所有约束条件(除变量的非负条件外)必须都是等式,约束条件右端常数项(right-hand-side)b_i必须全为非负值;
3)所有变量的取值必须全为非负值。
下面的模型即为线性规划问题的标准形式:
s.t.括号内的内容是约束条件和变量的非负约束条件,字母代表的含义分别是:
下面模型即为标准形式的展开型:
可行解与 最优解
若找到(x_1, x_2 ,..., x_n)的值满足所有约束条件,且每个变量的值非负,则(x_1, x_2 ,..., x_n)称为线性规划问题的可行解。使目标函数值达到最大值(或最小值)的可行解即为该问题最优解,求解线性规划问题的目标就是要找出目标函数的最优解。
如何将目标函数转化为标准型
因此,在用单纯形法求解前,需要将模型转化为标准形式。
这个过程包括四个部分的转换:
1. 目标函数的转换:
统一求极大值,若是求极小值,则可将下面的式子乘以(-1)。即:
转化为:
2. 变量的转换:
(1)对于已经是大于等于零的变量 x_j ≥ 0 不做变化;
(2)对于小于等于零的变量 x_j,取负号令其变为大于等于零的变量,即若 x_j ≤ 0,则 定义新变量x_j' = -x_j,x_j' ≥ 0;
(3)若 x_j 取值无约束,可令两个新的非负变量x_j', x_j'', 然后用x_j = x_j' - x_j''替换原问题中的x_j。
3. 约束条件的右端项常数的转换:
b_i < 0 时,只需将等式或不等式两端同乘(-1);
4. 约束条件的转换:
将所有不等式全部转换为等式:
对于“≤ ”型约束加入一个变量 x_s,x_s ≥ 0;
对于“≥ ”型约束则减去一个变量 x_s,x_s ≥ 0。
加到原约束条件中的变量,称为松弛变量,在实际问题中它表示未被充分利用的资源或缺少的资源,所以在引入模型后它们在目标函数中的系数均为零。
给定线性模型的标准形式,为了构造出初始基变量,约束条件还可能需要加上人工变量。人工变量最终必须等于0才能保持原问题性质不变。为保证人工变量为0,在目标函数中令其系数为M。M为无限大的正数,这是一个惩罚项,倘若人工变量不为零,则目标函数就永远达不到最优,所以必须将人工变量逐步从基变量中替换出去。如若到最终表中人工变量仍没有置换出去,那么这个问题就没有可行解,当然亦无最优解。
将上述的线性规划问题转化为标准形式,其结果应如下:
注:x_4, x_5是将自由变量x_3转化为非负变量而引入的新变量,x_6, x_7是松弛变量,x_8, x_9是人工变量。
矩阵式
因为表达方式简单,单纯形法的表示与定理的说明往往使用矩阵形式。上述标准形式的矩阵形式表示如下:
矩阵A如下式:
A为m×n矩阵。假设A的秩为m,即假设不存在冗余的约束条件,则m>n时,因为方程数量比变量数目多,必定有多个可行解,即可利用单纯形法来计算最优解。
单纯形法的算法步骤
1. 确定初始可行基和初始基可行解,
建立初始单纯形表;
2. 最优性检验 若在当前表的目标函数对应的行中,所有非基变量的系数非正,则可判断得到最优解,可停止计算。否则转入下一步;
3. 若单纯形表中1至m列构成单位矩阵,在j=m+1至n列中,若有某个对应x_k的系数列向量 P_k ≤ 0,则此问题是无界,停止计算。否则,转入下一步;
4. 挑选目标函数对应行中系数最大的非基变量作为进基变量。假设x_k为进基变量,按θ规则[1]计算,可确定x_l为出基变量,转下一步;
5. 以a_lk为主元素进行迭代(即用高斯消去法或称为旋转运算),把x_k所对应的列向量进行变换[2];
6. 重复2-5步,直到所有检验数非正后终止,得到最优解。
[1] θ规则
其中b_i是当前表中的右手项,a_ik即为在第i个约束中变量k的系数。
[2] x_k列变换
单纯形法举例
对于线性规划问题:
加入松弛变量,转化为标准形式得:
于是我们可以构造单纯形表,其中最后一行有星号的列为基变量。初始基可行解为(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).
/*
Author:Sun Jiaxuan
School:Huazhong University of Science and Technology
Compiler:TDM-GCC 4.7.1 64-bit Release
Time:2017-09-25 16:30
*/
#include <cstdio>
#include <climits>
#include <cstdlib>
#include <iostream>
#include <cstring>
#include <vector>
#include <fstream>
#include <set>
using namespace std;
vector<vector<double> > Matrix;//储存约束方程的系数矩阵
double Z;
set<int> P;//迭代器
size_t cn, bn;
bool Pivot(pair<size_t, size_t> &p)//返回0表示所有的非轴元素都小于0
{
int x = 0, y = 0;
double cmax = -INT_MAX;
vector<double> C = Matrix[0];
vector<double> B;
for( size_t i = 0 ; i < bn ; i++ )
{
B.push_back(Matrix[i][cn-1]);
}
for( size_t i = 0 ; i < C.size(); i++ )//在非轴元素中找最大的c
{
if( cmax < C[i] && P.find(i) == P.end())
{
cmax = C[i];
y = i;
}
}
if( cmax < 0 )
{
return 0;
}
double bmin = INT_MAX;//根据θ规则,求最小检验数,确定转出变量
for( size_t i = 1 ; i < bn ; i++ )
{
double tmp = B[i]/Matrix[i][y];
if( Matrix[i][y] != 0 && bmin > tmp )
{
bmin = tmp;
x = i;
}
}
p = make_pair(x, y);//存储转入转出变量
for( set<int>::iterator it = P.begin() ; it != P.end() ; it++)
{
if( Matrix[x][*it] != 0 )
{
P.erase(*it);
break;
}
}
P.insert(y);
return true;
}
void pnt()//输出原始系数矩阵
{
for( size_t i = 0 ; i < Matrix.size() ; i++ )
{
for( size_t j = 0 ; j < Matrix[0].size() ; j++ )
{
cout<<Matrix[i][j]<<"\t";
}
cout<<endl;
}
cout<<"result z:"<<-Matrix[0][cn-1]<<endl;
}
void Gaussian(pair<size_t, size_t> p)//高斯行变换
{
size_t x = p.first;
size_t y = p.second;
double norm = Matrix[x][y];
for( size_t i = 0 ; i < cn ; i++ )//主行归一化
{
Matrix[x][i] /= norm;
}
for( size_t i = 0 ; i < bn && i != x; i++ )
{
if( Matrix[i][y] != 0)
{
double tmpnorm = Matrix[i][y];
for( size_t j = 0 ; j < cn ; j++ )
{
Matrix[i][j] = Matrix[i][j] - tmpnorm * Matrix[x][j];
}
}
}
}
void solve()
{
pair<size_t, size_t> t;
while(1)
{
pnt();
if( Pivot(t) == 0 )//已得到最优解
{
return;
}
cout<<t.first<<" "<<t.second<<endl;
for( set<int>::iterator it = P.begin(); it != P.end() ; it++ )
{
cout<<*it<<" ";
}//输出迭代后的矩阵
cout<<endl;
Gaussian(t);
}
}
int main()
{
ifstream fin;
fin.open("sample.txt");//文件读入。在代码同一目录下创建sample的记事本文件,格式见后面批注
fin>>cn>>bn;
for( size_t i = 0 ; i < bn ; i++ )
{
vector<double> vectmp;
for( size_t j = 0 ; j < cn ; j++)
{
double tmp = 0;
fin>>tmp;
vectmp.push_back(tmp);
}
Matrix.push_back(vectmp);//初始矩阵存储
}
for( size_t i = 0 ; i < bn-1 ; i++ )
{
P.insert(cn-i-2);
}
solve();
}
/
//input:
///* Variables */
//var x1 >= 0;
//var x2 >= 0;
//var x3 >= 0;
///* Object function */
//maximize z: x1 + 14*x2 + 6*x3;
///* Constrains */
//s.t. con1: x1 + x2 + x3 <= 4;
//s.t. con2: x1 <= 2;
//s.t. con3: x3 <= 3;
//s.t. con4: 3*x2 + x3 <= 6;
//end;
/
//sample.txt文件示例:
//8 5
//1 14 6 0 0 0 0 0
//1 1 1 1 0 0 0 4
//1 0 0 0 1 0 0 2
//0 0 1 0 0 1 0 3
//0 3 1 0 0 0 1 6
/