Введение Процесс непрерывного микробиологического синтеза технологически реализуется в стационарных условиях, т. е. показатели процесса не изменяются во времени (концентрация биомассы, субстрата, продукта и т. п.). При реализации синтеза в процессе всегда имеют место возмущения (отклонения показателей процесса от стационарного значения). Если эти отклонения не приводят к нарушению технологического процесса и с течением времени значения показателей возвращаются к первоначальным, стационарное состояние полагается устойчивым. В противном случае для реализации процесса требуется внешнее управление. Так как малых возмущений практически избежать не удается, то, если при малых возмущениях наблюдается нарушение технологического процесса, приводящее к невозможности обеспечить его функциональное назначение, процесс считается неустойчивым в малом. Математический анализ такого рода устойчивости или неустойчивости в малом получил название «устойчивость первого приближения», оценка которой возможна по условиям Гурвица. Подробный математический анализ в доступной форме приведен в работах [1, 2]. Таким образом, настоящее сообщение будет базироваться на анализе устойчивости первого приближения. Ниже этот подход будет описан подробнее. Объект исследования Анализ устойчивости стационарных состояний технологического процесса возможен с использованием математической модели, воспроизводящей с достаточной точностью поведение показателей процесса. Процессы микробиологического синтеза распространены достаточно широко. Однако здесь речь идет о промышленно значимых процессах получения целевых продуктов [3], в частности о процессах получения пищевых кислот [4]. Рассматриваемые процессы отвечают следующим условиям: в процессе синтеза производится биомасса X, получается целевой продукт P, расходуется основной субстрат S. Сущность процесса определяется кинетикой роста биомассы, основным показателем которой является удельная скорость роста m, т. к. биомасса является продуцентом образования продукта. Соотношения для удельной скорости роста сформированы на основе неструктурированного подхода. Математически удельная скорость роста для рассматриваемого объекта является в общем случае нелинейной функцией X, S, P. Кроме того, в процессе ферментации удельная скорость роста в той или иной степени может быть ингибирована каждым из показателей X, S, P. В этом случае возможность ингибировании включается в математическое соотношение для m. Обращаясь к обзорам [5, 6], отметим, что для отдельных штаммов микроорганизмов при ферментации образуются, кроме целевого продукта, побочные, иногда представляющие практическую ценность. При этом таких продуктов образуется, как правило, незначительное количество [7]. С другой стороны, т. к. сырье является наиболее расходной статьей процесса, имеется тенденция [5, 6] использовать воспроизводимое сырье. Это отмечено в работах [8-10]. Особенность использования воспроизводимого сырья проявляется в том, что в этом сырье могут быть компоненты, которые в процессе организации технологии дают дополнительное количество основного субстрата S. Так, например, пшеничная мука (whole-wheat flour, wheat flour) в качестве сырья для ферментации содержит мальтозу, которая в процессе генерирует дополнительно основной субстрат - глюкозу, что повышает уровень эффективности использования сырьевых ресурсов. Схема такого процесса приведена на рисунке [10]. Схема процесса ферментации с получением дополнительного количества основного субстрата Учет этой составляющей, а также составляющей побочного продукта позволяет записать обобщенную математическую модель рассматриваемого объекта в виде системы дифференциальных уравнений, представляющих материальный баланс по четырем компонентам - S, X, P, M, где M есть концентрация побочного продукта или концентрация компонента, использование которого дает дополнительное количество основного субстрата. Система уравнений записывается для нестационарных условий процесса, организованного по непрерывной схеме: , (1) , (2) , (3) , (4) где F1, F2, F3, F4 - нелинейные функции показателей процесса (S, X, P, M); S, X, P - концентрации субстрата, биомассы, продукта соответственно; M - концентрация побочного продукта или компонента сырья, деградирующего в основной субстрат (в последнем случае M в соотношениях (2) и (3) отсутствует, соотношение (4) значительно упрощается (будет показано ниже)). Таким образом, объектом нашего исследования является процесс микробиологического синтеза, оформленный по непрерывному способу. Цель процесса - получение продукта в результате жизнедеятельности микроорганизмов. Математическая модель объекта представляет собой четыре уравнения материального баланса по компонентам: основному субстрату, биомассе, продукту, побочному продукту или компоненту сырья, восполняющему количество основного субстрата. Если в сырье отсутствует побочный продукт или компонент, дополнительно воспроизводящий основной субстрат, уравнения материального баланса записываются только по трем компонентам: S, X, P. С учетом вышеизложенного система уравнений математической модели в общем виде представлена уравнениями (1)-(4). Наиболее полно этому описанию объекта соответствует процесс получения молочной кислоты по непрерывному способу из пшеничной муки (wheat flour) [8]. Обращаясь к понятию «устойчивость первого приближения», запишем математические соотношения, определяющие это понятие. Система уравнений математической модели (1)-(4) преобразуется в систему линейных дифференциальных уравнений первого порядка с постоянными коэффициентами, записанных для малых отклонений от стационарных значений. Значения переменных представляются в следующем виде: S = Sст + d1; X = Xст + d2; P = Pст + d3; M = Mст + d4, где Sст, Xст, Pст, Mст - значения показателей для конкретного стационарного состояния; d1, d2, d3, d4 - малые отклонения от стационарных состояний. Нелинейные функции Fi в правой части системы (1)-(4) представлены рядом Тейлора по малым отклонениям с сохранением только двух первых членов ряда (в этом смысл первого приближения). Так, для функции F1 (1) имеем где частные производные вычисляются для конкретного стационарного состояния, т. е. они являются константами. В наших условиях F1ст = 0. В итоге, обозначив: ; ; ; и использовав эти обозначения для остальных уравнений системы (1)-(4), получим систему уравнений первого приближения: , (5) , (6) , (7) . (8) В соответствии с [1, 2], если какое-либо решение системы устойчиво, устойчива и вся система. Таким образом, система (5)-(8) служит основой для получения заключения об устойчивости первого приближения для любого рассматриваемого стационарного состояния. Цель настоящей работы - получить алгоритм, выполнение которого позволит оценить условия устойчивости первого приближения для объекта микробиологического синтеза, описанного выше. Описание алгоритма Алгоритм включает пять последовательно выполняемых этапов, что обеспечивает решение поставленной задачи, т. е. получение возможности оценить устойчивость первого приближения принятого стационарного состояния. Этап 1 является определяющим в формировании условий для математического обеспечения решения задачи. Записываются система уравнений математической модели (1)-(4) и соотношение для удельной скорости роста m с учетом возможного ингибирования по разным показателям процесса (уравнения получают с использованием источников [5, 6]). Записываются численные значения констант. Определяются пределы изменения входных показателей процесса ; (по возможностям сырьевого обеспечения). Определяются значения величины протока D : Dmin и Dmax. Значение Dmax вычисляется по условию Dmax < Dпред, где Dпред есть значение протока, при котором исходный субстрат концентрацией S0 полностью вымывается из аппарата. Так как по материальному балансу по X имеем уравнение для непрерывного процесса в виде D = m, значение Dпред будет равно m, в котором S = S0; X = 0; P = 0; M = 0. В результате имеется система уравнений математической модели, ограничения по входным показателям процесса, численные значения констант. На этапе 2 определяются стационарные состояния (одно или несколько), для которых необходимо оценить условия устойчивости. Стационарные состояния определяются заданием численных значений входных показателей S0, M0 и D, естественно, в заданных пределах. Для принятых стационарных состояний вычисляются показатели процесса Sст, Xст, Pст, Mст по решению системы линейных алгебраических уравнений вида На этапе 3 формируется система линейных дифференциальных уравнений первого порядка с постоянными коэффициентами на основе уравнений математической модели - система первого приближения. Оцениваются значения коэффициентов системы для каждого принятого стационарного состояния. Система уравнений первого приближения для рассматриваемого объекта имеет вид (5)-(8). Коэффициенты уравнений вычисляются разложением функций Fi в ряд Тейлора с использованием первых двух членов ряда, если иметь в виду, что в стационарном состоянии F1ст = F2ст = F3ст = F4ст = 0. Таким образом, коэффициенты ai, bi, ci, di вычисляются по следующим формулам: ; ; ; . Значения частных производных получаем с использованием значений Sст, Xст, Pст, Mст для каждого принятого стационарного состояния. Эти значения получены на этапе 2. В результате имеем систему линейных дифференциальных уравнений первого порядка с известными постоянными коэффициентами. Реализация этапа 4 позволяет вычислить элементы матрицы Гурвица, которая является основой оценки условий устойчивости. Матрица Гурвица для рассматриваемого объекта имеет вид . (9) Приведем два варианта вычисления элементов матрицы Гурвица P1, P2, P3, P4. Первый вариант реализуется путем преобразования системы четырех дифференциальных уравнений первого порядка с постоянными коэффициентами в одно дифференциальное уравнение четвертого порядка, т. е. система (6)-(9) преобразуется в уравнение . (10) Уравнение (10) записано относительно переменной d1. Отметим, что преобразование возможно по любой из переменных. Элементы матрицы Гурвица Pi вычисляются по формулам, приведенным в табл. 1. Соотношения табл. 1 получены в последовательности преобразования (5)-(8) в уравнение четвертого порядка (последовательность преобразования не приводится). Таблица 1 Соотношения для расчета коэффициентов P0, P1, P2, P3, P4 Соотношение ; ; ; M1 = b1a2 + c1a3 + d1a4; M2 = b1b2 + c1b3 + d1b4; M3 = b1c2 + c1c3 + d1c4; M4 = b1d2 + c1d3 + d1d4 N1 = M2a2 + M3a3 + M4a4; N2 = M2b2 + M3b3 + M4b4; N3 = M2c2 + M3c3 + M4c4; N4 = M2d2 + M3d3 + M4d4 L1 = N2a2 + N3a3 + N4a4; L2 = N2b2 + N3b3 + N4b4; L3 = N2c2 + N3c3 + N4c4; L4 = N2d2 + N3d3 + N4d4 K1 = M2 + b1a1; K2 = M1b1 - M2a1; K3 = M3b1 - M2c1; K4 = M4b1 - M2d1 R1 = M1b1 + N2; R2 = N1b1 - N2a1; R3 = N3b1 - N2c1; R4 = N4b1 - N2d1 S1 = N1b1 + L2; S2 = L1b1 - L2a1; S3 = L3b1 - L2c1; S4 = L4b1 - L2d1 G1 = R3b1 + K3a1b1; G2 = K3R1 - K1R3; G3 = K3R2 - K2R3; G4 = K3R4 - K4R3 F1 = S3b1 + M1K3b1; F2 = K3S1 - K1S3; F3 = K3S2 - K2S3; F4 = K3S4 - K4S3 P0 = K3G4b1; P1 = -(K3F4b1 + K3G4a1b1); P2 = -(G4F1 - G1F4); P3 = -(G4F2 - G2F4); P4 = -(G4F3 - G3F4) P1 = P1 / P0; P2 = P2 / P0; P3 = P3 / P0; P4 = P4 / P0; P0 = 1 Второй вариант реализуется следующим образом. Система уравнений (5)-(8) записывается в векторной форме: , где A - матрица коэффициентов уравнений системы (5)-(8): Используя собственные числа l и матрицу A, записываем характеристическое уравнение в виде Раскрывая определитель, получаем уравнение вида На этом завершается этап 4. Этап 5 является завершающим в описании алгоритма, который дает ответ на вопрос об устойчивости или неустойчивости рассматриваемых стационарных состояний. Заполняется матрица Гурвица (10) и вычисляются следующие определители: ; ; ; . (11) Необходимое и достаточное условие устойчивости . Если какое-либо из Ai < 0, стационарное состояние неустойчиво. Использование алгоритма приведено в следующем разделе. Анализ устойчивости стационарных состояний процесса получения молочной кислоты из пшеничной муки Исходные данные для моделирования в достаточно полном объеме представлены в работе [10]. Уравнения математической модели для нестационарных условий процесса по непрерывной технологии: , (12) , (13) , (14) . (15) Уравнение для удельной скорости роста записано в виде . (16) Учитывая, что KS << S (в результате оценок по эксперименту) соотношение (16) используется в виде . В соотношениях (12)-(16) m, mmax - удельная и максимальная скорость роста, ч-1; S, X, P, M - концентрация компонентов, г/л; S0, M0 - концентрация компонентов во входном потоке, г/л; D - величина протока, ч-1; YS, YP - стехиометрические коэффициенты, безразм.; Km, Pmax, n - константы. Численные значения констант: mmax = 0,28 ч-1; Pmax = 98,6 г/л; KS = 0,5 г/л; YS = 0,053; YP = 0,52; KM = 0,035; n = 3,0. Область задания S0, M0 и D: S0 = 50-120 г/л; M0 = 20-80 г/л; Dпред = 0,28 ч-1 (определено по условию материального баланса по X). К рассмотрению принято два стационарных состояния: первое - произвольное в допустимых пределах по S0, M0 и D; второе - для оптимальных условий, которые определены для максимальной продуктивности Qp = DP = maxQP по величине протока Dopt. На этом этап 1 завершен. Стационарное состояние 1. Значения показателей получены решением уравнений (12)-(15) по условию Условия на входе: S0 = 100 г/л; M0 = 50 г/л; D = 0,15 ч-1. Показатели стационарного состояния: Sст = 86,8731 г/л; Xст = 1,1971 г/л; Pст = 18,5205 г/л; Mст = 40,5405 г/л; Qp = 2,7781 г/(л · ч). Стационарное состояние 2. Получено по решению задачи оптимизации для тех же исходных значений S0 и M0. Оптимальное значение Dopt = 0,1181 ч-1. Значения показателей в стационарных условиях: Sст = 81,3675 г/л; Xст = 1,5932 г/л; Pст = 24,65 г/л; Mст = 38,5714 г/л; maxQP = 2,9118 г/( л · ч). Завершается этап 2. Далее вычисляются коэффициенты системы уравнений первого приближения (5)-(8) ai, bi, ci, di. Для вычисления коэффициентов составляем табл. 2. Стационарное состояние 1: a1 = -0,15; b1 = -2,8302; c1 = 0,1269; d1 = 0,035; a2 = 0,0; b2 = 0,0; c2 = -0,006727; d2 = 0,0; a3 = 0,0; b3 = 2,3207; c3 = -0,25407; d3 = 0,0; a4 = 0,0; b4 = 0,0; c4 = 0,0; d4 = -0,185. Стационарное состояние 2: a1 = -0,118125; b1 = -2,2288; c1 = 0,14405; d1 = 0,035; a2 = 0,0; b2 = 0,0; c2 = -0,0076347; d2 = 0,0; a3 = 0,0; b3 = 1,8276; c3 = -0,2369; d3 = 0,0; a4 = 0,0; b4 = 0,0; c4 = 0,0; d4 = -0,1531. Таблица 2 Соотношения для вычисления коэффициентов ai, bi, ci, di a1 = -D d1 = KM a2 = 0,0 b2 = 0,0 d2 = 0,0 a3 = 0,0 d3 = 0,0 a4 = 0,0 b4 = 0,0 c4 = 0,0 Завершен этап 3. Для вычисления элементов матрицы Гурвица использует второй вариант расчета. Для системы уравнений первого приближения вида (5)-(8) матрица A с учетом численных значений констант ai, bi, ci, di имеет вид . Характеристическое уравнение: или , где ; ; ; Таким образом, для первого стационарного состояния имеем: a1 = -0,15; c2 = -0,006727; b3 = 2,3207; c3 = -0,25407; d4 = -0,185; P0 =1,0; P1 = 0,58907; P2 = 0,12807; P3 = 0,01228; P4 = 0,0004332. Для второго: a1 = -0,118125; c2 = -0,0076347; b3 = 1,8276; c3 = -0,2369; d4 = -0,1531; P0 = 1,0; P1 = 0,508125; P2 = 0,096291; P3 = 0,008069; P4 = 0,000252. Завершен этап 4. Используя Pi, записываем матрицу Гурвица: По матрице Гурвица вычисляются определители (11): - для первого стационарного состояния: A1 = 0,58907; A2 = 0,0634; A3 = 0,000628; A4 = 2,72×10-7; - для второго стационарного состояния: A1 = 0,50125; A2 = 0,04086; A3 = 2,646×10-4; A4 = 6,6684×10-8. Поскольку определители для первого и второго состояния положительны, по необходимому и достаточному условию [1, 2] оба стационарных состояния устойчивы по первому приближению. Заключение В заключение отметим, что разработанный алгоритм может использоваться вне зависимости от сложности математического представления кинетических соотношений и числа уравнений математической модели. Единственным ограничением является то, что анализ проводится для условий первого приближения, т. е. в малой окрестности стационарного состояния, что в целом не позволяет делать заключение об устойчивости процесса при наличии больших возмущающих воздействий.