请选择 进入手机版 | 继续访问电脑版

Hello Mat

 找回密码
 立即注册

QQ登录

只需一步,快速开始

查看: 658|回复: 0

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

[复制链接]

691

主题

813

帖子

2万

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
29493
发表于 2019-3-4 23:51:18 | 显示全部楼层 |阅读模式
百度网盘链接: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。根据约束条件第一个等式可找到βjck存在的关系因而转化为求以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. figure('color',[1,1,1])
  83. plot(fitness_iter,'ro-','linewidth',2)
  84. xlabel('迭代次数');
  85. ylabel('适应度值');

  86. %% 结果分析
  87. constraint = 0;
  88. for i=1:nvar
  89.     constraint = constraint+i*abs(zbest(i));
  90. end
  91. if constraint<1
  92.     disp('满足约束条件 sum(k*c(k))<1')
  93. else
  94.     disp('满足约束条件 sum(k*c(k))<1')
  95. end
  96. fprintf('\n')
  97. [c, ceq] = nonline_constraint(beta_zbest, zbest);
  98. disp('等式约束关系 please check ceq .')
  99. disp(['非线性等式约束关系   ', num2str(ceq)]);

  100. save results.mat
复制代码



本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
程序:算法QQ  3283892722
群智能算法视频:X元,链接http://halcom.cn/forum.php?mod=forumdisplay&fid=73
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则


Python|Opencv|MATLAB|Halcom.cn  

GMT+8, 2019-6-17 09:16 , Processed in 0.101484 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表