% 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 point
% We will use k to record iteration times.
x = [1; 1; 1; 1];
u = [0; 0];
s = [1; 1; 1; 1];
normd = max([norm(A*x-b)/(norm(b)+1), norm(A'*u+s-p)/(norm(p)+1), abs(x'*s/4)]);
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)
  % Compute D
  D = diag(x./s);
  invS = diag(1./s);

  % Compute mu, ra, rb, rc
  mu = x'*s/4/(k+1);
  ra = -A*x+b;
  rb = -A'*u-s+p;
  rc = -x.*s+mu;

  % Compute the direction vector dx, ds, du
  du = inv(A*D*A')*(ra + A*D*rb - A*invS*rc);
  ds = rb - A'*du;
  dx = invS*rc - D*ds;

  % Compute ax, as, au
  negx = (dx<0);
  ax = min([1, beta*min(-x(negx)./dx(negx))]);
  negs = (ds<0);
  as = min([1, beta*min(-s(negs)./ds(negs))]);
  au = as;

  % Update x, s, u
  x = x + ax*dx;
  s = s + as*ds;
  u = u + au*du;

  % Update normd
  normd = max([norm(A*x-b)/(norm(b)+1), norm(A'*u+s-p)/(norm(p)+1), abs(x'*s/4)]);

  % 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 = [1.000, 1.000, 1.000, 1.000]

After 1 iterations:
        x = [3.343, 1.114, 1.971, 0.200]

After 2 iterations:
        x = [3.849, 0.288, 3.438, 0.013]

After 3 iterations:
        x = [3.990, 0.003, 4.012, 0.018]

After 4 iterations:
        x = [3.995, 0.006, 3.992, 0.003]

After 5 iterations:
        x = [3.999, 0.001, 3.998, 0.001]

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

Published with MATLAB® 7.6