% lecture 10 generalized least squares
% 27 Sept 13 -- bwb
%            -- for lecture on 30 Sept 2013
%
% illustrates two types of optimization-related problems
% Problem 1 -- least squares fit of a line to data
% Problem 2 -- minimizes sum of squares of 3 inlet flows to
%              meet an outlet flow constraint
%           -- also, a weighted sum of squares is solved
%
% -------------------------------------------------------------
% parameter fitting as a least squares problem
%
  a = 1;
  b = 1;
%
  t = 0:2:10;   % independent variable (creates row vector)
  t = t';       % creates column vector
%  add some noise
  rng('default')   % same random sequence each time the script is run
  noise = 0.5*randn(6,1);
  z = a + b.*t + noise;   % measured dependent variable
%
  figure(1)               % plot of data
  plot(t,z,'ko')          % black circles for data points
  xlabel('t')             % xaxis label
  ylabel('z')             % yaxis label
  title('Find best fit of z(i) = alpha + beta*t(i)')
%
  A = [ones(6,1) t];
  b = z;
  x = inv(A'*A)*A'*b          % least squares solution for alpha & beta
%
  alpha = x(1);
  beta  = x(2);
%
  zhat = alpha + beta.*t;   % model dependent variable
%
  figure(2)                 % compare data and model on same plot
  plot(t,z,'ko',t,zhat,'k')
  legend('data','model')
  xlabel('t')
  ylabel('z')
  title('Best fit of z(i) = alpha + beta*t(i)')
%
%  alternative to inv in MATLAB
%
  xcheck = (A'*A)\A'*b
%
%  use pinv -- based on SVD
%
  xpinv = pinv(A)*b
%
%  note that all three results yielded the same alpha and beta
%
%------------------------------------------------
% flow example
% outlet flow, F4, equals sum of inlet flows, F1+F2+F3
% 
% minimize the sum F1^2 + F2^2 + F3^2
%  s.t. F1+F2+F3 = F4
% that is, min x'x, s.t. Ax=b
%
  Aflow = [1 1 1]
  bflow = 10
  Aflow'*Aflow
  rank(Aflow'*Aflow)
%  Flowvec = (Aflow'*Aflow)\Aflow'*bflow  (*** not invertible! ***)
%
%  the alternative x = A'(AA')^-1*b is better!
%
  Flowv   = Aflow'*inv(Aflow*Aflow')*bflow
%
%  now, use the SVD based method (pinv)
%
  Flowvec = pinv(Aflow)*bflow   % pseudo-inverse uses SVD
  sum(Flowvec)
%
%  illustrate SVD to calculate the generalized inverse
%
  [U,S,V] = svd(Aflow)
%  using knowledge of dimensions for the next line
  FlowvecSVD = V(:,1)*(1/S(1,1))*U'*bflow
%
% ---------------------------------------
% weighted flow example
% valves sized such that valve 1 can handle 4 times as much flow as valve 2
%   which can handle 4 times as much flow as valve 3
% The following weight matrix is then used
  Wflow = [1 0 0; 0 4 0; 0 0 16]
% 
% FlowvecW = (Aflow'*Wflow*Aflow)\Aflow'*Wflow*bflow
  Winv = inv(Wflow);
  FlowvecW = Winv*Aflow'*inv(Aflow*Winv*Aflow')*bflow
  sum(FlowvecW)
