t s模糊辨识 matlab,T-S模糊辨识与广义预测控制(GPC)MATLAB源程序

T-S模糊辨识与广义预测控制(GPC)MATLAB源程序

clear all

close all

clear all

%获得样本数据

y(1)=1;y(2)=1;u(1)=0;r=16;N0=500;yr1=1;

for k=3:N0;

u(k-1)=5*sin(pi*(k-1)/100);

y(k)=(y(k-1)*y(k-2)*(y(k-1)-2.5)/(1+y(k-1)^2+y(k-2)^2))+u(k-1);

end

U1=zeros(4,N0);

for k=1:N0;

if y(k)<=-5

U1(4,k)=-0.1*y(k)-0.5;

%%%%%%%

U1(3,k)=0.2*y(k)+2;

elseif y(k)<=0

U1(3,k)=-0.2*y(k);

U1(2,k)=0.2*y(k)+1;

elseif y(k)<=5

U1(2,k)=-0.2*y(k)+1;

U1(1,k)=0.1*y(k);

else

U1(1,k)=0.1*y(k);

end

end

for i=1:4;

den=0.0001;

for k=1:N0

den=U1(i,k)+den;

end

C(i)=(U1(i,:)*y')/den;

end

for k=1:N0;

for i=1:4;

D(i,k)=abs(y(k)-C(i));

end

if (D(:,k)>0)

for i=1:4;

den=0.00001;

for j=1:4;

den=D(i,k)/D(j,k)+den;

end

U2(i,k)=1/den;

end

else

flag=(D(:,k)>0);

for i=1:4

if (flag(i)==1)

U2(i,k)=0;

else

U2(i,k)=1/sum(flag);

end

end

end

end

for i=1:4

den=0.00001;

for k=1:N0

den=U2(i,k)+den;

end

C(i)=(U2(i,:)*y')/den;

end

for k=1:N0

den=0.00001;

for i=1:4

den=U2(i,k)+den;

end

for i=1:4

Beta(i,k)=U2(i,k)/den;

end

end

for k=3:N0

fai(:,k-2)=[Beta(:,k-2)

Beta(:,k-2)*y(k-1)

Beta(:,k-2)*y(k-2)

Beta(:,k-2)*u(k-1)];

end

ft=fai';

w=zeros(N0-2,16);

w(:,1)=ft(:,1);

for k=1:N0-2

for h=2:16

for i=1:h-1

den(1,i)=0.001;

for k=1:N0-2

den(1,i)=den(1,i)+w(k,i);

end

alpha(h,i)=0;

for k=1:N0-2

alpha(h,i)=alpha(h,i)+w(k,i)*ft(k,h)/den(1,i);

end

end

s(1,h)=0;

for i=1:h-1

s(1,h)=s(1,h)+alpha(h,i)*w(k,i);

end

w(k,h)=ft(k,h)-s(1,h);

end

end

W=w';

at=alpha';

for i=1:r;

den=0.0001;

for k=1:N0-2

den=W(i,k)+den;

end

g(i)=(W(i,:)*y(1:N0-2)')/den;

end

theta(r)=g(r);

for i=(r-1):(-1):1;

theta(i)=g(i)-at(i,:)*theta';

end

x=[-15:0.1:15];

for k=1:301

if x(k)<=-10

UNB(k)=-0.1*x(k)-0.5;

UNS(k)=0;

UPS(k)=0;

UPB(k)=0;

elseif x(k)<=-5

UNB(k)=-0.1*x(k)-0.5;

UNS(k)=0.2*x(k)+2;

UPS(k)=0;

UPB(k)=0;

elseif x(k)<=0

UNB(k)=0;

UNS(k)=-0.2*x(k);

UPS(k)=0.2*x(k)+1;

UPB(k)=0;

elseif x(k)<=5

UNB(k)=0;

UNS(k)=0;

UPS(k)=-0.2*x(k)+1;

UPB(k)=0.1*x(k);

else

UNB(k)=0;

UNS(k)=0;

UPS(k)=0;

UPB(k)=0.1*x(k);

if UPB(k)>=1

UPB(k)=1;

end

end

end

hold on;

plot(x,UNB,'b');

plot(x,UNS,'k');

plot(x,UPS,'r');

plot(x,UPB,'g');

subplot(2,1,1):plot(x,UNB,'b',x,UNS,'k',x,UPS,'r',x,UPB,'g');

xlabel('隶属函数图');;ylabel('隶