Описание задачи

Следующий эскиз адаптирован из Ray and Liew [2].

Этот эскиз представляет собой балку, приваренную к подложке. Балка поддерживает нагрузку P на расстоянии L от подложки. Балку сваривают на подложку верхним и нижним сварными швами, каждый из которых имеет длину l и толщину h. Балка имеет прямоугольное сечение, ширину b и высоту t. Материал балки стальной.

Двумя целями являются стоимость изготовления балки и отклонение конца балки под приложенной нагрузкой P. Нагрузка P фиксируется на уровне 6000 фунтов и расстояние L фиксируется на уровне 14 дюймов.

Переменные проекты:

Стоимость изготовления балки пропорциональна количеству материала в балке, $$(l+L)tb$$, плюс количество материала в сварных швах, $$l{h}^2$$. Используя константы пропорциональности из цитируемых статей, первой целью является

Отклонение луча пропорционально P и обратно пропорционально $$b{t}^3$$. Снова, используя константы пропорциональности из цитируемых документов, второй целью является

$$F_2(x)=\frac{P}{x_4{x}_3^3}C$$, где $$C=\frac{4(14{)}^3}{30\times 10^6}\approx 3.6587\times 10^{-4}$$ и $\mathit{P}=6,000$.

Задача имеет несколько ограничений.

Aineq = [1,0,0,-1];
bineq = 0;

Сводные данные, напряжение сдвига на сварных швах имеет ограничение $$\tau (x)$$ <= 13600.

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 решатель держит меньше половины своих решений на лучшем фронте Парето, использует в два раза больше точек, чем раньше.

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 принимает только тысячи.

Сравнение решений

The 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

The gamultiobj решение лучше в самом правом фрагменте графика, в то время как paretosearch решение лучше в самом левом фрагменте. paretosearch использует много меньше вычислений функции, чтобы получить ее решение.

Обычно, когда задача не имеет нелинейных ограничений, paretosearch по крайней мере, так же точно, как gamultiobj. Однако получившиеся наборы Парето могут иметь несколько различные области значений. В этом случае наличие нелинейного ограничения вызывает paretosearch решение быть менее точным по части области значений.

Одно из главных преимуществ paretosearch это то, что обычно требуется гораздо меньше вычислений функции.

Начните с однообъектных решений

Чтобы помочь решателям найти лучшие решения, запустите их с точек, которые являются решениями для минимизации отдельных целевых функций. The 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 Hybrid Function.

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:

Целевые и нелинейные ограничения

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] Deb, Kalyanmoy, J. Sundar, Udaya Bhaskara Rao N, and Shamik Chaudhuri. Основанная на контрольной точке многозначная оптимизация с использованием эволюционных алгоритмов. International Journal of Computational Intelligence Research, том 2, № 3, 2006, стр. 273-286. Доступно в https://www.softcomputing.net/ijcir/vol2-issu3-paper4.pdf

[2] Рэй, Т. и К. М. Лиев. Метафора роя для мультибъективной оптимизации проекта. Оптимизация проектирования 34, 2002, стр. 141-153.