带非线性约束的智能算法求解

Halcom

百度网盘链接：https://pan.baidu.com/s/1gpGRTRP_nflufLwxp8vzVA   提取码：vmiw

具体链接在halcom.cn论坛，联系人QQ：3283892722
该论坛是一个学习交流平台，我会逐一的和大家分享学习。
欢迎大家录制视频，并提交给我，我来设置视频，你可在论坛进行打赏分享。
视频专用播放器：http://halcom.cn/forum.php?mod=viewthread&tid=258&extra=page%3D1

问题：

已知31个点的极坐标（rj，αj）。设有5个待求未知数ck(c1,c2,c3,c4,c5)，根据约束条件求极小值f下的ck。理论上f越接近0越好。同时，顺带也能求出与31个αj对应的31个未知数βj。根据约束条件第一个等式可找到βj与ck存在的关系，因而转化为求以ck为自变量的函数极小值问题。

附录：31个点坐标数据

极坐标形式：

（rj, αj）={{11.055, 0.},{11.055, 0.10472}, {11.055, 0.20944}, {11.055, 0.314159}, {11.055, 0.418879},{11.055, 0.523599}, {11.055, 0.628319}, {11.055, 0.733038}, {11.055, 0.837758},{11.055, 0.942478}, {11.055, 1.0472}, {11.055, 1.15192}, {11.055, 1.25664},{11.055, 1.36136}, {11.055, 1.46608}, {11.055, 1.5708}, {11.1159, 1.67552},{11.302, 1.78024}, {11.6239, 1.88496}, {12.1012, 1.98968}, {11.06, 2.0944},{9.4082, 2.19911}, {8.2645, 2.30383}, {7.4413, 2.40855}, {6.8355, 2.51327},{6.3855, 2.61799}, {6.0533, 2.72271}, {5.8146, 2.82743}, {5.6535, 2.93215},{5.5605, 3.03687}, {5.53, 3.14159}}

单独拿出来：

rj={11.055, 11.055,11.055, 11.055, 11.055, 11.055, 11.055, 11.055, 11.055, 11.055, 11.055, 11.055,11.055, 11.055, 11.055, 11.055, 11.1159, 11.302, 11.6239, 12.1012, 11.06,9.4082, 8.2645, 7.4413, 6.8355, 6.3855, 6.0533, 5.8146, 5.6535, 5.5605, 5.53}

αj={0., 0.10472,0.20944, 0.314159, 0.418879, 0.523599, 0.628319, 0.733038, 0.837758, 0.942478,1.0472, 1.15192, 1.25664, 1.36136, 1.46608, 1.5708, 1.67552, 1.78024, 1.88496,1.98968, 2.0944, 2.19911, 2.30383, 2.40855, 2.51327, 2.61799, 2.72271, 2.82743,2.93215, 3.03687, 3.14159}

fun2适应度函数如下：

1. function [betaj, s] = fun2(ck, popmin, popmax)
   
2. global rj aj nvar
   
3. % 上下限
   
4. % lb = -pi*ones(size(aj));
   
5. lb = aj-0.1;
   
6. lb(1) = 0;  lb(end) = pi;
   
7. % ub = 5*pi*ones(size(aj));
   
8. ub = aj+0.1;
   
9. ub(1) = 0;  ub(end) = pi;
   
10. constraint = 0;
    
11. for i=1:nvar
    
12.     constraint = constraint+i*abs(ck(i));
    
13. end
    
14. if(constraint<1)
    
15.     beta0 = rand( size(aj) );
    
16.     options = optimoptions( 'fmincon', 'Algorithm', 'sqp', 'Display', 'off' );
    
17.     [betaj, s] = fmincon( @(x)fun(x, ck), beta0, [],[],[],[], lb, ub, @(x)nonline_constraint(x, ck), options );
    
18. else
    
19.     betaj = zeros( size(aj) );
    
20.     s = 1e3;
    
21. end
    
22. s = abs(s);

复制代码

fun目标函数如下：

1. function y = fun(betaj, ck)
   
2. global rj aj nvar
   
3. y=0;
   
4. % 误差最小
   
5. for i=1:length(rj)
   
6.     y=y+rj(i)+11.055/(1+sum(ck))*(cos(aj(i)-betaj(i))+ck(1)*cos(aj(i)+1*betaj(i))+...
   
7.     ck(2)*cos(aj(i)+2*betaj(i))+ck(3)*cos(aj(i)+3*betaj(i))+ck(4)*cos(aj(i)+4*betaj(i))+...
   
8.     ck(5)*cos(aj(i)+5*betaj(i)));
   
9. end

复制代码

非线性月数函数如下：

1. function [c, ceq] = nonline_constraint(betaj, ck)
   
2. global rj aj nvar
   
3. c = [ ];      % 不等式约束
   
4. ceq=[ ];      % 等式约束
   
5. for i=1:length(aj)
   
6.     ceq(i) = sin(aj(i)-betaj(i))+ck(1)*cos(aj(i)+1*betaj(i))+...
   
7.     ck(2)*cos(aj(i)+2*betaj(i))+ck(3)*cos(aj(i)+3*betaj(i))+ck(4)*cos(aj(i)+4*betaj(i))+...
   
8.     ck(5)*cos(aj(i)+5*betaj(i));
   
9. end
   

复制代码

主程序如下：

1. clc,clear,close all
   
2. warning off;
   
3. parameters;
   
4. global rj aj nvar
   
5. rj = data(:,1);
   
6. aj = data(:,2);
   
7. lendata = length(rj);   % 点的数量
   
8. tic;
   
9. % PSO 参数
   
10. c1 = 1.4995;  
    
11. c2 = 1.4995;
    
12. Vmin = -1;
    
13. Vmax = 1;
    
14. maxiter = 100;  % 迭代次数
    
15. sizepop = 1000;  % 种群数量
    
16. popmin = -0.1;  popmax = 0.25; % x
    
17. nvar = 5;   % 5个未知量
    
18. %% 初始化种群
    
19. for i=1:sizepop
    
20.     pop(i,:) = popmin + (popmax-popmin).*rand(1,nvar);
    
21.     [betaj(i,:), fitness(i)] = fun2( pop(i,:), popmin, popmax );
    
22.     V(i,1:nvar) = Vmin + (Vmax-Vmin).*rand(1,nvar);
    
23. end
    
24. % 记录一组最优值
    
25. [bestfitness,bestindex]=min(fitness);
    
26. zbest=pop(bestindex,:);   % 全局最佳
    
27. beta_zbest = betaj(bestindex,:);   % 全局最佳
    
28. gbest=pop;                % 个体最佳
    
29. fitnessgbest=fitness;     % 个体最佳适应度值
    
30. fitnesszbest=bestfitness; % 全局最佳适应度值
    
31. wmax = 0.9;  wmin = 0.4;  % 权重
    
32. %% 迭代寻优
    
33. for i=1:maxiter
    
34.     disp(['当前迭代次数：', num2str(i)])
    
35.     for j=1:sizepop
    
36.         % 自适应权重1
    
37. %         w = wmin + (wmax-wmin)*(fitnessgbest(j)-min(fitness))/( max(fitness)-min(fitness) );
    
38.         % 自适应权重2
    
39. %         w = wmin - (wmax-wmin)*(fitnessgbest(j)-min(fitness))/( mean(fitness)-min(fitness) );
    
40.         % 自适应权重3
    
41.         if fitnessgbest(j)<=mean(fitness)
    
42.             w = wmin - (wmax-wmin)*(fitnessgbest(j)-min(fitness))/( mean(fitness)-min(fitness) );
    
43.         else
    
44.             w = wmax;
    
45.         end
    
46.         % 速度更新
    
47.         V(j,:) = w*V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
    
48.         V(j,:) = lb_ub(V(j,:), Vmin, Vmax);
    
49.         
    
50.         % 个体更新
    
51.         pop(j,:) = pop(j,:) + 0.5 * V(j,:);
    
52.         pop(j,:) = lb_ub(pop(j,:), popmin, popmax);
    
53.         
    
54.         % 适应度更新
    
55.         [betaj(j,:), fitness(j)] = fun2(pop(j,:), popmin, popmax);
    
56.         
    
57.         % 比较  个体间比较
    
58.         if fitness(j)<fitnessgbest(j)
    
59.             fitnessgbest(j) = fitness(j);
    
60.             gbest(j,:) = pop(j,:);
    
61.         end
    
62.         if fitness(j)<bestfitness
    
63.             bestfitness = fitness(j);
    
64.             zbest =  pop(j,:);
    
65.             beta_zbest = betaj(j,:);   % 全局最佳
    
66.         end
    
67.         
    
68.     end
    
69.    
    
70.     fitness_iter(i) = bestfitness;
    
71. end
    
72. %% 结果
    
73. toc ;
    
74. times = toc;
    
75. fprintf('\n')
    
76. disp(['计算时间 Time =    ',   num2str(times) ])
    
77. fprintf('\n')
    
78. disp(['最优解   ', num2str(zbest)])
    
79. disp(['最优解对应的适应度函数   ', num2str(bestfitness)])
    
80. disp(['最优解对应的beta   ', num2str(beta_zbest)])
    
81. fprintf('\n')
    
82. 
    
83. figure('color',[1,1,1])
    
84. plot(fitness_iter,'ro-','linewidth',2)
    
85. xlabel('迭代次数');
    
86. ylabel('适应度值');
    
87. 
    
88. %% 结果分析
    
89. constraint = 0;
    
90. for i=1:nvar
    
91.     constraint = constraint+i*abs(zbest(i));
    
92. end
    
93. if constraint<1
    
94.     disp('满足约束条件 sum(k*c(k))<1')
    
95. else
    
96.     disp('满足约束条件 sum(k*c(k))<1')
    
97. end
    
98. fprintf('\n')
    
99. [c, ceq] = nonline_constraint(beta_zbest, zbest);
    
100. disp('等式约束关系 please check ceq .')
     
101. disp(['非线性等式约束关系   ', num2str(ceq)]);
     
102. 
     
103. save results.mat

复制代码