联系我们contact us
400-123-4657全国服务热线:
- 地址:
- 广东省广州市天河区88号
- 邮箱:
- admin@youweb.com
- 手机:
- 13800000000
- 电话:
- 400-123-4657
【5. 智能优化算法案例】Matlab案例代码解析:智能优化算法篇添加时间:2024-07-22
记录于 2021-10-07。
- 差分进化算法(DE)
- 离散差分进化算法(DE)
- 差分进化算法实现指数拟合(DE-ExpFit)
- 梯度下降法曲线拟合(GradientDescent)
- 灰狼算法(GWO)
- 灰狼-布谷鸟算法(GWOCS)
- 粒子群算法(PSO)
- 离散粒子群算法(PSO)
- 粒子群算法实现指数拟合(PSO-ExpFit)
- 随机游走(RandWalk)
- 模拟退火算法实现指数拟合(SA)
- 模拟退火算法求解旅行商问题(SA-TSP)
- 樽海鞘算法(SSA)
- 状态转移算法(STA)
- 改进状态转移算法(STA)
- 鲸鱼算法(WOA)
基础篇:
https://zhuanlan.zhihu.com/p/408388753Demo 篇:
https://zhuanlan.zhihu.com/p/409300387基本流程:
初始化:
变异:
交叉:
选择:试验值与当前目标值比较,实验值小则更新对应的种群值
边界处理:超过边界的可以用范围内随机值代替或者进行边界吸收直接用边界值代替
改进:
自适应差分进化算法:
%%%%%%%%%%%%% 初始化 %%%%%%%%%%%%%
clear; clc;
% 种群数量
NP = 20;
% 变量维度
D = 2;
% 最大进化代数
G = 200;
% 初始变异算子
F0 = 0.5;
% 交叉算子
CR = 0.1;
% 上限
Xs = 10;
% 下限
Xx = -10;
% 阈值
yz = 1e-8;
%%%%%%%%%%%%% 赋初值 %%%%%%%%%%%%%
Ob = zeros(NP, 1);
Ob1 = zeros(NP, 1);
% 变异种群
v = zeros(D, NP);
% 选择种群
u = zeros(D, NP);
% 初始种群
x = rand(D, NP) * (Xs - Xx) + Xx;
%%%%%%%%%%%%% 计算目标函数 %%%%%%%%%%%%%
for m = 1 : NP
Ob(m) = ObjFun(x(:, m));
end
%%%%%%%%%%%%% 差分进化循环 %%%%%%%%%%%%%
trace = zeros(G, 1);
trace(1) = min(Ob);
for gen = 1 : G
%%%%%%%%%%%%% 变异操作 %%%%%%%%%%%%%
%%%%%%%%%%%%% 自适应变异算子 %%%%%%%%%
lamda = exp(1 - G / (G + 1 - gen));
F = F0 * 2^lamda;
%%%%%%%%%% r1、r2、r3和m互不相同 %%%%%%%%%%
for m = 1 : NP
r1 = randi([1, NP], 1);
while r1 == m
r1 = randi([1, NP], 1);
end
r2 = randi([1, NP], 1);
while r2 == m || r2 == r1
r2 = randi([1, NP], 1);
end
r3 = randi([1, NP], 1);
while r3 == m || r3 == r2 || r3 == r1
r3 = randi([1, NP], 1);
end
v(:, m) = x(:, r1) + F*(x(:, r2) - x(:, r3));
end
%%%%%%%%%%%%% 交叉操作 %%%%%%%%%%%%%
r = randi([1, D], 1);
for n = 1 : D
cr = rand;
if cr <= CR || n == r
u(n, :) = v(n, :);
else
u(n, :) = x(n, :);
end
end
%%%%%%%%%%%%% 边界条件的处理 %%%%%%%%%%%%%
for n = 1:D
for m = 1:NP
if u(n, m) < Xx || u(n, m) > Xs
u(n, m) = rand*(Xs - Xx) + Xx;
end
end
end
%%%%%%%%%%%%% 边界吸收 %%%%%%%%%%%%%
% for n=1 : D
% for m=1 : NP
% if u(n, m) < Xx
% u(n, m)=Xx;
% end
% if u(n, m) > Xs
% u(n, m)=Xs;
% end
% end
% end
%%%%%%%%%%%%% 选择操作 %%%%%%%%%%%%%
for m = 1 : NP
Ob1(m) = ObjFun(u(:, m));
end
for m = 1 : NP
if Ob1(m) < Ob(m)
x(:, m) = u(:, m);
end
end
for m = 1:NP
Ob(m) = ObjFun(x(:, m));
end
trace(gen + 1) = min(Ob);
if min(Ob) < yz
break;
end
end
[SortOb, Index] = sort(Ob);
x = x(:, Index);
% 最优变量
X = x(:, 1);
% 最优值
Y = min(Ob);
disp(['最优x:' num2str(X')]);
disp(['最优y:' num2str(Y)]);
subplot(1, 2, 1)
plot(trace, 'LineWidth', 2);
xlabel('迭代次数');
ylabel('适应度值');
title('适应度进化曲线');
len = 50;
xRange = linspace(Xx, Xs, len);
yRange = linspace(Xx, Xs, len);
[xMap, yMap] = meshgrid(xRange, yRange);
zMap = zeros(len);
for i = 1 : len
for j = 1 : len
zMap(i, j) = ObjFun([xMap(i, j), yMap(i, j)]);
end
end
subplot(1, 2, 2)
surf(xRange, yRange, zMap);
colorbar;
view(-45, -45);
shading interp
hold on
plot3(X(1), X(2), Y, 'o', 'MarkerFaceColor', 'r', 'MarkerSize', 10);
hold off
set(gcf, 'Position', [400, 200, 900, 400]);
function result = ObjFun(x)
x1 = x(1);
x2 = x(2);
t = x1^2 + x2^2;
result = 0.5 + (sin(sqrt(t))^2 - 0.5) / (1 + 0.001 * t)^2;
end
针对整数变量,用floor函数向下取整:
%%%%%%%%%%%%% 离散差分进化算法求函数极值 %%%%%%%%%%%%%
%%%%%%%%%%%%% 初始化 %%%%%%%%%%%%%
clear;clc;
% 种群数量
NP = 20;
% 变量维度
D = 2;
% 最大进化代数
G = 200;
% 初始编译算子
F = 0.5;
% 交叉算子
CR = 0.1;
% 上限
Xs = 100;
% 下限
Xx = -100;
%%%%%%%%%%%%% 赋初值 %%%%%%%%%%%%%
Ob = zeros(NP, 1);
Ob1 = zeros(NP, 1);
trace = zeros(G + 1, 1);
% 变异种群
v = zeros(D, NP);
% 选择种群
u = zeros(D, NP);
% 初始种群
x = randi([Xx, Xs], D, NP);
%%%%%%%%%%%%% 计算目标函数 %%%%%%%%%%%%%
for m = 1:NP
Ob(m) = ObjFunInt(x(:, m));
end
trace(1) = min(Ob);
%%%%%%%%%%%%% 差分进化循环 %%%%%%%%%%%%%
for gen = 1:G
%%%%%%%%%%%%% 变异操作 %%%%%%%%%%%%%
%%%%%%%%%% r1、r2、r3和m互不相同 %%%%%%%%%%
for m = 1:NP
r1 = randi([1, NP], 1);
while r1 == m
r1 = randi([1, NP], 1);
end
r2 = randi([1, NP], 1);
while r2 == m || r2 == r1
r2 = randi([1, NP], 1);
end
r3 = randi([1, NP], 1);
while r3 == m || r3 == r2 || r3 == r1
r3 = randi([1, NP], 1);
end
v(:, m) = floor(x(:, r1) + F*(x(:, r2) - x(:, r3)));
end
%%%%%%%%%%%%% 交叉操作 %%%%%%%%%%%%%
r = randi([1, D], 1);
for n = 1:D
cr = rand;
if cr <= CR || n == r
u(n, :) = v(n, :);
else
u(n, :) = x(n, :);
end
end
%%%%%%%%%%%%% 边界条件的处理 %%%%%%%%%%%%%
%%%%%%%%%%%%% 边界吸收 %%%%%%%%%%%%%
for n = 1:D
for m = 1:NP
if u(n, m) < Xx
u(n, m) = Xx;
end
if u(n, m) > Xs
u(n, m) = Xs;
end
end
end
%%%%%%%%%%%%% 选择操作 %%%%%%%%%%%%%
for m = 1:NP
Ob1(m) = ObjFunInt(u(:, m));
end
for m = 1:NP
if Ob1(m) > Ob(m)
x(:, m) = u(:, m);
end
end
for m = 1:NP
Ob(m) = ObjFunInt(x(:, m));
end
trace(gen + 1) = min(Ob);
end
[SortOb, Index] = sort(Ob);
x = x(:, Index);
% 最优变量
X = x(:, 1);
% 最优值
Y = min(Ob);
disp(['最优x:' num2str(X')]);
disp(['最优y:' num2str(Y)]);
subplot(1, 2, 1)
plot(trace, 'LineWidth', 2);
xlabel('迭代次数');
ylabel('适应度值');
title('适应度进化曲线');
len = 200;
xRange = linspace(Xx, Xs, len);
yRange = linspace(Xx, Xs, len);
[xMap, yMap] = meshgrid(xRange, yRange);
zMap = zeros(len);
for i = 1 : len
for j = 1 : len
zMap(i, j) = ObjFunInt([xMap(i, j), yMap(i, j)]);
end
end
subplot(1, 2, 2)
surf(xRange, yRange, zMap);
colorbar;
view(-45, 45);
shading interp
hold on
plot3(X(1), X(2), Y, 'o', 'MarkerFaceColor', 'r', 'MarkerSize', 10);
hold off
set(gcf, 'Position', [400, 200, 900, 400]);
function result = ObjFunInt(x)
result = -((x(1).^2 + x(2) - 1).^2 + (x(1) + x(2).^2 - 7).^2)/200 + 10;
end
clear;clc;
t = 0.2*(1:3000)';
data = 400*exp(-t/5) + 10*randn(size(t));
tic;
p = DE_ExpFit(t, data)
toc;
fit = p(1)*exp(-t/p(2)) + p(3);
plot(t, data, t, fit, 'LineWidth', 2)
legend('model', 'DE');
title('差分进化算法实现指数拟合');
function [x_opt, y_opt]= DE_ExpFit(t, Et)
%{
函数功能:差分进化算法实现指数拟合:y=a*exp(-x/b) + c
输入:
t:自变量;
Et:因变量;
输出:
x_opt:最优解;
y_opt:适应度(目标函数值);
%}
% 初始值
NP = 50; % 种群数量
D = 3; % 变量维度
G = 200; % 最大进化代数
F0 = 0.5; % 初始变异算子
CR = 0.9; % 交叉算子
max_sig = max(Et);
Xmin = [0.5*max_sig, -2, 0]; % 下限
Xmax = [1.5*max_sig, 4, 0.5*max_sig]; % 上限;
%%%%%%%%%%%%% 赋初值 %%%%%%%%%%%%%
Ob = zeros(NP, 1);
Ob1 = zeros(NP, 1);
v = zeros(NP, D); % 变异种群
u = zeros(NP, D); % 选择种群
x = rand(NP, D).*repmat(Xmax - Xmin, NP, 1) + repmat(Xmin, NP, 1);
%%%%%%%%%%%%% 计算目标函数 %%%%%%%%%%%%%
for m = 1 : NP
Ob(m) = Obj_Fit(x(m, :), t, Et);
end
y_opt = zeros(1, G);
%%%%%%%%%%%%% 差分进化循环 %%%%%%%%%%%%%
for gen = 1 : G
%%%%%%%%%%%%% 变异操作 %%%%%%%%%%%%%
%%%%%%%%%%%%% 自适应变异算子 %%%%%%%%%
lamda = exp(1 - G/(G + 1 - gen));
F = F0*2^lamda;
%%%%%%%%%% r1、r2、r3和m互不相同 %%%%%%%%%%
for m = 1 : NP
r1 = randi([1, NP], 1);
while r1 == m
r1 = randi([1, NP], 1);
end
r2 = randi([1, NP], 1);
while r2 == m || r2 == r1
r2 = randi([1, NP], 1);
end
r3 = randi([1, NP], 1);
while r3 == m || r3 == r2 || r3 == r1
r3 = randi([1, NP], 1);
end
v(m, :) = x(r1, :) + F*(x(r2, :) - x(r3, :));
end
%%%%%%%%%%%%% 交叉操作 %%%%%%%%%%%%%
r = randi([1, D], 1);
for n = 1:D
cr = rand;
if cr <= CR || n == r
u(:, n) = v(:, n);
else
u(:, n) = x(:, n);
end
end
%%%%%%%%%%%%% 边界条件的处理 %%%%%%%%%%%%%
% ①范围内随机值
for n = 1:NP
for m = 1:D
if u(n, m) < Xmin(m) || u(n, m) > Xmax(m)
u(n, m) = rand*(Xmax(m) - Xmin(m)) + Xmin(m);
end
end
end
% ②边界吸收
% for n=1 : NP
% for m=1 : D
% if u(n, m) < Xmin(m)
% u(n, m)=Xmin(m);
% end
% if u(n, m) > Xmax(m)
% u(n, m)=Xmax(m);
% end
% end
% end
%%%%%%%%%%%%% 选择操作 %%%%%%%%%%%%%
for m = 1 : NP
Ob1(m) = Obj_Fit(u(m, :), t, Et);
end
for m = 1 : NP
if Ob1(m) < Ob(m)
x(m, :) = u(m, :);
end
end
for m = 1 : NP
Ob(m) = Obj_Fit(x(m, :), t, Et);
end
y_opt(gen) = min(Ob);
end
[~, Index] = min(Ob);
x_opt = x(Index, :);
x_opt(2) = 10^(x_opt(2));
end
% 目标函数
function fitError = Obj_Fit(p0, t ,Et)
A = p0(1);
B = p0(2);
C = p0(3);
f = A*exp(-t/10^B) + C;
fitError = norm(Et - f);
end
clear; clc;
t = [ 0.1350, 0.2600, 0.3800, 0.5050, 0.6450, 0.8100, 1.0000, 1.2150, 1.4600, 1.7400, 2.0550, 2.4150, ...
2.8150, 3.2650, 3.7650, 4.3200, 4.9350, 5.6100, 6.3450, 7.1500, 8.0200, 8.9650, 9.9800, 11.0650, ...
12.2200, 13.4500, 14.7500, 16.1150, 17.5450, 19.0400, 20.5950, 22.1950, 23.8450, 25.5350, ...
27.2600, 29.2650, 32.2700, 36.2800, 41.4850, 47.8900, 55.8950, 65.5000, 77.5050, 92.3100, ...
108.7150, 128.3250, 151.9300, 179.9350, 212.9400, 252.9450]';
M = length(t);
% 构造模型
A =1800;
B = 15;
qMrl = 5;
AA = 50;
C = 10;
rng('default');
ratio = 0;
noise = ratio * randn(M, 1);
Et = A * exp(-t / B - 0.5 * qMrl * t.^2) + AA * exp(-t / B) + C + noise;
% 梯度下降法学习
% 初始参数
alpha = 1; % 学习率大收敛快,可能有震荡
iteMax = 1000000;
iniA = 1;
iniB = rand;
iniQMrl = rand;
iniAA = rand;
iniC = rand;
initW = [iniA; iniB; iniQMrl; iniAA; iniC];
err = 1e-6;
J = zeros(iteMax, 1);
G = zeros(size(initW));
e = 0.1;
convIter = 0;
for i = 1 : iteMax
gradA = 2 * iniA * exp(- 0.5 * iniQMrl^2 * t.^2 - t / iniB^2);
gradB = (2 * iniA^2 * t .* exp(-0.5 * iniQMrl^2 * t.^2 - t / iniB^2)) / iniB^3 + (2 * iniAA^2 * t .* exp(-t /iniB^2)) / iniB^3;
gradQMrl = - iniA^2 * iniQMrl * t.^2 .* exp(-0.5 * iniQMrl^2 * t.^2 - t / iniB^2);
gradAA = 2 * iniAA * exp(-t / iniB^2);
gradC = 2 * iniC * ones(M, 1);
grad = 1 / M * [gradA, gradB, gradQMrl, gradAA, gradC]' * (iniA^2 * exp(-t / iniB^2 - 0.5 * iniQMrl^2 * t.^2) + iniAA^2 * exp(-t / iniB^2) + iniC^2 - Et);
% Adagrad 法 x=x + yita * inv(G) * grad;
G = G + grad.^2;
newW = initW - alpha * diag(1 https://zhuanlan.zhihu.com/p/ sqrt(G + e)) * grad;
% newW=initW - alpha * grad;
if max(abs((newW - initW))) < err
convIter = convIter + 1;
if convIter == 10
J(i + 1 : end) = [];
disp(['梯度下降法迭代次数:', num2str(i)]);
result = newW';
result(1) = newW(1)^2;
result(2) = newW(2)^2;
result(3) = newW(3)^2;
result(4) = newW(4)^2;
result(5) = newW(5)^2;
disp(['梯度下降法迭代结果:', num2str(result)]);
break;
end
else
convIter = 0;
initW = newW;
iniA = newW(1);
iniB = newW(2);
iniQMrl = newW(3);
iniAA = newW(4);
iniC = newW(5);
J(i) = 1 / (2 * M) * norm(iniA^2 * exp(-t / iniB^2 - 0.5 * iniQMrl^2 * t.^2) + iniAA^2 * exp(-t / iniB^2) + iniC^2 - Et)^2;
end
end
fit = iniA^2 * exp(-t / iniB^2 - 0.5 * iniQMrl^2 * t.^2) + iniAA^2 * exp(-t / iniB^2) + iniC^2;
% 绘图
figure('Position', [50, 50, 900, 400]);
subplot(1, 2, 1)
loglog(J, 'LineWidth', 2)
legend(['alpha=', num2str(alpha)]);
subplot(1, 2, 2)
semilogx(t, Et, 'ok', 'MarkerFaceColor', 'r', 'MarkerSIze', 5, 'LineWidth', 1.5)
hold on
plot(t, fit, 'LineWidth', 1.5)
hold off
legend('Model', 'Fitting')
title('梯度下降法曲线拟合');
粒子状态更新
clear;clc;
% 群体例子个数
N = 100;
% 粒子维度
D = 20;
T = 200; % 最大迭代次数
c1 = 1.5; % 学习因子1
c2 = 1.5; % 学习因子2
wMax = 0.8; % 惯性权重最大值
wMin = 0.4; % 惯性权重最小值
xMax = 10; % 位置最大值
xMin = -10; % 位置最小值
vMax = 10; % 速度最大值
vMin = -10; % 速度最小值
%%%%%%%%%% 初始化种群个体(限定位置和速度) %%%%%%%%%
% 随机获得二进制编码的初始种群
x = randi([0, 1], N, D);
v = rand(N, D) * (vMax - vMin) + vMin;
%%%%%%%%%% 初始化个体最优位置和最优值 %%%%%%%%%
p = x;
pbest = ones(N, 1);
for i = 1:N
pbest(i) = ObjFunDiscrete(x(i, :), xMax, xMin);
end
%%%%%%%%%% 初始化全局最优位置和最优值 %%%%%%%%%
g = ones(1, D);
gbest = inf;
for i = 1:N
if pbest(i) < gbest
g = p(i, :);
gbest = pbest(i);
end
end
gb = ones(1, T);
vx = zeros(N, D);
%%%%%%%%%% 按照公式依次迭代直到满足精度或者迭代次数 %%%%%%%%%
for i = 1:T
for j = 1:N
%%%%%%%%%% 更新个体最优位置和最优值 %%%%%%%%%
if ObjFunDiscrete(x(j, :), xMax, xMin) < pbest(j)
p(j, :) = x(j, :);
pbest(j) = ObjFunDiscrete(x(j, :), xMax, xMin);
end
%%%%%%%%%% 更新全局最优位置和最优值 %%%%%%%%%
if pbest(j) < gbest
g = p(j, :);
gbest = pbest(j);
end
%%%%%%%%%% 动态计算惯性权重值 %%%%%%%%%
w = wMax - (wMax - wMin) * i / T;
%%%%%%%%%% 更新位置和速度值 %%%%%%%%%
v(j, :) = w * v(j, :) + c1 * rand * (p(j, :) - x(j, :)) + c2 * rand * (g - x(j, :));
x(j, :) = x(j, :) + v(j, :);
%%%%%%%%%% 边界条件处理 %%%%%%%%%
for ii = 1:D
if (v(j, ii) > vMax) || (v(j, ii) < vMin)
v(j, ii) = rand * (vMax - vMin) + vMin;
end
end
vx(j, :) = 1 https://zhuanlan.zhihu.com/p/ (1 + exp(-v(j, :)));
for jj = 1 : D
if vx(j, jj) > rand
x(j, jj) = 1;
else
x(j, jj) = 0;
end
end
end
%%%%%%%%%% 记录历代全局最优值 %%%%%%%%%
gb(i) = gbest;
end
m = 0;
for j = 1:D
m = m + g(j) * 2^(j - 1);
end
f1 = xMin + m * (xMax - xMin)/(2^D - 1);
disp(['最优个体:' num2str(f1)]);
disp(['最优值:' num2str(gb(end))]);
subplot(1, 2, 1);
plot(gb, 'LineWidth', 2);
xlabel('迭代次数');
ylabel('适应度值');
title('适应度进化曲线');
subplot(1, 2, 2)
xRange = linspace(xMin, xMax, 100);
yRange = xRange + 6 * sin(4 * xRange) + 9 * cos(5 * xRange);
plot(xRange, yRange, 'LineWidth', 2);
hold on;
plot(f1, gb(end),'o', 'MarkerFaceColor', 'r', 'MarkerSize', 10);
hold off;
set(gcf, 'Position', [400, 200, 900, 400]);
function result = ObjFunDiscrete(x, xMax, xMin)
%%%%%%%%%% 适应度函数 %%%%%%%%%%%
m = 0;
D = length(x);
for j = 1:D
m = m + x(j) * 2^(j - 1);
end
% 译码成十进制数
f = xMin + m * (xMax - xMin) / (2^D - 1);
fit = f + 6 * sin(4 * f) + 9 * cos(5 * f);
result = fit;
end
粒子群算法实现指数拟合(PSO-ExpFit)
clear;clc;
t = 0.2*(1:3000)';
data = 400*exp(-t/5) + 10*randn(size(t));
tic;
p = PSO_ExpFit(t, data)
toc;
fit = p(1)*exp(-t/p(2)) + p(3);
plot(t, data, t, fit, 'LineWidth', 2)
legend('model', 'PSO');
title('压缩因子粒子群算法实现指数拟合');
function [x_opt, y_opt]= PSO_ExpFit(t, Et)
%{
函数功能:压缩因子粒子群算法实现指数拟合:y=a*exp(-x/b) + c
输入:
t:自变量;
Et:因变量;
输出:
x_opt:最优解;
y_opt:适应度(目标函数值);
%}
% 初始值
D = 3; % 粒子维度,未知数个数;
IterMax = 200; % 最大迭代次数;
Vmin = -1; % 速度最小值;
Vmax = 1; % 速度最大值;
c1 = 2.05; % 学习因子1,认知部分;
c2 = 2.05; % 学习因子2,社会部分;
N = 200; % 种群个数;
max_sig = max(Et);
Xmin = [0.5*max_sig, -2, 0]; % 下限
Xmax = [1.5*max_sig, 4, 0.5*max_sig]; % 上限;
% 计算
phy = c1 + c2;
lamda = 2 / abs(2 - phy - sqrt(phy^2 - 4*phy)); % 压缩因子
% 初始化种群个体(限定位置和速度)
x = rand(N, D).*repmat(Xmax - Xmin, N, 1) + repmat(Xmin, N, 1);
v = rand(N, D)*(Vmax - Vmin) + Vmin;
% 初始化个体最优位置和最优值
p = x;
pbest = ones(1, N);
for i = 1 : N
pbest(i) = Obj_Fit(x(i, :), t ,Et);
end
% 初始化全局最优位置和最优值
g = ones(1, D);
gbest = inf;
for i = 1 : N
if pbest(i) < gbest
g = p(i, :);
gbest = pbest(i);
end
end
% 迭代直到满足精度或者迭代次数
fx = zeros(1, N);
y_opt = zeros(1, IterMax);
for i = 1 : IterMax
for j = 1 : N
% 更新个体最优位置和最优值
fx(j) = Obj_Fit(x(j, :), t ,Et);
if fx(j) < pbest(j)
p(j, :) = x(j, :);
pbest(j) = fx(j);
end
% 更新全局最优位置和最优值
if pbest(j) < gbest
g = p(j, :);
gbest = pbest(j);
end
% 更新位置和速度值
v(j, :) = lamda*v(j, :) + c1*rand*(p(j, :) - x(j, :)) + c2*rand*(g - x(j, :));
x(j, :) = x(j, :) + v(j, :);
% 边界条件处理
for ii = 1 : D
if (v(j, ii) > Vmax) || (v(j, ii) < Vmin)
v(j, ii) = rand*(Vmax - Vmin) + Vmin;
end
if (x(j, ii) > Xmax(ii)) || (x(j, ii) < Xmin(ii))
x(j, ii) = rand*(Xmax(ii) - Xmin(ii)) + Xmin(ii);
end
end
end
end
x_opt = g;
x_opt(2) = 10^(x_opt(2));
end
% 目标函数
function fitError = Obj_Fit(p0, t ,Et)
A = p0(1);
B = p0(2);
C = p0(3);
f = A*exp(-t/10^B) + C;
fitError = norm(Et - f);
end
clear;clc;
lb = [-10, -10];
ub = [10, 10];
x0 = rand(1, 2) .* (ub - lb) + lb;
rng('default');
[mx, minf] = OptRandWalk(x0, 10, 100, 10);
disp(['最优值:' num2str(mx)]);
disp(['最优值对应目标函数:' num2str(minf)]);
len = 100;
xRange = linspace(-10, 10, len);
yRange = linspace(-10, 10, len);
[xMap, yMap] = meshgrid(xRange, yRange);
zMap = zeros(len);
for i = 1 : len
for j = 1 : len
zMap(i, j) = func([xMap(i, j), yMap(i, j)]);
end
end
surf(xRange, yRange, zMap);
view(-45, -45);
shading interp
hold on
plot3(mx(1), mx(2), minf, 'o', 'MarkerFaceColor', 'r', 'MarkerSize', 10);
hold off
title('随机行走法求函数的极小值');
function[mx, minf]= OptRandWalk(x, lamda, N, n)
%{
函数功能:随机行走法求函数的极小值
x:初始值;
lamda:步长;
N:为了产生较好点的迭代次数;
n:单步循环行走次数,目的是尽可能走到全局最优点附近
mx:最优解;
minf:最优值。
%}
F = zeros(n, 1);
D = length(x);
X = zeros(n, D);
epsilon = 1e-6;
f1 = func(x);
while lamda >= epsilon
k = 1;
while(k <= N)
u = 2*rand(n, D) - 1;
for ii =1 : n
X(ii, :) = x + lamda*(u(ii, :) / norm(u(ii, :)));
F(ii) = func(X(ii, :));
end
[f11, kk] = min(F);
if f11 < f1
f1 = f11;
x = X(kk, :);
k = 1;
else
k = k + 1;
end
end
lamda = lamda / 2;
end
mx = X(kk, :);
minf = f1;
end
function f = func(x)
% f=-sin(sqrt((x(1) - 100).^2 + (x(2) - 100).^2 ) + exp(1)) https://zhuanlan.zhihu.com/p/ (sqrt((x(1) - 100).^2 + (x(2) - 100).^2 ) + exp(1)) - 1;
% f=sum(x.^2);
% f=3*cos(x(1)*x(2)) + x(1) + x(2)^2;
% f=4*x(1)^2-2.1*x(1)^4+x(1)^6/3+x(1)*x(2)-4*x(2)^2+4*x(2)^4;
% f=0.5 + (sin(sqrt(x(1)^2 + x(2)^2))^2 - 0.5) / (1 + 0.001*(x(1)^2 + x(2)^2))^2;
x1 = x(1);
x2 = x(2);
t = x1^2 + x2^2;
f = 0.5 + (sin(sqrt(t))^2 - 0.5) / (1 + 0.001 * t)^2;
end
clear;clc;
t = 0.2*(1:3000)';
data = 400*exp(-t/5) + 10*randn(size(t));
tic;
p = SAA_ExpFit(t, data)
toc;
fit = p(1)*exp(-t/p(2)) + p(3);
plot(t, data, t, fit, 'LineWidth', 2)
legend('model', 'SA');
title('模拟退火算法实现指数拟合');
function p0 = SAA_ExpFit(t, Et)
%{
函数功能:模拟退火算法实现指数拟合:y=a*exp(-x/b) + c
输入:
t:自变量;
Et:因变量;
输出:
x_opt:最优解;
y_opt:适应度(目标函数值);
%}
% initial condition
T0= 1000; % initial temperature
alpha= 0.99; % annealing factor
N_max= 3000; % max iterations
% range of varibles
max_sig = max(Et);
Xmin = [0.5*max_sig, -2, 0]; % 下限
Xmax = [1.5*max_sig, 4, 0.5*max_sig]; % 上限;
% initial value
var = rand(1,3).*(Xmax - Xmin) + Xmin;
fun = Obj_Fit(var, t, Et);
% iterative
for i= 1:N_max
Tk= T0*alpha^(i-1);
dvar = Tk*sign(rand(1,3)-0.5).*((1+1/Tk).^abs(2*rand(1,3)-1)-1);
var1 = var + dvar.*(Xmax - Xmin);
% accept it or not
if all(var1 >= Xmin) && all(var1 <= Xmax)
new_fun = Obj_Fit(var1, t, Et);
if (new_fun <= fun)
var = var1;
fun = new_fun;
else
if exp(-(new_fun-fun)/Tk)>rand
var = var1;
fun = new_fun;
end
end
end
end
% results
p0 = [var(1), 10^var(2), var(3)]';
end
function fitError = Obj_Fit(p0, t ,Et)
A1 = p0(1);
B1 = p0(2);
A0 = p0(3);
f = A1*exp(-t/10^(B1 + eps)) + A0;
fitError = norm(Et - f);
end
clear;clc;
ratio = 0.99;
t_start = 100;
t_end = 1;
Markov_length = 5000;
coordinates = load('CityCoordinates.txt');
amount = size(coordinates, 1); % 城市的数目
% 通过向量化的方法计算距离矩阵
coor_x_tmp1 = coordinates(:,1) * ones(1, amount);
coor_x_tmp2 = coor_x_tmp1';
coor_y_tmp1 = coordinates(:, 2) * ones(1, amount);
coor_y_tmp2 = coor_y_tmp1';
dist_matrix = sqrt((coor_x_tmp1 - coor_x_tmp2).^2 + (coor_y_tmp1 - coor_y_tmp2).^2);
sol_new = 1 : amount; % 产生初始解
% sol_new 是每次产生的新解;
% sol_current 是当前解;
% sol_best 是冷却中的最好解;
% E_current 是当前解对应的回路距离;
% E_new 是新解的回路距离;
% E_best 是最优解的
E_current = inf;
E_best = inf;
sol_current = sol_new;
sol_best = sol_new;
t = t_start;
while t >= t_end
% 当前温度下循环
for r = 1 : Markov_length
% 产生随机扰动
if (rand < 0.5) % 随机决定是进行两交换还是三交换
% 两交换
ind = randperm(amount, 2);
ind1 = ind(1);
ind2 = ind(2);
tmp1 = sol_new(ind1);
sol_new(ind1) = sol_new(ind2);
sol_new(ind2) = tmp1;
else
% 三交换
ind = randperm(amount, 3);
% 确保ind1 < ind2 < ind3
ind = sort(ind);
ind1 = ind(1);
ind2 = ind(2);
ind3 = ind(3);
tmplist1 = sol_new(ind1 + 1 : ind2 - 1);
sol_new(ind1 + 1 : ind1 + ind3 - ind2 + 1) = sol_new(ind2 : ind3);
sol_new(ind1 + ind3 - ind2 + 2 : ind3) = tmplist1;
end
% 计算目标函数值
E_new = 0;
for i = 1 : amount - 1
E_new = E_new + dist_matrix(sol_new(i), sol_new(i + 1));
end
% 再算上从最后一个城市到第一个城市的距离
E_new = E_new + dist_matrix(sol_new(amount), sol_new(1));
if E_new < E_current
E_current = E_new;
sol_current = sol_new;
if E_new < E_best
% 把冷却过程中最好的解保存下来
E_best = E_new;
sol_best = sol_new;
end
else
% 若新解的目标函数值大于当前解的,则仅以一定概率接受新解
if rand < exp(-(E_new - E_current)/t)
E_current = E_new;
sol_current = sol_new;
else
sol_new = sol_current;
end
end
end
t = t * ratio; % 降温
end
disp(['最优解为:', num2str(sol_best)]);
disp(['最短距离:', num2str(E_best)]);
x_path = [coordinates(sol_best, 1); coordinates(sol_best(1), 1)];
y_path = [coordinates(sol_best, 2); coordinates(sol_best(1), 2)];
plot(x_path, y_path, 'LineWidth', 2)
算法下载:
SSA: Salp Swarm Algorithmclear;clc;
% 搜索数量
SearchAgents_no = 30;
% 最大迭代次数
Max_iter = 200;
% 变量维度
dim = 2;
% 上限
ub = 10 * ones(1, dim);
% 下限
lb = -10 * ones(1, dim);
% 初始化位置
Positions = rand(SearchAgents_no, dim) .* repmat(ub - lb, SearchAgents_no, 1) + repmat(lb, SearchAgents_no, 1);
Convergence_curve = zeros(1, Max_iter);
SalpFitness = nan(1, Max_iter);
Sorted_salps = nan(SearchAgents_no, dim);
% 计算适应度值
for i = 1 : SearchAgents_no
SalpFitness(1, i) = ObjFun(Positions(i, :));
end
[sorted_salps_fitness, sorted_indexes] = sort(SalpFitness);
for i = 1:SearchAgents_no
Sorted_salps(i, :) = Positions(sorted_indexes(i), :);
end
FoodPosition = Sorted_salps(1, :);
FoodFitness = sorted_salps_fitness(1);
Convergence_curve(1) = FoodFitness;
% 主循环
for l = 2:Max_iter
c1 = 2*exp(-(4*l/Max_iter)^2);
for i = 1:SearchAgents_no
SalpPositions = Positions';
if i <= SearchAgents_no/2
for j = 1:dim
c2 = rand();
c3 = rand();
if c3 < 0.5
SalpPositions(j, i) = FoodPosition(j) + c1*((ub(j) - lb(j))*c2 + lb(j));
else
SalpPositions(j, i) = FoodPosition(j) - c1*((ub(j) - lb(j))*c2 + lb(j));
end
end
elseif i > SearchAgents_no/2
point1 = SalpPositions(: , i-1);
point2 = SalpPositions(: , i);
SalpPositions(:, i) = (point2 + point1) / 2;
end
Positions= SalpPositions';
end
% 更新
for i = 1:SearchAgents_no
% 边界处理:越界赋值为边界值
Flag4ub = Positions(i, :) > ub;
Flag4lb = Positions(i, :) < lb;
Positions(i, :) = (Positions(i, :).*(~(Flag4ub + Flag4lb))) + ub.*Flag4ub + lb.*Flag4lb;
SalpFitness(i) = ObjFun(Positions(i, :));
if SalpFitness(i) < FoodFitness
FoodPosition = Positions(i, :);
FoodFitness = SalpFitness(i);
end
end
Convergence_curve(l) = FoodFitness;
end
disp(['最优值:' num2str(FoodPosition)]);
disp(['最优值对应目标函数:' num2str(FoodFitness)]);
subplot(1, 2, 1)
plot(Convergence_curve, 'LineWidth', 2);
xlabel('迭代次数');
ylabel('适应度值');
title('适应度曲线');
len = 50;
xRange = linspace(lb(1), ub(1), len);
yRange = linspace(lb(2), ub(2), len);
[xMap, yMap] = meshgrid(xRange, yRange);
zMap = zeros(len);
for i = 1 : len
for j = 1 : len
zMap(i, j) = ObjFun([xMap(i, j), yMap(i, j)]);
end
end
subplot(1, 2, 2)
surf(xRange, yRange, zMap);
view(-45, -45);
shading interp
hold on
plot3(FoodPosition(1), FoodPosition(2), FoodFitness, 'o', 'MarkerFaceColor', 'r', 'MarkerSize', 10);
hold off
set(gcf, 'Position', [400, 200, 900, 400]);
function result = ObjFun(x)
x1 = x(1);
x2 = x(2);
t = x1^2 + x2^2;
result = 0.5 + (sin(sqrt(t))^2 - 0.5) / (1 + 0.001 * t)^2;
end
算法下载:
https://au.mathworks.com/matlabcentral/fileexchange/55667-the-whale-optimization-algorithmThe Whale Optimization Algorithm
The Whale Optimization Algorithmclear;clc;
% 鲸鱼数量
SearchAgents_no = 30;
% 最大迭代次数
Max_iter = 200;
% 变量维度
dim = 2;
% 上限
ub = 10 * ones(1, dim);
% 下限
lb = -10 * ones(1, dim);
%%%%% 初始化
Leader_pos = zeros(1, dim);
Leader_score = inf;
% 初始化位置
Positions = rand(SearchAgents_no, dim) .* repmat(ub - lb, SearchAgents_no, 1) + repmat(lb, SearchAgents_no, 1);
Convergence_curve = zeros(1, Max_iter);
% 主循环
for t = 1:Max_iter
for i = 1:SearchAgents_no
% 边界处理:越界赋值为边界值
Flag4ub = Positions(i, :) > ub;
Flag4lb = Positions(i, :) < lb;
Positions(i, :) = (Positions(i, :).*(~(Flag4ub + Flag4lb))) + ub.*Flag4ub + lb.*Flag4lb;
% 计算目标函数
fitness = ObjFun(Positions(i, :));
% 更新领头位置
if fitness < Leader_score
Leader_score = fitness;
Leader_pos = Positions(i, :);
end
end
% 线性递减
a = 2 - t * 2 / Max_iter;
a2 = -1 - t / Max_iter;
% 更新位置
for i = 1 : SearchAgents_no
r1 = rand();
r2 = rand();
A = 2*a*r1 - a;
C = 2*r2;
b = 1;
l = (a2-1)*rand+1;
p = rand();
for j = 1 : dim
if p < 0.5
if abs(A) >= 1
rand_leader_index = floor(SearchAgents_no*rand() + 1);
X_rand = Positions(rand_leader_index, :);
D_X_rand = abs(C*X_rand(j) - Positions(i, j));
Positions(i, j) = X_rand(j) - A*D_X_rand;
elseif abs(A) < 1
D_Leader = abs(C*Leader_pos(j) - Positions(i, j));
Positions(i, j) = Leader_pos(j) - A*D_Leader;
end
else
distance2Leader = abs(Leader_pos(j) - Positions(i, j));
Positions(i, j) = distance2Leader*exp(b.*l).*cos(l.*2*pi) + Leader_pos(j);
end
end
end
Convergence_curve(t) = Leader_score;
end
disp(['最优值:' num2str(Leader_pos)]);
disp(['最优值对应目标函数:' num2str(Leader_score)]);
subplot(1, 2, 1)
plot(Convergence_curve, 'LineWidth', 2);
xlabel('迭代次数');
ylabel('适应度值');
title('适应度曲线');
len = 50;
xRange = linspace(lb(1), ub(1), len);
yRange = linspace(lb(2), ub(2), len);
[xMap, yMap] = meshgrid(xRange, yRange);
zMap = zeros(len);
for i = 1 : len
for j = 1 : len
zMap(i, j) = ObjFun([xMap(i, j), yMap(i, j)]);
end
end
subplot(1, 2, 2)
surf(xRange, yRange, zMap);
view(-45, -45);
shading interp
hold on
plot3(Leader_pos(1), Leader_pos(2), Leader_score, 'o', 'MarkerFaceColor', 'r', 'MarkerSize', 10);
hold off
set(gcf, 'Position', [400, 200, 900, 400]);
function result = ObjFun(x)
x1 = x(1);
x2 = x(2);
t = x1^2 + x2^2;
result = 0.5 + (sin(sqrt(t))^2 - 0.5) / (1 + 0.001 * t)^2;
end