Оценивание скоростей реакций и неизмеряемых внешних возмущений биотехнологического процесса Текст научной статьи по специальности «Промышленные биотехнологии»

Вычислительные технологии

Том 14, № 2, 2009

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

А.Ю. Торглшов

Институт автоматики и процессов управления ДВО РАН, Владивосток, Россия

e-mail: [email protected]

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

Ключевые слова: биотехнологический процесс, неизмеряемые внешние возмущения, оценивание параметров, регуляризирующий оператор.

Введение

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

— расширенных фильтров Калмана [1];

— искусственных нейронных сетей (ИНС) [2];

— регуляризационно-фильтрационных наблюдателей [3, 4].

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 08-08-00004-а; 06-08-96014-р-восток-а) и Совета по грантам Президента РФ для поддержки молодых российских ученых (грант МК-2034.2008.8).

© ИВТ СО РАН, 2009.

от вышеперечисленных недостатков и поэтому наиболее применим в условиях неопределенности, т.е. когда:

— отсутствует точная информация о кинетических закономерностях процесса;

— присутствуют время задержки измерений концентраций реагирующих веществ и погрешность измерений технологических параметров;

— изменения состава субстратов являются неизвестными функциями во времени.

Алгоритм оценивания строится на основе инверсии нелинейных динамических систем и последующей фильтрации. Параметры регуляризации при этом находятся с помощью, например, метода ¿-кривой [5]. Однако регуляризационно-фильтрационный подход, предлагаемый авторами работ [3, 4], еще недостаточно проработан. В частности, не показана возможность использования нескольких параметров регуляризации в случае оценивания одновременно нескольких переменных. Также все расчеты проводились при довольно низком уровне помехи и неясно, как отразится увеличение шума на качестве оценивания. В данной работе некоторые перечисленные проблемы решены на примере двух моделей биореактора и предложен более универсальный алгоритм оценивания на основе регуляризационно-фильтрационного подхода.

1. Математическое моделирование биотехнологического процесса

В настоящей работе рассматривается БТП, протекающий в биореакторе, т. е. аппарате, в котором одни вещества, называемые субстратами, превращаются в другие, называемые продуктами, при участии живых организмов (бактерии, дрожжи, клеточные культуры, макроорганизмы) или ферментативных систем, выделенных из них. Входные и выходные переменные БТП приведены на блок-схеме (рис. 1). Важной особенностью биотехнологических процессов при участии живых организмов является то, что эти организмы работают аналогично катализатору превращения субстрата (питательные вещества, минеральные соли) в продукт. Если скорости процессов зависят только от концентрации одного субстрата, в то время как остальные присутствуют в избытке, то такой субстрат принято называть лимитирующим. Лимитирующим субстратом обычно становится главный источник питания — глюкоза, сахароза или какой-либо другой источник углеводов. Продуктом называют вещество (или группу веществ), которое выделяется организмами в окружающую среду в процессе жизнедеятельности. Продукт может быть как ценным (белки, лекарственные препараты, различные биологически активные вещества), так и не имеющим практической ценности (углекислый газ и другие продукты жизнедеятельности). Кроме того, помимо реакций образования продук-

ш

п

и [£>]

> Биореактор

Рис. 1. Переменные биотехнологического процесса: и — управляющее воздействие; ! — внешнее неизмеряемое ограниченное возмущение; п — шум измерений; у — измеряемый вектор выходных переменных

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

Биореактор — достаточно сложный объект наблюдения и управления. Дифференциальные уравнения, описывающие динамику биотехнологических процессов, сильно различаются в зависимости от особенностей каждого процесса. Даже в простейшем случае эти уравнения являются нелинейными. Кинетические закономерности биохимических реакций выражаются в зависимости удельных (на единицу биомассы) скоростей соответствующих реакций от концентраций. Классическим примером такой зависимости служит уравнение Моно

^тахБ

V

Ks + S'

(1)

где ^ — удельная скорость роста биомассы; ^тах — максимальная удельная скорость; Кб — субстратспецифичная константа; Б — концентрация лимитирующего субстрата. Более сложной моделью кинетики является модель Холдейна

V

VmaxS

Ks + S + S2/Kt

(2)

где К — константа ингибирования избытком субстрата. В данной модели учтено инги-бирование реакции субстратом при его большой концентрации. В случае если наблюдается ингибирование роста биомассы продуктом Р, то имеет место одна из зависимостей:

V

или

V

í VmaxS

\Ks + S,

VmaxS

K

P

Kp + Pt

\Ks + S

1 -

P

либо их комбинация

V

VmaxS

Ksi + S + S 2/Кг

1-

Pm

P

Pm

K

P

Kp + Pt

(3)

(4)

(5)

Здесь Ртах — концентрация продукта, при которой рост популяции прекращается; Кр — константа ингибирования; т, п — положительные числа. Очевидно, что все приведенные уравнения (1)-(5) — нелинейные относительно концентрации субстрата. Скорости роста биомассы М и биосинтеза V, характеризующие процесс в целом, являются произведением соответствующих удельных скоростей и концентрации биомассы X:

M = vX, V = vX.

(6) (7)

n

m

m

n

2. Режимы функционирования БТП и постановка задачи

В математическую модель БТП входит система уравнений материальных балансов из [6]:

dX/dt = vX - DX,

dP/dt = vX - DP, (8)

dS/dt = DSin - DS - YxvX - YpvX,

где

^ _ (кi + S + s2/KJ - P^x) (K + P) ' (9)

v = K^ • (io)

K S2 + S

Здесь D — скорость разбавления (может являться управляющим воздействием); Sin — концентрация подаваемого субстрата; Yx — коэффициент выхода биомассы; Yp — коэффициент выхода продукта; Ksi, Ks2 — субстратспецифичные константы. Приведенная система уравнений материального баланса (8) справедлива для широкого круга биотехнологических процессов (с общими скоростями реакций (6), (7)), в которых используется один вид микроорганизмов, культивируемый непрерывно на одном субстрате с целью получения единственного продукта-метаболита в биореакторе, режим работы которого близок к режиму идеального перемешивания. Анализ статических режимов системы, описанной уравнениями (8)-(10) путем имитационного моделирования при различных начальных условиях и скоростях разбавления D, показал, что система стремится к одному из двух типов стационарных режимов: при избыточном разбавлении и при недостаточной концентрации субстрата. При избыточном разбавлении субстрат усваивается частично, т. е. его концентрация внутри биореактора, а следовательно, и в выходном потоке далека от нулевой. При такой концентрации субстрата синтез продукта идет активно. При недостаточной концентрации субстрата он усваивается практически полностью, т. е. его концентрация как внутри биореактора, так и в выходном потоке близка к нулевой. Однако при понижении концентрации субстрата скорость роста биомассы падает в значительно меньшей степени, чем скорость биосинтеза. Следовательно, при таком режиме практически весь субстрат переходит в биомассу, а количество выделяющегося продукта мало. Другая модель БТП была предложена ранее [7] для описания биокинетики разложения перекиси водорода морским организмом. Система уравнений материального баланса в этом случае выглядит следующим образом:

dP/dt _ vFR - vSobF,

(11)

dS/dt _ -2 dP/dt,

где P — концентрация растворенного кислорода (продукт); S — концентрация перекиси (субстрат); v — удельная скорость образования кислорода; F — коэффициент, учитывающий число, форму и массу морского организма; R — коэффициент, характеризующий интенсивность перемешивания; vSob — собственная скорость поглощения кислорода морским организмом. Для удельной скорости образования кислорода v принимается кинетическая модель Холдейна, учитывающая ингибирование избытком субстрата:

_ _vmaxS__(12)

v _ KS + S + S2/Ki, ( )

где vmax — максимальная удельная скорость образования кислорода; Ks — субстрат-специфичная константа; Ki — константа ингибирования избытком субстрата. Скорости биосинтеза и роста биомассы являются параметрами системы и не могут быть измерены напрямую.

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

3. Оценивание неизвестных переменных БТП

Очевидно, что поставленная задача является обратной, т. е. когда по известным выходным сигналам необходимо определить входные сигналы. Выходные сигналы для системы (8) — X, Б и Р, а для системы (11) — Б и Р. Строго говоря, ^ и V являются здесь параметрами системы, однако мы будем рассматривать их как неизвестные входные функции времени. У такого подхода есть огромное преимущество — отсутствие необходимости знания зависимостей (9), (10), (12) для оценивания соответствующих сигналов. При этом требуется знание только системы уравнений материального баланса (8) и (11). Пусть Ш = [ш1'ш2---'шп1\т — приведенный вектор входа, включающий в себя параметры объекта (которые подлежат оцениванию) и действующие входные возмущения; У = [у1у2---уп]т — вектор выхода, тогда зависимость выхода от входа можно представить в виде операторного уравнения

У = АШ. (13)

С учетом воздействия шумов в преобразованном по Лапласу виде уравнение (13) представляется следующим выражением:

У(в) = Ор(з)Х (¿0 + N (в), (14)

где У (в) — измеренный (зашумленный) выход; Ор — передаточная функция объекта; N (в) — шум с неизвестными статистическими характеристиками; в — оператор Лапласа. При отсутствии ошибок измерений выходных сигналов можно восстановить входные сигналы, однако на практике всегда имеется некоторая шумовая составляющая регистрируемых сигналов, что делает решение неустойчивым, поэтому применяется ре-гуляризирующий оператор (фильтр) Gf (в). Регуляризированное решение будет иметь следующий вид:

Ж(в) = Gf (в^(в)У(в), (15)

где Ш(в) — оцениваемое значение входа; G~1 — инверсный оператор.

Распространим вышеизложенные выражения (14) и (15) на нелинейные системы нашего случая (8) и (11). Обращенная система уравнений (8) будет выглядеть следующим образом:

Д = ^ + А (16)

х

,> = ^ + о1, (17)

XX

Бп = — {¿Б/¿г + Ух йх/¿г + Ур ¿Р/¿г) + § + р + X ,

(18)

где IV = [¡1Б ¿У1; оценки производных находятся путем непосредственного дифференцирования отфильтрованных выходных переменных, т.е. ¿{¿X/¿г} = эХ = sXGf, ¿{¿Р/¿г} = эР = sPGf и ¿{¿¿/¿г} = эБ = sSGf. Для того чтобы оценить три входных сигнала, необходимо знать три выходных. В данном случае могут быть измерены сигналы Р, и X.

Обращенная система уравнений (11) примет вид

¿Р/¿г + и8оЪ Б = РЕ '

/¿г = -мР/¿г,

(19)

(20)

для которой достаточно измерение только одной переменной или Р.

Авторы работ [3, 4] предлагали использовать непрерывный регуляризирующий фильтр с передаточной функцией:

Gf (Э) 7 I \пг,

(э + а) '

где nf — порядок фильтра. В настоящей работе предлагается использовать его дискретный аналог, так как управление осуществляется с помощью компьютера в дискретные моменты времени. Для экспоненциального фильтра первого порядка (nf = 1) такой аналог мы легко найдем, применяя ^-преобразование:

Gf (г) = z{L-1 ^ (э) \п/=1)}

ах

г — е

-аТ'

(21)

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

Рис. 2. Ь-кривая для поиска оптимального параметра регуляризации применительно к системе (11); у = [Б; Р]Т

метода L-кривой. Описание этого метода можно найти в работе [5], а также в других работах ее автора. Согласно данному методу оптимальное значение параметра соответ--кривой с максимальной кривизной Ksoi, значение которой определялось

из условия

Ksol = sup {Qk} ,

где Ki =

yi - yi-1 yi+1 - yi

xi xi—1 Xi+1 xi

G Qr , Xi и y i — координаты точек L-кривой (абсцисса

и ордината соответственно); г — номер точки. Каждая точка ¿-кривой соответствует определенному значению параметра регуляризации а. Например, в случае системы (11) оцениванию подвергалась величина V. На рис. 2 представлена ¿-кривая, полученная путем моделирования процесса оценивания при варьировании параметра регуляризации а.

Заключение

В случае системы (8) оцениванию подвергались сразу три переменных: v и Sin на основе измерений S, P и X. Для сопоставимости полученных результатов с результатами работы [3] был установлен тот же уровень дисперсии белого шума, а скорость разбавления D изменялась по трапецеидальному закону. Входное возмущение Sin изменялось около значения 95 г/л ступенчато на величину 20%. Параметра регуляризации было два: один — для фильтра удельных скоростей / и v, другой — для фильтра Sin (обозначены соответственно как а\ и ^2).

Исходя из структуры инверсной системы (16)-(18), возможно независимое определение оптимальных значений а\ и а2. Численным методом найдено a°pt = 0.016, которое в дальнейшем фиксировалось при определении a2pt (рис. 3 и 4). При оптимальных параметрах регуляризации удельные скорости / и v оценивались с минимальной ошибкой (рис. 5). Имеют место запаздывание по фазе примерно на 7 ч и большие отклонения

зппппп

Рис. 3. L-кривая для поиска оптимального параметра регуляризации ai; y = [X; P; S]T;

w = [м; v ]T

Рис. 4. L-кривая для поиска оптимального параметра регуляризации a2

0,22 0,2 0,18

М,я

-1

0,16 0,14 0,12 0,1

А

/V * \ * А

/ / \ ^ ^ V1 к // г

1 1 \ \ /V*

V

0

50

100 Время, ч

150

0,12 0,1

V, ч

0,08 1

0,06 0,04 0,02 О

г\ ¡ А I М К 1 -■

У V 1. 11 № / V//

1 1 * *

(

11

200

50

100

Время, ч

150

200

Рис. 5. Значения ^ и V при оптимальных параметрах регуляризации: 1 — реальные; 2 — оцененные

20

18

16

14

12

(-

сС 10

К о

о

6

4

2

0

2 л

/ л

Д1

/ \

£ J ' \ 1 1 / .-» \

г к 1. \ > \ * Г /VI. \

/ \

/ // . -У-/ ^л?

50

100 Время, ч

150

200

160 140 120 100

^

^ 80 со

60 40 20 0

1

* ¿лА

V Л

V» V 1*1

11 7ч

50

100

Время, ч

150

200

Рис. 6. Значения переменных X, Р и 5 при оптимальных параметрах регуляризации: 1 реальные; 2 — оцененные

Рис. 7. Значения скоростей роста биомассы М и биосинтеза V: 1 — реальные; 2 — оцененные

Рис. 8. Значения внешнего возмущения Sin: 1 — реальные; 2 — оцененные

в первые 20 ч (шаг моделирования — 0,1 ч). При этом концентрации X, Р, £ уточняются (восстанавливаются из зашумленных данных) с большей точностью, особенно для £ (рис. 6). Скорости роста биомассы М и биосинтеза V, характеризующие процесс в целом, показаны на рис. 7. Очевидно, что ошибки велики только при быстрых скачкообразных изменениях скоростей реакций. Оценивание (рис. 8) было достигнуто на приемлемом уровне.

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

Список литературы

[1] Haykin S. Kalman Filtering and Neural Networks. N.Y.: John Wiley and Sons, 2001.

[2] Poznyak A., Sanchez E., Yu W. Dynamic Neural Networks for Nonlinear Control: Identification, State Estimation and Trajectory Tracking. New Jersey: World Scientific, 2001.

[3] Mhamdi A., Marquardt W. Estimation of reaction rates by nonlinear system inversion // Proc. ADCHEM'03. Hong-Kong, 2003. P. 171-176.

[4] Mhamdi A. State and Input Estimation Based on the Theory of Inverse Problems. Dusseldorf: VDI Verlag GmbH, 2004.

[5] Hansen C. Rank-deficient and Discrete Ill-posed Problems: Numerical Aspects of Linear Inversion. SIAM Monographs on mathematical modeling and computation. Philadelphia, 1998.

[6] Farza M., Busawon K., Hammouri H. Simple nonlinear observers for on-line estimation of kinetic rates in bioreactors // Automatica. 1998. Vol. 34, N 3. P. 301-318.

[7] Киллскин В.Е., Торгашов А.Ю., Артюков А.А. Разработка алгоритма оценивания биокинетических параметров процесса разложения перекиси в биореакторе // Сб. трудов 20-й Международной научной конференции "Математические методы в технике и технологиях". ММТТ-20. Ярославль, 2007. Т. 6, секция 12. С. 99-102.

[8] Леонов А.С., Яголл А.Г. Метод L-кривой всегда дает неустранимую систематическую ошибку // Вестник МГУ. Серия 3. Физика. Астрономия. 1997. № 6. С. 17-19.

Поступила в редакцию 23 июня 2008 г в переработанном виде — 1 октября 2008 г.