%
% Main Script for determining PO curve of PO tuning
%               writen by: Chriss Grimholt

st = 100;  % simulation time
rt = 1;    % step change time
dt = 1;    % disturbance change time
res = [];  % result matrix
ptau = []; % continious plotting matrix


% retrieving the process model

icase = 1;   % icase = 1 is Optimal tuning, icase = 2 is PO curve
first = 1;   % flagg for PO curve
teller = 0;

if icase == 1
    ms_temp = 1.59;
    simulation_runs = [0.01:0.01:0.99,1:0.1:15.9,16:1:50];
elseif icase == 2
    ms_temp = [1.59,1.4:0.01:2.5];
    simulation_runs = 1;
end

for ms_std = ms_temp
    
    
teller = teller +1;


for pmodel = simulation_runs;

[gp kc_s ti_s] = model(pmodel); % gets process model


stuning = [kc_s ti_s]; %SIMC tuning


pf = []; % temperary result storage

% simc tuning

gc = controller(1,kc_s,ti_s);   % PI-controller

    ms_s = ms(gp,kc_s,ti_s);    % SIMC ms

   [rv dv] = simtype(1);            % type = 1 gives setpoint change
   
   sim('system_model')
   
   pf_r_t = [iae(t,e) tv(u)];       % finds IAE and TV for setpoint
   
   [rv dv] = simtype(2);            % type = 2 gives disturbance change
   
   sim('system_model')
   
   pf_d_t = [iae(t,e) tv(u)];       % finds IAE and TV for disturbance
   
   pf_s = [pf_r_t, pf_d_t];
   
   
if icase  == 1
    ms_std = ms_s;
end

ti_sim=[0.01:0.005:3]';

[ti kc] = optimaltuning(gp, ms_std, pmodel, ti_sim);  % finds all kc values that 
                                              % gives ms=ms_std with given 
                                              % ti.

otuningc = [kc, ti];    % storing tuning

figure(1)

hold on
plot(ti(1:end-4),kc(1:end-4));  %plotting ti kc relation
hold off




for i = 1:size(otuningc,1);
    
   gc = controller(1,kc(i),ti(i));   % input: type, kc, ti
    
   [rv dv] = simtype(1);            % type = 1 gives setpoint change
   
   sim('system_model')
   
   pf_r_t = [iae(t,e) tv(u)];
   
   [rv dv] = simtype(2);            % type = 2 gives disturbance change
   
   sim('system_model')
   
   pf_d_t = [iae(t,e) tv(u)];
   
   pf = [pf; pf_r_t, pf_d_t];
   
end





switch icase
case 1
    [r_min r_i] = min(pf(:,1));
    [d_min d_i] = min(pf(:,3));
case 2
    
if ms_std == 1.59 && first == 1;
    [r_min r_i] = min(pf(:,1));
    [d_min d_i] = min(pf(:,3));
    first = 0;
else
    [r_minx r_i] = min(pf(:,1));
    [d_minx d_i] = min(pf(:,3));
end

    
end

rtuning = otuningc(r_i,:);
dtuning = otuningc(d_i,:);

% calculated the minimum value for the opbective function

pfw = [pf(:,1)/r_min, pf(:,3)/d_min];
pfw_s = [pf_s(:,1)/r_min, pf_s(:,3)/d_min];

pfws = pfw*[0.5 0.5]';
pfws_s = pfw_s*[0.5 0.5]';

[pfws_min pfws_i] = min(pfws);

otuning = otuningc(pfws_i,:);


% storing restults

res_temp = [pmodel,gp.outputd,rtuning,r_min,dtuning,d_min,otuning,pfws_min,stuning,pfws_s,(pfws_s-pfws_min)/pfws_min,ms_std];

res = [res;res_temp];

% printing results

switch icase    
case 1
figure(2)
ptau = [ptau;pmodel];
plot(ptau,res(:,8),ptau,res(:,11))
case 2
figure(2)
plot(res(2:end,16),res(2:end,11),ms_s,pfws_s,'o')    
end    
    
end

end

disp(res)


