MATHEMATICAL MODELING OF STABILITY OF BIOTECHNOLOGICAL PROCESS IN ITS OPTIMUM CONDITION
Abstract and keywords
Abstract (English):
Practical realization of microbiological synthesis in continuous operating conditions is connected with the stability of operation at steady state. Prediction of the parameters of the process, provided sustainable operation, is carried out using mathematical modeling, methodology of which includes assessment of steady states parameters and changes of the character for possible disturbances in the process. The methodology has been developed to estimate the stability of the stationary states of the process of anaerobic microbial synthesis in optimal operating conditions. These methodology based on the use of the linearized system of equations of unsteady state recorded in deviations. The relations are obtained for calculating the coefficients of the linearized equation. The system of equations is reduced to a single differential equation of the third order. Hurwitz matrix was composed and the necessary and sufficient conditions for stability were recorded. The system of equations has been simplified and possible variants of the solutions for different types of roots of the characteristic equation have been received. The analytical solutions for unsteady state have been received. The numerical example of this methodology used known kinetic relations for different processes of microbiological synthesis is presented in this report. The productivity of the desired product has been used as an optimality criterion. The results of modeling showed that the optimal condition is stable for these process parameters. The methodological approach can be recommended for the analysis of other kinetic schemes of the microbiological processes.

Keywords:
biotechnology, modeling, stability
Text
Введение Одно из направлений совершенствования процессов микробиологического синтеза связано с созданием оптимальных условий, определяющих максимальный выход целевых продуктов. Решение этой задачи не может быть реализовано только на базе экспериментальных исследований. Моделирование оптимальных условий осуществляется с использованием кинетических соотношений, вид которых все более усложняется [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 стрелками указано направление движения по фазовым траекториям. Для полученного устойчивого состояния фазовые траектории сходятся к точке, называемой «устойчивый узел». Заключение Разработанная методология оценки устойчивости стационарных состояний микробиологических процессов, реализуемых по непрерывной схеме, носит достаточно общий характер. Указанная оценка является особенно важной для процессов с нелинейными кинетическими соотношениями, для которых, как правило, затруднительно получать аналитические решения и почти невозможно прогнозировать устойчивость состояний экспериментально. Необходимо заметить, что рассмотренный числовой пример относится к процессу, протекающему в оптимальных условиях. В то же время использование методологии в полном объеме требует наличия кинетических данных в широком диапазоне изменения показателей процесса, при этом возникают повышенные требования к точности этих данных с тем, чтобы избежать неправильных выводов.
References

1. Gordeeva Yu. L. Modelirovanie nepreryvnogo processa biotehnologicheskogo polucheniya molochnoy kisloty / Yu. L. Gordeeva, Yu. A. Ivashkin, L. S. Gordeev // Teoreticheskie osnovy himicheskoy tehnologii. 2012. T. 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. Gordeeva Yu. L. Modelirovanie periodicheskogo processa mikrobiologicheskogo sinteza s nelineynoy kinetikoy rosta mikroorganizmov / Yu. L. Gordeeva, Yu. A. Ivashkin, L. S. Gordeev, M. B. Glebov. M.: RHTU im. D. I. Mendeleeva, 2011. 100 s.

4. Demidovich B. P. Lekcii po matematicheskoy teorii ustoychivosti / B. P. Demidovich. M.: Izd-vo MGU, 1998. 480 s.

5. Stepanov V. V. Kurs differencial'nyh uravneniy / V. V. Stepanov. M.: LKI, 2008. 472 s.

6. Fel'dbaum A. A. Teoreticheskie osnovy svyazi i upravleniya / A. A. Fel'dbaum, A. D. Dudykin, A. P. Manovcev, N. N. Mirolyubov. M.: Gos. izd-vo fiz.-mat. lit., 1963. 932 s.


Login or Create
* Forgot password?