% ==================== PRIPRAVA DAT K IDENTIFIKACI ====================
% nacteni namerenych vzorku
load('meas_data','*');   

% vytvoreni matice s vzorky pro jednotlive teploty
Ny = horzcat(meas_data(:,1),meas_data(:,2),meas_data(:,3),...
     meas_data(:,4),meas_data(:,5));  %Ny

% vytvoreni matice s vzorky pro jednotlive klapky
Nu = horzcat(meas_data(:,6),meas_data(:,7),meas_data(:,8),...
     meas_data(:,9),meas_data(:,10)); %Nu
 
% vytvoreni objektu tridy iddata s vstupnimi i vystupnimi signaly
%  s nastavenim vlastnosti
alldata = iddata(Ny,Nu,40);
alldata.Tstart = 0;
alldata.TimeUnit = 'sec';
alldata.Notes = 'split to 12 experiments, one exp. is 4.6 day';
alldata.InputName = {'K1','K2','K3','K4','K5'};
alldata.OutputName = {'T1','T2','T3','T4','T5'};
alldata.InputUnit = {'%','%','%','%','%'};
alldata.OutputUnit = {'°C','°C','°C','°C','°C'};

% rozdeleni namerenych dat do 13 casovych useku o delce 4.5 dne
dexp1 = alldata(1:10000,:,:);
dexp2 = alldata(10001:20000,:,:);
dexp3 = alldata(20001:30000,:,:);
dexp4 = alldata(30001:40000,:,:);
dexp5 = alldata(60001:70000,:,:);
dexp6 = alldata(70001:80000,:,:);
dexp7 = alldata(90001:100000,:,:);
dexp8 = alldata(100001:110000,:,:);
dexp9 = alldata(110001:120000,:,:);
dexp10 = alldata(120001:130000,:,:);
dexp11 = alldata(130001:140000,:,:);
dexp12 = alldata(140001:150000,:,:);
dexp13 = alldata(150001:160000,:,:);

% spojeni vsech casti do jednoho iddata objektu s vice experimenty
dmerge1 = merge(dexp1,dexp2,dexp3,dexp4,dexp5,dexp6,dexp7,dexp8,...
                dexp9,dexp10,dexp11,dexp12,dexp13);
                               
% vytvoreni validacnich dat
    dvalid1 = alldata(100001:125000,:,:);
    dvalid1.Tstart = 0;
    dvalid1.Notes = 'validation data';

% odstraneni stredni hodnoty                
ddetr1 = detrend(dmerge1,0);
dvdetr1 = detrend(dvalid1,0);

% prevzorkovani na periodu Ts = 120s
Q = 3;
P = 1;
dres = resample(ddetr1,P,Q);
dvres = resample(dvdetr1,P,Q);


% ==================== IDENTIFIKACE =====================================

% identifikace modelu pro teplotu T1
%  + prevody, validace, prechodove char.
dresT1 = dres(:,1,1:5);
dvresT1 = dvres(:,1,1:5);
n4s1T1 = n4sid(dresT1,1,'Focus','Sim','nk',[1 1 1 1 1],...
               'InitialState','backcast');
cn4s1T1 = d2c(n4s1T1);
czpkT1 = zpk(cn4s1T1,'m');
scrsz = get(0,'ScreenSize');
figure('Position',[10, 70, scrsz(3)-40, scrsz(4)-150]);
compare(dvresT1, n4s1T1,'r','InitialState','estimate');
figure('Position',[10, 70, scrsz(3)-40, scrsz(4)-150]);
step(czpkT1);

% identifikace modelu pro teplotu T2
%  + prevody, validace, prechodove char.
close all;
dresT2 = dres(:,2,1:5);
dvresT2 = dvres(:,2,1:5);
n4s1T2 = n4sid(dresT2,1,'Focus','Sim','nk',[1 1 1 1 1],...
               'InitialState','backcast');
cn4s1T2 = d2c(n4s1T2);
czpkT2 = zpk(cn4s1T2,'m');
figure('Position',[10, 70, scrsz(3)-40, scrsz(4)-150]);
compare(dvresT2, n4s1T2,'InitialState','estimate');
figure('Position',[10, 70, scrsz(3)-40, scrsz(4)-150]);
step(czpkT2);

% identifikace modelu pro teplotu T3
%  + prevody, validace, prechodove char.
close all;
dresT3 = dres(:,3,1:5);
dvresT3 = dvres(:,3,1:5);
n4s1T3 = n4sid(dresT3,1,'Focus','Sim','nk',[1 1 1 1 1],...
               'InitialState','backcast');
cn4s1T3 = d2c(n4s1T3);
czpkT3 = zpk(cn4s1T3,'m');
figure('Position',[10, 70, scrsz(3)-40, scrsz(4)-150]);
compare(dvresT3, n4s1T3,'InitialState','estimate');
figure('Position',[10, 70, scrsz(3)-40, scrsz(4)-150]);
step(czpkT3);

% identifikace modelu pro teplotu T4
%  + prevody, validace, prechodove char.
dresT4 = dres(:,4,1:5);
dvresT4 = dvres(:,4,1:5);
n4s1T4 = n4sid(dresT4,1,'Focus','Sim','nk',[1 1 1 1 1],...
               'InitialState','backcast','DisturbanceModel','None');
cn4s1T4 = d2c(n4s1T4);
czpkT4 = zpk(cn4s1T4,'m');
figure('Position',[10, 70, scrsz(3)-40, scrsz(4)-150]);
compare(dvresT4, n4s1T4,'InitialState','estimate');
figure('Position',[10, 70, scrsz(3)-40, scrsz(4)-150]);
step(czpkT4);

% identifikace modelu pro teplotu T5
%  + prevody, validace, prechodove char.
close all;
dresT5 = dres(:,5,1:5);
dvresT5 = dvres(:,5,1:5);
n4s1T5 = n4sid(dresT5,1,'Focus','Sim','nk',[1 1 1 1 1],...
               'InitialState','backcast','DisturbanceModel','None');
cn4s1T5 = d2c(n4s1T5);
czpkT5 = zpk(cn4s1T5,'m');
figure('Position',[10, 70, scrsz(3)-40, scrsz(4)-150]);
compare(dvresT5, n4s1T5,'InitialState','estimate');
figure('Position',[10, 70, scrsz(3)-40, scrsz(4)-150]);
step(czpkT5);
close all;

% vytvoreni matice prenosu
S = [czpkT1;czpkT2;czpkT3;czpkT4;czpkT5];
tFs = tf(S);

% dilci prenosy matice
s11 = tf(tFs.num{1,1},tFs.den{1,1});
s12 = tf(tFs.num{1,2},tFs.den{1,2});
s13 = tf(tFs.num{1,3},tFs.den{1,3});
s14 = tf(tFs.num{1,4},tFs.den{1,4});
s15 = tf(tFs.num{1,5},tFs.den{1,5});

s21 = tf(tFs.num{2,1},tFs.den{2,1});
s22 = tf(tFs.num{2,2},tFs.den{2,2});
s23 = tf(tFs.num{2,3},tFs.den{2,3});
s24 = tf(tFs.num{2,4},tFs.den{2,4});
s25 = tf(tFs.num{2,5},tFs.den{2,5});

s31 = tf(tFs.num{3,1},tFs.den{3,1});
s32 = tf(tFs.num{3,2},tFs.den{3,2});
s33 = tf(tFs.num{3,3},tFs.den{3,3});
s34 = tf(tFs.num{3,4},tFs.den{3,4});
s35 = tf(tFs.num{3,5},tFs.den{3,5});

s41 = tf(tFs.num{4,1},tFs.den{4,1});
s42 = tf(tFs.num{4,2},tFs.den{4,2});
s43 = tf(tFs.num{4,3},tFs.den{4,3});
s44 = tf(tFs.num{4,4},tFs.den{4,4});
s45 = tf(tFs.num{4,5},tFs.den{4,5});

s51 = tf(tFs.num{5,1},tFs.den{5,1});
s52 = tf(tFs.num{5,2},tFs.den{5,2});
s53 = tf(tFs.num{5,3},tFs.den{5,3});
s54 = tf(tFs.num{5,4},tFs.den{5,4});
s55 = tf(tFs.num{5,5},tFs.den{5,5});

% prenosova matice statickych zesileni (s=0)
S0 = [0.7    -0.126  0.243 0.281 -0.04;
      -0.013  0.344  0.037 0.24   0.203;
      -0.009  0.062  0.226 0.191  0.052;
      0.0169  0.029  0.078 0.182  0.063;
      0.03    0.023  0.052 0.067  0.079];

% cislo podminenosti pro s = 0
[U E V] = svd(S0);
cislo_podminenosti = max(max(E))/min(max(E));

% matice relativnich zesileni pro s = 0, RGA
St = (inv(S0))';
rgaS0 = S0.*St;


