title: Revised Simplex Method - Matlab Script
body: <p>I've been asked to write down a Matlab program in order to solve LPs using the <strong>Revised Simplex Method</strong>. </p>
  
  <p>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).</p>
  
  <p>The problem is only related to a part of the code, the one in the bottom of the script aiming at:</p>
  
  <p><strong>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.</strong></p>
  
  <p>Here's the code I wrote:</p>
  
  <pre><code>    %% Implementation of the revised Simplex. Solves a linear
      % programming problem of the form
      %
      %   min c'*x
      %   s.t. Ax  = b
      %         x &gt;= 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 &lt; 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 &gt; 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
  </code></pre>