%% Math 464 Linear Optimization (Spring 2023) %% Lecture 19, 03/21/2023 %% Lexicographic pivoting rule to avoid cycling, and %% the big-M simplex method %% Lexicographic pivoting rule %% We start with the same example where simplex method with default %% pivoting rules cycles back to the starting tableau >> format rat >> format compact >> T = [ 0 -3/4 20 -1/2 6 0 0 0; 0 1/4 -8 -1 9 1 0 0; 0 1/2 -12 -1/2 3 0 1 0; 1 0 0 1 6 0 0 1] T = 0 -3/4 20 -1/2 6 0 0 0 0 1/4 -8 -1 9 1 0 0 0 1/2 -12 -1/2 3 0 1 0 1 0 0 1 6 0 0 1 >> [m,n] = size(T) m = 4 n = 8 %% We still pick j with the most negative c'_j to enter >> % Iteration 1 >> j=2; % x1 enters %% We now use the lexicographic rule for selecting the leaving variable >> i=2; T(i,:) = T(i,:)/T(i,j) T = 0 -3/4 20 -1/2 6 0 0 0 0 1 -32 -4 36 4 0 0 0 1/2 -12 -1/2 3 0 1 0 1 0 0 1 6 0 0 1 >> i=3; T(i,:) = T(i,:)/T(i,j) T = 0 -3/4 20 -1/2 6 0 0 0 0 1 -32 -4 36 4 0 0 0 1 -24 -1 6 0 2 0 1 0 0 1 6 0 0 1 %% Row 2 is lexicographically smaller than Row 3 now, as -32 < -24 ([0 1..] are same at start) >> l=2; %% We do the pivoting operations. Notice we do not have to scale the %% pivot row now, as all candidate rows are already scaled in the %% previous step. >> i=1; T(i,:) = T(i,:) - T(i,j)*T(l,:) T = 0 0 -4 -7/2 33 3 0 0 0 1 -32 -4 36 4 0 0 0 1 -24 -1 6 0 2 0 1 0 0 1 6 0 0 1 >> i=3; T(i,:) = T(i,:) - T(i,j)*T(l,:) T = 0 0 -4 -7/2 33 3 0 0 0 1 -32 -4 36 4 0 0 0 0 8 3 -30 -4 2 0 1 0 0 1 6 0 0 1 >> % Iteration 2 >> j=3; % x2 enters %% There is only one candidate for leaving variable, so no ties %% We can just do the pivoting as usual >> l=3; >> i=1; T(i,:) = T(i,:) - T(i,j)*T(l,:) T = 0 0 0 -2 18 1 1 0 0 1 -32 -4 36 4 0 0 0 0 1 3/8 -15/4 -1/2 1/4 0 1 0 0 1 6 0 0 1 >> i=2; T(i,:) = T(i,:) - T(i,j)*T(l,:) T = 0 0 0 -2 18 1 1 0 0 1 0 8 -84 -12 8 0 0 0 1 3/8 -15/4 -1/2 1/4 0 1 0 0 1 6 0 0 1 >> % Iteration 3 >> j=4; % x3 enters %% All three rows are candidates; but Row 4 already has a 1 in the %% pivot column, so we scale Rows 2 and 3. >> i=2; T(i,:) = T(i,:)/T(i,j) T = 0 0 0 -2 18 1 1 0 0 1/8 0 1 -21/2 -3/2 1 0 0 0 1 3/8 -15/4 -1/2 1/4 0 1 0 0 1 6 0 0 1 >> i=3; T(i,:) = T(i,:)/T(i,j) T = 0 0 0 -2 18 1 1 0 0 1/8 0 1 -21/2 -3/2 1 0 0 0 8/3 1 -10 -4/3 2/3 0 1 0 0 1 6 0 0 1 >> l=3; % [0 0 8/3... ] is the lex min row >> i=1; T(i,:) = T(i,:) - T(i,j)*T(l,:) T = 0 0 16/3 0 -2 -5/3 7/3 0 0 1/8 0 1 -21/2 -3/2 1 0 0 0 8/3 1 -10 -4/3 2/3 0 1 0 0 1 6 0 0 1 >> i=2; T(i,:) = T(i,:) - T(i,j)*T(l,:) T = 0 0 16/3 0 -2 -5/3 7/3 0 0 1/8 -8/3 0 -1/2 -1/6 1/3 0 0 0 8/3 1 -10 -4/3 2/3 0 1 0 0 1 6 0 0 1 >> i=4; T(i,:) = T(i,:) - T(i,j)*T(l,:) T = 0 0 16/3 0 -2 -5/3 7/3 0 0 1/8 -8/3 0 -1/2 -1/6 1/3 0 0 0 8/3 1 -10 -4/3 2/3 0 1 0 -8/3 0 16 4/3 -2/3 1 >> % Iteration 4 >> j = 5; % x4 enters %% Notice while z did not change in Iterations 1-3, it will change in %% this iteration, as the only candidate for leaving variable will be %% set as 1/16 > 0! >> l=4; >> i=4; T(i,:) = T(i,:)/T(i,j) T = 0 0 16/3 0 -2 -5/3 7/3 0 0 1/8 -8/3 0 -1/2 -1/6 1/3 0 0 0 8/3 1 -10 -4/3 2/3 0 1/16 0 -1/6 0 1 1/12 -1/24 1/16 >> i=1; T(i,:) = T(i,:) - T(i,j)*T(l,:) T = 1/8 0 5 0 0 -3/2 9/4 1/8 0 1/8 -8/3 0 -1/2 -1/6 1/3 0 0 0 8/3 1 -10 -4/3 2/3 0 1/16 0 -1/6 0 1 1/12 -1/24 1/16 >> i=2; T(i,:) = T(i,:) - T(i,j)*T(l,:) T = 1/8 0 5 0 0 -3/2 9/4 1/8 1/32 1/8 -11/4 0 0 -1/8 5/16 1/32 0 0 8/3 1 -10 -4/3 2/3 0 1/16 0 -1/6 0 1 1/12 -1/24 1/16 >> i=3; T(i,:) = T(i,:) - T(i,j)*T(l,:) T = 1/8 0 5 0 0 -3/2 9/4 1/8 1/32 1/8 -11/4 0 0 -1/8 5/16 1/32 5/8 0 1 1 0 -1/2 1/4 5/8 1/16 0 -1/6 0 1 1/12 -1/24 1/16 >> % Iteration 5 >> j=6; % x5 enters %% Only one candidate for leaving variable. >> l=4; >> i=4; T(i,:) = T(i,:)/T(i,j) T = 1/8 0 5 0 0 -3/2 9/4 1/8 1/32 1/8 -11/4 0 0 -1/8 5/16 1/32 5/8 0 1 1 0 -1/2 1/4 5/8 3/4 0 -2 0 12 1 -1/2 3/4 >> i=1; T(i,:) = T(i,:) - T(i,j)*T(l,:) T = 5/4 0 2 0 18 0 3/2 5/4 1/32 1/8 -11/4 0 0 -1/8 5/16 1/32 5/8 0 1 1 0 -1/2 1/4 5/8 3/4 0 -2 0 12 1 -1/2 3/4 >> i=2; T(i,:) = T(i,:) - T(i,j)*T(l,:) T = 5/4 0 2 0 18 0 3/2 5/4 1/8 1/8 -3 0 3/2 0 1/4 1/8 5/8 0 1 1 0 -1/2 1/4 5/8 3/4 0 -2 0 12 1 -1/2 3/4 >> i=3; T(i,:) = T(i,:) - T(i,j)*T(l,:) T = 5/4 0 2 0 18 0 3/2 5/4 1/8 1/8 -3 0 3/2 0 1/4 1/8 1 0 0 1 6 0 0 1 3/4 0 -2 0 12 1 -1/2 3/4 %% Notice we have an optimal solution, but have to scale Row-1 so as to %% get the x1 column basic in that row. >> i=2; T(i,:) = T(i,:)/T(i,2) T = 5/4 0 2 0 18 0 3/2 5/4 1 1 -24 0 12 0 2 1 1 0 0 1 6 0 0 1 3/4 0 -2 0 12 1 -1/2 3/4 >> % x1 = 1, x3= 1, x5 = 3/4 >> % Big-M method on favorite example >> open LPExample.m >> A = [1 1 -1 0 0; 3 1 0 -1 0; 3 2 0 0 1]; b = [2 4 10]'; c = [2 1 0 0 0]'; [m,n]=size(A); >> A = [A eye(m)] A = 1 1 -1 0 0 1 0 0 3 1 0 -1 0 0 1 0 3 2 0 0 1 0 0 1 %% Choosing M=10,000 should work well in this instance. >> M = 10000; >> c = [c; M*ones(m,1)] c = 2 1 0 0 0 10000 10000 10000 %% We choose the 3 artifical variables in the starting bfs >> Bind = [6 7 8]; >> Binv = eye(m); >> cB = c(Bind); >> cP = (c' - cB'*Binv*A)'; x = zeros(n,1); xB = Binv*b; x(Bind) = xB; >> cP cP = -69998 -39999 10000 10000 -10000 0 0 0 %% starting tableau >> T = [ [-cB'*xB cP']; [xB Binv*A]]; >> T T = -160000 -69998 -39999 10000 10000 -10000 0 0 0 2 1 1 -1 0 0 1 0 0 4 3 1 0 -1 0 0 1 0 10 3 2 0 0 1 0 0 1 %% We reinitialize m,n to size of T now! >> [m,n] = size(T) m = 4 n = 9 %% We use the default pivoting rule: most negative c'_j enters, and %% smallest l among candidate x_B(l)'s leaves; but we will not face %% any ties for the leaving variable in this problem >> % Iteration 1 >> j=2; % x1 enters >> inds = find (T(2:m,j) > 0); [theta,l] = min( T(inds+1,1)./T(inds+1,j) ); l=inds(l)+1; [theta,l] ans = 4/3 3 % l=3 ==> x7 leaves >> T(l,:) = T(l,:)/T(l,j); for i=setdiff(1:m,l); T(i,:) = T(i,:) - T(i,j)*T(l,:);end; T T = -200008/3 0 -49999/3 10000 -39998/3 -10000 0 69998/3 0 2/3 0 2/3 -1 1/3 0 1 -1/3 0 4/3 1 1/3 0 -1/3 0 0 1/3 0 6 0 1 0 1 1 0 -1 1 >> % Iteration 2 >> j=3; % x2 enters >> inds = find (T(2:m,j) > 0); [theta,l] = min( T(inds+1,1)./T(inds+1,j) ); l=inds(l)+1; [theta,l] ans = 1 2 % l=2 ==> x6 leaves >> T(l,:) = T(l,:)/T(l,j); for i=setdiff(1:m,l); T(i,:) = T(i,:) - T(i,j)*T(l,:);end; T T = -50003 0 0 -29999/2 -9999/2 -10000 49999/2 29999/2 0 1 0 1 -3/2 1/2 0 3/2 -1/2 0 1 1 0 1/2 -1/2 0 -1/2 1/2 0 5 0 0 3/2 1/2 1 -3/2 -1/2 1 %% *** At this point, we could pivot x5 into the basis and get to the %% optimal tableau! But that would be cheating in some sense, as %% we *KNOW* that the optimal basis has {x1,x2,x5}! Let's just %% continue with the default pivot rule instead. %% *** >> % Iteration 3 >> j=4; % x3 enters >> inds = find (T(2:m,j) > 0); [theta,l] = min( T(inds+1,1)./T(inds+1,j) ); l=inds(l)+1; [theta,l] ans = 2 3 % l=3 ==> x1 leaves >> T(l,:) = T(l,:)/T(l,j); for i=setdiff(1:m,l); T(i,:) = T(i,:) - T(i,j)*T(l,:);end; T T = -20004 29999 0 0 -19999 -10000 10000 29999 0 4 3 1 0 -1 0 0 1 0 2 2 0 1 -1 0 -1 1 0 2 -3 0 0 2 1 0 -2 1 >> % Iteration 4 >> j=5; % x4 enters >> inds = find (T(2:m,j) > 0); [theta,l] = min( T(inds+1,1)./T(inds+1,j) ); l=inds(l)+1; [theta,l] ans = 1 4 % l=4 ==> x8 leaves >> T(l,:) = T(l,:)/T(l,j); for i=setdiff(1:m,l); T(i,:) = T(i,:) - T(i,j)*T(l,:);end; T T = -5 1/2 0 0 0 -1/2 10000 10000 19999/2 5 3/2 1 0 0 1/2 0 0 1/2 3 1/2 0 1 0 1/2 -1 0 1/2 1 -3/2 0 0 1 1/2 0 -1 1/2 >> % Iteration 5 >> j=6; % x5 enters >> inds = find (T(2:m,j) > 0); [theta,l] = min( T(inds+1,1)./T(inds+1,j) ); l=inds(l)+1; [theta,l] ans = 2 4 % l=4 ==> x4 leaves >> T(l,:) = T(l,:)/T(l,j); for i=setdiff(1:m,l); T(i,:) = T(i,:) - T(i,j)*T(l,:);end; T T = -4 -1 0 0 1 0 10000 9999 10000 4 3 1 0 -1 0 0 1 0 2 2 0 1 -1 0 -1 1 0 2 -3 0 0 2 1 0 -2 1 >> % Iteration 6 >> j=2; % x1 enters >> inds = find (T(2:m,j) > 0); [theta,l] = min( T(inds+1,1)./T(inds+1,j) ); l=inds(l)+1; [theta,l] ans = 1 3 % l=3 ==> x3 leaves >> T(l,:) = T(l,:)/T(l,j); for i=setdiff(1:m,l); T(i,:) = T(i,:) - T(i,j)*T(l,:);end; T T = -3 0 0 1/2 1/2 0 19999/2 19999/2 10000 1 0 1 -3/2 1/2 0 3/2 -1/2 0 1 1 0 1/2 -1/2 0 -1/2 1/2 0 5 0 0 3/2 1/2 1 -3/2 -1/2 1 >> % This tableau is optimal! >> % z* = 3 for x1=1, x2=1, and x5=5. >> % Note that all three artificial variables x6,x7,x8 are zero here.