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 