%模糊控制實例系統為20/(1.6s^2+4.4s+1)
%被控系統建模
num=20; %系統的分子
den=[1.6 4.4 1];%系統的分母
[a1,b,c,d]=tf2ss(num,den);%傳遞函數轉換到狀態空間
x=[0;0];
%系統參數
T=0.01;h=T;N=1000;
R=ones(1,N);%參考輸入
%定義輸入和輸出變量及隸屬度函數
a=newfis('Simple');
a=addvar(a,'input','e',[-4 4]);
a=addmf(a,'input',1,'NB','trimf',[-4,-4,-2]);
a=addmf(a,'input',1,'NS','trimf',[-4,-2,0]);
a=addmf(a,'input',1,'ZO','trimf',[-2,0,2]);
a=addmf(a,'input',1,'PS','trimf',[0,2,4]);
a=addmf(a,'input',1,'PB','trimf',[2,4,4]);
a=addvar(a,'input','de',[-4 4]);
a=addmf(a,'input',2,'NB','trimf',[-4,-4,-2]);
a=addmf(a,'input',2,'NS','trimf',[-4,-2,0]);
a=addmf(a,'input',2,'ZO','trimf',[-2,0,2]);
a=addmf(a,'input',2,'PS','trimf',[0,2,4]);
a=addmf(a,'input',2,'PB','trimf',[2,4,4]);
a=addvar(a,'output','u',[-4 4]);
a=addmf(a,'output',1,'NB','trimf',[-4,-4,-2]);
a=addmf(a,'output',1,'NS','trimf',[-4,-2,0]);
a=addmf(a,'output',1,'ZO','trimf',[-2,0,2]);
a=addmf(a,'output',1,'PS','trimf',[0,2,4]);
a=addmf(a,'output',1,'PB','trimf',[2,4,4]);
%模糊規則矩陣
rr=[5  5  4  4  3
    5  4  4  3  3
    4  4  3  3  2
    4  3  3  2  2
    3  3  2  2  1];
r1=zeros(prod(size(rr)),3);k=1;
for i=1:size(rr,1)
    for j=1:size(rr,2)
        r1(k,:)=[i,j,rr(i,j)];
        k=k+1;
    end
end
[r,s]=size(r1);
r2=ones(r,2);
rulelist=[r1,r2];
a=addrule(a,rulelist);
%采用模糊控制器的二階系統仿真
e=0;de=0;
ke=30;kd=5;ku=1;
for k=1:N
    %輸入變量變換至論域
    e1=ke*e;
    de1=kd*de;
    if e1>=4
        e1=4;
    elseif e1<=-4
        e1=-4;
    end
    if de1>=4
        de1=4;
    elseif de1<=-4
        de1=-4;
    end
    %模糊推理，計算出被控對象的控制輸入
    in=[e1 de1];
    u=ku*evalfis(in,a);
    uu(1,k)=u;
    if k<=ku
        u=0;
    else
        u=uu(1,k-ku);
    end
   
    %控制作用于被控系統，計算系統輸出；利用龍格—庫塔法進行系統仿真
    k0=a1*x+b*u;
    k1=a1*(x+h*k0/2)+b*u;
    k2=a1*(x+h*k1/2)+b*u;
    k3=a1*(x+h*k2)+b*u;
    x=x+(k0+2*k1+2*k2+k3)*h/6;
    y=c*x+d*u;
    yy(1,k)=y;
    %計算系統輸出誤差及誤差導數
    e1=e;e=y-R(1,k);
    de=(e-e1)/T;
end
%典型二階環節的模糊控制輸出曲線
kk=[1:N]*T;
figure(1);
plot(kk,R,'k',kk,yy,'r');grid on
xlabel('時間（秒）');ylabel('輸出');