МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ УСТОЙЧИВОСТИ БИОТЕХНОЛОГИЧЕСКОГО ПРОЦЕССА В ЕГО ОПТИМАЛЬНОМ СОСТОЯНИИ
Аннотация и ключевые слова
Аннотация (русский):
Практическая реализация процессов микробиологического синтеза в непрерывных условиях функционирования связана с устойчивостью функционирования в стационарном состоянии. Прогнозирование показателей процесса, обеспечивающих устойчивое функционирование, осуществляется с использованием математического моделирования, методология которого включает оценку показателей стационарного состояния и исследование характера изменения этих показателей при возможных возмущениях в процессе. Разработана методология оценки устойчивости стационарных состояний анаэробного процесса микробиологического синтеза в оптимальных условиях функционирования. Методология базируется на использовании линеаризованной системы уравнений нестационарного состояния, записанной в отклонениях. Получены соотношения для расчета коэффициентов линеаризованных уравнений, система уравнений сведена к одному дифференциальному уравнению третьего порядка, составлена матрица Гурвица и записаны необходимые и достаточные условия устойчивости. Однако для оптимального стационарного состояния получено дифференциальное уравнение второго порядка, для которого необходимые и достаточные условия устойчивости отвечают требованию положительности коэффициентов характеристического уравнения. Система уравнений в отклонениях упрощена с учетом оптимальности, рассмотрены возможные варианты расчетов, ориентированные на вид корней характеристического уравнения. Приведены аналитические решения для нестационарных условий в окрестности стационарного состояния. Рассмотрен пример численной реализации методологии в окрестности стационарного состояния с использованием известных кинетических соотношений, характерных для многих процессов микробиологического синтеза. В качестве критерия оптимальности принята продуктивность по целевому продукту. Расчеты для численных значений кинетических констант показали, что оптимальное состояние является устойчивым. Разработанный методологический подход может быть рекомендован для анализа других кинетических схем микробиологических процессов.

Ключевые слова:
биотехнология, моделирование, устойчивость
Текст
Введение Одно из направлений совершенствования процессов микробиологического синтеза связано с созданием оптимальных условий, определяющих максимальный выход целевых продуктов. Решение этой задачи не может быть реализовано только на базе экспериментальных исследований. Моделирование оптимальных условий осуществляется с использованием кинетических соотношений, вид которых все более усложняется [1]. В то же время обеспечение устойчивого функционирования процессов требует использования систем управления, поддерживающих необходимый режим работы ферментера при наличии возможных возмущений, приводящих к отклонению показателей процесса от регламентных. Понятие устойчивости является чрезвычайно важным в технологическом обеспечении микробиологических процессов, организованных по непрерывной схеме. В случае устойчивого стационарного состояния существуют возможности так называемого саморегулирования, т. е. восстановления регламентных значений показателей процесса после возмущения без вмешательства системы управления. Для оценки условий устойчивости необходимо использовать математический аппарат, достаточно хорошо разработанный теоретически. В настоящем сообщении сформирована методология оценки условий устойчивости для микробиологического процесса, функционирующего в оптимальных условиях. Уравнения математической модели процесса микробиологического синтеза в нестационарном состоянии в непрерывных условиях функционирования имеют следующий вид [2, 3]: (1) (2) (3) где D = Q/V, ч-1; V - объем заполнения реактора, м3; Q - объемная скорость поступающего потока м3/ч; μ - удельная скорость роста биомассы, ч-1; YX/S - стехиометрический коэффициент, г/г; X, S, P - концентрация биомассы, субстрата и продукта соответственно на выходе из реактора, г/л; Sf - концентрация субстрата в потоке, поступающем в реактор, г/л; α, β - константы. (4) где μmax - максимальная удельная скорость роста, ч-1; Pm - константа насыщения продукта, г/л; Ki - константа ингибирования субстрата, г/л; Km - константа насыщения субстрата, г/л. Стационарные условия процесса: (5) Оптимальное стационарное состояние соответствует максимальной продуктивности Qp по целевому продукту P: QP = D · P= max. Формулы для вычисления показателей процесса в оптимальных условиях приведены в [3]: (6) (7) Popt = Pm / 2, (8) . (9) (10) Условия устойчивости Для получения дальнейших решений используем понятие «устойчивость по Ляпунову» [4, 5]. Решение уравнений (1)-(3) для стационарных условий (5) при t = t0 назовем невозмущенным движением [4, 5]. Возмущенным движением назовем решение уравнений (1)-(3) при t > t0, которое возникает при начальных условиях, отличных от тех, которые были при t = t0. Приведем следующую формулировку понятия «устойчивость по Ляпунову». Стационарное состояние (или невозмущенное движение) устойчиво по Ляпунову, если для любого значения e > 0 найдется такое значение d > 0, при котором разность координат невозмущенного и возмущенного движений будет меньше e при любом значении t: t0 £ t < + ¥. В противном случае невозмущенное движение неустойчиво. Координаты возмущенного движения удобно представить в виде суммы координат невозмущенного движения и приращения, являющегося функцией времени t. Так, для системы (1)-(3) координаты возмущенного движения будут: X = Xст + d1(t); S = Sст + d2(t); P = Pст + d3(t). Качественная теория дифференциальных уравнений дает возможность оценить устойчивость по Ляпунову, не прибегая к необходимости интегрировать систему (1)-(3) [4-6]. Решение вопроса об устойчивости в математике ограничивалось в основном изучением устойчивости в первом приближении [4-6], суть которого заключалась в преобразовании системы (1)-(3) в систему относительно приращений di(t). Как показано в [6], решение вопроса об устойчивости исходной системы (1)-(3) по первому приближению не всегда дает ответ об устойчивости системы в целом. Основное заключение состоит в том, что для многих систем нелинейных дифференциальных уравнений существует ограниченная область значений, для которой имеет место устойчивость по Ляпунову. Вне этой области исходная система неустойчива. В дальнейшем системы с подобными условиями получили название систем, устойчивых «в малом» [6]. Оценка границ устойчивости «в малом» является сложной задачей, решение которой зависит от вида и характера нелинейности в уравнениях модели. Некоторые примеры такой оценки приведены в [5, 6], но, к сожалению, только для системы двух дифференциальных уравнений. Для оценки условий устойчивости используем два понятия. Первое - исследуется устойчивость в «малом», т. е. возмущения в процессе есть отклонения показателей от стационарных значений на малую величину. Это означает, что правые части уравнений (1)-(3) могут быть вычислены использованием только двух членов разложения в ряд Тейлора. Малые отклонения показателей X, S и P от стационарного состояния обозначим δ1, δ2, δ3 соответственно. Второе - устойчивость системы уравнений может быть оценена только по одному показателю решения системы (например, по δ1, или δ2, или δ3) [4]. Линеаризованная система уравнений (1)-(3) в отклонениях имеет вид (11) Формулы расчета коэффициентов уравнений (11) приведены в табл. 1, где (12) Частные производные в уравнениях (12) вычисляются для стационарных состояний. Показатели стационарных состояний вычисляются в результате решения (1)-(4) при выполнении (5). В общем случае для оценки условий устойчивости достаточно получить из системы (11) дифференциальное уравнение относительно одной переменной и использовать качественную теорию дифференциальных уравнений [4] для оценки устойчивости системы. Преобразование системы (11) к одному дифференциальному уравнению относительно δ1 дает (13) где , (14) , (15) , (16) (17) Учитывая, что в табл. 1 a1 = 0, расчетные соотношения (14)-(17) упрощаются. Таблица 1 Расчетные соотношения для вычисления коэффициентов уравнений (11) Обозначение коэффициента Расчетная формула a1 0 b1 c1 a2 a1 0 b2 c2 a3 b3 c3 Характеристическое уравнение для (13) имеет вид (18) Необходимое и достаточное условие устойчивости по Гурвицу определяется следующим образом: - матрица Гурвица: - для характеристического уравнения (18) матрица Гурвица будет: Необходимое и достаточное условие устойчивости: (19) Соотношение (19) справедливо для оценки устойчивости любого стационарного состояния системы (1)-(3) с использованием линеаризованных уравнений (11), для оптимального стационарного состояния система уравнений (11) упрощается. Обращаясь к соотношению (6) для оптимальных условий и к табл. 1, замечаем: (20) Система уравнений (11) примет вид (21) , (22) . (23) Уравнения (21) и (23) не содержат δ2, т. е. эти уравнения преобразуются в одно дифференциальное уравнение второго порядка независимо от уравнения (22): (24) или, с учетом ранее принятых обозначений, характеристическое уравнение (25) Необходимое и достаточное условие устойчивости для уравнения второго порядка Численные значения определяются значениями констант кинетического соотношения и показателями процесса в стационарном состоянии. Переходные процессы в окрестности стационарного состояния Используя систему уравнений (21)-(23) с соотношениями для коэффициентов, можно оценить характер переходных процессов в окрестности стационарного состояния и, что наиболее важно, их длительность. Указанная оценка необходима, т. к. в переходном режиме нарушается регламентное значение показателей процесса, что, вообще говоря, нежелательно. Нетрудно видеть, что система (21)-(23) может быть достаточно просто проинтегрирована, в результате чего получаем аналитическое решение для переходных процессов в окрестности стационарного состояния. Корни характеристического уравнения (25): (26) Вид решения (24) определяется видом корней (26): А. Корни r1 и r2 действительные и разные: Решение: (27) Б. Корни r1 и r2 действительные и равные: Решение: В. Корни r1 и r2 комплексные, сопряженные: . Решение: Рассмотрим решение (21)-(23) для трех вариантов типа корней (26). Вариант А: По уравнению (27) вычисляем : По уравнению (21) вычисляем δ3: (28) Соотношения (27) и (28) подставляем в (22) и получаем решение для δ2 интегрированием дифференциального уравнения первого порядка методом Бернулли: , где Таким образом, решение (21)-(23) получаем в виде . Константы интегрирования С1, С2, С3 вычисляем, используя начальные условия по методу Крамера. Начальные условия: при , где δi0 - возмущение по X (i = 1), S (i = 2), P (i = 3) соответственно. Константа с номером k вычисляется по соотношению (29) где G - определитель матрицы системы: Gk - определитель G, в котором столбец с номером k заменен на столбец h: (30) Решение (21)-(23) для вариантов Б и В получено по вышеописанной схеме. Вариант Б: (31) где Ck вычисляется по (29), где h - по (30). Вариант В: (32) где . (33) Ck вычисляется по (29), где h - по (30). Результаты численного моделирования Для численного моделирования использованы данные кинетических исследований анаэробного процесса микробиологического синтеза (табл. 2 [2, 3]). Таблица 2 Численные значения констант μ, ч-1 Pm, г/л Km, г/л Ki, г/л YX/S, г/г α, г/г β, ч-1 0,48 50 1,2 22 0,4 2,2 0,2 В основу расчета положена система линеаризованных уравнений (21)-(23), по которой получены дифференциальное уравнение (24) и характеристическое уравнение (25). Численные значения коэффициентов характеристического уравнения: P0 = 1,0 > 0; P1 = -c1 = 0,0478 и P2 = -c1a3 = 0,0268. Таким образом, Pi > 0; i = 0, 1, 2, что является необходимым и достаточным условием устойчивости. Коэффициенты системы уравнений (21)-(23) вычисляются по соотношениям в табл. 1. Данные оптимального стационарного состояния (6)-(10) представлены в табл. 3. Таблица 3 Численные значения показателей процесса в оптимальных условиях Sfopt, г/л Dopt, ч-1 Sopt, г/л Xopt, г/л Popt, г/л QP, г/(л∙ч) 23,4 0,1636 5,138 7,304 25,0 4,09 Значения коэффициентов уравнений (21)-(23): c1 = -0,0478; a2 = -0,409; b2 = -0,1636; c2 = 0,1195; a3 = 0,5599; c3 = -0,2687. Корни характеристического уравнения (26): r1,2 = -0,13437 ± 0,0933·i. Таким образом, получаем расчетную схему (вариант В) - корни комплексные, сопряженные: r1 = α0 + iβ0; r2 = α0 - iβ0, где α0 = -0,1344; β0 = 0,0933. Зависимости δ1(t), δ2(t), δ3(t) вычисляются по уравнениям (32). Константы интегрирования C1, C2 и C3 находятся по правилу Крамера по определителю G (33) с использованием (29) и (31): где r3 = b2 = -0,1636; δ10, δ20, δ30 - начальные значения δ1, δ2, δ3 при t = 0; ; . Результаты вычислений по (32) показаны на рис. 1-3: рис. 1 - зависимости δ10, δ20, δ30 при δ10 = ± 0,01; δ20 = δ30 = 0,0; рис. 2 - зависимости δ10, δ20, δ30 при δ20 = 0,01; δ10 = δ30 = 0,0; рис. 3 - зависимости δ10, δ20, δ30 при δ30 = ± 0,01; δ20 = δ10 = 0,0. Рис. 1. Переходный процесс: при δ10 = 0,01; δ20 = δ30 = 0,0 - линии δ1, δ2, δ3; при δ10 = -0,01; δ20 = δ30 = 0,0 - линии δ1*, δ2*, δ3* Рис. 2. Переходный процесс: при = 0,01; = = 0,0 - линия δ2; при = -0,01; = = 0,0 - линия δ2* Рис. 3. Переходный процесс: при = 0,01; = = 0,0 - линии δ1, δ2, δ3; при = -0,01; = =0,0 - линии δ1*, δ2*, δ3* На рис. 4 показаны фазовые траектории в пространстве переменных δ1, δ2, δ3. Рис. 4. Фазовые траектории в переходном процессе: 1 - = 0,01; = = 0,0; 1* - = -0,01; = = 0,0; 2 - = 0,01; = = 0,0; 2* - = -0,01; = = 0,0; 3 - = 0,01; = = 0,0; 3* - = -0,01; = = 0,0 На рис. 4 стрелками указано направление движения по фазовым траекториям. Для полученного устойчивого состояния фазовые траектории сходятся к точке, называемой «устойчивый узел». Заключение Разработанная методология оценки устойчивости стационарных состояний микробиологических процессов, реализуемых по непрерывной схеме, носит достаточно общий характер. Указанная оценка является особенно важной для процессов с нелинейными кинетическими соотношениями, для которых, как правило, затруднительно получать аналитические решения и почти невозможно прогнозировать устойчивость состояний экспериментально. Необходимо заметить, что рассмотренный числовой пример относится к процессу, протекающему в оптимальных условиях. В то же время использование методологии в полном объеме требует наличия кинетических данных в широком диапазоне изменения показателей процесса, при этом возникают повышенные требования к точности этих данных с тем, чтобы избежать неправильных выводов.
Список литературы

1. Гордеева Ю. Л. Моделирование непрерывного процесса биотехнологического получения молочной кислоты / Ю. Л. Гордеева, Ю. А. Ивашкин, Л. С. Гордеев // Теоретические основы химической технологии. 2012. Т. 46, № 3. C. 324-328.

2. Kumar G. P. Periodic operation of a bioreactor with input multiplicities / G. P. Kumar, J. V. K. Sastry, M. Chidambaram // The Canadian Journal of Chemical Engineering. 1993. Vol. 71, no. 5. P. 766-770.

3. Гордеева Ю. Л. Моделирование периодического процесса микробиологического синтеза с нелинейной кинетикой роста микроорганизмов / Ю. Л. Гордеева, Ю. А. Ивашкин, Л. С. Гордеев, М. Б. Глебов. М.: РХТУ им. Д. И. Менделеева, 2011. 100 с.

4. Демидович Б. П. Лекции по математической теории устойчивости / Б. П. Демидович. М.: Изд-во МГУ, 1998. 480 с.

5. Степанов В. В. Курс дифференциальных уравнений / В. В. Степанов. М.: ЛКИ, 2008. 472 с.

6. Фельдбаум А. А. Теоретические основы связи и управления / А. А. Фельдбаум, А. Д. Дудыкин, А. П. Мановцев, Н. Н. Миролюбов. М.: Гос. изд-во физ.-мат. лит., 1963. 932 с.