蚁群算法ACO优化的控制系统参数逼近算法

Halcom

蚁群算法ACO优化的控制系统参数逼近算法
参考链接：链接：https://pan.baidu.com/s/1ow0aKl3jcyvZrsFQAzghXA 提取码：ix7u

1. clc,clear,close all
   
2. warning off
   
3. feature jit off
   
4. load('U.mat')
   
5. load('Y.mat')
   
6. % zbest_tf = importfile1('zbest.mat');
   
7. % 对Y进行抽样
   
8. Y1 = Y(1:1000:end);
   
9. % ACO 参数
   
10. maxiter = 100;  % 迭代次数
    
11. sizepop = 30;   % 种群数量
    
12. popmin = [0.1,14,0.1]; % x_min
    
13. popmax = [100,20,100]; % x_max
    
14. Rou = 1.0;      % 信息素增量强度
    
15. P0 = 0.01;      % 转移概率   
    
16. nvar = 3;       % nvar个自变量
    
17. % 初始化种群
    
18. for i=1:sizepop
    
19.     pop(i,:) = popmin + (popmax-popmin).*rand(1,nvar)*0.5;
    
20.     fitness(i) = BJ_Fun4( pop(i,:), Y1 );
    
21. end
    
22. % 记录一组最优值
    
23. [bestfitness,bestindex]=min(fitness);
    
24. zbest=pop(bestindex,:);   % 全局最佳
    
25. gbest=pop;                % 个体最佳
    
26. fitnessgbest=fitness;     % 个体最佳适应度值
    
27. fitnesszbest=bestfitness; % 全局最佳适应度值
    
28. %% 迭代寻优
    
29. for i=1:maxiter
    
30.     disp(['当前迭代次数：', num2str(i)])
    
31.     lamda = 1/i;       % 信息素挥发因子
    
32.     [bestfit,flag] =min(fitness);
    
33.     for j=1:sizepop
    
34.         Pt(j) = (fitness(j)-bestfit)./bestfit;
    
35.     end
    
36.     for j=1:sizepop
    
37.         % 转移概率值
    
38.         if Pt(j)<P0
    
39.             pop(j,:) = pop(j,:) + (2*rand(1, nvar)-1)*lamda/2;
    
40.         else
    
41.             pop(j,:) = pop(j,:) + (popmax-popmin).*(rand(1, nvar));
    
42.         end   
    
43.         
    
44.         % 越界处理
    
45.         for k=1:nvar
    
46.             if pop(j,k)>popmax(k)
    
47.                 pop(j,k)=popmax(k);
    
48.             end
    
49.             if pop(j,k)<popmin(k)
    
50.                 pop(j,k)=popmin(k);
    
51.             end
    
52.         end
    
53.         
    
54.         % 适应度值
    
55.         fitness1 = BJ_Fun4(pop(j,:), Y1);
    
56.         fitness2 = BJ_Fun4(gbest(j,:), Y1);
    
57.         if fitness1<fitness2
    
58.             gbest(j,:) = pop(j,:);
    
59.         else
    
60.             pop(j,:) = gbest(j,:);
    
61.         end
    
62.         
    
63.         fitness(j) = (1-Rou)*fitness(j) + BJ_Fun4(pop(j,:), Y1);
    
64.         
    
65.         if fitness(j) < fitnesszbest
    
66.             fitnesszbest = fitness(j);
    
67.             fitnessgbest(j) = BJ_Fun4(pop(j,:), Y1);
    
68.         end
    
69.     end
    
70.     [fitness_iter(i),index]= min(fitnessgbest);
    
71. end
    
72. zbest = pop(index,:);  % 最优个体
    
73. save zbest.mat zbest
    
74. %% 显示
    
75. disp('最优解')
    
76. disp(['最优解： ', num2str( zbest )])
    
77. fprintf('\n')
    
78. disp(['最优解对应的目标函数值： ', num2str( min(fitness_iter) )])
    
79. disp('传递函数');
    
80. sys = tf( zbest(1), [1, zbest(2), zbest(1)] )
    
81. %% 画图显示
    
82. figure('color',[1,1,1])
    
83. plot(fitness_iter,'ro-','linewidth',2)
    
84. grid on
    
85. xlabel('迭代次数');
    
86. ylabel('适应度函数值');
    
87. %% 响应曲线
    
88. Ygap = 1000;
    
89. [rint_BFO, yout_BFO] = BJ_Fun_response4(zbest, Y(1:Ygap:end));
    
90. figure;hold on;
    
91. plot(rint_BFO.*max(Y(1:Ygap:end)),'r-','linewidth',2);
    
92. plot(yout_BFO.*max(Y(1:Ygap:end)),'b.-','linewidth',2);
    
93. legend('阶跃控制信号','ACO算法优化的信号逼近')
    
94. axis tight;
    
95. hold off
    
96. 
    
97. %% 响应曲线
    
98. Amp= 5000;
    
99. Ygap = 1;
    
100. [rint_BFO2, yout_BFO2] = BJ_Fun_response5( Y(1:Ygap:end), zbest);
     
101. figure;hold on;
     
102. plot(rint_BFO2.*Amp,'k-','linewidth',2);
     
103. plot(Y(1:Ygap:end),'r-','linewidth',2);
     
104. plot(yout_BFO2*Amp,'b.-','linewidth',2);
     
105. legend('阶跃控制信号','实际采集信号','ACO算法优化的信号逼近')
     
106. axis tight;
     
107. hold off

复制代码

适应度函数如下：

1. function fitness = BJ_Fun4(xx, Y1)
   
2. ts=1e-2;
   
3. aa=xx(1); bb=xx(2); cc=xx(3);
   
4. % sys = tf( [1.6], [1,1.5,1.6] );
   
5. sys = tf( aa, [1, bb, aa] );
   
6. % dsys = c2d(sys, ts, 'z');
   
7. % [num,den] = tfdata( dsys, 'v' );
   
8. for k=1:length(Y1)
   
9.     timef(k) = k*ts;
   
10.     r(k) = 1.0;
    
11. %     r(k) = Y1(k)./(max(Y1)+eps);
    
12. %     r(k) = stepfun(timef(k),0);
    
13.     u(k) = Y1(k)./(max(Y1)+eps);
    
14. end
    
15. rout =  lsim(sys, r, timef);
    
16. yout =  lsim(sys, u, timef);
    
17. error1 = r - rout';
    
18. error2 = u - yout';
    
19. fitness1 = sum( abs(error1) ) + 100*max( abs( error1 ) );
    
20. fitness2 = sum( abs(error2) ) + 0*max( abs( error2 ) );
    
21. fitness = fitness1+fitness2;

复制代码

响应度曲线：

1. function [r, yout] = BJ_Fun_response4(xx, Y1)
   
2. ts=1e-0;
   
3. aa=xx(1); bb=xx(2); cc=xx(3);
   
4. % sys = tf( [1.6], [1,1.5,1.6] );
   
5. sys = tf( aa, [1, bb, aa] );
   
6. % dsys = c2d(sys, ts, 'z');
   
7. % [num,den] = tfdata( dsys, 'v' );
   
8. for k=1:length(Y1)
   
9.     timef(k) = k*ts;
   
10. %     r(k) = 1.0;
    
11.     r(k) = Y1(k)./(max(Y1)+eps);
    
12. %     r(k) = stepfun(timef(k),0);
    
13. end
    
14. yout =  lsim(sys, r, timef);
    
15. % error = r - yout';
    
16. % fitness = sum( abs(error) ) + 300*min( abs( error ) );

复制代码