|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
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 编辑 ] |
|