Оптимизационные модели в проблеме загрязнения атмосферы

Рисунок 2.1 - Переход к локальным системам координат

Координаты фиксированной и локальной систем связаны соотношением:

(2.3)

Тогда в базовой системе координат имеем:

,(2.4)

где

- координаты источника в базовой системе координат .

2.2 Выбор методов решения поставленной оптимизационной задачи

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

2.2.1 Описание адаптивного метода случайного поиска

Метод случайных направлений - один из методов случайного поиска. Для подобных методов общим является применение случайных чисел.

Задается начальная точка . Каждая последующая точка находится по формуле

,

где - величина шага; - случайный вектор единичной длины, определяющий направление поиска; - номер итерации. На текущей итерации при помощи генерирования случайных векторов получаются точки, лежащие на гиперсфере радиуса с центром в точке . Если значение функции в полученной точке не меньше, чем в центре, шаг считается неудачным (точки при поиске из , при поиске из ).

Если число неудачных шагов из текущей точки достигает некоторого числа М, дальнейший поиск продолжается из той же точки, но с меньшим шагом до тех пор, пока он не станет меньше заранее заданной величины . Если же значение функции в полученной точке меньше, чем в центре, шаг считается удачным и в найденном направлении делается увеличенный шаг, играющий роль ускоряющего шага. Если при этом значение функции снова меньше, чем в центре, направление считается удачным и дальнейший поиск продолжается из этой точки (точки при поиске из , при поиске из ) /11/. Если же значение функции стало не меньше, чем в центре, направление считается неудачным и поиск продолжается из старого центра (в точке при поиске из функция меньше, чем в , а в точке уже не меньше, поэтому направление неудачное).

Алгоритм адаптивного метода случайного поиска

1. Задать начальную точку , коэффициенты расширения и сжатия , М - максимальное число неудачно выполненных испытаний на текущей итерации, - начальную величину шага, - минимальную величину шага, - максимальное число итераций. Положить .

2. Получить случайный вектор

,

где - случайная величина, равномерно распределенная на интервале .

3. Вычислить

.

4. Проверить выполнение условий.

4а. Если ,

шаг удачный. Положить

.

Определить, является ли текущее направление удачным:

если ,

направление поиска удачное. Положить

, ,

и проверить условие окончания. Если , положить и перейти к п. 2. Если , поиск завершить:

;

если ,

направление поиска неудачное, перейти к п. 5

4б. Если ,

шаг неудачный и перейти к п. 5.

5. Оценить число неудачных шагов из текущей точки.

5а. Если , следует положить и перейти к п. 2;

5б. Если , проверить условие окончания:

если , процесс закончить:

, ;

если , положить , и перейти к п. 2 /11/.

2.2.2 Описание метода конфигураций

Метод конфигураций (метод Хука-Дживса) представляет собой комбинацию исследующего поиска с циклическим изменением переменных и ускоряющего поиска по образцу. Исследующий поиск ориентирован на выявление локального поведения целевой функции и определение направления ее убывания вдоль «оврагов». Полученная информация используется при поиске по образцу при движении вдоль «оврагов».

Исследующий поиск начинается в некоторой начальной точке , называемой старым базисом. В качестве множества направлений поиска выбирается множество координатных направлений. Задается величина шага, которая может быть различной для разных координатных направлений и переменной в процессе поиска. Фиксируется первое координатное направление и делается шаг в сторону увеличения соответствующей переменной. Если значение функции в пробной точке меньше значения функции в исходной точке, шаг считается удачным. В противном случае необходимо вернуться в исходную точку и сделать шаг в противоположном направлении с последующей проверкой поведения функции. После перебора всех координат исследующий поиск завершается. Полученная точка называется новым базисом. Если исследующий поиск с данной величиной шага неудачен, то она уменьшается и процедура продолжается. Поиск заканчивается, когда текущая величина шага станет меньше некоторой величины /11/.

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

Обозначим через координатные направления:

.

При поиске по направлению меняется только переменная , а остальные переменные остаются зафиксированными.

Алгоритм метода конфигураций

1. Задать начальную точку , число для остановки алгоритма, начальные величины шагов по координатным направлениям

,

ускоряющий множитель , коэффициент уменьшения шага . Положить

.

2. Осуществить исследующий поиск по выбранному координатному направлению.

2а. Если

, то есть ,

шаг считается удачным. В этом случае положить

и перейти к п. 3.

2б. Если в пункте 2а шаг неудачен, то делается шаг в противоположном направлении. Если

, то есть ,

шаг считается удачным. В этом случае следует положить

и перейти к п. 3.

2в. Если в пунктах 2а и 2б шаги неудачны, положить

.

3. Проверить условия.

3а. Если , то положить и перейти к п. 2 (продолжить исследующий поиск по оставшимся направлениям).

3б. Если , проверить успешность исследующего поиска:

если , перейти в п. 4;

если , перейти в п. 5.

4. Провести поиск по образцу. Положить

,

и перейти к п. 2.

5. Проверить условие окончания.

5а. Если все , то поиск закончить:

.

5б. Для тех , для которых , уменьшить величину шага

.

Положить

и перейти к п. 2.

В алгоритме можно использовать одинаковую величину шага по координатным направлениям, то есть вместо применять /11/.

2.3 Определение максимальной концентрации загрязненности атмосферы промышленного региона

В настоящее время затруднительно найти реальные показатели по выбросам вредных веществ, поэтому изучив литературу и осуществив поиск в сети Internet были подобранны исходные данные близкие к реальным.

Пусть в промышленном регионе имеется три предприятия: . Географическое расположение этих промышленных объектов представлено на рисунке 2.2. Координаты точек на графике:



Поиск максимума проводится в заданной области :

(2.5)

Все необходимые исходные данные каждого предприятия для решения поставленной оптимизационной задачи представлены в таблице 2.1.

Таблица 2.1 - Исходные данные для предприятий

Предприятия	Координаты предприятий, м	Выбросы загрязняющих веществ, г/с	Ско-рость ветра, м/с	Углы, оС	Высота подъемов струи, м
				Азота оксид	Ацетон	Хлор
	4000	7500	1	0,023389	0,674001	0,0031	9	45	180
	1500	4000	1	0,056042	0,000299	0,000001	9	45	250
	5000	1000	1	0,000016	0,000094	0,000006	9	45	200

2.4 Результаты расчетов и их анализ

С помощью описанных методов был произведен расчет максимальной концентрации загрязняющих веществ в атмосфере, в ходе которого необходимо было выявить точку максимального загрязнения. Результаты расчетов приведены в таблице 2.2. Также был проведен расчет при скорости ветра 12 и 20 м/с; угол 60оС и при расположении предприятий: А (1000,3000), В (1000,1000), С (2000,2000). Результаты расчетов приведены в таблицах 2.3 и 2.4.

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

Из таблиц 2.3 и 2.4 можно сделать вывод, что при изменении скорости ветра изменяется значение максимальной концентрации, а ее точка остается прежней.

Результаты полученные с помощью разработанного программного пакета совпадают с расчетами проведенными в системе MathCad, что показано на рисунках 2.3 - 2.8.

Таблица 2.2 - Расчет концентрации оксида азота-(1), ацетона-(2), хлора-(3)

	Адаптивный метод случайного поиска	Метод конфигураций
	Максимальная концентрация,	Искомая точка, м	Кол-во итераций	Максимальная концентрация,	Искомая точка, м	Кол-во итераций
1		894,79	3395,39	200		894,79	3395,39	158
2		641,22	4141,22	244		641,22	4141,22	159
3		641,22	4641,22	233		641,22	4641,22	154

Таблица 2.3 - Расчет концентрации оксида азота-(1), ацетона-(2), хлора-(3)

	Адаптивный метод случайного поиска	Метод конфигураций
	Максимальная концентрация,	Искомая точка, м	Кол-во итераций	Максимальная концентрация,	Искомая точка, м	Кол-во итераций
1		893,11	814,82	230		893,11	814,82	147
2		641,22	2378,58	235		641,22	2378,58	171
3		641,22	2378,58	220		641,22	2378,58	170

Таблица 2.4 - Расчет концентрации оксида азота-(1), ацетона-(2), хлора-(3)

	Адаптивный метод случайного поиска	Метод конфигураций
	Максимальная концентрация,	Искомая точка, м	Кол-во итераций	Максимальная концентрация,	Искомая точка, м	Кол-во итераций
1		893,11	814,82	236		893,11	814,82	148
2		641,22	2378,58	229		641,22	2378,58	172
3		641,22	2378,58	206		641,22	2378,58	171



Рисунок 2.3 - Графическая иллюстрация распределения концентрации оксида азота: а) Линии постоянной концентрации б) Поверхности уровня концентрации



Рисунок 2.4 - Графическая иллюстрация распределения концентрации ацетона а) Линии постоянной концентрации б) Поверхности уровня концентрации

Рисунок 2.5 - Графическая иллюстрация распределения концентрации хлора: а) Линии постоянной концентрации б) Поверхности уровня концентрации

Рисунок 2.6 - Графическая иллюстрация распределения концентрации оксида азота при расположении предприятий: А (1000, 3000), В (1000, 1000), С (2000, 2000) а) Линии постоянной концентрацииб) Поверхности уровня концентрации



Рисунок 2.7 - Графическая иллюстрация распределения концентрации ацетона при расположении предприятий: А (1000, 3000), В (1000, 1000), С (2000, 2000): а) Линии постоянной концентрацииб) Поверхности уровня концентрации



Рисунок 2.8 - Графическая иллюстрация распределения концентрации хлора при расположении предприятий:А (1000, 3000), В (1000, 1000), С (2000, 2000): а) Линии постоянной концентрацииб) Поверхности уровня концентрации

2.5 Программная реализация расчета максимальной концентрации загрязненности промышленного региона

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

Рисунок 2.6 -Главная форма программы

В правой части формы пользователю предоставляется возможность задать начальную точку (стартовую точку) и высоту наблюдения. Далее пользователь выбирает метод, которым он будет решать поставленную оптимизационную задачу (рисунок 2.7).



Рисунок 2.7 - Форма выбора метода решения задачи с заданием начальной точки и высоты наблюдения

После выбора метода решения оптимизационной задачи становится активной кнопка «Данные». При ее нажатии появляется окно, в котором пользователь выбирает файл исходных данных «input» (рисунок 2.8 - 2.9).

Рисунок 2.8 - Окно загрузки исходных данных «input»

Файл с исходными данными «input» представлен на рисунке 2.9.



Рисунок 2.9 - Файл с исходными данными «input»

После выбора файла с исходными данными становится активной кнопка «Расчет» (рисунок 2.10), при нажатии которой выводится окно с результатами расчета значения максимальной концентрации загрязненности атмосферы и координаты точки максимальной загрязненности (рисунок 2.11), также результаты выводятся в файл «output», в котором еще и записывается количество итераций.

Рисунок 2.10 - Форма с активной кнопкой «Расчет»



Рисунок 2.11 - Окно с результатами расчетов программы

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

3. МНОГОКРИТЕРИАЛЬНАЯ ОПТИМИЗАЦИЯ ПРЕДЕЛЬНО ДОПУСТИМЫХ ВЫБРОСОВ ЗАГРЯЗНЕННОСТИ АТМОСФЕРЫ ПРОМЫШЛЕННОГО РЕГИОНА ОТ НЕСКОЛЬКИХ ИСТОЧНИКОВ

3.1 Постановка задачи

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

(3.1)

Решением поставленной оптимизационной многомерной задачи будет совокупность переменных

.

Если в (3.1) обозначить

, , ,

то формула примет вид:

,(3.2)

где ,(3.3)

.

Получаем задачу оптимизации:

(3.4)

Для частного случая :

(3.5)

Тогда общая задача векторной оптимизации примет вид:

(3.6)

Получили задачу многокритериальной многомерной оптимизации.

3.2 Выбор методов решения поставленной многокритериальной оптимизационной задачи

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

3.2.1 Метода главного критерия для решения задачи многокритериальной оптимизации

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

Пусть

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

по важности:

(3.7)

и выберем в качестве главного наиболее важный критерий. Означим его через

.

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

(3.8)

К допустимому множеству альтернатив добавим условие (3.8). Расширенное, таким образом, допустимое множество будем обозначать . Таким образом, вместо многокритериальной задачи решается однокритериальная по известным методам решения.

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

/12/.

3.2.2 Метод последовательных уступок для решения задачи многокритериальной оптимизации

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

(3.9)

По наиболее важному критерию определяется наибольшее значение и назначается величина допустимого снижения (уступки). Рассматриваются в дальнейшем только те альтернативы (решения), значения частичных критериев которых входят в назначенные допуски.

Метод «уступок» позволяет использование на практике вместо оптимальных так называемые субоптимальные (близкие к оптимальным) решения /12/. Данные решения незначительно отличаются от оптимальных по главному критерию, но данная «уступка» - отклонение от оптимального значения - даст возможность учесть дополнительные ограничения по другим, вспомогательным критериям.

Максимизация критериев записывается в виде

,(3.10)

где

допустимое множество решений (альтернатив.)

Пусть

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

(3.11)

Затем находится максимальное значение второго критерия

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

,(3.12)

где ;

- максимальное значение по -му критерию.

Аналогично для критериев назначается величина допустимой уступки

расчет концентрация загрязнение атмосфера выброс

и находится максимальное значения вплоть и для критерия :

,(3.13)

где ;

- максимальное значение по -му критерию;

- величина допустимой «уступки» для -го критерия.

При этом получим сужающуюся последовательность допустимых множеств решений (альтернатив):

,

где .

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

3.3 Расчеты и анализ полученных результатов

С помощью описанных методов многокритериальной оптимизации были найдены значения мощностей выбросов загрязняющих веществ каждым предприятием, при которых суммарная концентрация вредного вещества каждого типа не превышает предельно допустимых концентраций (ПДК). Исходные данные для расчета были взяты из таблицы 2.1 главы 2 дипломной работы.

Значения ПДК для рассматриваемых загрязняющих /18/ веществ представлены в таблице 3.1.

Таблица 3.1 - Предельно допустимые концентрации загрязняющих веществ в атмосфере

Наименование загрязняющего вещества	ПДК,
Азота оксид	0,0004
Ацетон	0,00001
Хлор	0,0001

Задача расчета выбросов загрязняющих веществ в атмосферу сводилась к задаче многокритериальной оптимизации и решалась методами главного критерия и последовательных уступок. В частности, в методе главного критерия, поочередно, каждый критерий назначался главным, при этом на остальные накладывались ограничения и проводился расчет. При пользовании же методом последовательных уступок для проведения расчета критерии упорядочивались в порядке убывания важности, после чего выбирались уступки и проводился расчет. Результаты проделанной работы приведены в таблицах 3.2 и 3.3.

Таблица 3.2 - Расчет максимальной мощности выбросов концентрации оксида азота, хлора и ацетона методом главного критерия

№ ГК	Ограничения на остальные функции	Максимально возможное значение функции выброса по главному критерию,	Координаты точки выброса
Расчет максимальной мощности выбросов концентрации Оксида азота
1				6898,79	6440,86
2				4902,99	4968,42
3				2173,47	4949,06
Расчет максимальной мощности выбросов концентрации Ацетона
1				5861,66	5566,07
2				4963,66	5004,16
3				2378,23	4871,64
Расчет максимальной мощности выбросов концентрации Хлора
1				5701,13	5749,57
2				5564,88	6130,47
3				2026,21	4532,03

Таблица 3.3 - Расчет максимальной мощности выбросов концентрации оксида азота, хлора и ацетона методом последовательных уступок

№ ВК	Уступки на остальные функции	Максимально возможные значеня мощности выброса каждым предприятием,	Координаты точки выброса
Расчет максимальной мощности выбросов концентрации Оксида азота
1,2,3						6901,63	6456,76
2,1,3						4962,87	4974,11
3,1,2						2201,58	5003,94
Расчет максимальной мощности выбросов концентрации Ацетона
1,2,3						5893,63	5498,12
2,1,3						4982,91	4994,95
3,1,2						2392,56	4863,97
Расчет максимальной мощности выбросов концентрации Хлора
1,2,3						5693,74	5740,78
2,1,3						5549,87	6139,35
3,1,2						2003,54	4525,87

Таблица 3.4 - Сравнение полученных результатов двумя методами: методом главного критерия и методом последовательных уступок

	Максимально возможные значения функции выбросов,	Координаты точки выброса
	МГК	МПУ	МГК	МПУ
Расчет максимальной мощности выбросов концентрации Оксида азота
1			6898,79	6440,86	6901,63	6456,76
2			4902,99	4968,42	4962,87	4974,11
3			2173,47	4949,06	2201,58	5003,94
Расчет максимальной мощности выбросов концентрации Ацетона
1			5861,66	5566,07	5893,63	5498,12
2			4963,66	5004,16	4982,91	4994,95
3			2378,23	4871,64	2392,56	4863,97
Расчет максимальной мощности выбросов концентрации Хлора
1			5701,13	5749,57	5693,74	5740,78
2			5564,88	6130,47	5549,87	6139,35
3			2026,21	4532,03	2003,54	4525,87

Из таблицы 3.4 видно, что оба метода многокритериальной оптимизации дают почти одинаковые результаты, связано это с тем, что ограничения в методе главного критерия и уступки в методе последовательных уступок выбирались согласовано. А небольшая разница в полученных результатах объясняется недостаточно точным согласованием ограничений и уступок и конечно использованием разных методов.

Значения максимальных мощностей выбросов при данных значениях ПДК, приведенные в таблицах 3.2 и 3.3 принадлежат пространству оптимальных решений поставленной задачи. Выбор в качестве решения одного из вариантов зависит от того, насколько максимальная концентрация при этом решении приближается к ПДК. Так для оксида азота это

,

для ацетона: ,

для хлора: .

3.4 Программная реализация расчета значения мощностей выбросов загрязняющих веществ промышленного региона

Созданный программный пакет решает задачу многокритериальной многомерной оптимизации, не накладывая ограничения на количество источников выбросов (рисунок 3.1).

Рисунок 3.1 - Главная форма программы

В правой части формы пользователю предоставляется возможность задать начальную точку (стартовую точку), высоту наблюдения и загрузить файл с исходными данными, для этого необходимо нажать на кнопку «Данные». При ее нажатии появляется окно, в котором выбирается файл исходных данных «input» (рисунок 3.2).

Рисунок 3.2 - Окно загрузки исходных данных «input»

алее пользователь выбирает метод, которым он будет решать поставленную многокритериальную многомерную оптимизационную задачу. Если пользователь выбрал метод главного критерия, то ему необходимо задать номер главного критерия и ограничения на другие функции (рисунок 3.3). Если пользователь выбрал метод последовательных уступок, то ему необходимо задать величину уступки (рисунок 3.4).



Рисунок 3.3 - Форма с выбором метода главного критерия

Рисунок 3.4 - Форма с выбором метода последовательных уступок

После выбора файла с исходными данными и вводом необходимых условий для каждого многокритериального метода пользователь выбирает метод, которым он будет решать полученную однокритериальную задачу (рисунок 3.5).



Рисунок 3.5 - Форма с активной кнопкой «Расчет»

Далее становится активной кнопка «Расчет» (рисунок 3.5), при нажатии которой результаты расчетов программы выводятся в файл «output».

ЗАКЛЮЧЕНИЕ

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

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

В разделе безопасности жизнедеятельности приведен сравнительный анализ практических методик расчета концентраций вредных примесей в атмосфере.

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

СПИСОК ЛИТЕРАТУРЫ

Атмосферная турбулентность и моделирование распространения примесей, // Под ред. Ф.Т.М. Ньистадт и Х. Ван Доп; Л: Гидрометеоиздат, 1985.

Безопасность жизнедеятельности в условиях чрезвычайных ситуаций / Под ред. Н.Г. Евтушенко, А.П. Кузьмина. -- Курган, 1994.

Берлянд М.Е. и др. Об усовершенствовании методов расчета загрязнения атмосферы. // Труды ГГО, 1987, вып. 511. с. 3-23.

Берлянд М.Е., Прогноз и регулирование загрязнения в атмосфере, Л: Гидрометеоиздат, 1985.

Берлянд М.Е., Современные проблемы атмосферной диффузии и загрязнения атмосферы, Л: Гидрометеоиздат, 1975.

Бузало Н.С. Математическое моделирование переноса примеси в мезометеорологическом пограничном слое атмосферы // Дисс. работа на соискание уч. ст. канд. наук, 2003.

Гаврилов А.С. Программный комплекс ZONE - интеллектуальная геоинформационная система для управления качеством атмосферы города. // Информационный бюллетень ГИС ассоциации N1(13) 1998.

Гигиенические нормативы ГН 2.1.6.1338-03 "Предельно допустимые концентрации (ПДК) загрязняющих веществ в атмосферном воздухе населенных мест" (утв. Главным Государственным санитарным врачом РФ 31 мая 2003 г.)

Гигиенические нормативы ГН 2.1.6.1338-03 "Предельно допустимые концентрации (ПДК) загрязняющих веществ в атмосферном воздухе населенных мест" (утв. Главным Государственным санитарным врачом РФ 31 мая 2003 г.)

Замай С.С. Якубайлик О.Э. Модели оценки и прогноза загрязнения атмосферы промышленными выбросами в информационно-аналитической системе природоохранных служб крупного города. Учебное пособие. Красноярск, 1998.

Корнеенко В.П., Методы оптимизации. - М.: Высшая школа, 2007 г. - 664 с.

Марчук Г.И., Математическое моделирование в проблеме окружающей среды, М.: Наука, 1971. - 350 с.

Методика расчета концентраций в атмосферном воздухе вредных веществ, содержащихся в выбросах предприятий (ОНД-86). Л.: Гидрометеоиздат, 1987. - 93 с.

Методика расчета концентраций в атмосферном воздухе вредных веществ, содержащихся в выбросах предприятий (ОНД-86). Л.: Гидрометеоиздат, 1987. - 93 с.

Миргородский В. Способы, средства и особенности ликвидации химически опасных аварий.// Мир и безопасность. - №6. - 2000.

Никифоров А.Н., Методы оптимизации: Учеб. Пособие / Юж.-Рос. Гос. Техн. Ун-т. - Новочеркасск: ЮРГТУ, 2006.

СДЯВ и защита от них. / Под ред. В.А. Владимирова. - М.: Воениздат, 1989.

Тищенко Н.Ф., Охрана атмосферного воздуха. Расчет содержания вредных веществ и их распределения. - М.: Наука, 1991.

ПРИЛОЖЕНИЕ

using System;

using System.Collections.Generic;

using System.Linq;

using System.Windows.Forms;

namespace Optimization

{

static class Program

{

/// <summary>

/// The main entry point for the application.

/// </summary>

[STAThread]

static void Main()

{

Application.EnableVisualStyles();

Application.SetCompatibleTextRenderingDefault(false);

Application.Run(new Form1());

}

}

}

using System;

using System.Collections.Generic;

using System.ComponentModel;

using System.Data;

using System.Drawing;

using System.Text;

using System.Windows.Forms;

using System.IO;

using System.Text.RegularExpressions;

namespace Optimization

{

class random_search

{

const int n = 250; //максимальное число итераций

const int m = 450; //максимальное чило испытаний

static double alpha1 = 1.618;

static double beta = 0.618;

static double[] t = new double[n + 1];//вектор шага на каждой итерации

static double R = 1e-15; //величина минимального шага шага

static int dim = 2;

public static double[] ARS_PU(int dimIV, double[] VV, double[][] xyz, double[] alpha, double[] HH,

double[] QQ, double zz, double shag0, double[] x0, int NCF, double[] c_k, double[] UST)

{

double fx_extr;

double Xmax = 0; double Ymax = 0;

for (int i = 0; i < dimIV; i++)

{

if (Xmax < xyz[i][0])

{

Xmax = xyz[i][0];

}

if (Ymax < xyz[i][1])

{

Ymax = xyz[i][1];

}

}

t[0] = shag0; ;// шаг на нулевой итерации

double[] x = new double[dim + 2];

double[] xrez = new double[dim + 2];

x[0] = x0[0]; x[1] = x0[1];

double[] y = new double[dim]; double[] z = new double[dim];

double[] ksi = new double[dim]; double[] tvec = new double[dim];

bool priznak = false;//priznak = true если шаг стал <= R

double fx = Function(dimIV, VV, xyz, alpha, HH, QQ, x, zz, NCF);

double fy; double fz; int num_it = 0; int num_isp = 0;

int pr1 = 0;

for (int k = 0; k < n; k++)//счетчик итераций

{

for (int j = 1; j <= m; j++)//счетчик испытаний

{

ksi = GetRandomVector(dim);// случайный вектор

tvec = MultNV(t[k], VdivNorm(ksi, dim), dim);

y = Add(x, tvec, dim); // точка в случайном направлении

if (Cond_PU(dimIV, VV, xyz, alpha, HH, QQ, y, zz, NCF, c_k, UST))

{

pr1 = 1;

xrez[0] = y[0]; xrez[1] = y[1];

}

else

{

pr1 = 0;

}

fx = Function(dimIV, VV, xyz, alpha, HH, QQ, x, zz, NCF);

fy = Function(dimIV, VV, xyz, alpha, HH, QQ, y, zz, NCF);

if (fy > fx) // проверка на удачность шага //<

{

tvec = MultNV(alpha1, Sub(y, x, dim), dim);//alpa*(y-x)

z = Add(x, tvec, dim);// точка, расположенная дальше в том же направлении

if (Cond_PU(dimIV, VV, xyz, alpha, HH, QQ, z, zz, NCF, c_k, UST))

{

pr1 = 1;

xrez[0] = z[0]; xrez[1] = z[1];

}

else

{

pr1 = 0;

}

fz = Function(dimIV, VV, xyz, alpha, HH, QQ, z, zz, NCF);

if (fz > fx)// проверка удачности направления//<

{

for (int i = 0; i < dim; i++)//направление удачное x = z

{

if (pr1 == 1)

{

x[i] = xrez[i];

}

else

{

x[i] = z[i];

}

}

t[k + 1] = alpha1 * t[k];//шаг для след итерации делаем больше

break;//прекращаем испытания

}

}

else

{

t[k + 1] = t[k];

}

if (j == m)//все испытания неудачные

if (t[k] > R)//можно ли уменьшать шаг

{

t[k] = t[k] * beta;// уменьшили шаг

k--;//надо пройти итерацию снова уменьшив шаг

}

else

{

priznak = true;

num_it = k; num_isp = j;

}

}

if (priznak) break;// шаг нельзя уменьшать, прекращаем итер

}

fx_extr = fx;

if (Math.Abs(fx_extr) <= 0)

{

MessageBox.Show("Решение не полученно!");

}

xrez[dim] = num_it;

xrez[dim + 1] = num_isp;

return xrez;

}

public static double[] HukJivs_PU(int dimIV, double[] VV, double[][] xyz, double[] alpha, double[] HH,

double[] QQ, double zz, double shag0, double[] x0, int NCF, double c_k, double ust_k)

{

double f1 = 0; double f0 = 0; double lam = 0.5;

double[] y1 = new double[dim + 1];

double[] yt = new double[dim];

double[] x1 = new double[dim];

double DD = 100;

bool priznak = false;

for (int ii = 0; ii < 2; ii++)

{

y1[ii] = x0[ii];

yt[ii] = x0[ii];

}

for (int k = 0; k < n; k++)//счетчик итераций

{

f0 = Function(dimIV, VV, xyz, alpha, HH, QQ, y1, zz, NCF); //1

for (int i = 0; i < 2; i++)//счетчик по базису

{

yt[i] = yt[i] + DD;//2

if (Function(dimIV, VV, xyz, alpha, HH, QQ, yt, zz, NCF) >= c_k - ust_k)

{

f1 = Function(dimIV, VV, xyz, alpha, HH, QQ, yt, zz, NCF);//2

if (f1 > f0)//2a

{

for (int ii = 0; ii < 2; ii++)

{

y1[ii] = yt[ii];//2a

}

priznak = true;

yt[i] = yt[i] - DD;

}

else

{

yt[i] = yt[i] - 2 * DD;//2b

if (Function(dimIV, VV, xyz, alpha, HH, QQ, yt, zz, NCF) >= c_k - ust_k)

{

f1 = Function(dimIV, VV, xyz, alpha, HH, QQ, yt, zz, NCF);//2b

if (f1 > f0)//2b

{

for (int ii = 0; ii < 2; ii++)

{

y1[ii] = yt[ii];//2b

}

priznak = true;

yt[i] = yt[i] + DD;

}

}

}

}

if (priznak == false)//2в

{

for (int ii = 0; ii < 2; ii++)

{ y1[ii] = yt[ii]; }

}

for (int ii = 0; ii < 2; ii++)//???

{ yt[ii] = y1[ii]; }

}

if (Function(dimIV, VV, xyz, alpha, HH, QQ, y1, zz, NCF) >= c_k - ust_k)

{

f1 = Function(dimIV, VV, xyz, alpha, HH, QQ, y1, zz, NCF);//3b

if (f1 > f0)//3b

{

for (int ii = 0; ii < 2; ii++)//4

{

x1[ii] = y1[ii];

y1[ii] = x1[ii] + lam * (x1[ii] - yt[ii]);

}

}

else//3b

{

if (DD <= R)//5

{

break;//5a

}

else

{

DD = DD / alpha1;//5b

}

}

}

/**/

y1[dim] = k;

}

return y1;

}

public static double[] ARS_GK(int dimIV, double[] VV, double[][] xyz, double[] alpha, double[] HH,

double[] QQ, double zz, double shag0, double[] x0, int NCF, double[] OGR)

{

double fx_extr;

double Xmax = 0; double Ymax = 0;

for (int i = 0; i < dimIV; i++)

{

if (Xmax < xyz[i][0])

{

Xmax = xyz[i][0];

}

if (Ymax < xyz[i][1])

{

Ymax = xyz[i][1];

}

}

t[0] = shag0; ;// шаг на нулевой итерации

double[] x = new double[dim + 2];

double[] xrez = new double[dim + 2];

x[0] = x0[0]; x[1] = x0[1];

double[] y = new double[dim]; double[] z = new double[dim];

double[] ksi = new double[dim]; double[] tvec = new double[dim];

bool priznak = false;//priznak = true если шаг стал <= R

double fx = Function(dimIV, VV, xyz, alpha, HH, QQ, x, zz, NCF);

double fy; double fz; int num_it = 0; int num_isp = 0;

int pr1 = 0;

for (int k = 0; k < n; k++)//счетчик итераций

{

for (int j = 1; j <= m; j++)//счетчик испытаний

{

ksi = GetRandomVector(dim);// случайный вектор

tvec = MultNV(t[k], VdivNorm(ksi, dim), dim);

y = Add(x, tvec, dim); // точка в случайном направлении

fx = Function(dimIV, VV, xyz, alpha, HH, QQ, x, zz, NCF);

fy = Function(dimIV, VV, xyz, alpha, HH, QQ, y, zz, NCF);

if (fy > fx) // проверка на удачность шага //<

{

if ((Cond(dimIV, VV, xyz, alpha, HH, QQ, y, zz, NCF, OGR)) && (fy > fx))

{

pr1 = 1;

xrez[0] = y[0]; xrez[1] = y[1];

}

else

{

pr1 = 0;

}

tvec = MultNV(alpha1, Sub(y, x, dim), dim);//alpa*(y-x)

z = Add(x, tvec, dim);// точка, расположенная дальше в том же направлении

fz = Function(dimIV, VV, xyz, alpha, HH, QQ, z, zz, NCF);

if (fz > fx)// проверка удачности направления//<

{

if ((Cond(dimIV, VV, xyz, alpha, HH, QQ, z, zz, NCF, OGR)) && (fz > fx))

{

pr1 = 1;

xrez[0] = z[0]; xrez[1] = z[1];

}

else

{

pr1 = 0;

}

for (int i = 0; i < dim; i++)//направление удачное x = z

{

if (pr1 == 1)

{

x[i] = xrez[i];

}

else

{

x[i] = z[i];

}

}

t[k + 1] = alpha1 * t[k];//шаг для след итерации делаем больше

break;//прекращаем испытания

}

}

else

{

t[k + 1] = t[k];

}

if (j == m)//все испытания неудачные

if (t[k] > R)//можно ли уменьшать шаг

{

t[k] = t[k] * beta;// уменьшили шаг

k--;//надо пройти итерацию снова уменьшив шаг

}

else

{

priznak = true;

num_it = k; num_isp = j;

}

}

if (priznak) break;// шаг нельзя уменьшать, прекращаем итер

}

fx_extr = fx;

if (Math.Abs(fx_extr) <= 0)

{

MessageBox.Show("Решение не полученно!");

}

xrez[dim] = num_it;

xrez[dim + 1] = num_isp;

return xrez;

}

public static double[] HukJivs_GK(int dimIV, double[] VV, double[][] xyz, double[] alpha, double[] HH,

double[] QQ, double zz, double shag0, double[] x0, int NCF, double[] OGR)

{

double f1 = 0; double f0 = 0; double lam = 0.5;

double[] y1 = new double[dim + 1];

double[] xrez = new double[dim + 1];

double[] yt = new double[dim];

double[] x1 = new double[dim];

int pr1 = 0;

double DD = 100;

bool priznak = false;

for (int ii = 0; ii < 2; ii++)

{ y1[ii] = x0[ii]; yt[ii] = x0[ii]; }

for (int k = 0; k < n; k++)//счетчик итераций

{

f0 = Function(dimIV, VV, xyz, alpha, HH, QQ, y1, zz, NCF); //1

for (int i = 0; i < 2; i++)//счетчик по базису

{

yt[i] = yt[i] + DD;//2

f1 = Function(dimIV, VV, xyz, alpha, HH, QQ, yt, zz, NCF);//2

if (f1 > f0)//2a

{

if ((Cond(dimIV, VV, xyz, alpha, HH, QQ, yt, zz, NCF, OGR)) && (f1 > f0))

{

pr1 = 1;

xrez[0] = yt[0]; xrez[1] = yt[1];

}

else

{

pr1 = 0;

}

for (int ii = 0; ii < 2; ii++)

{

y1[ii] = yt[ii];//2a

}

priznak = true;

yt[i] = yt[i] - DD;

}

else

{

yt[i] = yt[i] - 2 * DD;//2b

f1 = Function(dimIV, VV, xyz, alpha, HH, QQ, yt, zz, NCF);//2b

if (f1 > f0)//2b

{

if ((Cond(dimIV, VV, xyz, alpha, HH, QQ, yt, zz, NCF, OGR)) && (f1 > f0))

{

pr1 = 1;

xrez[0] = yt[0]; xrez[1] = yt[1];

}

else

{

pr1 = 0;

}

for (int ii = 0; ii < 2; ii++)

{

y1[ii] = yt[ii];//2b

}

priznak = true;

yt[i] = yt[i] + DD;

}

}

if (priznak == false)//2в

{

for (int ii = 0; ii < 2; ii++)

{

y1[ii] = yt[ii];

}

}

for (int ii = 0; ii < 2; ii++)//???

{

yt[ii] = y1[ii];

}

}

f1 = Function(dimIV, VV, xyz, alpha, HH, QQ, y1, zz, NCF);//3b

if (f1 > f0)//3b

{

if ((Cond(dimIV, VV, xyz, alpha, HH, QQ, y1, zz, NCF, OGR)) && (f1 > f0))

{

pr1 = 1;

xrez[0] = y1[0]; xrez[1] = y1[1];

}

else

{

pr1 = 0;

}

for (int ii = 0; ii < 2; ii++)//4

{

x1[ii] = y1[ii];

y1[ii] = x1[ii] + lam * (x1[ii] - yt[ii]);

}

}

else//3b

{

if (DD <= R)//5

{

break;//5a

}

else

{

DD = DD / alpha1;//5b

}

}

/**/

xrez[dim] = k;

}

return xrez;

}

public static double[] xyz_conv(double[] x, double z, double[] xyz, double al)

{

double[] xyz2 = new double[3];

xyz2[0] = (x[0] - xyz[0]) * Math.Cos(al) + (x[1] - xyz[1]) * Math.Sin(al);

xyz2[1] = -(x[0] - xyz[0]) * Math.Sin(al) + (x[1] - xyz[1]) * Math.Cos(al);

xyz2[2] = z - xyz[2];

return xyz2;

}

public static double Function(int dimIV, double[] VV, double[][] xyz, double[] alpha,

double[] HH, double[] QQ, double[] xx, double zz, int NCF)

{

double s = 0;

double[,] DD = new double[dimIV, 2];//дисперсии

double[] xConv;

int i = NCF;

double dx = Math.Abs(xx[0]);

DD[i, 0] = (0.22 * dx) / Math.Sqrt(1 + 0.0001 * dx);

DD[i, 1] = (0.2 * dx);

xConv = xyz_conv(xx, zz, xyz[i], alpha[i]);

s = (QQ[i] / (2 * Math.PI * VV[i] * DD[i, 0] * DD[i, 1])) *

Math.Exp(-xConv[1] * xConv[1] / (2 * DD[i, 0] * DD[i, 0])) *

(Math.Exp(-(xConv[2] - HH[i]) * (xConv[2] - HH[i]) / (2 * DD[i, 1] * DD[i, 1])) +

Math.Exp(-(xConv[2] + HH[i]) * (xConv[2] + HH[i]) / (2 * DD[i, 1] * DD[i, 1])));

return s;

}

public static double[] FuncArr(int dimIV, double[] VV, double[][] xyz, double[] alpha,

double[] HH, double[] QQ, double[][] xx, double zz)

{

double[] s = new double[dimIV];

double[,] DD = new double[dimIV, 2];//дисперсии

double[] xConv;

for (int i = 0; i < dimIV; i++)

{

double dx = Math.Abs(xx[i][0]);

DD[i, 0] = (0.22 * dx) / Math.Sqrt(1 + 0.0001 * dx);

DD[i, 1] = (0.2 * dx);

xConv = xyz_conv(xx[i], zz, xyz[i], alpha[i]);

s[i] = (QQ[i] / (2 * Math.PI * VV[i] * DD[i, 0] * DD[i, 1])) *

Math.Exp(-xConv[1] * xConv[1] / (2 * DD[i, 0] * DD[i, 0])) *

(Math.Exp(-(xConv[2] - HH[i]) * (xConv[2] - HH[i]) / (2 * DD[i, 1] * DD[i, 1])) +

Math.Exp(-(xConv[2] + HH[i]) * (xConv[2] + HH[i]) / (2 * DD[i, 1] * DD[i, 1])));

}

return s;

}

public static double[] Add(double[] x, double[] y, int dim)// сложение векторов

{

double[] z = new double[dim];

for (int i = 0; i < dim; i++)

z[i] = x[i] + y[i];

return z;

}

public static double[] Sub(double[] x, double[] y, int dim)// вычитание векторов

{

double[] z = new double[dim];

for (int i = 0; i < dim; i++)

z[i] = x[i] - y[i];

return z;

}

public static double[] MultNV(double a, double[] x, int dim)// умножение вектора на число

{

double[] z = new double[dim];

for (int i = 0; i < dim; i++)

z[i] = a * x[i];

return z;

}

public static double[] VdivNorm(double[] x, int dim)// вектор деленный на свою норму

{

double[] z = new double[dim];

double norm = 0.0;

for (int i = 0; i < dim; i++)

{

norm = norm + x[i] * x[i];

}

norm = Math.Sqrt(norm);

z = MultNV(1.0 / norm, x, dim);

return z;

}

public static double[] GetRandomVector(int dim)// вектор случайных величин от -1 до +1

{

Random random = new Random();

double[] ksi = new double[dim];

for (int i = 0; i < dim; i++)

ksi[i] = 2.0 * random.NextDouble() - 1.0;

return ksi;

}

public static bool Cond(int dimIV, double[] VV, double[][] xyz, double[] alpha, double[] HH,

double[] QQ, double[] xt, double zz, int NCF, double[] OGR)

{

bool PR = true;

for (int i = 0; i < dimIV; i++)

{

if (i != NCF)

{

if (Function(dimIV, VV, xyz, alpha, HH, QQ, xt, zz, i) < OGR[i])

{

PR = false;

break;

}

}

}

return PR;

}

public static bool Cond_PU(int dimIV, double[] VV, double[][] xyz, double[] alpha, double[] HH,

double[] QQ, double[] xt, double zz, int NCF, double[] c_k, double[] UST)

{

bool PR = true;

for (int i = 0; i < NCF; i++)

{

if (Function(dimIV, VV, xyz, alpha, HH, QQ, xt, zz, i) < c_k[i] - UST[i])

{

PR = false;

break;

}

}

return PR;

}

}

}

Размещено на Allbest.ru

...