% RBF02.m
clear all;
close all;
kp=0.03;
ki=0.01;
kd=0.03;

yout_1=0;
u_1=0;
e_1=0;e_2=0;
ts=0.001;

for k=1:1:1000
    time(k)=k*ts;
    % reference input
    rin(k)=1.0*sign(sin(2*pi*k*ts));
    % nonlinear plant
    yout(k)=(-0.2*yout_1+0.5+u_1)/(1.5+yout_1^2);
    % error
    e(k)=rin(k)-yout(k);
    % PID
    xc(1)=e(k)-e_1;
    xc(2)=e(k);
    xc(3)=e(k)-2*e_1+e_2;
    du(k)=kp*xc(1)+ki*xc(2)+kd*xc(3);
    u(k)=u_1+du(k);
    % update yout(k-1),u(k-1)
    yout_1=yout(k);
    u_1=u(k);
    % update e(k-1),e(k-2)
    e_1=e(k);
    e_2=e_1;    
end
figure(1);
plot(time,rin,'b',time,yout,'r')
xlabel('time(s)');ylabel('rin,yout');