function [bestever,bestP,conv_curv,total_time]=SAEA_RFS_origfPB(lb,ub,ndims,npop,maxFEs,fobj,initPos,initFit) %% SAEA-RFS Adapted by Pedro Bento November 2023 (retrieved from the repository https://github.com/h-hg/DDEA/blob/master/README.md) %_____________start function_______________________________% warning off tic; flag = 'cubic'; cycle = 0; ts_num = 200; % the number of initial training samples K = 5; % old 20! the number of sub-problems sub_size = 20; %old 100! maximum dimension of all sub-problems nIter = 5; % manter! maximum iterations of sub-problem optimization lu = [lb*ones(1,ndims);ub*ones(1,ndims)]; % Initial population Parent = initPos; % generate initial archive FES = ts_num; %start the FEcounter with the number of started positions for the archive % the initial population goes to the archive, meaning initialize the remainder archive=initPos; init_l=ts_num-npop; XRRmin = repmat(lu(1, :), init_l, 1); XRRmax = repmat(lu(2, :), init_l ,1); archive = [archive; XRRmin + (XRRmax - XRRmin) .* lhsdesign(init_l ,ndims)]; % an extra column is used to store the fitness of each position stored in the archive archive(1:npop, ndims+1) = initFit; for i=npop+1:ts_num archive(i, ndims+1) = fobj(archive(i,1:ndims)); end bestever = min(archive(:,ndims+1)); bestP = inf; array=[linspace(1,ts_num,ts_num)', bestever.*ones(ts_num,1)]; % main loop while (FES < maxFEs) %% form sub-problem and training ts = archive; [lambda, gamma, dv , sp ] = RFS(sub_size,ts(:, 1:ndims),... ts(:, end), flag, K); %% Optimization [~, index] = min(archive(:, end)); bestX = archive(index, 1:ndims); sub_eval_min = []; for i = 1:K sub_Parent = Parent(:,dv{i}); % the population of sub-problem fit_Parent = RBF_eval(sub_Parent,ts(sp{i}, dv{i}), ... lambda{i}, gamma{i}, flag); subBound = lu(:, dv{i}); % use DE as a basic optimizer for each sub-problem for iter = 1:nIter sub_child = DE(sub_Parent, bestX(1,dv{i}), subBound); % Model evaluation fit_child = RBF_eval(sub_child, ts(sp{i}, dv{i}), ... lambda{i}, gamma{i}, flag); fit_pop = [fit_Parent;fit_child]; [~,y] = sort(fit_pop); com_pop = [sub_Parent;sub_child]; % Environment selection sub_Parent = com_pop(y(1:npop),:); fit_Parent = fit_pop(y(1:npop)); end % save the best solution and its estimated value of sub-problem sub_eval_min(i,1) = fit_Parent(1); sub_eval_min(i,2:length(dv{i})+1) = sub_Parent(1,:); % update the problem population Parent(:,dv{i}) = sub_Parent; end %% Determine the solution to be evaluated using real funcion solution = bestX; [min_model,sign] =min(sub_eval_min(:,1)); temp = RBF_eval(Parent(1,dv{sign}),ts(sp{sign},dv{sign}),lambda{sign},gamma{sign},flag); if temp <= min_model solution(:,dv{sign}) = Parent(1,dv{sign}); else solution(:,dv{sign}) = sub_eval_min(sign,2:length(dv{sign})+1); end %% search difference diff = setdiff(solution,archive(:,1:ndims),'rows'); for i=1:size(diff,1) diff(i,ndims+1) = fobj(diff(i,1:ndims)); end FES = FES + 1; archive = [archive;diff]; cycle = cycle+1; % update the global solution % bestever = min(min(archive(:,ndims+1)),bestever); %old [bestever_pot,index] = min(archive(:,ndims+1)); if (bestever_pot<bestever) bestever=bestever_pot; bestP=archive(index,ndims+1); else bestP=bestX; end array = [array;FES,bestever]; end total_time=toc; conv_curv=array(:,2); end %_____________end main function_______________________________% %% Secondary Functions %-------------------------------------------------------------------------- function [lambda,gamma,c,r]=RFS(sub_size,S,Y,flag,K) % %Input: %sub-size - maximum dimension of all sub-problems %S - Sample site matrix (m x d), m=number of sample sites, d=dimension %Y - Objective function values corresponding to points in S %flag - string determining what type of RBF model should be used %K - the number of sub-problems % %Output: %lambda, gamma - vectors with RBF model parameters %c - the desicion variables of each sub-problem %r - the samples used for training models of each sub-problem [n,m]=size(S); %% RF training i = 1; while i<=K M = randperm(m,randperm(sub_size,1)); % form sub-problem d = length(M); L = randperm(n,2*d); S1 = S(L,M); % training set Y1 = Y(L,:); [lambda{i},gamma{i}] = RBF(S1,Y1,flag); % train a model for each sub-problem c{i} = M; r{i} = L; if sum(isnan(lambda{i}))==0 i=i+1; end end end %function %-------------------------------------------------------------------------- function [lambda, gamma]=RBF(S,Y,flag) %function [lambda, gamma, rho]=RBF(S,Y,flag) % For Gaussian basis function %computes the parameters of the radial basis function interpolant %Copyright (c) 2012 by Juliane Mueller % % This file is part of the surrogate model module toolbox. % %Author information %Juliane Mueller %Tampere University of Technology, Finland %juliane.mueller2901@gmail.com % %Input: %S - Sample site matrix (m x d), m=number of sample sites, d=dimension %Y - Objective function values corresponding to points in S %flag - string determining what type of RBF model should be used % %Output: %lambda, gamma - vectors with RBF model parameters [m,n]=size(S); P=[S,ones(m,1)]; R=zeros(m,m); R = pdist2(S,S,'euclidean'); %rho=2; %Sd=real(sqrt(S.^2*ones(size(S'))+ones(size(S))*(S').^2-2*S*(S'))); % For Gaussian basis function %rho=max(max(Sd))/(n*m)^(1/n); % For Gaussian basis function if strcmp(flag,'cubic') %cubic RBF Phi= R.^3; elseif strcmp(flag,'TPS') %thin plate spline RBF R(R==0)=1; Phi=R.^2.*log(R); elseif strcmp(flag,'linear') %linear RBF Phi=R; elseif strcmp(flag,'Gaussian') Phi=exp(-R.^2/rho^2); elseif strcmp(flag,'multiquad') Phi=sqrt((R.^2 + rho^2)); elseif strcmp(flag,'invmultiquad') Phi=1./sqrt(R.^2 + rho^2); end A=[Phi,P;P', zeros(n+1,n+1)]; RHS=[Y;zeros(n+1,1)]; params=A\RHS; lambda=params(1:m); gamma=params(m+1:end); end %function %-------------------------------------------------------------------------- function Yest=RBF_eval(X,S,lambda,gamma,flag) %function Yest=RBF_eval(X,S,lambda,gamma,rho,flag) % For Gaussian basis function %Copyright (c) 2012 by Juliane Mueller % % This file is part of the surrogate model module toolbox. % %Author information %Juliane Mueller %Tampere University of Technology, Finland %juliane.mueller2901@gmail.com % %input: %X are points where function values should be calculated %S are points where the function values are known %lambda parameter vector %gamma contains the parameters of the optional polynomial tail %flag is a string indicating which RBF to be used %output: %the estimated function values at the points in X %rho=2; [mX,nX]=size(X); %dimensions of the points where function value should be predicted [mS,nS]=size(S); %dimesnions of sample sites taken so far if nX~=nS %check that both matrices are of the right shape X=X'; [mX,nX]=size(X); end R=zeros(mX,mS); %compute pairwise distances of points in X and S R = pdist2(X,S,'euclidean'); if strcmp(flag,'cubic') Phi= R.^3; elseif strcmp(flag,'TPS') R(R==0)=1; Phi=R.^2.*log(R); elseif strcmp(flag,'linear') Phi=R; elseif strcmp(flag,'Gaussian') Phi=exp(-R.^2/rho^2); elseif strcmp(flag,'multiquad') Phi=sqrt((R.^2 + rho^2)); elseif strcmp(flag,'invmultiquad') Phi=1./sqrt(R.^2 + rho^2); end Yest1=Phi*lambda; %first part of response surface Yest2=[X,ones(mX,1)]*gamma; %optional polynomial tail Yest=Yest1+Yest2; %predicted function value end %function %-------------------------------------------------------------------------- function U=DE(X,bestX,Bound) %变异向量用函数mutation(X,bestX,F,mutationStrategy) %交叉向量用函数crossover(X,V,CR,crossStrategy) %mutation %mutationStrategy=1：DE/rand/1, %mutationStrategy=2：DE/best/1, %mutationStrategy=3：DE/rand-to-best/1, %mutationStrategy=4：DE/best/2, %mutationStrategy=5：DE/rand/2. %crossover %crossStrategy=1:binomial crossover %crossStrategy=2:Exponential crossover F=0.8;%scaling factor 缩放因子 CR=1;%crossover rate 交叉概率 mutationStrategy=2;%变异策略 crossStrategy=1;%交叉策略 % while Generation<maxIteration V=mutation(X,bestX,F,mutationStrategy); U=crossover(X,V,CR,crossStrategy,Bound); % end end %-------------------------------------------------------------------------- function V=mutation(X,bestX,F,mutationStrategy) [NP,D]=size(X); for i=1:NP %在[1 NP]中产生nrandI个互不相等的随机数，且与i皆不相等 nrandI=5; r=randi([1,NP],1,nrandI); for j=1:nrandI equalr(j)=sum(r==r(j)); % 计数：产生的随机数中相同的个数 end equali=sum(r==i); % 行号与产生的随机数相同的次数 equalval=sum(equalr)+equali; % 相同的值 while(equalval>nrandI) r=randi([1,NP],1,nrandI); for j=1:nrandI equalr(j)=sum(r==r(j)); end equali=sum(r==i); equalval=sum(equalr)+equali; end switch mutationStrategy case 1 %mutationStrategy=1:DE/rand/1; V(i,:)=X(r(1),:)+F*(X(r(2),:)-X(r(3),:)); case 2 %mutationStrategy=2:DE/best/1; V(i,:)=bestX+F*(X(r(1),:)-X(r(2),:)); case 3 %mutationStrategy=3:DE/rand-to-best/1; V(i,:)=X(i,:)+F*(bestX-X(i,:))+F*(X(r(1),:)-X(r(2),:)); case 4 %mutationStrategy=4:DE/best/2; V(i,:)=bestX+F*(X(r(1),:)-X(r(2),:))+F*(X(r(3),:)-X(r(4),:)); case 5 %mutationStrategy=5:DE/rand/2; V(i,:)=X(r(1),:)+F*(X(r(2),:)-X(r(3),:))+F*(X(r(4),:)-X(r(5),:)); otherwise error('没有所指定的变异策略，请重新设定mutationStrategy的值'); end end end %-------------------------------------------------------------------------- function U=crossover(X,V,CR,crossStrategy,Bound) [NP,Dim]=size(X); switch crossStrategy %crossStrategy=1:binomial crossover case 1 for i=1:NP jRand=floor(rand*Dim);%由于jRand要在[0,1)*Dim中取值，故而用floor for j=1:Dim if rand<CR||j==jRand U(i,j)=V(i,j); else U(i,j)=X(i,j); end end end %crossStrategy=2:Exponential crossover case 2 for i=1:NP j=floor(rand*Dim);%由于j在[0,1)*Dim中取值，故而用floor L=0; U=X; U(i,j)=V(i,j); j=mod(j+1,D); L=L+1; while(rand<CR&&L<Dim) U(i,j)=V(i,j); j=mod(j+1,D); L=L+1; end end otherwise error('没有所指定的交叉策略，请重新设定crossStrategy的值'); end %% 越界重置 for i=1:NP for j=1:Dim while (U(i,j)>Bound(2,j)) || (U(i,j)<Bound(1, j)) U(i,j)=Bound(1,j)+rand*(Bound(2,j) - Bound(1,j)); end end end end %--------------------------------------------------------------------------