数学规划模型（一）

一.概述

        我们先来了解一下什么是数学规划？

        数学规划是运筹学的一个分支，其用来研究：在给定的条件下（约束条件），如何按照某一衡量指标（目标函数）来寻求规划，管理工作中的最优方案。简单的说，也就是求目标函数在一定的约束条件下的极值问题。

        再来看数学规划的一般形式。

        用数学式子将其表示出来，x代表决策变量，一般会有多个自变量；f(x)代表目标函数；下面的g(x)代表的是约束条件，约束条件可分为不等式约束、等式约束、整数约束。在一个题目中，可能会有一种约束条件，也可能会有多种约束条件。

        对于不同的问题，可将数学规划分类。

        在这些类型中，线性规划是条件要求最苛刻的，需要目标函数和约束条件均是决策变量的线性表达式，但是求解起来会比较方便，如在MATLAB中使用单纯形法就可以求解。而非线性规划的条件相对比较宽松，只要目标函数或约束条件中有一个是决策变量的非线性表达式，但是因为目标函数可能会有多种非线性表达形式，约束条件也是，所以求非线性规划会出现很多种情况，因此求解非线性规划非常困难，在后面我们可以用蒙特卡洛模拟来求解非线性规划，也可以使用很多智能算法，如模拟退火、粒子群等。而整数规划可以简单的理解为约束条件中含有整数约束的数学规划问题，而0-1规划是整数规划的特例。

二.线性规划问题的求解

        我们在MATLAB中求解线性规划，首先需要将问题转化为MATLAB能识别的标准形式。

        MATLAB中线性规划的标准型，只能用来求解目标函数的最小值，但若是问题要求我们求解目标函数的最大值怎么办？我们可以将目标函数前加上负号，于是我们求解的负的目标函数的最小值即可转化为正的目标函数的最大值。

        目标函数可看成是两个向量的内积，其中向量c是目标函数的系数向量，x是决策变量的向量，向量的行数均取决于决策变量的个数。

        对于问题的约束条件要求我们转换成上述形式，同样的，若不等式约束是>=，我们则需要将不等式约束转换成<=的形式，方法是在不等式两边同时加上负号，这样不等式号就会改变；等式约束不必多说，上下界约束也可以当成不等式约束，而MATLAB为了求解起来方便而加入了上下界约束。注意，约束条件可能只对部分决策变量有约束。

        我们以第一题为例，向量c是目标函数的系数向量，x是决策变量的向量，向量的行数均取决于决策变量的个数；矩阵A的行数取决于不等式约束条件的个数，列数取决于决策变量的个数，是不等式约束条件的系数；向量b是不等式约束条件的上界值；1b向量为下届约束；而本题由于没有上界约束条件，因此ub向量的值均为正无穷。

        在MATLAB中我们可以很简单的求解线性规划问题。

        对应上面几道练习题，我们可以很容易的得出代码。

        我们来看具体的代码。

%% Matlab求解线性规划
% [x fval] = linprog(c, A, b, Aeq, beq, lb,ub, x0)  
% c是目标函数的系数向量，A是不等式约束Ax<=b的系数矩阵，b是不等式约束Ax<=b的常数项
% Aeq是等式约束Aeq x=beq的系数矩阵，beq是等式约束Aeq x=beq的常数项
% lb是X的下限，ub是X的上限，X是向量[x1,x2,...xn]' , 即决策变量。
% 迭代的初始值为x0（一般不用给）

%% 例题1
c = [-5 -4 -6]';  % 加单引号表示转置
% c = [-5 -4 -6];  % 写成行向量也是可以的，不过不推荐，我们按照标准型来写看起来比较正规
A = [1 -1 1;
        3 2 4;
        3 2 0];
b = [20 42 30]';   
lb = [0 0 0]'; 
[x fval] = linprog(c, A, b, [], [], lb)  % ub我们直接不写，则意味着没有上界的约束
% x =
%          0
%    15.0000
%     3.0000
% 
% fval =
%    -78

%% 例题2
c = [0.04 0.15 0.1 0.125]';  
A = [-0.03 -0.3 0 -0.15;
        0.14 0 0 0.07];
b = [-32 42]';
Aeq = [0.05 0 0.2 0.1];
beq = 24;
lb = [0 0 0 0]';
[x fval] = linprog(c, A, b, Aeq, beq, lb)
% x =
%          0
%   106.6667
%   120.0000
%          0
% 
% fval =
%     28

% 这个题可能有多个解，即有多个x可以使得目标函数的最小值为28（不同的Matlab版本可能得到的x的值不同，但最后的最小值一定是28）
% 例如我们更改一个限定条件：令x1要大于0（注意Matlab中线性规划的标准型要求的不等式约束的符号是小于等于0）
% x1 >0  等价于  -x1 < 0，那么给定 -x1 <= -0.1 (根据实际问题可以给一个略小于0的数-0.1)，这样能将小于号转换为小于等于号，满足Matlab的标准型
c = [0.04 0.15 0.1 0.125]';  
A = [-0.03 -0.3 0 -0.15;
        0.14 0 0 0.07
        -1 0 0 0];
b = [-32 42 -0.1]';
Aeq = [0.05 0 0.2 0.1];
beq = 24;
lb = [0 0 0 0]';
[x fval] = linprog(c, A, b, Aeq, beq, lb)
% x =
%     0.1000
%   106.6567
%   119.9750
%          0
%
% fval =
%    28.0000

%% 例题3
c = [-2 -3 5]';
A = [-2 5 -1;
          1 3 1];
b = [-10 12];
Aeq = ones(1,3);
beq = 7;
lb = zeros(3,1);
[x fval] = linprog(c, A, b, Aeq, beq, lb)
fval = -fval % 注意这个fval要取负号（原来是求最大值，我们添加负号变成了最小值问题）
% x =
%     6.4286
%     0.5714
%          0
% fval =
%   -14.5714
% fval =
%    14.5714

%% 多个解的情况
% 例如 ： min z = x1 + x2   s.t.  x1 + x2 >= 10
c = [1 1]';   
A = [-1 -1];
b = -10;
[x fval] = linprog(c, A, b)   % Aeq, beq, lb和ub我们都没写，意味着没有等式约束和上下界约束
% x有多个解时，Matlab会给我们返回其中的一个解

%% 不存在解的情况
% 例如 ： min z = x1 + x2   s.t.  x1 + x2 = 10 、 x1 + 2*x2 <= 8、 x1 >=0 ，x2 >=0 
c = [1 1]'; 
A = [1 2];
b = 8;
Aeq = [1 1];
beq = 10;
lb = [0 0]';
[x fval] = linprog(c, A, b, Aeq, beq, lb)  % Linprog stopped because no point satisfies the constraints.（没有任何一个点满足约束条件）

三.线性规划的典型例题

(1)例题1:生产决策问题

该题可参考解答列出目标函数和约束条件，直接给出代码供理解。

%% 生产决策问题
format long g   %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
% (1) 系数向量
c = zeros(9,1); % 初始化目标函数的系数向量全为0
c(1) = 1.25 -0.25 -300/6000*5;  % x1前面的系数是c1
c(2) = 1.25 -0.25 -321/10000*7;
c(3) = -250 / 4000 * 6;
c(4)  = -783/7000*4;
c(5) = -200/4000 * 7;
c(6) = -300/6000*10;
c(7) = -321 / 10000 * 9;
c(8) = 2-0.35-250/4000*8;
c(9) = 2.8-0.5-321/10000*12-783/7000*11;
c = -c;  % 我们求的是最大值，所以这里需要改变符号
% (2) 不等式约束
A = zeros(5,9);
A(1,1) = 5;  A(1,6) = 10;
A(2,2) = 7;  A(2,7) = 9; A(2,9) = 12;
A(3,3) = 6;  A(3,8) = 8;
A(4,4) = 4;  A(4,9) = 11;
A(5,5) = 7;  
b = [6000 10000 4000 7000 4000]';
% (3) 等式约束
Aeq = [1 1 -1 -1 -1 0 0 0 0;
            0 0 0 0 0 1 1 -1 0];
beq = [0 0]';
%（4）上下界
lb = zeros(9,1);

% 进行求解
[x fval] = linprog(c, A, b, Aeq, beq, lb)
fval = -fval
% fval =
%           1146.56650246305
%  注意，本题应该是一个整数规划的例子，我们在后面的整数规划部分再来重新求解。

        另外注意：本题实质上应该是个整数规划，我们的件数不应该出现小数，我们在后面部分再来求解。

(2)例题2:投料问题

        上面是我们针对该题作出的解设，目标函数即是料场j往工地i运送的水泥的吨数Xij乘上料场j与工地i的距离，再进行求和。我们重点理解约束条件，首先我们有不等式约束Xij>=0，很容易理解，因为运送水泥的吨数肯定是个正数；其次我们有第一个等式约束，两个料场送往每个工地的水泥的吨数之和应等于每个工地水泥的日需求量；第二个不等式约束是每个料场送往各个工地的水泥的吨数之和应小于料场的日储量20吨。

        接着我们要把该问题转换成MATLAB的标准型才能进行求解。

        我们的决策变量一共有12个，其中X1表示料场1送往工地1的水泥的吨数，X2则表示料场1送往工地2的水泥的吨数，以此类推，X7表示料场2送往工地1的水泥的吨数...目标函数的系数向量c是距离是已知的；将不等式约束和等式约束也转换成标准型。

        借助代码可以更好的理解本题的求解过程。

%% 投料问题
clear,clc
format long g   %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
% (1) 系数向量
a=[1.25  8.75  0.5  5.75  3  7.25];  % 工地的横坐标
b=[1.25  0.75  4.75	5  6.5  7.25];   % 工地的纵坐标
x = [5  2];  % 料场的横坐标
y = [1  7];  % 料场的纵坐标
c = [];  % 初始化用来保存工地和料场距离的向量 (这个向量就是我们的系数向量）
for  j =1:2
    for i = 1:6
        c = [c;  sqrt( (a(i)-x(j))^2 + (b(i)-y(j))^2)];  % 每循环一次就在c的末尾插入新的元素
    end
end
% (2) 不等式约束
A =zeros(2,12);
A(1,1:6) = 1;
A(2,7:12) = 1;
b = [20,20]';
% (3) 等式约束
Aeq = zeros(6,12);  
for i = 1:6
    Aeq(i,i) = 1;  Aeq(i,i+6) = 1;
end
% Aeq = [eye(6),eye(6)]  % 两个单位矩阵横着拼起来
beq = [3 5 4 7 6 11]';  % 每个工地的日需求量
%（4）上下界
lb = zeros(12,1);

% 进行求解
[x fval] = linprog(c, A, b, Aeq, beq, lb)
x = reshape(x,6,2)  % 将x变为6行2列便于观察（reshape函数是按照列的顺序进行转换的，也就是第一列读完，读第二列，即x1对应x_1,1，x2对应x_2,1）

% fval =
%           135.281541790676

        使用Lingo这个软件求解会更加方便，但其实用MATLAB也不是很复杂。

        下次我们将继续学习线性规划模型的其他分类情况。