% aby wywolac te funkcje musimy podac PID i model np.:
% obiekt=[K,T1,T2,tau,czas symulacji] - w odpowiednie miejsca podac wartosc
% PID = [P,I,D]
% w command window wpisac np "out=PIDOptimFun([P,I,D],[K,T1,T2,tau,czas])"
function [FC, p, cr, u, I2] = PIDOptimFun(PID,obiekt)

% dane naszego regulatora
P=PID(1);
I=PID(2);
D=PID(3);

% dane naszego modelu
K=obiekt(1);
tau=obiekt(4);
T1=obiekt(2);
T2=obiekt(3);
czasSym=obiekt(5);



% symulacja naszego modelu w simulink
options = simset('SrcWorkspace', 'current');
simOut = sim('SI',czasSym ,options);

% wartosc marginesu delta wynosi od 1% do 5, my zakladamy 5, bo czemu nie:)
d=5;

% przeregulowanie zgodnie ze wzorem z prezentacji
p=abs(((max(simOut.y)-simOut.y(end))/simOut.y(end))*100);

% uchyb jest roznica miedzy wartoscia osiagnieta przez obiek po regulacji a
% sygnalem zadanym. Abs, bo mozliwe jest otrzymanie wartosci ujemnej, czego
% nie chcemy
u=abs(simOut.y(end) - simOut.u(end));


% czas regulacj wyniga z marginesu delta. Sam czas regulacji to czas po
% ktorego uplywie wartosc odpowiedzi ukladu nie przekracza odchylenia od
% wartosci zadanej (wynikajacej z delty)
for i=length(simOut.t):-1:1
    if simOut.y(i)>(1+d/100)*simOut.u(end) || simOut.y(i)<(1-d/100)*simOut.u(end)
       cr=simOut.t(i);
       break;
    end
end

% kryterium całkowe, potrzebne do funkcji celu - w tym przypadku calka
% kwadratu bledu
I2=trapz(simOut.t, simOut.e.^2);

% funkcja celu jak wynik funkcji
FC = 0.5*p+0.1*cr+0.2*u+0.2*I2;
end