%This is the main file to solve Self Optimizing Control Problems
% To find full H with structural constraints and
% as combinations of optimal measurement subsets

%Written by: Ramprasad Yelchuru, PhD student, NTNU, Trondheim.
%Date : 15th Nov 2009 

%The data required to execute these codes are 
% Process gains from u to y : Gy
% Process disturbance gains from d to y : Gdy
% Hessian of cost function w.r.t u : Juu
% Hessian of cost function w.r.t u and d : Jud
% The magnitude of disturbances : Wd
% The magnitude of noise on measurement y : Wn


load 'EvaporatorData.mat' %To load the data for Evaporator case study
% load 'OriginalData.mat' %To load the data for Distillation case study
[ny,nu]=size(Gy)

for i = 2:10 % i can be any value between nu to ny
    
    % M used in bigM formulation in MIQP
    % The user should choose appropriate bigM for each problem
    bigM = 200;  %bigM = 2 (distillation case study), bigM = 200 (evaporator case study)
    
    %Structual constraints to select measurements from different sections
    StructCon =[zeros(1,ny*nu) 0 1 0 1 1 0 1 0 0 1;...
        zeros(1,ny*nu) 1 0 1 0 0 1 0 1 1 0]; %To 
    Con = [floor(i/2);i-floor(i/2)];
    
    %Below cplex_miqp_full_structure function is to solve the MIQP problem
    [HFS,lossFS,HGyFS,nmesFS,ctS]= cplex_miqp_full_structure(Gy,Gdy,Juu,Jud,Wd,Wn,i,bigM,StructCon,Con);
    
    ssetFS{i,1}=nonzeros(nmesFS)';  % The i optimal measurement subsets with structural constraints 
    
    LFS(i) = lossFS; % The optimal loss with i measurement subsets with structural constraints 
    
    ctFS(i) = ctS;  % The computational time taken to solve MIQP with structural constraints 
    
    nmesFS'  %To display the i optimal measurements with structural constraints 
end
no = nu:i;

% To plot the optimal loss function with i optimal measurements with
% structural constraints
figure(11)
plot(no,LFS(nu:i),'k-s','LineWidth',2,'MarkerSize',8);
xlabel('No. of measurements');
ylabel('Loss');
