💥1 概述
該算法旨在對評估成本高的黑盒函數進行全局多目标優化。例如，該算法已應用于同時優化生命周期評估（LCA）和化學過程模拟的成本[2]。但是，該算法也可以應用于其他黑盒功能，例如CFD模拟。它基于貝葉斯優化方法，該方法建構高斯過程代理模型以加速優化。此外，該算法可以在每次疊代中識别幾個有希望的點（批量順序模式）。這允許并行評估多個模拟。
📚2 運作結果
%% Initialize option structure
it = 1;
Opt = set_option_structure(Opt,X,Y);%% Write initial text file for TSEMO_log
create_log_file(X,Y,Opt,f,lb,ub)for i = 1:ceil(Opt.maxeval/Opt.NoOfBachSequential)
tic;
%% Scale Variables
[Xnew,Ynew] = ScaleVariables(X,Y,lb,ub,Opt) ;
%% Training of GP
for j = 1:Opt.Gen.NoOfGPs
Opt.GP(j).hyp = TrainingOfGP(Xnew,Ynew(:,j),Opt.GP(j));
end
%% Draw samples of GPs
for j = 1:Opt.Gen.NoOfGPs
Opt.Sample(j).f = posterior_sample(Xnew,Ynew(:,j),Opt.GP(j));
end
%% Determine pareto front of function samples
[Sample_pareto,Sample_xpareto,Sample_nadir] = Find_sample_pareto(Opt,i);
Opt.warmstart_pareto = Sample_xpareto;
%% Sampling point at maximum hypervolume improvement
[index,hv_imp] = hypervolume_improvement_index(Ynew,Sample_nadir,Sample_pareto,Opt);
xNew = Sample_xpareto(index,:);
Xnew = [Xnew;xNew];
for j = 1:Opt.Gen.NoOfInputDim
xnewtrue(:,j) = xNew(:,j)*(ub(j)-lb(j)) + lb(j);
end
for l = 1 : size(xnewtrue,1)
ytrue(l,:) = f(xnewtrue(l,:));
end
X = [X;xnewtrue];
Y = [Y;ytrue];
Ynew = zeros(size(Y,1),Opt.Gen.NoOfGPs);
for j = 1:Opt.Gen.NoOfGPs
Ynew(:,j) = (Y(:,j) - mean(Y(:,j)))/std(Y(:,j));
end
front = paretofront(Y);
Xpareto = X(front,:);
Ypareto = Y(front,:);
%% Update log each iterations
update_log_file(it,hv_imp,toc,xnewtrue,ytrue,Opt,Y,ub,lb)
%% Display
if it == 1
fprintf('%10s %10s %10s \n','Iteration', 'HypImp', 'Time(s)');
end
fprintf('%10d %10.4g %10.3g \n', it, hv_imp, toc);
it = it+1;
%% Determine final hyperparameter values for analysis
if i == ceil(Opt.maxeval/Opt.NoOfBachSequential)
for j = 1:Opt.Gen.NoOfGPs
Opt.GP(j).hyp = TrainingOfGP(Xnew,Ynew(:,j),Opt.GP(j));
end
hypf = zeros(Opt.Gen.NoOfInputDim+2,Opt.Gen.NoOfGPs);
for j = 1:Opt.Gen.NoOfGPs
covhyp = exp(Opt.GP(j).hyp.cov);
hypf(:,j) = [covhyp(1:Opt.Gen.NoOfInputDim).*(ub-lb)';covhyp(end)*std(Y(:,j));exp(Opt.GP(j).hyp.lik)*std(Y(:,j))];
end
%% Obtain Pareto front from spectral Gaussian process model
for j = 1:Opt.Gen.NoOfGPs
[Opt.Mean(j).f,Opt.Mean(j).varf] = mean_sample(Xnew,Ynew(:,j),Opt.GP(j));
end
[Mean_pareto,Mean_xpareto] = Find_mean_pareto(Opt);
for j = 1:Opt.Gen.NoOfInputDim
XParetoGP(:,j) = Mean_xpareto(:,j)*(ub(j)-lb(j)) + lb(j);
end
for j = 1:Opt.Gen.NoOfGPs
YParetoGP(:,j) = Mean_pareto(:,j)*std(Y(:,j)) + mean(Y(:,j));
end
for j = 1:Opt.Gen.NoOfGPs
for k = 1:Opt.pop
YParetoGPstd(k,j) = sqrt(Opt.Mean(j).varf(Mean_xpareto(k,:)))*std(Y(:,j));
end
end
%% Update log with final results
final_log_update(Xpareto,Ypareto,X,Y,XParetoGP,YParetoGP,hypf,Opt)
end
end
returnfunction create_log_file(X,Y,Opt,f,lb,ub)
try
function_name = func2str(f);
string1 = '';
string2 = {};
string3 = '';
string4 = '';
string5 = {};
string6 = '';
string7 = '';
for i = 1:size(X,2)
string1 = strcat(string1,'%8.4f ');
string2 = {string2{:},strcat('x',num2str(i))};
string3 = strcat(string3,'%+8s ');
end
for i = 1:size(Y,2)
string4 = strcat(string4,'%8.4f ');
string5 = {string5{:},strcat('f',num2str(i))};
string6 = strcat(string6,'%+8s ');
string7 = strcat(string7,'%8d ');
end
TSEMO_log = fopen( 'TSEMO_log.txt', 'w');
fprintf(TSEMO_log,'\n %s %s \n', 'TSEMO log file created on',date);
fprintf(TSEMO_log,'\n %s','This file shows the initial specifications of TSEMO and logs the output.');
fprintf(TSEMO_log,'\n %s \n', 'License information');
fprintf(TSEMO_log,'\n %s \n', 'BSD 2-Clause License');
fprintf(TSEMO_log,'\n %s', 'Copyright (c) 2017, Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin');
fprintf(TSEMO_log,'\n %s \n', 'All rights reserved.');
fprintf(TSEMO_log,'\n %s', 'Redistribution and use in source and binary forms, with or without');
fprintf(TSEMO_log,'\n %s \n', 'modification, are permitted provided that the following conditions are met:');
fprintf(TSEMO_log,'\n %s ', '*Redistributions of source code must retain the above copyright notice, this');
fprintf(TSEMO_log,'\n %s \n', ' list of conditions and the following disclaimer.');
fprintf(TSEMO_log,'\n %s ', '*Redistributions in binary form must reproduce the above copyright notice,');
fprintf(TSEMO_log,'\n %s ', ' this list of conditions and the following disclaimer in the
fprintf(TSEMO_log,strcat('%s',string7,'\n'), ' Number of spectral sampling points: ',Opt.GP(1:size(Y,2)).nSpectralpoints);
fprintf(TSEMO_log,strcat('%s',string7), ' Type of matern function: ',Opt.GP(1:size(Y,2)).matern);
fprintf(TSEMO_log,strcat('\n','%s',string7), ' Direct evaluations per input dimension: ',Opt.GP(1:size(Y,2)).fun_eval);
function update_log_file(it,hv_imp,toc,xnewtrue,ytrue,Opt,Y,ub,lb)
try
string1 = '';
string2 = {};
string3 = '';
string4 = '';
string5 = {};
string6 = '';
string7 = '';
string8 = {};
for i = 1:Opt.Gen.NoOfInputDim
string1 = strcat(string1,'%8.4f ');
string2 = {string2{:},strcat('x',num2str(i))};
string3 = strcat(string3,'%+8s ');
string8 = {string8{:},strcat('lambda',num2str(i))};
end
string8 = {string8{:},'sigmaf'};
string8 = {string8{:},'sigman'};
for i = 1:Opt.Gen.NoOfGPs
string4 = strcat(string4,'%8.4f ');
string5 = {string5{:},strcat('f',num2str(i))};
string6 = strcat(string6,'%+8s ');
string7 = strcat(string7,'%8d ');
hypcov = exp(Opt.GP(i).hyp.cov);
hypmat(:,i) = [hypcov(1:end-1).*(1./(ub-lb))';hypcov(end)*std(Y(:,i));exp(Opt.GP(i).hyp.lik)*std(Y(:,i))];
end