关于(微粒群)pso优化pid

hunan

小女子我在做毕业设计时,遇到了个大难题,是关于pso优化pid的,有哪位高手可以指点下么? 万分感谢...



具体题目:

利用标准PSO算法（带惯性权重的PSO算法）、带收缩因子的PSO算_\带交叉因子的改进PSO算法   分别实现PID控制器的控制参数的优

化与计算机仿真









[ 本帖最后由 hunan 于 2007-5-19 18:22 编辑 ]

风花雪月

有什么具体问题？

hunan

就是要个matlab的程序 优化一个pid的函数

hunan

你qq多少类 ..加下我吧 39383685

风花雪月

1. clear;
   
2. x0=[2.813,1.719,1.151];%x=[Kp,Ki,Kd]
   
3. c1=1.495;c2=1.495;n=3;group=50;Dmax=100;
   
4. % Vmm=[0.4 1.2;0.35 1.05;0.05 0.15];
   
5. Vmm=[-2.5 2.5;-2.5 2.5;-2.5 2.5];%可以为负吗
   
6. 
   
7. Xmm=[0 5;0 5;0 5];
   
8. X=zeros(n,group);V=zeros(n,group);
   
9. % omega=zeros(1,group);
   
10. % e=0.001;
    
11. 
    
12. X(1,:)=x0(1)+5*rand(1,group);V(1,:)=Vmm(1,1)+5*rand(1,group);
    
13. X(2,:)=x0(2)+5*rand(1,group);V(2,:)=Vmm(2,1)+5*rand(1,group);
    
14. X(3,:)=x0(3)+5*rand(1,group);V(3,:)=Vmm(3,1)+5*rand(1,group);
    
15. P_p=X;globe=zeros(n,Dmax);
    
16. 
    
17. %%%评价每个粒子适应值，寻找出 P_globle
    
18. for j=1:group
    
19.     xx=X(:,j)';
    
20. %     fz(j)=J_ITAE(xx);
    
21.     fz(j)=J_ITAE_discrete(xx);
    
22. end
    
23. [P_g,I]=min(fz);%P_g  1*1 ?
    
24. globe(:,1)=X(:,I);
    
25. 
    
26. for itrtn=1:Dmax
    
27.     omega=0.739;
    
28.     r1=rand(1);r2=rand(1);
    
29.     for i=1:group
    
30.         V(:,i)=omega*V(:,i)+c1*r1*(P_p(:,i)-X(:,i))+c2*r2*(globe(itrtn)-X(:,i));
    
31.     end
    
32.     X=X+V;%先虑出位置，如果先是速度虑出，则位置虑出的受速度影响。
    
33.     for i=1:group
    
34.     %compare each dimension of V
    
35.         for row=1:n
    
36.             if X(row,i)>Xmm(row,2)
    
37.                 X(row,i)=Xmm(row,2);
    
38.             elseif X(row,i)<Xmm(row,1)
    
39.                 X(row,i)=Xmm(row,1);
    
40.             else
    
41.             end
    
42.         end
    
43.     end
    
44.     for i=1:group
    
45.         %compare each dimension of V
    
46.         for row=1:n
    
47.             if V(row,i)>Vmm(row,2)
    
48.                 V(row,i)=Vmm(row,2);
    
49.             elseif V(row,i)<Vmm(row,1)
    
50.                 V(row,i)=Vmm(row,1);
    
51.             else
    
52.             end
    
53.         end
    
54.     end
    
55.    
    
56.     for k=1:group
    
57.         xx=X(:,j)';
    
58. %         fz1(k)=J_ITAE(xx);
    
59.         fz1(j)=J_ITAE_discrete(xx);
    
60.         if fz1(k)<fz(k)
    
61.             P_p(:,k)=X(:,k);
    
62.             fz(k)=fz1(k);
    
63.         end
    
64.         if fz(k)<P_g
    
65.             P_g=fz(k);
    
66.             
    
67.         end
    
68.     end
    
69.     [P_g1 I]=min(fz);
    
70.     Ess(itrtn)=P_g1;
    
71. %     globe=P_p(:,I);
    
72.     globe(:,itrtn+1)=X(:,I)
    
73. %     if P_g1<e   
    
74. %         break;
    
75. %     end
    
76. end
    
77. kp=globe(1,100);ti=globe(2,100)/globe(1,100);td=globe(3,100)/globe(2,100);T=0.5;
    
78.    globe_0=[1.719,2.813,1.151];
    
79. kp_0=globe_0(1);ti_0=globe_0(2)/globe_0(1);td_0=globe_0(3)/globe_0(2);
    
80. 
    
81. numpid=[td*ti,ti,1];denpid=[0,ti/kp,0];G_pid=tf(numpid,denpid);
    
82.   numpid_0=[td_0*ti_0,ti_0,1];denpid_0=[0,ti_0/kp_0,0];G_pid_0=tf(numpid_0,denpid_0);
    
83. 
    
84. [numz,denz]=pade(T,4);G_exp=tf(numz,denz);
    
85. 
    
86. numd=([0,0,1]);dend=([1,2,1]);G_d=tf(numd,dend);
    
87. 
    
88. 
    
89. G_c=feedback(G_d*G_pid,G_exp);step(G_c,'b');hold on
    
90.     G_c_0=feedback(G_d*G_pid_0,G_exp);step(G_c_0,'r') ;   

复制代码

1. function  q=J_ITAE(x)%(x,ht)
   
2. % axis([0,40,1,1.2]);
   
3. Kp=x(1);Ki=x(2);Kd=x(3);
   
4. Ti=Kp/Ki;Td=Kd/Kp;
   
5. T=0.5
   
6. 
   
7.     numpid=[Kp*Td*Ti,Kp*(Ti+Td),Kp];denpid=[Td*Ti,Ti,0];
   
8.     [numz,denz]=pade(T,4);
   
9.     numd=([0,0,1]);dend=([1,2,1]);
   
10. %     num=conv(conv(numpid,numd),denz);xyj
    
11. %     num=conv(conv(numpid,numd),numz); jsx1
    
12.     num=conv(conv(numpid,numd),denz);%jsx2
    
13.     den1=conv(conv(denpid,dend),denz);
    
14.     den2=conv(conv(numpid,numd),numz);
    
15.     den=den1+den2;
    
16.    
    
17. %     t=0:0.1:50;xyj
    
18.     t=0:0.1:100;
    
19. %     ii=find(t>=T);
    
20. %     [y,x]=step(num,den,t);
    
21. %     y=[zeros(ii(1)-1,1);y((ii(1)+1):length(t))];
    
22. 
    
23. %     y(1:length(t)-ii(1)+1)];
    
24.    
    
25. %     if (ht==1) plot(t,y,'-');
    
26. %     end
    
27. %     if (ht==2) plot(t,y,'--');
    
28. %     end
    
29.     q=0;tt=0;
    
30.     for j=1:501
    
31.         q=q+abs(1-y(j))*tt*0.1;
    
32.         tt=tt+0.1;
    
33.     end
    
34. end

复制代码

风花雪月

原帖由 hunan 于 2007-5-29 21:29 发表

                               
登录/注册后可看大图

你qq多少类 ..加下我吧 39383685

不好意思，我的QQ纯用于闲聊
另，今后matlab问题请到matlab讨论区去问

guanglwp

我做的毕业设计也是粒子群优化PID参数的，我找了一个粒子群的源程序，结合了PID程序，但是不能运行，请大侠帮忙给我改改，万分感谢

程序如下：
size=30;
code=3;
minx(1)=zeros(1);
maxx(1)=20*ones(1);
minx(2)=zeros(1);
maxx(2)= 1.0*ones(1);
minx(3)=zeros(1);
maxx(3)=1.0*ones(1);

   
c1=1.4962;               %学习因子1
c2=1.4962;               %学习因子2
w=0.7298;                %惯性权重
MaxDT=200;               %最大迭代次数
D=3;                     %搜索空间维数（未知数个数）
N=40;                    %初始化群数个体数目      
eps=10;                  %设置精度
%――――初始化种群的个体（可以在这里设置位置和速度的范围）―――

for  i=1:N
    for j=1:D
X(i,j)=randn;     %随机初始化位置
V(i,j)=randn;     %随机初始化速度
end        
end
%―――先计算各个粒子的适应度，并初始化pi和pg―――
for i=1:N
    p(i)=pid_gaf(x(i,:));
    y(i,:)=x(i,:);
                                                                                
end                                                                           
pg=x(1,:);               %pg为全局最优
for  i=2:N
    if  pid_gaf(x(I,:),D)<pid_gaf(pg,D)        
   pg=x(i,:);
    end
end            
%―――进入主要循环，按照公式依次迭代，直到满足精度要求―――

for  t=1:maxDT
for  i=1:N
  v(i,:)=w*v(i,:)+c1*rand(y(i,:)-x(i,:)) +c2*rand*(pg-x(i,: ));
x(i,: )=x(i,:)+v(i,:);   
if   pid_gaf(x(i,:),D)<p(i)
     p(i)=pid_gaf(x(i,:),D);
     y(i,: )=x(i,:);
end

if p(i)<fitness(pg,D)
    pg=y(i,:);
    end
pbest(t)=fitness(pg,D);
end
%最后给出计算结果。
Figure(1);
Xlabel('time');
Ylabel('best j');

Figure(2);
Polt(timef,rin,'r',timef,yout,'b');
Xlabel('Time(s)');
Ylabel('rin,yout');



%子程序
function[kpid, bsj]=pid_gaf(x(i,:))
global rin yout timef

ts=0.001;       %采样时间
sys=tf(400,[1,50,0]);        %传递函数
dsys=c2d(sys,ls,'z'); %将传递函数进行Z变换
[num,den]=tfdata(dsys,'v');


rin=1.0;
u_1=0.0;u_2=0.0;
y_1=0.0;y_2=0.0;
x=[0,0,0];
b=0;
error_1=0;
tu=1;
s=0;
p=100;
  

for k=1:1:p
    timef(k)=k*ts;
    r(k)=rin;
  

u(k)=kpidi(i)*x(1)+kpidi(2)*x(2)+kpid(3)*x(3);


if  u(k)>=10
     u(k)=10;
end

if  u(k)<=-10
     u(k)=-10;
   
end

yout(k)=-den(2)*y_1-den(3)*y_2+num(2)*u_1+num(3)*u_2;
error(k)=r(k)-yout(k);
%-----return of PID parameters--------

u_2=u_1;u_1=u(k);
y_2=y_1;y_1=yout(k);

x(1)=error(k);                       %calculating p
x(2)=(error(k)-error(k)_1)/ts;       %calculating D
x(3)=x(3)+error(k)*ts;               %calculating  I

error_2=error_1;
error_1=error(k);
if  s==0
    if  yout(k)>0.95&yout(k)<1.05
tu=timef(k)
   s=1;
end
end
end
for   i=1:1:P
ji(i)=abs(error(i))
b=b+ji(i);
if  i>1
error(i)=yout(i)-yout(i-1);
if      erry(i)<0
b=b+100*abs(erry(i));
   end
end
end

风花雪月

你给出的这个程序错误无数，根本没办法修改，没有什么参考价值

zzfox2002

可是你给的程序肯定无法运行。。。。