www.pudn.com > stf.rar > OBSE5.M
>-----------------<<<Structure parameter>>>
kp=100; ti=0.4; td=0.1;
dt=0.2;
n=3; > state variable
m=1; > output variable
np=1; > input variable
nv=2; > process noise
ne=1; > measurement noise
mm=400;
>-----------------<<<System variable>>>
x1d=zeros(mm,1);
x=zeros(n,mm+1);
y=zeros(m,mm+1);
u=zeros(np,mm+1);
e=zeros(ne,mm+1);
v=zeros(nv,mm+1);
ep=zeros(m,mm+1);
>-----------------<<<Filter variable>>>
x_=zeros(n,mm+1);
vv0=zeros(m);
ff1=zeros(n);
hh=zeros(m,n);
rf=ones(n,1);
r1=zeros(m,1);
lmd=eye(n);
kk=zeros(n,m);
x1=zeros(n,1);
gama=zeros(n,nv);
beta=1;
rrr=zeros(1,mm+1);
mb1=zeros(mm,1);
threshold=(15)*ones(mm,1);
vv_=zeros(mm);
>-----------------<<<Teperaray variable>>>
lsnn1=zeros(n); lsnn2=zeros(n); lsnn3=zeros(n);
lsnn4=zeros(n); lsnn0=zeros(n);
lsnm1=zeros(n,m); lsnm2=zeros(n,m); lsnm3=zeros(n,m);
lsmn1=zeros(m,n); lsmn2=zeros(m,n); lsmn3=zeros(m,n);
lsmm1=zeros(m,m); lsmm2=zeros(m,m); lsmm3=zeros(m,m);
ls1m1=zeros(1,m); lsm11=zeros(1,m);
ls111=zeros(1,1);
u1=zeros(mm,np); y1=zeros(mm,1); y2=zeros(mm,1);
xx1=zeros(mm,1); xx3=zeros(mm,1);
xx2=zeros(mm,1); ; xx1_=zeros(mm,1); xx2_=zeros(mm,1); xx3_=zeros(mm,1);
xx_=zeros(n,1);
pp2=zeros(n); rr1=zeros(mm,1); rr2=zeros(mm,1);
b1=zeros(1,mm+1);
b2=zeros(1,mm+1);
>_________________<<<Filter parameter>>>
ra=0.95;
rf(n,n)=1.;
>_________________<<< FDD Algorithm Variable>>>
mb1=zeros(mm,1);
seta=zeros(mm,1);
>__________________<<< Parameters of the FDD Algorithm>>>
n3=5;
seta0=100; ro2=1*1;
n0=20; n1=10; n2=10;
>-----------------<<<Parameters of the system>>>
a1=0.0;
a2=0.0;
a3=0.08;
qq=zeros(nv);
qq(1,1)=a1*a1;
qq(2,2)=a2*a2;
rr=zeros(ne);
rr(1,1)=a3*a3;
>____________________<<< Initial value of the plant>>>
x(1,1)=0.2; x(2,1)=400; x(3,1)=100;
b1(1,1)=0;
b2(1,1)=0;
> -------------------<<<Initial value of the filter>>>
x_(1,1)=0.1; x_(2,1)=350; x_(3,1)=80;
pp=(100)*eye(n); > P(0|0)
>-----------------<<<Generation of random variable>>>
q11=100; v11=100; k011=exp(13.4); er11=5360;
dh11=-17835.821; rho11=1000; cp11=0.239;
ua11=11950; caf11=1; tf11=400;
for i=1: mm;
if i <= 200
x1d(i,1)=0.2;
else
x1d(i,1)=0.3;
end;
ep(1,i)=-(x1d(i,1)-x_(1,i));
if i==1
u(1,i)=kp*ep(1,i)+419;
elseif i==2
u(1,i)=kp*(ep(1,i)*(1+dt/ti)+ td*(ep(1,i)-ep(1,i-1))/dt)+419 ;
else
la0=kp*(1+dt/ti+td/dt);
la1=kp*(1+2*td/dt);
la2=kp*td/dt;
u(1,i)=u(1,i-1)+la0*ep(1,i)-la1*ep(1,i-1)+la2*ep(1,i-2);
end;
>u(1,i)=419;
if i <= 99;
b1(1,i+1)=0;
elseif (i >=100) &amt; (i <=109);
b1(1,i+1)=-4;
else
b1(1,i+1)=0;
end;
if i<= 199;
b2(1,i+1)=0;
elseif (i >=200) &amt; (i < 299);
b2(1,i+1)=+5/100;
else
b2(1,i+1)= 0;
end;
x(1,i+1)=x(1,i)+ dt*(x(3,i)/v11)*(caf11-x(1,i)) - dt*k011*exp(-er11/x(2,i))*x(1,i);
x(2,i+1)=x(2,i)+dt*(x(3,i)/v11)*(tf11-x(2,i))+ dt*(-dh11)*(k011/(rho11*cp11))*exp(-er11/x(2,i))*x(1,i) + dt*(ua11/(v11*rho11*cp11))*(u(1,i)-x(2,i));
x(3,i+1)=x(3,i)+b1(1,i+1)+b2(1,i+1);
y(1,i+1)=x(1,i+1)*x(2,i+1);
>[i u(1,i) x(1,i) x(2,i) x(3,i)]
>------------filtering begin:-------------------------------------
ff1(1,1)=1+dt*(-x_(3,i)/v11- k011*exp(-er11/x_(2,i))); > A(k)
ff1(1,2)=-dt*k011*x_(1,i)*exp(-er11/x_(2,i))*er11/(x_(2,i)*x_(2,i));
ff1(1,3)=dt*(caf11-x_(1,i))/v11;
ff1(2,1)=dt*(-dh11)*k011*exp(-er11/x_(2,i))/(rho11*cp11);
ff1(2,2)=1- dt*(x_(3,i)/v11 + dh11*k011*x_(1,i)*exp(-er11/x_(2,i))*er11/(rho11*cp11*x_(2,i)*x_(2,i)) + ua11/(v11*rho11*cp11) );
ff1(2,3)= dt*(tf11-x_(2,i))/v11;
ff1(3,3)=1;
x1(1,1)=x_(1,i)+dt*(x_(3,i)*(caf11-x_(1,i))/v11 - k011*exp(-er11/x_(2,i))*x_(1,i));
x1(2,1)=x_(2,i)+ dt*(x_(3,i)*(tf11-x_(2,i))/v11- dh11*k011*exp(-er11/x_(2,i))*x_(1,i)/(rho11*cp11) + ua11*(u(i)- x_(2,i))/(v11*rho11*cp11) );
x1(3,1)=x_(3,i);
hh(1,1)=x1(2,1); hh(1,2)=x1(1,1);
> f(x(k|k),..)
r1(1,1)=y(1,i+1)-x1(1,1)*x1(2,1);
if i==1
ls1m1=r1'; vv0=r1*ls1m1;
> lmd=tdfnm(n,m,qq,rr,vv0,ff1,pp,hh,rf,gama,beta);
rr1(i,1)=lmd(3,3);
[xx_,pp2,kk,pp1,vv_]=filternm(n,qq,rr,ff1,r1,pp,lmd,x1,hh,gama);
pp=pp2;
for j1=1:n; x_(j1,i+1)=xx_(j1,1); end;
> [i+1 x_(1,i+1) x_(2,i+1) x_(3,i+1) lmd(1,1)]
elseif i>= 2
ls1m1=r1'; lsmm1=r1*ls1m1; lsmm2=(ra)*vv0; lsmm3=lsmm2+lsmm1;
vv0=(1/(1+ra))*lsmm3;
> lmd=tdfnm(n,m,qq,rr,vv0,ff1,pp,hh,rf,gama,beta);
rr1(i,1)=lmd(3,3);
[xx_,pp2,kk,pp1,vv_]=filternm(n,qq,rr,ff1,r1,pp,lmd,x1,hh,gama);
pp=pp2;
for j1=1:n; x_(j1,i+1)=xx_(j1,1); end;
> [i+1 x_(1,i+1) x_(2,i+1) x_(3,i+1) lmd(1,1)]
end;
ccc=0;
for j1=1 : n;
ccc=ccc+(x(j1,i+1)-x_(j1,i+1))*(x(j1,i+1)-x_(j1,i+1));
end;
rr2(i,1)=sqrt(ccc);
end;
for j1=1:mm;
xx1(j1,1)=x(1,j1);
xx2(j1,1)=x(2,j1);
xx1_(j1,1)=x_(1,j1);
xx2_(j1,1)=x_(2,j1);
xx3(j1,1)=x(3,j1);
xx3_(j1,1)=x_(3,j1);
u1(j1,1)=u(1,j1);
y1(j1,1)=y(1,j1);
end;
j1=(1:mm)';
subplot(221);
plot(j1, [x1d xx1 xx1_]);
xlabel(' steps (a) ');
ylabel('CONCENTRATION');
subplot(222);
plot(j1, [xx2 xx2_ ] );
xlabel(' steps (b) ');
ylabel('TEMPERATURE');
subplot(223);
plot(j1, [xx3 xx3_]);
xlabel(' steps (c) ');
ylabel('FLOWRATE');
Subplot(224);
plot(j1, rr1);
xlabel(' steps (d) ');
ylabel(' FADING FACTOR ');
>print c:obse1 -dmeta
clear *.*;
clear;
>function lmd=tdfnm(n,m,qq,rr,vv0,ff1,pp,hh,rf,gama,beta)
>lsnn1=ff1'; lsnn2=pp*lsnn1; lsnn1=ff1*lsnn2; lsnm1=hh';
>lsnn2=lsnm1*hh; lsnn3=lsnn1*lsnn2; >> M(k+1)
>lsnnv1=gama*qq; lsnvn1=gama'; lsnn1=lsnnv1*lsnvn1;
>lsmn1=hh*lsnn1; lsnm1=hh'; lsmm1=lsmn1*lsnm1;
>lsmm3=(beta)*rr;
>lsmm2=lsmm1+ lsmm3;
>lsmm1=vv0-lsmm2; >> N(k+1)
>llm=0;
>for j1=1:n; llm=llm+rf(j1,1)*lsnn3(j1,j1); end;
>lln=0;
>for j1=1:m; lln=lln+lsmm1(j1,j1); end;
>c=lln/llm;
>for j1=1:n; b=rf(j1,1)*c;
> if b <= 1
> lmd(j1,j1)=1;
> else lmd(j1,j1)=b;
> end;
>end;
>function [xx_,pp2,kk,pp1,vv_]=filternm(n,qq,rr,ff1,r1,pp,lmd,x1,hh,gama)
>lsnn1=ff1'; lsnn0=pp*lsnn1; lsnn1=ff1*lsnn0; lsnn2=lmd*lsnn1;
>lsnnv1=gama*qq; lsnvn1=gama'; lsnn1=lsnnv1*lsnvn1;
>lsnn3=lsnn2+lsnn1;
>pp1=lsnn3; > P(k+1|k)
>lsnm1=hh'; lsnm2=lsnn3*lsnm1;
>lsmm1=hh*lsnm2; lsmm2=lsmm1+rr; lsmm1=inv(lsmm2);
> vv_=lsmm1;
>kk=lsnm2*lsmm1; >> k(k+1)
>lsnn1=eye(n);
>lsnn2=kk*hh; lsnn4=lsnn1-lsnn2;
>pp2=lsnn4*lsnn3;
>lsnn1=kk*r1; xx_=x1+lsnn1;
>function mb1=mb(mm,n0,n1, seta, seta0, ro2 )
>for j1=n0+n1+1:mm;
>c1=0;
>for j2=1:n1; c1=c1+seta(1,j1-j2,); end;
>c11=c1/n1;
>c2=0; c3=0;
>for j2=1:n1;
>c2=c2+(seta(1,j1-j2)-seta0)*(seta(1,j1-j2)-seta0);
>c3=c3+(seta(1,j1-j2)-c11)*(seta(1,j1-j2)-c11);
>end;
>c21=c2/(n1-1);
>c31=c3/(n1-1);
>mb1(j1,1)=c21/ro2 - log(c31/ro2) -1;
>[i j1 mb1(j1,1)]
>end;