Сравнение 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. The 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);
The 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
