✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab仿真内容点击👇
⛄ 内容介绍
针对工业对永磁同步电机调速系统的更高调速精度,更快响应速度这些要求,该文提出了一种新的永磁同步电机控制策略,即利用粒子群算法对模糊PI控制器(Fuzzy PI)的2个参数因子kp,ki进行全局优化,充分发挥了粒子群算法的快速性.利用Matlab工具进行仿真验证,观察控制系统的一阶动态响应.结果表明,系统具有很强的鲁棒性,能够很好地跟踪负载变化,动态响应快,速度跟随准确.
⛄ 部分代码
% PSO Algorithm Radoslaw Wierzbowski
%%
%Motor parameters initialization
clear;
clc;
Pn=4e3; %Nominal power
Un=280; %Nominal voltage
In=17.5; %Nominal current
wn=1980*pi/30; %Angular velocity
Rn=1.69; %Armatures resistance
Ln=6.6e-3; %Armature inductance
J=0.02*2; %Inertia (double catalogue value)
Tn=19.3; %Nominal torque
%Fluxes calculated from diffrent equations basing on nominal values:
psi_E=(Un-Rn*In)/wn; %from armature equation (Un=Rn*In+wn*psi_E)
psi_M=Tn/In; %from torque equation (Tn=psi_M*In)
f_sw=4000; %Converters switching frequency
T_sw=1/f_sw;
Tdelay=0.5*T_sw; %Approximated delay
V_dc=1.2*Un;
Kconv=V_dc; %Converters amp
c=0.1*Tn/wn; %friction coefficient
ki=1/In; %Measuring amp - current
kw=1/wn; %Measuring amp - velocity
ku=1/Un; %Measuring amp - voltage
Tsim=1e-4; %Simulation max step size
%Calculation of parameters for PI current controller (modulus optimum)
Ks=Kconv/Rn/In;
T1=Ln/Rn; %Dominant time constant
Tsum=Tdelay; %Sum of small time constants
Kri=1/(2*Ks*Tsum);
Tr=T1;
%Upper boundaries calculation (symmetrical optimum)
Jprop=10*J;
Tr2=8*Tsum;
Kr2=Jprop/(8*psi_M*kw*(1/ki)*(2*Tsum)^2);
Kpmax=Tr2*Kr2;
Kimax=Kr2;
%%
%Swarm initialization:
n = 20; %Number of particles
dim =2; %Dimensions
steps=120; %Number of steps
c1=2.05; %Cognitive factor
c2=2.05; %Exploration factor
fi=c1+c2;
X=2/abs((2-fi-((fi^2)-(4*fi))^(1/2))); %Constriction factor
rng('shuffle'); %Seed initialization
Evap=1.5; %Evaporation factor
velocity_clamping = 10; %Velocity clamping
diversity_limit = 1; %Diversity limit
psi1=rand(2,1);
psi2=rand(2,1); %Initial random coefficients
swarmdir(1)=1;
swarmdir(2)=1; %Initial directions
position=(50*rand(dim,n)); %Position initialization
speed=10*rand(dim,n); %Velocity initialization
change1=35; %Step of first inertia change
change2=80; %Step of second inertia change
%%
%Initial values
Pbest=position; %Matrix of personal best positions
ObjFun=0*rand(n,1); %Matrix of objective function values
step=1;
%Calculation of initial Kp and Ki and ObjFun values:
for i=1:n
Kp=Pbest(1,i);
Ki=Pbest(2,i);
ObjFun(i)=optimize(Kp,Ki); %Optimize functions is in another file
end
bestlocalObjFun=ObjFun; %Matrix of best local objective function values
[globalObjFun,g]=min(bestlocalObjFun); %Best global objective function value
Gbest = Pbest(:,g) ; %Best global position
speed(:,i)=X*(speed(:,i)+(fi/2)*swarmdir'.*(psi1.*(Pbest(:,i)-position(:,i)))+(fi/2)*swarmdir'.*(psi2.*(Gbest-position(:,i)))); %tablica predkosci
position(:,i)=position(:,i)+speed(:,i);
%
%%
%Initialization done
%Main loop:
for step=1:steps
fprintf('Step %i \n', step);
fprintf('Found optimum: \n Kp= %0.3f \n Ki= %0.3f \n', Gbest(1,1),Gbest(2,1));
%Inertia change conditions
if step==change1
J=J*5;
end
if step==change2
J=J*0.2;
end
%Asynchronous update rule
for i=1:n
%Kp and Ki assignment and simulation
Kp=position(1,i);
Ki=position(2,i);
ObjFun(i)=optimize(Kp,Ki);
%Evaporation conditions
if 0.01+ObjFun(i)>=bestlocalObjFun(i)*Evap
bestlocalObjFun(i)=Evap*bestlocalObjFun(i);
else
bestlocalObjFun(i)=ObjFun(i)+0.01;
Pbest(:,i)=position(:,i);
end
%Gbest calculation
[globalObjFun,g]=min(bestlocalObjFun);
Gbest=Pbest(:,g);
%Swarm diversity calculation
swarm_div(1) = 0.5*(max(position(1,:))-min(position(1,:)));
swarm_div(2) = 0.5*(max(position(2,:))-min(position(2,:)));
%Direction change in case of too high diversity
for d=1:2,
if swarm_div(d) < diversity_limit,
swarmdir(d) = -1;
else
swarmdir(d) = 1;
end
end
%Position and speed calculation
psi1=rand(2,1);
psi2=rand(2,1);
speed(:,i)=X*(speed(:,i)+(fi/2)*swarmdir'.*(psi1.*(Pbest(:,i)-position(:,i)))+(fi/2)*swarmdir'.*(psi2.*(Gbest-position(:,i))));
speed(1,i)= min(max(-velocity_clamping,speed(1,i)),velocity_clamping);
speed(2,i)= min(max(-velocity_clamping,speed(2,i)),velocity_clamping);
position(:,i)=position(:,i)+speed(:,i);
%Walls on Kp=0, Ki=0, Kpmax and Kimax
if position(1,i)<0
position(1,i) = abs(position(1,i));
speed(1,i)=abs(speed(1,i));
end
if position(2,i)<0
position(2,i) = abs(position(2,i));
speed(2,i) = abs(speed(2,i));
end
if position(1,i)>Kpmax
position(1,i) = position(1,i)-(2*(position(1,i)-Kpmax));
speed(1,i)=-speed(1,i);
end
if position(2,i)>Kimax
position(2,i) = position(2,i)-(2*(position(2,i)-Kimax));
speed(2,i)=-speed(2,i);
end
%History logging for plotting purposes
History(1,step,:)=position(1,:);
History(2,step,:)=position(2,:);
speedHistory(1,step,:)=speed(1,:);
speedHistory(2,step,:)=speed(2,:);
GbestHistory(step,:)=Gbest;
Gbestvalue(step,:)=globalObjFun;
%figure(1);
%for i=1:n
% subplot(2,1,1);
% scatter(log10(position(2,i)),log10(position(1,i)));
% title('Pozycje czastek');
% xlabel('log10(Ki)'); ylabel('log10(Kp)');
% axis([0 3 0 3]);
% hold on
%end
% hold off;
end
end
Kp=Gbest(1,1);
Ki=Gbest(2,1);
fprintf('Optimum found: \n Kp= %0.3f \n Ki= %0.3f \n', Gbest(1,1),Gbest(2,1));
figure();
subplot(3,1,1);
scatter(1:steps, GbestHistory(:,1));
hold;
plot([change1-1, change1-1], [1, max(GbestHistory(:,1))], 'LineWidth', 1, 'Color', 'blue');
plot([change2-1, change2-1], [1, max(GbestHistory(:,1))], 'LineWidth', 1, 'Color', 'blue');
xlabel('iteration');
ylabel('Kp');
grid;
subplot(3,1,2);
scatter(1:steps, GbestHistory(:,2));
hold;
plot([change1-1, change1-1], [1, max(GbestHistory(:,2))], 'LineWidth', 1, 'Color', 'blue');
plot([change2-1, change2-1], [1, max(GbestHistory(:,2))], 'LineWidth', 1, 'Color', 'blue');
xlabel('iteration');
ylabel('Ki');
grid;
subplot(3,1,3);
scatter(1:steps, Gbestvalue);
hold;
plot([change1-1, change1-1], [1, max(Gbestvalue)], 'LineWidth', 1, 'Color', 'blue');
plot([change2-1, change2-1], [1, max(Gbestvalue)], 'LineWidth', 1, 'Color', 'blue');
xlabel('iteration');
ylabel('Gbest Value');
grid;
⛄ 运行结果
⛄ 参考文献
[1]陈旭, 宋正强. 基于粒子群优化的永磁同步电机PI控制[J]. 中国高新技术企业, 2008(22):2.
