В этом примере показано, как исследовать компромисс между силой и стоимостью луча. Несколько публикаций используют этот пример в качестве тестовой задачи для различных многоцелевых алгоритмов, включая Деб и др. [1] и Рэй и Лив [2].
Для видео обзора этого примера смотрите Множества Парето для Многоцелевой Оптимизации.
Следующий эскиз адаптируется от Рэя и Лива [2].
Этот эскиз представляет луч, сваренный на подложку. Луч поддерживает загрузку P на расстоянии L от подложки. Луч сварен на подложку с верхними и более низкими сварками, каждой длиной l и толщиной h. Луч имеет прямоугольное поперечное сечение, ширина b и высота t. Материал луча является сталью.
Эти две цели являются стоимостью производства луча и отклонением конца луча при прикладной загрузке P. Загрузка P фиксируется на уровне 6 000 фунтов, и расстояние L фиксируется на уровне 14 дюймов.
Переменные проекта:
x (1) = h, толщина сварок
x (2) = l, продолжительность сварок
x (3) = t, высота луча
x (4) = b, ширина луча
Стоимость производства луча пропорциональна на сумму материала в луче, , плюс сумма материала в сварках, . Используя коэффициенты пропорциональности из процитированных бумаг, первая цель
Отклонение луча пропорционально P и обратно пропорционально . Снова, с помощью коэффициентов пропорциональности из процитированных бумаг, вторая цель
, где и .
Проблема имеет несколько ограничений.
Толщина сварки не может превысить ширину луча. В символах, x (1) <= x (4). В синтаксисе тулбокса:
Aineq = [1,0,0,-1]; bineq = 0;
Напряжение сдвига на сварках не может превысить 13 600 фунтов на квадратный дюйм. Чтобы вычислить напряжение сдвига, сначала вычислите предварительные выражения:
Таким образом, напряжение сдвига на сварках имеет ограничение <= 13600.
Нормальное напряжение на сварках не может превысить 30 000 фунтов на квадратный дюйм. Нормальное напряжение .
Проваливающаяся грузоподъемность в вертикальном направлении должна превысить прикладную загрузку 6 000 фунтов. Используя значения модуля Янга psi и psi, проваливающееся ограничение загрузки . Численно, это становится неравенством .
Границы на переменных 0.125 <=x (1) <= 5, 0.1 <= x (2) <= 10, 0.1 <= x (3) <= 10, и 0.125 <= x (4) <= 5. В синтаксисе тулбокса:
lb = [0.125,0.1,0.1,0.125]; ub = [5,10,10,5];
Целевые функции появляются в конце этого примера в функциональном objval(x)
. Нелинейные ограничения появляются в конце этого примера в функциональном nonlcon(x)
.
paretosearch
РешениеМожно оптимизировать эту проблему несколькими способами:
Установите максимальное отклонение и найдите одно-объективную минимальную стоимость производства по проектам, которые удовлетворяют максимальному ограничению отклонения.
Установите максимальную стоимость производства и найдите одно-объективное минимальное отклонение по проектам, которые удовлетворяют ограничению стоимости производства.
Решите многоцелевую задачу, визуализируя компромисс между этими двумя целями.
Чтобы использовать многоцелевой подход, который дает больше информации о проблеме, устанавливает целевую функцию и нелинейную ограничительную функцию.
fun = @objval; nlcon = @nonlcon;
Решите задачу с помощью paretosearch
с 'psplotparetof'
функция plot. Чтобы уменьшать сумму диагностической информации об отображении, установите Display
опция к 'off'
.
opts_ps = optimoptions('paretosearch','Display','off','PlotFcn','psplotparetof'); rng default % For reproducibility [x_ps1,fval_ps1,~,psoutput1] = paretosearch(fun,4,Aineq,bineq,[],[],lb,ub,nlcon,opts_ps);
disp("Total Function Count: " + psoutput1.funccount);
Total Function Count: 1320
Для более сглаженной передней стороны Парето попытайтесь использовать больше точек.
npts = 160; % The default is 60
opts_ps.ParetoSetSize = npts;
[x_ps2,fval_ps2,~,psoutput2] = paretosearch(fun,4,Aineq,bineq,[],[],lb,ub,nlcon,opts_ps);
disp("Total Function Count: " + psoutput2.funccount);
Total Function Count: 4408
Это решение похоже на более сглаженную кривую, но оно имеет меньшую степень Объективных 2. Решатель принимает в три раза больше функциональных оценок при использовании 160 точек Парето вместо 60.
gamultiobj
РешениеЧтобы видеть, имеет ли решатель значение, попробуйте gamultiobj
решатель на проблеме. Установите эквивалентные опции как в предыдущем решении. Поскольку gamultiobj
решатель сохраняет меньше чем половину своих решений на лучшей передней стороне Парето, используйте в два раза больше точек, чем прежде.
opts_ga = optimoptions('gamultiobj','Display','off','PlotFcn','gaplotpareto','PopulationSize',2*npts); [x_ga1,fval_ga1,~,gaoutput1] = gamultiobj(fun,4,Aineq,bineq,[],[],lb,ub,nlcon,opts_ga);
disp("Total Function Count: " + gaoutput1.funccount);
Total Function Count: 33921
gamultiobj
берет десятки тысяч функциональных оценок, тогда как paretosearch
берет только тысячи.
gamultiobj
решение, кажется, отличается от paretosearch
решение, несмотря на то, что это затрудняет, чтобы сказать, потому что нанесенные на график шкалы отличаются. Постройте эти два решения на том же графике, с помощью той же шкалы.
fps2 = sortrows(fval_ps2,1,'ascend'); figure hold on plot(fps2(:,1),fps2(:,2),'r-') fga = sortrows(fval_ga1,1,'ascend'); plot(fga(:,1),fga(:,2),'b--') xlim([0,40]) ylim([0,1e-2]) legend('paretosearch','gamultiobj') xlabel 'Cost' ylabel 'Deflection' hold off
gamultiobj
решение лучше в самом правом фрагменте графика, тогда как paretosearch
решение лучше в крайнем левом фрагменте. paretosearch
использование много меньше функциональных оценок, чтобы получить его решение.
Как правило, когда проблема не имеет никаких нелинейных ограничений, paretosearch
по крайней мере, так же точно как gamultiobj
. Однако получившиеся Множества Парето могут иметь несколько различные области значений. В этом случае присутствие нелинейного ограничения вызывает paretosearch
решение быть менее точным по части области значений.
Одно из основных преимуществ paretosearch
это, обычно многим требуется меньше функциональных оценок.
Чтобы помочь решателям найти лучшие решения, запустите их с точек, которые являются решениями минимизации отдельных целевых функций. pickindex
функция возвращает одну цель в objval
функция. Используйте fmincon
найти одно-объективный optima. Затем используйте те решения в качестве начальных точек для многоцелевого поиска.
x0 = zeros(2,4); x0f = (lb + ub)/2; opts_fmc = optimoptions('fmincon','Display','off','MaxFunctionEvaluations',1e4); x0(1,:) = fmincon(@(x)pickindex(x,1),x0f,Aineq,bineq,[],[],lb,ub,@nonlcon,opts_fmc); x0(2,:) = fmincon(@(x)pickindex(x,2),x0f,Aineq,bineq,[],[],lb,ub,@nonlcon,opts_fmc);
Исследуйте одно-объективный optima.
objval(x0(1,:))
ans = 1×2
2.3810 0.0158
objval(x0(2,:))
ans = 1×2
76.7253 0.0004
Минимальная стоимость 2.381, который дает отклонение 0,158. Минимальное отклонение 0.0004, который имеет стоимость 76,7253. Нанесенные на график кривые довольно круты около концов их областей значений, означая, что вы получаете намного меньше отклонения, если вы берете стоимость немного выше ее минимума, или намного менее стоимости, если вы берете отклонение немного выше его минимума.
Запустите paretosearch
из одно-объективных решений. Поскольку вы построите решения позже тот же график, удалите paretosearch
функция plot.
opts_ps.InitialPoints = x0;
opts_ps.PlotFcn = [];
[x_psx0,fval_ps1x0,~,psoutput1x0] = paretosearch(fun,4,Aineq,bineq,[],[],lb,ub,nlcon,opts_ps);
disp("Total Function Count: " + psoutput1x0.funccount);
Total Function Count: 5024
Запустите ga
от тех же начальных точек, и удаляют его функцию построения графика.
opts_ga.InitialPopulationMatrix = x0;
opts_ga.PlotFcn = [];
[~,fval_ga,~,gaoutput] = gamultiobj(fun,4,Aineq,bineq,[],[],lb,ub,nlcon,opts_ga);
disp("Total Function Count: " + gaoutput.funccount);
Total Function Count: 46721
Постройте решения на тех же осях.
fps = sortrows(fval_ps1x0,1,'ascend'); figure hold on plot(fps(:,1),fps(:,2),'r-') fga = sortrows(fval_ga,1,'ascend'); plot(fga(:,1),fga(:,2),'b--') xlim([0,40]) ylim([0,1e-2]) legend('paretosearch','gamultiobj') xlabel 'Cost' ylabel 'Deflection' hold off
Путем запуска с одно-объективных решений, gamultiobj
решение немного лучше, чем paretosearch
решение в нанесенной на график области значений. Однако gamultiobj
берет почти в десять раз больше функциональных оценок, чтобы достигнуть его решения.
gamultiobj
может вызвать гибридный функциональный fgoalattain
автоматически попытаться достигнуть более точного решения. Смотрите, улучшает ли гибридная функция решение.
opts_ga.HybridFcn = 'fgoalattain'; [xgah,fval_gah,~,gaoutputh] = gamultiobj(fun,4,Aineq,bineq,[],[],lb,ub,nlcon,opts_ga); disp("Total Function Count: " + gaoutputh.funccount);
Total Function Count: 45512
fgah = sortrows(fval_gah,1,'ascend'); figure hold on plot(fps(:,1),fps(:,2),'r-') plot(fga(:,1),fga(:,2),'b--') plot(fgah(:,1),fgah(:,2),'g-') xlim([0,40]) ylim([0,1e-2]) legend('paretosearch','gamultiobj','gamultiobj/fgoalattain') xlabel 'Cost' ylabel 'Deflection' hold off
Гибридная функция обеспечивает небольшое улучшение на gamultiobj
решение, в основном в крайней левой части графика.
fgoalattain
Вручную от paretosearch
Точки решенияНесмотря на то, что paretosearch
не имеет никакой встроенной гибридной функции, можно улучшить paretosearch
решение путем выполнения fgoalattain
от paretosearch
точки решения. Создайте цель и веса для fgoalattain
при помощи той же настройки для fgoalattain
как описано в gamultiobj Гибридной Функции.
Fmax = max(fval_ps1x0); nobj = numel(Fmax); Fmin = min(fval_ps1x0); w = sum((Fmax - fval_ps1x0)./(1 + Fmax - Fmin),2); p = w.*((Fmax - fval_ps1x0)./(1 + Fmax - Fmin)); xnew = zeros(size(x_psx0)); nsol = size(xnew,1); fvalnew = zeros(nsol,nobj); opts_fg = optimoptions('fgoalattain','Display','off'); nfv = 0; for ii = 1:nsol [xnew(ii,:),fvalnew(ii,:),~,~,output] = fgoalattain(fun,x_psx0(ii,:),fval_ps1x0(ii,:),p(ii,:),... Aineq,bineq,[],[],lb,ub,nlcon,opts_fg); nfv = nfv + output.funcCount; end disp("fgoalattain Function Count: " + nfv)
fgoalattain Function Count: 13350
fnew = sortrows(fvalnew,1,'ascend'); figure hold on plot(fps(:,1),fps(:,2),'r-') plot(fga(:,1),fga(:,2),'b--') plot(fgah(:,1),fgah(:,2),'g-') plot(fnew(:,1),fnew(:,2),'k.-') xlim([0,40]) ylim([0,1e-2]) legend('paretosearch','gamultiobj','gamultiobj/fgoalattain','paretosearch/fgoalattain') xlabel 'Cost' ylabel 'Deflection'
Комбинация paretosearch
и fgoalattain
создает самую точную переднюю сторону Парето. Увеличьте масштаб, чтобы видеть.
xlim([3.64 13.69])
ylim([0.00121 0.00442])
hold off
Даже с дополнительным fgoalattain
расчеты, общая функция значит комбинацию, меньше, половина функции значат gamultiobj
одно только решение.
fprintf("Total function count for gamultiobj alone is %d.\n" + ... "For paretosearch and fgoalattain together it is %d.\n",... gaoutput.funccount,nfv + psoutput1x0.funccount)
Total function count for gamultiobj alone is 46721. For paretosearch and fgoalattain together it is 18374.
Нанесенные на график точки показывают оптимальные значения в функциональном пространстве. Чтобы определить, какие параметры достигают этих значений функции, найдите размер луча и размер сварки для того, чтобы понять конкретную мысль стоимости/отклонения. Например, график paretosearch
сопровождаемый fgoalattain
показывает точки со стоимостью приблизительно 6 и отклонением приблизительно 3.5e–3. Определите размеры луча и сварки, которые достигают этих точек.
whichgood = find(fvalnew(:,1) <= 6 & fvalnew(:,2) <= 3.5e-3); goodpoints = table(xnew(whichgood,:),fvalnew(whichgood,:),'VariableNames',{'Parameters' 'Objectives'})
goodpoints=4×2 table
Parameters Objectives
________________________________________ ___________________
0.63457 1.5187 10 0.67261 5.6973 0.0032637
0.61635 1.5708 10 0.63165 5.391 0.0034754
0.63228 1.5251 10 0.6674 5.6584 0.0032892
0.65076 1.4751 10 0.70999 5.976 0.0030919
Четыре набора параметров достигают стоимости меньше чем 6 и отклонения меньше, чем 3.5e–3:
Сварите толщину немного более чем 0,6
Сварите длину приблизительно 1,5
Излучите высоту 10 (верхняя граница)
Излучите ширину между 0,63 и 0.71
function [Cineq,Ceq] = nonlcon(x) sigma = 5.04e5 ./ (x(:,3).^2 .* x(:,4)); P_c = 64746.022*(1 - 0.028236*x(:,3)).*x(:,3).*x(:,4).^3; tp = 6e3./sqrt(2)./(x(:,1).*x(:,2)); tpp = 6e3./sqrt(2) .* (14+0.5*x(:,2)).*sqrt(0.25*(x(:,2).^2 + (x(:,1) + x(:,3)).^2)) ./ (x(:,1).*x(:,2).*(x(:,2).^2 / 12 + 0.25*(x(:,1) + x(:,3)).^2)); tau = sqrt(tp.^2 + tpp.^2 + (x(:,2).*tp.*tpp)./sqrt(0.25*(x(:,2).^2 + (x(:,1) + x(:,3)).^2))); Cineq = [tau - 13600,sigma - 3e4,6e3 - P_c]; Ceq = []; end function F = objval(x) f1 = 1.10471*x(:,1).^2.*x(:,2) + 0.04811*x(:,3).*x(:,4).*(14.0+x(:,2)); f2 = 2.1952./(x(:,3).^3 .* x(:,4)); F = [f1,f2]; end function z = pickindex(x,k) z = objval(x); % evaluate both objectives z = z(k); % return objective k end
[1] Деб, Kalyanmoy, Дж. Сундар, Удая Бхэскара Рао Н и Самик Каудури. Контрольная точка Основанная Многоцелевая Оптимизация Используя Эволюционные Алгоритмы. Международный журнал Вычислительного Интеллектуального Исследования, Издания 2, № 3, 2006, стр 273–286. Доступный в https://www.softcomputing.net/ijcir/vol2-issu3-paper4.pdf
[2] Излучите, T. и К. М. Лив. Метафора Роя для Многоцелевой Оптимизации Проекта. Техническая Оптимизация 34, 2002, pp.141–153.