%DRNN Tunning PID Controler for coupling plant
clear all;
close all;

u1_1=0.0;u1_2=0.0;u1_3=0.0;
u2_1=0.0;u2_2=0.0;u2_3=0.0;
y1_1=0;y2_1=0;

wd1=rands(7,1);
wo1=rands(7,1);
wi1=rands(3,7);
wd1_1=wd1;wo1_1=wo1;wi1_1=wi1; 

x1=zeros(7,1);
x1_1=x1;

wd2=rands(7,1);
wo2=rands(7,1);
wi2=rands(3,7);
wd2_1=wd2;wo2_1=wo2;wi2_1=wi2; 

xitei1=0.40;xited1=0.40;xiteo1=0.40;
xitei2=0.40;xited2=0.40;xiteo2=0.40;
alfa1=0.04;alfa2=0.04;

x2=zeros(7,1);
x2_1=x2;

error1_1=0;error1_2=0;
error2_1=0;error2_2=0;

kp1=0.30;ki1=0.150;kd1=0.2;
kp2=0.30;ki2=0.150;kd2=0.2;

kp1_1=kp1;kd1_1=kd1;ki1_1=ki1;  
kp2_1=kp2;kd2_1=kd2;ki2_1=ki2;  

xitekp1=0.5;xitekd1=0.3;xiteki1=0.001;
xitekp2=0.5;xitekd2=0.3;xiteki2=0.0001;

ei1=0;ei2=0;
ts=1;
for k=1:1:1500
time(k)=k*ts;

%Step Signal
R=[1;0];
%R=[0;1];

%Coupling Plant
yout1(k)=1.0/(1+y1_1)^2*(0.8*y1_1+u1_2+0.2*u2_3);
yout2(k)=1.0/(1+y2_1)^2*(0.9*y2_1+0.3*u1_3+u2_2);

In1=[u1_1,yout1(k),1]';
for j=1:1:7
    si(j)=In1'*wi1(:,j)+wd1(j)*x1(j);
end

for j=1:1:7
    x1(j)=(1-exp(-si(j)))/(1+exp(-si(j)));
end

Pi=0*x1;
for j=1:1:7
    Pi(j)=wo1(j)*(1+x1(j))*(1-x1(j))*x1_1(j);
end

Qi=0*wi1;
for j=1:1:7
   for i=1:1:3
       Qi(i,j)=wo1(j)*(1+x1(j))*(1-x1(j))*In1(i);
   end
end   

ym=0;
for j=1:1:7
    ym=ym+x1(j)*wo1(j);
end
ymout1(k)=ym;

wo1=wo1+xiteo1*(yout1(k)-ymout1(k))*x1+alfa1*(wo1-wo1_1);
wd1=wd1+xited1*(yout1(k)-ymout1(k))*Pi+alfa1*(wd1-wd1_1); 
for j=1:1:7
   if abs(wd1(j))>1
      wd1(j)=0.5*sign(wd1(j));
   end
end   
wi1=wi1+xitei1*(yout1(k)-ymout1(k))*Qi+alfa1*(wi1-wi1_1);   
  
yu=0;
for j=1:1:7
   yu=yu+wo1(j)*wi1(1,j)*(1+x1(j))*(1-x1(j));
end   
dyout1(k)=yu;

error1(k)=R(1)-yout1(k);
  
xc1(1)=error1(k);             %Calculating P
xc1(2)=error1(k)-error1_1;    %Calculating D
ei1=ei1+error1(k)*ts;
xc1(3)=ei1;                   %Calculating I
      
kp1(k)=kp1_1+xitekp1*error1(k)*xc1(1);
kd1(k)=kd1_1+xitekd1*error1(k)*xc1(2);
ki1(k)=ki1_1+xiteki1*error1(k)*xc1(3);  

kp1(k)=kp1_1+xitekp1*error1(k)*xc1(1)*dyout1(k);
kd1(k)=kd1_1+xitekd1*error1(k)*xc1(2)*dyout1(k);
ki1(k)=ki1_1+xiteki1*error1(k)*xc1(3)*dyout1(k);  

if kp1(k)<0
   kp1(k)=0;
end
if kd1(k)<0
   kd1(k)=0;
end
if ki1(k)<0
   ki1(k)=0;
end
   
u1(k)=kp1(k)*xc1(1)+kd1(k)*xc1(2)+ki1(k)*xc1(3); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In2=[u2_1,yout2(k),1]';

for j=1:1:7
   si(j)=In1'*wi2(:,j)+wd2(j)*x2(j);
end

for j=1:1:7
   x2(j)=(1-exp(-si(j)))/(1+exp(-si(j)));
end

Pi=0*x2;
for j=1:1:7
   Pi(j)=wo2(j)*(1+x2(j))*(1-x2(j))*x2_1(j);
end

Qi=0*wi2;
for j=1:1:7
   for i=1:1:3
      Qi(i,j)=wo2(j)*(1+x2(j))*(1-x2(j))*In2(i);
   end
end   

ym=0;
for j=1:1:7
   ym=ym+x2(j)*wo2(j);
end
ymout2(k)=ym;

wo2=wo2+xiteo2*(yout2(k)-ymout2(k))*x2+alfa2*(wo2-wo2_1);
wd2=wd2+xited2*(yout2(k)-ymout2(k))*Pi+alfa2*(wd2-wd2_1);
  for j=1:1:7
     if abs(wd2(j))>1
        wd2(j)=0.5*sign(wd2(j));
     end
  end   
wi2=wi2+xitei2*(yout2(k)-ymout2(k))*Qi+alfa2*(wi2-wi2_1);   
  
yu=0;
for j=1:1:7
    yu=yu+wo2(j)*wi2(1,j)*(1+x2(j))*(1-x2(j));
end   
dyout2(k)=yu;

error2(k)=R(2)-yout2(k);

xc2(1)=error2(k);             %Calculating P
xc2(2)=error2(k)-error2_1;    %Calculating D
ei2=ei2+error2(k)*ts;
xc2(3)=ei2;                   %Calculating I
      
kp2(k)=kp2_1+xitekp2*error2(k)*xc2(1)*dyout2(k);
kd2(k)=kd2_1+xitekd2*error2(k)*xc2(2)*dyout2(k);
ki2(k)=ki2_1+xiteki2*error2(k)*xc2(3)*dyout2(k);  

if kp2(k)<0
   kp2(k)=0;
end
if kd2(k)<0
   kd2(k)=0;
end
if ki2(k)<0
   ki2(k)=0;
end
   
u2(k)=kp2(k)*xc2(1)+kd2(k)*xc2(2)+ki2(k)*xc2(3);

%Return of PID parameters
error1_2=error1_1;
error1_1=error1(k);

error2_2=error2_1;
error2_1=error2(k);

wd1_1=wd1;wo1_1=wo1;wi1_1=wi1; 
wd2_1=wd2;wo2_1=wo2;wi2_1=wi2; 

u1_4=u1_3;u1_3=u1_2;u1_2=u1_1;u1_1=u1(k);
u2_4=u2_3;u2_3=u2_2;u2_2=u2_1;u2_1=u2(k);

y1_1=yout1(k);
y2_1=yout2(k);
end
figure(1);
plot(time,dyout1,'r');
xlabel('time(s)');ylabel('Jacobian');
hold on;
plot(time,dyout2,'b');
xlabel('time(s)');ylabel('Jacobian');
figure(2);
plot(time,R(1),'b',time,yout1,'r');
hold on;
plot(time,R(2),'r',time,yout2,'b');
xlabel('time(s)');ylabel('rin,yout');
figure(3);
plot(time,kp1,'r',time,ki1,'b',time,kd1,'k');
xlabel('time(s)');ylabel('kp1 ki1 kd1');
figure(4);
plot(time,kp2,'r',time,ki2,'b',time,kd2,'k');
xlabel('time(s)');ylabel('kp2 ki2 kd2');