【优化求解】基于自适应风驱动算法求解多目标优化问题含Matlab源码

【优化求解】基于自适应风驱动算法求解多目标优化问题含Matlab源码

1 简介

2 部分代码

function [rec] = purecmaes_wdo(counteval,rec,npop,pressure)

%   coeff -- CMAES optimized coefficients returned to WDO.

%   counteval -- Iteration counter.

%   rec -- Record of prior values used in CMAES.

%   npop -- number of population members from WDO, each member gets their

%          own set of coefficients determined by the CMAES.

%   pressure -- pressure(cost function) computed by WDO for the set of

%               coefficients that CMEAS picked last iteration

%---------------------------------------------------------

%   This script is modified version of the publicly available purecmaes.m

%   Original "purecmaes.m" can be accessed from:

%   URL: .m

%

%   This function takes in pressure values from the WDO and computes new

%   set of coefficients for the WDO and returns them. In simpler terms, it

%   optimizes the WDO coefficients of "alpha, g, RT, c"

%

%   This function is not optimized, only proof-of-concept code to

%   illustrate that the WDO coefficients can be optimized by CMAES, 

%   which in turn creates and Adaptive WDO algorithm.

%---------------------------------------------------------

if counteval==2    %initialization only happens when the CMAES called for the first time.

    rec.N = 5;  % problem dimension is fixed to 4: alp, g, RT, c 

    an = rand(rec.N,1);  % objective variables initial point

    rec.sigma = 0.5;        % coordinate wise standard deviation (step size)

  

    % Strategy parameter setting: Selection  

%     lambda = 4+floor(3*log(N));  % population size, offspring number

    rec.lambda = npop; %this is defined by the wdo population size.

    rec.mu = rec.lambda/2;               % number of parents/points for recombination

    rec.weights = log(rec.mu+1/2)-log(1:rec.mu)'; % muXone array for weighted recombination

    rec.mu = floor(rec.mu);        

    rec.weights = rec.weights/sum(rec.weights);     % normalize recombination weights array

    rec.mueff=sum(rec.weights)^2/sum(rec.weights.^2); % variance-effectiveness of sum w_i x_i

    % Strategy parameter setting: Adaptation

    rec = (4 + rec.mueff/rec.N) / (rec.N+4 + 2*rec.mueff/rec.N); % time constant for cumulation for C

    rec.cs = (rec.mueff+2) / (rec.N+rec.mueff+5);  % t-const for cumulation for sigma control

    rec.c1 = 2 / ((rec.N+1.3)^2+rec.mueff);    % learning rate for rank-one update of C

    u = min(1-rec.c1, 2 * (rec.mueff-2+1/rec.mueff) / ((rec.N+2)^2+rec.mueff));  % and for rank-mu update

    rec.damps = 1 + 2*max(0, sqrt((rec.mueff-1)/(rec.N+1))-1) + rec.cs; % damping for sigma 

                                                      % usually close to 1

    % Initialize dynamic (internal) strategy parameters and constants

    rec.pc = zeros(rec.N,1); 

    rec.ps = zeros(rec.N,1);   % pc, ps: evolution paths for C and sigma

    rec.B = eye(rec.N,rec.N);                       % B defines the coordinate system

    rec.D = ones(rec.N,1);                      % diagonal D defines the scaling

    rec.C = rec.B * diag(rec.D.^2) * rec.B';            % covariance matrix C

    rec.invsqrtC = rec.B * diag(rec.D.^-1) * rec.B';    % C^-1/2 

    rec.eigeneval = 0;                      % track update of B and D

    rec.chiN=rec.N^0.5*(1-1/(4*rec.N)+1/(21*rec.N^2));  % expectation of  ||N(0,I)|| == norm(randn(N,1))

end

    % get the fitness from WDO pressure:

    rec.arfitness = pressure';

    

    % Sort by fitness and compute weighted mean into xmean

    [rec.arfitness, rec.arindex] = sort(rec.arfitness);  % minimization

    ld = an;

    an = rec.arx(:,rec.arindex(1:rec.mu)) * rec.weights;  % recombination, new mean value

    

    % Cumulation: Update evolution paths

    rec.ps = (1-rec.cs) * rec.ps ... 

          + sqrt(rec.cs*(2-rec.cs)*rec.mueff) * rec.invsqrtC * (ld) / rec.sigma; 

    rec.hsig = sum(rec.ps.^2)/(1-(1-rec.cs)^(2*counteval/rec.lambda))/rec.N < 2 + 4/(rec.N+1);

    rec.pc = (1-rec) * rec.pc ...

          + rec.hsig * sqrt(rec*(2-rec)*rec.mueff) * (ld) / rec.sigma; 

    % Adapt covariance matrix C

    rec.artmp = (1/rec.sigma) * (rec.arx(:,rec.arindex(1:rec.mu)) - ld,1,rec.mu));  % mu difference vectors

    rec.C = (u) * rec.C ...                   % regard old matrix  

         + rec.c1 * (rec.pc * rec.pc' ...                % plus rank one update

                 + (1-rec.hsig) * rec*(2-rec) * rec.C) ... % minor correction if hsig==0

         + u * rec.artmp * diag(rec.weights) * rec.artmp'; % plus rank mu update 

    % Adapt step size sigma

    rec.sigma = rec.sigma * exp((rec.cs/rec.damps)*(norm(rec.ps)/rec.chiN - 1)); 

    

    % Update B and D from C

    if counteval - rec.eigeneval > rec.lambda/(rec.c1&#u)/rec.N/10  % to achieve O(N^2)

      rec.eigeneval = counteval;

      rec.C = triu(rec.C) + triu(rec.C,1)'; % enforce symmetry

      [rec.B,rec.D] = eig(rec.C);           % eigen decomposition, B==normalized eigenvectors

      rec.D = sqrt(diag(rec.D));        % D contains standard deviations now

      rec.invsqrtC = rec.B * diag(rec.D.^-1) * rec.B';

    end

    

    % Generate and evaluate lambda offspring

    for k=1:rec.lambda,

        rec.arx(:,k) = an + rec.sigma * rec.B * (rec.D .* randn(rec.N,1)); % m + sig * Normal(0,C) 

    end

end

% end-of-CMAES.

3 仿真结果

4 参考文献

[1]陈彬彬, 曹中清, 余胜威. 基于风驱动优化算法WDO的PID参数优化[J]. 计算机工程与应用, 2016, 52(014):250-253.

[2]寻琦, 赵惠玲, and 王肇平. "一种基于自适应风驱动优化算法的天线阵方向图综合方法.". 

​博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。

部分理论引用网络文献，若有侵权联系博主删除。