Документация

Описание проблемы

Следующий эскиз адаптируется от Рэя и Лива [2].

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

Эти две цели являются стоимостью производства луча и отклонением конца луча при прикладной загрузке P. Загрузка P фиксируется на уровне 6 000 фунтов, и расстояние 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(\mathrm{}14{)}^3}{\mathrm{}30\times 10^6}\approx 3.\mathrm{}\mathrm{}\mathrm{}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: 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.

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:

Объективные и нелинейные ограничения

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.