CMAES的Matlab实现

% CMA-ES: Evolution Strategy with Covariance Matrix Adaptation for % nonlinear function minimization. % % This code is an excerpt from cmaes.m and implements the key parts % of the algorithm. It is intendend to be used for READING and % UNDERSTANDING the basic flow and all details of the CMA *algorithm*. % Computational efficiency is sometimes disregarded. % -------------------- Initialization -------------------------------- % User defined input parameters (need to be edited) clear all;clc; target_function = @(x)(x(1)-1).^2+0.3*sin(x(1))+(x(2)-1).^2+0.3*sin(x(2)); x_o = [3:-0.01:0]; y_o = zeros(length(x_o),length(x_o)); for index_i = 1:length(x_o) for index_j = 1:length(x_o) y_o(index_i,index_j) = target_function([x_o(index_i) x_o(index_j)]); end end subplot(1,2,1); mesh(x_o,x_o,y_o); subplot(1,2,2); mesh(x_o,x_o,y_o); % subplot(2,2,2); % contourf(x_o,x_o,y_o,20); set(gcf,'outerposition',get(0,'screensize')); noise = 0.1; N = 2; % number of objective variables/problem dimension xmean = [2.5;2.5]; % objective variables initial point sigma = 0.05; % coordinate wise standard deviation (step-size) stopeval = 20; % Strategy parameter setting: Selection lambda = 4+floor(3*log(N)); % population size, offspring number lambda = 20; mu = lambda/2; % lambda=12; mu=3; weights = ones(mu,1); would be (3_I,12)-ES weights = log(mu+1/2)-log(1:mu)'; % muXone recombination weights mu = floor(mu); % number of parents/points for recombination weights = weights/sum(weights); % normalize recombination weights array mueff=sum(weights)^2/sum(weights.^2); % variance-effective size of mu % Strategy parameter setting: Adaptation cc = (4+mueff/N) / (N+4 + 2*mueff/N); % time constant for cumulation for C cs = (mueff+2)/(N+mueff+5); % t-const for cumulation for sigma control c1 = 2 / ((N+1.3)^2+mueff); % learning rate for rank-one update of C cmu = 2 * (mueff-2+1/mueff) / ((N+2)^2+2*mueff/2); % and for rank-mu update damps = 1 + 2*max(0, sqrt((mueff-1)/(N+1))-1) + cs; % damping for sigma % Initialize dynamic (internal) strategy parameters and constants pc = zeros(N,1); ps = zeros(N,1); % evolution paths for C and sigma B = eye(N); % B defines the coordinate system D = eye(N); % diagonal matrix D defines the scaling C = B*D*(B*D)'; % covariance matrix eigeneval = 0; % B and D updated at counteval == 0 chiN=N^0.5*(1-1/(4*N)+1/(21*N^2)); % expectation of Opt_Record = []; n=0; pause(1); % -------------------- Generation Loop -------------------------------- counteval = 0; % the next 40 lines contain the 20 lines of interesting code while counteval < stopeval % Generate and evaluate lambda offspring for k=1:lambda arz(:,k) = randn(N,1); % standard normally distributed vector arx(:,k) = xmean + sigma * (B*D * arz(:,k)); % add mutation % Eq. 40 arfitness(k) = target_function(arx(:,k)) + noise*randn; % objective function call end %% 绘制响应面和采集函数 subplot(1,2,2); mesh(x_o,x_o,y_o); hold on plot3(arx(2,:),arx(1,:),arfitness,'k.','MarkerSize',20); plot3(xmean(2),xmean(1),target_function(xmean),'r.','MarkerSize',30); axis([0 3 0 3 0 9]) title(['Iteration = ' num2str(counteval+1)]) hold off n = n+1; saveas(gcf,['baysian_' num2str(n) '.bmp']) % Sort by fitness and compute weighted mean into xmean [arfitnesss, arindex] = sort(arfitness); % minimization xmean = arx(:,arindex(1:mu))*weights; % recombination % Eq. 42 zmean = arz(:,arindex(1:mu))*weights; % == D?-1*B’*(xmean-xold)/sigma % Cumulation: Update evolution paths ps = (1-cs)*ps + (sqrt(cs*(2-cs)*mueff)) * (B * zmean); % Eq. 43 hsig = norm(ps)/sqrt(1-(1-cs)^(2*counteval/lambda))/chiN < 1.4+2/(N+1); pc = (1-cc)*pc + hsig * sqrt(cc*(2-cc)*mueff) * (B*D*zmean); % Eq. 45 % Adapt covariance matrix C C = (1-c1-cmu) * C ... % regard old matrix % Eq. 47 + c1 * (pc*pc' ... % plus rank one update + (1-hsig) * cc*(2-cc) * C) ... % minor correction + cmu ... % plus rank mu update * (B*D*arz(:,arindex(1:mu))) ... * diag(weights) * (B*D*arz(:,arindex(1:mu)))'; % Adapt step-size sigma sigma = sigma * exp((cs/damps)*(norm(ps)/chiN - 1)); % Eq. 44 % Update B and D from C C=triu(C)+triu(C,1)'; % enforce symmetry [B,D] = eig(C); % eigen decomposition, B==normalized eigenvectors D = diag(sqrt(diag(D))); % D contains standard deviations now Opt_Record = [Opt_Record xmean]; counteval = counteval + 1; pause(1.1); end % while, end generation loop