%
% OptIdent_2_v7   Identifikace dynam.systemu OBEC.BUDICI FUNKCE
%				Minimalizace kvadrat.funkcionalu identifikacni odchylky
%                                                     JaJ, TUL, FM, MTI


clear all;	close all;	hold off;  clc;
global a1 a2 Tmax k a1v a2v Jv J
k=1;

disp(' ')
disp(' ')
disp('****************************************************************')
disp('    ITERACNI IDENTIFIKACE DYNAMICKEHO SYSTEMU')
disp('               s obecnym budicim signalem')
disp(' ')
disp('    Minimalizace kvadrat.funkcionalu identifikacni odchylky ')
disp('****************************************************************')

% Blok vstupnich parametru pro rizeni vypoctu
%**********************************************************************
Tmax=20;			% zadana delka simulace (a kriterialniho intervalu)
%dT=0.01;			% tvrde zadany max.simulacni krok
a2=0.1; a1=0.5;		% vychozi odhad parametru na zacatku optimalizace
%**********************************************************************

x0=[a2 a1];			% naplneni vektoru parametru na pocatku optimalizace
sim('Ident');	    % spousti z aktualniho adresare simulaci "Ident" 
                    % s pocatecnim nastavenim parametru
format short;
disp('   Vychozi odhad :')

tf(2,[a2 a1 1])     % vystup vychozi prenos.fce na obrazovku

disp('  a1, a2 ... koef.aprox.přenosu,J ... hodnota kriteria'), a1, a2, J
									
Ident		        % otevre okno se simulacnim schematem 

figure(1);			% otevre okno 1 a vykresli ys(t) s vychozim nastavenim 
plot(Time,ym,'g',Time,ys,'b','LineWidth',3); hold on;
title({'\bf Vyvoj odhadu reakce modelu'});
legend('\bf Namereny prubeh y_m','\bf Simulovany prubeh y_s');
plot([0 Tmax],[0 0],'k','LineWidth',1);
xlabel('\bf cas [s]'); ylabel('\bf y_m,  y_s');

figure(2);			% otevre okno 2 vizualizace parametru a vykresli osy
plot([-0.1 0.6],[0 0],'k',[0 0],[-0.1 0.6],'k','LineWidth',2); hold on; 
xlabel('\bf osa a_1'); ylabel('\bf osa a_2'); grid on; 
plot(a1,a2,'r x','LineWidth',3,'MarkerSize',12);
title({'\bf Vyvoj odhadu parametru a_1, a_2'});

a1v(k)=a1; a2v(k)=a2; Jv(k)=J;
disp('****************************************************************')
pause               % stiskni "Enter" pro pokracovani vypoctu

%**********************************************************************
disp('    Zacatek iterace')
disp(' ')
                    % nastav.parametru optimalizace(4.par.=max.pocet kroku)
OPTIONS=optimset('LargeScale','off','MaxIter',25,'Display','iter');	

x=fminsearch('KritIdent',x0,OPTIONS);  % optimalizace simplexovou metodou 
									   % (1.par.= minimalizovana ucel.funkce
                                       % je urcena pomoci"KritIdent"
                                       % (nutno volat z aktualniho adresare!)
%**********************************************************************
disp('     Konec iterace')
disp('****************************************************************')
pause               % stiskni "Enter" pro pokracovani vypoctu

a1opt=x(1);	a2opt=x(2);
a1=a1opt;	a2=a2opt;	sim('Ident');

disp('   Vysledek optimalizace :')
disp(' ')
disp('     a1opt., a2opt...koef.aprox.prenosu, J...opt.hodnota kriteria'),a1opt, a2opt, J
disp(' ')
tf(2,[a2 a1 1])     % vystup vysledne prenos.fce na obrazovku
disp(' ')
figure(1);			% otevre okno a vykresli prech.char. s konecnym nastav.
plot(Time,ys,'r','LineWidth',3); 

figure(2); plot(a1,a2,'r o','LineWidth',3,'MarkerSize',8);
disp('****************************************************************')
disp('     Konec prikladu iteracni identifikace')
disp('                                               JaJ, TUL, FM, MTI')
disp('****************************************************************')
disp(' ')