Задача минимальной высоты стека

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

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

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

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 концентрацию $0.000125g{m}^{-3}$. Это - установленный законом максимум, при котором концентрация диоксида серы не должна превышать концентрацию в области R. На графике ясно видно, что концентрация SO2 превышает максимальную для начальной высоты стека.

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

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

Как fseminf Работы

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

таким, что

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

$Aeq*x=beq$ (Линейные ограничения равенства)

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

$ceq(x)=0$ (Нелинейные Ограничения Равенства)

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

и

$K_j(x,w)<=0$, где $w\in I_j$ для $j=1,...,n_{inf}$ (Нелинейные полунепрерывные ограничения)

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

Функция fseminf решает общую полубесконечную задачу, начиная с начального значения, $x_0$и использование итерационной процедуры для получения оптимального решения, $x_{opt}$.

Ключевым компонентом алгоритма является обработка «полубесконечных» ограничений, $K_j$. В $x_{opt}$ требуется, чтобы $K_j$ должны быть допустимыми при каждом значении $w$ в интервале $I_j$. Это ограничение может быть упрощено, принимая во внимание все локальные максимумы $K_j$ по отношению к $w$ в интервале $I_j$. Исходное ограничение эквивалентно требованию, чтобы значение $K_j$ на каждом из вышеуказанных локальных максимумов возможно.

fseminf вычисляет приближение ко всем локальным максимальным значениям каждого полубесконечного ограничения, $K_j$. Для этого fseminf сначала вычисляет каждое полубесконечное ограничение по mesh из $w$ значения. Простая схема дифференцирования затем используется, чтобы вычислить все локальные максимальные значения $K_j$ из вычисленного полуинфинитного ограничения.

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

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

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

Прежде чем мы сможем позвонить fseminf чтобы решить задачу, нам нужно написать функцию, чтобы вычислить нелинейные ограничения в этой задаче. Необходимо реализовать ограничение, заключающееся в том, что концентрация диоксида серы на уровне земли не должна превышать $0.000125g{m}^{-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 оценивает «полунепрерывные» ограничения на mesh как часть общего вычисления этих ограничений. Когда ваша функция ограничения вызывается fseminf, интервал mesh, который вы должны использовать, передается вашей функции. Fseminf первоначально вызовет функцию ограничения с интервалом между сетками, s, установлено на NaN. Это позволяет вам инициализировать размер сетки для оценки ограничений. Здесь у нас есть одно «бесконечное» ограничение в двух «бесконечных» переменных. Это означает, что нам нужно инициализировать размер сетки в матрицу 1 на 2, в этом случае s = [1000 4000].

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

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

3. Вычислим ограничения над mesh

После определения mesh ограничения можно вычислить по ней. Затем эти ограничения возвращаются к 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 в этой точке была ниже $0.000125g{m}^{-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. Определите mesh, которая будет использоваться для оценки ограничений

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

3. Вычислим ограничения над mesh

Остальная часть кода вычисляет концентрацию SO2 в каждой точке mesh скорость/направление ветра. Затем эти ограничения возвращаются к 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™ пользователя в главе Использование решателей Optimization Toolbox в разделе Ограниченная нелинейная оптимизация: fseminf Формулировка задачи и алгоритм.