% doing revised simplex with bounds "manually".
% problem is: min 3x_1 + 4x_2 - 5x_3 + 7x_4
% subject to x_1 - x_2 + 2x_3 = 6, -x_1 + 2x_2 + x_4 = 5,
% 1 <= x_1 <= 2,
% 2 <= x_2 <= 6,
% -1 <= x_3 <= 4,
% 1 <= x_4
%
p=[3 4 -5 7]'; A=[1 -1 2 0; -1 2 0 1]; b=[6 5]';
% choose initial nonbasic set to be {1,2},
% with x_1 at upper bound and x_2 at lower bound
N=[1 -2]; x(1)=2; x(2)=2;
% The nonbasic set is thus {3,4}, and we can calculate the
% corresponding values of x(3) and x(4) from the equality constraints
[L,U]=lu(A(:,B));
B=[3 4];
x(B)=U\(L\(b-A(:,abs(N))*x(abs(N))))
x =
2
2
3
3
% now compute reduced cost vector
[L,U]=lu(A(:,B));
u=L'\(U'\p(B)); c=p(abs(N))-A(:,abs(N))'*u
c =
12.5000
-12.5000
% BOTH nonbasic elements are eligible pivots. Since c_1>0, objective will
% decrease as x_1 is decreased below its upper bound, and since c_2<0, the
% objective will also decrease when x_2 is increased away from its lower bound.
% We set s=1 to try decreasing x_1.
s=1;
% calculate d, the negative of column 1 of the tableau.
d=U\(L\A(:,abs(N(s))))
d =
0.5000
-1.0000
% if x1=2-lambda, x3=3+.5*lambda and x4=3-lambda. the first of these
% three variables to hit a bound is actually x_1, which hits its lower
% bound when lambda=1. Let's take this step in x_1 and recalculate
% the basic variables, and make the relevant changes to N
N(1)=-1;
lambda=1; x(1)=2-lambda;
x(B)=U\(L\(b-A(:,abs(N))*x(abs(N))))
x =
1.0000
2.0000
3.5000
2.0000
% Now to the next step. Since B and abs(N) have not changed, c and u are the same.
c
c =
12.5000
-12.5000
% now the only eligible pivot column is s=2, since we can decrease the
% objective by allowing x_2 to increase from its lower bound, since
% c_N(2)<0.
s=2;
% calculate the negative of this tableau column
d=U\(L\(A(:,abs(N(s)))));
d
d =
-0.5000
2.0000
% when x_2=2+lambda, basic variables become x_3=3.5+.5*lambda and
% x4=2-2*lambda. The first variable to hit a bound is x_2, when
% lambda=.5. So set r=2, since B(2)=4;
r=2;
lambda=0.5;
x(2)=x(2)+lambda;
x(B)=x(B)-lambda*d;
x
x =
1.0000
2.5000
3.7500
1.0000
% update B and N
B(2)=2; N(2)=-4;
B
B =
3 2
N
N =
-1 -4
% start another iteration.
[L,U]=lu(A(:,B));
u=L'\(U'\p(B));
c=p(abs(N))-A(:,abs(N))'*u
c =
6.2500
6.2500
% Current x is optimal! Both nonbasic variables are at their lower
% bound, and in both cases the objective will increase if we increase
% them away from their lower bound. calculate optimal objective
z=p'*x
z =
1.2500
diary off