% Main script for SIMC/I-SIMC comparison
% input process model number (Exx) for pmodel
%                       writen by: Chriss Grimholt
%

st = 200;  % simulation time
rt = 1;    % step change time
dt = 1;    % disturbance change time

res = [];

% getting process model

for pmodel = 1
   
    

[gp kp tp theta type] = model(pmodel);

% tuning factor

tc = theta;

% controller settings

[gc_s kc_s ti_s] = stuning(kp,tp,theta,tc,type);

[gc_i kc_i ti_i] = ituning(kp,tp,theta,tc,type);

ms_s = ms(gp, kc_s, ti_s);


ms_i = ms(gp, kc_i, ti_i);



% finding PO refrence

[kc_os, ti_os] = optimaltuning(gp, ms_s, type);
[kc_oi, ti_oi] = optimaltuning(gp, ms_i, type);

plot(ti_os(1:end-4),kc_os(1:end-4))

%finding refrence IEA
pf_os = [];
pf_oi = [];

for n = 1:2
for i = 1:size(ti_os,1);
    
    switch n
        case 1 
            gc = controller(1,kc_os(i),ti_os(i));
        case 2
            gc = controller(1,kc_oi(i),ti_oi(i));
    end

    [rv dv] = simtype(1);            % type = 1 gives setpoint change
   
    sim('system_model')
   
    pf_r = [iae(t,e), tv(u)];
   
    [rv dv] = simtype(2);            % type = 2 gives disturbance change
   
    sim('system_model')
   
    pf_d = [iae(t,e), tv(u)];
    n
    i/size(ti_os,1)
    switch n
        
        case 1
            pf_os = [pf_os; pf_r, pf_d];
        case 2
            pf_oi = [pf_oi; pf_r, pf_d];
    end
    
end
end

[iae_rs iae_rs_i] = min(pf_os(:,1));
[iae_ds iae_ds_i] = min(pf_os(:,3));

check1 = [ti_os(iae_rs_i),ti_os(iae_ds_i)];
check1a = [kc_os(iae_rs_i),kc_os(iae_ds_i)];


ref_s = [iae_rs iae_ds];

[iae_ri iae_ri_i] = min(pf_oi(:,1));
[iae_di iae_di_i] = min(pf_oi(:,3));

check2 = [ti_oi(iae_ri_i),ti_oi(iae_di_i)];

ref_i = [iae_ri iae_di];

Jos_temp = [pf_os(:,1)/iae_rs, pf_os(:,3)/iae_ds]*[0.5 0.5]';
[Jos jos_i]= min(Jos_temp);

Joi_temp = [pf_oi(:,1)./iae_ri, pf_oi(:,3)./iae_di]*[0.5 0.5]';
[Joi joi_i]= min(Joi_temp);

% Simulating responces;

for i = 1:2
    
    switch i
        
        case 1
            gc = gc_s;
        case 2
            gc = gc_i;
    end


    [rv dv] = simtype(1);            % type = 1 gives setpoint change
   
    sim('system_model')
   
    pf_r = [iae(t,e), tv(u)];
   
    [rv dv] = simtype(2);            % type = 2 gives disturbance change
   
    sim('system_model')
   
    pf_d = [iae(t,e), tv(u)];
    

    switch i
        
        case 1
            pf_s = [pf_r, pf_d];
        case 2
            pf_i = [pf_r, pf_d];
    end
    
end


% calculating the perfomrance

Js = 0.5*(pf_s(1)/iae_rs + pf_s(3)/iae_ds);
Ji = 0.5*(pf_i(1)/iae_ri + pf_i(3)/iae_di);


JsJ = Js/Jos;
JiJ = Ji/Joi;

% storing results

res_temp = [pmodel,tp/theta, type, kc_s, ti_s, ms_s, pf_s, JsJ, check1
            pmodel,tp/theta, type, kc_i, ti_i, ms_i, pf_i, JiJ, check2];
        
res=[res;res_temp];


end
res



