 %pfc_ramp
 clear
 close
%horizon
tf=2500;tech=1;w=1:1:tf;
%model process level1
u=ones(1,tf)*50; cons=u; sp=u;sm=u;
dc=u*0;e=u;Cons=u;pert=u*0;
uu=u*0;c1=u*0;c2=c1;uuu=u*0;
%-------------
tm=80; am=exp(-tech/tm);bm=1-am;
tp=80; ap=exp(-tech/tp);bp=1-ap; 
km=1;
kp=1;
% coincidence horizon
h2=10;
h1=2; 
%rd=input('rd (10)'); %process delay 
rd=10*1;%------------------------------------------------------RETARD----------------
trbf=40;la=exp(-tech*3/trbf);
%coincidence  free mode
lah1=1-la^h1;lah2=1-la^h2;bmh1=1-am^h1;bmh2=1-am^h2;
%slope of setpoint
p=0.0217;%pente
vl=4*p;
H=1;
% set point design with 2 docking curves for connection------------------------------
%---

TEMPINIT=50;TEMPFIN=20;%temp init and final
TRINIT=650;
TRFIN= TRINIT+ ((TEMPINIT-TEMPFIN)/p); %date init and final
TAC1=400;%  soft start
DT=TRINIT-TAC1;% docking to the ramp
TAC2=TRFIN-DT;% smooth landing
%
DAC=TEMPINIT-p*(TAC2-TRINIT)-TEMPFIN;
tau2=DAC/p;
K1=p/(2*(TRINIT+DT-TAC1));
K2=-p/(2*(TAC2-TRFIN -DT));

for ii = 2:1:tf;
   if ii< TRINIT; cons(ii)=TEMPINIT;end
   if ii>=TRINIT; cons(ii)=TEMPINIT-p*(ii-TRINIT);end
   if ii>=TRFIN;  cons(ii)=TEMPINIT-p*(-TRINIT+TRFIN);end
   %---
   if ii<TAC1; C(ii)=TEMPINIT;end
   if ii>=TAC1 & ii<TRINIT+DT;C(ii)=TEMPINIT-K1*(ii-TAC1)^2;end
   if ii>=TRINIT+DT;C(ii)=TEMPINIT-p*(ii-TRINIT );end
   if ii>=TAC2; C(ii)=TEMPFIN +K2*(ii-TRFIN-DT)^2;end
   if ii>=TRFIN+DT; C(ii)=TEMPFIN;end
   dc(ii)=C(ii)-C(ii-1);
end

Cons=C;
figure(1);
plot(w,cons, 'r',w,Cons,'k');
grid;zoom;
axis ([1 tf 15 55])
%---------------------------------------------------- end of set point proG

 
%break
%response to base functions yn :ramp z :step
zc1=zeros(1, h2+2);zc2=zc1;
for ii=2:1:h2+2;%  2 methods
   t(ii)=1;tt(ii)=(ii-1);
   zc1(ii)=zc1(ii-1)*am+bm*t(ii-1);
   zc2(ii)=zc2(ii-1)*am+bm*tt(ii-1);
end
%--- may simpler procedur
   y1=h1*km-km*(1-am^h1)/(1-am);
   z1=km*(1-am^h1);  
   y2=h2*km-km*(1-am^h2)/(1-am);
   z2=km*(1-am^h2);
%----
eo=e;spo=sp;smo=sm;
lh=1-exp(-tech*3/120);

%real time
for ii=2+rd:1:tf-50;
   if ii>250; pert(ii)=0;else pert(ii)=0; end
   if ii>400; pert(ii)=0;end
 %simulation processus
 sp(ii)=sp(ii-1)*ap+bp*(e(ii-1-rd)+pert(ii))*kp;%---------------------
 %control
   sm(ii)=sm(ii-1)*am+bm*e(ii-1)*km;
   spest=sp(ii)+sm(ii)-sm(ii-rd);%estimation of predicted CV
   %incremental targets at h1 et h2
   %retard-----------------------------------------rd--------------------?
   d1=(Cons(ii+h1)-Cons(ii))+(Cons(ii)-spest+dc(ii)*rd)*lah1+sm(ii)*bmh1;% 
   d2=(Cons(ii+h2)-Cons(ii))+(Cons(ii)-spest+dc(ii)*rd)*lah2+sm(ii)*bmh2;%
   e(ii)=(y2*d1-y1*d2)/(y2*z1-y1*z2);%coef step
   %constraint
   if e(ii)>e(ii-1)+vl  ; e(ii)=e(ii-1)+vl;end
   %constraint
   if e(ii)<e(ii-1)-vl ; e(ii)=e(ii-1)-vl;end  

end
TCENT=floor(TAC2);
ER=Cons(TCENT)-sp(TCENT)
ww=2*w;
E=input(' E   ');% error E to evaluate how PFC follows up a ramp with no lag !
figure(2);
plot(ww,e,'b',ww,Cons,'k',ww,pert*2+10,'g',ww,sp+E,'r',ww,spo*0,'g');
hold on
grid on
zoom on
axis([1 2*tf 15 55])

%<<<<<<<<<<<<<FOLLOW THE SET POINT ZITH NO DELAY . unique property of PFC /////

title ( ' SETPOINT k / Tmass r / TECalo b  / EXOthermicity g / no error' )
text (50, 51 , 'd°')
text(1210, 17 , 'minutes')

