clc
close all;
clear all;
tc=newfis('tc','sugeno');
tc=addvar(tc,'input','theta',[-90 90]);
tc=addmf(tc,'input',1,'about 0','trimf',[-90 0 90]);
tc=addmf(tc,'input',1,'about 90','trimf',[0 90 90]);
tc=addmf(tc,'input',1,'about-90','trimf',[-90 -90 0]);

tc=addvar(tc,'input','omega',[-5 5]);
tc=addmf(tc,'input',2,'NG','gaussmf',[1.8 -5]);
tc=addmf(tc,'input',2,'ZR','gaussmf',[1.8 0]);
tc=addmf(tc,'input',2,'PO','gaussmf',[1.8 5]);

tc=addvar(tc,'output','u',[-2 2]);
tc=addmf(tc,'output',1,'No.1','linear',[126.07 26.3 0]);
tc=addmf(tc,'output',1,'No.2','linear',[2722.3 883.4 0]);

rulelist=[1 0 1 1 1;2 0 2 1 1;3 0 2 1 1];
tc=addrule(tc,rulelist);


model=newfis('model','sugeno');
model=addvar(model,'input','theta',[-90 90]);
model=addmf(model,'input',1,'about 0','trimf',[-90 0 90]);
model=addmf(model,'input',1,'about 90','trimf',[0 90 90]);
model=addmf(model,'input',1,'about-90','trimf',[-90 -90 0]);

model=addvar(model,'input','omega',[-5 5]);
model=addmf(model,'input',2,'NG','gaussmf',[1.8 -5]);
model=addmf(model,'input',2,'ZR','gaussmf',[1.8 0]);
model=addmf(model,'input',2,'PO','gaussmf',[1.8 5]);

model=addvar(model,'input','u',[-5 5]);
model=addmf(model,'input',3,'any','gaussmf',[1.5 -5]);

model=addvar(model,'output','d_theta',[-2 2]);
model=addmf(model,'output',1,'No.1','linear',[0 1 0 0]);
model=addmf(model,'output',1,'No.2','linear',[0 1 0 0]);

model=addvar(model,'output','d_omega',[-2 2]);
model=addmf(model,'output',2,'No.1','linear',[17.2941 0 -0.1765 0]);
model=addmf(model,'output',2,'No.2','linear',[14.4706 0 -0.1765 0]);
rulelist1=[1 0 0 1 1 1 1;2 0 0 2 2 1 1;3 0 0 2 2 1 1];
model=addrule(model,rulelist1);

showrule(tc)
showrule(model)

N=1000;h=0.005; %N为迭代次数，h为步长
x=[1.309;0];    %设定系统的初值
y=zeros(2,N);
u=zeros(1,N);
for k=2:N               %R-K迭代法
u(k)=evalfis([x(1),x(2)],tc);
k0=evalfis([x(1),x(2),u(k)],model)';
x1=x+h*k0/2;
k1=evalfis([x1(1),x1(2),u(k)],model)';
x1=x+h*k1/2;
k2=evalfis([x1(1),x1(2),u(k)],model)';
x1=x+h*k2;
k3=evalfis([x1(1),x1(2),u(k)],model)';
x=x+(k0+2*k1+2*k2+k3)*h/6;
y(:,k)=x;
end
kk=[1:N]*h;
figure
plot(kk,y),grid on
xlabel('x'),ylabel('系统状态')
axis([0 5 -2 1.5])

figure
plot(kk,u),grid on
xlabel('u(k)'),ylabel('模糊控制器输出')
axis([0 5 0 250])