This is the Affine Scaling Algorithm outlined by:

Linear and nonlinear programming with Maple: an interactive, applications ...

Paul E. Fishback,Paul F. Fishback

http://books.google.se/books?id=qwRLD9oXvQAC&pg=PA68&lpg=PA68&dq=%22Linear+and+Nonlinear+Programming+with+Maple%22+Section+2.6:+Interior+Point+Algorithm+I&source=bl&ots=gvpyCxR-eA&sig=zN00EKB-syghqS0zK2xuEa9rdXk&hl=sv&ei=nLGuToDuBOT04QSy1-n4Dg&sa=X&oi=book_result&ct=result&resnum=3&ved=0CCsQ6AEwAg#v=onepage&q=Fuel%20LP%20Pro&f=false

restart:

with(plots):

with(Optimization):

Ob := 4*x[1]+3*x[2]:

con1 := x[1] <= 8:

con2 := 2*x[1]+2*x[2] <= 28:

con3 := 3*x[1]+2*x[2] <= 32:

con4 := x[1] >= 0:

con5 := x[2] >= 0:

con := [con1, con2, con3, con4, con5]:

evalf(simplex[maximize](Ob, con), 3);

Maximize(Ob, con, iterationlimit = 1000);

constraintsA := {con1, con2, con3, con4, con5};

inequal(constraintsA, x[1] = 0 .. 15, x[2] = 0 .. 20, labels = ["x[1]", "x[2]"], labelfont = [TIMES, ROMAN, 14]);

** {x[1] = 4., x[2] = 10.}**

** [46., [x[1] = 4.00000000000000178, x[2] = 9.9999999999999982]]**

** constraintsA := {0 <= x[1], 0 <= x[2], x[1] <= 8, 2 x[1] + 2 x[2] <= 28, 3 x[1] + 2 x[2] <= 32}**

restart:

with(plots):

with(LinearAlgebra):

c := Vector[row]([4, 3, 0, 0, 0]):

X := Vector[row]([x[1], x[2], x[3], x[4], x[5]]):

c.X:

A := Matrix(3, 5, [1, 0, 1, 0, 0, 2, 2, 0, 1, 0, 3, 2, 0, 0, 1]):

b := `<,>`(8, 28, 32):

N := 15:

lambda := .75:

x := array(0 .. N):

x[0] := `<,>`(3, 4, 5, 14, 15):

for i from 0 to N-1 do

d := DiagonalMatrix(convert(x[i], list));

`~A` := A.d;

`~c` := c.d;

P := IdentityMatrix(5)-Transpose(`~A`).MatrixInverse(`~A`.Transpose(`~A`)).`~A`;

cp := P.Transpose(`~c`);

alpha := 1/abs(min(seq(cp[j], j = 1 .. 5)));

x[i+1] := d.(`<,>`(1, 1, 1, 1, 1)+alpha*lambda.cp)

end do:

interface(rtablesize = 20):

Matrix([convert([0, init, seq(it[i], i = 1 .. N-1)], Vector), Transpose(Matrix([convert([x1, x2, slack1, slack2, slack3, OB], Vector), seq(`<,>`(x[i], 4*x[i][1]+3*x[i][2]), i = 0 .. N-1)]))]);

constraintsA := {xx[1] >= 0, xx[2] >= 0, xx[1] <= 8, 2*xx[1]+2*xx[2] <= 28, 3*xx[1]+2*xx[2] <= 32}:

Ap := Array(1 .. 1, 1 .. 2):

Ap[1, 1] := display({inequal(constraintsA, xx[1] = 0 .. 15, xx[2] = 0 .. 20, labels = ["x[1]", "x[2]"], labelfont = [TIMES, ROMAN, 14]), pointplot([seq([x[i][1], x[i][2]], i = 0 .. N-1)], color = red, thickness = 3, connect = true), pointplot([seq([x[i][1], x[i][2]], i = 0 .. N-1)], color = red, thickness = 3, symbol = solidcircle, symbolsize = 14)}, title = "Constraint", titlefont = [TIMES, ROMAN, 16]):

Ap[1, 2] := display({pointplot([seq([i, 4*x[i][1]+3*x[i][2]], i = 0 .. N-1)], connect = true), pointplot([seq([i, 4*x[i][1]+3*x[i][2]], i = 0 .. N-1)], color = black, thickness = 3, symbol = solidcircle, symbolsize = 14)}, title = "Objective Function", titlefont = [TIMES, ROMAN, 16]):

display(Ap)