%CMAC and PID Concurrent Control
function CMAC=CMAC_PID()
clear all;
close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%system%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ts=0.001;                  %采样时间越短，控制精度越难
sys=tf(1770,[1,60,1770]);
dsys=c2d(sys,ts,'z');
[num,den]=tfdata(dsys,'v');
%%%%%%%%%%%%%%%%%%%%initial vallue and set the parameter of system %%%%%%%%
alfa=0.04;   xite=0.10;
N=100;C=5;
w=zeros(N+C,1);
tmax=1;
kp=25; ki=0.0; kd=0.28;
%%%%%%%%%%%%%%%%%%%%%%loop needed value%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w_1=w;w_2=w;y_1=0;y_2=0;y_3=0;
u_1=0.0;u_2=0.0;u_3=0.0;
x=[0,0,0]';error_1=0;

%%%%%%%%%%%%%%%%%%%%%%%%主程序数据准备%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%概念映射，将输入空间rin<【-5 5】进行量化，概念映射
A=0.50;   Smin=-A; Smax=A;
dvi=(Smax-Smin)/(N-1);
for i=1:1:C                %C size
    v(i)=Smin;
end
for i=C+1:1:C+N            %N size
    v(i)=v(i-1)+dvi;
end
for i=N+C+1:1:N+2*C        %C size
    v(i)=Smax;
end
%%%%%%%%%%%%%%%%%%%%%%%%%主程序——仿真%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 

for k=1:1:tmax/ts
%%%%%%%%%%%%%%%%%%%%%be controled model%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
rin=0.5*sin(2*2*pi*k*ts);          % rin(k)=A*sign(sin(2*2*pi*k*ts));   跟踪难度大于正弦
yout=-den(2)*y_1-den(3)*y_2+num(2)*u_1+num(3)*u_2;
error=rin-yout;  
%%%%%%%%&&&&&&&&&&&&&&CMAC神经网络的输出 %%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%% 
for i=1:1:N+C
if rin>=v(i)&rin<=v(i+C);   a(i)=1;else;  a(i)=0;end
end
un=a*w;
%%%%%%%%%%%%%%%%%%%%%%%%%PID输出及总得控制输出%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%CMAC的实际映射
up=kp*x(1)+kd*x(2)+ki*x(3);                           %PID
u=up+un;                                              %如果不加UN为仅在PID作用下控制
%if k==150 ;   u(k)=u(k)+5.0;end  %disturbance
if u>=10;   u=10;end;if u<=-10;   u=-10; end           %输出阀值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%CMAC net learing%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
d_w=a'*xite*(u-un)/C;             %performance index is u-un
w=w_1+ d_w+alfa*(w_1-w_2);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%parameter update%%%%%%%%%%%%%%%%%%%%%%%%%
 w_3=w_2;w_2=w_1;w_1=w;u_2=u_1;u_1=u;y_2=y_1;y_1=yout;
x(1)=error;                % Calculating P
x(2)=(error-error_1)/ts;   % Calculating D
x(3)=x(3)+error*ts;
error_2=error_1;error_1=error;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%output save%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 time(k)=ts*k;                  %record the time
 RR(k)=rin; YY(k)=yout; UN(k)=un;U(k)=u;
end
figure(1)
plot(time,RR,'b',time,YY,'r');grid on
%figure(2)
%plot(time,UN,'b',time,U,'r');grid on
end