% setup the linear programming problem
%        min  f = p' * x
% subject to  A * x = b
%             x >= 0
p = [-5; -1; 0; 0];
A = [2,3,1,0; 2,1,0,1];
b = [12, 8];

% Set the identity matrix and a parameter beta
Id = eye(4);
beta = 0.99;

% Step 0, set the initial interior point
% We will use k to record iteration times.
% normd is the norm of direction vector d,
% which will be used in the stoping criteria.
x = [1/2; 2; 5; 5];
normd = 1;
k = 0;
fprintf('Starting interior point: \n');
fprintf('        x = [%5.3f, %5.3f, %5.3f, %5.3f]\n\n', x(1), x(2), x(3), x(4));

% path records the coordinates of points in the iteration.
% This information will be used to draw the graph later.
path = [x(1), x(2)];

% Now we start the iteration.
% For this problem, the stopting criteria is normd <= 1e-3.
while (normd >1e-3)
  % T defines the linear transformation x = T*y,
  % where y = [1;1;1;1]
  T = diag(x);
  y = [1;1;1;1];

  % Compute the transformed linear programming problem
  %        min  f = pk' * y
  % subject to  Ak * y = b
  %             y >= 0
  pk = T*p;
  Ak = A*T;

  % Compute the direction vector d, for the transformed problem.
  d = (Id - Ak'*inv(Ak*Ak')*Ak) * (-pk);
  normd = norm(d);

  % Compute alpha, which is the maximum value that ensures
  % y + alpha * d >= 0
  % Here negd contains all negative entries of vector d.
  negd = d(d<0);
  alpha = min(-1./negd);

  % Update y, and then x.
  ynew = y + beta*alpha*d;
  x = T*ynew;

  % Record the new point in path
  path = [path; x(1), x(2)];

  % update k and print result.
  k = k+1;
  fprintf('After %d iterations:\n', k);
  fprintf('        x = [%5.3f, %5.3f, %5.3f, %5.3f]\n\n', x(1), x(2), x(3), x(4));
end

% Draw the graph
figure(1);
clf;
hold on
% Draw the feasible region, which is a polygon.
plot([0,4,3,0,0], [0,0,2,4,0], '-', 'LineWidth', 2);
% Draw the path.
plot(path(:,1), path(:,2), '-*r');
hold off;
axis equal;
axis([0 4 0 4]);

Starting interior point: 
        x = [0.500, 2.000, 5.000, 5.000]

After 1 iterations:
        x = [1.534, 2.960, 0.050, 1.971]

After 2 iterations:
        x = [2.998, 1.985, 0.050, 0.020]

After 3 iterations:
        x = [3.020, 1.961, 0.079, 0.000]

After 4 iterations:
        x = [3.990, 0.020, 3.961, 0.000]

After 5 iterations:
        x = [4.000, 0.000, 4.000, 0.000]

After 6 iterations:
        x = [4.000, 0.000, 4.000, 0.000]

Published with MATLAB® 7.6