clc
clear
%被控系统建模
num=20;
den=[8 6 1];
[a1,b,c,d]=tf2ss(num,den);%传递函数转换到状态空间
x=[0 0]';
%系统参数
T=0.01;
h=T;
N=1000;
%输入给定为1
R=ones(1,N);
%定义输入和输出变量及隶属度函数
a=newfis('Simple');
a=addvar(a,'input','e',[-6,6]); %误差 e
a=addmf(a,'input',1,'NB','trimf',[-6,-6,-4]);
a=addmf(a,'input',1,'NM','trimf',[-6,-4,-2]);
a=addmf(a,'input',1,'NS','trimf',[-4,-2,0]);
a=addmf(a,'input',1,'Z','trimf',[-2,0,2]);
a=addmf(a,'input',1,'PS','trimf',[0,2,4]);
a=addmf(a,'input',1,'PM','trimf',[2,4,6]);
a=addmf(a,'input',1,'PB','trimf',[4,6,6]);
a=addvar(a,'input','de',[-6,6]); %误差变化率 de
a=addmf(a,'input',2,'NB','trimf',[-6,-6,-4]);
a=addmf(a,'input',2,'NM','trimf',[-6,-4,-2]);
a=addmf(a,'input',2,'NS','trimf',[-4,-2,0]);
a=addmf(a,'input',2,'Z','trimf',[-2,0,2]);
a=addmf(a,'input',2,'PS','trimf',[0,2,4]);
a=addmf(a,'input',2,'PM','trimf',[2,4,6]);
a=addmf(a,'input',2,'PB','trimf',[4,6,6]);
a=addvar(a,'output','u',[-6,6]); %控制量 u
a=addmf(a,'output',1,'NB','trimf',[-6,-6,-4]);
a=addmf(a,'output',1,'NM','trimf',[-6,-4,-2]);
a=addmf(a,'output',1,'NS','trimf',[-4,-2,0]);
a=addmf(a,'output',1,'Z','trimf',[-2,0,2]);
a=addmf(a,'output',1,'PS','trimf',[0,2,4]);
a=addmf(a,'output',1,'PM','trimf',[2,4,6]);
a=addmf(a,'output',1,'PB','trimf',[4,6,6]);
%模糊规则矩阵
rulelist=[1 1 7 1 1;
1 2 7 1 1;
1 3 6 1 1;
1 4 6 1 1;
1 5 5 1 1;
1 6 5 1 1;
1 7 4 1 1;
2 1 7 1 1;
2 2 6 1 1;
2 3 6 1 1;
2 4 5 1 1;
2 5 5 1 1;
2 6 4 1 1;
2 7 3 1 1;
3 1 6 1 1;
3 2 6 1 1;
3 3 5 1 1;
3 4 5 1 1;
3 5 4 1 1;
3 6 3 1 1;
3 7 3 1 1;
4 1 6 1 1;
4 2 5 1 1;
4 3 5 1 1;
4 4 4 1 1;
4 5 3 1 1;
4 6 3 1 1;
4 7 2 1 1;
5 1 5 1 1;
5 2 5 1 1;
5 3 4 1 1;
5 4 3 1 1;
5 5 3 1 1;
5 6 2 1 1;
5 7 2 1 1;
6 1 5 1 1;
6 2 4 1 1;
6 3 3 1 1;
6 4 3 1 1;
6 5 2 1 1;
6 6 2 1 1;
6 7 1 1 1;
7 1 4 1 1;
7 2 3 1 1;
7 3 3 1 1;
7 4 2 1 1;
7 5 2 1 1;
7 6 1 1 1;
7 7 1 1 1];
a=addrule(a,rulelist);
%采用模糊控制器的二阶系统仿真
e=0;de=0;
ke=10;kd=15;ku=1;
for k=1:N
%输入变量变换至论域
e1=ke*e;
de1=kd*de;
if e1>=6
e1=6;
elseif e1<=-6
e1=-6;
end
if de1>=6
de1=6;
elseif de1<=-6
de1=-6;
end
%模糊推理，计算出被控对象的控制输入
in=[e1 de1];
u=ku*evalfis(in,a);
uu(1,k)=u;
%控制作用于被控系统，计算系统输出；利用龙格—库塔法进行系统仿真
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('输出');
