% Solves: Maximize c_1^Tx_1 + c_2^Tx_2 + d^Tw
%         s.t.     A_1x_1 + A_2x_2 + Dw = b
%                  B_1x_1               = d_1
%                           B_2x_2      = d_2
%                  x_1 >= 0, x_2 >= 0, w >= 0
% using the Dantzig-Wolfe decomposition scheme
% The idea is to modify the revised simplex code "revised_simplex.m"
% on the course webpage to handle the block angular structure in the LP

% Basic idea:
% (1) Reformulate the block angular LP into another LP having fewer
% constraints but MANY MANY more variables using the FINITE BASIS theorem
% (Theorem 16.2) on page 244 of Chvatal.
% (2) Solve the reformulated LP using the revised simplex method with
%  column generation.


% Kartik's MATLAB code for MA/OR/IE 505
% March 30, 2006

% Key in problem parameters
% This is the LP (26.8) on page 435 of Chvatal.
% This is the block angular system that we considered in class.

c_1 = [3 4 0 0]';
c_2 = [4 7 0 0]';
d = [0 0]';
A_1 = [1 2 0 0; 1 1 0 0];
A_2 = [2 3 0 0; 1 1 0 0];
D = [1 0; 0 1];
B_1 = [2 3 1 0; 1 1 0 1];
B_2 = [1 3 1 0; 2 1 0 1];
b = [225 114]';
d_1 = [120 50]';
d_2 = [150 100]';

% Set the tolerance parameter
eps = 1e-6;


%%%%%%%%%%%%%%%%%%%
% The extreme points of
% P_1 = {x_1| B_1x_1 = d_1, x_1 >= 0}
% are (0 0 120 50), (0 40 0 0), (30 20 0 0), and (50 0 20 0) respectively
% The extreme points of
% P_2 = {x_2| B_2x_2 = d_2, x_2 >= 0}
% are (0 0 150 100), (50 0 100 0), (0 50 0 50), and (30 40 0 0)
% respectively.
% The feasible regions P_1 and P_2 are bounded, so we do not have any
% extreme rays
% Let us start with the 1st extreme point (origin in the space of original variables) of P_1 and
% the 1st extreme point (origin in the space of the original variables) of P_2 and the columns of D in the initial basis
% Use this information to construct A_B and c_B

A_B = [A_1*[0 0 120 50]' A_2*[0 0 150 100]' D; 1 0 0 0; 0 1 0 0];
c_B = [c_1'*[0 0 120 50]'; c_2'*[0 0 150 100]'; d];

btilde = [b; 1; 1];

% Construct the first basic feasible solution for the refomulated problem
x_B = A_B\btilde;

iter = 0;
x1 = [0 0 120 50]';
x2 = [0 0 150 100]';

while 1 == 1,
iter = iter + 1;

% Update lower bounds as discussed in class
% The value of a basic feasible solution in every iteration provides this
% bound
if iter == 1,
  lb(iter) = c_B'*x_B;
else
  lb(iter) = max(lb(iter-1),c_B'*x_B);
end  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 2 :- Solve A_B^Ty = c_B
% Set up the objective function for the two subproblems
% Subproblem i is
% theta_i = Min (c_i-A_i^Ty)^Tx_i 
%           s.t. B_ix_i = d_i, x_i >= 0
% Reduced costs are s_i = theta_i - z_i

y = A_B'\c_B;
z = y(3:4);
y = y(1:2);
s_1 = c_1-A_1'*y;
s_2 = c_2-A_2'*y;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the subproblems using the revised simplex code on the course
% webpage.
% We will use sensitivity analysis to solve the subproblems as soon as
% possible
% Use the previous primal optimal solution to the subproblem as the
% starting point

[obj1,x1,u1] = revised_simplex_nondeg(s_1,B_1,d_1,1e-6,x1);
[obj2,x2,u2] = revised_simplex_nondeg(s_2,B_2,d_2,1e-6,x2);
% We should be solving these two subproblems in parallel!

% Update the upper bounds using the subproblem information
if (iter == 1),
  ub(iter) = b'*y + d_1'*u1 + d_2'*u2;
else
  ub(iter) = min(ub(iter-1),b'*y + d_1'*u1 + d_2'*u2);
end  

sN = [obj1-z(1); obj2-z(2); d-D'*y];

[sNmax,k] = max(sN);

if sNmax <= eps,
   fprintf('We are done\n');
   fprintf('Number of iterations is %d\n',iter);
   fprintf('Optimal objective value is %f\n',c_B'*x_B);
   return;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 3 :- Solve A_Bd = a_{N(k)}
% Find theta = Min_{i=1,...,m|d_i > 0} xB(i)/d(i)
% Let theta = xB(l)/d(l)
% x_{B(l)} is the leaving basic variable
% Also check for unboundedness if d <= 0

if (k == 1),
  aN = [A_1*x1; 1; 0];
  temp = c_1'*x1;
elseif (k == 2),
  aN = [A_2*x2; 0; 1];
  temp = c_2'*x2;
elseif (k ==3),  
  aN = [D(:,1); 0; 0]';
  temp = d(1);
else
  aN = [D(:,2); 0; 0];
  temp = d(2);
end  

dcol = A_B\aN;
zz = find(dcol > eps)';

if (isempty(zz))
  error('System is unbounded\n');
end  

[theta,ii] = min(x_B(zz)./dcol(zz));
l= zz(ii(1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Update x_B, A_B, and c_B

A_B(:,l) = aN;
x_B = x_B - theta*dcol;
x_B(l) = theta;
c_B(l) = temp;

end; % while

% One can construct a solution for the original problem from the values of x_B
% but I have not done it in this demo.