声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 911|回复: 2

[编程技巧] 出错代码:跪求 粒子群优化PID参数的程序 那里出错了? 在线等~!!! 运行不了~!

[复制链接]
发表于 2010-4-7 20:25 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
clear;
x0=[2.813,1.719,1.151];%x=[Kp,Ki,Kd]
c1=1.495;c2=1.495;n=3;group=50;Dmax=100;
% Vmm=[0.4 1.2;0.35 1.05;0.05 0.15];
Vmm=[-2.5 2.5;-2.5 2.5;-2.5 2.5];%可以为负吗

Xmm=[0 5;0 5;0 5];
X=zeros(n,group);V=zeros(n,group);
% omega=zeros(1,group);
% e=0.001;

X(1,:)=x0(1)+5*rand(1,group);V(1,:)=Vmm(1,1)+5*rand(1,group);
X(2,:)=x0(2)+5*rand(1,group);V(2,:)=Vmm(2,1)+5*rand(1,group);
X(3,:)=x0(3)+5*rand(1,group);V(3,:)=Vmm(3,1)+5*rand(1,group);
P_p=X;globe=zeros(n,Dmax);

%%%评价每个粒子适应值,寻找出 P_globle
for j=1:group
    xx=X(:,j)';
%     fz(j)=J_ITAE(xx);
    fz(j)=J_ITAE_discrete(xx);
end
[P_g,I]=min(fz);%P_g  1*1 ?
globe(:,1)=X(:,I);

for itrtn=1:Dmax
    omega=0.739;
    r1=rand(1);r2=rand(1);
    for i=1:group
        V(:,i)=omega*V(:,i)+c1*r1*(P_p(:,i)-X(:,i))+c2*r2*(globe(itrtn)-X(:,i));
    end
    X=X+V;%先虑出位置,如果先是速度虑出,则位置虑出的受速度影响。
    for i=1:group
    %compare each dimension of V
        for row=1:n
            if X(row,i)>Xmm(row,2)
                X(row,i)=Xmm(row,2);
            elseif X(row,i)<Xmm(row,1)
                X(row,i)=Xmm(row,1);
            else
            end
        end
    end
    for i=1:group
        %compare each dimension of V
        for row=1:n
            if V(row,i)>Vmm(row,2)
                V(row,i)=Vmm(row,2);
            elseif V(row,i)<Vmm(row,1)
                V(row,i)=Vmm(row,1);
            else
            end
        end
    end
   
    for k=1:group
        xx=X(:,j)';
%         fz1(k)=J_ITAE(xx);
        fz1(j)=J_ITAE_discrete(xx);
        if fz1(k)<fz(k)
            P_p(:,k)=X(:,k);
            fz(k)=fz1(k);
        end
        if fz(k)<P_g
            P_g=fz(k);
            
        end
    end
    [P_g1 I]=min(fz);
    Ess(itrtn)=P_g1;
%     globe=P_p(:,I);
    globe(:,itrtn+1)=X(:,I)
%     if P_g1<e   
%         break;
%     end
end
kp=globe(1,100);ti=globe(2,100)/globe(1,100);td=globe(3,100)/globe(2,100);T=0.5;
   globe_0=[1.719,2.813,1.151];
kp_0=globe_0(1);ti_0=globe_0(2)/globe_0(1);td_0=globe_0(3)/globe_0(2);

numpid=[td*ti,ti,1];denpid=[0,ti/kp,0];G_pid=tf(numpid,denpid);
  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);

[numz,denz]=pade(T,4);G_exp=tf(numz,denz);

numd=([0,0,1]);dend=([1,2,1]);G_d=tf(numd,dend);


G_c=feedback(G_d*G_pid,G_exp);step(G_c,'b');hold on
    G_c_0=feedback(G_d*G_pid_0,G_exp);step(G_c_0,'r') ;   
----------------------------------------------------------------------------------------------------------------------------




function  q=J_ITAE(x)%(x,ht)
% axis([0,40,1,1.2]);
Kp=x(1);Ki=x(2);Kd=x(3);
Ti=Kp/Ki;Td=Kd/Kp;
T=0.5

    numpid=[Kp*Td*Ti,Kp*(Ti+Td),Kp];denpid=[Td*Ti,Ti,0];
    [numz,denz]=pade(T,4);
    numd=([0,0,1]);dend=([1,2,1]);
%     num=conv(conv(numpid,numd),denz);xyj
%     num=conv(conv(numpid,numd),numz); jsx1
    num=conv(conv(numpid,numd),denz);%jsx2
    den1=conv(conv(denpid,dend),denz);
    den2=conv(conv(numpid,numd),numz);
    den=den1+den2;
   
%     t=0:0.1:50;xyj
    t=0:0.1:100;
%     ii=find(t>=T);
%     [y,x]=step(num,den,t);
%     y=[zeros(ii(1)-1,1);y((ii(1)+1):length(t))];

%     y(1:length(t)-ii(1)+1)];
   
%     if (ht==1) plot(t,y,'-');
%     end
%     if (ht==2) plot(t,y,'--');
%     end
    q=0;tt=0;
    for j=1:501
        q=q+abs(1-y(j))*tt*0.1;
        tt=tt+0.1;
    end
end

[ 本帖最后由 jerry_wn 于 2010-4-7 21:19 编辑 ]
回复
分享到:

使用道具 举报

发表于 2010-4-7 21:11 | 显示全部楼层
怎今天大家都忘记了本版规则
6) 求助完整格式:出错代码和出错提示
发表于 2010-4-8 10:27 | 显示全部楼层
什么错误?                       .
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-13 19:12 , Processed in 0.072494 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表