авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 |   ...   | 2 | 3 || 5 | 6 |   ...   | 8 |

«Д.С. ДВОРЕЦКИЙ, С.И. ДВОРЕЦКИЙ, Г.М. ОСТРОВСКИЙ НОВЫЕ ПОДХОДЫ К ПРОЕКТИРОВАНИЮ ХИМИКО-ТЕХНОЛОГИЧЕСКИХ ПРОЦЕССОВ, АППАРАТОВ И СИСТЕМ В УСЛОВИЯХ ...»

-- [ Страница 4 ] --

Для устранения трудоемкой операции вычисления многомерного интеграла при решении задачи ДЗИП1 будем использовать формули ровку двухэтапной задачи оптимизации, в которой минимизируется верхняя граница целевой функции. Введем семейство областей, обладающих свойством Pr{ } =.

Тогда эта задача может быть записана следующим образом:

C= min u, (2.84) a, d, z ( ),u, max C (a, d, z (), ) u, (2.85) 1 (a, d ) 0. (2.86) Обозначим решение этой задачи через a, d, z (), u,. По скольку ограничение C (a, d, z (), ) u 0 удовлетворяется в ка ждой точке области, то вероятность удовлетворения этого нера венства равна. Это означает, что с вероятностью целевая функ ция будет меньше, чем u.

Рассмотрим задачу получения верхней оценки величины u. Для этого выберем какую-либо одну область 1 из семейства и решим задачу (2.84) – (2.86) при фиксированной области = min u, a, d, z ( ),u max C (a, d, z (), ) u, 1 (a, d ) 0.

Пусть мы получили оптимальное значение u1. Поскольку взята произвольная область из семейства, то имеет место соотношение u1 u. Таким образом, с вероятностью целевая функция будет меньше, чем u1. Поскольку при = 1 семейство состоит из одной области, то в этом случае мы получаем точное решение. Отсюда ясно, что чем ближе к единице, тем ближе u1 к u.

Если все параметры i являются независимыми и имеют нор мальное распределение, то область является многомерным прямо угольником вида = { : iN k i i i iN + ki i, I = 1,..., n }, где N номинальная точка и величина k i определяются следую 1 / n щим образом: Pr { iN ki i i iN + k i i } =.

Если параметры i не являются независимыми, но имеют нор мальное распределение, то область является многомерным эллип соидом.

Структура задачи (2.84) – (2.86) близка к структуре ДЗИП1 и не требует вычисления многомерных материалов, но трудоемкость ее решения будет намного меньше трудоемкости решения ДЗИП1.

Перейдем к рассмотрению метода решения ДЗИП2. Обобщим развитый в работе [47] метод разбиений и границ (алгоритм 2.5) для решения задачи ДЗИП2.

Рассмотрим алгоритм вычисления верхней границы оптимального значения целевой функции ДЗИП2 – решения задачи полубесконечно го программирования вида:

i r C (a, d, zi, 1,i, 2,r ), C2,( ) = min U (2.87) a, d, z i, z l iI rQ max g j (a, d, z i, 1,i, 2 ) 0, j = 1,..., m, i I1, (2.88) 2 max g j (a, d, z l, 1, 2 ) 0, l = 1,..., N, j = 1,..., m. (2.89) 11, ( ) l 2 { } Пусть l( ) = 1, 2 : 1 1 ;

2 l2. Введем множество крити l ческих точек S 2, p ) для каждой подобласти l :

( l { } S 2, p ) = 1,l, s, 2,l, s : (1,l, s, 2,l, s ) l, s I 2, p ), ( ( l l где I 2, p ) – множество номеров критических точек, принадлежащих по ( l добласти l( ), p номер итерации в алгоритме внешней аппроксима ции решения задачи (2.87) – (2.89), номер итерации в алгоритме ре шения ДЗИП2. При решении задачи (2.87) – (2.89) номер итерации фиксирован, а число подобластей и их размеры остаются постоянными.

Алгоритм вычисления верхней границы оптимального значения целевой функции ДЗИП2 (решения задачи (2.87) – (2.89)) имеет сле дующий вид.

Алгоритм 2. Шаг 1. Положим p = 1, зададим множество подобластей l( ) (l = 1,..., N ) z i,( 0), z l,( 0), a (0), d (0) и начальные значения (i I1, l = 1,..., N ) соответствующих переменных, множества I1, Q1 и S 2,1), достаточно малое число 0.

( l Шаг 2. Решаем задачу i r C (a, d, z i, 1,i, 2,r ), min a, d, z i, z l iI rQ g j ( a, d, z i, 1,i, 2,i,q ) 0, i I1, 2,i, q S2, j = 1,..., m, i g j (a, d, z l, 1,l, s, 2,l, s ) 0, j = 1,..., m, (1,l, s, 2,l, s ) S 2, p ), l 1,..., N.

( l [ ] Пусть a, d, z i, z l – решение этой задачи.

Шаг 3. Для каждой точки L( ) решаем задачу max g j ( a, d, z 1,i, 1,i, 2 ), j = 1,..., m. (2.90) 2 Шаг 4. Для каждой подобласти l (l = 1,..., N ) решаем задачу max g j (a, d, z l, 1, 2 ), j = 1,..., m. (2.91) l 2 ( ) Обозначим через l, j = 1,l, j, 2,l, j решение задачи (2.91) и через решение задачи (2.90).

2,l, j Шаг 5. Проверяем выполнение условий g j (a, d, z 1,i, 1,i, 2,i, j ), j = 1,..., m, i I1, (2.92) g j (a, d, z l, 1,l, j, 2,l, j ), j = 1,..., m, l = 1,..., N. (2.93) Если неравенства (2.92), (2.93) удовлетворяются, то решение най дено, переходим к шагу 9;

в противном случае – к шагу 6.

Шаг 6. Сформируем множества { } R i = 2,i, j : g j (a, d, z i, 1,i, 2,i, j ), = { ) }.

: g j (a, d, z l, 1,l, j, 2,l, j Rl l, j Множества R i и R l содержат точки, в которых ограничения на рушаются.

Шаг 7. Сформируем новые множества критических точек { U R l }, S 2i = {S 2i U R i }.

S 2, p +1) = S 2, p ) ( ( l l Шаг 8. Положим p : = p + 1 и переходим к шагу 2.

Шаг 9. Сформируем множество R ( ), являющееся объединением всех множеств R l, l = 1,..., N :

{ } R ( ) = 1,l, j, l = 1,..., N, j = 1,..., m.

Так как множество R ( ) соответствует последней итерации, то точки этого множества соответствуют активным ограничениям, т.е.

S AP = R ( ). Множество S AP будет использоваться в качестве S 2 ) в () ( ()..

задаче вычисления нижней границы:

i rC (a, d, z1,i, 1,i, 2,r ), С2,( ) = min L (2.94) a, d, z i, z r iI rQ max g j (a, d, z i, 1,i, 2 ) 0, j = i,..., m, i I1, (2.95) 2 max g j (a, d, z 1, s, 1, s, 2 ) 0, 1, s S 2 ), j = 1,..., m.

( (2.96) 2 Приведем далее алгоритм решения ДЗИП2, реализующий метод разбиений и границ. Для реализации метода разбиений и границ необ ходимо иметь правило разбиения области 1 и правило выбора множества критических точек в задаче (2.94) – (2.96) вычисления ниж ней границы. Аналогично методу разбиения и границ для решения ДЗИП1 мы будем использовать следующее эвристическое правило разбиения области 1 : будут разбиваться только те подобласти 1,( ) (l = 1,..., N ), для которых ограничения (2.89) в задаче (2.87) – l (2.89) вычисления верхней границы будут активны в точке решения, т.е. если выполняются следующие условия:

max g j (a ( ), d ( ), z l,( ), ) = 0, j J 1, ( ) l [ ] где a ( ), d ( ), z l,( ) – решение задачи (2.87) – (2.89).

В качестве множества S 2 ) в задаче (2.94) – (2.96) необходимо ( использовать множество S AP активных точек задачи (2.87) – (2.89) ().

вычисления верхней границы.

Введем множество L( ) подобластей 1,( ) следующим образом:

i { } L( ) = 1,( ) : r (1,( ) ), i i где некоторое положительное число, которое будет настраиваться во время работы алгоритма решения ДЗИП2.

Алгоритм 2. Шаг 1. Положим = 1 и зададим начальное разбиение области 1 на подобласти 1,(1), (l = 1,..., N1 ), множество аппроксимационных l точек 1,i, 2,r, i I1, r Q1, начальное множество критических точек S 20), начальные значения z i,( 0), z l,( 0), a ( 0), d ( 0) (i I1, l = 1,..., N1 ) ( 1 0, 2 0, 1 0, соответствующих переменных, величины 2 0, ( 2 1, 1 2 ), где 1 и 2 малые величины. Положим С2,( 0) = c, С 2,( 0) = c, где c достаточно большое число (c C2 ).

U L Шаг 2. Вычислим верхнюю границу С2,( ) решением задачи U с использованием алгоритма Пусть (2.87) – (2.89) 2.6.

[ ] ( ) ( ) i,( ) i,( ) (i I1, l = 1,..., N ) решение этой задачи.

a,d,z,z { } Шаг 3. Определим множество Q ( ) = 1,( ) : l I Q ) подобластей ( l 1,( ), которым соответствуют активные ограничения в задаче i (2.87) – (2.89) U,l (a ( ), d ( ) ) = 0, 1,( ) Q ( ).

2 l Шаг 4. Если Q ( ) пустое множество, то решение ДЗИП2 найде но, алгоритм останавливает свою работу.

Шаг 5. Если выполняется условие С2,( 1) C 2,( ) 1 C 2,( ), U U U то решение ДЗИП2 найдено, алгоритм останавливает свою работу.

В противном случае проверяем выполнение условия С2,( 1) С2,( ) 2 C2,( ), U U U и если оно нарушается, то переходим к шагу 8.

С2,( ), решая задачу L Шаг 6. Определяем нижнюю границу (2.94) – (2.96). Здесь мы полагаем S 2 ) = S AP.

( ().

Шаг 7. Если выполняется условие С2,( ) С2,( ) 1 С2,( ), U L U то решение ДЗИП2 найдено, алгоритм останавливает свою работу.

Шаг 8. Если выполняется условие r (1,( ) ) 2, i = 1,..., N, то i решение ДЗИП 2 найдено, алгоритм заканчивает свою работу.

Шаг 9. Если выполняется условие r ( l( ) ) 1, l = 1,..., N, то переходим к шагу 11.

Шаг 10. Формируем множество L( ) и находим множество V ( ) подобластей 1,( +1), принадлежащих одновременно множествам L( ) i и Q ( ), т.е. V ( ) = L( ) I Q ( ).

Разбиваем каждую подобласть 1,( ) V ( ) на две подобласти i 1,( ) = 1,( +1) U 1,( +1), образуем новое множество ( +1) подобла i i i 1 ( ), заменяя каждую подобласть стей из старого множества 1,( ) (1,( ) V ( ) ) новыми подобластями i( +1) и i( +1), и переходим i i 1 к шагу 12.

Шаг 11. Положим 1 : = 1 / 2 и переходим к шагу 9.

Шаг 12. Положим : = + 1 и переходим к шагу 2.

В работе [31] доказано, что: 1) решение ДЗИП2 найдено, если множество Q ( ) пусто;

2) алгоритм 2.7 сходится. Конечно, это утвер ждение верно в локальном смысле, т.е. метод разбиений и границ дает по крайней мере локальный минимум ДЗИП2.

2.4. АЛГОРИТМЫ РЕШЕНИЯ ОДНО- И ДВУХСТАДИЙНЫХ ЗАДАЧ ИНТЕГРИРОВАННОГО ПРОЕКТИРОВАНИЯ ХТС С «МЯГКИМИ» ОГРАНИЧЕНИЯМИ При формулировании ОЗИП для случая 3 предположим, что мы имеем полную информацию относительно функции распределения вероятностей для.

В качестве функции I (a d, z ) мы будем использовать среднее C ( a, d, z, ), значение первоначальной целевой функции т.е.

I ( a, d, z ) = M {C ( a, d, z, )}, либо верхнюю границу исходного крите рия С (a, d, z, ) с доверительной вероятностью 0. В этом случае одностадийная задача интегрированного проектирования имеет вид min M {C (a, d, z, )}, a,d, z { } P()d j, Вер g j (a, d, z, ) 0 = j = 1,..., m, j { } j = : g j (a, d, z, ) 0,, где P() функция плотности вероятности;

j заданная вероятность выполнения j-го ограничения.

Можно записать другую формулировку ОЗИП, в которой в каче стве критерия будет использоваться его верхняя граница, которая не может быть нарушена с заданной вероятностью 0. Перепишем задачу ОЗИП для случая 3 в виде min u, a, d, z,u Вер { C (a, d, z, ) u 0} 0, { } Вер g j (a, d, z, ) 0 j, j = 1,..., m.

Здесь u новая переменная, ограничивающая целевую функцию C (a, d, z, ), каждое ограничение должно удовлетворяться с вероятно стью не меньшей, чем j.

[a, d ] решение этой задачи. Тогда конструкция, z, u Пусть ХТС ( a, d ) и режим z гарантируют, что в течение всего этапа функционирования ХТС целевая функция C (a, d, z ) будет меньше, чем u с вероятностью 0. Используя ту же целевую функцию, можно свести задачу с жесткими ограничениями (2.50) – (2.52) к следующей задаче с жесткими и мягкими ограничениями:

min u, a, d, z,u Вер { C (a, d, z, ) u 0} 0, max g j (a, d, z, ) 0, j = 1,..., m.

Мы рассмотрим две формулировки ОЗИП для случая 3, которые часто встречаются на практике.

1. В стохастическом программировании обычно предполагают, что функции распределения вероятностей известны. Однако на прак тике это предположение не выполняется, и функции распределения вероятностей зачастую неизвестны. В этом случае принимают допу щение о равномерном (равновероятном) законе распределения вероят ностей неопределенных параметров и формулировка ОЗИП имеет сле дующий вид: требуется определить тип a аппаратурного оформления ХТС, векторы конструктивных параметров d* технологического обо рудования и режимных (управляющих) переменных z такие, что С (a, d, z ) = min M {C (a, d, z, y (a, d, z, ), )} (2.97) a,d, z при связях в форме уравнений математической модели химико технологического объекта y = (a, d, z, ) (2.98) и ограничениях { } Вер g j (a, d, z, ) y j зад y j 0 j, j J. (2.99) Перепишем задачу (2.97) – (2.99) в терминах А-задач стохастиче ского программирования: требуется определить тип a A аппаратур ного оформления ХТС, векторы конструктивных d* и управляющих z* переменных, а также m-мерный вектор постоянных величин * = (1*, 2*, …, m*) такие, что С (a, d, z ) = min min i C (a, d, z, i ) | g j (a, d, z, ) j, j J, a, d, z i = I (2.100) [ ] { } = | j Вер g j (a, d, z, ) 0 j, i = 1 ;

I1 – множество аппрокси где i – весовые коэффициенты, iI мационных точек в области.

Идея такого подхода, в сущности, очень проста. Поясним ее на примере одномерной задачи стохастического программирования с од ним ограничением g (a, d, z, ) 0. На рисунке 2.5, а заштрихована недопустимая область ограничения. Пусть соотношение между целе вой функцией M {C (a, d, z, )} и g (a, d, z, ) такое, как показано на рис. 2.5, б. Следует заметить, что такое соотношение (кроме, конечно, экзотических случаев) в оптимизационных задачах химической техно логии бывает всегда, т.е. наиболее предпочтительные значения целе вой функции лежат в недопустимой области (см. рис. 2.5), поскольку в противном случае ограничение было бы неактивным и его не следова ло бы учитывать. В этом случае решение z' традиционной задачи оп тимизации достигается при g (a, d, z, ) = 0. Очевидно, при реализа ции этого решения z' значения g (a, d, z, ) будут иметь случайный разброс вследствие наличия случайной величины. На рисунке 2.5, в показан этот разброс, который может имитироваться на математиче ской модели (a, d, z, ).

В зависимости от вхождения случайной величины в функцию g (a, d, z, ) закон распределения этой функции может изменяться.

Следовательно, эта вероятность может быть как меньше, так и больше 0,5. Таким образом, при решении традиционной задачи при = мы даже не знаем, какова вероятность нарушения технологических огра ничений.

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

Идея А-задач стохастического программирования заключается в следующем: исходное ограничение задачи заменяется на ограничение вида g (a, d, z, ), где 0, т.е. исходное ограничение как бы уже сточается (рис. 2.5, г). После этого решается детерминированная зада ча оптимизации с новыми ограничениями { } min C (a, d, z, ) | g (a, d, z, ) 1.

z При этом решение задачи z будет соответствовать тому, что технологические ограничения g (a, d, z, ) будут равными 1 (см.

рис. 2.5, г). Соответственно, вероятность нарушения ограничения а) M {C (a, d, z, )} g (a, d, z, ) б) в) P ( g ) dg P(g) ( S1 = = 1 ) (1 зад ) g (a, d, z ', ) = г) P(g) S (1 зад ) ( S 2 = 2 ) ( S1 = 1 ) g ( a, d, z '', ) = 1 0 д) P(g) S 3 = 3 = (1 зад ) g ( a, d, z *, ) = * Рис. 2.5. Геометрическая иллюстрация идеи решения А-задачи стохастической оптимизации уменьшается по сравнению с 1, т.е. 2 1, а значение целевой функции возрастает (см. рис. 2.5, б). Таким образом, мы приблизились к опти мальному решению задачи, которое изображено на рис. 2.5, д. Отметим, что при выполнении этой процедуры мы не вычисляли вероятность вы полнения (нарушения) ограничения на каждом шаге поиска z*. Вычисле ние Вер [ g (a, d, z, ) 0] производится в оптимальной точке u* для то го, чтобы проверить выполнение условия Вер [ g (a, d, z, ) 0] зад.

В том случае, если эти условия не выполняются, выбирается новое число 2 1 0 и вновь решается детерминированная задача оптими зации с ограничением g (a, d, z, ) 2. Процедура продолжается до тех пор, пока не будет найдено такое *, при котором технологическое ограничение g (a, d, z, ) 0 выполняется с заданной вероятностью, т.е. Вер [ g (a, d, z, ) 0] зад или 3 1 – зад.

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

В соответствии с методом А-задач стохастической оптимизации нами разработан следующий алгоритм решения задачи (2.100) [10].

Алгоритм 2. Шаг 1. Полагаем число итераций = 1, задаем значения вектора гарантированной вероятности j, j = 1,..., m и точности решения { } задачи ОЗИП, множество аппроксимационных точек S1 = i : i I1, ( ) и начальные при начальное значение вектора = 10), ( 0),..., ( 0) (0) ( m ( 0) ( 0) ( 0) ближения a, d, z.

Шаг 2. Методом последовательного квадратичного программи рования решаем задачу НЛП i C (a, d, z, y(a, d, z, i )) I (a, d, z ) = min (2.101) a,d, z iI при связях y = ( a, d, z, i ) и ограничениях g j (a, d, z, y (a, d, z, i )) (j ), (j ) 0, j = 1,..., m, i I1. (2.102) Шаг 3. В точке (a ( 0 ), d ( 0 ), z ( 0 ) ), которая является решением задачи (2.101)–(2.102), вычисляются вероятности выполнения ограни чений g j (a, d, z, ) y j зад y j 0 с использованием генератора слу чайных чисел i, i I1 с равномерным законом распределения и ма тематической модели y = (a, d, z, i ), i I и проверяется выполнение условий { } Вер g j (a, d, z ) 0 j, j = 1,..., m.

Шаг 4. Если вероятностные ограничения не выполняются, т.е.

( ), включается алгоритм входа в допустимую область. Про стейшим алгоритмом такого типа является уменьшение j() для нару шенных ограничений. Далее число увеличивается на 1, т.е. : = + и следует переход к шагу 2.

Шаг 5. Если вероятностные ограничения выполняются, то век тор * находим из решения внешней А-задачи оптимизации I (a, d, z ) = min I (a, d, z ). (2.103) В общем случае задача (2.103) может быть решена подходящим методом нелинейного программирования. Однако нами применялcя простейший алгоритм коррекции вектора путем увеличения его компонентов на величину j = ( ) (Вер [ g j (•) 0] j ), где () – шаг коррекции на -й итерации, подбираемый опытным пу тем. Поиск * прекращается, если j для j становится меньше зара нее заданного малого числа (точность поиска *).

Вычисление вероятностных интегралов производится стандарт ными методами (HSS, Монте-Карло).

2. Пусть = (1, 2). Для совокупности случайных параметров 1 функции распределения вероятностей известны, в то время как для совокупности случайных параметров 2 2 функции распределе ния вероятностей неизвестны, где 1 и 2 – области неопределенности параметров 1 и 2 соответственно. В этом случае формулировка ОЗИП имеет вид { } min max M 1 C (a, d, z, 1, 2 ), a,d, z { } min Вер 1 g j (a, d, z, 1, 2 ) 0 j, j = 1,..., m, 2 где { } P ( ) d, Вер 1 g j (a, d, z, 1, 2 ) 0 = 1 j (2 ) { } j ( 2 ) = 1 : g j (a, d, z, 1, 2 ) 0, 1 1.

Можно также сформулировать одноэтапную задачу оптимизации с усредненными ограничениями:

min M { C (a, d, y (), z, )}, a,d, z y = (a, d, z, ), { } M g j (a, d, y (), z, ) 0, j = 1,..., m, где y () многомерная функция неопределенных параметров.

Сравнительный анализ результатов решения одноэтапных задач интегрированного проектирования с жесткими и усредненными огра ничениями показывает, что { } M g j (a, d, y (), z, ) max ( g j (a, d, z, )), и оптимальное значение целевой функции ОЗИП с усредненными огра ничениями лучше (меньше), чем оптимальное значение целевой функ ции ОЗИП (2.50) – (2.52) с жесткими ограничениями (теорема П2 [31]).

Предположим, что мы решили ОЗИП и получили решение [ a, d, z ]. Чтобы реализовать это решение, мы должны поддержи вать выполнение условия z = z* на стадии функционирования ХТС. Это означает, что мы не можем настраивать управляющие переменные на стадии функционирования. Ясно, что постановка и решение ОЗИП на стадии проектирования приводит к неэкономичной конструкции ХТС, так как не используется настройка управляющих переменных. Пред положим теперь, что мы используем ДЗИП1 для решения задачи с же сткими и мягкими ограничениями. Это означает, что задача С (a, d, ) = min C (a, d, z, ), z g j (a, d, z, ) 0, j = 1,..., m, решается в каждой точке области. Отсюда следует, что все ограни чения g j (a, d, z, ) 0, j = 1,..., m, будут удовлетворяться в каждой точке области с вероятностью i = 1. Использование ДЗИП3 даст более экономичную конструкцию ХТС, так как мягкие ограничения могут удовлетворяться с вероятностью i 1.

Далее мы дадим строгую формулировку двухстадийной задачи интегрированного проектирования в случае мягких ограничений. Мы будем предполагать, что функции распределения вероятностей извест ны, все ограничения являются мягкими и должны быть удовлетворены с вероятностью зад. В качестве целевой функции этой задачи мы бу дем использовать среднее (ожидаемое) значение некоторой меры, ха рактеризующей работу ХТС на стадии функционирования. В этом случае формулировка ДЗИП для случая 3 (ДЗИП3) имеет следующий вид:

( ) I1 = min I11) (a, d ) + I1( 2) (a, d ), ( (2.104) d P() d зад, (2.105) C C (a, d, ) P() d ;

I11) (a, d ) = (a, d, )P () d ;

I1 2 ) ( a, d ) = ( ( где 1 \ C (a, d, ) решение min C (a, d, z, ) | g j (a, d, z, ) 0, задачи z j = 1,..., m ;

1 = { : h(a, d, ) 0} ;

\ 1 теоретико-множественная h(a, d, ) разность;

функция определяется формулой h(a, d, ) = min max g j (a, d, z, ).

jJ z [ ] решение этой задачи. В противоположность Пусть a, d 1, z1 () ДЗИП1 задача (min C (a, d, z, ) | g j (a, d, z, ) 0, j = 1,..., m) решается z только в точках области 1. Так как в каждой точке области выполняет ся условие h(a, d, ) 0, то задача (min C (a, d, z, ) | g j (a, d, z, ) 0, z j = 1,..., m) имеет решение в каждой точке области 1, и в каждой точке области 1 можно найти такие значения z (), что все ограни чения будут выполняться. Вероятность попадания в область больше или равна зад. Ясно, что в области 1 выполняется условие max min max g j (a, d, z, ) 0. Если зад 1, то 1 и в пределе 1 jJ z при зад = 1 (если существуют векторы a и d такие, что 1 (a, d ) 0 ) получаем I12 (a, d ) = 0. При этом задача (2.104), (2.105) преобразуется в двухэтапную задачу интегрированного проектирования с жесткими ограничениями (2.10), (2.11) (ДЗИП1).

Для реализации найденного оптимального режима на стадии функционирования можно использовать следующий подход. Согласно задаче (2.104), (2.105) задача min C (a, d, z, ) | g j (a, d, z, ) 0, z j = 1,..., m должна решаться только внутри области { } 1 = : h(a, d 1, ) 0 (см. соотношение 1 = { : h(a, d, ) 0} ).

Будем использовать следующее правило: запомним уравнение поверх ности области 1 и если на этапе функционирования измеренное зна чение будет находиться внутри области 1 (т.е. 1 ), то мы должны решать задачу min C (a, d, z, ) | g j (a, d, z, ) 0, j = 1,..., m, в z ( 1 ) – задачу ( min C (a, d, z, ) ). В этом случае противном случае z мы должны запоминать только одну многомерную функцию, хотя это тоже может потребовать огромного объема памяти.

В связи с этим введем штрафную функцию m p( g j (a, d, z, )) С (a, d, z,, ) = C (a, d, z, ) +, (2.106) j = где 0, если g j 0, p( g j ) = 2 (2.107) g j, если g j 0.

Здесь p( g j (a, d, z, )) штрафной член и штрафной коэффи циент. Таким образом, для всех значений мы будем использовать следующую внутреннюю задачу оптимизации:

С (a, d,, ) = min С (a, d, z,, ). (2.108) z Ясно, что оптимальное значение критерия этой задачи зависит от a, d,,. Пусть z (a, d,, ) решение этой задачи. Переменные a, d, должны быть выбраны таким образом, чтобы ограничения удовлетворялись с вероятностью зад и среднее значение величины С (a, d,, ) принимало минимальное значение. Поэтому формули ровка ДЗИП3 будет иметь вид { } I = min M C (a, d,, ), (2.109) a, d, P()d зад, (2.110) { } где 2 = : g j (a, d, z (a, d,, ), ) 0.

Неравенство (2.110) гарантирует удовлетворение ограничений с вероятностью не меньшей зад. В определении области 2 использу ются значения z ( a, d,, ) управляющих переменных, полученных решением внутренней задачи оптимизации (2.108). В задаче (2.109) штрафной коэффициент используется как дополнительная поиско вая переменная.

Рассмотрим две формулировки ДЗИП с усредненными мягкими ограничениями. Первая формулировка имеет вид I = min M { C (a, d, z (), )}, (2.111) a, d, z () { } M g j (a, d, z (), ) 0, j = 1,..., m.

Заметим, что здесь опять поисковые переменные z () являются многомерными функциями.

Во второй формулировке мы будем использовать задачу (2.108) в качестве внутренней задачи оптимизации. Оптимальное значение кри терия этой задачи зависит от a, d,,. Переменные a, d, должны быть выбраны таким образом, чтобы ограничения удовлетворялись в среднем на этапе функционирования. В этом случае ДЗИП с усреднен ными мягкими ограничениями будет иметь вид { } I1 = min M C (a, d,, ) (2.112) a, d, { } M g j (a, d, z (a, d,, ), ) 0, j = 1,..., m.

Следует заметить, что решение z (a, d,, ) задачи (2.108) ис пользуется в ограничениях задачи (2.112). Используя задачу (2.108), мы можем преобразовать задачу (2.112) к виду I1 = min min C (a, d, z,, ) P () d, (2.113) a,d z { } M g j (a, d, z (a, d, z (a, d, z (a, d, )), ) 0, j = 1,..., m. (2.114) Изменим порядок операций интегрирования и минимизации в формуле (2.113):

C (a, d, z(), ) P() d, I 2 = min (2.115) a,d, z () M { g j (a, d, z (a, d,, ), )} 0, j = 1,..., m.

Задача (2.115) не эквивалентна задаче (2.113), поскольку пере менные z в критерии (2.113), соответствующие разным значениям, не независимы (они связаны условиями (2.114)): I 2 I1. Это означает, что решение задачи (2.115) определяет только верхнюю границу опти мального значения задачи (2.113).

Следующая процедура последовательных приближений позволяет получить решение задачи (2.113).

Алгоритм 2. Шаг 1. Положим = 1. Зададим начальное приближение z ( 0 ) ( ) для управляющих переменных и малую величину 1 0.

Шаг 2. Решаем следующую задачу:

I1 ) = min min C (a, d, z,, ) P () d, ( (2.116) a,d z { } M g j (a, d, z ( 1) (), ) 0, j = 1,..., m. (2.117) Пусть a ( ), d ( ), z ( ) () решение этой задачи, а z ( 1) () есть значение вектора управляющих переменных, полученное на предыду щей итерации.

Поскольку ограничения (2.117) связывают только конструктив ные переменные a ( ), d ( ), то в отличие от задачи (2.113) управляю щие переменные z задачи (2.116), соответствующие разным значениям, независимы. Поэтому мы можем изменить порядок операций интег рирования и минимизации в (2.116):

I 2 ) = min C (a, d, z,, ) P() d, ( a,d, z () { } M g j (a, d, z ( 1) (), ) 0, j = 1,..., m.

Шаг 3. Если I1 ) I1 1) 1 I1 ), то решение задачи (2.113) ( ( ( найдено, и алгоритм останавливает свою работу.

Шаг 4. Положим = + 1 и переходим к шагу 2.

Далее рассмотрим методы решения ДЗИП3, когда каждое j-е ограничение должно выполняться с заданной вероятностью j [152].

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

C (a, d, z(), ) P() d, С1 = min a,d, z ( ) Pr{g j (d, z (), ) 0} j, j = 1,..., m + p, p число g m+l hl (a, d, z ()) 0, l = 1,..., p ;

где ограничений hl (a, d, z ()) функции, образующие область H, которой принадле жат переменные d, z, т.е. H = {d, z : hl (a, d, z ) 0, l = 1,..., p}.

Будем искать приближенное решение сформулированной задачи в предположении, что управляющие переменные z () аппроксимиру ются линейными функциями ~1 () = b 0 + b11 +... + b n, где bi век n z тор размерности nz = dim z, компоненты которого будем обозначать через bk, k = 1,..., nz.

i Поисковыми переменными теперь будут компоненты векторов bi.

Подставив ~1 () в исходную задачу, получим z С1 = min C (a, d, b, ) P() d, a, d,b Pr{g j (d, b, ) 0} j, j = 1,..., m + p, где n C (a, d, b, ) f 0 (a, d, b 0 + b11 +... + b n, ), n g j (a, d, b, ) f j (a, d, b0 + b11 +... + b n, ), j = 1,..., m + p.

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

Рассмотрим теперь случай, когда в качестве целевой функции двухстадийной задачи используется верхняя граница исходной целе вой функции с доверительной вероятностью 0. В этом случае двух стадийная задача с мягкими ограничениями имеет вид min u, a,d, z ( ),u Pr{C (a, d, z (), ) u 0} 0, Pr{g j (a, d, z (), ) 0} j, j = 1,..., m + p.

Опять будем искать приближенное решение задачи в предполо жении, что управляющие переменные z () являются линейными функциями параметров. В этом случае ДЗИП3 имеет следующий вид:

min u, a,d,b,u Pr{C (a, d, b, ) u 0} 0, Pr{g j (a, d, b, ) 0} j, j = 1,..., m + p.

Вид задачи такой же, как у одноэтапной вероятностной задачи с мягкими ограничениями.

2.5. АЛГОРИТМЫ РЕШЕНИЯ ОДНО- И ДВУХСТАДИЙНЫХ ЗАДАЧ ОПТИМИЗАЦИИ СО СМЕШАННЫМИ ОГРАНИЧЕНИЯМИ При формулировании ДЗИП для случая 4 (ДЗИП4) предположим, что: 1) ограничения с номерами j J1 = (1,..., m1 ) являются мягкими и должны удовлетворяться с вероятностью j, а ограничения с номера ми j J 2 = (m1 + 1,..., m) – жесткими, J J1 I J 2 = (1,..., m). В этом случае формулировка ДЗИП4 будет иметь вид [31, 48] min I (a, d ) = I1 (a, d ) + I 2 (a, d ), (2.118) a,d Вер [ 1 ] зад, (2.119) 1 (a, d ;

J 2 ) 0, (2.120) где, область 1 определяется как 1 = : min max g j (a, d, z, ) 0, jJ z (min C (a, d, z, ) | g j (a, d, z, ) 0, j J ) P() d;

I1 ( a, d ) = z [ min C ( a, d, z, ) + A max{ max g j ( a, d, z, ), 0} I 2 (a, d ) = jJ z \ | g j ( a, d, z, ) 0, j J ] P () d;

1 (a, d ;

J 2 ) = max min max g j ( a, d, z, ), jJ z A штрафной коэффициент;

Вер [ 1 ] зад – величина стохасти ческой гибкости [151].

Здесь 1 – множество тех значений, для которых могут быть выполнены ограничения задачи, и Вер [ 1 ] зад. В критерии оп тимизации для каждого 1 переменная z должна выбираться из условия минимума C (a, d, z, ) при выполнении ограничений h(a, d, ) = min max g j (a, d, z, ) 0, а при 1 из условия мини jJ z мизации функции, учитывающей величину C (a, d, z, ) и штраф за нарушение ограничений g j (a, d, z, ) 0.

Ограничение (2.119) гарантирует удовлетворение всех ограниче ний с вероятностью зад. Кроме того, дополнительное ограничение (2.120) гарантирует жесткое удовлетворение ограничений с номерами j J 2 для всех точек области. Если существует такие a, d, что 1 (a, d ;

J 2 ) 0, зад 1, 1, то при и в пределе при зад = 1, I 2 (a, d ) = 0 задача (2.118) – (2.120) превращается в ДЗИП1.

Пусть в качестве целевой функции двухстадийной задачи исполь зуется верхняя граница исходной целевой функции с доверительной вероятностью 0. В этом случае формулировка ДЗИП4 может быть представлена и в другом виде (см. п. 2.1):

I = min u, (2.121) a, d,u, z ( ) Вер { g 0 = С (a, d, z (), ) u 0} 0, (2.122) { } Вер g j (a, d, z (), ) 0 j, j J1, (2.123) 1 (a, d ) = max min max g j (a, d, z, ) 0. (2.124) jJ z В задаче (2.121) – (2.124) u скалярная переменная (аналог кон структивных переменных);

Вер {•} вероятность выполнения ограни {•} ;

g 0, g j функции g 0 (a, d, z (), ) = чения ограничений;

= С (a, d, z (), ) – критерий интегрированного проектирования;

g j (a, d, z, ) y j зад y j 0, j = 1,..., m функции ограничений;

y j = (a, d, z, ), j = 1,..., m;

0, j заданные значения вероятности выполнения ограничений;

1 (d ) функция гибкости ХТС.

Введем обозначения g j (a, d, z, ) u, j = 0;

g j ( a, d, u, z, ) = g j (a, d, z, ), j J1 ;

{ } и множество S ( ) = i : i I ( ) накопления точек с индексами iI ( ), в которых нарушаются ограничения (2.122) – (2.124), причем во множестве точек S1( ) будут накапливаться точки, в которых нару шаются жесткие ограничения, а во множестве S 2 ) – точки, в которых ( нарушаются мягкие ограничения. Кроме того, мы будем использовать вспомогательную задачу нелинейного программирования (В):

I = min u;

d,u, z i g j (d, u, z i, i ) 0, j J1, i I ( ) ;

(В) ( ) g j (d, z, ) 0, j J 2, i I i i.

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

Алгоритм 2.10 [49] Шаг 1. Принимаем = 1, задаем начальные множества S ( 0), I ( 0) из условия наилучшей аппроксимации функций z (), число номеров n то чек i, i I ( 0) и начальные приближения a ( 0), d ( 0), u ( 0), z i,( 0), i I (0).

Шаг 2. Решаем вспомогательную задачу (В) I = min u;

d,u, z i g j (d, u, z i, i ) 0, j J1, i I ( 1) ;

g j (d, z i, i ) 0, j J 2, i I ( 1).

И пусть a ( ), d ( ), u ( ), z ( ) есть решение этой задачи.

Шаг 3. Вычисляем ( ) 1 (d ( ) ) = max min max g j a ( ), d ( ), z, (2.125) jJ z с использованием алгоритма внешней аппроксимации. Обозначим че рез ( ) решение задачи (2.125) и проверяем выполнение условия 1 (d ( ) ) 0. (2.126) в точке решения ( ) задачи (2.125). Если условие (2.126) не выполня ется, то переходим к шагу 4, в противном случае – к шагу 5.

Шаг 4. Дополним множество точек S1( ), в которых нарушаются ограничения (2.126), точкой ( ), т.е.

S1( ) = S1( 1) U ( ), I1( ) = I1 1) U (n + 1);

n := n + 1.

( Шаг 5. Проверяем выполнение мягких (вероятностных) ограни чений { } Вер g j (a ( ), d ( ), z (), ) 0 j, j J1. (2.127) На данном шаге мы не имеем функций z = z (), а известны толь ко значения этих функций в дискретных точках i, i I ( 1). Поэтому эти точки будем использовать для аппроксимации функций z = z ().

Если условие (2.126) выполняется, а условие (2.127) не выполня ется, то переходим к шагу 6.

Если условия (2.126), (2.127) выполняются, то решение найдено a = a ( ), d = d ( ), z = z i,( ) и алгоритм останавливает свою работу.

Шаг 6. Вычисляем ( ) 2 (a ( ), d ( ) ) = max min max g j a ( ), d ( ), u ( ), z, (2.128) jJ z с использованием алгоритма внешней аппроксимации где J1 = (0, 1, 2, …, m1). Обозначим через ( ) решение задачи (2.127) и дополним точ кой ( ) множество точек S2k ), в которых нарушаются мягкие ограни ( чения (2.128), т.е.

S 2 ) = S 2 1) U ( ), ( ( I 2 ) = I 2 1) U (n + 1) ;

n := n + 1.

( ( S ( ) = S1 ) U S 2(), ( Шаг 7. Сформируем множества I ( ) = I1( ) U I 2(), положим := + 1 и перейти к шагу 2.

Дадим некоторые пояснения алгоритму.

На шаге 5 осуществляется многомерная интерполяция с помощью функций z = z () по известным дискретным точкам i, z i, i I ( ).

Это можно сделать с помощью многомерных кубических сплайнов или с использованием процедуры приближенной аппроксимации, суть ко торой заключается в следующем. При реализации имитационной мо дели для каждого полученного случайного значения в качестве со ответствующего z () берем значение z l ( l ), l I ( ), которое соот ветствует точке i, наиболее близкой к точке, т.е.

n ( j ij ) 2, i I () = I1() U I 2(), r i (, i ) = n = dim, j = ) = min r i (, l ) i = arg min r i (, i ) z = z i.

iI ( k ) iI ( k ) Фактически в описанной процедуре мы используем здесь кусоч но-постоянную аппроксимацию функций z = z().

На шаге 6 неравенство 2 (a ( ), d ( ) ) 0 означает, что мягкие ог раничения выполняются с вероятностью 1. Поэтому, если не выполня ется условие (2.127), то заведомо 2 (a ( ), d ( ) ) 0 и, следовательно, мы получим точку (k ), в которой нарушаются мягкие ограничения.

При использовании дополнительной переменной u проводим масштабирование поисковых переменных, чтобы диапазоны их изме нения были примерно одинаковы.

2.6. ДИНАМИЧЕСКАЯ ОПТИМИЗАЦИЯ ХТС В УСЛОВИЯХ НЕОПРЕДЕЛЕННОСТИ При интегрированном проектировании ХТС осуществляется со вместное проектирование ХТС и системы автоматического управления [4, 54]. Дополнительное усложнение вносится наличием неопределен ных параметров в математических моделях [55]. Это требует рассмот рения проблемы динамической оптимизации ХТС в условиях неопре деленности. По аналогии с проблемой статической оптимизации тре буется ввести понятие динамической гибкости, одноэтапной и двух этапной задач оптимизации в условиях неопределенности.

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

Сформулируем задачу оптимальной стабилизации для класса ра зомкнутых систем управления: требуется найти управление u * (t ) U, доставляющее минимум функционалу качества вида:

( )( ) I (u * (t ), ) = min y (t k ) y зад (t k ), F y (t k ) y зад (t k ) + uU t ( (y(t ) yзад (t )), Q(t ) (y(t ) yзад (t )) + ] u (t ), G (t ) u (t ) )dt (2.129) 1k + 2t при связях в форме уравнений математической модели динамики не линейного химического процесса х(t ) = f ( х(t ), u (t ), );

& х(t0 ) = х0 ;

(2.130) и ограничениях на качественные показатели переходного процесса в системе автоматического управления.

Здесь х(t ), хзад (t ) – п-мерные векторы текущего и заданного состояния (программы изменения) процесса, соответственно;

u (t ) – т-мерный вектор управления;

f (•) – нелинейная по х и u вектор функция;

F, Q(t ) – положительно полуопределенные матрицы (n n) ;

G (t ) – симметричная положительно определенная матрица (m m) ;

• – скалярное произведение векторов.

Для решения задачи (2.129), (2.130) нами применялся метод «по следовательной итерации», суть которого состоит в замене исходной нелинейной задачи сходящейся последовательностью линейных [56].

Каждая линейная задача последовательности получается путем линеа ризации нелинейной вектор-функции f (•) в окрестности траектории состояния ХТП и управления, полученных при решении предыдущей линейной задачи. В первом приближении функция f (•) линеаризует ся в окрестности траектории x(t ) = xзад (t ), u (t ) = u (0) для задачи опти мальной стабилизации. В этом случае система линеаризованных урав нений имеет вид x ( +1) (t ) = A( ) x ( +1) (t ) + B ( ) (t )u ( +1) (t ) + h ( ) (t );

& x ( +1) (t 0 ) = y 0 ;

(2.131) f ( x, u, ) A( ) (t ) = ;

x x = x ( ) (t ) u =u ( ) ( t ) где f ( x, u, ) B ( ) (t ) = ;

u x = x( ) (t ) u =u ( ) ( t ) h ( ) (t ) = f ( x ( ), u ( ), ) A( ) (t ) x ( ) (t ) B ( ) (t )u ( ) (t ).

Задача (2.129), (2.131) линейна по переменным x((t)+1), u ((t)+1) и ее решение определяется известным соотношением для оптимального управления [50]. Последовательность линейных задач решается до тех пор, пока при некотором = q выполняется неравенство x ( q ) x ( q 1).

При этом вектор управления u (q ) принимается в качестве реше ния задачи (2.129), (2.131), т.е. u* = u ( q ). Заметим, что сходимость ите раций в сильной степени зависит от удачного выбора начального при ближения u ( 0).

Решение задачи синтеза оптимального управления в замкнутой системе может быть получено на базе метода АКОР по критерию обобщенной работы А.А. Красовского [18]. В соответствии с этим ме тодом для процесса, описываемого уравнениями r ij ( x, t )u j xi + f i ( x, ) = (i = 1, n), & j = оптимальными в смысле минимума функционала u 2 + u 2 оп tk t m 1k j j I = V3 x(t k ) + Q( x, t )dt + dt, k 2t j =1 j t0 являются управления V n u j = u j оп = k 2 k j ( x, t ), j xk k = где V = V ( x, t ) – решение уравнения V V n fi = Q t i =1 xi при граничном условии Vt =tk = Vз, f j, ij, Q, Vз – заданные непрерыв ные функции, k 2 0 – заданные коэффициенты.

j На решении уравнений свободного движения x M + f ( x M,, t ) = & (2.132) левая часть дифференциального уравнения обращается в полную про & изводную по времени: V = Q.

tk Отсюда следует: V ( x M (t k )) V ( x M (t 0 )) = Q ( x M (t )) dt.

t По условию для терминальной задачи V ( xM (t k )) = Vз ( xM (t k )). Та tk ким образом, V ( xM (t k )) = Vз ( x M (t k )) + Q( xM, t ) dt.

t Допустим, что текущее время и интервал оптимизации разбиты на достаточно короткие циклы длиной t ц. Начало очередного цикла с точностью до t ц совпадает с текущим моментом t. Предположим, что в начале каждого цикла система контроля и оценивания реального управляемого процесса определяет вектор состояния x(t ) и задает его в качестве начального значения в модель (2.132) свободного движения, обеспечивая в начале каждого цикла равенство x M (t ) = x(t ). Таким образом, интегрируя уравнения (2.132) свободного движения на ин тервале от t до tk, можно вычислить tk V ( x(t )) = Vз ( xM (t k )) + Q( xM, t )) dt. (2.133) t Однако оптимальные управления рассчитываются по формулам:

V n u j оп = k 2 k j, j = 1, r, (2.134) j xk k = V и конечной целью является вычисление частных производных.

x Точнее, как видно из (2.134), требуется определить r скалярных про V V x,..., x на векторы 1 j,..., nj.

изведений вектора градиента 1 n Обычно число управлений r меньше размерности пространства со стояний n, и выгодно сразу определять проекции вектора градиента ( ) V на 1 j,..., nj, а не на координатные оси т.е.. Применим для xk вычисления компонент и проекций градиента схему правой разности.

В результате получим выражение для расчета оптимальных управле ний в виде k 2 tk Vз ( x M (t k )) + Q dt j u j oп = j xM (t ) = x + j t tk Vз ( xM (t k )) + Q dt, j = 1, r, (2.135) xM = x ( t ) t где j – вектор (столбец) с компонентами 1 j,..., nj ;

– норма это го вектора;

– малая действительная величина.

Заданная функция Vз и квадратура в квадратных скобках вычис ляются на траекториях свободного движения объекта (2.135), возбуж даемого начальными условиями, которые для первой скобки соответ ствуют вектору x(t ) + j, для второй скобки – x(t ). Для определения значений всех r управлений согласно (2.138) проводим r + «запуск» прогнозирующей модели (2.135). Модель свободного движе ния объекта можно заставить работать в ускоренном времени, вводя масштаб по времени = t /, где = const 1. Тогда уравнения про гнозирующей модели имеют вид dxM + f ( xM,, ) = 0.

d Темп интегрирования, характеризуемый величиной, должен быть таков, чтобы за каждый цикл t ц осуществлялось достаточное число «прогонок» свободного движения на интервале t k t, необхо V димое для численного определения частных производных. В этом xi случае 2 =t 2 / Q( xM, ) d.

V ( x()) = Vз ( x M ( k )) + =t / Сформулированные согласно (2.135) управления подаются на объект x + f ( x,, t ) = u & и остаются неизменными в течение определенного цикла t ц.

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

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

Рассмотрим соответствующий вариант алгоритма оптимального управления с прогнозирующей моделью. Уравнения управляемого процесса записываются в виде xi + f i ( x1,..., xn, y1,..., yr, ) = 0, i = 1, n, & y j = u j, j = 1, r ;

& где y = ( y1,..., yn ) – вектор органов управления;

u = (u1,..., u r ) – вектор управления. Таким образом, в данном случае осуществляется управле ние скоростями перемещения органов управления.

Отмечая, что свободное движение воспроизводится прогнози рующей моделью в ускоренном времени, записываем dxM + f ( x M, y M, ) = 0, d dy M = 0.

d В начале каждого цикла t ц переменные состояния процесса вводятся в прогнозирующую модель. В данной задаче имеем расши ренный вектор состояния ( x, y ). Для численного определения частных V производных осуществляем варьирование начальных условий по y j y M в каждом запуске прогнозирующей модели.

Для преобразования данного алгоритма в алгоритм нетерминаль ного управления необходимо осуществить переход к скользящему ин тервалу оптимизации, при котором t k = t + T, где T – заданная длина интервала оптимизации. При этом определение tk осуществля ется из условия достижения, например выхода целевых продуктов за данного (максимального) значения.

Наиболее трудоемкой операцией в алгоритме с прогнозирующей моделью является численное интегрирование уравнений свободного движения, выполняемое в каждом цикле r + 1 раз. Предположим, что для численного интегрирования с достаточной точностью уравнений свободного движения (2.132) на начальном интервале оптимизации t k t необходимо M операций. Интервал интегрирования с каждым тактом сокращается (при фиксированном моменте времени tk ) и сред нее число операций однократного интегрирования (при достаточно большом числе циклов) будет равным 0,5M. Обозначим общее число циклов, на которые разбит интервал оптимизации (t k t ), через nц.

Тогда общее число операций, необходимых для решения задачи синте за оптимального управления, выражается формулой ~ nц (r + 1) M.

Выше были рассмотрены алгоритмы решения задач оптимальной стабилизации режимов и синтеза оптимального управления ХТП при фиксированном значении вектора неопределенных параметров. Да лее предполагается, что эффективность функционирования описанного выше алгоритма управления будет учитываться с учетом действия случайных возмущений. При этом возможна коррекция оптималь ного управления для тех областей изменения неопределенных пара метров оптимального управления, для которых эффективность функционирования системы управления ХТП окажется недостаточно эффективной. В этом случае необходимо формулировать и решать за дачу проектирования системы и алгоритма управления в условиях не определенности. Для таких систем в работе [57] рассматривается по становка задачи оптимизации в виде t k I (u ) = M Q ( x(t, ), u (t )) dt min (2.136) uU t0 при связях x = f ( x, u (t ), ), t [t0, t k ], & x(0) = x0 (2.137) и ограничениях { } Bep g j ( x(t, ), u (t ) j, j J, (2.138) где x E n – вектор состояния ХТП;

E L – вектор неопределенных (случайных) параметров, имеющих известную функцию плотности распределения P() с некоторой областью определения ;

f (•) – п-мерная вектор-функция, имеющая li -ю непрерывную производную по yi, i = 1, n, и m j -ю производную по j, j = 1, L;

u (t ) – вектор функция управления из некоторого допустимого множества U ;

j – заданные действительные числа;

– заданные вероятности;

Q(•), g (•) – достаточно гладкие функции своих аргументов.

Будем считать, что f (•) аналитична по x в E n и по вектору па раметров в области. Тогда решение x (t, u (t ), ) также будет ана литично в и может быть представлено сходящимся рядом Тейлора:

x ( q ) (t, u (t ), ) n x ( q ) (t, u (t ), ) = x ( q ) (t, u (t ), 2 ) + ( i ) + 2(i ) i = 1 x (t, u (t ), ) L 2 (q) + (i ) ( j ) +...;

q = 1, n, (2.139) (i ) ( j) i, j = где x (q ) – q-я координата вектора решения системы дифференциальных уравнений состояния ХТП (2.137);

( ) – случайный вектор отклонений от номинального (среднего) значения вектора неопреде ленных параметров = P() d ;

обозначим производные решения x(t, u (t ), ) по координатам вектора через x ( qi)) и предположим, ( ( j ) что подынтегральная функция Q( x(t, ), u (t )) имеет все непрерывные частные производные по x порядка n m. Если теперь воспользоваться свойством линейности оператора математического ожидания, то мож но получить для функционала I (u ) новое выражение { } tk {Q( x(t, ) + Q( x(t, )) x (i ) n L ~ I (u ) = M ( j ) + x (i ) ( j) i =1 j = t { } 1 2Q( x(t, )) x ( q ) x (i ) L n L n + M ( j ) ( r ) + x (i ) x ( q ) ( r ) ( j ) r =1 i =1 j =1 q = { } } 1 Q( x(t, )) 2 x (i ) L n L + M ( j ) ( r ) +... + R(Q ( m ) ) dt.

x ( i ) ( j ) ( r ) r =1 i =1 j = (2.140) { } Члены вида M ( (1) ) n1... ( ( L) ) nL, где nr, r = 1, L – некото рые степени, являются центральными моментами соответствующего порядка и легко считаются, поскольку предполагается известной функция P() – плотность распределения случайного вектора не определенных параметров. Однако, если случайная величина (i ) при каждом значении других случайных величин ( j ), j = 1,..., i 1, i + 1,..., L имеет симметричное распределение, то { } член M ( (1) ) n1... ( ( L ) ) nL равен нулю, как только (i ) входит в него в нечетной степени. Так, например, если вектор имеет сим метричное относительно своего математического ожидания распреде ление (нормальное, равномерное и т.п.), то ненулевыми будут лишь { } такие M ( (1) ) n1... ( ( L ) ) nL, в которых все степени nr четные.

Если теперь ограничиться конечным числом членов в разложении ~ (2.140), то мы получим детерминированный функционал I (u ).Таким образом можно перейти от стохастической задачи к детерминирован ~ ной задаче минимизации функционала I (u ) с уравнениями связи (2.137).


Далее в терминах теории А-задач стохастической оптимизации исходную задачу (2.139) – (2.141) с учетом (2.139), (2.140) можно за писать в виде ~ ~ I (u A* ) = min {min I ( x(t, ), u ), x = f ( x, u, ), g j ( x(t, ), u, t ) A j (t ), & A uU (2.141) j J, t [t 0, t k ]}, где { } [ ] = j, Bep g j ( x(t, ), u A, t ) 0, j J, t [t 0, t k ].

Для определения вектора A* (t ) можно построить минимизирую щую последовательность А-задач оптимизации: {AS } : m( AS ) m( A* ), ~ где m( A) = min I (u ).

uU Алгоритм ее построения сводится к следующему.

Шаг 1. Задается начальное приближение при = 0 для вектор функции, A(t ) – некоторые кусочно-постоянные функции ~ ( ) ~ () A j (ti ), j J, A j (t ), A j (t ) = A j (ti ) аппроксимирующие при t [t i, t i +1 ], i = 0, k 1.

Шаг 2. Методами АКОР или Крылова–Черноусько [58] решается задача оптимального управления, записанная в фигурных скобках (2.141).

Шаг 3. При найденном управлении u A (t ) производится проверка выполнения вероятностных ограничений (2.138) на каждом отрезке t [ti, ti +1 ], i = 0, k 1.

Шаг 4. Если вероятностные ограничения на каком-либо проме ~ жутке (t i, t i +1 ) не выполняются, т.е. A ( ) (t ), то для нарушенных ~ ограничений проводится уменьшение значений A ( ) (ti ) при t [t i, t i +1 ], i = + 1, и следует переход к шагу 2. В противном случае переходят к шагу 5.

~ Шаг 5. Кусочно-постоянные вектор-функции A ( ) (t ), j J, и j соответствующие им управления u A( ) находятся из решения задачи:

I (u A* ) = min {•}.

~ (2.142) A В общем случае задача (2.142) при условии кусочно-постоянной аппроксимации функции A(t ) может быть решена эффективным ме тодом последовательного квадратичного программирования.

Таким образом, все описанные выше алгоритмы оптимальной стабилизации и оптимального управления для детерминированных и стохастических объектов управления будут использоваться нами при интегрированном проектировании ХТП, аппаратов и систем.

3. ПРИМЕРЫ ИНТЕГРИРОВАННОГО ПРОЕКТИРОВАНИЯ АППАРАТУРНО-ТЕХНОЛОГИЧЕСКОГО ОФОРМЛЕНИЯ ГИБКИХ ХИМИКО-ТЕХНОЛОГИЧЕСКИХ ПРОЦЕССОВ И СИСТЕМ 3.1. ЭНЕРГО- И РЕСУРСОСБЕРЕГАЮЩЕЕ НЕПРЕРЫВНОЕ ПРОИЗВОДСТВО АЗОПИГМЕНТОВ Азопигменты получают при последовательном проведении процес сов тонкого органического синтеза (диазотирования и азосочетания), фильтрования, сушки и диспергирования пигмента. Процессы тонкого органического синтеза – химические реакции диазотирования и азосоче тания осуществляются в действующих производствах периодического действия в емкостных аппаратах большого объема с механическими перемешивающими устройствами, в которых невозможно обеспечить достаточную степень однородности полей концентраций и температуры в аппарате, а следовательно, и требуемое качество органических полу продуктов и красителей [59].

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

3.1.1. МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ НЕПРЕРЫВНОГО ПРОЦЕССА ДИАЗОТИРОВАНИЯ АРОМАТИЧЕСКИХ АМИНОВ Обзор литературных данных и результатов экспериментальных ис следований, проведенных в работах [60 – 67], позволил установить пе речень наиболее вероятных реакций, протекающих при синтезе азокра сителей, кинетические уравнения и константы химических процессов диазотирования (табл. 3.1) и азосочетания, математические модели этих процессов, осуществляемых в многоступенчатых проточных аппаратах смешения. В меньшей степени были исследованы непрерывные процес сы диазотирования и азосочетания, осуществляемые в турбулентных аппаратах трубчатого и комбинированного типов. В последние годы интенсивно развивается новый энерго- и ресурсосберегающий непре рывный способ получения целого ряда химических продуктов с исполь зованием малогабаритных трубчатых турбулентных аппаратов [68 – 71].

На рисунках 3.1, 3.2 представлены трубчатые турбулентные аппа раты, позволяющие осуществлять химические процессы в высокотур булентных потоках.

III II III II III II III II IV I а) D B dk / dd /2 IV I A C б) Рис. 3.1. Малогабаритный турбулентный трубчатый реактор:

а – цилиндрического типа;

Рис. 3.2. Турбулентный б – диффузор конфузорного типа: трубчатый реактор dd – диаметр диффузора, комбинированного типа:

I, II, III, IV – потоки веществ dk – диаметр конфузора 3.1. Перечень наиболее вероятных реакций, протекающих при диазотировании ароматических аминов нитритом натрия, и уравнения кинетики процесса диазотирования Диазотирование (реакция экзетермическая) растворение твердой фазы амина в среде соляной кислоты [ArNH2]s ArNH W образование диазотирующего агента (HNO2) NaNO2 + HCl W = HNO2 + NaCl целевая реакция диазотирования ArN2Cl + 2H2O + q (Дж) ArNH2 + HNO2 + HCl W разложение азотистой кислоты 3HNO2 W 2NO + HNO3 + H2O образование диазосмол ArN2Cl + HNO2 W разложение диазосоединения с образованием диазосмол ArN2Cl W образование диазоаминосоединений ArN2Cl + ArNH2 W 6 Ar2N3H + HCl Продолжение табл. 3. Уравнения кинетики и кинетические константы процесса диазотирования Предэкспоненциальный Энергия Порядок Кинетическое уравнение множитель, активации, реакции (м3)n–1/мольn–1/с Дж/моль 3,75105 46, W2 = k2 [ArNH2] [HNO2] W3 = k3[HNO 2 ]4 /PNO 7,171021/(9,81104)2 119, 0,32105 63, W4 = k4 [ArN2Cl] [HNO2] 1,101010 87, W5 = k5 [ArN2Cl] О б о з н а ч е н и я: А – ArNH2 (амин);

АК – HNO2;

СК – HCl;

– NO;

D – ArN2Cl (диазосоединение);

= (1, 2);

PNO – парциальное давление нитрозных газов;

r – радиус частицы амина;

l – текущая координата длины реактора;

cA – равновесная концентрация амина;

A, M A – плотность и * молекулярная масса амина;

( r, l ) – ненормированная плотность распре деления числа N частиц амина по размеру, т.е. = dN / dr ;

Gl, GS, GN – расходы жидкой и твердой фаз суспензии амина и нитрита натрия;

S тр – площадь поперечного сечения трубы реактора;

c p – теплоемкость;

– скорость потока;

K – коэффициент теплопередачи;

q – тепловой эффект реакции;

индексы: (0) – вход в реактор;

х – хладагент;

(i) – точка распреде ленного по длине ввода реагентов в реактор.

Такие аппараты весьма эффективны не только для проведения химических процессов, но и для интенсивного смешения компонентов с целью гомогенизации жидкости. Они успешно используются в каче стве смесителей (для получения целевого продукта), предреакторов (для приготовления исходной реакционной смеси), малогабаритных реакторов (для проведения быстрых химических реакций), а также осуществления массообменных процессов [71 – 74].

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

1) смешения дозируемых в аппарат исходных реагентов в камере сме шения, установленной на входе в аппарат;

2) химического взаимодей ствия ароматического амина с диазотирующим агентом (азотистой ки слотой) в трубчатой части аппарата и в камерах смешения (в аппаратах комбинированного типа);

3) изменения условий протекания элементар ных процессов химического взаимодействия при наличии в аппарате устройства турбулизации потока (например, диффузор-конфузора).

Рассмотрим основные этапы построения математической модели процесса диазотирования, осуществляемого в турбулентных аппаратах трубчатого типа в среде соляной кислоты. Исходное сырье – ароматический амин (например, 3-нитро-4-амонотолуол) – непрерыв но подается в реактор в виде солянокислой суспензии. Одновременно в реактор вводится водный раствор нитрита натрия, подача которого может распределяться по длине трубчатой части реактора или по сек циям комбинированного реактора (см. рис. 3.1, 3.2).

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

Допустимый диаметр трубы реактора по условию неосаждения агрегатов твердой фазы амина определяется как d d кр, где d кр – первый критический диаметр трубы реактора рассчитывает ся по формуле Re µc d кр =, (3.1) Koc c где µ c – средний динамический коэффициент вязкости суспензии;

oc – скорость осаждения агрегатов амина;

c – средняя плотность суспензии;

K = п / oc, п – скорость потока.

Средняя плотность суспензии при заданной плотности твердой фазы T и массовая доля кислоты в водном растворе X CK рассчиты ваются по формулам:

1 c = = ;

;

X T 1 X T X CK 1 X CK + + T CK B сA ) M A ( сCKM CK X CK = = (3.2) ;

, (1 ) T где X CK, X T – массовые доли концентрированной соляной кислоты в растворе и твердой фазы амина в суспензии;

T, CK, B – плотности твердой фазы амина, концентрированной соляной кислоты и воды со ответственно, cA0), cCK = 2,5cA – концентрация амина и соляной ки ( (0) слоты в суспензии на входе в реактор.

Из условий надежной работы клапанов и насосов, питающих ре актор, концентрация амина в суспензии на входе в реактор не должна превышать 800 моль/м3. Отсюда максимальная объемная доля твердой фазы амина в суспензии будет равна = 8,4%, и этому значению соответствует следующая формула для расчета вязкости суспензии:


µ c = µ (1 + 2,5), (3.3) где µ c =, µ, µ CK, µ B – средние динамические коэф X CK 1 X CK + µ CK µB фициенты вязкости жидкой фазы, соляной кислоты и воды соответст венно.

Среднюю скорость свободного осаждения твердой фазы можно рассчитать по формуле d T ( T ) max oc = 5,46. (3.4) После подстановки формул (3.2) – (3.4) в выражение (3.1) полу чим, что d кр зависит от концентрации амина в суспензии на входе в реактор, т.е.

d кр = f (cA ) ).

( (3.5) Вторым условием работоспособности реактора диазотирования является обеспечение турбулентного режима течения суспензии амина – Re = 10 4. Определяемый с учетом этого условия диаметр будем нахо дить из следующего выражения:

4Gl(0) c d кр =, Re µc где d кр – второй критический диаметр труб реактора;

Gl(0) – расход жидкой фазы на входе в реактор.

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

Q 10 Gl(0) =, cA )Tэфф K D K Az М П ( где K D, K Az, Tэфф – выходы по стадиям диазотирования, азосочетания и эффективный фонд рабочего времени реакторной установки синтеза азопигментов соответственно;

М П – молекулярный вес пигмента.

Диаметр d кр зависит от концентрации амина в питании реактора и производительности d кр = f (cA ), Q), ( (3.6) где Q – производительность реактора по пигменту.

На рисунке 3.3, а представлены графики, соответствующие урав нениям (3.5) и (3.6), из которых видно, что при производительности меньше 1000 т/год диаметр трубы реактора определяется условием турбулентности.

Задаваясь максимально возможным диаметром трубы реактора (40 мм), можно определить максимальную концентрацию амина с A0) ( в питании реактора. Минимальную концентрацию амина можно опре делить по точке перегиба кривой для каждой производительности Q (см. рис. 3.3, а). Соответственно такими концентрациями будут:

– для производительности 500 т/год cA ) [130, 210] моль/м3;

( – для производительности 1000 т/год cA ) [280, 410] моль/м3.

( На рисунке 3.3, б представлены графики зависимости среднего времени пребывания реакционной смеси от концентрации ами на с A0). Их анализ показывает, что при производительности больше ( 1000 т/год и концентрации амина не более 100 моль/м3 время пребыва ния будет расти, так как определяющим условием при этом является неосаждение твердой фазы амина (рис. 3.3, а). При концентрациях больше 100 моль/м3 и Q 1000 т/год диаметр реактора определяется только из условий турбулентности реакционного потока и существен но зависит от изменения концентрации амина, что вызывает большее изменение объема реактора, нежели расхода суспензии на входе в ре актор. При Q 1000 т/год среднее время пребывания при всех концен трациях падает. Только при очень больших концентрациях (700…800) моль/м3 влияние на время пребывания расхода суспензии и объема реактора становится равнозначным.

На стадии функционирования ХТС диаметр трубы реактора по стоянен и при необходимости увеличения производительности Q можно наращивать его длину. Этого результата можно добиться и тех нологическим приемом, повышая концентрацию твердой фазы амина с целью снижения расхода солянокислой суспензии амина в питании реактора. Для увеличения концентрации амина без нарушения условий работоспособности реактора необходимо выбрать такой диаметр тру бы (см. рис. 3.3, а), которому будет соответствовать максимально воз можный интервал допустимых концентраций амина для выбранной производительности. Такому условию соответствует диаметр трубы, равный 40 мм.

/ 10, c d кр, м 0, 0, 0,43 1 0, 4 0,15 0,01 350 500 650 350 500 650 cA0), моль/м ( cA0), ( моль/м б) а) Рис. 3.3. Зависимости критических диаметров трубы реактора dкр (а) и среднего времени пребывания реакционной смеси (б) от концентрации амина на входе в реактор:

а – 1 – d кр = f (cA ) );

2 – 5 – dкр = f (cA ), Q);

2 – Q = 2000 т/год;

1 ( 1 ( 3 – Q = 1500 т/год;

4 – Q = 1000 т/год;

5 – Q = 500 т/год;

б – 1 – Q = 2000 т/год;

2 – Q = 1500 т/год;

3 – Q = 1000 т/год;

4 – Q = 500 т/год Вторым важным достоинством применения реактора с макси мальным диаметром (40 мм) является возможность его использования в малотоннажных производствах, так как при больших диаметрах не будет выполняться условие турбулентности реакционной смеси.

На рисунке 3.4 представлен интервал допустимых значений кон центраций амина на входе в реактор, по которому определяются ми нимально и максимально возможные значения концентраций амина при различных производительностях реакторной установки по пиг менту. На стадии эксплуатации реактора концентрация амина будет изменяться относительно номинального значения, ориентировочно, в пределах ±2,5%. В связи с этим на рис. 3.4 изображена допустимая область значений концентраций амина, ограниченная кривыми cA0max = f (Q) / 0,975 и cA0max = f (Q) / 1,025 (допустимая область за () () штрихована).

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

При меньших производительностях целесообразно использовать реак тор с мешалкой периодического действия или реактор смешения типа «царга-тарелка» [75, 76].

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

Схема материальных потоков в одной секции трубчатого реактора приведена на рис. 3.5.

cA ), моль/м ( cA0) ( max cA0min () 869 1246 1623 115 Q, т/год Рис. 3.4. Допустимая область концентрации амина на входе в реактор l Gx ( L), Tx ( L) cA ), cCK, Gl( 0), (0 ( 0) Хладагент GS0), (0, r ), T ( 0) ( Солянокислая r c ( L), ( L, r ) суспензия амина Раствор диазосоединения cD ), cAK, c0), (0 (0) ( l1 l cN0), GN0), TN0) ( ( ( Gx (0), Tx (0) Водный раствор нитрита натрия Рис. 3.5. Схема материальных потоков в трубчатом реакторе диазотирования:

(r ) – гранулометрический состав;

x – хладагент;

l – жидкая фаза;

S – твердая фаза;

, – нитрозные газы и диазосмолы соответственно При составлении уравнений математической модели примем сле дующие допущения: 1) гидродинамические режимы течения реакци онной смеси и хладагента близки к режиму идеального вытеснения;

2) реактор является объектом с распределенными только по длине l координатами, скорость потока (l ) = const;

3) потери тепла в окру жающую среду пренебрежимо малы;

4) теплофизические характери стики принимаются постоянными в рабочем диапазоне температур;

5) кристаллы амина будем считать сферическими частицами;

6) реак ция диазотирования происходит в растворе.

С учетом принятых обозначений и допущений составим уравне ния материального покомпонентного баланса на участке трубы ( l1, l 2 ) за промежуток времени ( t1, t2 ) по растворенному амину:

t [cA (l2, t ) cA (l1, t )]dt + Gl t t2 l +S W2 (cA, cAK, T ) (l, r )W1 (cA, T ) dr dl dt = t1 l1 MA l = S [с(l, t 2 ) с(l, t1 )]dl.

l Пользуясь теоремой о среднем, получим равенство Gt [cA (l 2, t ) cA (l1, t )] t + t = t + S W2 (cA, cAK, T ) (l, r )W1 (cA, T ) dr lt = MA l = l t =t = S [c(l, t 2 ) c(l, t1 )] l, l =l которое при помощи теоремы о конечных приращениях можно преоб разовать к виду c (l, t ) tl + Gl A l l = l t = t + S W2 (cA, cAK, T ) (l, r )W1 (cA, T ) dr lt = MA l =l t =t c(l, t ) =S lt.

t l = l t = t Переходя к пределу при l1, l 2 l и t1, t 2 t, получим уравнение cA c + W2 (cA, cAK, T ) (l, r ) W1 (cA, T ) dr = ;

(3.7) l t MA cA (l, 0) = cA (l );

cA (0, t ) = cA (t ).

(0) l Аналогичным образом можно получить уравнения динамики трубчатого реактора и для других компонентов реакционной смеси:

по азотистой кислоте (АК) cAK c + W2 (cA, cAK, T ) + W3 (cAK, T ) + W4 (cAK, c D, T ) = AK ;

(3.8) l t cAK (l, 0) = cAK (l );

cAK (0, t ) = c N G N0) (t ) / Gl(0) (t );

l ( по диазосоединению (D) c D c W2 (cA, cAK, T ) + W4 (cAK, c D, T ) + W5 (c D, T ) = D ;

(3.9) l t c D (l, 0) = c D (l );

c D (0, t ) = c D ) (t );

( по продуктам разложения (, ) c c W3 (cAK, T ) = ;

c (l, 0) = c (l );

c (0, t ) = c0) (t );

0 ( l t (3.10) c c W4 (cAK, c D, T ) W5 (c D, T ) = ;

(3.11) l t c (l, 0) = c (l );

c (0, t ) = c0) (t ).

0 ( Составим уравнение динамики для фракции частиц амина, харак теризующихся размером от r до r + dr на участке трубы (l1, l2 ) за время (t1, t2 ) :

t [N (l2, t )(l2, t, r )dr N (l1, t )(l1, t, r )dr ]dt + t [ ] tl 1 + N (l, t ) (l, t, r )W1( r ) (l, t, r ) (l, t, r + dr ) W1( r ) (l, t, r + dr ) dl dt = t l l [N (l, t 2 )(l, t 2, r )dr N (l, t1 )(l, t1, r ) dr ]dl, = l которое с использованием приведенной выше техники можно преобра зовать к уравнению вида [ ] (l, t, r ) (l, t, r ) (l, t, r )W1( r ) (l, t, r ) = ;

(3.12) l r t (0, t, r ) = (0) (t, r );

(l, 0, r ) = 0 (l, r ).

Получим теперь уравнения динамики теплообмена в трубчатом реакторе по реакционной смеси:

t c p [T (l 2, t ) T (l1, t )]d + t t2 l [ qSW2 (l, t ) + K1D[T (l, t ) Tx (l, t )]]dldt = + t1 l l = c p S [T (l, t 2 ) T (l, t1 )]dl ;

l по хладагенту (х):

t2 t 2 l [Tx (l2, t ) Tx (l1, t )]dt [K1D[T (l, t ) Tx (l, t )]]dldt = c x x Gx p t1 t1 l l = c x x S p [T (l, t 2 ) Tx (l, t1 )]dl.

p l Проводя рассуждения, аналогичные предыдущим, получим урав нения:

T (l, t ) T (l, t ) qSW2 (l, t ) + K1D [T (l, t ) Tx (l, t )] = c p S c p Gl ;

(3.13) l t T (l, 0) = T 0 (l );

T (0, t ) = T (0) (t );

Tx (l, t ) T (l, t ) K1D [T (l, t ) Tx (l, t )] = c x x S p x c x x Gx ;

(3.14) p p l t Tx (l, 0) = T0 (l );

Tx (0, t ) = Tx( L ) (t ).

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

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

c T = 0;

= 0;

= 0.

t t t Математическая модель статики процесса диазотирования, осу ществляемого в трубчатом реакторе, представляет собой систему не линейных обыкновенных дифференциальных уравнений. К выходным переменным процесса диазотирования помимо концентраций, темпе ратуры и расходов потоков относятся: производительность реакторной установки – Q ( a, d, z, ) = с D Glвых M D 28 1000 т пигмента/год;

вых с D Glвых вых выход диазосоединения – K D (a, d, z, ) = 100 98% ;

кон [с ] G вх вх Аs l вых центрация азотистой кислоты 7 моль/м3 ) 14 моль/м cАК (a, d, z, на выходе из реактора;

количество диазосмол в диазорастворе – с Glвых вых П ( a, d, z, ) = 100% 0,5% ;

количество нитрозных газов в []сА s Glвх вх с Glвых вых диазорастворе – П (a, d, z, ) = 100% 1,0%.

вх вх c N GN Наибольшую сложность при решении системы дифференциаль ных уравнений, описывающих статические режимы непрерывного процесса диазотирования, представляет уравнение, описывающее гра нулометрический состав твердой фазы амина в l-м сечении трубчатого реактора:

[ ], (l, r ) = (l, r )W1( r ) (l, r ) l r (0, r ) = (0) (r ).

(3.15) В случае линейного уравнения кинетики растворения частицы амина ) = Ar exp ( E1 / RT )(cA c A / A.

dr * (3.16) dt Решение уравнения (3.15) может быть получено методом харак теристик в аналитическом виде [77], где A, – кинетические констан ты;

cA, cA – равновесная и текущая концентрации амина;

A – плот * ность амина.

Запишем решение уравнения (3.16) в виде 1+ l A exp ( E1 / RT )(cA cA ) * ~ r (l ) = f (r0, l ) = r0 +1 (1 + ) 0 dl, A откуда можно рассчитать начальный радиус r0 частицы по формуле 1+ l A exp( E1 / RT ) (cA cA ) * ~ r0 = f1 (r, l ) = r +1 + (1 + ) 0 dl.

A В этом случае решение уравнения (3.15) с начальным условием может быть записано в следующем виде [77]:

[ ] l W ( r ) ~ ~ ~ (r, l ) = (0) ( f1 (r, l )) exp r ( l, f1 (r, l )) d l.

(3.17) 0 r В случае нелинейного уравнения кинетики растворения частицы, например, в виде dm = * (cA cA ) S, * dt где – эффективный коэффициент массоотдачи;

S – поверхность час тицы, необходимо использовать численный алгоритм решения уравне ния (3.15). Аппроксимируя дифференциальные уравнения в частных производных (3.15) конечной системой дифференциальных уравнений в обычных производных с использованием конечно-разностной схемы первого порядка, получим i 1 W d i r +r (r ) = W1( r ) i i 1, cA, cAK i (ri, cA, cAK ) i, (3.18) ri r dl 1 (0) = (0) (ri, l ), ri – шаг сетки.

Полученная система обыкновенных дифференциальных уравне ний (3.18) одновременно с другими уравнениями модели может быть решена каким-либо численным методом. При этом могут возникнуть сложности, поскольку в начальной фазе процесса диазотирования ско рость растворения твердой фазы и скорость реакции диазотирования различаются на несколько порядков, т.е. система дифференциальных уравнений математической модели процесса диазотирования является жесткой. В этом случае явные методы Рунге–Кутта исключаются из рассмотрения, и для решения системы дифференциальных уравнений модели статики процесса диазотирования мы применяли два метода:

неявный метод трапеций [78] и метод Дормана–Принса 5-го порядка точности с автоматическим выбором шага [79], которые дают вполне сопоставимые результаты и обеспечивают получение решения с за данной точностью.

Математическая модель динамики процесса диазотирования, осуществляемого в турбулентном трубчатом аппарате, включает нели нейные дифференциальные уравнения с частными производными (3.7) – (3.14), для решения которых использовали неявный конечно разностный метод [80].

Для математического описания процесса диазотирования, осуще ствляемого в турбулентном трубчатом реакторе комбинированного типа, необходимо к уравнениям (3.7) – (3.14) добавить уравнения, опи сывающие протекание процесса диазотирования в камере смешения.

В неустановившемся режиме текущий радиус частицы r кроме начального значения r0 зависит также от возраста частицы и теку щего времени t. Кинетика растворения частицы описывается квазили нейными дифференциальным уравнением с частными производными первого порядка:

r r + = W1 (r, cCK, cA, T ), (3.19) t где W1 – скорость растворения частицы;

cCK, cA – концентрации со ляной кислоты и амина соответственно;

T – температура.

Начальные условия определим следующим образом:

t 0 = 0, = 0, r ( 0 ) = r0 ( 0 ). (3.20) Решение уравнений (3.19) при заданных начальных условиях (3.20) имеет вид r ( +1) (, t ) r0 +1 ( 0 ) t ~* ~~ A exp ( E1 / RT ( t )(cA (cCK, T ) cA ( t ))d t, = +1 +1 где t ~* ~ ~ ~~ W1 (r, cCK, cA, T ) = A exp( E1 / RT ( t )(cA (cCK ( t ), T ( t )) cA ( t ))d t.

Учитывая, что 0 = t t 0, получим 0 = t + t 0 и r ( +1) (, t ) r0 +1 ( t ) t ~* ~ ~ ~~ = A exp ( E1 / RT ( t )(cA (cCK ( t ), T ( t )) cA ( t )) d t.

+1 +1 Плотность распределения массы частиц в каждый момент време ни описывается дифференциальными уравнениями, получаемыми из материального баланса для фракции частиц:

(r, t ) [ W1 ] ( 0) ( 0) = + n (r, t ) (r, t ) / (t ), (3.21) t r где (t ) – среднее время пребывания частиц в модуле-реакторе.

Решение уравнения (3.21) при тех же начальных условиях (3.20) будет следующим:

(r, t ) = exp [F (0, t, r, cCK, c A, T ) { 0 [(0, t, r, cCK, cA, T )] }] + t [ ] [ ] + Gs 0) (t, t, r, cCK, c A, T ), t exp F (0, t, r, cCK, cA, T ), cCK, T d t, ~ ~ ~ ~ ( (3.22) где t {W1r [(t, t, r, cCK, cA, T ), cCK, cA, T ]+ 1 / ( t )}d t, ~ ~ ~ F (0, t, r, cCK, cA, T ) = W1r = dW1 / dr = Ar ( +1) exp ( E1 / RT )(cA (cCK, T ) сA ).

* Подставляя зависимости (3.22) в уравнения покомпонентного ма териального и теплового балансов, получим rmax (r, t )W1(r, cCK, cA, T )dr /V + cA Gi cA = / V W2 (cA, cAK, T ) cAGl ;

( 0) ( 0) & с A (0) = c A 0 ;

(3.23) cAK = cAKGl( 0) V + cN0)GN0) / V (W2 (cA, cAK, T ) + ( 0) ( ( & + W3 (cAK, T ) + W4 (cA, cD, T ) cAKGl V ;

сAK (0) = cAK 0 ;

(3.24) cCK = cCK Gl(0) V c N0)G N0) / V W2 (cA, cAK, T ) cCK Gl V ;

( 0) ( ( & сCK (0) = cCK 0 ;

(3.25) cD = cD ) Gl(0) V + (W2 (cA, cAK, T ) W4 (cAK, c D, T ) W5 (cD, T )) c DGl V ;

( & с D (0) = c D 0 ;

(3.26) c = c0) Gl( 0) / V + W4 (c AK, c D, T ) + W5 (c D, T ) c Gl / V ;

( & с (0) = c 0 ;

(3.27) & T = Gl( 0)T (0) / V + cv GN0)TN0) / cvV GlT / V + N( ( W2 (cA, cAK, T )h2 W3 (cAK, T )h + + KF (T Tx ) / cvV ;

T (0) = T0 ;

(3.28) cvV cvV & Tx = Gx0)Tx(0) / V p GxTx / V p + KF (T Tx ) / cv V p ;

( x Tx (0) = Tx 0. (3.29) Система нелинейных дифференциальных уравнений (3.22) – (3.29) представляет собой математическую модель динамики процесса диа зотирования, осуществляемого в камере смешения, которую мы будем использовать в дальнейшем для получения линеаризованных описаний динамики в окрестности номинального статического режима и провер ки их адекватности в различных диапазонах изменения возмущающих воздействий.

Линеаризуем уравнения (3.23) – (3.29) в окрестности номинально го установившегося режима, обозначаемого *, и получим x = Fx + Gu, & (3.30) r где x T = {, c, T, Tx }, u T = {Gl(0), G N0), Gx0) }.

( ( c, T, Tx, Gl(0), GN0), Gx0) обо ( ( Векторы с компонентами значают отклонения от стационарных значений, т.е. ci = ci ci*, T = T T *, Tx = Tx Tx*, Gl(0) = Gl(0) Gl(0)*, GN0) = GN0) GN0)*, ( ( ( Gx = Gx0) Gx0)*.

( ( Элементы матриц F, G системы дифференциальных уравнений (3.30) имеет следующий вид:

rmax W1 (r, cCK, cA, T f11 = * * * )dr / V ;

rmax W f12 = Gl(0)* GN0)* V+ * (r, t ) (r, cCK, cA, T * )dr / V k2 (T * )cAK ;

( * * * V cA f13 = k2 (T * ) cA ;

* rmax W f14 = * (r, cCK, cA, T * )dr / V ;

f15 = 0;

f16 = 0;

* * cCK rmax W1 k (T * ) * * f17 = * (r, cCK, cA, T * )dr V * * cA cAK ;

T T f18 = 0;

f 21 = 0;

f 22 = k 2 (T ) c AK ;

f 23 = Gl( 0) / V G N0) / V k 2 (T * )cAK 4k 3 (T * )cAK k 4 (T * )c * ;

( * * D f 24 = 0;

f 25 = k 4 (T )cAK ;

f 26 = 0;

* k (T * ) *4 k 4 (T * ) * * k 2 (T * ) * * f 27 = cA cAK 3 cAK cAK c D ;

T T T f 28 = 0;

f 31 = 0;

f 32 = k 2 (T * ) cAK ;

* f 33 = k 2 (T * ) cAK Gl(0)* / V G N0)* / V ;

( * k 2 (T * ) * * f 34 = f 35 = f 36 = 0;

f 37 = c A c AK ;

T f 38 = f 41 = 0;

f 42 = k 2 (T * ) cAK ;

f 43 = k 2 (T * ) cA k 4 (T * ) c * ;

f 44 = 0;

* * D f 45 = Gl(0)* / V k 4 (T * ) c AK k 5 (T * );

f 46 = 0;

* k (T * ) * * k 5 (T * ) * k 2 (T * ) * * f 47 = cA c AK 4 cAK c D cD ;

T T T f 48 = f 51 = f 52 = 0;

f 53 = k4 (T * )c* ;

f 54 = 0;

D f 55 = k 4 (T * ) cAK + k 5 (T * );

f 56 = Gl(0)* / V G N0)* / V ;

( * k4 (T * ) * * k5 (T * ) * f57 = cAK cD + cD ;

f 58 = 0;

T T f 61 = 0;

f 62 = k 2 (T * ) cAK h2 / cvV ;

f 63 = k 2 (T * ) c A h2 / cvV ;

* * f 64 = f 65 = f 66 = 0;

k 2 (T * ) * * f 67 = Gl(0 ) / V G N0) / V + cA cAK h2 / cvV KF / cvV ;

( T f 68 = KF / 2cvV ;

f 71 = f 72 = f 73 = f 74 = f 75 = f 76 = 0;

f 77 = KF / cv V p ;

f 78 = Gx0) / V p KF / 2cv V p ;

x ( x g11 = (cA0)* cA ) / V ;

g12 = cA / V ;

g13 = 0;

( * * g 21 = (cAK cAK ) / V ;



Pages:     | 1 |   ...   | 2 | 3 || 5 | 6 |   ...   | 8 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.