pidsimu.m

蚁群算法优化ＰＩＤ参数．．只运行ant.m这个文件就行了

function [y]=PIDsimu(kp,ti,td)
a(1,1)=1;b(1,1)=1;c(1,1)=1;d(1,1)=2;cp(1,1)=1.0;
a(1,2)=1;b(1,2)=1;c(1,2)=2;d(1,2)=2;cp(1,2)=2.0;
a(1,3)=1;b(1,3)=1;c(1,3)=3;d(1,3)=5;cp(1,3)=3.0;
a(2,1)=1;b(2,1)=2;c(2,1)=2;d(2,1)=0;cp(2,1)=2.0;
a(2,2)=1;b(2,2)=2;c(2,2)=3;d(2,2)=0;cp(2,2)=3.0;
a(2,3)=1;b(2,3)=2;c(2,3)=4;d(2,3)=0;cp(2,3)=4.0;


h=0.1;T=1;
%kp=[5,5;5,5;5,5];ti=[1,1;1,1;1,1];td=[0.4821,0.5;0.4,0.5;0.4,0.5];
%lam1=1;lam2=1;
for i=1:3
    for j=1:2002
        u(i,j)=0;
    end
end
for i=1:2
   for j=1:2002
       e(i,j)=0;y(i,j)=0;ef(i,j)=0;%e为控制器输入，ef为输出偏差；
   end
end
for i=1:2
   for j=1:2002
       u1(i,j)=0;u2(i,j)=0;u3(i,j)=0;
   end
end
for i=1:2
    for j=1:3
     %xm1(i,j)=0;xm2(i,j)=0;
     x(i,j)=0;x1(i,j)=0;
 end
end
             for j=1:3               
             for i=1:d(1,j)/h
                ys1(j,i)=0;
                ysm1(j,i)=0;
               end
           end
              for j=1:3               
             for i=1:d(2,j)/h
                ys1(j,i)=0;
                ysm1(j,i)=0;
               end
           end
%龙格库塔求控制作用
for m=1:200/h
    %for n=1:1:2
        if m==1
            for n=1:2
            e(n,m)=1;
            end
        else
            for n=1:2
               e(n,m)=1-y(n,m-1);
            end%求e1 and e2
        end
    %end
    for i=1:2
           u1(i,m+2)=u1(i,m+1)+kp(1,i)*(e(i,m+2)-e(i,m+1))+T/ti(1,i)*e(i,m+2)+td(1,i)/T*(e(i,m+2)-2*e(i,m+1)+e(i,m));
           u2(i,m+2)=u2(i,m+1)+kp(2,i)*(e(i,m+2)-e(i,m+1))+T/ti(2,i)*e(i,m+2)+td(2,i)/T*(e(i,m+2)-2*e(i,m+1)+e(i,m));
           u3(i,m+2)=u3(i,m+1)+kp(3,i)*(e(i,m+2)-e(i,m+1))+T/ti(3,i)*e(i,m+2)+td(3,i)/T*(e(i,m+2)-2*e(i,m+1)+e(i,m));
           %ut=u(m+2);
           % u(i,m)=u(i,m)+f4(i,j)*xu(i,j)+f2(i,j)*e(j,m);%u1/2=g1*e1+g2*e2
    end
    for i=1:2
        u(1,m+2)=u(1,m+2)+u1(i,m+2);u(2,m+2)=u(2,m+2)+u2(i,m+2);u(3,m+2)=u(3,m+2)+u3(i,m+2);
    end
    %龙格库塔求输出          
    %真实对象状态
  for i=1:2
      for j=1:3
            for k=1:1
                k1(i,j)=h*(-b(i,j)*x(i,j)+u(j,m+2));
                k2(i,j)=h*(-b(i,j)*(x(i,j)+0.5*k1(i,j))+u(j,m+2));%为什么是0.5？
                k3(i,j)=h*(-b(i,j)*(x(i,j)+0.5*k2(i,j))+u(j,m+2));%为什么是0.5？
                k4(i,j)=h*(-b(i,j)*(x(i,j)+k3(i,j))+u(j,m+2));
                x(i,j)=x(i,j)+(k1(i,j)+2*k2(i,j)+2*k3(i,j)+k4(i,j))/6;
            end
        end
    end
            %switch i
            %case '1'

            for j=1:3
                ys1(j,m+d(1,j)/h)=cp(1,j)*x(1,j);

                ys2(j,m+d(2,j)/h)=cp(2,j)*x(2,j);
            end
            
                y(1,m)=ys1(1,m)+ys1(2,m)+ys1(3,m);
                y(2,m)=ys2(1,m)+ys2(2,m)+ys2(3,m);
end

%plot(y(1,:),'r');
%hold on;
%plot(y(2,:));