Revised Simplex Method - Matlab Script

Question

I've been asked to write down a Matlab program in order to solve LPs using the Revised Simplex Method.

The code I wrote runs without problems with input data although I've realised it doesn't solve the problem properly, as it does not update the inverse of the basis B (the real core idea of the abovementioned method).

The problem is only related to a part of the code, the one in the bottom of the script aiming at:

Computing the new inverse basis B^-1 by performing elementary row operations on [B^-1 u] (pivot row index is l_out). The vector u is transformed into a unit vector with u(l_out) = 1 and u(i) = 0 for other i.

Here's the code I wrote:

    %% Implementation of the revised Simplex. Solves a linear
    % programming problem of the form
    %
    %   min c'*x
    %   s.t. Ax  = b
    %         x >= 0   
    %
    % The function input parameters are the following:
    %     A: The constraint matrix 
    %     b: The rhs vector 
    %     c: The vector of cost coefficients 
    %     C: The indices of the basic variables corresponding to an
    %        initial basic feasible solution
    % 
    % The function returns:
    %     x_opt: Decision variable values at the optimal solution  
    %     f_opt: Objective function value at the optimal solution
    %
    % Usage: [x_opt, f_opt] = S12345X(A,b,c,C)
    %  NOTE: Replace 12345X with your own student number 
    %        and rename the file accordingly

    function [x_opt, f_opt] = SXXXXX(A,b,c,C)
    %% Initialization phase 
    % Initialize the vector of decision variables
    x = zeros(length(c),1);

    % Create the initial Basis matrix, compute its inverse and    
    % compute the inital basic feasible solution
    B=A(:,C);
    invB = inv(B);
    x(C) = invB*b;


    %% Iteration phase
    n_max = 10;        % At most n_max iterations
    for n = 1:n_max    % Main loop

        % Compute the vector of reduced costs c_r 

        c_B = c(C);         % Basic variable costs
        p = (c_B'*invB)';   % Dual variables
        c_r = c' - p'*A;    % Vector of reduced costs

        % Check if the solution is optimal. If optimal, use 
        % 'return' to break from the function, e.g.

        J = find(c_r < 0); % Find indices with negative reduced costs

        if (isempty(J))
            f_opt = c'*x;
            x_opt = x;
            return;
        end

        % Choose the entering variable
        j_in = J(1);

        % Compute the vector u (i.e., the reverse of the basic directions) 

        u = invB*A(:,j_in);
        I = find(u > 0);

        if (isempty(I))
           f_opt = -inf;  % Optimal objective function cost = -inf
           x_opt = [];        % Produce empty vector []
           return         % Break from the function
        end

        % Compute the optimal step length theta

        theta = min(x(C(I))./u(I));

        L = find(x(C)./u == theta); % Find all indices with ratio theta

        % Select the exiting variable

        l_out = L(1);

        % Move to the adjacent solution 

        x(C) = x(C) - theta*u;
        % Value of the entering variable is theta
        x(j_in) = theta;


        % Update the set of basic indices C 

        C(l_out) = j_in;

        % Compute the new inverse basis B^-1 by performing elementary row 
        % operations on [B^-1 u] (pivot row index is l_out). The vector u is trans-
        % formed into a unit vector with u(l_out) = 1 and u(i) = 0 for
        % other i.

        M=horzcat(invB,u);

        [f g]=size(M);
        R(l_out,:)=M(l_out,:)/M(l_out,j_in); % Copy row l_out, normalizing M(l_out,j_in) to 1
        u(l_out)=1;
        for k = 1:f % For all matrix rows
           if (k ~= l_out) % Other then l_out
           u(k)=0;
           R(k,:)=M(k,:)-M(k,j_in)*R(l_out,:); % Set them equal to the original matrix Minus a multiple of normalized row l_out, making R(k,j_in)=0
           end
        end

        invM=horzcat(u,invB);

        % Check if too many iterations are performed (increase n_max to 
        % allow more iterations)
        if(n == n_max)
            fprintf('Max number of iterations performed!\n\n');
            return
        end
    end  % End for (the main iteration loop)
end  % End function

%% Example 3.5 from the book (A small test problem)
% Data in standard form:
% A = [1 2 2 1 0 0;
%     2 1 2 0 1 0;
%     2 2 1 0 0 1];
% b = [20 20 20]';
% c = [-10 -12 -12 0 0 0]';
% C = [4 5 6];           % Indices of the basic variables of 
%                        % the initial basic feasible solution
% 
% The optimal solution
% x_opt = [4 4 4 0 0 0]' % Optimal decision variable values
% f_opt = -136           % Optimal objective function cost

Solution

Ok, after a lot of hrs spent on the intensive use of printmat and disp to understand what was happening inside the code from a mathematical point of view I realized it was a problem with the index j_in and normalization in case of dividing by zero therefore I managed to solve the issue as follows. Now it runs perfectly. Cheers.

%% Implementation of the revised Simplex. Solves a linear
% programming problem of the form
%
%   min c'*x
%   s.t. Ax  = b
%         x >= 0
%
% The function input parameters are the following:
%     A: The constraint matrix 
%     b: The rhs vector 
%     c: The vector of cost coefficients 
%     C: The indices of the basic variables corresponding to an
%        initial basic feasible solution
% 
% The function returns:
%     x_opt: Decision variable values at the optimal solution  
%     f_opt: Objective function value at the optimal solution
%
% Usage: [x_opt, f_opt] = S12345X(A,b,c,C)
%  NOTE: Replace 12345X with your own student number 
%        and rename the file accordingly

function [x_opt, f_opt] = S472366(A,b,c,C)
    %% Initialization phase 
    % Initialize the vector of decision variables
    x = zeros(length(c),1);

    % Create the initial Basis matrix, compute its inverse and    
    % compute the inital basic feasible solution
    B=A(:,C);
    invB = inv(B);
    x(C) = invB*b;


    %% Iteration phase
    n_max = 10;        % At most n_max iterations
    for n = 1:n_max    % Main loop

        % Compute the vector of reduced costs c_r 

        c_B = c(C);         % Basic variable costs
        p = (c_B'*invB)';   % Dual variables
        c_r = c' - p'*A;    % Vector of reduced costs

        % Check if the solution is optimal. If optimal, use 
        % 'return' to break from the function, e.g.

        J = find(c_r < 0); % Find indices with negative reduced costs

        if (isempty(J))
            f_opt = c'*x;
            x_opt = x;
            return;
        end

        % Choose the entering variable
        j_in = J(1);

        % Compute the vector u (i.e., the reverse of the basic directions) 

        u = invB*A(:,j_in);
        I = find(u > 0);

        if (isempty(I))
           f_opt = -inf;  % Optimal objective function cost = -inf
           x_opt = [];        % Produce empty vector []
           return         % Break from the function
        end

        % Compute the optimal step length theta

        theta = min(x(C(I))./u(I));

        L = find(x(C)./u == theta); % Find all indices with ratio theta

        % Select the exiting variable

        l_out = L(1);

        % Move to the adjacent solution 

        x(C) = x(C) - theta*u;
        % Value of the entering variable is theta
        x(j_in) = theta;


        % Update the set of basic indices C 

        C(l_out) = j_in;

        % Compute the new inverse basis B^-1 by performing elementary row 
        % operations on [B^-1 u] (pivot row index is l_out). The vector u is trans-
        % formed into a unit vector with u(l_out) = 1 and u(i) = 0 for
        % other i.

        M=horzcat(u, invB);
        [f g]=size(M);
        if (theta~=0)
        M(l_out,:)=M(l_out,:)/M(l_out,1); % Copy row l_out, normalizing M(l_out,1) to 1
        end
        for k = 1:f % For all matrix rows
           if (k ~= l_out) % Other then l_out
           M(k,:)=M(k,:)-M(k,1)*M(l_out,:); % Set them equal to the original matrix Minus a multiple of normalized row l_out, making R(k,j_in)=0
        end
        end
        invB=M(1:3,2:end);


        % Check if too many iterations are performed (increase n_max to 
        % allow more iterations)
        if(n == n_max)
            fprintf('Max number of iterations performed!\n\n');
            return
        end
    end  % End for (the main iteration loop)
end  % End function

%% Example 3.5 from the book (A small test problem)
% Data in standard form:
% A = [1 2 2 1 0 0;
%     2 1 2 0 1 0;
%     2 2 1 0 0 1];
% b = [20 20 20]';
% c = [-10 -12 -12 0 0 0]';
% C = [4 5 6];           % Indices of the basic variables of 
%                        % the initial basic feasible solution
% 
% The optimal solution
% x_opt = [4 4 4 0 0 0]' % Optimal decision variable values
% f_opt = -136           % Optimal objective function cost