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

Следующий эскиз адаптируется от Рэя и Лива [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(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' функция 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:

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

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.