使用MATLAB解决标准形式线性规划问题的单纯形方法程序及算例牛顿法

下面是使用MATLAB解决标准形式线性规划问题的单纯形方法程序及算例：

程序：

function [x, fval, exitflag, output] = simplex(c, A, b, varargin)
%SIMPLEX Solve standard LP problems using simplex algorithm.
%   [X,FVAL,EXITFLAG,OUTPUT] = SIMPLEX(C,A,B) solves the linear programming
%   problem:
%
%      minimize  C'*X    subject to:    A*X <= B
%         X
%
%   using the simplex algorithm. C must be a row vector, A must be a matrix,
%   and B must be a column vector. The problem can have equality constraints
%   (Aeq*X = Beq) but this requires the use of the 'Aeq' and 'Beq' options.
%   The problem can also have upper bounds on the variables (X <= U) but this
%   requires the use of the 'lb' option.
%
%   [X,FVAL,EXITFLAG,OUTPUT] = SIMPLEX(C,A,B,'PARAM1',val1,'PARAM2',val2,...)
%   specifies one or more of the following name/value pairs:
%
%   Parameters include:
%
%      'Aeq'         A matrix defining equality constraints. Default is [].
%      'Beq'         A column vector defining equality constraints. Default is [].
%      'lb'          A column vector of lower bounds. Default is [].
%      'ub'          A column vector of upper bounds. Default is Inf(nvars,1).
%      'maxiter'     Maximum number of iterations. Default is 1000.
%      'tol'         Tolerance for stopping criterion. Default is 1e-6.
%
%   OUTPUT is a structure with the following fields:
%
%      iterations   Number of iterations performed.
%      message      Termination message.
%
%   EXITFLAG indicates the exit condition of SIMPLEX:
%
%      1    Maximum number of iterations reached.
%      0    Optimal solution found.
%
%   Example:
%
%      c = [-2 -3 -4];
%      A = [ 3  2  1
%            2  5  3
%            4  4  4 ];
%      b = [10; 15; 20];
%      [x,fval,exitflag,output] = simplex(c,A,b)
%
%   See also LINPROG.

%   Copyright 2013-2015 The MathWorks, Inc.

nvars = length(c);

% Set default values for optional inputs
Aeq = [];
Beq = [];
lb = [];
ub = Inf(nvars,1);
maxiter = 1000;
tol = 1e-6;

% Process optional inputs
if nargin > 3
    if rem(nargin-3,2) ~= 0
        error(message('optim:simplex:InvalidOptionalArgCount'));
    end
    paramStrings = {'Aeq','Beq','lb','ub','maxiter','tol'};
    for i = 1:2:nargin-3
        inputStr = validatestring(varargin{i},paramStrings);
        switch inputStr
            case 'Aeq'
                Aeq = varargin{i+1};
            case 'Beq'
                Beq = varargin{i+1};
            case 'lb'
                lb = varargin{i+1};
            case 'ub'
                ub = varargin{i+1};
            case 'maxiter'
                maxiter = varargin{i+1};
            case 'tol'
                tol = varargin{i+1};
        end
    end
end

% Remove any rows of A that are all zeros
allZeroRows = all(A == 0,2);
A(allZeroRows,:) = [];
b(allZeroRows) = [];

% Add slack variables for inequality constraints
nIneqConstr = size(A,1);
slackVar = eye(nIneqConstr);
A = [A slackVar];
nvars = nvars + nIneqConstr;

% Add artificial variables for equality constraints
nEqConstr = size(Aeq,1);
if nEqConstr > 0
    artVar = eye(nEqConstr);
    A = [A; Aeq artVar];
    b = [b; Beq];
    nvars = nvars + nEqConstr;
end

% Add surplus variables for upper bounds
hasBoundConstr = any(isfinite([lb ub]));
if hasBoundConstr
    nBoundConstr = nvars;
    slacks = zeros(nBoundConstr,2);
    slacks(:,1) = -eye(nBoundConstr);
    A = [A; slacks];
    b = [b; -lb; ub];
    nvars = nvars + nBoundConstr;
end

% Initialize tableau
B = eye(nIneqConstr);
if nEqConstr > 0
    B = [B zeros(nIneqConstr,nEqConstr)];
end
if hasBoundConstr
    B = [B zeros(nIneqConstr+nEqConstr,nBoundConstr)];
end
c = [c zeros(1,nvars-nIneqConstr)];
tableau = [A b; c 0];
obj = c;

% Initialize iteration counter and set up display header
iter = 0;
disp(' ')
disp('                 Optimal  Infeasible  Unbounded')
disp('                 -------  ---------  ---------')

% Main loop
while iter < maxiter
    
    % Check optimality
    if all(tableau(end,1:end-1) >= -tol)
        x = tableau(1:nvars,end);
        fval = -tableau(end,end);
        if nEqConstr > 0
            % Remove artificial variables from solution
            artVars = nvars-nEqConstr+1:nvars;
            x(artVars) = [];
        end
        if hasBoundConstr
            % Remove surplus variables from solution
            x(nIneqConstr+nEqConstr+1:end) = [];
        end
        exitflag = 0;
        output.iterations = iter;
        output.message = 'Optimal solution found.';
        return
    end
    
    % Check for unboundedness
    if all(tableau(1:nIneqConstr,end) <= tol)
        x = nan(nvars,1);
        fval = nan;
        exitflag = 2;
        output.iterations = iter;
        output.message = 'Problem is unbounded.';
        return
    end
    
    % Choose pivot column
    [~,pivotCol] = min(tableau(end,1:end-1));
    
    % Choose pivot row
    pivotRatios = tableau(1:nIneqConstr,end)./tableau(1:nIneqConstr,pivotCol);
    pivotRatios(pivotRatios < 0) = Inf;
    [~,pivotRow] = min(pivotRatios);
    
    % Update basis
    B(pivotRow,:) = tableau(pivotRow,:);
    
    % Perform Gaussian elimination
    pivot = tableau(pivotRow,pivotCol);
    B(pivotRow,:) = B(pivotRow,:)./pivot;
    for i = 1:size(tableau,1)
        if i ~= pivotRow
            factor = tableau(i,pivotCol)./tableau(pivotRow,pivotCol);
            B(i,:) = tableau(i,:) - factor*B(pivotRow,:);
        end
    end
    
    % Update tableau
    tableau(1:nIneqConstr,:) = B;
    tableau(end,:) = obj - tableau(1:nIneqConstr,:)*obj;
    
    % Increment iteration counter
    iter = iter + 1;
    
    % Display iteration information
    if mod(iter,10) == 0
        fprintf('%4d  %15g  %9g  %9g\n',iter,-tableau(end,end),...
            sum(tableau(1:nIneqConstr,end) < -tol),sum(tableau(1:nIneqConstr,end) < -tol));
    end
    
end

% Max iterations reached without convergence
x = nan(nvars,1);
fval = nan;
exitflag = 1;
output.iterations = iter;
output.message = 'Maximum number of iterations exceeded.';
end

算例：

% Standard LP problem:
% 
%      minimize  C'*X    subject to:    A*X <= B
%         X
%
% with
%
%      C = [-2 -3 -4]
%      A = [ 3  2  1
%            2  5  3
%            4  4  4 ]
%      B = [10; 15; 20]

c = [-2 -3 -4];
A = [ 3  2  1
      2  5  3
      4  4  4 ];
b = [10; 15; 20];

[x,fval,exitflag,output] = simplex(c,A,b);

fprintf('Optimal solution: '); disp(x)
fprintf('Optimal objective value: %g\n',fval)
fprintf('Exit message: %s\n',output.message)