% y_nom=[109 62 77.67 73.16 58.44 3785 1800 1189]';
% 
% yd_opt= [109 61.98 77.68 73.16 58.44 3785 1800 1188  ;
%          109 85.86 62.00 62.00 58.41 3800 1799 1159  ;
%          109 62.00 85.21 63.00 58.45 3800 1799 1169  ;
%          109 62.00 62.00 85.86 58.41 3800 1799 1159]';
% 
% d=[ 0.1 0.1 0.1 0.1];
% 
% F = (yd_opt-repmat(y_nom,1,4))./repmat(d,length(y_nom),1);
% 
% H=null(F');
% H=H'


%================== SCALED CALCULATION ================================= 

% clc
% 
% y_max=[40 40 40 40 20 1000 1000 1000]';
% 
% y=y_nom./y_max;
% y_o=yd_opt./repmat(y_max,1,4);
% 
% d=[ 0.1 0.1 0.1 0.1];
% 
% 
% F = (y_o-repmat(y,1,4))./repmat(d,length(y),1);
% 
% H=null(F');
% H=H'


%================ NOMINAL VALUE SCALING ==================================

% clc
% 
% y_nom=[109 62 77.67 73.16 58.44 3785 1800 1189]';
% 
% y_opt= [109 61.98 77.68 73.16 58.44 3785 1800 1188  ;
%          109 85.86 62.00 62.00 58.41 3800 1799 1159  ;
%          109 62.00 85.21 63.00 58.45 3800 1799 1169  ;
%          109 62.00 62.00 85.86 58.41 3800 1799 1159]';
% 
% dy_max=[40 40 40 40 20 500 500 500]';
% 
% %y=y_nom./y_max;
% %y_o=yd_opt./repmat(y_max,1,4);
% 
% d=[ 0.01 0.01 0.01 0.01];
% 
% dy_opt = (y_opt-repmat(y_nom,1,4))./repmat(dy_max,1,4);
% 
% F = dy_opt./repmat(d,length(y_nom),1);
% 
% H=null(F');
% H=H'

%============ NO SCALING ==================================================

% y_nom=[98.585250352942225  78.504968998544442  81.693071578026647  78.504970176203173 54.855571356582168 3550 989.92 1600]';
% 
% y_opt= [98.585231890763637  78.505634474258841  81.691834146606467  78.505633519811752 54.864710377134038 3550 993.20 1600;
%        98.585371776347799  78.999972500134561  81.701208595312536  78.001228685264806 54.855571349067638 3550 989.91 1600;
%        98.585082958820337  78.510993735230727  81.681851925846431  78.510993696054214 54.859101600518350 3550 991.19 1600;
%        98.585372076153917  78.001190449564248  81.701231466368682  78.999986237008727 54.855571356582168 3550 989.92 1600]';
% 
%     d=[ 0.01 0.01 0.01 0.01];
%    
%    dy_opt=(y_opt-repmat(y_nom,1,4));
%    F = dy_opt./repmat(d,8,1);
%    H=null(F');
%    H=H'
%    
%    c1=H(1,:)
%    c2=H(2,:)
%    c3=H(3,:)
%    c4=H(4,:)
%    
%    null=c1*y_nom
%    
%===================== NEW CASE ===========================================
% Using initial values:
%x0=[ 100  50  70 50 1920    82     410    1550    440    660    750 340    510   1550    440    660    5770  1302  2240   -3.52  23.52 ];


% y_nom= [109.0000   81.9553   63.0727   81.9553  54.2464  3300    1055.5    1600]';
% 
% y_opt=[ 109.0000   81.9553   63.0727   81.9553  54.2455  3300    1055.2    1600 ;
%         109.0000   88.2364   69.0001   69.7469  54.1975  3300    1037.3    1600 ;
%         108.8560   81.4038   64.9367   81.4038  54.1731  3300    1028.1    1600 ;
%         109.0000   69.7469   69.0001   88.2364  54.1975  3300    1037.3    1600 ]';


% Nominal values for initial va colose to optima
 y_nom2=[109.0000  81.5473   82.3443   63.0918  54.2464  3300    1055.5    1600]';

% Using optimal initial values / near optimal 
y_nom= [109.0000   85.5045   77.5209   63.9579  54.2464  3300    1055.5    1600]';


% Optimal values for disturbence implementation
y_opt=[ 109.0000   85.5045   77.5208   63.9580  54.2455  3300    1055.2    1600 ;       % d1
        109.0000   88.2364   69.7468   69.0002  54.1975  3300    1037.3    1600 ;       % d2
        109.0000   69.7469   88.2364   69.0000  54.2163  3300    1044.3    1600 ;       % d3
        109.0000   69.0000   69.7469   88.2364  54.1975  3300    1037.3    1600 ]';     % d4


        d=[ 0.07 0.07 0.07 0.07];           % Disturbance                   
   
   dy_opt=(y_opt-repmat(y_nom,1,4));
   F = dy_opt./repmat(d,8,1)                % Sensitivity Matrix
   H=null(F')                              % Nullspace of F
   H=H';                                    % H-atrix
   
% Self-optimizing parameters from H-matrix

   c1=H(1,:)                            
   c2=H(2,:)
   c3=H(3,:)
   c4=H(4,:)
   
   null=c2*y_nom
   
%================= Ørten ==================================================

% y_nom= [  109.0000   84.3287   79.2181   63.4366  54.2464 3300    1055.5    1600]';
% 
% y_opt=[ 109.0000   84.3304   79.2158   63.4371  54.2455  3300    1055.2    1600 ;
%         109.0000   88.2364   69.7470   68.9999  54.1975  3300    1037.3    1600 ;
%         108.8560   81.4038   64.9367   81.4038  54.1731  3300    1028.1    1600 ;
%         108.9781   85.0042   78.0949   63.9999  54.1553  3300    1021.5    1600 ]';
% 
% 
% 
% d=[ 0.07 0.07 0.07 0.07];
%    
%    dy_opt=(y_opt-repmat(y_nom,1,4));
%    F = dy_opt./repmat(d,8,1);
%    H=null(F');
%    H=H'
%    
%    c1=H(1,:)
%    c2=H(2,:)
%    c3=H(3,:)
%    c4=H(4,:)
%    
%    null=c1*y_nom

%================ Case 3 ===================================================
%Using intitial values 
%x0=[ 109  82 63.1 82 200  4.8 24   1450  330   500  493   390     586  1250  330    1700  3300  1055.5  1600 34.2  54 ];


% y_nom= [  109.0000   81.5473   82.3443   63.0918  54.2464 3300    1055.5    1600]';
% 
% y_opt=[ 109.0000   81.5572   82.3348   63.0913  54.2455  3300    1055.2    1600 ;
%         109.0000   88.1393   70.8440   68.0000  54.1973  3300    1037.2    1600 ;
%         109.0000   69.7469   88.2364   69.0000  54.2163  3300    1044.3    1600 ;
%         109.0000   69.0000   69.7469   88.2364  54.1975  3300    1037.3    1600 ]';
% 
% 
% 
% d=[ 0.07 0.07 0.07 0.07];
% 
%    dy_opt=(y_opt-repmat(y_nom,1,4));
%    F = dy_opt./repmat(d,8,1);
%    H=null(F');
%    H=H'
%    
%    c1=H(1,:)                %c-atrix saved in optimizing_par3
%    c2=H(2,:)
%    c3=H(3,:)
%    c4=H(4,:)
%    
%    null=c1*y_nom