带非线性约束的智能算法求解
百度网盘链接：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适应度函数如下：
function = fun2(ck, popmin, popmax)
global rj aj nvar
% 上下限
% lb = -pi*ones(size(aj));
lb = aj-0.1;
lb(1) = 0;lb(end) = pi;
% ub = 5*pi*ones(size(aj));
ub = aj+0.1;
ub(1) = 0;ub(end) = pi;
constraint = 0;
for i=1:nvar
constraint = constraint+i*abs(ck(i));
end
if(constraint<1)
beta0 = rand( size(aj) );
options = optimoptions( 'fmincon', 'Algorithm', 'sqp', 'Display', 'off' );
= fmincon( @(x)fun(x, ck), beta0, [],[],[],[], lb, ub, @(x)nonline_constraint(x, ck), options );
else
betaj = zeros( size(aj) );
s = 1e3;
end
s = abs(s);fun目标函数如下：
function y = fun(betaj, ck)
global rj aj nvar
y=0;
% 误差最小
for i=1:length(rj)
y=y+rj(i)+11.055/(1+sum(ck))*(cos(aj(i)-betaj(i))+ck(1)*cos(aj(i)+1*betaj(i))+...
ck(2)*cos(aj(i)+2*betaj(i))+ck(3)*cos(aj(i)+3*betaj(i))+ck(4)*cos(aj(i)+4*betaj(i))+...
ck(5)*cos(aj(i)+5*betaj(i)));
end非线性月数函数如下：
function = nonline_constraint(betaj, ck)
global rj aj nvar
c = [ ]; % 不等式约束
ceq=[ ]; % 等式约束
for i=1:length(aj)
ceq(i) = sin(aj(i)-betaj(i))+ck(1)*cos(aj(i)+1*betaj(i))+...
ck(2)*cos(aj(i)+2*betaj(i))+ck(3)*cos(aj(i)+3*betaj(i))+ck(4)*cos(aj(i)+4*betaj(i))+...
ck(5)*cos(aj(i)+5*betaj(i));
end
主程序如下：clc,clear,close all
warning off;
parameters;
global rj aj nvar
rj = data(:,1);
aj = data(:,2);
lendata = length(rj); % 点的数量
tic;
% PSO 参数
c1 = 1.4995;
c2 = 1.4995;
Vmin = -1;
Vmax = 1;
maxiter = 100;% 迭代次数
sizepop = 1000;% 种群数量
popmin = -0.1;popmax = 0.25; % x
nvar = 5; % 5个未知量
%% 初始化种群
for i=1:sizepop
pop(i,:) = popmin + (popmax-popmin).*rand(1,nvar);
= fun2( pop(i,:), popmin, popmax );
V(i,1:nvar) = Vmin + (Vmax-Vmin).*rand(1,nvar);
end
% 记录一组最优值
=min(fitness);
zbest=pop(bestindex,:); % 全局最佳
beta_zbest = betaj(bestindex,:); % 全局最佳
gbest=pop; % 个体最佳
fitnessgbest=fitness; % 个体最佳适应度值
fitnesszbest=bestfitness; % 全局最佳适应度值
wmax = 0.9;wmin = 0.4;% 权重
%% 迭代寻优
for i=1:maxiter
disp(['当前迭代次数：', num2str(i)])
for j=1:sizepop
% 自适应权重1
% w = wmin + (wmax-wmin)*(fitnessgbest(j)-min(fitness))/( max(fitness)-min(fitness) );
% 自适应权重2
% w = wmin - (wmax-wmin)*(fitnessgbest(j)-min(fitness))/( mean(fitness)-min(fitness) );
% 自适应权重3
if fitnessgbest(j)<=mean(fitness)
w = wmin - (wmax-wmin)*(fitnessgbest(j)-min(fitness))/( mean(fitness)-min(fitness) );
else
w = wmax;
end
% 速度更新
V(j,:) = w*V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
V(j,:) = lb_ub(V(j,:), Vmin, Vmax);
% 个体更新
pop(j,:) = pop(j,:) + 0.5 * V(j,:);
pop(j,:) = lb_ub(pop(j,:), popmin, popmax);
% 适应度更新
= fun2(pop(j,:), popmin, popmax);
% 比较个体间比较
if fitness(j)<fitnessgbest(j)
fitnessgbest(j) = fitness(j);
gbest(j,:) = pop(j,:);
end
if fitness(j)<bestfitness
bestfitness = fitness(j);
zbest =pop(j,:);
beta_zbest = betaj(j,:); % 全局最佳
end
end
fitness_iter(i) = bestfitness;
end
%% 结果
toc ;
times = toc;
fprintf('\n')
disp(['计算时间 Time = ', num2str(times) ])
fprintf('\n')
disp(['最优解 ', num2str(zbest)])
disp(['最优解对应的适应度函数 ', num2str(bestfitness)])
disp(['最优解对应的beta ', num2str(beta_zbest)])
fprintf('\n')
figure('color',)
plot(fitness_iter,'ro-','linewidth',2)
xlabel('迭代次数');
ylabel('适应度值');
%% 结果分析
constraint = 0;
for i=1:nvar
constraint = constraint+i*abs(zbest(i));
end
if constraint<1
disp('满足约束条件 sum(k*c(k))<1')
else
disp('满足约束条件 sum(k*c(k))<1')
end
fprintf('\n')
= nonline_constraint(beta_zbest, zbest);
disp('等式约束关系 please check ceq .')
disp(['非线性等式约束关系 ', num2str(ceq)]);
save results.mat
页:
[1]