【控制】粒子群算法优化SO 调谐 PI 控制器，用于可变惯量 BLDC 电机的速度控制附matlab代码.zip

% 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;