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

%Written by: Ramprasad Yelchuru, PhD student, NTNU, Trondheim.
%Date : 13th 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)
    
    %Below cplex_miqp_full_final function is to solve the MIQP problem 
    [HF,lossF,HGyF,nmesF,ct]= cplex_miqp_full_final(Gy,Gdy,Juu,Jud,Wd,Wn,i,bigM);
    
    sset{i,1}=nonzeros(nmesF)'; % The i optimal measurement subsets  
    
    LF(i) = lossF; % The optimal loss with i measurement subsets
    
    ctF(i) = ct;  % The computational time taken to solve MIQP 
    
    nmesF' %To display the i optimal measurements
end
no = 2:i;

% To plot the optimal loss function with i optimal measurements
figure(10)
plot(no,LF(nu:i),'r-*','LineWidth',2,'MarkerSize',8);
xlabel('No. of measurements');
ylabel('Loss');
