clear;

// tu można zamieścić podstawienia i obliczenia pomocnicze:

// które pozwolą w efekcie wyznaczyć parametry modeli transmitancyjnych:
//k1 = 
//T1 = 
//k2 = 
//T2 = 

xdel(winsid());
clc;

LAB_FS = 3;
PLOT_STEP = 100;

d = fscanfMat('dane.txt');
t = d(1:PLOT_STEP:$, 1)/1000;
u = d(1:PLOT_STEP:$, 4);
y = d(1:PLOT_STEP:$, 2);

// przebieg zarejestrowany
subplot(2, 1, 1);
plot(t, y);
title('Wyjście - zarejestrowany przebieg', 'fontsize', 3);
a = min(y), b = max(y), r = b - a;
ax = gca();
ax.data_bounds = [0, a - 0.05 * r; max(t), b + 0.05 * r];
ax.tight_limits = 'on';
h = get('hdl');
h.children(1).thickness = 2;
xlabel('t [s] - skala 1:300', 'fontsize', LAB_FS);
ylabel('glukoza [mg/dL]', 'fontsize', LAB_FS);
xgrid;
subplot(2, 1, 2);
plot(t, u);
title('Sterowanie - zarejestrowany przebieg', 'fontsize', 3);
a = min(u), b = max(u), r = b - a;
ax = gca();
ax.data_bounds = [0, a - 0.05 * r; max(t), b + 0.05 * r];
ax.tight_limits = 'on';
h = get('hdl');
h.children(1).thickness = 2;
h.children(1).foreground = 5;
xlabel('t [s] - skala 1:300', 'fontsize', LAB_FS);
ylabel('insulina [%]', 'fontsize', LAB_FS);
xgrid;

ue = u($);
ii = find(u == ue);
t = t(ii);
t = t - t(1);
y = y(ii);
y = -y + y(1);
t = t';
y = y';
du = u($) - u(1);

scf();
plot(t, y);
title('Wydzielony i odwrócony przebieg odpowiedzi skokowej', 'fontsize', 3);
h = get('hdl');
h.children(1).thickness = 2;
xlabel('t [s] - skala 1:300', 'fontsize', LAB_FS);
ylabel('glukoza [mg/dL]', 'fontsize', LAB_FS);
xgrid;

def = exists(['k1', 'T1', 'k2', 'T2'])
if sum(def) < 4 then
    error('Aby wykonać dalszą część skryptu, zdefiniuj zmienne: k1, T1, k2, T2');
end

// FUNKCJE DOPSAOWANIA
// inercja I-rzędu
function e = G1(p, z),
    y = z(1), t = z(2), k = p(1), T = p(2);
    ym = du*(k-k*%e^(-t/T));
    e = y - ym;
endfunction

// inercja II-rzędu, T1 = T2
function e = G2(p, z),
    y = z(1), t = z(2), k = p(1), T = p(2);
    ym = du*(-(k*t*%e^(-t/T))/T-k*%e^(-t/T)+k);
    e = y - ym;
endfunction

// inercja II-rzędu, T1 <> T2
function e = G3(p, z),
    y = z(1), t = z(2), k = p(1), T1 = p(2), T2 = p(3);
    ym = du*(-(T2*k*%e^(-t/T2))/(T2-T1)+(T1*k*%e^(-t/T1))/(T2-T1)+k);
    e = y - ym;
endfunction

// model 1
G1a =  syslin('c', k1 / (T1 * %s + 1));
y1a = du * csim('step', t, G1a);
er1a = sum((y - y1a).^2);
[p1, er1b] = datafit(G1, [y; t], [k1; T1]);
k1b = p1(1), T1b = p1(2);
G1b =  syslin('c', k1b / (T1b * %s + 1));
y1b = du * csim('step', t, G1b);
er1b = sum((y - y1b).^2);
scf();
plot(t, y, t, y1a, t, y1b);
h = get('hdl');
h.children(3).thickness = 2;
h = legend('eksperyment', ..
sprintf('model t10-t90 (k = %.2f, T = %.1f, er = %.1f)', k1, T1, er1a), ..
sprintf('model datafit (k = %.2f, T = %.1f, er = %.1f)', k1b, T1b, er1b), opt=4);
h.font_size=3;
xlabel('t [s] - skala 1:300', 'fontsize', LAB_FS);
ylabel('glukoza [mg/dL]', 'fontsize', LAB_FS);
title(['Dopasowanie do modelu inercyjnego I-rzędu:' ...
'$G(s)=\frac{k}{Ts+1}$'], 'fontsize', 3);
xgrid;

// model 2
G2a = syslin('c', k2 / (T2 * %s + 1)^2);
y2a = du * csim('step', t, G2a);
er2a = sum((y - y2a).^2);
[p2, er2b] = datafit(G2, [y; t], [k2; T2]);
k2b = p2(1), T2b = p2(2);
G2b =  syslin('c', k2b / (T2b * %s + 1)^2);
y2b = du * csim('step', t, G2b);
er2b = sum((y - y2b).^2);
scf();
plot(t, y, t, y2a, t, y2b);
h = get('hdl');
h.children(3).thickness = 2;
h = legend('eksperyment', ..
sprintf('model t10-t90 (k = %.2f, T = %.1f, er = %.1f)', k2, T2, er2a), ..
sprintf('model datafit (k = %.2f, T = %.1f, er = %.1f)', k2b, T2b, er2b), opt=4);
h.font_size=3;
xlabel('t [s] - skala 1:300', 'fontsize', LAB_FS);
ylabel('glukoza [mg/dL]', 'fontsize', LAB_FS);
title(['Dopasowanie do modelu inercyjnego II-rzędu' ...
'$G(s)=\frac{k}{(Ts+1)^2}$'], 'fontsize', 3);
xgrid;

// model 3
p = datafit(G3, [y; t], [k2; 0.9 * T2; 1.1 * T2]);
k3 = p(1), Ta3 = p(2), Tb3 = p(3);
G3 = syslin('c', k3, (Ta3 * %s + 1) * (Tb3 * %s + 1));
y3 = du * csim('step', t, G3);
er3 = sum((y - y3).^2);
scf();
plot(t, y , t, y3);
h = get('hdl');
h.children(2).thickness = 2;
h.children(1).foreground = 5;
h = legend('eksperyment', ..
sprintf('model datafit (k = %.2f, T1 = %.1f, T2 = %.1f, er = %.1f)', ..
k3, Ta3, Tb3, er3), opt=4);
h.font_size=3;
xlabel('t [s] - skala 1:300', 'fontsize', LAB_FS);
ylabel('glukoza [mg/dL]', 'fontsize', LAB_FS);
title(['Dopasowanie do modelu inercyjnego II-rzędu' ...
'$G(s)=\frac{k}{(T_1s+1)(T_2s+1)}$'], 'fontsize', 3);
xgrid;