《多尺度量子谐振子优化算法》中的代码
解释很有创意的一种优化算法。
01 MQHOA算法(i-MQHOA算法)基本结构:
Initialize: k,
σ
m
i
n
\sigma_{min}
σmin,
σ
s
\sigma_{s}
σs = max - min
Randomly generate
x
i
(
i
=
1
,
⋯
 
,
k
)
x_i(i= 1, \cdots, k)
xi(i=1,⋯,k) in domain
Caculate the standard deviation
σ
k
′
\sigma_k^{\prime}
σk′ for all
x
i
x_i
xi
D
o
Do
Do
D
o
\quad Do
Do
D
o
\quad \quad Do
Do
∀
x
i
,
\qquad \quad \forall x_i,
∀xi, generate
x
i
′
∼
N
(
x
i
,
σ
s
2
)
x_i^{\prime} \sim N(x_i, \sigma _{s}^{2})
xi′∼N(xi,σs2)
∀
x
i
\qquad \quad \forall x_i
∀xi and
x
i
′
,
x_i^{\prime},
xi′, if
f
(
x
i
′
)
<
f
(
x
i
)
f(x_i^{\prime}) < f(x_i)
f(xi′)<f(xi) then
x
i
=
x
i
′
x_i = x_i^{\prime}
xi=xi′
\qquad \quad
Calculate the standard deviation
σ
k
\sigma_k
σk
Δ
σ
k
=
∣
σ
k
−
σ
k
′
∣
\qquad \quad \Delta \sigma_k = | \sigma_k - \sigma_k^{\prime} |
Δσk=∣σk−σk′∣
\qquad \quad
Update
σ
k
′
:
σ
k
′
←
σ
k
\sigma_k^{\prime}:\sigma_k^{\prime} \leftarrow \sigma_k
σk′:σk′←σk
\qquad
While(
Δ
σ
k
>
σ
s
\Delta \sigma_k > \sigma_s
Δσk>σs)
\qquad
Update the worst solution
x
w
o
r
s
t
=
x
m
e
a
n
x^{worst} = x^{mean}
xworst=xmean
\quad
While(
σ
k
>
σ
s
\sigma_k > \sigma_s
σk>σs)
\quad
Update
σ
s
:
σ
s
←
σ
s
2
\sigma_s : \sigma_s \leftarrow \cfrac{\sigma_s}{2}
σs:σs←2σs
While(
σ
s
>
σ
m
i
n
\sigma_s > \sigma_{min}
σs>σmin)
Output
x
b
e
s
t
x^{best}
xbest and
f
(
x
b
e
s
t
)
f(x^{best})
f(xbest)
02 matlab demo
func.m
function y = func(x, DIM)
% ------------------------------------------
% Griewank function% The global minima: f(x) = 0
% Position(0)
% ------------------------------------------
n = DIM;
fr = 4000;
s = 0;
p = 1;
for j = 1:n; s = s + x(j)^2; end
for j = 1:n; p = p*cos(x(j)/sqrt(j));end
y = s/fr-p+1;
end
MQHOA.m
% ---------------------------------------------------
% MQHOA algorithm for sigle model objective function
% Matlab Code by Wang Peng (Dec.21, 2015)
% This program must include func.m
% This is global minimum program.
% All rights reserved by Parallel Computing Lab.
% ---------------------------------------------------
clear all;
feature jit off; % shut down the accelerator
format long;
global DIM;
% -------参数定义区域开始
repeat = 10; % 定义重复计算次数
DIM = 300; % objective function's dimentions
ACCURACY = 0.000001; % 定义计算精度,sigma_min
optimalNum = 15; % k值的大小
minDomain = -10; maxDomain = 10; % objective function's domain
% -------参数定义区域结束
% -------多次重复计算所需参数
gmeanV = zeros(1, repeat);
gbestV = zeros(1, repeat);
gfangcha = zeros(1, repeat);
% -------多次重复计算
for kk = 1: repeat
% --参数初始化区域开始
funcV = zeros(1, optimalNum); % 定义保存k个采样点的函数值,初始化为0
samplePos = zeros(DIM, optimalNum); % 采样点矩阵
sigma = maxDomain - minDomain; % 初始尺度为目标函数定义域的大小
stdPre = zeros(1, DIM); % 存储每一个维度上的标准差
stdNow = zeros(1, DIM);
% -- 计算k个采样位置的目标函数值
% 定义k个DIM维采样点的坐标,并初始化
optimalSolution = unifrnd(minDomain, maxDomain, DIM, optimalNum);
stdPre = std(optimalSolution, 1, 2); % 计算初始化k个采样点的标准差,按行求
w = 0; % function evolution times
for k = 1: optimalNum % 求最优解函数值
funcV(k) = func(optimalSolution(:,k), DIM);
w = w + 1;
end
while 1 % M iteration begin
while 1 % QHO iteration begin
while 1 %稳定性收敛迭代开始
for i = 1: optimalNum
% 采用Box-Muller方法生成新的正态分布随机采样点
theat = 2*pi*rand(DIM, 1);
R = sqrt(-2.0*log(rand(DIM, 1)));
gaosiRand = R.*cos(theat);
samplePos(:, i) = optimalSolution(:, i) + sigma*gaosiRand;
% 采用Box-Muller方法生成新的正态分布采样点
% 处理越界采样点
for i1 = 1: DIM
if samplePos(i1, i) > maxDomain
samplePos(i1, i) = maxDomain;
end
if samplePos(i1, i) < minDomain
samplePos(i1, i) = minDomain;
end
end
sampleValue = func(samplePos(:, i), DIM); % 求第i个采样点的函数值
w = w + 1;
if sampleValue < funcV(i) % 如果采样点值小于当前点函数值,则替换
funcV(i) = sampleValue;
optimalSolution(:, i) = samplePos(:, i);
end
end
stdNow = std(optimalSolution, 1, 2); % 新解标准差
c = stdPre - stdNow; % 两次最优解标准差的差值
stdPre = stdNow;
if max(abs(c)) <= (sigma) % 如果所有维中的标准差最大的值都小于Sigma,则表明在该能级采样已稳定
break;
end
end % 稳定性收敛迭代结束
% -- 以下代码将k个采样点中均值的位置合并到最大值的位置中,实现能级的下降
meanV = mean(funcV);
meanPos = mean(optimalSolution, 2); % 取得平均坐标
[v_max, index_max] = max(funcV); % 取得最大值的序号index_max
optimalSolution(:, index_max) = meanPos; % 用平均坐标替换最大值对应坐标
funcV(index_max) = func(meanPos, DIM);
w = w + 1;
stdPre = std(optimalSolution, 1, 2); % 新解标准差
% -- 能级下降结束
if max(stdPre) < sigma
break;
end
end % QHO iteration end
sigma = sigma/2.0;
if sigma <= ACCURACY
break;
end
end % M iteration end
[global_min, index] = min(funcV);
gmeanV(kk) = mean(funcV);
gbestV(kk) = min(funcV);
fprintf('No. %d run, DIM=%d, Global minimum function value is %e. function evolution = %d\n', kk, DIM, global_min, w);
end % MQHOA end
meanV = mean(gmeanV);
bestV = min(gbestV);
fangcha = std(gbestV);
fprintf('Repeat=%d, MeanValue=%d, BestValue=%d, Std=%d\n', repeat, meanV, bestV, fangcha);
运行结果:
No. 1 run, DIM=300, Global minimum function value is 7.284663e-10. function evolution = 381435
No. 2 run, DIM=300, Global minimum function value is 5.056802e-10. function evolution = 389947
No. 3 run, DIM=300, Global minimum function value is 7.396041e-03. function evolution = 381926
No. 4 run, DIM=300, Global minimum function value is 1.232099e-02. function evolution = 358647
No. 5 run, DIM=300, Global minimum function value is 4.499122e-10. function evolution = 392409
No. 6 run, DIM=300, Global minimum function value is 7.396041e-03. function evolution = 389146
No. 7 run, DIM=300, Global minimum function value is 6.205952e-10. function evolution = 388832
No. 8 run, DIM=300, Global minimum function value is 5.278187e-10. function evolution = 381695
No. 9 run, DIM=300, Global minimum function value is 4.817348e-10. function evolution = 399043
No. 10 run, DIM=300, Global minimum function value is 4.652333e-10. function evolution = 384442
Repeat=10, MeanValue=2.711308e-03, BestValue=4.499122e-10, Std=4.566761e-03