В этом примере показано, как проанализировать компромисс между силой и стоимостью балки. В нескольких публикациях этот пример используется в качестве проблемы тестирования для различных мультиобъективных алгоритмов, включая Deb et al. [1] и Ray and Liew [2].
Видеообзор этого примера см. в разделе Аппараты Парето для многообъективной оптимизации.
Следующий эскиз адаптирован из Ray and Liew [2].
Этот эскиз представляет балку, приваренную к подложке. Луч поддерживает нагрузку Р на расстоянии L от подложки. Балку сваривают на подложку с верхним и нижним сварными швами, каждый длиной l и толщиной h. Балка имеет прямоугольное сечение, ширину b и высоту T. Материал балки - сталь.
Двумя целями являются стоимость изготовления балки и отклонение конца балки под приложенной нагрузкой P. Нагрузка P зафиксирована на уровне 6000 фунтов, а расстояние L зафиксировано на уровне 14 дюймов.
Конструктивные переменные:
x (1) = h, толщина сварных швов
x (2) = l, длина сварных швов
x (3) = t, высота балки
x (4) = b, ширина балки
Стоимость изготовления балки пропорциональна количеству материала в балке, tb, плюс количество материала в lh2. Используя константы пропорциональности из цитируемых документов, первой целью является
0,04811x3x4 (14 + x2).
Отклонение луча пропорционально P и обратно пропорционально . Опять же, используя константы пропорциональности из цитируемых документов, вторая цель
Px4x33C, 330×106≈3.6587×10-4 P = 6000.
Проблема имеет несколько ограничений.
Толщина сварного шва не может превышать ширину балки. В символах x (1) < = x (4). В синтаксисе панели инструментов:
Aineq = [1,0,0,-1]; bineq = 0;
Срезное напряжение (на сварных швах не может превышать 13600 psi. Чтобы рассчитать напряжение сдвига, сначала вычислите предварительные выражения:
+ x3) 2
(x1 + x3) 2
2τ1τ2x2R.
Таким образом, напряжение сдвига на сварных швах имеет ограничение (< = 13600.
Нормальное напряжение (на сварные швы не может превышать 30000 фунт/кв. дюйм. Нормальное напряжение P6Lx4x32≤30×103.
Нагрузка на изгиб в вертикальном направлении должна превышать приложенную нагрузку 6 000 фунт. Используя значения модуля 106 psi 12 × 106 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' функция графика. Для уменьшения объема диагностической информации установите 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: 1870
Для более гладкого фронта Парето попробуйте использовать больше точек.
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: 6254
Это решение выглядит как более гладкая кривая, но оно имеет меньшую степень цели 2. Решатель принимает в три раза больше оценок функций при использовании 160 точек Парето вместо 60.
gamultiobj РешениеЧтобы увидеть, имеет ли решатель значение, попробуйте gamultiobj решатель проблемы. Задайте эквивалентные параметры, как в предыдущем решении. Потому что gamultiobj решатель хранит меньше половины своих решений на лучшем фронте Pareto, используя в два раза больше точек, чем раньше.
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: 38401
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 найти оптимум с одной целью. Затем используйте эти решения в качестве начальных точек для многообъективного поиска.
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);
Изучите оптимум с одной целью.
objval(x0(1,:))
ans = 1×2
2.3810 0.0158
objval(x0(2,:))
ans = 1×2
76.7188 0.0004
Минимальная стоимость - 2,381, что даёт прогиб 0,158. Минимальный прогиб составляет 0,0004, что имеет стоимость 76,7253. Кривые на графике довольно крутые вблизи концов их диапазонов, что означает, что вы получаете гораздо меньше прогиба, если вы принимаете стоимость немного выше его минимума, или гораздо меньше, если вы берете прогиб немного выше его минимума.
Начать paretosearch из одноцелевых решений. Поскольку решения будут выведены на график позже на том же графике, удалите paretosearch функция графика.
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: 4839
Начать 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: 37441
Постройте графики решений на тех же осях.
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: 57478
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: 14049
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 37441. For paretosearch and fgoalattain together it is 18888.
Точки на графике показывают наилучшие значения в пространстве функций. Чтобы определить, какие параметры достигают этих значений функции, найдите размер балки и размер сварного шва, чтобы получить конкретную точку затрат/отклонения. Например, график 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.67262 5.6974 0.0032637
0.61635 1.5708 10 0.63165 5.391 0.0034753
0.63228 1.5251 10 0.6674 5.6584 0.0032892
0.65077 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] Деб, Калянмой, Дж. Сундар, Удая Бхаскара Рао Н и Шамик Чаудхури. Многоцелевая оптимизация на основе опорных точек с использованием эволюционных алгоритмов. Международный журнал исследований вычислительной разведки, том 2, № 3, 2006, стр. 273-286. Доступно по адресу https://www.softcomputing.net/ijcir/vol2-issu3-paper4.pdf
[2] Луч, T. и К. М. Лив. Метафора роя для многообъективной оптимизации конструкции. Инженерная оптимизация 34, 2002, стр.141-153.