CC BY
52
УДК 66.01
Т. В. Лаптева, Н. Н. Зиятдинов, Д. Д. Первухин,
Г. М. Островский
НИЖНЯЯ ОЦЕНКА ОДНОЭТАПНОЙ ЗАДАЧИ ОПТИМАЛЬНОГО ПРОЕКТИРОВАНИЯ С ВЕРОЯТНОСТНЫМИ ОГРАНИЧЕНИЯМИ
Ключевые слова: оптимизация, оптимальное проектирование, одноэтапная задача, вероятностные
ограничения.
Учет неопределенности в задачах оптимального проектирования часто приводит к появлению операторов математического ожидания и вероятности, требующих многократного расчета многомерных интегралов на каждой итерации оптимизационной процедуры. Существующие процедуры численного интегрирования для достижения высокой точности вычислений требуют большого числа узловых точек, из-за чего растут временные затраты на получение решения. В статье рассматривается подход, сводящий одноэтапную задачу оптимального проектирования с вероятностными ограничениями и математическим ожиданием в критерии к последовательности задач нелинейного программирования с линейными ограничениями. Эффективность подхода демонстрируется на основе решения задачи оптимального проектирования системы реактор-теплообменник.
Keywords: optimization, optimal design, one-stage problem, chance constraints.
Taking into account uncertainty issue in optimal design problems frequently leads to expectation and probability operators which demand multiple calculation of multidimensional integrals on each iteration of optimization procedure. Up to date procedures for numerical integration of high precision require large amount of nodal points - a reason of long time to acquire solution. This paper reviews approach that reduces chance constrained one-stage problem for optimal design which has expectation in criteria to a sequence of NLP-problems with linear constraints. Efficiency of approach is demonstrated on optimal design problem for "reactor - heat exchanger" system.
Задача проектирования оптимальной химико-технологической системы (ХТС) решается в условиях некоторой неточности исходной технологической, физико-химической, экономической информации, что приводит к появлению неточностей в математических моделях. Неточность в исходную информацию также вносят изменения во внешних условиях функционирования и изменения внутренних характеристик ХТС. Решение такой задачи может проводиться в некоторой средней точке по неопределенным параметрам [1], а полученные параметры затем изменяются согласно опыту проектировщика. Полученная ХТС должна быть работоспособной в процессе функционирования несмотря на неточности в исходной информации. Условие работоспособности предусматривает выполнение предъявляемых к ХТС требований, которые в задачу оптимального проектирования включаются в виде ограничений. На ограничения может налагаться требование их безусловного выполнения, либо выполнение с заданной вероятностью.
В зависимости от выделения этапов проектирования и функционирования в жизни ХТС задача проектирования оптимальной ХТС ставится в одноэтапной (ОЭЗО) либо двухэтапной (ДЭЗО) постановках [2]. Будем рассматривать одноэтапную постановку задачи оптимального проектирования ХТС при наличии вероятностных или мягких ограничений.
В качестве критерия будем использовать математическое ожидание E[f(x, 0)] функции, характеризующей эффективность работы ХТС. Задача имеет вид [2]:
f * = min E [f (x,0)] (1)
x
Pr{g (x,0) < 0} >a, (2)
Здесь x - вектор, включающий конструктивные и управляющие параметры, которые в одноэтапной постановке задачи равноправны, 0 - вектор неопределенных параметров, f(x,0)
- критерий оптимальности, д^ (х,в) - ограничения на ХТС. Относительно неопределенных
параметров известно, что их значения принадлежат некоторой области Т
Пусть Т - область изменения параметров в и О] = {в: д(х,в) < 0, ве Т} . Тогда
Е[1(х,в)] = | 1 (х,в)р(в)бв , Рг{д](х,в) < 0} = ^р(в)с1в, где р(в) - плотность вероятности
Т 01
параметров в.
Сложность прямого решения задачи (1) заключается в необходимости многократного расчета (на каждой итерации процедуры оптимизации) многомерных интегралов. Существующие процедуры численного интегрирования с помощью квадратур (кубатур) [3] и метода Монте-Карло (и его производных - метода латинского гиперкуба, метода проб Хаммерслея) [4] требуют большого числа узловых точек, количество которых растет экспоненциально с ростом размерности пространства неопределенных параметров. В результате для достижения удовлетворительной точности могут потребоваться колоссальные временные затраты. В связи с этим весьма важным становится развитие методов решения одноэтапных задач оптимизации.
Мы будем предполагать независимость неопределённых параметров, а также, что их распределение подчинено закону нормального распределения. Пусть область неопределённости является многомерным прямоугольником вида
Т = {в: в^ <в! <ви, / = 1,..., п}. Аппроксимация критерия. Сначала рассмотрим проблему вычисления математического ожидания функции 1 (х,в) относительно вектора случайных параметров в. Для этого мы должны вычислить многомерный интеграл
Е[Г (х,в),Т ] = Е [1 (х,в)] = 11 (х,в)р(в)бв (3)
Т
Предполагая дифференцируемость 1 (х,в) по параметрам в, разложим в ряд Тэйлора функцию 1 (х,в) (по параметрам в) в окрестности некоторой точки в1 области Т, учитывая
члены только первого порядка малости. Обозначим через 1 (х,в,в!) линейную часть этого разложения
Г (х,вв) = Г (хв) + £ ^в1 (в-в!). /=1 в
Величина Еар[1(х,в),Т] = |¿(х,в,в')р(в)бв будет являться некоторым приближением к
Т
(3) на области Т. Для улучшения качества аппроксимации мы будем дробить исходную область Т на подобласти Т,. Причем ^ Т! = Т. На каждой итерации будем дробить ту
I
подобласть Т3, в которой качество аппроксимации функции 1 (х,в) соответствующей функцией 1 (х,в,вэ) наихудшее.
Тогда на некоторой итерации Еар[1 (х,в),Т] = ^Еар[1 (х,в),Т], где к - номер итерации.
I=1
Мы будем использовать два варианта аппроксимации критерия в подобласти:
1. кусочно-линейными в подобластях Т функциями
пв Ы ( х в1 )
ЕаРV(х,в);Т,] = а (х,в!) + (Е[в,Т] - а в);
/=1 ов1
2. кусочно-постоянными в подобластях Т; функциями
ЕаР[1 (х,в);Т,] = а,1 (х,в),
где а, = \р(в)с1в, Е[в,;Г,] = $в!Р(в)с1в.
Т Т,
Пользуясь независимостью параметров в,- и их нормальным распределением Крамер [5], имеем
а = [ф(ви ) - Фв )]-[ф(Ю) - ф(Ю:')];
9й
E [ 9 T ] = [Ф( ~UJ ) - Ф( ~LJ )] - j в; Р( в ; )d 9 . . . [Ф( 1 ) - Ф( в^ 1 )] ,
где = -Ев ¡], к = ^ -Ев ¡].
СТ ст
Аппроксимация ограничений. Перейдем к решению задачи (1). Будем предполагать, что функции д1 (х, в) являются дифференцируемыми по переменным х, в и выпуклыми по в .
ПустьЕ[в ¡], V[в,] - математическое ожидание и дисперсия (ст, )2, (ст,- = V[в /] ), случайной величины .
Подход к решению задачи (1) основывается на сведении ограничений (2) к детерминированному виду. Пусть вначале д1 (х, в) является линейной функцией по
параметрам в, и переменным у : gj(х,в) = в1х1 +... + впхп - Ь, в> 0, / = 1,..., п .
Здесь величина Ь также является случайной величиной, имеющей математическое
ожидание Е[Ь] и дисперсию V[Ь] = (стЬ)2, (стЬ = V[Ь] ). Интервал изменения ТЬ величины Ь
имеет вид ТЬ = {Ь: Ь- < Ь < Ьи}, где Ь- = Е[Ь] - кЬстЬ, Ьи = Е[Ь] + кЬстЬ.
Вероятностное ограничение (2) эквивалентно следующему детерминированному ограничению [6]
Г n
2 E [ в ,\x , +Ф\а).
V'=1
< E[b], (4)
d¿
2v [ в ¡]x 2 + V [b]
¡=1 ,
где Ф(а) - функция нормированного, нормального распределения, Ф-1 (а) - функция, обратная к Ф(а), (Ф(Ф_1(а)) = а)
1 ^ Г ¿
Ф(л) = ^= jexp -z-
<2жi L 2
В общем случае, когда функции g¡ (x, в) являются нелинейными, выпуклыми по в ,
разложим их в ряд Тэйлора в некоторой точке в1, ограничившись членами первого порядка малости. Мы будем использовать линейную часть разложения. Обозначим её через
- , , Л dg¡ (x, в1) ,
g¡ (x, в, в1) = gj (x, в1) + 2 дв (в; - в!). (5)
¡=1 дв;
Можно показать, что в данном случае по аналогии с линейным случаем можно свести ограничения к детерминированному виду (4) [2]. Тогда задача (1) примет вид
f = min E [f ( x,9)] (6)
xeX
n dg ■
2(в!■J -Е[в;])-в-gJ(x,9!■J)
'=1 ' ->Ф-1(а} ),j = 1,..., m. (7)
2v в ] ^
tí ¡ дв!
L,l
При замене функции gj(х,в) в ограничениях (1) их приближениями (5) задача (6) фактически принимает вид
f = min E[f (х,в)];
x
Pr{gj (х,в,в') < 0} >aj, j = 1,..., m. Введем область Qj (в') = {в: gj(х,в,в') < 0;öe T} . Для выпуклых по в функций выполняется условие [7]
gj(х,в) > gj(х,в') + ]Г в
. " dg,(х,в') , ')+Е (вI-в')
, " дд, (х,в) , при любых х . Отсюда, если д] (х,в) < 0, то и д](х,в ) + ^—--(0! - в,) < 0.
I=1 дв,
Следовательно, Пу сП, (в), у = 1,..., т . Отсюда Рг{в е (в)} > Рг{в е Пу}. Тогда, сравнивая задачи (1) и (6)
1 < 1',
т.е. решение задачи (6) даёт нижнюю оценку решения задачи (1) в случае, если функции ду(х,в) являются выпуклыми по в.
Уточнение оценки. Для уточнения получаемой оценки служит итерационная процедура, она основывается на следующих операциях.
Для улучшения аппроксимации ограничений (1) для каждого ограничения будут последовательно добавляться точки вк'', у = 1,..., т, в которых будет проводиться
формирование функций дУ(х,в,вк''), у = 1,..., т. На основе этих функций будут
конструироваться ограничения (7) и добавляться к ограничениям, построенным на предыдущих итерациях.
Для улучшения аппроксимации целевой функции (3) будет проводиться разбиение области Т на подобласти Тр.
Пример. Рассмотрим технологическую систему, состоящую из реактора и теплообменника (рис. 1) [8] с рециклом. В реакторе объема V протекает экзотермическая первого порядка реакция вида А ———> В. Рецикл (с расходом Г1) используется для управления температурой Т1 в реакторе. Противоточный теплообменник служит для охлаждения рециклического потока ^, используя холодную воду с расходом ^. Процесс служит для получения целевого продукта В .
Рис. 1 - Процесс, состоящий из реактора и теплообменника: 1 - реактор, 2 -теплообменник, 3 - компрессор
Математическая модель реактора
Го(ОАо - ОА1 У СА0 = Ук о ехр(-Е/ЯТ! )Ом, (-АНГ (ОЛо - ОЛ1VОЛо = Г>рср(Т - То) + 0НЕ.
Математическая модель теплообменника
ОНЕ = ^р (Т1 - Т2 ) = Р*ср(Tw2 - ^ ),
О = А и [(Т1 - 2 ) - (Т2 - 7*1 )]
НЕ~ 1 \п{(7 - т*2 )/(Т2 - 7,1 )У где Го, То, ОАо - расход (м3/ч), температура (К) входного потока и концентрация вещества А (кмоль/м3) во входном потоке; ко - константа скорости реакции, (м3/(кмоль*ч)); V, Т1, ОА1 -объём реактора (м3), температура в реакторе (К) и концентрация вещества А (кмоль/м3) в выходном потоке; (-АН) - теплота реакции (кДж/кмоль); ^,Т2 - расход (м3/ч) и температура рецикла (К); Ср, Срт - теплоёмкости рециклического потока и охлаждающей воды (кДж/(кг*К)); р, р* - плотность смеси в реакторе и потока холодной воды (кг/ м3); Тш, Т*2, Г* - входная и выходная температура (К), и расход охлаждающей воды (м3/ч); А( - поверхность теплообмена теплообменника (м2); и - полный коэффициент теплопередачи ( кДж/(м2 -ч^)).
Здесь ко = ко / ОАо, где ко - действительная константа скорости реакции.
Таблица 1 - Параметры модели
РСр 167,4 кДж(м3*К) и 1635 кДж(м2*ч*К)
4,19о кДж(м3*К) То 333 К
Го 45,36 м3/ч Т\н 1 3оо К
ОА о 32,о4 кмоль/м3 Е / Я 5бо К
ко 9,81 м /(кмоль*ч) АН -232бо кДж/кмоль
Мы предполагаем, что на этапе проектирования мы не можем получить точные значения параметров 9 = {Го,То,Тт,ко,и}, однако нам известны диапазоны изменения их значений. Область неопределенности будем задавать в виде Т = 0 :9?(1 - уб,) <9, <9^(1 + уд У),/ = 1,...,5}, где 9^ - номинальные значения неопределенных параметров, б - величина, задающая отклонение 0 от 9^, у - параметр, определяющий размер область неопределенности. Номинальные значения неопределенных параметров представлены в таблице 1. Отклонения составляют соответственно 1о%, 2%, 3%, 1о%, 1о%. Будем предполагать, что неопределенные параметры являются независимыми с нормально распределенными значениями.
Ограничения задачи
(Оао - Оа1 )/ Оао < о,9; (8)
7 -Т2 > о; (9)
Т\м 2 - Тт > о ; (1о)
311 < Т2 < 389; (11)
т2 - тт > 11,1; (12)
311 < Т < 389; (13)
7 - Т* 2 > 11,1; (14)
3о1 < Т*2 < 355. (15)
Целевая функция имеет вид
= 691,2 -V0,7 + 873 • At0'6 +1,76 • Fw + 7,056 • F1.
Целевая функция включает капитальные затраты на реактор и теплообменник (первые два члена) и эксплуатационные затраты на охлаждающую воду и перекачку рециклового потока (последние два члена).
Конструктивными параметрами системы являются объем реактора V и поверхность теплообмена At теплообменника. Управляющими переменными выбраны температура в реакторе
T1, температура охлаждающей воды ^2. Будем решать задачу проектирования оптимальной ХТС в постановке одноэтапной задачи оптимизации с вероятностными ограничениями (1), (2).
Ограничения (8)-(12) зависят от неопределенных параметров, ограничения (13)-(15) от неопределенных параметров не зависят. Будем предполагать, что каждое из ограничений (8)-(12) должно выполняться с заданной вероятностью. Среди этих ограничений следует выделить ограничение (9), т.к. его нарушение приводит к тому, что значение расхода рециклового потока F1 становится отрицательным, что противоречит физическому смыслу. Поэтому будем требовать, чтобы это ограничение выполнялось с вероятностью не ниже 0,95. Мы рассмотрели случаи, когда вероятность выполнения ограничений (8), (10)-(12) составила не менее 0,5; 0,75; 0,95.
Таблица 2 - Результаты решения задачи оптимизации
У а Номер итерации
1 2 3 4 5 6
гк 9783,6 9791,2 9801,4 9814,8 9821,6 9821,5
0,5 t 0,078 0,203 0,36 0,641 0,906 1,281
А 0,07% 0,1% 0,13% 0,07% 0,001%
9850,8 9858,6 9868,7 9882,1 9889,0 9888,9
1 0,75 t 0,109 0,218 0,406 0,609 0,875 1,171
А 0,08% 0,1% 0,14% 0,07% 0,001%
гк 9947,8 9955,8 9965,9 9979,4 9986,3 9986,1
0,95 t 0,094 0,203 0,36 0,563 0,813 1,188
А 0,08% 0,1% 0,14% 0,07% 0,002%
fk 9783,6 9795,3 9812,1 9834,8 9839,3 9839,1
0,5 t 0,093 0,265 0,453 0,656 0,921 1,203
А 0,12% 0,17% 0,22% 0,05% 0,002%
fk 9867,6 9879,7 9885,8 9900,7 9916,5 9916,7
1,25 0,75 t 1,109 0,296 0,468 0,687 0,953 1,328
А 0,12% 0,16% 0,14% 0,16% 0,02%
fk 9989,1 10001,6 10018,1 10023,7 10038,2 10060,6
0,95 t 0,093 0,234 0,39 0,593 0,921 1,234
А 0,12% 0,17% 0,05% 0,15% 0,22%
гк 9783,6 9800,1 9818,1 9849,3 9856,4 9856,1
0,5 t 0,093 0,265 0,437 0,781 1,046 1,406
А 0,17% 0,18% 0,32% 0,07% 0,005%
гк 9884,5 9901,6 9909,4 9930,7 9937,7 9937,4
1,5 0,75 t 0,094 0,235 0,453 0,672 0,938 1,266
А 0,16% 0,08% 0,21% 0,07% 0,003%
гк 10030,5 10048,3 10056,0 10077,4 10084,4 10084,1
0,95 t 0,093 0,234 0,421 0,703 0,968 1,312
А 0,18% 0,08% 0,21% 0,07% 0,003%
Для решения задачи используем предложенный подход, сводящий задачу (1), (2) к задаче (6), (7). Для формулирования задачи каждое из вероятностных ограничений (8)-(12) было приведено к виду (7). Результаты решения задачи приведены в таблице 2.
В таблице 2 -к - значение критерия задачи (6), (7), t - время расчета (в секундах), а -вероятность выполнения ограничений (8), (1о)-(12), при этом вероятность выполнения
ограничения (9) равна о,95, А = ——.
-к
Как видно из приведенных результатов к 5-ой итерации было получено решение. Было проведено вычисление методом Монте-Карло критерия оптимальности (1) в оптимальной точке, найденной решением задачи (6), (7). Результаты расчета приведены в таблице 3.
Таблица 3 - Результаты расчета критерия в точке оптимальности
У а fMK fLE А
1 0,5 9794,1 9821,5 0,28%
1 0,75 9861,4 9888,9 0,28%
1 0,95 9958,4 9986,1 0,28%
Где -МК - значение критерия оптимальности, рассчитанное методом Монте-Карло, -
- - -
Л1 или I I I
'MK LE fMK
Время расчета методом Монте-Карло составило не менее 40 минут. Таким образом, приведенные результаты показывают эффективность предложенного подхода в виде значительной экономии вычислительных затрат на решение задач оптимального проектирования.
Литература
1. Зиятдинов, Н.Н. Поиск энергосберегающих режимов работы установки разделения изоамилен-изопреновой фракции производства изопрена / Н.Н. Зиятдинов, Д.А. Рыжов, Т.В. Лаптева, В.А. Курбатов // Вестник Казан. технол. ун-та. - 2009. - № 6. - С. 249-258.
2. Островский, Г.М. Одноэтапная задача с мягкими ограничениями / Г.М. Островский, Н.Н. Зиятдинов, Т.В. Лаптева, Д.Д. Первухин // Теоретические основы химической технологии. -2008. - т. 43, № 4. - С 441-451.
3. Бахвалов, Н.С. Численные методы / Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. - М.: Бином. Лаборатория знаний, 2007. - 696 с.
4. Diwekar, U.M. Efficient Sampling Technique for Optimization under Uncertainty / U.M. Diwekar, JR. Kalagnanam // AIChE J. - 1997. - Vol.43. - P. 440.
5. Крамер Г. Математические методы статистики. - М.: РХД, 2003. - 648 с.
6. Liu, B. Theory and Practice of Uncertain Programming, 3rd Edition - UTLAB, 2009.
7. Базара, М. Нелинейное программирование: теория и алгоритмы / М. Базара, К. Шетти, пер.англ. Т.Д. Березневой и В.А. Березнева, под ред. Д.Б. Юдина. - М.: Мир, 1982. - 583 с.
8. Halemane, K.P. Optimal Process Design under Uncertainty / K.P. Halemane, I.E. Grossmann // AIChE J. -1983. - Vol.29. - Р. 425-433.
© Т. В. Лаптева - канд. техн. наук, проф. каф. системотехники КГТУ, tanlapteva@yandex.ru; Н. Н. Зиятдинов - д-р техн. наук, проф., зав. каф. системотехники КГТУ, nnziat@yandex.ru; Д. Д. Первухин - асп. той же кафедры; Г. М. Островский - д-р техн. наук, проф. каф. системотехники КГТУ, ostralex@yandex.ru.
