% pid_pfc_limited
clear
close
tech=1; %sampling
tf=400;w=1:1:tf;
u=zeros(1, tf);
e=u; sm=u; sp=u;
tau=30; 
am=exp(-tech/tau); bm=1-am;
ap=exp(-tech/tau); bp=1-ap;
r=20;%delay
km=1;kp=1;
cons=100;
trbf=60;
lh=1-exp(-tech*3/trbf);
%pid
yp=u; ye=u;intg=u;
Kp=2;
Tr=15;
emax=input(' emax 110 ');
for ii= 2+r:1:tf;
   %process
   sp(ii)=sp(ii-1)*ap+bp*e(ii-1-r)*kp;
   %----------------------------------------------PFC
   %model independant
   sm(ii)=sm(ii-1)*am+bm*e(ii-1)*km;
   %control
   spred= sp(ii)+ sm(ii)-sm(ii-r);
   d=( cons-spred)*lh+sm(ii)*bm;
   e(ii)= d/(km*bm);
   if e(ii)>emax; e(ii)=emax;end
   %-------------------------------------------------PI
   %pid (with no dealy otherwise explosion r=0)
   yp(ii)=yp(ii-1)*ap+bp*ye(ii-1);
   er=cons-yp(ii);
   intg(ii)=intg(ii-1)+er;
   ye(ii)=Kp*er+(1/Tr)*intg(ii);
   if ye(ii)>emax; ye(ii)=emax;end
end
   figure (1)
   plot( w,e, 'k',w, sp,'r')
   grid
   zoom
  axis([0 tf 0 160])
  title ( ' PFC constraint   e k / sp r')
   figure (2)
   plot( w,ye, 'k',w, yp,'r')
   grid
   zoom
  axis([0 tf 0 160])
  title ( ' PID contrainte   ye k / yp r')
  