线性规划问题的目标函数灵敏度分析

某工厂生产A\B\C三种型号的产品，三种型号产品收益分别为2、1、3元，其生产受下列条件限制：人工<=1,材料<=3，每种产品的人工消耗分别为1/3、1/3、1/3；每种产品的用料消耗分别为1/3、4/3、7/3，试决定该生产的最优方案。

按照单纯形法，计算如下：

clc
clear all
close all

syms x1 x2 x3 y1 y2 y3 y4 y5;
 
% 根据例题构造线性规划（LP）的标准方程
C = [0 0 -2 -3 -1];
A = [1 0 1/3 1/3 1/3;0 1 1/3 4/3 7/3];
b = [1 ;3];
X = [y1 y2 x1 x2 x3].';

resultTable = lpSolve(C,X,A,b);

求解过程以及结果：

构造LP的单纯形表
[  0, 0, y1, y2,  x1,  x2,  x3]
[  0, 0,  0,  0,  -2,  -3,  -1]
[ y1, 1,  1,  0, 1/3, 1/3, 1/3]
[ y2, 3,  0,  1, 1/3, 4/3, 7/3]
 
更新单纯形表
[  0, 0, y1, y2, x1, x2, x3]
[  0, 6,  6,  0,  0, -1,  1]
[ x1, 3,  3,  0,  1,  1,  1]
[ y2, 2, -1,  1,  0,  1,  2]
 
更新单纯形表
[  0, 0, y1, y2, x1, x2, x3]
[  0, 8,  5,  1,  0,  0,  3]
[ x1, 1,  4, -1,  1,  0, -1]
[ x2, 2, -1,  1,  0,  1,  2]

通过求解结果，可以看出，x1=1,x2=3,x3=0时候，可以获得最大收益S=8。

当产品的收益在一定范围内波动时候，可能会影响最优方案的制定，问每种产品价格在什么范围内波动时候，该方案仍然是最优方案？这就是典型的目标函数灵敏度分析问题，求解答案如下：

基变量 x1 的价格区间 
[0.750000,3.000000]
基变量 x2 的价格区间 
[2.000000,8.000000]
非基变量 x3 的价格不超过 
4.000000

以上即为保证投产方案的稳健区间，只要产品收益在该范围内，则仍然为最优方案。

LP求解程序如下：

function resultTable = lpSolve(C,X,A,b)
% 目标函数系数C，变量X，约束条件系数A，约束项右侧常数b,约束条件个数m
% 传入参数来自于标准形式的LP方程

%% 求解过程
% 将标准方程中的系数矩阵A分解成基B与非基N,将X分解为基变量XB与非基变量XN
m = size(A,1);
B = A(1:m,1:m);
N = A(1:m,m+1:end);
CB = C(1,1:m);
CN = C(1,m+1:end);
XB = X(1:m,:);
XN = X(m+1:end,:);

invB = pinv(B);
% XB = invB*b - invB*N*XN

% 令XN = zeros（3，1）.'构造方程的一个基本解Xbasic
% invBxb >= 0时候，称B是Lp的可行基，此时Xbasic是Lp的基本可行解
invBxb = invB*b;
Xbasic = [invBxb ; zeros(3,1)];
%{
% 验证C*X 是否等于 CB*invB*b + (CN - CB*invB*N)*XN
C*X
CB*invB*b + (CN - CB*invB*N)*XN
%}

% 验证基本可行解是否是LP的最优解,当creterion>=0时候，Xbasic是最优解，B是最优基
creterion = C - CB*invB*A;

% 将LP转写称为典式形式y00\Ym\invBxN\invBxb,Y0n是LP的检验数（相对成本系数）
% minS= y00 + y0n*XN
% s.t : XB = invBxb - invBxN*XN
% X>=0
y00 = CB*invB*b;
Y0n = CN - CB*invB*N;
invBxN = invB*N;

% 构造方程的单纯形表
disp('构造LP的单纯形表')
TableLP = [0 0 X.';0 -y00 zeros(1,m) Y0n;XB invBxb eye(m) invBxN];
disp(TableLP)
while(min(Y0n) < 0)
    % 对TableLP中第二行非基部分进行判断，如果有负数元素，则认为当前解不是最优解
    % 然后判断invBxN的元素与0的关系，全<=0时候没有最优解，当存在>0的元素时候，执行
    % 换基操作，循环迭代直至最优解（为更具普适性，可以改用bland规则进行判断换基）
    vectorBuffer = TableLP(2,2+1:end);
    goleSensTable2 = TableLP;
    for i = 1:m
        mNonBasic = find(goleSensTable2(1,:) == TableLP(i+2,1));
        vectorBuffer(mNonBasic - 2) = 0;
    end
    qPosition = find(vectorBuffer < 0,1);
    
    % 取q所在列为主列，判断主列元素全部<=0是否成立，若成立则方程没有最优解，计算结束
    if(max(TableLP(3:end,2 + qPosition)) < 0)
        disp('LP方程没有最优解');
        break;
    else
        % 若不成立继续换基操作
        yi0 = TableLP(3:end,2);
        yiq = TableLP(3:end,2+qPosition);
        theta = min(yi0(yiq>0)./yiq(yiq>0));
        Jp = find(yi0./yiq == theta,1);
        
        % 取Jp所在行为主行，并以ypq为主元继续进行换基操作,xq为进基，xp为出基
        ypq = TableLP(2+Jp,2+qPosition);
        TableLP(2 + Jp,1) = TableLP(1,2 + qPosition);
        TableLP(2+Jp,2:end) = TableLP(2+Jp,2:end)/ypq;
        for i = 1:m+1
            if(i ~= Jp + 1)
                TableLP(i + 1,2:end) = TableLP(i + 1,2:end) - TableLP(Jp + 2,2:end)*...
                    TableLP(i + 1,2 +qPosition);
            end
        end
    end
    disp('更新单纯形表')
    disp(TableLP)
    
    % 更新检验数Y0n
    goleSensTable1 = TableLP;
    for i = 1:m
    mNonBasic = find(goleSensTable1(1,:) == TableLP(i+2,1));
    goleSensTable1(:,mNonBasic) = [];
    end 
    Y0n = goleSensTable1(2,3:end);
end

%% 打印结果
fprintf('最优方案：\n');
for i = 1:numel(X)-m
    if(ismember(XN(i),TableLP(3:end,1)))
        Jp = find(TableLP(:,1) == XN(i));
        fprintf('%s = %f \n',XN(i),TableLP(Jp,2));
    else
        fprintf('%s = 0\n',XN(i));
    end
end
maxS = TableLP(2,2);
fprintf('最优目标函数值：%f\n\n',maxS);


%% 目标函数系数的灵敏度分析
% 分离非基系数矩阵
goleSensTable = TableLP;
for i = 1:m
    mNonBasic = find(goleSensTable(1,:) == TableLP(i+2,1));
    goleSensTable(:,mNonBasic) = [];
end
goleSensTable = goleSensTable(:,3:end);
% 判断基变量的灵敏度区间
for i = 1:numel(X)-m
    if(ismember(XN(i),TableLP(3:end,1)))
        fprintf('基变量 %s 的价格区间 \n',XN(i));
        Jp = find(TableLP(:,1) == XN(i));
        vector1 = goleSensTable(2,:);
        vector2 = goleSensTable(Jp,:);
        Dmax = vector1(vector2>0)./vector2(vector2>0)*-1;
        Dmin = vector1(vector2<0)./vector2(vector2<0)*-1;
        if(isempty(max(Dmax)))
            xMin = -Inf;
        else
            xMin = -C(X == XN(i)) + max(Dmax);
        end
        
        if(isempty(min(Dmin)))
            xMax = Inf;
        else
            xMax = -C(X == XN(i)) + min(Dmin);
        end
        fprintf('[%f,%f]\n',xMin,xMax);
    end
end

% 判断非基变量灵敏度区间
for i = 1:numel(X)-m
    if(~ismember(XN(i),TableLP(:,1)))
        fprintf('非基变量 %s 的价格不超过 \n',XN(i));
        xNonMax = -C(i + m) + TableLP(2,2+m+i);
        fprintf('%f\n',xNonMax);
    end
end

%% 返回结果
resultTable = TableLP;
end

完整的程序输出如下：

构造LP的单纯形表
[  0, 0, y1, y2,  x1,  x2,  x3]
[  0, 0,  0,  0,  -2,  -3,  -1]
[ y1, 1,  1,  0, 1/3, 1/3, 1/3]
[ y2, 3,  0,  1, 1/3, 4/3, 7/3]
 
更新单纯形表
[  0, 0, y1, y2, x1, x2, x3]
[  0, 6,  6,  0,  0, -1,  1]
[ x1, 3,  3,  0,  1,  1,  1]
[ y2, 2, -1,  1,  0,  1,  2]
 
更新单纯形表
[  0, 0, y1, y2, x1, x2, x3]
[  0, 8,  5,  1,  0,  0,  3]
[ x1, 1,  4, -1,  1,  0, -1]
[ x2, 2, -1,  1,  0,  1,  2]
 
最优方案：
x1 = 1.000000 
x2 = 2.000000 
x3 = 0
最优目标函数值：8.000000

基变量 x1 的价格区间 
[0.750000,3.000000]
基变量 x2 的价格区间 
[2.000000,8.000000]
非基变量 x3 的价格不超过 
4.000000

以上即为目标函数灵敏度分析的实现方法。