Проблема минимальной высоты стека

Рассмотрим область 20 км на 20 км, R, в которой должны быть размещены десять дымовых труб. Эти дымовые трубы выделяют в атмосферу несколько загрязняющих веществ, одним из которых является диоксид серы. Расположения штабелей по оси X и Y являются фиксированными, но высота штабелей может изменяться.

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

Сначала постройте график дымовых труб на их начальной высоте. Обратите внимание, что мы увеличили масштаб субрегиона R 4km на 4km, который содержит стопки дымохода.

h0 = [210;210;180;180;150;150;120;120;90;90];
plotChimneyStacks(h0, 'Chimney Stack Initial Height');

Эта проблема связана с двумя параметрами окружающей среды: скоростью и направлением ветра. Позже в этом примере мы позволим этим параметрам изменяться, но для первой проблемы мы установим для этих параметров типовые значения.

% Wind direction in radians
theta0 = 3.996;   
% Wind speed in m/s
U0 = 5.64;

Теперь давайте построим график концентрации диоксида серы (SO2) на уровне земли по всей области R (помните, что график дымовых труб был на меньшей области). Концентрация SO2 была рассчитана с заданными исходными высотами дымовых труб.

Мы видим, что концентрация SO2 варьируется по интересующей области. Заметим две особенности графика диоксида серы:

В очень простых терминах первый признак обусловлен преобладающим ветром, который дует SO2 к верхнему левому углу плоскости (x, y) в этом примере. Второй фактор обусловлен тем, что SO2 транспортируется на землю посредством диффузии. Это более медленный процесс по сравнению с преобладающим ветром и, таким образом, SO2 достигает уровня земли только в верхнем левом углу интересующей области.

Для более подробного обсуждения атмосферной дисперсии из дымовых труб см. ссылку, приведенную во введении.

Розовая плоскость указывает концентрацию SO2 $\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{}^{\mathrm{0,000125гм-3}}$. Это установленный максимум, для которого концентрация диоксида серы не должна превышать в области R. Из графика ясно видно, что концентрация SO2 превышает максимум для начальной высоты дымовой трубы.

Проверка файла MATLAB concSulfurDioxide чтобы увидеть, как рассчитывается концентрация диоксида серы.

plotSulfurDioxide(h0, theta0, U0, ...
    'Sulfur Dioxide Concentration at Initial Stack Height');

Как fseminf Работы

Прежде чем решить проблему минимальной высоты стека, мы наметим, как fseminf решает полулегкую задачу. Общая задача полубезграничного программирования может быть заявлена как:

такой, что

$\mathrm{Ax}<$= b (Линейные ограничения неравенства)

$\mathrm{Aeq}*x=$ beq (ограничения линейного равенства)

$c(x)\mathrm{<}$ = 0 (нелинейные ограничения неравенства)

$\mathrm{ceq}(x)\mathrm{}$= 0 (ограничения нелинейного равенства)

$l<=x$ < = u (Границы)

и

${}_{\mathrm{Kj}}(x,w)\mathrm{}$< = 0, ${}_{\mathrm{где}}$ w∈Ij $\mathrm{}дляj={}_{1,.}$.., ninf (Нелинейные полубесконечные ограничения)

Этот алгоритм позволяет задать ограничения для нелинейной задачи оптимизации, которые должны удовлетворяться в интервалах вспомогательной переменной. $$Обратите внимание, что для fseminf, эта переменная ограничена 1 или 2 мерными для каждого полупредельного ограничения.

Функция fseminf решает общую полубесконечную задачу, начиная с начального значения ${}_{\mathrm{x0}}$и используя итерационную процедуру для получения оптимального решения ${}_{\mathrm{xopt}}$.

Ключевая составляющая алгоритма - обработка «полубесконечных» ограничений, ${}_{\mathrm{Kj}}$. При ${}_{\mathrm{xopt}}$ требуется, чтобы ${}_{\mathrm{Kj}}$ был осуществим при каждом значении $w$ в интервале ${}_{\mathrm{Ij}}$. Это ограничение может быть упрощено путем учета всех локальных максимумов ${}_{\mathrm{Kj}}$ относительно $w$ в интервале ${}_{\mathrm{Ij}}$. Исходное ограничение эквивалентно требованию, чтобы значение ${}_{\mathrm{Kj}}$ в каждом из вышеупомянутых локальных максимумов было осуществимым.

fseminf вычисляет приближение ко всем локальным максимальным значениям каждого полубесконечного ограничения, ${}_{\mathrm{Kj}}$. Для этого fseminf сначала вычисляет каждое полулегкое ограничение по сетке из $$значений w. Затем используется простая разностная схема для вычисления всех локальных максимальных значений ${}_{\mathrm{Kj}}$ из вычисленного полубесконечного ограничения.

Как будет показано ниже, эта сетка создается в функции ограничения. Интервал, который следует использовать для каждой $$координаты w сетки, предоставляется функции ограничения fseminf.

При каждой итерации алгоритма выполняются следующие шаги:

Запись функции нелинейных ограничений

Прежде чем мы сможем позвонить fseminf для решения задачи необходимо записать функцию для оценки нелинейных ограничений в этой задаче. Ограничение, которое должно быть реализовано, состоит в том, что концентрация приземного диоксида серы не должна превышать $\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{}^{\mathrm{0,000125гм-3}}$ в каждой точке области R.

Это полулегкое ограничение, и реализация функции ограничения объясняется в этом разделе. Для проблемы минимальной высоты стека мы реализовали ограничение в файле MATLAB airPollutionCon.

type airPollutionCon.m

function [c, ceq, K, s] = airPollutionCon(h, s, theta, U)
%AIRPOLLUTIONCON Constraint function for air pollution demo
% 
%   [C, CEQ, K, S] = AIRPOLLUTIONCON(H, S, THETA, U) calculates the
%   constraints for the air pollution Optimization Toolbox (TM) demo. This
%   function first creates a grid of (X, Y) points using the supplied grid
%   spacing, S. The following constraint is then calculated over each point
%   of the grid:
%
%   Sulfur Dioxide concentration at the specified wind direction, THETA and
%   wind speed U <= 1.25e-4 g/m^3
%
%   See also AIRPOLLUTION

%   Copyright 2008 The MathWorks, Inc.

% Initial sampling interval
if nargin < 2 || isnan(s(1,1))
    s = [1000 4000]; 
end

% Define the grid that the "infinite" constraints will be evaluated over
w1x = -20000:s(1,1):20000;
w1y = -20000:s(1,2):20000;
[t1,t2] = meshgrid(w1x,w1y);

% Maximum allowed sulphur dioxide
maxsul = 1.25e-4; 

% Calculate the constraint over the grid
K = concSulfurDioxide(t1, t2, h, theta, U) - maxsul;

% Rescale constraint to make it 0(1)
K = 1e4*K;

% No finite constraints
c = [];
ceq = [];

Эта функция иллюстрирует общую структуру функции ограничения для задачи полулегкого программирования. В частности, функция ограничения для fseminf может быть разбит на три части:

1. Определение начального размера сетки для оценки зависимости

Напомним, что fseminf оценивает «полубесконечные» ограничения по сетке как часть общего расчета этих ограничений. Когда функция ограничения вызывается fseminf, интервал между сетками, который следует использовать, предоставляется функции. Fseminf сначала вызовет функцию ограничения с интервалом сетки, s, установите значение NaN. Это позволяет инициализировать размер сетки для вычисления ограничения. Здесь мы имеем одно «бесконечное» ограничение в двух «бесконечных» переменных. Это означает, что нам нужно инициализировать размер сетки для матрицы 1 на 2, в данном случае, s = [1000 4000].

2. Определение сетки, которая будет использоваться для оценки зависимости

Необходимо создать сетку, которая будет использоваться для оценки ограничения. Три строки кода после комментария «Определить сетку, по которой будут вычисляться» бесконечные «ограничения» в airPollutionCon может быть модифицирован для большинства 2-d задач полусреднего программирования.

3. Расчет ограничений по сетке

После определения сетки над ней можно рассчитать ограничения. Эти ограничения затем возвращаются в fseminf из вышеприведенной функции ограничения.

Заметим, что в этой проблеме мы также масштабировали ограничения так, чтобы они варьировались в масштабе, который ближе к масштабу целевой функции. Это помогает fseminf во избежание проблем масштабирования, связанных с целями и ограничениями, которые различаются в различных масштабах.

Решение проблемы оптимизации

Теперь мы можем позвонить fseminf для решения проблемы. Все дымовые трубы должны иметь высоту не менее 10 м, и мы используем начальную высоту трубы, указанную ранее. Обратите внимание, что третий входной аргумент для fseminf ниже (1) указывает, что существует только одно полулегкое ограничение.

lb = 10*ones(size(h0));
[hsopt, sumh, exitflag] = fseminf(@(h)sum(h), h0, 1, ...
    @(h,s) airPollutionCon(h,s,theta0,U0), [], [], [], [], lb);

Local minimum possible. Constraints satisfied.

fseminf stopped because the predicted change in the objective function
is less than the value of the function tolerance and constraints 
are satisfied to within the value of the constraint tolerance.

fprintf('\nMinimum computed cumulative height of chimney stacks : %7.2f m\n', sumh);

Minimum computed cumulative height of chimney stacks : 3667.19 m

Минимальная совокупная высота, вычисленная по fseminf значительно выше начальной общей высоты дымовых труб. Ниже в примере будет показано, как изменяется минимальная кумулятивная высота при добавлении к проблеме неопределенности параметров. Пока давайте построим графики дымовых труб на их оптимальной высоте.

Проверка файла MATLAB plotChimneyStacks для просмотра способа создания графика.

plotChimneyStacks(hsopt, 'Chimney Stack Optimal Height');

Проверка результатов оптимизации

Напомним, что fseminf определяет, что полусредненное ограничение выполняется везде, гарантируя, что дискретизированные максимумы ограничения находятся ниже заданной границы. Мы можем убедиться в том, что полулегкое ограничение выполняется везде, построив график концентрации диоксида серы на уровне земли для оптимальной высоты пакета.

Следует отметить, что концентрация диоксида серы принимает максимально возможное значение в левом верхнем углу плоскости (x, y), т.е. при x = -20000m, y = 20000m. Эта точка отмечена синей точкой на рисунке ниже и подтверждена путем вычисления концентрации диоксида серы в этой точке.

Проверка файла MATLAB plotSulfurDioxide чтобы увидеть, как создавались графики.

titleStr = 'Optimal Sulfur Dioxide Concentration and its maximum (blue)';
xMaxSD = [-20000 20000];
plotSulfurDioxide(hsopt, theta0, U0, titleStr, xMaxSD);

SO2Max = concSulfurDioxide(-20000, 20000, hsopt, theta0, U0);
fprintf('Sulfur Dioxide Concentration at x = -20000m, y = 20000m : %e g/m^3\n', SO2Max);

Sulfur Dioxide Concentration at x = -20000m, y = 20000m : 1.250000e-04 g/m^3

Учет неопределенности в факторах окружающей среды

Концентрация диоксида серы зависит от нескольких факторов окружающей среды, которые были сохранены при фиксированных значениях в вышеупомянутой проблеме. Двумя факторами окружающей среды являются скорость ветра и направление ветра. Более подробное обсуждение всех параметров проблемы см. в ссылке, приведенной во введении.

Мы можем исследовать изменение поведения системы в отношении скорости и направления ветра. В этом разделе примера мы хотим убедиться, что пределы диоксида серы удовлетворяются, даже если направление ветра изменяется от 3,82 рад до 4,18 рад, а средняя скорость ветра изменяется от 5 до 6,2 м/с.

Для обеспечения того, чтобы концентрация диоксида серы не превышала предела в области R. Это ограничение должно быть выполнимым для всех пар скорости и направления ветра.

Такое ограничение будет иметь четыре «бесконечные» переменные (скорость и направление ветра и координаты x-y земли). Однако любое полупредельное ограничение, предоставляемое fseminf может иметь не более двух «бесконечных» переменных.

Для реализации этого ограничения в подходящей форме для fseminf, мы напоминаем концентрацию SO2 при оптимальной высоте стопки в предыдущей проблеме. В частности, концентрация SO2 принимает максимально возможное значение при x = -20000m, y = 20000m. Чтобы уменьшить число «бесконечных» переменных, мы будем считать, что концентрация SO2 также примет свое максимальное значение в этой точке, когда присутствует неопределенность. Затем мы требуем, чтобы концентрация SO2 в этой точке была ниже $\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{}^{\mathrm{0,000125gm-3}}$ для всех пар скорости и направления ветра.

Это означает, что «бесконечными» переменными для этой задачи являются скорость и направление ветра. Чтобы увидеть, как было реализовано это ограничение, проверьте файл MATLAB. uncertainAirPollutionCon.

type uncertainAirPollutionCon.m

function [c, ceq, K, s] = uncertainAirPollutionCon(h, s)
%UNCERTAINAIRPOLLUTIONCON Constraint function for air pollution demo
% 
%   [C, CEQ, K, S] = UNCERTAINAIRPOLLUTIONCON(H, S) calculates the
%   constraints for the fseminf Optimization Toolbox (TM) demo. This
%   function first creates a grid of wind speed/direction points using the
%   supplied grid spacing, S. The following constraint is then calculated
%   over each point of the grid:
%
%   Sulfur Dioxide concentration at x = -20000m, y = 20000m <= 1.25e-4
%   g/m^3
%
%   See also AIRPOLLUTIONCON, AIRPOLLUTION

%   Copyright 2008 The MathWorks, Inc.

% Maximum allowed sulphur dioxide
maxsul = 1.25e-4; 

% Initial sampling interval
if nargin < 2 || isnan(s(1,1))
    s = [0.02 0.04]; 
end

% Define the grid that the "infinite" constraints will be evaluated over
w1x = 3.82:s(1,1):4.18; % Wind direction
w1y = 5.0:s(1,2):6.2;   % Wind speed
[t1,t2] = meshgrid(w1x,w1y);

% We assume the maximum SO2 concentration is at [x, y] = [-20000, 20000]
% for all wind speed/direction pairs. We evaluate the SO2 constraint over
% the [theta, U] grid at this point.
K = concSulfurDioxide(-20000, 20000, h, t1, t2) - maxsul;
 
% Rescale constraint to make it 0(1)
K = 1e4*K;

% No finite constraints
c = [];
ceq = [];

Эту функцию ограничения можно разделить на те же три секции, что и ранее:

1. Определение начального размера сетки для оценки зависимости

Код, следующий за комментарием «Начальный интервал выборки», инициализирует размер сетки.

2. Определение сетки, которая будет использоваться для оценки зависимости

Следующий раздел кода создает сетку (теперь по скорости и направлению ветра), используя конструкцию, аналогичную той, которая использовалась в начальной задаче.

3. Расчет ограничений по сетке

Остальная часть кода вычисляет концентрацию SO2 в каждой точке сетки скорость/направление ветра. Эти ограничения затем возвращаются в fseminf из вышеприведенной функции ограничения.

Теперь мы можем позвонить fseminf для решения проблемы высоты стека с учетом неопределенности факторов окружающей среды.

[hsopt2, sumh2, exitflag2] = fseminf(@(h)sum(h), h0, 1, ...
    @uncertainAirPollutionCon, [], [], [], [], lb);

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the  value of the constraint tolerance.

fprintf('\nMinimal computed cumulative height of chimney stacks with uncertainty: %7.2f m\n', sumh2);

Minimal computed cumulative height of chimney stacks with uncertainty: 3812.14 m

Теперь мы можем посмотреть на разницу между минимальной вычисленной совокупной высотой стека для проблемы с неопределенностью параметров и без нее. Вы должны видеть, что минимальная кумулятивная высота увеличивается при добавлении неопределенности к проблеме. Это ожидаемое увеличение высоты позволяет SO2 концентрации оставаться ниже установленного максимума для всех пар скорость/направление ветра в указанном диапазоне.

Мы можем проверить, чтобы концентрация диоксида серы не превышала предела по интересующей области посредством осмотра участка диоксида серы. Для данной (x, y) точки строим график максимальной концентрации SO2 для скорости и направления ветра в заданных диапазонах. Обратите внимание, что мы увеличили масштаб в левом верхнем углу плоскости X-Y.

titleStr = 'Optimal Sulfur Dioxide Concentration under Uncertainty';
thetaRange = 3.82:0.02:4.18;
URange = 5:0.2:6.2;
XRange = [-20000,-15000];
YRange = [15000,20000];
plotSulfurDioxideUncertain(hsopt2, thetaRange, URange, XRange, YRange, titleStr);

Наконец, мы планируем дымовые трубы на их оптимальной высоте, когда существует неопределенность в определении проблемы.

plotChimneyStacks(hsopt2, 'Chimney Stack Optimal Height under Uncertainty');

Существует множество вариантов алгоритма полунейминитального программирования, fseminf. Дополнительные сведения см. в Руководстве пользователя Optimization Toolbox™ в главе Использование решателей панели инструментов оптимизации в разделе Ограниченная нелинейная оптимизация: формулирование и алгоритм проблем fseminf.