Минимальная проблема высоты трубы

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

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

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

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) по всему региону Р (помните, что график стеков дымохода был по меньшей области). Концентрация SO2 была вычислена с набором стеков дымохода к их начальным высотам.

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

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

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

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

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

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

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

Прежде чем мы сможем вызвать fseminf чтобы решить задачу, мы должны записать функцию, чтобы оценить нелинейные ограничения в этой проблеме. Ограничение, которое будет реализовано, состоит в том, что концентрация Двуокиси серы уровня земли не должна превышать $0.000125g{m}^{-3}$ в каждой точке в области Р.

Это - полубесконечное ограничение, и реализация ограничительной функции объяснена в этом разделе. Для минимальной проблемы высоты трубы мы реализовали ограничение в файле 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. Задайте начальный размер mesh для ограничительной оценки

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

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

Должна быть создана mesh, которая будет использоваться для ограничительной оценки. Эти три строки кода после комментария "Задают сетку, что "бесконечные" ограничения будут оценены по" в airPollutionCon может быть изменен для большинства 2-х полубесконечных проблем программирования.

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) плоскость, i.e. в x =-20000m, y = 20000 м. Эта точка отмечена синей точкой на рисунке ниже и проверена путем вычисления концентрации двуокиси серы в этой точке.

Исследуйте файл 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 м/с.

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

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

Реализовывать это ограничение в подходящей форме для fseminf, мы вспоминаем концентрацию SO2 на оптимальной высоте трубы в предыдущей проблеме. В частности, концентрация SO2 принимает свое максимальное возможное значение в x =-20000m, y = 20000 м. Чтобы сократить количество "бесконечных" переменных, мы примем, что концентрация 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. Задайте начальный размер mesh для ограничительной оценки

Код после комментария "Начальный интервал выборки" инициализирует размер mesh.

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 и Алгоритм.