paretosearch и gamultiobjЭтот пример показывает, как создать набор точек на передней стороне Парето, использующей и paretosearch и gamultiobj. Целевая функция имеет две цели и двумерную контрольную переменную x. mymulti3 целевой функции доступен на вашем сеансе MATLAB®, когда вы нажимаете кнопку, чтобы отредактировать или попробовать этот пример. Также скопируйте код mymulti3 в свой сеанс. Для скорости вычисления векторизована функция.
type mymulti3function f = mymulti3(x) % f(:,1) = x(:,1).^4 + x(:,2).^4 + x(:,1).*x(:,2) - (x(:,1).^2).*(x(:,2).^2) - 9*x(:,1).^2; f(:,2) = x(:,2).^4 + x(:,1).^4 + x(:,1).*x(:,2) - (x(:,1).^2).*(x(:,2).^2) + 3*x(:,2).^3;
Найдите Множества Парето для целевых функций с помощью paretosearch и gamultiobj. Установите опцию UseVectorized на true для добавленной скорости. Включайте функцию построения графика, чтобы визуализировать Множество Парето.
rng default nvars = 2; opts = optimoptions(@gamultiobj,'UseVectorized',true,'PlotFcn','gaplotpareto'); [xga,fvalga,~,gaoutput] = gamultiobj(@(x)mymulti3(x),nvars,[],[],[],[],[],[],[],opts);
Optimization terminated: average change in the spread of Pareto solutions less than options.FunctionTolerance.
optsp = optimoptions('paretosearch','UseVectorized',true,'PlotFcn',{'psplotparetof' 'psplotparetox'}); [xp,fvalp,~,psoutput] = paretosearch(@(x)mymulti3(x),nvars,[],[],[],[],[],[],[],optsp);
Pareto set found that satisfies the constraints. Optimization completed because the relative change in the volume of the Pareto set is less than 'options.ParetoSetChangeTolerance' and constraints are satisfied to within 'options.ConstraintTolerance'.
Вычислите теоретически точные места на передней стороне Парето при помощи mymulti4. Функция mymulti4 доступна на вашем сеансе MATLAB®, когда вы нажимаете кнопку, чтобы отредактировать или попробовать этот пример.
type mymulti4function mout = mymulti4(x)
%
gg = [4*x(1)^3+x(2)-2*x(1)*(x(2)^2) - 18*x(1);
x(1)+4*x(2)^3-2*(x(1)^2)*x(2)];
gf = gg + [18*x(1);9*x(2)^2];
mout = gf(1)*gg(2) - gf(2)*gg(1);
Функция mymulti4 оценивает градиенты этих двух целевых функций. Затем, для области значений значений x(2), используйте fzero, чтобы определить местоположение точки, где градиенты точно параллельны, который является где вывод mout = 0.
a = [fzero(@(t)mymulti4([t,-3.15]),[2,3]),-3.15]; for jj = linspace(-3.125,-1.89,50) a = [a;[fzero(@(t)mymulti4([t,jj]),[2,3]),jj]]; end figure plot(fvalp(:,1),fvalp(:,2),'bo'); hold on fs = mymulti3(a); plot(fvalga(:,1),fvalga(:,2),'r*'); plot(fs(:,1),fs(:,2),'k.') legend('Paretosearch','Gamultiobj','True') xlabel('Objective 1') ylabel('Objective 2') hold off
gamultiobj находит точки с немного более широким распространением на пробеле целевой функции. Постройте решения на пробеле переменной решения, наряду с теоретической оптимальной кривой Парето и контурным графиком этих двух целевых функций.
[x,y] = meshgrid(1.9:.01:3.1,-3.2:.01:-1.8); mydata = mymulti3([x(:),y(:)]); myff = sqrt(mydata(:,1) + 39);% Spaces the contours better mygg = sqrt(mydata(:,2) + 28);% Spaces the contours better myff = reshape(myff,size(x)); mygg = reshape(mygg,size(x)); figure; hold on contour(x,y,mygg,50) contour(x,y,myff,50) plot(xp(:,1),xp(:,2),'bo') plot(xga(:,1),xga(:,2),'r*') plot(a(:,1),a(:,2),'-k') xlabel('x(1)') ylabel('x(2)') hold off
В отличие от решения paretosearch, решение gamultiobj имеет точки в экстремальных концах области значений на пробеле целевой функции. Однако решение paretosearch имеет больше точек, которые ближе к истинному решению и на пробеле целевой функции и на пробеле переменной решения. Число точек на передней стороне Парето отличается для каждого решателя, когда вы используете опции по умолчанию.
Что происходит, если решение вашей проблемы имеет контрольные переменные, которые являются большими? Исследуйте этот случай путем сдвига проблемных переменных. Для неограниченной проблемы может перестать работать gamultiobj, в то время как paretosearch более устойчив к таким сдвигам.
Для более легкого сравнения задайте 35 точек на передней стороне Парето для каждого решателя.
shift = [20,-30];
fun = @(x)mymulti3(x+shift);
opts.PopulationSize = 100; % opts.ParetoFraction = 35
[xgash,fvalgash,~,gashoutput] = gamultiobj(fun,nvars,[],[],[],[],[],[],opts);
Optimization terminated: average change in the spread of Pareto solutions less than options.FunctionTolerance.
gamultiobj не удается найти полезное Множество Парето.
optsp.PlotFcn = []; optsp.ParetoSetSize = 35; [xpsh,fvalpsh,~,pshoutput] = paretosearch(fun,nvars,[],[],[],[],[],[],[],optsp);
Pareto set found that satisfies the constraints. Optimization completed because the relative change in the volume of the Pareto set is less than 'options.ParetoSetChangeTolerance' and constraints are satisfied to within 'options.ConstraintTolerance'.
figure plot(fvalpsh(:,1),fvalpsh(:,2),'bo'); hold on plot(fs(:,1),fs(:,2),'k.') legend('Paretosearch','True') xlabel('Objective 1') ylabel('Objective 2') hold off
paretosearch находит распространение точек решения равномерно почти целой возможной областью значений.
Добавление границ, даже довольно свободных единиц, помогает и gamultiobj и paretosearch найти соответствующие решения. Установите нижние границы -50 в каждом компоненте и верхние границы 50.
opts.PlotFcn = []; optsp.PlotFcn = []; lb = [-50,-50]; ub = -lb; [xgash,fvalgash,~,gashoutput] = gamultiobj(fun,nvars,[],[],[],[],lb,ub,opts);
Optimization terminated: average change in the spread of Pareto solutions less than options.FunctionTolerance.
[xpsh2,fvalpsh2,~,pshoutput2] = paretosearch(fun,nvars,[],[],[],[],lb,ub,[],optsp);
Pareto set found that satisfies the constraints. Optimization completed because the relative change in the volume of the Pareto set is less than 'options.ParetoSetChangeTolerance' and constraints are satisfied to within 'options.ConstraintTolerance'.
figure plot(fvalpsh2(:,1),fvalpsh2(:,2),'bo'); hold on plot(fvalgash(:,1),fvalgash(:,2),'r*'); plot(fs(:,1),fs(:,2),'k.') legend('Paretosearch','Gamultiobj','True') xlabel('Objective 1') ylabel('Objective 2') hold off
В этом случае оба решателя находят хорошие решения.
paretosearch с решения gamultiobjПолучите подобную область значений решений от решателей путем запуска paretosearch с решения gamultiobj.
optsp.InitialPoints = xgash; [xpsh3,fvalpsh3,~,pshoutput3] = paretosearch(fun,nvars,[],[],[],[],lb,ub,[],optsp);
Pareto set found that satisfies the constraints. Optimization completed because the relative change in the volume of the Pareto set is less than 'options.ParetoSetChangeTolerance' and constraints are satisfied to within 'options.ConstraintTolerance'.
figure plot(fvalpsh3(:,1),fvalpsh3(:,2),'bo'); hold on plot(fvalgash(:,1),fvalgash(:,2),'r*'); plot(fs(:,1),fs(:,2),'k.') legend('Paretosearch','Gamultiobj','True') xlabel('Objective 1') ylabel('Objective 2') hold off
Теперь решение paretosearch подобно решению gamultiobj, несмотря на то, что некоторые точки решения отличаются.
Другой способ получить хорошее решение состоит в том, чтобы начать с точек, которые минимизируют каждую целевую функцию отдельно.
От многоцелевой функции создайте одно целевую функцию, которая выбирает каждую цель в свою очередь. Используйте переключенную функцию от предыдущего раздела. Поскольку вы даете точки хорошего начала решателям, вы не должны задавать границы.
nobj = 2; % Number of objectives x0 = -shift; % Initial point for single-objective minimization uncmin = cell(nobj,1); % Cell array to hold the single-objective minima allfuns = zeros(nobj,2); % Hold the objective function values eflag = zeros(nobj,1); fopts = optimoptions('patternsearch','Display','off'); % Use an appropriate solver here for i = 1:nobj indi = zeros(nobj,1); % Choose the objective to minimize indi(i) = 1; funi = @(x)dot(fun(x),indi); [uncmin{i},~,eflag(i)] = patternsearch(funi,x0,[],[],[],[],[],[],[],fopts); % Minimize objective i allfuns(i,:) = fun(uncmin{i}); end uncmin = cell2mat(uncmin); % Matrix of start points
Запустите paretosearch с одно-объективных минимальных точек и обратите внимание, что он имеет полный спектр в своих решениях. paretosearch добавляет, что случайная начальная буква указывает на предоставленные единицы в порядке иметь генеральную совокупность, по крайней мере, людей options.ParetoSetSize. Точно так же gamultiobj добавляет случайные точки в предоставленные единицы, чтобы получить генеральную совокупность, по крайней мере, людей (options.PopulationSize)*(options.ParetoFraction).
optsp = optimoptions(optsp,'InitialPoints',uncmin);
[xpinit,fvalpinit,~,outputpinit] = paretosearch(fun,nvars,[],[],[],[],[],[],[],optsp);Pareto set found that satisfies the constraints. Optimization completed because the relative change in the volume of the Pareto set is less than 'options.ParetoSetChangeTolerance' and constraints are satisfied to within 'options.ConstraintTolerance'.
Теперь решите проблему с помощью gamultiobj, начинающего с начальных точек.
opts = optimoptions(opts,'InitialPopulationMatrix',uncmin);
[xgash2,fvalgash2,~,gashoutput2] = gamultiobj(fun,nvars,[],[],[],[],[],[],opts);Optimization terminated: average change in the spread of Pareto solutions less than options.FunctionTolerance.
figure plot(fvalpinit(:,1),fvalpinit(:,2),'bo'); hold on plot(fvalgash2(:,1),fvalgash2(:,2),'r*'); plot(fs(:,1),fs(:,2),'k.') plot(allfuns(:,1),allfuns(:,2),'gs','MarkerSize',12) legend('Paretosearch','Gamultiobj','True','Start Points') xlabel('Objective 1') ylabel('Objective 2') hold off
Оба решателя заполняют переднюю сторону Парето между экстремальными точками с довольно точными и хорошо распределенными решениями.
Просмотрите конечные пункты на пробеле переменной решения.
figure; hold on xx = x - shift(1); yy = y - shift(2); contour(xx,yy,mygg,50) contour(xx,yy,myff,50) plot(xpinit(:,1),xpinit(:,2),'bo') plot(xgash2(:,1),xgash2(:,2),'r*') ashift = a - shift; plot(ashift(:,1),ashift(:,2),'-k') plot(uncmin(:,1),uncmin(:,2),'gs','MarkerSize',12); xlabel('x(1)') ylabel('x(2)') hold off
