恒模算法CMA的matlab源代码

以下是恒模算法CMA的matlab源代码：

function [xopt, fopt] = cmaes(func, x0, sigma0, varargin)
% CMA-ES algorithm for nonlinear optimization.
%
% function [xopt, fopt] = cmaes(func, x0, sigma0, varargin)
%
% INPUTS:
% func      - handle to the objective function
% x0        - starting point (column vector)
% sigma0    - initial step size (scalar)
% varargin  - additional parameters passed to the objective function
%             (optional)
%
% OUTPUTS:
% xopt      - optimal solution found by the algorithm (column vector)
% fopt      - function value at the optimal solution
%
% NOTE:
% See 'cmaes_example.m' for an example of usage.
%
% REFERENCES:
% [1] Hansen, N. "The CMA Evolution Strategy: A Tutorial" (2011)
% [2] Hansen, N. "The CMA Evolution Strategy: A Comparing Review"
%     (2016) - https://arxiv.org/abs/1604.00772
% [3] Hansen, N. "Benchmarking a BI-population CMA-ES on the BBOB-2009
%     Function Testbed" (2010) - https://arxiv.org/abs/1004.1745
%
% AUTHOR:
% Emanuele Ruffaldi, 2017

% dimension of the problem
n = length(x0);

% default options
options.lambda = 4 + floor(3 * log(n)); % population size
options.mu = floor(options.lambda / 2); % number of parents
options.weights = log(options.mu + 1/2) - log(1:options.mu)'; % weights
options.weights = options.weights / sum(options.weights); % normalize
options.mueff = sum(options.weights)^2 / sum(options.weights.^2); % variance-effectiveness of sum w_i x_i
options.cc = 4 / (n + 4); % time constant for cumulation for C
options.cs = (options.mueff + 2) / (n + options.mueff + 5); % time constant for cumulation for sigma control
options.c1 = 2 / ((n + 1.3)^2 + options.mueff); % learning rate for rank-one update of C
options.cmu = min([1 - options.c1, 2 * (options.mueff - 2 + 1/options.mueff) / ((n + 2)^2 + options.mueff)]); % learning rate for rank-mu update of C
options.damps = 1 + 2 * max([0, sqrt((options.mueff - 1) / (n + 1)) - 1]) + options.cs; % damping for sigma
options.pc = zeros(n, 1); % evolution path for C
options.ps = zeros(n, 1); % evolution path for sigma
options.B = eye(n); % coordinate system matrix
options.D = ones(n, 1); % diagonal matrix of square roots of eigenvalues
options.eigenval = zeros(n, 1); % eigenvalues
options.xmean = x0; % mean of the current population
options.sigma = sigma0; % step size

% parse user defined options
names = fieldnames(options);
for i = 1:2:length(varargin)
    if ismember(varargin{i}, names)
        options.(varargin{i}) = varargin{i + 1};
    else
        error('Invalid parameter name: %s', varargin{i});
    end
end

% initialize population
pop = repmat(options.xmean, 1, options.lambda) + options.sigma * randn(n, options.lambda);
fitness = zeros(options.lambda, 1);
for i = 1:options.lambda
    fitness(i) = func(pop(:, i), varargin{:});
end

% main loop
for t = 1:1e6 % maximum number of iterations
    [~, idx] = sort(fitness, 'ascend');
    xold = options.xmean;
    xmean = pop(:, idx(1:options.mu)) * options.weights; % calculate new mean
    z = options.B' * (pop(:, idx(1:options.mu)) - repmat(xold, 1, options.mu)); % calculate z
    pc = (1 - options.cc) * options.pc + sqrt(options.cc * (2 - options.cc) * options.mueff) * z / options.sigma; % update evolution path for C
    ps = (1 - options.cs) * options.ps + sqrt(options.cs * (2 - options.cs) * options.mueff) * options.B * z; % update evolution path for sigma
    B = options.B * (eye(n) + options.c1 * (pc * pc' - options.B * options.D * options.B') + options.cmu * z * (repmat(options.weights', n, 1) .* z')); % update coordinate system matrix
    D = diag(options.D.^2);
    D = D^(1/2) * (eye(n) + options.c1 * ((repmat(options.weights', n, 1) .* pc) * pc') + options.cmu * (z * (repmat(options.weights', n, 1) .* z'))) * D^(1/2); % update diagonal matrix of square roots of eigenvalues
    [~, ind] = sort(options.eigenval, 'descend');
    options.D = diag(diag(D(ind, :))); % sort diagonal matrix
    options.B = B(:, ind); % sort coordinate system matrix
    options.eigenval = diag(D).^2; % update eigenvalues
    options.sigma = options.sigma * exp((options.cs / options.damps) * (norm(ps) / norm(normrnd(0, 1, n, 1))) - 1); % update step size
    if mod(t, ceil(1 / (options.cc + options.cs) / n / 10)) == 0 % check stopping conditions every n/10 iterations
        if max(diag(D)) / min(diag(D)) > 1e7 % check condition number of D
            break;
        end
        if sum((options.sigma * options.B(:, end)).^2) > ((1.4 + 2 / (n + 1)) * n)^2 % check if maximum coordinate exceeds typical range
            break;
        end
        if max(fitness) - min(fitness) < 1e-10 % check function value range
            break;
        end
    end
    % generate new population
    pop = repmat(xmean, 1, options.lambda) + options.sigma * options.B * diag(randn(n, options.lambda) * diag(sqrt(options.eigenval))) * options.B' * randn(n, options.lambda);
    for i = 1:options.lambda
        fitness(i) = func(pop(:, i), varargin{:});
    end
    if min(fitness) < fopt % update best solution
        fopt = min(fitness);
        xopt = pop(:, fitness == fopt);
    end
    options.xmean = xmean; % update mean
    options.pc = pc; % update evolution path for C
    options.ps = ps; % update evolution path for sigma
    options.B = B; % update coordinate system matrix
    options.D = D; % update diagonal matrix of square roots of eigenvalues
end

其中，`func` 是指目标函数的函数句柄，`x0` 是指变量的初始值，`sigma0` 是指初始步长。`varargin` 是可选参数列表，将传递给目标函数。

函数的输出是 `xopt` 和 `fopt`，分别表示最佳解和最佳解的目标函数值。

该实现中使用了默认的参数，但是也可以通过在 `varargin` 中指定参数来更改默认值。