求解标准形式线性规划问题的单纯形方法matlab程序

以下是一个标准形式线性规划问题的单纯形法MATLAB程序示例：

function [x, fval, exitflag] = simplex(A, b, c)
% SIMPLEX Solve standard LP problem using simplex method.
%   [X, FVAL, EXITFLAG] = SIMPLEX(A, B, C) solves the linear programming
%   problem:
%
%      min    c'x
%      s.t.   Ax = b
%             x >= 0
%
%   using the simplex method. A must be a matrix of size m-by-n, where m is
%   the number of constraints and n is the number of variables. B must be a
%   column vector of length m. C must be a column vector of length n. X is
%   the solution, FVAL is the optimal objective function value, and
%   EXITFLAG is 1 if the method converged, -1 if the problem is
%   infeasible, or -2 if the problem is unbounded.

[m, n] = size(A);
basis = n+1:m+n; % initial basis

while true
    % Phase 1: find basic feasible solution
    [B, N, bvec, cvec, x] = init(A, b, c, basis);
    if isempty(B)
        x = zeros(n, 1);
        fval = 0;
        exitflag = -1; % infeasible
        return
    elseif all(x >= 0)
        break
    end
    
    % Phase 2: optimize
    [B, N, bvec, cvec, x, exitflag] = optimize(A, b, c, B, N, bvec, cvec, x);
    if exitflag == -2
        fval = -Inf;
        return
    end
end

fval = c' * x;

end

function [B, N, bvec, cvec, x] = init(A, b, c, basis)
% INIT Find basic feasible solution.
%   [B, N, BVEC, CVEC, X] = INIT(A, B, C, BASIS) finds a basic feasible
%   solution to the linear programming problem:
%
%      min    c'x
%      s.t.   Ax = b
%             x >= 0
%
%   given an initial basis. A must be a matrix of size m-by-n, where m is
%   the number of constraints and n is the number of variables. B must be a
%   column vector of length m. C must be a column vector of length n. BASIS
%   is a vector of size m containing the initial basis indices. B is the
%   current basis indices, N is the nonbasis indices, BVEC is the current
%   basic variable values, CVEC is the current nonbasic variable values,
%   and X is the current solution.

[m, n] = size(A);
B = basis;
N = setdiff(1:n, B);
cb = c(B);
cn = c(N);

T = [A(:,B) eye(m)];
z = [cb' zeros(1, m)];
cvec = [cb' cn'];
bvec = b;

while true
    % compute reduced costs
    xn = T \ A(:,N);
    rn = cn' - cb' * xn;
    if all(rn >= 0)
        x = T \ b;
        x = [x(1:length(B)); zeros(length(N), 1)];
        break
    end
    
    % find entering variable
    [~, j] = min(rn);
    aj = A(:,N(j));
    if all(aj <= 0)
        x = [];
        break
    end
    
    % compute ratios
    xB = T \ b;
    xj = xB - (T \ aj) * (xB(N(j)) / (aj' * (T \ aj)));
    xj(xj < 0) = 0;
    if all(xj == 0)
        x = [];
        break
    end
    
    % find leaving variable
    [~, i] = min(xB ./ xj);
    B(i) = N(j);
    N(j) = basis(i);
    
    % update tableau
    ej = zeros(n, 1);
    ej(B(i)) = 1;
    T = T - (T(:,i) / T(i,i)) * ej';
    T(i,:) = T(i,:) / T(i,i);
    z = z - (z(i) / T(i,i)) * ej;
    cb(i) = c(B(i));
    cn(j) = c(N(j));
    cvec = [cb' cn'];
end

end

function [B, N, bvec, cvec, x, exitflag] = optimize(A, b, c, B, N, bvec, cvec, x)
% OPTIMIZE Optimize using simplex method.
%   [B, N, BVEC, CVEC, X, EXITFLAG] = OPTIMIZE(A, B, C, B, N, BVEC, CVEC, X)
%   optimizes the linear programming problem:
%
%      min    c'x
%      s.t.   Ax = b
%             x >= 0
%
%   given an initial basis and feasible solution. A must be a matrix of size
%   m-by-n, where m is the number of constraints and n is the number of
%   variables. B must be a vector of size m containing the current basis
%   indices. N must be a vector of size n-m containing the current nonbasis
%   indices. BVEC must be a vector of size m containing the current basic
%   variable values. CVEC must be a vector of size n-m containing the current
%   nonbasic variable values. X must be the current feasible solution.
%   EXITFLAG is 1 if the method converged, -2 if the problem is unbounded.

[m, n] = size(A);
cb = c(B);
cn = c(N);

while true
    % compute reduced costs
    xn = A(:,N) - A(:,B) * (bvec \ A(:,N));
    rn = cn' - cb' * xn;
    if all(rn >= 0)
        exitflag = 1;
        break
    end
    
    % find entering variable
    [~, j] = min(rn);
    aj = A(:,N(j));
    if all(aj <= 0)
        exitflag = -2;
        break
    end
    
    % compute ratios
    xB = bvec \ A(:,B);
    xj = xB - (bvec \ aj) * (xB(j) / (aj' * (bvec \ aj)));
    xj(xj < 0) = 0;
    if all(xj == 0)
        exitflag = -2;
        break
    end
    
    % find leaving variable
    [~, i] = min(x ./ xj);
    B(i) = N(j);
    N(j) = basis(i);
    
    % update tableau
    ej = zeros(n, 1);
    ej(B(i)) = 1;
    T = [A(:,B) bvec; aj' 0];
    T = T / T(i,i);
    T(:,i) = -T(:,i) / T(i,i);
    T(i,:) = T(i,:) / T(i,i);
    T(i,i) = 1 / T(i,i);
    bvec = T(1:m,end);
    cb(i) = c(B(i));
    cn(j) = c(N(j));
    cvec = [cb' cn'];
    x = T(1:m,end);
end

end

使用方法示例：

% Define problem data
A = [1 2 1; 2 1 1];
b = [4; 5];
c = [1; 1; 1];

% Solve LP problem
[x, fval, exitflag] = simplex(A, b, c);