 

 
 %ISO_IDENT_ELEM
 clear
 close
 


%Elementary Global identification for simplicity tutorial sake !

%------

%Simulation of Process to be identified
% H(p)= kp. exp(-r.p)/(1+Tp.p) : transfer function
% r: pure time delay: here to be determined experimentaly 
%  but an external optimisation loop can be easily added 
% Kp ; gain et Tp : time constant are identified with the iso-distance domain

clear all
close all

%ELABORATION OF A TEST DATA FILE

% Process definition

kp=1;             % gain 
Tp= 50   ;        % time constant
tech =1  ;        % sampling period (arbitrary unit) 
ap=exp(-tech/Tp); bp=1-ap; % discrete model


% Generating a file supposedly coming from the real world, for simulation !
% variables
tf= 5000;        % observation horizon
w=1:1:tf;        % time
MV=ones(1,tf);   % manipulted variable : process known input
CV=ones(1,tf);   % process output
b=zeros(1,tf);   % noise on output 
bmv=zeros(1,tf); % noise on input
Tint=1300;       % test signal elaboration
sm=CV;           % model output
r= 10 ;          % pure delay expressed in samples
kb=2   ;         % additive noise on input
%-------------------------------------------
for ii=2+1+r:1:tf;  %ii=time
   b(ii)=b(ii-1)*.9+.1*(rand(1)-0.5)*kb;
   bmv(ii)=bmv(ii-1)+.8*.2*(rand(1)-0.5)*.1;
   % noise
   if ii< 801;MV(ii)=2+bmv(ii);CV(ii)=.5+b(ii); else;
      TT= Tint*0.99985^ii;
      MV(ii)= sign(sin(2*pi*ii/TT))+2;  %  steps of test-----------choice
      %MV(ii)=100;
   CV(ii)=CV(ii-1)*ap+bp*(MV(ii-1-r)*kp+b(ii)+.5); 
   end%  process output
end

%file 
save fichtest w CV MV
%------
%-------------------------------------------------
% Loading a file from real world ( here simulated)------------------------------------START
load fichtest
nfich='fichtest';
%change of name to make the procedure universal whatever the name of the physical measurements !
e=MV;sp=CV;sm=sp;
figure(1)%-----------------------------------------
plot(w,MV,'k',w,CV,'r',w,0,'b')
grid;zoom;
axis([1 tf 0 6])
text(200,1.1,'CV')
text(200,2.2,'MV')
title('ROUGH FILE ')
disp (' pause ')
pause

%break
%-----------------------------------------------FILTERING : essential
% Parallel filtering of course to avoid bias !

% HIGH PASS
Tph=500;% time constant
afh=exp(-tech/Tph);bfh=1-afh;
ufh=filter([1 -1],[1 -afh],e-e(1));
spfh=filter([1 -1],[1 -afh],sp-sp(1));

% LOW PASS 
Tpb=10;%time constant filter
afb=exp(-tech/Tpb);bfb=1-afb;
uf=filter([bfb],[1 -afb],ufh-ufh(1))+ufh(1);
spf=filter([bfb],[1 -afb],spfh-spfh(1))+spfh(1);
% restitution
e=uf;sp=spf;
figure(2)%-------------------------------------------
plot(w,e,'k',w,sp,'r')
grid;zoom;
axis([1 tf -3 3])
text(200,-0.7,'CV')
text(200,0.2,'MV')
title(' FILE ')

%break
%------------------------------------------------
%OBSERVATION of the physaical data
%selection of the identification horizon
tinit=input('time tinit(2500)= ');% time of start
tfin= input('time  tend(4950)= ');% time of end

%Computation of the norm of the CV and of the looked for level of ISO-DISTANCE
VPP=0;NIV =0;             % initialisation
niv=input('level of iso (0.05)=');  % ISO lEVEL
rm=input('delay (r=10)');    %pure time delay of model
for ii=2+r:1:tf
if ii>tinit & ii<tfin;VPP=VPP+CV(ii)*CV(ii);NIV =NIV +niv*niv;end
end

%------------------------------------------------
% Screening the parametric space for differnt times constant of the model
%N=input('N (80)'); %  number of time constants explored
N=70;
Cmin=10^10;         % max criterium for initialidsation
%the parametric space will be scanned from taumin to taumax
taumin = input('taumin(= 30)');
taumax = input('taumax(= 80)');

%Identification Global----------------------------------------------------------------ID

for j =1:1:N      
  taum(j)=taumin+(taumax-taumin)*(j-1)/(N-1);   % model time constant
  am=exp(-tech/taum(j));bm=1-am;                % discretisation model   
  sm(ii)=sp(ii);
  VMM=0;VMP=0;
  for ii=2+r:1:tf
     if ii<tinit;sm(ii)=sp(ii);else;
     sm(ii)=sm(ii-1)*am+bm*e(ii-1-rm);end       % THE GAIN IS =  ATTENTION ! 
     if ii>tinit & ii<tfin;
     VMM=VMM+sm(ii)*sm(ii); VMP=VMP+sm(ii)*sp(ii);end% projections
  end
  Vmm(j)=VMM;Vmp(j)=VMP;         % storing in mmeory
  K(j)=VMP/VMM;                  % optimal local gain for this time constant
  C=K(j)*K(j)*VMM-2*K(j)*VMP+VPP;% criterium  local for this time constant
  if C<Cmin; Cmin=C;kopt=K(j); tauopt=taum(j);Copt=C;jopt=j;end; % selection of the optimal critrium
end
for j=1:1:N      %  iso_distance computation for a given prescribed value + NIV
   det = Vmp(j)*Vmp(j)-Vmm(j)*(VPP-Copt-NIV);  % discrminant of second order eqaution giving 2 gains k1 et k2
  if det>0;k1(j)=(Vmp(j)-sqrt(det))/Vmm(j);k2(j)=(Vmp(j)+sqrt(det))/Vmm(j);
  else; k1(j)=Vmp(j)/Vmm(j);k2(j)=k1(j);end    % plotting the valley       
end
%-----------------------------------------------------------------------------------
kopt
tauopt
Crit=Cmin*100/(tfin+1-tinit)/VPP
am=exp(-tech/tauopt);bm=1-am;  % optima model


%model opt simulation
for ii=2+r:1:tf;
   if ii<tinit | ii>tfin;sm(ii)=sp(ii);else;
      sm(ii)=sm(ii-1)*am+bm*e(ii-1-r)*kopt;end;
end

figure(3)
plot(w, sp,'r',w,sm,'k',w,e,'b')
grid,zoom,
title ('IDENTIFICATION / SM k / SP r / E b')
disp (' pause ')
pause

figure(4)
plot(taum,k1,'k',taum,k2,'k',taum,K,'b',tauopt,kopt,'r*', Tp, kp,'k*');
hold on
grid on
zoom on
text(tauopt,k2(jopt)+0.00002,['rm=', num2str(floor(rm))],'FontSize',12)
xlabel('Taum, Tauopt r*'),ylabel('k1,k2,K, Kopt *')
title(' PARAMETRIC SPACE / ISO DISTANCE /')

% NOTE
% in professionnal use less parameters are to be introduced
