%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solution to Problem 3:
% Sensitivity Analysis:
% When a Little Means a Lot
%
% Generate three linear programming problems:
%       min c^T x
% subject to
%       Ax <= b
%        x >= 0
% where A, b, and c are defined below.
%
% Solve each problem, recording the Lagrange multipliers (dual variables).
% Also, for each problem, solve  nprob  problems formed from replacing
% A by A+E, where the elements of E are uniformly distributed in 
% [-tau,tau].
%
% exlinpro.m     Dianne O'Leary 05/2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize tau, nprob, and a figure-counter.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tau = .001;
jfig = 0;
nprob = 100;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loop through the three examples.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for iexample=1:3,

   if (iexample == 1)
      A = [1,1];
      b = 1;
      c = -[3; 1];
   elseif (iexample == 2)
      A = [1,1];
      b = 1;
      c = -[1.0005;0.9995];
   else
      A = [.01,5];
      b = .01;
      c = -[1;0];
   end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the problem as given.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   [xtrue,fval,exitflag,output,lambda] = linprog(c,A,b,[],[],[0;0],[Inf;Inf]);
   disp(sprintf('Lagrange multipliers for the three constraints of Problem %d',...
          iexample))
   disp(sprintf(' %f', lambda.ineqlin, lambda.lower))
   lamsav(iexample,:) = [lambda.ineqlin,lambda.lower'];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the perturbed problems.
% Save the solutions as rows of the array samp.
% Save the function values in f.
% Save the norm of the change in x in deltax.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   for i=1:nprob,
      Ap = A + tau*randn(1,2);
      Ap = A + 2*tau*(rand(1,2)-.5);
      [x,f(i)] = linprog(c,Ap,b,[],[],[0;0],[Inf;Inf]);
      samp(i,:) = x';
      deltax(i) = norm(x-xtrue);
   end
   xchange(iexample) = max(deltax);
   fchange(iexample) = max(abs(f-fval));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Plot all of the solutions for perturbations of this example.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   jfig = jfig + 1;
   figure(jfig)
   plot(samp(:,1),samp(:,2),'*')
   title(sprintf('LP Example %d: perturbed solutions',iexample))
   xlabel('x_1')
   ylabel('x_2')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Plot all of the function values for perturbations of this example.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   jfig = jfig + 1;
   figure(jfig)
   plot(f,zeros(1,nprob),'*')
   title(sprintf('LP Example %d: perturbed function values',iexample))
   xlabel('c^Tx')

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Print some summary information.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' ')
disp('Summary:')
disp(sprintf('Perturbation size:                %f',tau))
for iexample=1:3,
  disp(sprintf('Example %d Lagrange multipliers:',iexample)) 
  disp(sprintf('            for Ax <=b:           %f',lamsav(iexample,1)))
  disp(sprintf('            for lower bounds:     %f %f',lamsav(iexample,2:3)))
  disp(sprintf('            Maximal change in x:  %f',xchange(iexample)))
  disp(sprintf('            Maximal change in f:  %f',fchange(iexample)))
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Save the pictures.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
print -depsc lpx1.eps
figure(2)
print -depsc lpf1.eps
figure(3)
print -depsc lpx2.eps
figure(4)
print -depsc lpf2.eps
figure(5)
print -depsc lpx3.eps
figure(6)
print -depsc lpf3.eps