On a several methods of parametric identification in biokinetics

Full Text

Используемые сокращения и обозначения

МНК – метод наименьших квадратов;

ОДУ – обыкновенные дифференциальные уравнения;

ПраЧ – правая часть (системы дифференциальных уравнений);

a_j и F_j(t) – соответственно коэффициенты и функции, задающие линейную модель (Σ[a_j·F_j(t)]);

В – скорость разбавления (в хемостате);

С_i – текущая концентрация i-го компонента системы (i = 1, 2, …, K);

С_0i – начальная концентрация i-го компонента системы (i = 1, 2, …, K);

D – функция, задающая расстояние между двумя числами;

F – некоторая функция;

G – линеаризующее преобразование;

H(k) – невязка (общего вида), выражающая отклонение экспериментальных значений С_i(t_m) от рассчитанных φ_i(k,t_m);

К – количество компонентов в системе;

k – вектор параметров математической модели (его элементы: k₀, k₁, …, k_j);

N_i – количество экспериментально измеренных (в моменты t_m) значений концентраций для i-го компонента;

S(k) – сумма квадратов отклонений экспериментальных значений С_i(t_m) от рассчитанных φ_i(k,t_m);

S₀ – концентрация субстрата, поступающего в хемостат;

t – время;

t_m – моменты времени, в которые производятся измерения значений С_i (m = 1, 2, …, N_i);

v_i – масштабирующий параметр для i-го компонента;

v_im – масштабирующий параметр для i-го компонента в момент времени t_m;

W_m – «комбинированный» весовой коэффициент (учитывающий все различные причины «взвешивания»);

w_m – частный весовой коэффициент (учитывающий, например, только точность экспериментальных данных);

х – любая независимая переменная (абсцисса);

у – любая зависимая переменная (ордината);

ΔС_i – абсолютная погрешность измерений;

μ – удельная скорость роста микроорганизмов;

σ_m – среднеквадратичное отклонение;

φ_i(t), φ_i(k,t) – рассчитанное по математической модели значение концентрации i-го компонента.

При изучении наук примеры полезнее правил.

Исаак Ньютон¹

ВВЕДЕНИЕ

Кинетикой (и термодинамикой) все чаще и чаще пользуются для выяснения характеристик систем химических и физических процессов, протекающих в организмах. Область исследований, в которых применяют эти дисциплины для решения биохимических проблем, очень обширна и простирается от органической химии до биологии [Bray, White, 1957: p. 5]. В различных областях биологии, включая экологию, достаточно хорошо разработаны формальные методы построения биокинетических математических моделей, чаще всего представляющих собой системы дифференциальных или функционально-дифференциальных уравнений [Romanovskiy et al., 1975; Murray, 2002; Glagolev, Smagin, 2005; Glagolev et al., 2023]. Относительно этих систем можно поставить две принципиально разные задачи. Прямая задача состоит в том, чтобы найти решение системы при заданных значениях параметров (и, возможно, заданных дополнительных условиях, каковыми могут являться, например, начальные или краевые условия в случае систем дифференциальных уравнений). Одним из важнейших типов обратных задач является задача идентификации параметров математической модели («коэффициентная обратная задача»): найти такие численные значения параметров, при которых решение системы будет в некотором смысле наилучшим образом соответствовать имеющимся экспериментальным данным.

Обычно те из систем уравнений (математических моделей), которые являются относительно реалистичными и представляют интерес для практики, оказываются столь сложными, что не имеют аналитического решения. К счастью, были разработаны эффективные численные методы решения прямых задач (см., например [Johnson, 1980: Сh. 3, 4, 7; Engeln-Mullges, Uhlig, 1996: p. 13-154, 423-508; Glagolev et al., 2018; Davis, 2024: p. 173-214, 462-562]). А вот универсальных эффективных методов решения обратных задач, по-видимому, до сих пор не существует. И для надежного решения каждой конкретной задачи приходится подбирать свой метод. Целью настоящей лекции является на конкретных примерах познакомить студентов с некоторыми часто используемыми методами параметрической идентификации математических моделей биокинетики. Впрочем, задачи параметрической идентификации возникают не только в биокинетике, но и практически во всех тех областях науки вообще (и биологии в частности), где используется математическое моделирование (например, [Bard, 1974; Seinfeld, Lapidus, 1974; Ljung, 1987; Mamikhin et al., 2025]). Поэтому рассматриваемые ниже методы могут оказаться пролезными не только для биокинетиков, но для гораздо более широкого круга студентов и специалистов.

ПОСТАНОВКА ЗАДАЧИ ИДЕНТИФИКАЦИИ ПАРАМЕТРОВ МАТЕМАТИЧЕСКОЙ МОДЕЛИ

Простейшая постановка задачи параметрической идентификации

Пусть в системе имеется К различных компонентов, которые могут взаимодействовать между собой, в результате чего с течением времени концентрации (С_i, i = 1, 2, …, K) данных компонентов изменяются. И пусть имеется математическая модель, позволяющая рассчитать эти концентрации в произвольный момент времени t. Обозначим рассчитанные по модели концентрации через φ_i(k,t),² где k – вектор параметров модели (его элементы: k₀, k₁, …, k_j). Например, в качестве системы можно рассмотреть микроорганизмы (концентрация которых – С₁), потребляющие некоторый питательный субстрат (С₂). В результате взаимодействия микробов с субстратом (это взаимодействие выражается в потреблении субстрата) концентрация микробов с течением времени будет расти, а концентрация субстрата – падать (если, конечно, в систему не поступают новые порции субстрата). Математические модели, описывающие этот процесс (см., например, [Pirt, 1975; Romanovskiy et al., 1975: p. 122-140; Glagolev, 2021]), обычно содержат 3-4 параметра: максимальную удельную скорость роста организмов, экономический коэффициент и константу полунасыщения для их роста на данном субстрате, а также, возможно, удельную скорость отмирания.

Если из опыта известны наборы данных, содержащих концентрации компонентов С_i в различные моменты времени t_m, то возможный подход к определению неизвестных констант базируется на минимизации функции S, выражающей сумму квадратов отклонений экспериментальных значений С_i(t_m) от рассчитанных φ_i(k,t_m) [Korobov, Ochkov, 2009: p. 133]:

$S\left(\mathbf{k}\right)=\sum _{i=1}^K\sum _{m=1}^{N_i}{\left[C_i\right(t_m)-{\varphi }_i(\mathbf{k},t_m\left)\right]}^2$, (1)

где N_i – количество экспериментально измеренных значений С_i для i-го компонента (если в рассматриваемой задаче компонент всего один, то мы будем писать не «С₁» и «N₁», а просто «С» и «N»); k – вектор искомых констант. Именно путем варьирования этих констант (оптимизирующих переменных [Gartman, Klushin, 2020: p. 286]) осуществляется минимизация функции S – «невязки между экспериментальными и расчетными данными».

Основное требование к искомым параметрам модели состоит в том, чтобы их изменение оказывало заметное влияние на S (в особенности в окрестностях возможных экстремальных значений). В противном случае найденные оптимальные значения не будут однозначно характеризовать оптимальное решение. Это может привести к получению физически необоснованных и плохо определенных решений, что сведет на нет решение оптимизационной задачи [Gartman, Klushin, 2020: p. 286]. Для читателей, знакомых с разделением задач на «хорошо обусловленные» (well-conditioned) и «плохо обусловленные» (ill-conditioned) [Gerald, Wheatley, 1994: p. 117, 148-154, 172-173; Engeln-Mullges, Uhlig, 1996: р. 8; Glagolev, Sabrekov, 2019], конечно, очевидно, что слабое влияние k на S будет характеризовать именно последние. Подчеркнем, что в данной лекции мы предполагаем задачу минимизации S по k хорошо обусловленной, а к плохо обусловленным задачам обратимся позднее.

Различные типы невязок между расчетными и экспериментальными данными

Однако функция S, задаваемая формулой (1), – не единственная возможная невязка. Более того, иногда она является далеко не лучшей невязкой для минимизации. Пусть, например, ищутся константы для некоторого процесса, в котором 1-й компонент распадается на 2-й и 3-й, но 3-й очень быстро превращается в 4-й, в отличие от 2-го, который ни во что более не превращается

$C_1\begin{array}{c}{}_{\to C_2}\\ {}^{\to C_3\to C_4}\end{array}$

Понятно, что концентрации 3-го компонента могут быть существенно меньше, чем 2-го (т.е. С₃ < < С₂). Поэтому вклад невязки (в суммарную функцию S) между экспериментальными данными и расчетными значениями для 3-го компонента будет существенно меньше, чем для 2-го. Но тогда при минимизации S будет найден такой набор параметров k, при котором экспериментальные данные для 2-го компонента будут воспроизводиться теоретическими кривыми хорошо, а для 3-го – относительно плохо. Чтобы избежать этого, можно использовать функцию

$H\left(\mathbf{k}\right)=\sum _{i=1}^K\sum _{m=1}^{N_i}{\left\{\right[C_i\left(t_m\right)-{\varphi }_i(\mathbf{k},t_m)]/v_i\}}^2$,

где v_i представляет собой некоторый масштабирующий параметр для i-го компонента. Часто 1/v_i² называют «весовым коэффициентом» или просто «весом» (подробнее о весовых коэффициентах будет сказано ниже в соответствующем разделе).

Весовые коэффициенты в этих формулах определяют важность включения в соответствующие критерии рассогласования между экспериментальными данными и переменными, рассчитанными по модели. Чем больше весовой коэффициент, тем более важным считается соответствующее рассогласование (измерение), которое должно влиять на результаты решения идентификационной задачи. Теоретические аспекты решаемой задачи идентификации тесно связаны с математической статистикой и регрессионным анализом, что позволяет более обоснованно определять весовые коэффициенты, в частности, путем использования обратной величины дисперсии экспериментальных измерений (или какой-либо иной характеристики точности измерений; например, можно принять v_i = ΔС_i, где ΔС_i – абсолютная погрешность измерений [Alton, 2011]). Но, вообще говоря, существуют различные методы оценки весовых коэффициентов: эмпирические, статистические и т.п., которые могут выбираться исходя из физического смысла решаемой задачи [Gartman, Klushin, 2020: p. 285]. Таким образом, в общем виде следовало бы писать не v_i, а v_im, чтобы иметь возможность задавать «веса» не только для i-го компонента вообще, но и для каждого m-го экспериментально измеренного значения этого компонента (например, чтобы учесть, что некоторые значения одного и того же i-го компонента могут быть измерены довольно точно, а другие – весьма приближенно).³

В качестве масштабирующего параметра можно взять какую-либо характерную концентрацию i-го компонента. В частности, если v_im = С_i(t_m), то Н(k) будет представлять собой сумму квадратов относительных ошибок (тогда как S – сумма квадратов абсолютных ошибок). Впрочем, относительные ошибки – это тоже не всегда лучший вариант. Действительно, на практике нас обычно не интересует точное совпадение эксперимента и расчета для исчезающе малых концентраций, но именно для них относительные ошибки будут, скорее всего, наибольшими; следовательно, минимизация Н приведет к такому набору k, при котором отклонение (в абсолютных единицах) теоретических кривых от очень малых концентраций будет мало, а от больших – велико⁴.

Наконец, кроме квадратов разностей можно использовать и другие меры расстояния между измеренными и экспериментальными данными, поэтому в общем виде будем записывать невязку следующим образом:

$H\left(\mathbf{k}\right)=\sum _{i=1}^K\sum _{m=1}^{N_i}D\left(\right[C_i\left(t_m\right)-{\varphi }_i(\mathbf{k},t_m)]/v_{im})$. (2)

Здесь D – некоторая подходящая функция для вычисления вышеуказанного расстояния; например, D(х) = х², или D(х) = |х|, или др. Выбор вида функции D, равно как и параметров v_iт, определяется конкретной задачей. И, в частности, статистическими характеристиками экспериментальных данных, используемых в этой задаче (например, если погрешности лучше удовлетворяют двойному экспоненциальному распределению⁵, нежели гауссовскому, то A. Richardson, D. Hollinger рекомендовали в этом случае использовать в (2) не квадрат, а абсолютное значение [Alton, 2011]).

Впрочем, попытки любую возможную невязку описать единой формулой вряд ли имеют большой смысл, ибо, скорее всего, найдется невязка какого-то экзотического вида, которая не удовлетворяет предложенной общей формуле. А если единая формула все-таки будет найдена, то вид ее окажется столь общим, что практического значения иметь не будет. Поэтому гораздо интереснее и важнее вопрос о том, какие конкретные невязки и как часто использовались на практике. Т.Н. Гартман и Д.В. Клушин [Gartman, Klushin, 2020: p. 284-285] указывают 3 основных величины для оценки рассогласования расчетных и экспериментальных данных, которые они называют «критериями»⁶:

• критерий метода наименьших квадратов (МНК), задаваемый формулой (1);
• критерий метода наименьших модулей, задаваемый формулой (2) при D(х) = |х|;
• критерий Чебышева (минимаксный критерий): в частном случае одного компонента (т.е. при K = 1) минимизируется max(|C(t_m) - φ(t_m)|/v_m), где 1 ≤ m ≤ N.

Эти авторы считают, что наиболее употребляемым является критерий МНК, за ним следует минимаксный метод Чебышева и в последнюю очередь – критерий наименьших модулей.

МЕТОДЫ РЕШЕНИЯ ЗАДАЧИ ИДЕНТИФИКАЦИИ ПАРАМЕТРОВ

Типы реализации вычислений параметров моделей

Выше (во «Введении») было заявлено, что пока не существует универсальных эффективных методов решения обратных задач. Однако может показаться, что сказанное ранее (в разделе «Простейшая постановка задачи параметрической идентификации») противоречит этому. Действительно, казалось бы, достаточно записать для конкретной задачи параметрической идентификации невязку (1) и найти такие параметры k₀, k₁, …, k_j, которые обеспечивают минимум невязки, как задача будет решена. Однако никакого противоречия здесь нет, поскольку в общем случае найти этот минимум совсем не просто. Главная трудность состоит в том, что в реальных задачах функция S(k) может иметь несколько локальных минимумов. А универсальных эффективных методов поиска глобального минимума функции при наличии у нее множества локальных минимумов пока не существует. В связи с отсутствием универсальности методов решения выделяют несколько типов задач параметрической идентификации (для работы с которыми применяются те или иные специфические вычислительные методы). В зависимости от того, могут ли функции φ_i(t) быть представлены в аналитическом виде, выделяют два основных типа реализации вычислений параметров.

I. Первый из них применим к ситуациям, когда математическая модель интегрируется в квадратурах⁷. Тогда функции φ_i(t) задаются конкретными аналитическими выражениями, с помощью которых можно осуществить аппроксимацию набора кинетических данных. В свою очередь, здесь можно выделить два способа определения параметров аппроксимации [Korobov, Ochkov, 2009: p. 133].

I.1. Первый и наиболее простой способ можно осуществить, если уравнение зависимости φ_i(t) допускает запись в форме, указывающей на возможность линеаризации кинетических данных. Тогда выбирают соответствующие линеаризующие координаты для кинетической кривой и строят график прямой линии. Корректно построенный график должен удовлетворять критерию минимальной суммы квадратов отклонений, заданной формулой (1). Преимущество этого способа решения состоит в том, что он легко реализуем: практически все математические пакеты предусматривают наличие встроенных функций линейной регрессии. Кроме того, вычисления по этому способу не требуют предварительного определения начальных приближений для искомых параметров [Korobov, Ochkov, 2009: p. 133].

I.2. Но если уравнение линеаризовать не представляется возможным (или это очень сложно), то используются различные нелинейные подходы.

II. Более сложно решить обратную задачу, если математическая модель не интегрируется в квадратурах, т.е. аналитических выражений для функций, выражающих зависимости свойств системы от времени, в этом случае не существует. Однако и в таких ситуациях для идентификации параметров применимы различные подходы [Korobov, Ochkov, 2009: p. 134], которые мы рассмотрим ниже на примерах. Но начнем с простейших задач типа (I.1).

Определение параметров после линеаризующего преобразования модели

Линеаризующие преобразования

Из сказанного выше (об основной трудности минимизации невязки – многоэкстремальности) очевидно, что особое значение приобретают методы решения задачи параметрической идентификации, приводящие к минимизации функции с одним минимумом. В частности, невязка будет иметь лишь один минимум, если математическая модель φ_i(k,t) линейна по k, и в этом случае задача идентификации параметров очень просто решается при помощи, например, линейного метода наименьших квадратов.

К сожалению, анализ литературы [Bray, White, 1957; Romanovskiy et al., 1975; Murray, 2002; Glagolev, 2010] показывает, что линейные модели в биокинетике практически не встречаются. Однако сферу применения линейного МНК можно существенно расширить, поскольку многие часто используемые нелинейные формулы путем некоторых трансформаций преобразуются к линейному выражению [Panikov et al., 1993: p. 39-40; Engeln-Mullges, Uhlig, 1996: p. 215-217; Ryan, 1997: p. 71-79; Davis, 2024: p. 318-321, 325-326, 332-333, 337-339]. Те модели, которые могут быть преобразованы к линейным, Н. Дрейпер и Г. Смит [Draper, Smith, 1981: p. 278] называют «внутренне линейными» моделями (intrinsically linear model⁸), в противоположность «внутренне нелинейным» (intrinsically nonlinear model).

В качестве примера рассмотрим внутренне линейную функцию, очень часто применяющуюся в биокинетике и экологии. Это простейшая математическая модель, выраженная линейным дифференциальным уравнением 1-го порядка [Bray, White, 1957: ch. 6; Murray, 2002: §1.1; Glagolev, 2004; Davis, 2024: p. 340]:

dC/dt = -k₁·C. (3)

Данное уравнение легко интегрируется в квадратурах, и если выполняется начальное условие C(0) = C₀, то решением будет функция C = C₀·е^-k1·t. Очевидно, что функция эта нелинейна (нелинейна и по t, и по k). Пусть мы имеем экспериментальные данные – значения C, измеренные в некоторые моменты времени t_m (отсчитываемые от начального момента времени, т.е. времени начала эксперимента, которое будем считать нулевым): С(t₁), С(t₂), …, С(t_m), …, С(t_N). Существует очень простое преобразование, позволяющее записать экспоненциальную функцию в линейной форме. Это преобразование – логарифмирование. Действительно, после логарифмирования получаем:

ln(C) = ln(C₀) - k₁·t. (4)

Введем новый параметр а₀ = ln(C₀) и новую переменную у = ln(C). Тогда мы будем иметь линейную модель: у = а₀ - k₁·t. Поскольку значения С нам были известны, то и значения у известны, они получаются логарифмированием: у(t₁) = ln(С(t₁)), у(t₂) = ln(С(t₂)), …, у(t_m) = ln(С(t_m)), …, у(t_N) = ln(С(t_N)). Для линейной модели очень просто определить параметры а₀ и k₁ линейным методом наименьших квадратов [Ryan, 1997, p. 72]. Но здесь возникает новая проблема. Дело в том, что найденные этим методом параметры обращают в минимум сумму квадратов отклонений значений преобразованных величин у(t₁), у(t₂), …, у(t_N) от расчетных значений, а не сумму квадратов отклонений измеренных величин С(t₁), С(t₂), …, С(t_N) от соответствующих расчетных [Rumshiskiy, 1971: p. 84].

Весовые коэффициенты

К счастью, возникшая проблема имеет относительно простое решение. Чтобы вычислить C₀ и k₁ при помощи линейного метода наименьших квадратов (примененного к внутренне линейным моделям), нужно использовать специальные веса, учитывающие, какое конкретное преобразование было применено к нелинейной формуле. Однако, кроме того, исходные данные могли быть неравноточными, следовательно, результирующие («комбинированные») веса должны учитывать не только вид преобразования, но и веса, обусловленные погрешностями исходных экспериментальных данных (обозначим их через w_m).

Если распределение ошибок ординат подчиняется нормальному (гауссову) закону, то можно использовать w_m = 1/σ_m², где σ_m – среднеквадратичное отклонение для С(t_m). Этот метод называется инструментальным взвешиванием. Если σ_m неизвестны, то может оказаться полезным применение так называемого «статистического взвешивания»: w_m = 1/С(t_m) [Johnson, 1980: ch. 5]. Наконец, Engeln-Mullges and Uhlig [1996: p. 190] допускают использование w_i = 1/С_i². Когда истинные значения σ_i нам не известны, то такой способ вычисления весов имеет право на существование, поскольку в практике измерений часто оказывается, что σ_i ~ С_i.

Пусть исходная (нелинейная) формула имела вид С = f(k₀, k₁, …, t), где k_j – искомые параметры, C и t – соответственно зависимая и независимая переменные. И пусть G – некоторое преобразование (функция), которое, будучи применено к рассматриваемой нелинейной формуле, превращает ее в линейное выражение: G(y) = G(f(k₀, k₁, …, t)) = Σ[a_j·F_j(t)], где F_j(t) – некоторые функции, не содержащие искомых параметров; a_j – некоторые новые параметры (выражающиеся – возможно, нелинейно – через k₀, k₁, …). Тогда формула «комбинированных» весов (W_m) такова: W_m = w_m·(dG/dC)^-2|_m [Engeln-Mullges, Uhlig, 1996, p. 216]. Например, для рассмотренной выше линеаризации логарифмированием будем иметь:

W_m = w_m·{d[ln(C)]/dC}^-2|_m = w_m·(1/C)^-2|_m = w_m·C_m².

Метод приближенного линеаризующего преобразования

Обычно в литературе метод определения параметров с помощью линеаризации данных излагается именно в приложении к уравнениям реакции 1-го порядка (3), (4) или к другим простым формулам, которые с помощью некоторого преобразования можно привести в точности к линейному выражению (см., например, [Draper, Smith, 1981: § 5.2; Engeln-Mullges, Uhlig, 1996: p. 215-216; Korobov, Ochkov, 2009: p. 134-137; Davis, 2024: p. 318-321, 325, 332-333, 337-339]). Но чтобы читатель понимал большую универсальность данного метода, мы немного отойдем от указанной традиции и рассмотрим чуть более сложный пример, следуя Васильевой и др. [Vasil’eva et al., 1995]. Эти авторы рассмотрели потребление субстрата (4-хлоранилина) растущей биомассой микроорганизмов. Вообще говоря, в данном случае мы будем иметь систему как минимум из двух уравнений: одно – для описания динамики роста биомассы, а другое – для потребления субстрата:

dC₁/dt = μ(C₂)·C₁, dС₂/dt = -μ(C₂)·C₁/k₁, (5)

где C₁ и C₂ – соответственно концентрации микроорганизмов и питательного субстрата (в данном случае – хлоранилина); μ – удельная скорость роста (зависящая от C₂, а в некоторых моделях микробиологической кинетики – еще и от C₁); k₁ – экономический коэффициент. В зависимости от конкретного вида функции μ эти уравнения могут или не могут быть проинтегрированы в квадратурах. Для простоты Васильева и др. [Vasil’eva et al., 1995] рассматривают простейший случай постоянства удельной скорости роста: μ(C₂) = k₂ = const. Конечно, в реальности эта скорость не может быть всегда постоянной и, по крайней мере, должна обратиться в 0 при C₂ = 0, ибо очевидно, что микроорганизмы не могут расти, если у них совсем не будет еды. Но из микробиологической кинетики хорошо известно (см., например, [Pirt, 1975: § 2.6-2.7; Atkinson, Mavituna, 1983: p. 207; Romanovskiy et al., 1975: p. 124; Panikov et al., 1992]), что в довольно широком диапазоне значений C₂, действительно, μ остается практически постоянной – почти на уровне максимальной удельной скорости роста. Итак, если принять данное упрощение, то первое уравнение системы уже никак не связано со вторым и легко интегрируется при начальном условии C₁(0) = C₀₁:

C₁ = C₀₁·ехр(k₂·t).

Полученное выражение мы можем подставить во второе уравнение, что даст

dC₂/dt = -k₂·[C₀₁·ехр(k₂·t)]/k₁.

Проинтегрируем это уравнение при начальном условии С(0) = С₀₂:

С₂ = С₀₂ - (С₀₁/k₁)·[ехр(k₂·t) - 1]. (6)

Для относительно больших значений t, когда ехр(k₂·t) >> 1, будет выполняться следующее приближенное соотношение:

С₀₂ - С₂ ≈ (С₀₁/k₁)·ехр(k₂·t).

Оно легко линеаризуется в полулогарифмических координатах:

ln(С₀₂ - С₂) ≈ ln(С₀₁/k₁) + k₂·t. (7)

Т.е. если мы обозначим y = ln(С₀₂ - С₂) и а₁ = ln(С₀₁/k₁), то в координатах (t, y) получим прямую линию, задаваемую уравнением

y ≈ a₁ + k₂·t. (8)

Но, обрабатывая реальные экспериментальные данные, следует помнить о тех приближениях, которые мы сделали, чтобы получить (7). Во-первых, мы приняли μ = const, а это условие не может выполняться при малых С₂, поскольку микробам не будет хватать еды, чтобы поддерживать такую же большую удельную скорость роста, которая была в условиях изобилия пищи. Во-вторых, (7) получается из (6) только когда ехр(k₂·t) >> 1. В принципе, зная из микробиологии типичные значения максимальной удельной скорости роста (k₂) для микроорганизмов, можно приближенно оценить: при каких значениях t формулу (7) применять еще нельзя, а при каких – уже можно. Но тут можно сильно ошибиться (поскольку в зависимости от вида микробов и конкретного субстрата их максимальные удельные скорости роста могут различаться на порядки – ср., например, [Pirt, 1975: § 13.2; Lee, 1981: p. 159; Atkinson, Mavituna, 1983: p. 207; Dorofeev et al., 1992]). А кроме того, в этом нет смысла, поскольку область допустимого применения соотношения (7) будет видна из самого графика: экспериментальные данные, полученные при малых t, не лягут на прямую (из-за невыполнения второго допущения), равно как не лягут на нее и точки, полученные при больших t (здесь отклонение от прямой произойдет по вине первого допущения).

Рассмотрим вышеизложенное на численном примере, материал для которого возьмем все из той же статьи [Vasil’eva et al., 1995]. В Табл. 1 приведены экспериментальные данные⁹, которые мы будем использовать. Эти же данные (за исключением точки для t = 0) показаны на Рис. 1a в полулогарифмических координатах. Из рисунка видно, что полный набор данных действительно аппроксимируется прямой довольно плохо. Причем как раз (как мы теоретически предположили выше) именно из-за того, что от прямой резко отклоняются несколько самых правых и левых точек. А вот если отбросить по три крайние точки справа и слева, оставив только девять «центральных» точек, то прямая пройдет через них вполне удовлетворительно¹⁰. Итак, для дальнейшей работы оставляем только эти точки.

 

Таблица 1. Динамика потребления 4-хлоранилина культурой Paracoccus denitrificans.

Table 1. Dynamics of 4-chloroaniline consumption by a microbial culture Paracoccus denitrificans.

Время (t, часы)	Концентрация (С2, мг/л)		Время  (t, часы)	Концентрация (С2, мг/л)
0	20.0		69.0	14.9
16.6	19.3		75.4	11.2
34.0	19.1		80.0	7.6
39.5	18.8		82.8	5.6
46.9	19.0		85.5	3.3
51.5	18.5		89.2	1.7
59.8	18.3		92.9	1.1
65.3	16.5		98.4	0.1

 

Примечание: Полный набор экспериментальных данных представлен совокупностью точек обоих типов («●» и «о»); ________ - прямая, построенная МНК только по «о»; __ __ __ - прямая, построенная МНК по всем точкам. Note: The complete set of experimental data is repre-sented by points of both types («●» and «о»); the straight lines were constructed using least-squares fitting ________ - based on «о» points only; __ __ __ - based on all points.

Примечание: ● – экспериментальные данные; Кривые построены ▬▬▬ - по формуле (6) (по точкам «←»); _ _ _ _ _ - по формуле (6) (по точкам «↑»); _________ - по формуле (11). Note: ● – experimental data; The curves were constructed ▬▬▬ - by using formula (6) (from points «←»); _ _ _ _ _ - by using formula (6) (from points «↑»); _________ - by using formula (11).

Рис. 1. Моделирование динамики потребления 4-хлоранилина по данным [Vasil’eva et al., 1995].

Fig. 1. Modeling of the 4-chloroaniline consumption dynamics with data from [Vasil’eva et al., 1995].

 

Параметры формулы (8) в MATLAB можно найти при помощи функции polyfit. Подробно она описана в интерактивной справочной системе MATLAB и в литературе (см., например, [Chen et al., 1999: §4.4]), а сейчас мы лишь продемонстрируем использование этой функции для решения нашей конкретной задачи. Введем абсциссы и ординаты экспериментальных точек, а потом вызовем polyfit (задав в ней третий аргумент равным 1, что соответствует поиску коэффициентов полинома 1-го порядка, т.е. линейной функции):

>> t=[46.9 51.5 59.8 65.3 69 75.4 80.0 82.8 85.5];

>> y=[-.02 .38 .53 1.26 1.63 2.17 2.52 2.67 2.82];

>> P=polyfit(t,y,1),

мы сразу получим результат:

P = 0.0770 -3.7218

Он означает, что в (8) a₁ = -3.7218, k₂ = 0.0770. Поскольку на основании (7) a₁ = ln(С₀₁/k₁), то С₀₁/k₁ = ехр(a₁) = ехр(-3.7218) ≈ 0.0242, и, если нам известно С₀₁, мы можем найти k₁.

В заключение этого примера оценим, выполняется ли в данном случае допущение о том, что прямая строится на таком интервале времени, где единицей в (6) уже можно пренебречь (и потому, собственно говоря, получается именно прямая). Как мы помним, это допущение состоит в том, что ехр(k₂·t) >> 1. Крайняя левая точка, которую мы используем для построения прямой, соответствует значению t = 46.9 час. (см. Табл. 1 и Рис. 1a). Для нее имеем: ехр(0.077·46.9) ≈ 37. Понятно, что для следующих точек, расположенных правее, условие тем более будет выполняться.

Способ подстановки

Простейший случай

Когда математическая модель представляет собой явную аналитическую зависимость С(t) = φ(k,t), метод подстановки идейно весьма прост. Пусть, как и раньше, мы имеем N пар экспериментальных данных: {t₁, С(t₁)}; {t₂, С(t₂)}; …; {t_m, С(t_m)}; …; {t_N, С(t_N)}. Тогда для определения параметров k = {k₀, k₁, …, k_j} каждая из пар данных подставляется в уравнение модели, в результате чего мы будем иметь систему уравнений:

С(t₁) = φ(k₀, k₁, …, k_j, t₁), С(t₂) = φ(k₀, k₁, …, k_j, t₂), …, (9)

где t_i и С(t_i) – известные (из эксперимента) числа, а k₀, k₁, …, k_j – неизвестные величины, которые и следует найти, решая данную систему уравнений. Если нам нужно определить j неизвестных констант, то необходимо иметь не менее j пар экспериментальных данных (т.е. не менее j уравнений). Но обычно данных гораздо больше, поэтому оказывается возможным определить несколько наборов констант, которые затем усредняются. Справедливости ради заметим, что в таком простейшем виде этот способ потерял свою актуальность с появлением математических пакетов [Korobov, Ochkov, 2009: p. 134], способных осуществлять нелинейную регрессию. Тем не менее для полноты картины мы приведем пример применения этого метода, опять воспользовавшись данными Табл. 1 и уравнением (6). Однако это уравнение запишем в виде С₂ = С₀₂ - k₃·[ехр(k₂·t) - 1], где k₃ = С₀₁/k₁, поскольку, не зная значение С₀₁, мы можем определить только численное значение отношения С₀₁/k₁, но не сам коэффициент k₁. Итак, наше уравнение содержит два неизвестных параметра (k₃ и k₂), следовательно, нужно создать систему из двух уравнений, т.е. взять две пары экспериментальных данных. Но какие конкретно пары взять?

Из сказанного выше вроде бы следует, что взять можно какие угодно пары; главное перебрать все возможные пары – найти различные k₃ и различные k₂ для всех возможных пар {t_m, С₂(t_m)}, {t_n, С₂(t_n)} (m ≠ n) и усреднить все полученные значения k₃ и k₂. Однако эта рекомендация в некоторых случаях может оказаться плохой. Например, если экспериментальных точек не очень много и среди них есть хотя бы пара значений, измеренных с достаточно большими погрешностями, то k₃ и k₂, найденные по этим значениям, могут оказаться настолько «плохими», что при усреднении сильно «испортят» всю картину. На практике иногда оказывается более полезным выбрать лишь j (т.е. в нашем примере – две) наиболее характерные точки и записать систему уравнений только для них. На Рис. 1 показаны полученные в эксперименте значения из работы [Vasil’eva et al., 1995]. Напомним, что уравнение (6) выведено из (5) в предположении того, что значения С₂ достаточно велики – только тогда можно считать, что μ(С₂) = сonst = k₂. Если проанализировать функцию (6) стандартными средствами математического анализа, то можно легко установить, что это выпуклая функция (т.е. d²С₂/dt² < 0). Но из Рис. 1b ясно видно, что примерно при t = 86 находится точка перегиба, правее которой выпуклость сменяется вогнутостью. Следовательно, три крайние справа точки (при t = 89.2, 92.9, 98.4) нельзя использовать в системе уравнений для определения k₃ и k₂, поскольку для этих точек уравнение (6) вообще не выполняется. Далее заметим, что кривая (6) в любом случае пройдет через крайнюю левую точку (при t = 0): С₂ = С₀₂ - k₃·[ехр(k₂·0) - 1] = С₀₂ - k₃·(1 - 1) = С₀₂. В принципе любые другие точки (т.е. кроме крайней первой и трех последних) можно выбирать для формирования нашей системы уравнений. Но представим, что мы выбрали 3-ю и 4-ю точки (т.е. {34, 19.1} и {39.5, 18.8}, на Рис. 1b они указаны вертикальными стрелками). Если записать систему уравнений для этой пары, то решения (а это будут k₃ = 0.304735 и k₂ = 0.040429) обеспечат прохождение кривой С₂(t) точно через точки {34, 19.1}, {39.5, 18.8} и… и далее кривая пойдет как Бог на душу положит (см. Рис. 1b). Интуитивно понятно, что если мы хотим обеспечить аппроксимацию экспериментальных данных в целом, то нужно, чтобы аппроксимирующая кривая прошла через крайние точки и точку посередине (имеется в виду середина оси ординат!). Но через крайнюю левую точку кривая, как только что было показано чуть выше, проходит «автоматически», значит, нужно выбрать еще четвертую точку справа (т.е. {85.5, 3.3} – именно она является «крайней», в которой уравнение (6) еще выполняется) и, например, точку {75.4, 11.2}. Тогда система уравнений метода подстановки для определения k₃ и k₂ будет такой:

3.3 = 20 - k₃·[ехр(k₂·85.5) - 1], 11.2 = 20 - k₃·[ехр(k₂·75.4) - 1].

Эту систему можно численно решить при помощи различных программных средств¹¹. Ее решение таково: k₃ = 0.076631, k₂ = 0.063026. Как видно из Рис. 1b, при таких значениях параметров расчётная кривая вполне удовлетворительно описывает экспериментальные данные в той области, где справедливо уравнение (6).

Обобщение на неявные функции

В предыдущем разделе мы предполагали, что уравнение математической модели может быть записано в виде С(t) = φ(k,t). Но поскольку исходно кинетическая модель обычно представляет собой дифференциальное уравнение (или систему таких уравнений), как, например, (3) или (5), то нельзя ожидать, что решение этих уравнений всегда будет иметь столь простой вид. Какие еще могут быть формы решений (если, конечно, решение в принципе может быть записано в аналитическом виде)? С логической точки зрения есть еще две возможности. Во-первых, решение дифференциального уравнения даст функцию t от С, из которой в явном виде выразить С (как некоторую функцию t) невозможно. И, во-вторых, возможен более общий случай – решением дифференциального уравнения будет некоторое выражение, связывающее С и t, но при этом из него невозможно ни С выразить через t в явном виде, ни, наоборот, t через С: F(С,t,k) = 0. Однако поскольку систему уравнений метода подстановки мы все равно решаем численно, а не аналитически, то такое усложнение связи С с t не играет вообще никакой роли – система метода подстановки (9) будет иметь вид

F(С(t₁), t₁, k₀, k₁, …, k_j) = 0, F(С(t₂), t₂, k₀, k₁, …, k_j) = 0, …, (10)

и решаться теми же самыми численными методами, что и (9), т.к. решить-то ее нам нужно относительно k₀, …, k_j, а не относительно С или t. Чтобы сказанное было яснее, приведем конкретный пример, опять обратившись к системе (5).

Для описания зависимости μ(C₂) в микробиологической кинетике очень часто используется уравнение Моно (см., например, [Pirt, 1975: § 2.6; Atkinson, Mavituna, 1983: p. 204-210; Romanovskiy et al., 1975: p. 123; Lokshina et al., 2019]): μ(C₂) = k₂·C₂/(k₄ + C₂), где k₄ – «константа полунасыщения». Система (5) с μ(C₂), заданной по уравнению Моно, имеет следующее аналитическое решение:

$\frac{k_4+C_{02}+k_3}{k_2\cdot (C_{02}+k_3)}\cdot \ln \left(\frac{C_{02}-C_2}{k_3}+1\right)-\frac{k_4}{k_2\cdot (C_{02}+k_3)}\cdot \ln \left(\frac{C_2}{C_{02}}\right)=t$. (11)

Данное выражение задает однозначную связь C₂ с t (заданному значению t соответствует некоторое значение C₂), т.е. C₂ можно рассматривать как функцию от t.

Но это – неявная функция, поскольку мы не можем выразить из (11) C₂ в виде С₂(t) = φ(k₂, k₃, k₄, t) [Grossman, Turner, 1974: App. B2]. Однако, повторим, такой сложный вид связи C₂ с t никак не влияет на возможность применения метода подстановки. Запишем уравнения этого метода в данном случае. Прежде всего отметим, что поскольку теперь математическая модель содержит 3 искомых параметра (k₂, k₃ и k₄), то мы должны записать 3 уравнения, для чего, следовательно, нам понадобятся 3 экспериментальные точки. Обратите внимание: интегрируя систему (5) с μ(C₂), заданной по уравнению Моно, мы не использовали никаких дополнительных ограничительных предположений, следовательно, теперь мы можем использовать любые экспериментальные точки, в том числе и три крайние справа, которые мы были вынуждены отбрасывать ранее. Более того, интуитивно понятно, что для определения k₄ как раз нужно использовать какую-либо из этих правых точек¹². Итак, выберем следующие 3 точки: {51.5, 18.5}, {85.5, 3.3} и {98.4, 0.1}. Тогда система уравнений метода подстановки для определения k₂, k₃ и k₄ будет такой:

$\left\{\begin{array}{l}\frac{k_4+20+k_3}{k_2\cdot (20+k_3)}\cdot \ln \left(\frac{20-18.5}{k_3}+1\right)-\frac{k_4}{k_2\cdot (20+k_3)}\cdot \ln \left(\frac{18.5}{20}\right)=51.5\\ \frac{k_4+20+k_3}{k_2\cdot (20+k_3)}\cdot \ln \left(\frac{20-3.3}{k_3}+1\right)-\frac{k_4}{k_2\cdot (20+k_3)}\cdot \ln \left(\frac{3.3}{20}\right)=85.5\\ \frac{k_4+20+k_3}{k_2\cdot (20+k_3)}\cdot \ln \left(\frac{20-0.1}{k_3}+1\right)-\frac{k_4}{k_2\cdot (20+k_3)}\cdot \ln \left(\frac{0.1}{20}\right)=98.4\end{array}\right.$

Ее решение: k₄ ≈ 6.975717, k₃ ≈ 0.020932, k₂ ≈ 0.112743. Как видно из Рис. 1b, при таких значениях параметров расчётная кривая вполне удовлетворительно описывает все экспериментальные данные, включая и крайние правые точки.

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

Если дифференциальные уравнения кинетической модели не могут быть проинтегрированы аналитически, то мы не имеем явной зависимости С(t) = φ(k,t) или хотя бы неявного задания функции С(t) в виде F(С(t), t, k₀, k₁, …, k_j) = 0. В этом случае систематически осуществить способ подстановки мы не можем. Но если модель такова, что кривые имеют экстремумы, то в точках экстремумов производные равны нулю¹³, следовательно, в этих точках можно в уравнения (уже не содержащие производных!) подставить экспериментальные данные. Конечно, таких ситуаций почти всегда будет меньше, чем коэффициентов, которые необходимо определить, поэтому обычно все коэффициенты этим способом определить мы не сможем. Однако даже определение одного коэффициента или хотя бы связи между некоторыми коэффициентами может существенно облегчить решение общей обратной задачи (например, если всего надо определить 2 коэффициента и 1 мы вычислили по способу подстановки в точке экстремума, то остается определить еще лишь 1 коэффициент, что легко сделать хотя бы простейшим подбором). Рассмотрим пример не слишком простой модели.

Ю.М. Романовский с соавт. [Romanovskiy et al., 1975: р. 161-162] привели систему уравнений роста биомассы в проточном культиваторе с учетом инерционности процесса синтеза рибосом:

 dС₁/dt = k₂·С₂·С₁·С₃/[(k₀+С₂)·k₃] - В·С₁, (12)

 dС₂/dt = - k₂·С₂·С₁·С₃/[(k₀+С₂)·k₃·k₅] - В·(S₀ - С₂), (13)

dС₃/dt = [k₃·(k₀+С₂)/(k₁+С₂) - С₃]/k₄, (14)

где В – скорость разбавления; k₀, k₁, k₂, k₃, k₄ и k₅ – различные кинетические параметры¹⁴ модели (в частности, k₅ – экономический коэффициент, показывающий, какая часть поглощенного субстрата идет на построение биомассы; k₂ – максимальная удельная скорость роста); С₁ и С₂ – концентрации биомассы микробов и питательного субстрата соответственно; С₃ – количество рибосом, отнесенное к единице биомассы (концентрация рибосомальной РНК в клетке); S₀ – концентрация субстрата, поступающего в культиватор; t – время (независимая переменная). Параметры B и S₀ мы задаем (а если их менять невозможно, то в любом случае для данного конкретного ферментера они известны). Понятно, что определить нам нужно коэффициенты k₀, k₁, k₂, k₃, k₄ и k₅ по измеряемой в эксперименте динамике зависимых переменных С₁, С₂ и С₃.

Экспериментальные данные возьмем из [Ugodchikov, 1975: р. 194] – мы приводим их в Табл. 2. Важно отметить, что эти экспериментальные данные получены для частного случая периодической культуры. Следовательно, для их описания в уравнениях (12) и (13) нужно положить В = 0. Но нам это только на руку, ведь тогда уравнения становятся еще проще!

 

Таблица 2. Динамика роста культуры Escherichia coli (шт. М-17) на глюкозо-солевой среде [Ugodchikov, 1975: р. 194].

Table 2. Growth dynamics of an Escherichia coli (strain M-17) culture on a glucose-mineral medium [Ugodchikov, 1975: p. 194].

Время (t, час)	Концентрация
	биомассы (С1, мг/л)	субстрата (С2, мг/л)	РНК (С3, гРНК/гБиомассы)
0.0	53	130	0.040
0.5	59	130	0.041
1.0	65	130	0.057
1.5	80	130	0.068
2.0	106	126	0.089
2.5	150	120	0.107
3.0	213	109	0.127
3.5	314	90	0.140
4.0	479	63	0.140
4.5	711	36	0.134
5.0	899	0	0.122
6.0	947	0	0.116

 

Легко заметить, что уравнение (14) описывает кривую с экстремумом:

dС₃/dt = 0 при [k₃·(k₀+С₂)/(k₁+С₂) - С₃]/k₄ = 0 => k₃·(k₀+С₂^е)/(k₁+С₂^е) = С₃^е, (15)

где верхним индексом «е» мы обозначили значения переменных в точке экстремума. Из табл. 2 видим, что экстремальная (максимальная) концентрация рибосомальной РНК достигается, вероятно, где-то в интервале от 3.5 до 4 час. Кажется естественным в качестве С₃^е и С₂^е взять средние арифметические значения:

С₃^е = [С₃(3.5) + С₃(4)]/2 = (0.14 + 0.14)/2 = 0.14, С₂^е = [С₂(3.5) + С₂(4)]/2 = (90 + 63)/2 = 76.5 мг/л.

Подставляя эти численные значения в (15), имеем соотношение

k₃·(k₀+76.5)/(k₁+76.5) = 0.14, (16)

задающее связь между k₃, k₀ и k₁. Как видим, определить какой-либо коэффициент нам пока не удалось, но если в дальнейшем мы сможем вычислить, например, k₀ и k₁, то из (16) без труда определим k₃.

Исключение времени как независимой переменной

Прием исключения времени как независимой переменной можно эффективно применять при решении некоторых обратных задач [Korobov, Ochkov, 2009: p. 141]. Но сразу подчеркнем, что прием этот – весьма частный. Впрочем, именно в кинетических задачах его оказывается возможным применить не так уж редко, поскольку в них обычно выполняются законы сохранения для тех или иных компонентов, а такие законы, как правило, и представляют собой форму уравнений с исключенным временем. С другой стороны, следует отметить, что метод применим не только в обратных задачах кинетики. Более широко можно говорить про метод исключения независимой переменной вообще (а не только времени). Тем не менее мы рассмотрим его именно в приложении к кинетической задаче и исключать будем именно время, которое в таких задачах как раз и является независимой переменной.

Суть описываемого метода состоит в следующем. Чаще всего кинетическую систему стараются записать в виде автономных дифференциальных уравнений, т.е. таких, в правую часть (ПраЧ) которых независимая переменная (время) не входит. В этом случае, разделив одно уравнение системы на другое, можно получить уравнение, в которое время не будет входить вообще¹⁵. Полученное уравнение может оказаться существенно проще для интегрирования, чем исходные, поэтому есть смысл определять значения параметров для него. Правда, обычно таким образом можно определить далеко не все параметры исходной системы, поскольку лишь отношения некоторых уравнений порождают уравнение, существенно более простое, чем исходные. Но в любом случае идентификация хоть какого-то количества параметров упрощает дальнейшую идентификацию остальных. Покажем на конкретном примере, как действует описанный метод. В качестве примера опять обратимся к системе (12)-(14).

Очевидно, что если разделить уравнение (12) на (13), то получим гораздо более простое уравнение¹⁶ с разделяющимися переменными (С₁ и С₂), которое очень легко интегрируется в квадратурах:

dС₁/dС₂ = -k₅ => С₁ - С₀₁ = k₅·(С₀₂ - С₂), (17)

здесь интегрирование было проведено в границах от начальных концентраций (С₀₁ и С₀₂ в момент времени t = 0) до текущих (С₁ и С₂ в произвольный момент времени t > 0). Из полученного соотношения видно, что если построить зависимость у = (C₁‑C₀₁) от х = (C₀₂-C₂), то экспериментальные данные должны лечь на прямую, наклон которой составит k₅ (если, конечно, эти данные действительно удовлетворяют нашей математической модели). На рисунке 2 выполнено такое построение. Из рисунка видно, что данные измерений действительно прекрасно (r² ≈ 0.996) ложатся на прямую. При помощи линейной регрессии можно легко определить, что k₅ ≈ 6.73.

 

Рис. 2. Определение экономического коэффициента в координатах (C₀₂-C₂), (C₁‑C₀₁).

Fig. 2. Determination of the yield constant in the coordinates (C₀₂-C₂), (C₁‑C₀₁).

 

Итак, при помощи метода исключения независимой переменной нам удалось из первых двух уравнений математической модели (12)-(14) получить весьма простое дифференциальное уравнение, интегральной кривой которого является прямая. И по наклону этой прямой мы идентифицировали численное значение одного из параметров модели – экономического коэффициента. К сожалению, исключение времени путем деления уравнений (12) или (13) на уравнение (14) приводит к довольно сложным выражениям. Так что каких-либо простых путей для определения данным методом значений других параметров модели не просматривается (напомним, что нам надо еще идентифицировать k₀, k₁, k₂, k₃ и k₄). Поэтому перейдем к рассмотрению других методов.

Упрощение модели в предельных случаях

Строго говоря, этот метод тоже весьма частный. Но опыт показывает, что почти всегда можно найти ситуацию, когда сложная модель может быть хорошо аппроксимирована простой моделью, допускающей аналитическое интегрирование в такой форме, к которой можно подобрать несложное линеаризующее преобразование. Мы нечто подобное уже делали выше в разделе «Метод приближенного линеаризующего преобразования». Ведь уравнение (6) мы не смогли линеаризовать и рассматривали предельный случай больших t, поскольку при этом становится возможной очень простая линеаризация в полулогарифмических координатах. Но в указанном разделе мы аналитически проинтегрировали дифференциальное уравнение, получили в качестве его решения нелинейную зависимость (6) и в дальнейшем работали с приближением к этому аналитическому решению. Иначе говоря, приближение в нашем анализе возникало после аналитического интегрирования. А теперь на примере модели (12)-(14) обсудим возможность приближенного интегрирования более сложных уравнений, т.е. воспользуемся приближением до интегрирования.

Обычно относительно сложные системы уравнений оказывается возможным аналитически проинтегрировать в предельных случаях очень малых или очень больших значений переменных. Пусть С₂ – велико. Самые большие значения концентрации субстрата, очевидно, соответствуют начальному периоду роста биомассы, когда микробы только-только начали потреблять субстрат. Если выполняется условие

k₀ < < С₂ >> k₁,

то k₀ + С₂ ≈ С₂ и k₁ + С₂ ≈ С₂, следовательно, уравнение (12) приближенно можно заменить более простым уравнением:

dС₁/dt ≈ k₂·С₁·С₃/k₃.

Кажется, что проинтегрировать это уравнение нельзя, поскольку оно содержит две независимых переменных – С₁ и С₃. Но ведь и уравнение (14) при больших С₂ упрощается:

dС₃/dt ≈ (k₃ - С₃)/k₄.

Это уравнение легко интегрируется в квадратурах: С₃ ≈ k₃ - (k₃-С₀₃)·exp(-t/k₄). Следовательно, на начальном участке динамической кривой биомасса должна описываться вот таким простым дифференциальным уравнением:

dС₁/dt ≈ k₂·С₁·[1 - (1-С₀₃/k₃)·exp(-t/k₄)],

легко интегрируемым в квадратурах.

К сожалению, описываемый метод обладает некоторыми недостатками, в связи с чем его применение крайне ограничено. Во-первых, не всегда понятно, где проходит тот рубеж, после которого упрощенной моделью пользоваться уже нельзя. А во-вторых, в область упрощенной модели часто попадает слишком мало экспериментальных данных, поэтому надежная идентификация параметра невозможна. И если первый недостаток можно легко преодолеть (это мы уже обсуждали выше в разделе «Метод приближенного линеаризующего преобразования»; в частности, см. Рис. 1a), то со вторым недостатком какими-то математическими средствами содержательно бороться нельзя и для надежной идентификации следует лишь наращивать объем экспериментальной информации.

Представление модели в виде ряда Тейлора

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

На Рис. 3 в графическом виде представлена часть экспериментальных данных из Табл. 2. Легко заметить, что доля РНК (Рис. 3б) в промежутке времени примерно с 1.5 до 3.5 час. возрастает практически линейно. Разложим функцию R(t) в ряд Тейлора в окрестности точки с абсциссой, являющейся центром этого интервала (т.е. 2 час.), и отбросим все члены ряда, кроме двух первых:

С₃(t) ≈ С₃(2) + (t-2)·dС₃/dt|_t=2.

 

Рис. 3. Динамика биомассы (а) и РНК (б) в культуре Е. coli.

Fig. 3. Dynamics of biomass (a) and RNA (b) in an E. coli culture.

 

Выражение для dС₃/dt нам известно – это ур. (14), но чтобы получить его значение при t = 2, нужно подставить в (14) численные значения С₂(2) и С₃(2). Из Табл. 2 видим, что С₂(2) = 126 мг/л и С₃(2) = 0.089, поэтому

С₃(t) ≈ 0.089 + (t-2)·[k₃·(k₀+126)/(k₁+126) - 0.089]/k₄.

С другой стороны, мы можем провести по 5 экспериментальным точкам (как на Рис. 3б – в интервале t от 1.5 час. до 3.5 час.) эмпирическую прямую. Будем искать прямую следующего вида: С₃(t) ≈ 0.089 + (t-2)·а, где а – неизвестный коэффициент, который, собственно говоря, и требуется определить. Наилучшим образом (в смысле наименьших квадратов) экспериментальные данные аппроксимируются такой прямой:

 С₃(t) = 0.089 + 0.036·(t-2),

следовательно,

[k₃·(k₀+126)/(k₁+126) - 0.089]/k₄ = 0.036. (18)

Подчеркнем, что представление С₃(t) в виде ряда Тейлора пока не дало нам возможности определить численное значение какого-то конкретного коэффициента, но позволило найти дополнительную связь между коэффициентами. Понятно, что каждая такая связь фактически сокращает количество неизвестных коэффициентов на один. Действительно, предположим, что в дальнейшем мы сможем найти каким-либо способом коэффициенты k₀, k₁ и k₃. Тогда из (18) мы сможем вычислить k₄.

Теперь найдем еще одну связь между коэффициентами, обратившись на сей раз к функции С₁(t). Легко заметить, что биомасса бактерий (Рис. 3а) в промежутке времени примерно с 2 до 5 час. возрастает практически экспоненциально (обратите внимание: на этом рисунке масштаб по оси ординат – логарифмический, а ось абсцисс – обычная, следовательно, экспонента в этих осях будет представлена прямой линией). Разложим функцию ln(С₁) в ряд Тейлора (по t) в окрестности точки с абсциссой, являющейся центром этого интервала (т.е. 3.5 час.), и отбросим все члены ряда кроме двух первых:

ln(С₁(t)) ≈ ln(X(3.5)) + (t-3.5)·dln(С₁)/dt|_t=3.5.

Однако откуда же мы возьмем выражение для dln(С₁)/dt? Прежде всего воспользуемся известным правилом: dln(С₁)/dt = С₁^-1·dС₁/dt. А выражение для dС₁/dt нам известно – это уравнение (12), следовательно (не забываем, что В = 0!):

dln(С₁)/dt = k₂·С₂·С₃/[(k₀+С₂)·k₃].

Но чтобы получить его значение при t = 3.5, нужно подставить численные значения С₁(3.5), С₂(3.5) и С₃(3.5). Из таблицы 2 видим, что С₁(3.5) = 314, С₂(3.5) = 90 мг/л и С₃(3.5) = 0.14, поэтому

ln(С₁) ≈ 5.749 + (t-3.5)·k₂·12.6/[(k₀+90)·k₃].

С другой стороны, мы можем провести по 7 экспериментальным точкам (как на Рис. 3а – в интервале t от 2 до 5 час.) эмпирическую экспоненту. Будем искать экспоненту следующего вида: ln(С₁) ≈ 5.749 + (t-3.5)·а. Наилучшим образом (в смысле наименьших квадратов) экспериментальные данные аппроксимируются вот такой экспонентой:

ln(С₁) ≈ 5.749 + (t-3.5)·0.7383,

следовательно,

k₂·12.6/[(k₀+90)·k₃] = 0.7383 < => k₂/[(k₀+90)·k₃] = 0.0586. (19)

Итак, мы нашли еще одну связь между коэффициентами.

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

Линеаризация с численным интегрированием кинетических данных

У читателя могло сложиться мнение о том, что способ определения кинетических параметров, основанный на линеаризации данных, имеет ограниченное применение. Покажем, однако, что круг задач, решаемых этим способом, достаточно широк и может включать те из них, где математическая модель представляется в виде нелинейных систем дифференциальных уравнений [Korobov, Ochkov, 2009: p. 142]. Хотя, конечно, речь не идет о том, что этим методом можно определять параметры любых нелинейных систем дифференциальных уравнений. Вообще говоря, к сожалению, для работы этим методом нет каких-то универсальных алгоритмов (или, по крайней мере, авторам такие методы не известны). Как сейчас станет ясно читателю, возможность успешного использования метода в той или иной ситуации зависит от уровня знаний и смекалки исследователя, пытающегося применить его. Обратимся опять к нашему примеру (12)-(14).

Очевидно, что уравнение (13) – нелинейное. Мы сейчас произведем некоторые тождественные преобразования этого уравнения, поскольку для успешного применения излагаемого метода будет лучше, если все переменные имеют примерно один порядок. В этой связи представляется удобным, если все переменные будут находиться в интервале от 0 до 1, чего можно добиться, нормируя их на максимальные значения:

(1/t_mах)·d(С₂/S_mах)/d(t/t_mах) = -k₂·(С₂/S_mах)·Х_mах·(С₁/Х_mах)·R_mах·(С₃/R_mах)/[(k₀+S_mах·{С₂/S_mах})·k₃·k₅],

где S_mах = C₀₂, поскольку очевидно, что в отсутствие подачи питательного субстрата (а мы помним, что B = 0) его концентрация в ферментере с течением времени может только падать из-за потребления микробами; Х_mах = C₀₁ + k₅·C₀₂, как это следует из (17) при C₂ = 0, т.е. при полном исчерпании субстрата; t_mах = 6 час. и R_mах = 0.14, как это видно из табл. 2 (хотя R можно было и не нормировать, поскольку, будучи долей РНК, R уже находится в нужных нам границах). Для удобства введем новые (нормированные) безразмерные переменные r = C₃/R_mах, s = C₂/S_mах, τ = t/t_mах, х = C₁/Х_mах и два параметра

K₁ = t_mах·k₂·Х_mах·R_mах/(k₃·k₅·S_mах), K₂ = k₀/S_mах (20)

(очевидно, что K₁ и K₂ также являются безразмерными). Заметим, что используемое здесь значение экономического коэффициента мы уже определили выше (в разделе «Исключение времени как независимой переменной»): k₅ ≈ 6.73.

В новых переменных наше уравнение выглядит так:

ds/dτ = -K₁·s·x·r/(K₂ + s),

по сравнению с (13) качественно ничего, разумеется, не изменилось – уравнение осталось нелинейным (с тем же самым типом нелинейности, что и раньше). Но, домножив обе части последнего уравнения на (k₂ + s), получим

(k₂ + s)·ds/dτ = -k₁·s·x·r.

Если бы s, x и r были известными функциями τ, то полученное выражение можно было бы проинтегрировать следующим образом:

$\underset{1}{\overset{s}{\int }}(K_2+s)\text{d}s=-K_1\cdot \underset{0}{\overset{\text{τ}}{\int }}s\left(\text{τ}\right)\cdot x\left(\text{τ}\right)\cdot r\left(\text{τ}\right)\cdot \text{dτ}$;

здесь мы сразу учли, чему равен нижний предел интегрирования левого интеграла: s(0) = C₂(0)/S_mах = 1. Но несмотря на то, что нам не известен аналитический вид функций s(τ), x(τ) и r(τ), мы для заданного τ можем найти значение интеграла в правой части уравнения численно – ведь значения s, x и r для некоторых τ нам известны из эксперимента. Тогда, если взять множество значений τ, для каждого из них подсчитать интеграл и аппроксимировать полученные значения интеграла какой-то аналитической зависимостью, то фактически мы будем иметь некоторую функцию (обозначим ее z₁) от τ. А интеграл в левой части уравнения легко берется аналитически. В результате вместо последнего уравнения мы будем иметь:

K₂·(s-1) + 0.5·(s²-1) = -K₁·z₁(τ).

Акцентируем внимание на последнем уравнении, поскольку именно его мы будем использовать для определения неизвестных констант. Для этого нужно иметь информацию, относящуюся к текущим значениям s, x и r. Эти данные следует линеаризовать в несколько необычных координатах [Korobov, Ochkov, 2009: p. 144]. Нетрудно заметить, что если ввести переменные у = 0.5·(1-s²) и z₂ = s-1, то у можно рассматривать как линейную функцию z₁ и z₂:

у = K₁·z₁ + K₂·z₂ (21)

и, следовательно, для определения параметров K₁, K₂ можно применить процедуру множественной линейной регрессии. Вычислить у и z₂ весьма просто. Но как быть с интегралом

$z_1=\underset{0}{\overset{\text{τ}}{\int }}s\left(\text{τ}\right)\cdot x\left(\text{τ}\right)\cdot r\left(\text{τ}\right)\cdot \text{dτ}$?

Численное интегрирование становится возможным после проведения предварительной интерполяции модифицированных данных [Korobov, Ochkov, 2009: p. 144]. Впрочем, в данном конкретном случае интегрирование без интерполяции почти не отличается от интегрирования после интерполяции, поэтому, чтобы не запутывать читателя, интерполяцию проводить мы не будем. Обратимся к конкретным вычислениям. Значения исходных нормированных переменных τ, s, x, r, а также зависимой переменной у и линеаризующих независимых переменных z₁, z₂ сведены в Табл. 3. Чтобы найти численные значения параметров K₁ и K_2, нужно осуществить линейную регрессию у по z₁, z₂. Это можно сделать, например, в Excel при помощи средства «Регрессия» из надстройки «Пакет анализа» («Анализ данных»)¹⁷. В качестве входных данных следует взять данные из столбцов «у», «z₁» и «z₂» Табл. 3 (при этом у – объясняемая переменная, а z₁ и z₂ – объясняющие) и провести регрессию с условием «Константа-ноль» – в этом режиме регрессия будет проведена именно для формулы (21), тогда как по умолчанию проводилась бы регрессия более общего вида: у = K₁·z₁ + K₂·z₂ + K₃.

 

Таблица 3. Нормирование и линеаризация данных таблицы 2.

Table 3. Normalization and linearization of the data from Table 2.

τ	х	s	r	у = 0.5·(1-s2)	$z_1={\int }_0^{\tau }s\left(\tau \right)·x\left(\tau \right)·r\left(\tau \right)d\tau$	z2 = s - 1
0	0.057	1	0.282	0	0	0
0.083	0.064	1	0.290	0	0.00144	0
0.167	0.070	1	0.406	0	0.00339	0
0.25	0.086	1	0.481	0	0.00631	0
0.333	0.114	0.969	0.632	0.030	0.01095	-0.031
0.417	0.162	0.923	0.760	0.074	0.01859	-0.077
0.5	0.230	0.838	0.904	0.148	0.03057	-0.162
0.583	0.338	0.692	0.999	0.260	0.04756	-0.308
0.667	0.516	0.485	1	0.383	0.06774	-0.515
0.75	0.766	0.277	0.957	0.462	0.08663	-0.723
0.833	0.969	0	0.872	0.500	0.09509	-1
1	1.021	0	0.828	0.500	0.09509	-1

 

Множественная линейная регрессия дает результат, который неопытному глазу может показаться чрезвычайно странным. На Рис. 4 мы привели копию экрана нашего компьютера с результатами линейной регрессии, выполненной в Ехсеl. Удивительно здесь то, что для коэффициента К₂ («Переменная Х 2») получено отрицательное значение (-0.04). Но из (20) очевидно, что должно быть К₂ > 0 (поскольку по биологическому смыслу k₀ > 0 и S_mах > 0). Однако если проанализировать значения в столбцах строки «Переменная Х 2», то все становится кристально ясным. Действительно, число -0.04 получено с огромной погрешностью, а интервал с 95%-й вероятностью, содержащий истинное значение K₂, простирается от ‑0.17 до +0.10, т.е. положительные значения K₂ вполне возможны. Это свидетельствует о плохой обусловленности задачи построения регрессии в виде (21) по экспериментальным данным Табл. 3: они оказываются в некотором смысле «проще», чем такая «сложная» зависимость, как (21). Точнее говоря, «информативность» конкретно этих экспериментальных данных недостаточно велика, чтобы можно было идентифицировать параметры «сложной» модели. Поэтому упростим математическое выражение с тем, чтобы его сложность стала соответствовать информативности данных.

 

Рис. 4. Результат множественной регрессии у по z₁ и z₂ (данные Табл. 3) по формуле (21).

Fig. 4. The result of multiple regression of y on z₁ and z₂ (data from Table 3) using formula (21).

 

Сделать это можно совершенно формально. Раз мы получили, что K₂ значимо не отличается от нуля, то просто в (21) нужно подставить K₂ = 0, в результате чего получим более простое уравнение регрессии: у = K₁·z₁. Если провести линейную регрессию для этого выражения (например, опять в среде Excel) – см. Рис. 5, то получим значение K₁ = 5.30, причем интервал «95%-й вероятности» не слишком велик: от 5.08 до 5.52, т.е. неопределенность K₁ составляет всего лишь ±4%. Используя (20), легко вычислить отношение

k₂/k₃ = К₁·k₅·S_mах/(t_mах·Х_mах·R_mах) = 5.30·6.73·130/(6·928·0.14) ≈ 5.95 1/час. (22)

 

Рис. 5. Линейная регрессия «зависимой переменной» у = 0.5·(1-s²) по «независимой переменной» $z_1=\underset{0}{\overset{\text{τ}}{\int }}s\left(\text{τ}\right)\cdot x\left(\text{τ}\right)\cdot r\left(\text{τ}\right)d\text{τ}$

Fig. 5. Linear regression of the "dependent variable" у = 0.5·(1-s²) against the "independent variable" $z_1=\underset{0}{\overset{\text{τ}}{\int }}s\left(\text{τ}\right)\cdot x\left(\text{τ}\right)\cdot r\left(\text{τ}\right)d\text{τ}$

 

Кстати, теперь, зная это отношение, можно при помощи (19) оценить k₀: k₀ = k₂/(k₃·0.0586) ‑ 90 ≈ 5.95/0.0586 - 90 = 11.5 мг/л. Следовательно, К₂ = k₀/S_mах = 11.5/130 ≈ 0.0885, что как раз оказывается в интервале, найденном нами выше, когда мы проводили множественную линейную регрессию для (21).

Одномерная нелинейная минимизация невязки

Оглянемся на то, что нам удалось сделать в области определения численных значений параметров модели (12)-(14). Всего нужно было идентифицировать 6 параметров: k₀, k₁, k₂, k₃, k₄ и k₅. Мы смогли довольно точно определить k₅ = 6.73 и – с грехом пополам – k₀ = 11.5 мг/л. Кроме того, нам удалось найти 3 соотношения, связывающие все параметры, а именно речь идет о формулах (16), (18) и (22). Если подставить в них найденное значение k₀, то будем иметь систему из 3 уравнений относительно 4 переменных (k₁, k₂, k₃ и k₄). Перепишем ее в следующем виде:

k₃ = 0.14·(k₁+76.5)/88, k₂ = 5.95·k₃, k₄ = [137.5·k₃/(k₁+126) - 0.089]/0.036. (23)

Понятно, что однозначно решить такую систему нельзя. Но если бы мы знали численное значение k₁, то

• первое уравнение немедленно дало бы нам значение k₃;
• подставив полученное k₃ во второе уравнение, мы получили бы значение k₂;
• а подстановка значений k₁ и k₃ в третье уравнение дала бы k₄.

Казалось бы, к чему эти рассуждения, если k₁ мы все равно не знаем. Однако не будем забывать о наличии у нас набора экспериментальных данных для C₁, C₂ и C₃. Очевидно, что «правильное» значение k₁ должно обеспечить прохождение динамических кривых C₁(t), C₂(t) и C₃(t), рассчитанных по модели (12)-(14), вблизи экспериментальных данных. Иначе говоря, нас устроит вовсе не всякое значение k₁, а лишь такое, при котором отклонение расчетных кривых от экспериментальных данных будет минимальным. Как уже говорилось выше (см. раздел «Подходы к решению»), количественно охарактеризовать степень отклонения теоретических кривых от экспериментальных данных можно при помощи функции (2). В формуле (2) рассчитанные значения обозначены функциями φ_i(k,t). В данном случае это будут решения системы дифференциальных уравнений (12)-(14). Например, φ₁(k,t) = C₁(t), φ₂(k,t) = C₂(t) и φ₃(k,t) = C₃(t), где k = k₁. То, как можно решить систему ОДУ численными методами, описано в обширной литературе – см., например, [Gerald, Wheatley, 1994: p. 393-534; Engeln-Mullges, Uhlig, 1996: p. 423-506; Korobov, Ochkov, 2009: p. 84-131; Davis, 2024: p. 462-548] (и отдельно упомянем несколько пособий, рассматривающих численное решение дифференциальных уравнений именно средствами MATLAB: [Chen et al., 1999: §7.1-7.2; Shampine et al., 2003; Glagolev et al., 2018]). Поэтому здесь мы не будем подробно на этом останавливаться, а просто приведем текст программы, решающей систему (12)-(14) при заданном значении k₁ – см. функцию Ugodchikov3 в Приложении 2 (там k₁ обозначено как Ks).

Вызов этой функции с параметром, представляющим собой значение k₁, приводит к выводу на экран графика, содержащего экспериментальные и расчётные данные, а также четырех числовых значений: k₃, k₂, k₄ (в программе они обозначены соответственно Rm, mu, Т) и невязки (HEB). В данном случае невязка рассчитывается по формуле (2), где D(х) = |х|, а в качестве w_i берутся средние значения соответствующих экспериментальных данных (таким образом, w₁ – это среднее всех экспериментально измеренных концентраций биомассы, т.е. среднее значение по 2-му столбцу таблицы 2; w₂ – среднее экспериментально измеренных концентраций субстрата, т.е. среднее по 3-му столбцу таблицы; w₃ – среднее концентраций РНК, т.е. среднее по 4-му столбцу). Например, если мы выполним расчет для k₁ = 18:

>> HEB=Ugodchikov3(18),

то на экране появится следующий результат:

Rm = 0.1503 mu = 0.8945 T = 1.5154 HEB = 3.2270

На Рис. 6 показаны графики, соответствующие двум вариантам расчетов. Из рисунка 6б видно, что при k₁ = 5 мг/л расчетная кривая C₃(t) качественно не соответствует экспериментальным данным. Действительно, в эксперименте мы получаем максимальное значение C₃ при t = 3.5÷4, после чего наблюдается некоторое снижение доли рибосом (по времени соответствующее исчерпанию субстрата в среде). Расчет же дает монотонно возрастающую кривую, причем при исчерпании субстрата вместо снижения наблюдается резкий рост C₃. Немудрено, что численное значение невязки в этом случае будет относительно велико. Картина на рисунке 6а (полученная в результате расчета с k₁ = 19 мг/л) гораздо лучше соответствует экспериментальным данным: понятно, что и невязка в этом случае будет меньше. Проведя множество расчетов с самыми разными значениями k₁, можно построить кривую зависимости невязки от k₁ (см. Рис. 7).

 

Рис. 6. Результаты численного решения системы ОДУ (12)-(14) при k₁ = 19 (a) и k₁ = 5 (б).

Примечание:

Концентрация биомассы (мг/л): о – экспериментальные данные; ^________ – расчет по модели.

Концентрация субстрата (мг/л): х – экспериментальные данные; ^{. . . . . . . . .} – расчет по модели.

Доля РНК (‰): * – экспериментальные данные; ^{__ __ __} – расчет по модели.

Fig. 6. Results of the numerical solution of the ODE system (12)-(14) for k₁=19 (a) and k₁=5 (б).

Note:

Biomass concentration (mg/L): о – experimental data; ^________ – model calculation.

Substrate concentration (mg/L): х – experimental data; ^{. . . . . . . . .} – model calculation.

RNA fraction (‰): * – experimental data; ^{__ __ __} – model calculation.

 

Рис. 7. Подбор численного значения Ks (т.е. k₁) путем минимизации невязки между решением системы (12)-(14) и экспериментальными данными.

Fig. 7. Estimation of the numerical value of Ks (i.е. k₁) by minimizing the residual between the solution of the system (12)-(14) and the experimental data.

 

Как это часто бывает в реальных (а не чисто учебных) задачах, результат опять оказался несколько неожиданным. Существует довольно четкий минимум невязки (НЕВ ≈ 3.227) при k₁ ≈ 18 мг/л. При меньших значениях k₁ невязка быстро возрастает (причина этого была вскрыта только что выше – при малых k₁ расчетная кривая C₃(t) перестает не только количественно, но и качественно соответствовать экспериментальным данным). А вот при увеличении k₁ невязка сначала немного увеличивается, достигая максимума примерно при 78÷81 мг/л, а потом очень медленно снижается. Важно отметить, что при k₁ ≈ 722 мг/л невязка становится такой же, какой она была при k₁ ≈ 18 мг/л, однако при дальнейшем увеличении k₁ она продолжает падать, хотя все более и более медленно. Так что же – наилучшее значение k₁ находится где-то в области очень больших значений (или, может быть, оно вообще бесконечно)?

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

К сожалению, очень медленное изменение невязки (с изменением k₁) наводит на мысль, что задача идентификации k₁ описанным методом в этой области становится плохо обусловленной или даже вообще некорректной. Более строгое рассмотрение вопроса о степени обусловленности может опираться на следующие рассуждения. Каждая экспериментальная точка получена нами с какой-то погрешностью. Поэтому даже если бы наша математическая модель позволяла абсолютно правильно рассчитать значения концентраций биомассы и субстрата, а также доли РНК, которые были в реальности у культуры E. coli в рассматриваемом эксперименте, невязка между таким расчетом и экспериментально измеренными значениями была бы ненулевой, ибо измерения в любом случае отягощены погрешностями. Для относительно простых невязок и законов вероятностных распределений погрешностей можно аналитически получить оценку того, какой была бы невязка в рассматриваемом идеальном случае. А вообще для любых распределений ее можно получить средствами статистического моделирования (см., например, [Brandt, 1999: Ch. 4; Chen et al., 1999: Ch. 6; Glagolev, Filippov, 2011: §1.5.3;¹⁸ Davis, 2024: §8.3]). Покажем на простом примере, как это сделать.

Пусть погрешности измерения распределены по известному закону. Для определенности будем рассматривать нормальные распределения со средними me_Х, me_S, me_R и стандартными отклонениями σ_Х, σ_S, σ_R (соответственно, для независимых переменных C₁, C₂ и C₃). В Ugodchikov3 мы использовали следующую невязку:

$H=\sum _{m=2}^{12}\left(\frac{\left|C_1\right(t_m)-{\varphi }_1(t_m\left)\right|}{w_1}+\frac{\left|C_2\right(t_m)-{\varphi }_2(t_m\left)\right|}{w_2}+\frac{\left|C_3\right(t_m)-{\varphi }_3(t_m\left)\right|}{w_3}\right)$,

где через С₁, С₂ и С₃ обозначены экспериментальные данные (концентрации биомассы, субстрата и доли РНК), а через φ₁, φ₂ и φ₃ – соответствующие расчетные величины, полученные численным интегрированием модели (12)-(14); m – номер строки в Табл. 2; w₁, w₂ и w₃ – средние значения С₁, С₂ и С₃ (т.е. средние значения по соответствующим столбцам Табл. 2). Для простоты сосредоточим наше внимание пока только на концентрации биомассы. В любой момент времени t_m ее измеренное значение С₁(t_m) будет представлять собой сумму истинного значения этой концентрации, которую мы обозначим через χ(t_m), и случайной погрешности ΔХ (удовлетворяющей нормальному распределению со средним me_Х и стандартным отклонением σ_Х). Очевидно, что при абсолютно правильной работе модели

С₁(t_m) - φ₁(t_m) = χ(t_m) + ΔХ - χ(t_m) = ΔХ,

т.е. невязка все равно не будет нулевой из-за погрешностей измерений биомассы. Разумеется, аналогичные рассуждения относятся к любым измерениям, следовательно, С₂(t_m) - φ₂(t_m) = ΔS и С₃(t_m) - φ₃(t_m) = ΔR, где ΔS и ΔR – соответствующие случайные погрешности измерения С₂ и С₃. Поэтому, чтобы узнать, какой будет невязка в идеальном случае – при абсолютно правильной работе модели, надо просто в формуле для невязки заменить С₁(t_m) - φ₁(t_m) на случайные величины¹⁹ ΔХ; С₂(t_m) - φ₂(t_m) – заменить на ΔS, а С₃(t_m) - φ₃(t_m) – на ΔR. При помощи генератора случайных чисел для распределений с заданными параметрами (в данном случае для трех нормальных распределений: me_Х, σ_Х; me_S, σ_S; me_R, σ_R) следует получить конкретные реализации ΔХ, ΔS и ΔR (в данном случае – 11 реализаций для ΔХ, 11 – для ΔS и 11 – ΔR). В Приложении 1 приведена написанная на языке MATLAB функция stat_mod.

Пусть оказалось, что измерения концентрации биомассы характеризовались средней погрешностью 24 мг/л, причем стандартное отклонение составляло 17 мг/л; при измерениях концентрации субстрата средняя погрешность и ее стандартное отклонение составляли соответственно 10 мг/л и 7 мг/л; а при измерениях доли РНК – соответственно 0.01 и 0.007. Если мы хотим методом статистического моделирования определить среднее значение невязки, обусловленное такими погрешностями экспериментальных данных, то вызвать функцию stat_mod можно следующим образом (используем 10 000 случайных реализаций невязки):

>> [HEBmin,HEBavg,HEBmax,HEB]=stat_mod([24 10 0.01],[17 7 0.007],10000,0.9);

здесь последнее число указывает, что рассчитывается такой доверительный интервал, в который среднее значение попадет с 90%-й вероятностью. Понятно, что статистическое моделирование при каждом запуске программы дает немного различающиеся результаты²⁰. Мы получили следующий результат:

HEBmin = 3.2148 HEBavg = 3.2204 HEBmax = 3.2261

Таким образом, если измерения проведены с описанной выше точностью, то невязка, большая чем 3.22, может быть обусловлена погрешностью модели (а к меньшим невязкам вообще нет смысла стремиться, поскольку при этом мы, скорее всего, пытаемся точно описать не столько особенности исследуемой биологической системы, сколько шум).

Однако значение невязки, примерно соответствующее вышеуказанному «пороговому» значению, достигается на Рис. 6 в двух местах: в достаточно «остром» минимуме (k₁ ≈ 18 мг/л) и существенно правее – при k₁ ≈ 731 мг/л. Какое же значение K_S выбрать? К сожалению, этот вопрос относится, скорее, к теории некорректных задач, которая в данной лекции не рассматривается.

СОПОСТАВЛЕНИЕ РЕЗУЛЬТАТОВ РАЗЛИЧНЫХ МЕТОДОВ

Итак, мы рассмотрели целый ряд методов идентификации параметров математических моделей, но не затрагивали вопрос о погрешности полученных оценок. Эта погрешность может быть оценена при помощи как обычных методов математической статистики, так и нетрадиционных методов, с интенсивным использованием ЭВМ, типа уже упоминавшихся выше методов бутстрэпа, Монте-Карло, складного ножа и др. (см., например, [Efron, 1979; Ryan, 1997: p. 116, 126, 272-275, 290-291, 315-316, 430; Brandt, 1999: §9.8, 9.13; Glagolev et al., 2017: р. 59]), которые из-за ограниченности объема статьи мы не будем здесь демонстрировать. К сожалению, вопросам погрешности получаемых параметров в работах по моделированию часто не уделяется должного внимания. Поэтому в заключение мы считаем нужным сказать об этом несколько слов.

Выше в разделах «Метод приближенного линеаризующего преобразования», «Простейший случай» и «Обобщение на неявные функции» тремя различными методами для одних и тех же экспериментальных данных искались такие параметры роста микробов на 4-хлоранилине, как максимальная удельная скорость роста (k₂) и отношение начальной концентрации микроорганизмов к экономическому коэффициенту их роста на данном субстрате (k₃). Соответственно, эти методы дали следующие значения: k₂ ≈ 0.077, 0.063, 0.113; k₃ = 0.024, 0.077, 0.021. Как видим, значения, полученные разными методами, довольно сильно отличаются друг от друга. Возможным объяснением могло бы быть то, что какого-то метода (или каких-то методов) получены неверные значения. Но поскольку мы старались каждый раз применять конкретный метод в той области значений, где он действительно применим, то это объяснение, по-видимому, должно быть отброшено. Остается предположить следующее: погрешность полученных значений достаточно велика – такова, что с учетом этой погрешности они не противоречат друг другу. Действительно, пусть, например, погрешность определения коэффициента k₂ для каждого метода составляет ±29%, тогда k₂ ≈ 0.077 ± 0.22, 0.063 ± 0.18, 0.113 ± 0.33, и если истинное значение k₂ составляло около 0.081, то это удовлетворяет всем этим оценкам.

АВТОМАТИЗАЦИЯ РЕШЕНИЯ ЗАДАЧ МИНИМИЗАЦИИ В MATLAB

Стандартные решатели MATLAB для оптимизационных задач

Выбор стандартных решателей оптимизационных задач зависит от типа решаемой задачи, ее физического смысла, рельефа целевой функции (критерия оптимальности), оптимизирующих переменных и наличия ограничений. В MATLAB реализованы следующие стандартные решатели: fminbnd, polyfit, lsqcurvefit, fminsearch, fmincon, первый из которых может применяться исключительно при решении задач одномерной оптимизации; второй и третий – при решения задач многомерной параметрической идентификации для процессов, зависящих от одной влияющей на их состояние переменной²¹ (при этом в случае решателя polyfit аппроксимирующая функция представляет собой многочлен произвольной степени, которую надо задать для решения); а два последних – для решения задач как одномерной, так и многомерной оптимизации [Gartman, Klushin, 2020: p. 286-288]. К сожалению, в краткой статье у нас нет возможности описать подробно все вышеперечисленные MATLAB-решатели. Мы более или менее подробно рассмотрим применение лишь одного из них (fminbnd), обсудим вопрос о том, зачем в нашей задаче – идентификации параметров – нужна многомерная минимизация, а далее отошлем читателя к богатой литературе.

Каждый из перечисленных решателей может включать в себя один или несколько алгоритмов оптимизации, которые при корректном задании начального приближения значений оптимизирующих переменных позволяют итеративным путем найти искомые оптимальные величины оптимизирующих переменных [Gartman, Klushin, 2020: p. 287]. Большинство вычислительных функций, в т.ч. fminbnd и fminsearch, написаны на языке программирования MATLAB. Они имеют открытый код и расположены в подкаталоге \toolbox\matlab\funfun [Anufriyev et al., 2005: p. 261-262].

Выбор стандартного решателя MATLAB может быть реализован на компьютере путем задания следующих основных параметров:

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

В некоторых случаях (например, при решении задач нелинейного программирования²²) необходимо дополнительно задать:

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

Наконец, можно задать параметры стандартной функции optimset или optimoptions. В частности, с их помощью возможно задание точности вычислений, отличающейся от задаваемой по умолчанию. А для контроля за ходом вычислительного процесса можно затребовать вывод в командное окно MATLAB последовательных итераций со значениями оптимизирующих переменных и целевых функций [Gartman, Klushin, 2020: p. 287].

Решение задач одномерной оптимизации при помощи fminbnd

Основным решателем для задач одномерной оптимизации является fminbnd (хотя могут применяться fminsearch и fmincon, предназначенные для решения задач многомерной оптимизации) [Gartman, Klushin, 2020: p. 288]. Теперь для этого решателя нам нужно задать конкретную целевую функцию.

Ранее мы уже запрограммировали целевую функцию, рассчитывающую невязку при заданном значении k₁ – это Ugodchikov3. В принципе для решателя fminbnd ее можно взять точно в таком виде, как она использовалась нами выше (т.е. точно в таком виде, как в разделе «Приложение 2…»). Однако в строке этой функции, помеченной комментарием из 7 восклицательных знаков, значения Rm, mu и T не только вычисляются, но и выводятся на дисплей. Поскольку на каждой итерации fminbnd будет предполагать в качестве минимального какое-то значение k₁ и заставлять функцию Ugodchikov3 вычислять невязку при этом значении, то на экран ЭВМ будет выведено много в общем-то не слишком интересной для нас информации – значения Rm, mu и T, вычисленные для всех опробованных значений k₁ (оптимальным из которых будет только одно!). Чтобы избежать этого, достаточно переписать указанную строку в Ugodchikov3, отключив вывод на монитор при помощи точек с запятой, т.е. таким образом:

Rm=0.14*(Ks+76.5)/88; mu=5.95*Rm; T=(137.5*Rm/(Ks+126)-0.089)/0.036;

Наконец, если задаться каким-либо интервалом k₁, на котором мы хотим найти минимум невязки (пусть это будет, например, 0.1 ≤ k₁ ≤ 100), то решение задачи оптимизации займет лишь 1 строку (прямо в командном окне MATLAB – безо всякого программирования):

>> k1_opt=fminbnd('Ugodchikov3',0.1,100)

и довольно быстро (в пределах секунды) на экране появится ответ:

k1_opt = 17.3252,

который совершенно понятен безо всяких объяснений: минимальное значение невязки достигается при значении k₁ = 17.3252. Однако приведенный пример представляет собой самый простой способ вызова fminbnd. Существуют и другие, которые мы сейчас проиллюстрируем на том же самом примере.

В выражение для решателя fminbnd может быть включена стандартная функция optimset [Gartman, Klushin, 2020: p. 288], параметры которой описаны в Табл. 4. Но можно поступить иначе – задать параметры в управляющей структуре (которую мы будем, как и в справочной системе MATLAB, называть options, хотя имя может быть произвольным). Тогда перед вызовом вычислительных функций следует, воспользовавшись функцией optimset, сформировать переменную²³ options в соответствии с характером требуемого контроля. В общем случае входные аргументы optimset задаются попарно [Anufriyev et al., 2005: p. 258-259]:

options=optimset(…,'вид_контроля','значение',…).

 

Таблица 4. Параметры optimset [Anufriyev et al., 2005: p. 259-260].

Table 4. Optimset parameters [Anufriyev et al., 2005: p. 259-260].

Вид контроля	Значение	Результат	Примечание
Display	'off'	Информация о процессе вычислений не выводится
	'iter'	Выводится информация о каждом шаге процесса вычислений
	'final'	Выводится информация только о завершении вычислений
	'notify'	Выводится предупреждение, если процесс не сходится	Значение по умолчанию
MaxFunEvals	Целые числа	Максимальное количество вызовов исследуемой функции	Только для fminbnd и fminsearch
MaxIter		Максимальное количество итераций вычиcлительного процесса
TolFun	Вещественные положительные числа	Точность по функции для останова вычислений
TolX		Точность по аргументу функции для останова вычислений

 

Приступим к формированию структуры options на примере минимизации функции одной переменной при помощи fminbnd. Получим сначала информацию о работе алгоритма минимизации на его последнем шаге. Для этого вызовем функцию optimset с одним выходным аргументом options и двумя входными: 'Display' и 'final', а затем укажем options в дополнительном четвертом входном аргументе MATLAB-функции fminbnd (кроме того, установим формат вывода long, поскольку в дальнейшем нам понадобятся все значащие цифры для анализа результата) [Anufriyev et al., 2005: p. 258]. Все это можно сделать прямо в командном окне MATLAB:

>> options=optimset('Display','final');

>> format long, k1_opt=fminbnd('Ugodchikov3',0.1,100,options)

Ответ MATLAB, появляющийся на экране, будет таким:

Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004

k1_opt = 17.32521144373114.

Как видим, кроме координаты точки локального минимума выводится информация об успешном завершении вычислительного процесса и точности, с которой (по умолчанию) найдено решение. Для изменения точности следует заново сформировать структуру options, указав еще одну пару входных аргументов. В следующем примере мы задаем не только опцию 'Display', но и точность 10^-9 по аргументу при помощи параметра TolX [Anufriyev et al., 2005: p. 259]. Обратите внимание, что для 'Display' мы на этот раз выбрали другое значение. Итак, введем в командном окне MATLAB:

>> options=optimset('Display','iter','TolX',1e-9);

>> k1_opt=fminbnd('Ugodchikov3',0.1,100,options).

Получаем (для экономии места приводим начало и конец «распечатки»):

Func-count x f(x) Procedure

 3 38.2584 3.37439 initial

 4 61.8416 3.43065 golden

 5 23.6832 3.29657 golden

 6 14.6752 3.32481 golden

 7 23.5419 3.29519 parabolic

 8 20.2604 3.25705 parabolic

 9 18.127 3.22887 golden

 …              …          …               …

 34 17.3252 3.21645 parabolic

 35 17.3252 3.21645 golden

Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-009

k1_opt = 17.32519999149728.

Значение, полученное с точность 10^-9, можно считать «истинным» по сравнению с тем, которое получено с точностью 10^-4 (по умолчанию fminbnd использует точность 10^-4 [Anufriyev et al., 2005: p. 259]), поэтому погрешность последнего легко вычислить:

17.32521144373114 - 17.32519999149728 = 0.00001145223386.

Следовательно, при первом вычислении точность была ~10^-5, а не 10^-4, как выведено в сообщении. Информация, возвращаемая функцией fminbnd, содержит гарантированную точность, в то время как реальная точность может быть лучше [Anufriyev et al., 2005: p. 259].

Пользователям, имеющим представление о численных методах, могут оказаться полезны сведения о ходе вычислений, выводимые на экран при значении параметра 'Display', равном 'iter' (как мы это задали в последнем примере). Появилась таблица, каждая строка которой соответствует одной итерации и содержит информацию о том, какой раз вызывалась исследуемая функция – «Func-count», текущее приближение – «x» и значение функции от него – «f(x)», а также метод, применяемый на данной итерации, – «Procedure» [Anufriyev et al., 2005: p. 260]. Может показаться, что последние 8 итераций (с 28-й по 35-ю включительно) не имели смысла, ибо ни значение оптимизируемой переменной (17.3252), ни значение минимизируемой функции (3.21645) не менялись. Однако это постоянство только кажущееся. Просто у оптимизируемой переменной мы видим лишь 4 десятичных знака, а точность задали на уровне 10^-9. На итерациях с 27-й по 35-ю оптимизируемая переменная меняется, но меняется в 5-м знаке; потом – в 6-м и т.д.

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

>> options=optimset('TolX',1e-9);

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

>> options=optimset(options,'TolFun',1e-7,'Display','iter').

При этом можно создать несколько структур с различными именами (к примеру, options1 и options2) для разных вариантов управления вычислительным алгоритмом [Anufriyev et al., 2005: p. 261].

Для просмотра всех текущих опций можно вывести структуру options в командное окно или открыть ее в редакторе массивов двойным щелчком по строке с options в окне браузера переменных рабочей среды [Anufriyev et al., 2005: p. 261]. Проиллюстрируем только первую возможность. Если в точности были выполнены все вышеприведенные примеры манипуляций с optimset, то введя в командной строке

>> options,

получим внушительный список параметров, в котором будут заданы только те, которые мы определили выше при помощи optimset (т.е. Display: 'iter', TolFun: 1.0000e-007 и TolX: 1.0000e-009).

Получение значения отдельного параметра производится при помощи функции optimget, входными аргументами которой являются имя структуры и название требуемого параметра [Anufriyev et al., 2005: p. 261]. Например, если мы хотим узнать, какая задана точность функции для останова вычислений, то следует ввести в командной строке

>> optimget(options,'TolFun')

и незамедлительно будет получен ответ:

ans = 1.0000e-007.

Если данному параметру не было присвоено значение при генерации структуры, то возвращается пустой массив [Anufriyev et al., 2005: p. 261]. Например, выше мы не задавали максимальное количество итераций²⁴. Поэтому при обращении к соответствующему параметру

>> optimget(options,'MaxIter')

в качестве ответа будет получен пустой массив:

ans = [ ].

Вызов функции optimset без входных аргументов

>> optimset

приводит к отображению в командном окне названий всего множества параметров вместе с их допустимыми значениями. В фигурных скобках указаны значения, используемые в вычислительном алгоритме по умолчанию [Anufriyev et al., 2005: p. 261].

В заключение данного раздела опишем способы вызова fminbnd с некоторыми дополнительными выходными параметрами. Еще раз решим задачу нашего примера, но обратимся к решателю вот так:

>>[k1_opt,HEB_min,Exitflag]=fminbnd('Ugodchikov3',0.1,100,options).

Результат будет таким:

Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004

k1_opt = 17.3252 HEB_min = 3.2165 Exitflag = 1.

Exitflag – числовая характеристика вычислительного процесса: Exitflag > 0 – решение найдено с требуемой точностью, Exitflag = 0 – достигнуто максимальное число итераций [Gartman, Klushin, 2020: p. 293], Exitflag < 0 – отсутствие сходимости вычислительного процесса. И, наконец, обратившись к fminbnd таким образом:

>>[k1_opt,HEB_min,Exitflag,output]=fminbnd('Ugodchikov3',0.1,100,options),

получим следующий результат:

Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004

Ks_opt = 17.3252 HEB_min = 3.2165 Exitflag = 1

output = iterations: 24

 funcCount: 27

 algorithm: 'golden section search, parabolic interpolation'

 message: [1x112 char],

откуда становится понятно, что структура output содержит в своих четырех полях информацию о процессе оптимизации: в output.algorithm – использованный алгоритм; в output.funcCount и output.iterations – соответственно количество вычислений целевой функции и количество итераций, понадобившиеся для поиска минимума с необходимой точностью; в output.message – сообщение о том, чем закончилось применение решателя (в данном случае это «Optimization terminated: the current…»).

На этом мы завершим описание fminbnd и optimset, но заметим, что в разделе «Mathematics: Function Functions: Minimizing Functions and Finding Zeros» справочной системы MATLAB размещена подробная информация о минимизации функций и, кроме того, приведены ссылки на литературу с описанием вычислительных методов, реализованных в MATLAB [Anufriyev et al., 2005: p. 261]²⁵.

Решение задач многомерной оптимизации при помощи fminsearch

Однако, увлекшись вопросами программной реализации, мы несколько отвлеклись от существа задачи поиска параметров математической модели. Окинем взглядом решенную нами выше задачу микробиологической кинетики еще раз. При помощи различных частных методов мы свели ее к нескольким линейным задачам поиска тех или иных параметров и к нелинейной задаче поиска параметра k₁. Линейные задачи решались при помощи чрезвычайно простой процедуры линейной регрессии (реализованной даже в электронной таблице Excel), а нелинейная – MATLAB-решателем fminbnd. Понятно, что любой коэффициент был получен нами с некоторой погрешностью. Иначе говоря, мы не можем быть уверены, что, например, экономический коэффициент (k₅) равен в точности 6.73 (напомним, что такое значение было получено нами выше в разделе «Исключение времени как независимой переменной»). Но значение k₁ = 17.3252 получено (в разделе «Решение задач одномерной оптимизации при помощи fminbnd») именно в предположении того, что k₅ = 6.73. А если k₅ = 6.8²⁶? Разумеется, нам не составит никакого труда определить k₁, найдя минимум невязки и в случае k₅ = 6.8. Но тут возникает важный вопрос: будет ли минимальное значение невязки таким же, каким оно было для k₁ = 17.3252 (и, соответственно, для k₅ = 6.73)? Это ниоткуда не следует! Таким образом, вполне может оказаться, что при k₅ = 6.8 будет найдено такое значение k₁ (разумеется, опять из условия минимальной невязки), что минимальная невязка будет меньше, чем она была для k₅ = 6.73 и, соответственно, k₁ = 17.3252.

Если читатель думает, что это все какие-то отвлеченные теоретизирования, что такая ситуация может встретиться редко, то он жестоко ошибается. Давайте практически проделаем все то, что описано выше. Прежде всего, еще раз решим нашу нелинейную задачу при k₅ = 6.73, но вызовем решатель fminbnd таким образом, чтобы увидеть значение минимальной невязки при оптимальном значении k₁:

>> [Ks_opt, HEB_min]=fminbnd('Ugodchikov3',0.1,100,options)

Результат не заставит себя ждать:

Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004

Ks_opt = 17.3252 HEB_min = 3.2165.

А теперь изменим в функции Ugodchikov3 (см. функцию в приложении 2) значение Y: вместо Y=6.73 зададим Y=6.8 и вновь запустим процедуру минимизации. Результат будет уже другим:

Ks_opt = 17.3906 HEB_min = 3.1754.

Для более «плохого» (с точки зрения соответствующей задачи линейной регрессии) значения k₅ = 6.8 удалось найти настолько «хорошее» значение k₁ (= 17.3906), что общая невязка оказалась не больше, а меньше! Иначе говоря, параметры k₅ = 6.8 и k₁ = 17.3906 (все остальные параметры остались неизменными) лучше удовлетворяют экспериментальным данным, чем k₅ = 6.73 и k₁ = 17.3252.

Но ведь мы фактически ткнули пальцем в небо – произвольно выбрали значение k₅ = 6.8. Быть может, для k₅ = 6.77 невязка была бы еще меньше? Или, скажем, для k₅ = 6.83… Очевидно, что среди всех (имеющих биологический смысл) значений k₅ найдется такое, при котором невязка будет «самой-самой минимальной». Итак, возникает задача уже не одномерной, а двумерной оптимизации: подобрать такие значения k₁ и k₅, чтобы невязка была минимальна. Как уже говорилось выше (в разделе «Стандартные решатели MATLAB для оптимизационных задач»), алгоритмы, выполняющие такую оптимизацию, имеются. Более того, существование универсальных решателей, способных найти минимум для функции не только двух, но и трех и вообще многих переменных, наводит на следующую мысль: нельзя ли минимизировать невязку, подбирая значения сразу для всех 6 параметров: k₀, k₁, …, k₅?

Разработано множество решателей задач многомерной оптимизации, которые используются в среде MATLAB. Часть из них поставляется со стандартным пакетом, а остальные (и их большинство) приходится дополнительно приобретать у фирм-дистрибьюторов или скачивать из Интернета [Gartman, Klushin, 2020: p. 294]. Кроме того, значительное число оптимизационных программ (как общего назначения, так и уже сопряженных именно с задачей идентификации параметров) опубликовано в литературе на языке MATLAB и на других языках программирования (см., например, [Johnson, 1980: § 5.4; Bunday, 1984; Brandt, 1999: § 9.6, 9.11, 9.14, 10.8, 10.11, 10.13-10.15, 10.19; Mathews, Fink, 1999: § 8.1]). Они различаются используемыми в решателе алгоритмами оптимизации, возможностями варьирования различных алгоритмов и их параметрами, а также применимостью к решению разнотипных задач [Gartman, Klushin, 2020: p. 294]. Здесь мы рассмотрим многомерную оптимизацию только при помощи стандартной функции MATLAB fminsearch.

При использовании решателя fminsearch достаточно задать начальные значения оптимизирующих переменных в виде вектора, компоненты которого будут изменяться в соответствии с алгоритмом оптимизации, используемым в решателе²⁷ [Gartman, Klushin, 2020: p. 295]. Но, разумеется, теперь (в отличие от одномерного случая) и целевая функция должна воспринимать независимую переменную как вектор.

Выше в наших примерах мы везде рассматривали ситуацию, когда функция Ugodchikov3 получает значение k₁ и вычисляет (при этом заданном k₁) значение невязки между экспериментальными данными и модельными кривыми. Напомним, что все другие коэффициенты модели (k₀, k₂, k₃, k₄ и k₅) при помощи разнообразных частных способов идентификации параметров были уже либо вычислены (k₅ = 6.73; k₀ = 11.5), либо связаны с k₁ однозначными зависимостями (23), по которым вычислялись в Ugodchikov3, как только туда поступало значение k₁. Теперь нам нужно переписать функцию Ugodchikov3 таким образом, чтобы на входе она воспринимала не один лишь параметр (k₁), а все 6 коэффициентов модели. Эту новую функцию мы назовем Ugodchikov2.

Собственно говоря, непосредственно на входе Ugodchikov2 ничего менять не нужно, просто теперь входной параметр coef будем рассматривать (и потому – будем задавать!) не как единственное число, а как вектор, содержащий 6 элементов: значения k₀, k₁, …, k₅. Вообще, изменения, преобразующие Ugodchikov3 в Ugodchikov2, очень незначительны: нам всего-то нужно убрать блок расчета коэффициентов модели, т.е. 4 строки, начиная с комментария «% КОЭФФИЦИЕНТЫ МОДЕЛИ». Кроме того, теперь мы предусмотрели, что в процессе поиска минимума решатель может захотеть проверить отрицательные значения коэффициентов, но т.к. они лишены биологического смысла, то делать этого нельзя. Поэтому перед вызовом функции ode23tb, осуществляющей численное интегрирование дифференциальных уравнений, мы добавили строку, обеспечивающую неотрицательность коэффициентов. Наконец, поскольку в дальнейшем предполагается многократно вызывать Ugodchikov2 из какого-либо решателя оптимизационных задач, мы посчитали, что попытки построить график при каждом таком вызове могут тормозить вычислительный процесс, поэтому «закомментарили» строки, обеспечивающие графический вывод (подчеркнем: лишь «закомментарили», но не стерли их, поскольку когда-нибудь построение графика может понадобиться – и тогда достаточно будет «раскомментарить» эти строки обратно). Ну и, разумеется, мы частично изменили комментарии в «шапке» (см. текст программы Ugodchikov2 в приложении 3). Теперь, когда у нас появилась программа, вычисляющая невязку после получения 6 коэффициентов модели, никакого труда не составит подать имя этой программы и начальные значения коэффициентов на вход функции fminsearch, а она уже автоматически подберет такие их значения, которые обеспечат минимум невязки и, следовательно, наилучшее совпадение кривых, рассчитанных по модели, с экспериментальными данными. В качестве начального приближения к значениям коэффициентов будем использовать те, что были получены (при помощи частных способов) ранее:

k₃ = 0.15, k₂ = 0.89, k₄ = 1.5, k₅ = 6.73, k₀ = 11.5, k₁ = 17.3252. (24)

Итак, в командном окне MATLAB вызовем функцию fminsearch следующим образом:

>>opti=optimset('Display','iter');

>>x_opt=fminsearch('Ugodchikov2',[.15 .89 1.5 6.73 11.5 17.3252],opti)

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

Iteration Func-count min f(x) Procedure

 0 1 3.29005

 1 7 3.03938 initial simplex

 2 8 3.03938 reflect

 … … … …

 4 12 2.94759 contract inside

 … … … …

 8 17 2.93738 reflect

 … … … …

 16 31 2.87947 expand

 … … … …

 32 51 2.54836 expand

 … … … …

 64 101 1.91947 reflect

 … … … …

 128 207 1.87345 contract inside

 … … … …

 212 353 1.87261 contract inside

Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-004

x_opt = 0.1973 1.1790 4.2399 6.8769 -1.2888 5.8503.

На первый взгляд, результат может показаться странным: значение k₀ – отрицательное! Однако вспомним, что в теле функции Ugodchikov2 мы принудительно делали все коэффициенты неотрицательными (при помощи строки, помеченной комментарием: «По биологическому смыслу параметры неотрицательны»). Иначе говоря, когда fminsearch генерирует значение k₀ = -1.2888 и передает его в Ugodchikov2, там система дифференциальных уравнений решается все равно для значения k₀ = +1.2888. Во-первых, это ясно видно из текста программы, а во-вторых, читатель может проверить это самым прямым образом, выполнив расчет модельных кривых при помощи команд

>> HEB = Ugodchikov2([0.1973 1.179 4.2399 6.8769 -1.2888 5.8503])

и

>> HEB = Ugodchikov2([0.1973 1.179 4.2399 6.8769 1.2888 5.8503])

(разумеется, предварительно «раскомментарив» графический вывод, т.е. 3 строки, начиная с plot…). В обоих случаях будет получен совершенно одинаковый ответ:

HEB = 1.8732

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

И в заключение данного раздела нам хотелось бы кратко коснуться вопроса, относящегося не столько к технике решения оптимизационной задачи, сколько к сути получаемых результатов. Выше (в разделе «Уровень невязки, обусловленный погрешностью экспериментальных данных») мы показали, что если измерения концентрации биомассы, субстрата и доли РНК характеризовались соответственно такими средними абсолютными погрешностями (± их стандартные отклонения): 24 ± 17 мг/л; 10 ± 7 мг/л; 0.010 ± 0.007, то невязка, меньшая чем примерно 3.2, может быть обусловлена неточностью самих экспериментальных данных. Если отнести эти средние погрешности к средним значениям измеряемых показателей, то получим следующие относительные погрешности: 24·100%/340 ≈ 7%, 10·100%/89 ≈ 11% и 0.010·100%/0.098 ≈ 10%. Такой уровень относительной погрешности представляется достаточно типичным для измерений в биологических системах, особенно если учитывать вариабельность этих систем, а не только погрешность измерительного прибора.

Приблизительно такой же уровень невязки (между расчетными кривыми и экспериментальными данными) обеспечивает математическая модель (12)-(14) с численными значениями коэффициентов (24). Но теперь нам удалось найти другой набор численных значений тех же самых коэффициентов, который обеспечивает невязку ≈1.9. Следовательно, раз мы считали значения (24) обеспечивающими адекватность модели, то уж тем более таковыми следует признать значения из нового набора. Более того, можно вполне обоснованно предположить, что на самом деле мы имеем не два дискретных набора коэффициентов, а бесконечное множество. Действительно, посмотрим, например, на k₃: в (24) мы имели значение k₃ ≈ 0.15 (при этом НЕВ ≈ 3.2), а теперь – k₃ = R_m ≈ 0.20 (НЕВ ≈ 1.9); вероятно, для любого значения k₃ из интервала от 0.15 до 0.20 можно найти такой набор всех остальных параметров (k₀, k₁, k₂, k₄ и k₅), что невязка окажется где-то между 1.9 и 3.2, т.е. модель с таким набором параметров также можно признать адекватной экспериментальным данным. Понятно, что все то же самое можно сказать относительно любого из коэффициентов модели. Но обратите внимание, что величины «интервалов адекватности» для разных параметров существенно различаются. Если в качестве характерных значений этих параметров взять середины интервалов, то мы получим, что относительная погрешность (будем обозначать ее значком «δ») δk₃ ≈ (0.2-0.15)·100%/0.175 ≈ 29%, а, например, δk₄ ≈ (4.2-1.5)·100%/[(4.2+1.5)/2] ≈ 95%. Подчеркнем: модель (12)-(14) и со значением k₄ ≈ 1.5 час, и с k₄ ≈ 4.2 час (и, скорее всего, с любым значением из интервала 1.5¸4.2), в общем, удовлетворяет имеющимся экспериментальным данным. Понятно, что такая низкая точность определения коэффициента модели вряд ли может удовлетворить современного взыскательного исследователя (хотя, конечно, все зависит от поставленной задачи).

Вспомним, что максимальный уровень невязки (≈ 3.2) соответствовал характерной погрешности экспериментальных данных 7-11%. А для погрешности k₄ мы получаем (из этих относительно точных данных!) – 95%!!! Это свидетельствует либо о плохой обусловленности использованного метода решения обратной кинетической задачи, либо о плохой обусловленности самой этой задачи. Но, как уже было сказано выше, понятия хорошей/плохой обусловленности (и корректности/некорректности) составят тему другой лекции.

БЛАГОДАРНОСТИ

Работа выполнена в рамках государственного задания МГУ им. М.В. Ломоносова.

 

ПРИЛОЖЕНИЕ 1: MATLAB-ФУНКЦИЯ stat_mod

APPENDIX 1: MATLAB FUNCTION stat_mod

function [HEBmin,HEBavg,HEBmax,HEB]=stat_mod(me,sigma,N,P)

%**************************************************************************

% СТАТИСТИЧЕСКОЕ МОДЕЛИРОВАНИЕ НЕВЯЗКИ, ПОРОЖДАЕМОЙ *

% НОРМАЛЬНО РАСПРЕДЕЛЕННОЙ ПОГРЕШНОСТЬЮ ЭКСПЕРИМЕНТАЛЬНЫХ ДАННЫХ *

%**************************************************************************

% ВЫЗОВ: HEB=stat_mod(me,sigma,N,P) *

% *

% ВХОДНЫЕ ПАРАМЕТРЫ: *

% Фактически, экспериментальные данные также являются входным параметром, *

% но ради удобства они задаются непосредственно в теле данной функции *

% (в матрице Y, где каждый столбец соответствует данным одного типа, *

% например, если в разные моменты времени измерялись концентрации био-*

% массы, субстрата и рибосом, то в 1-м столбце матрицы можно размес- *

% тить данные по биомассе, во 2-м - по субстрату и в 3-м - по рибо- *

% сомам; тогда каждая строка матрицы Y будет соответствовать конкрет- *

% ному времени измерений); *

% me, sigma - векторы-строки, соответственно, средних значений абсолютных *

% погрешностей и стандартных отклонений (каждому столбцу матрицы Y со-*

% ответствует свой элемент вектора me и свой элемент вектора sigma); *

% N - количество реализаций вычислений невязки; *

% P - доверительная вероятность для расчета доверительного интервала, со- *

% держащего среднее значение невязки. %

% *

% ВЫХОДНЫЕ ПАРАМЕТРЫ: *

% HEBmin,HEBmax - соответственно, левая и правая границы доверительного *

% интервала, содержащего среднее арифметическое значение невязки, *

% обусловленной погрешностью экспериментальных данных; *

% HEB - вектор, содержащий все значения невязок, полученных в процессе *

% сатистического моделирования; *

% HEBavg - среднее арифметическое значение, вычисленное по элементам HEB. *

%**************************************************************************

% НЕОБХОДИМАЯ ВНЕШНЯЯ ФУНКЦИЯ (из Statistics Toolbox): *

% tinv - Inverse of Student's T cumulative distribution function. *

%**************************************************************************

% ПРИМЕЧАНИЕ: Вычисление доверительного интервала производится в предпо- *

% ложении нормального распределения значений невязки. Поскольку это пред- *

% положение в реальности может не выполняться, к величине доверительного *

% интервала следует относиться с осторожностью. *

%**************************************************************************

% ЛИТЕРАТУРА: Румшиский Л.З. 1971. Математическая обработка результатов *

% эксперимента. М.: Наука. С. 29. *

%**************************************************************************

% ПРОГРАММИСТ: Глаголев М.В. 16/02/2020 *

%**************************************************************************

Y =[59 130 0.041

 65 130 0.057

 80 130 0.068

 106 126 0.089

 150 120 0.107

 213 109 0.127

 314 90 0.140

 479 63 0.140

 711 36 0.134

 899 0 0.122

 947 0 0.116];

[m,n] = size(Y); % Количество строк (m) и столбцов (n) матрицы Ye

YAvg = mean(Y); % Характерные значения переменных

for k=1:N

 y=randn(m,n); % Генератор нормально распределенных случайных чисел

 HEB(k)=0;

 for i=1:m for j=1:n HEB(k)=HEB(k)+abs(me(j)+sigma(j)*y(i,j))/YAvg(j);

 end %for j...

 end %for i...

end %for k...

% Расчет доверительного интервала в предположении нормального распределения

HEBavg=mean(HEB); x=tinv(1-0.5*(1-P),N-1)*std(HEB)/N^0.5;

HEBmin=HEBavg-x; HEBmax=HEBavg+x;

 

ПРИЛОЖЕНИЕ 2: MATLAB-ФУНКЦИЯ Ugodchikov3

APPENDIX 2: MATLAB FUNCTION Ugodchikov3

function HEB=Ugodchikov3(coef)

%**************************************************************************

% ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ *

% СТРУКТУРИРОВАННОЙ МОДЕЛИ ПЕРИОДИЧЕСКОЙ КУЛЬТУРЫ МИКРООРГАНИЗМОВ: *

% dX/dt = m*S*X*R/[(Kr+S)*Rm]; *

% dS/dt = -m*S*X*R/[(Kr+S)*Rm*Y]; *

% dR/dt = [Rm*(Kr+S)/(Ks+S) - R]/T; *

%и расчет невязки между интегральными кривыми и экспериментальными данными*

%**************************************************************************

% ВЫЗОВ: HEB=Ugodchikov3(coef) *

% *

% ВХОДНОЙ ПАРАМЕТР: coef - коэффициент Ks в вышеприведенной системе урав- *

% нений (другие коэффициенты определены из экспериментальных данных в *

% блоке КОЭФФИЦИЕНТЫ МОДЕЛИ). *

% *

% ВЫХОДНОЙ ПАРАМЕТР: HEB - невязка между интегральными кривыми и экспери- *

% ментальными данными, представляющая собой сумму нормированных квадратов *

% отклонений (нормирование осуществляется путем деления квадратов отклоне-*

% ний на среднее значение соответствующей переменной, вычисленное по экс- *

% периментальным данным). *

% *

% ПРОГРАММИСТ: Глаголев М.В. 10/02/2019. (Версия 1.0) *

%**************************************************************************

% ПРИМЕР: при вызове *

% *

% HEB=Ugodchikov3(5) *

% *

% должен быть выдан результат: Rm=0.1297 mu=0.7715 T=1.3081 HEB=4.5821 *

% и выведен график, содержащий экспериментальные и расчетные данные (на *

% графике X и S приведены в [мг/л], а R – в [мгРНК/гБиомассы]). *

%**************************************************************************

% ЭКСПЕРИМЕНТАЛЬНЫЕ ДАННЫЕ:

% X[мг/л] S[мг/л] R[гРНК/гБиомассы] t(Абсцисса, [час])

Ye=[53 130 0.040 0

 59 130 0.041 0.5

 65 130 0.057 1

 80 130 0.068 1.5

 106 126 0.089 2

 150 120 0.107 2.5

 213 109 0.127 3

 314 90 0.140 3.5

 479 63 0.140 4

 711 36 0.134 4.5

 899 0 0.122 5

 947 0 0.116 6];

[m,n] = size(Ye); % Количество строк (m) и столбцов (n) матрицы Ye

time = Ye(:,n)'; % Моменты времени вывода решения

% КОЭФФИЦИЕНТЫ МОДЕЛИ

Ks=coef; Y=6.73; Kr=11.5;

Rm=0.14*(Ks+76.5)/88, mu=5.95*Rm, T=(137.5*Rm/(Ks+126)-0.089)/0.036 % !!!!!!!

cAbs(1)=Rm; cAbs(2)=mu; cAbs(3)=T; cAbs(4)=Y; cAbs(5)=Kr; cAbs(6)=Ks;

Y0 = [Ye(1,1) % Начальные значения: концентрации биомассы

 Ye(1,2) % лимитирующего субстрата

 Ye(1,3)]; % PHK

YeAvg = mean(Ye); % Характерные значения переменных

[T,Y] = ode23tb(@Monod,time,Y0,[],cAbs); % Обращаемся к решателю ОДУ

% Вычисление невязки между экспериментальными данными и расчетом по модели

HEB=0; for i=2:m for j=1:n-1 HEB=HEB+abs(Y(i,j)-Ye(i,j))/YeAvg(j); end, end

plot(T,Y(:,1),'k',T,Y(:,2),':k',T,1000*Y(:,3),'--k',T,...

 Ye(:,1),'ok',T,Ye(:,2),'xk',T,1000*Ye(:,3),'*k');

xlabel('Time, h'), ylabel('Concentration') %Оси координат

%**************************************************************************

% ФУНКЦИЯ, ЗАДАЮЩАЯ ПРАВУЮ ЧАСТЬ СИСТЕМЫ УРАВНЕНИЙ *

%**************************************************************************

function dy=Monod(t,z,cAbs)

Rm=cAbs(1); mu=cAbs(2); T=cAbs(3); Y=cAbs(4); Kr=cAbs(5); Ks=cAbs(6);

X=z(1); S=z(2); R=z(3); % Обозначаем переменные для удобства

dy=[mu*S*X*R/(Kr+S)/Rm % Уравнения динамики концентраций: биомассы

 -mu*S*X*R/(Kr+S)/Rm/Y % лимитирующего субстрата

 (Rm*(Kr+S)/(Ks+S) - R)/T]; % PHK

 

ПРИЛОЖЕНИЕ 3: MATLAB-ФУНКЦИЯ Ugodchikov2

APPENDIX 3: MATLAB FUNCTION Ugodchikov2

function HEB=Ugodchikov2(coef)

%**************************************************************************

% ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ *

% СТРУКТУРИРОВАННОЙ МОДЕЛИ ПЕРИОДИЧЕСКОЙ КУЛЬТУРЫ МИКРООРГАНИЗМОВ: *

% dX/dt = m*S*X*R/[(Kr+S)*Rm]; *

% dS/dt = -m*S*X*R/[(Kr+S)*Rm*Y]; *

% dR/dt = [Rm*(Kr+S)/(Ks+S) - R]/T; *

%и расчет невязки между интегральными кривыми и экспериментальными данными*

%**************************************************************************

% ВЫЗОВ: HEB=Ugodchikov2(coef) *

% *

% ВХОДНОЙ ПАРАМЕТР: coef - вектор коэффициентов [Rm m T Y Kr Ks] в выше- *

% приведенной системе уравнений. *

% *

% ВЫХОДНОЙ ПАРАМЕТР: HEB - невязка между интегральными кривыми и экспери- *

% ментальными данными, представляющая собой сумму нормированных квадратов *

% отклонений (нормирование осуществляется путем деления квадратов отклоне-*

% ний на среднее значение соответствующей переменной, вычисленное по экс- *

% периментальным данным). *

% *

% ПРОГРАММИСТ: Глаголев М.В. 06/04/2020. (Версия 1.0) *

%**************************************************************************

% ПРИМЕР: при вызове *

% *

% HEB=Ugodchikov3([0.1297 0.7715 1.3081 6.73 11.5 5]) *

% *

% должен быть выдан результат: HEB=4.5848. *

%**************************************************************************

% ЭКСПЕРИМЕНТАЛЬНЫЕ ДАННЫЕ:

% X[мг/л] S[мг/л] R[гРНК/гБиомассы] t(Абсцисса, [час])

Ye=[53 130 0.040 0

 59 130 0.041 0.5

 65 130 0.057 1

 80 130 0.068 1.5

 106 126 0.089 2

 150 120 0.107 2.5

 213 109 0.127 3

 314 90 0.140 3.5

 479 63 0.140 4

 711 36 0.134 4.5

 899 0 0.122 5

 947 0 0.116 6];

[m,n] = size(Ye); % Количество строк (m) и столбцов (n) матрицы Ye

time = Ye(:,n)'; % Моменты времени вывода решения

Y0 = [Ye(1,1) % Начальные значения: концентрации биомассы

 Ye(1,2) % лимитирующего субстрата

 Ye(1,3)]; % PHK

YeAvg = mean(Ye); % Характерные значения переменных

cAbs = abs(coef); % По биологическому смыслу параметры неотрицательны

[T,Y] = ode23tb(@Monod,time,Y0,[],cAbs); % Обращаемся к решателю ОДУ

% Вычисление невязки между экспериментальными данными и расчетом по модели

HEB=0; for i=2:m for j=1:n-1 HEB=HEB+abs(Y(i,j)-Ye(i,j))/YeAvg(j); end, end

%plot(T,Y(:,1),'k',T,Y(:,2),':k',T,1000*Y(:,3),'--k',T,...

 % Ye(:,1),'ok',T,Ye(:,2),'xk',T,1000*Ye(:,3),'*k');

%xlabel('Time, h'), ylabel('Concentration') %Оси координат

%**************************************************************************

% ФУНКЦИЯ, ЗАДАЮЩАЯ ПРАВУЮ ЧАСТЬ СИСТЕМЫ УРАВНЕНИЙ *

%**************************************************************************

function dy=Monod(t,z,cAbs)

Rm=cAbs(1); mu=cAbs(2); T=cAbs(3); Y=cAbs(4); Kr=cAbs(5); Ks=cAbs(6);

X=z(1); S=z(2); R=z(3); % Обозначаем переменные для удобства

dy=[mu*S*X*R/(Kr+S)/Rm % Уравнения динамики концентраций: биомассы

 -mu*S*X*R/(Kr+S)/Rm/Y % лимитирующего субстрата

 (Rm*(Kr+S)/(Ks+S) - R)/T]; % PHK

 

¹ Цитируется по [Ivanitskiy, 2001: p. 108]).

² Поскольку в дальнейшем будет важна зависимость φ_i не столько от k, сколько от t, мы иногда будем опускать k, подразумевая, конечно, что φ_i(t) = φ_i(k,t).

³ Когда веса для всех экcпериментальных точек i-го компонента одинаковы, вполне возможно использовать лишь индекс i; а когда, напротив, веса различаются для экспериментальных точек i-го компонента, но не различаются между компонентами (или вообще существует только один компонент), – лишь индекс m.

⁴ Кроме того, если окажется, что С_i(t_m) = 0 для каких-то i и m, то мы вообще «сядем в лужу», поскольку деление на 0 не определено.

⁵ В литературе это распределение («double exponential distribution») встречается также под множеством других названий: «билатеральное экспоненциальное распределение», «двойное показательное распределение», «двустороннее показательное распределение», «показательное распределение с двумя хвостами», «1-й закон Лапласа», «1-й закон ошибок Пуассона», «1-й закон распределения ошибок Лапласа» [Johnson et al., 1994: § 19. 10; Johnson et al., 1995: § 24.1].

⁶ Согласно определению, «критерий» (от греч. kritérion – средство для суждения) – признак, на основании которого производится оценка, определение или классификация чего-либо; мерило оценки [Prokhorov, 1983, p. 654]. В узком смысле критерий должен включать в себя некоторое условие (например, в случае классификации оно указывает, почему рассматриваемый объект относится или не относится к данному классу). Следовательно, просто какое-то число критерием быть не может. Но, видимо, иногда «критерий» понимается более широко – как некоторый показатель (например, число), который следует подвергнуть сравнению с заданным показателем (или заданными показателями) такого же типа, чтобы произвести оценку изучаемого объекта или классифицировать его. Также отметим, что все свои «критерии» Гартман и Клушин [Gartman, Klushin, 2020: p. 285] записывают только для частного случая К = 1.

⁷ «Квадратура» имеет несколько значений (см., например, [Stepanov, 1958: p. 8; Matveyev, 1963: р. 8; Prokhorov, 1983: р. 562; Prokhorov, 1988: p. 264]), в частности «квадратурой» называют операцию взятия неопределенного интеграла [Stepanov, 1958: p. 8]. Видимо, именно в этом смысле понимают данный термин Коробов и Очков [Korobov, Ochkov, 2009: p. 133], когда говорят, что если система дифференциальных уравнений (математическая модель) интегрируется в квадратурах, то «функции φ_i(t)… задаются конкретными аналитическими выражениями, с помощью которых можно осуществить аппроксимацию опытного набора кинетических данных». Напомним, что аналитическое выражение, формула – это выражение, определяющее совокупность действий, которые нужно проделать в определенном порядке над значением аргумента и константами, чтобы получить значение функции [Prokhorov, 1988: p. 71-72]. Поэтому, например, выражение $y={(t^2+1)}^{\frac{1}{3}}$ является аналитическим, а $y=\int {(t^2+1)}^{\frac{1}{3}}dt$ – не является. Но, вообще говоря, в теории дифференциальных уравнений принято более широкое толкование термина «интегрируемость в квадратурах».

Дифференциальное уравнение считается проинтегрированным в квадратурах, если его общее решение получено в явной или неявной форме, которая может содержать еще не взятые интегралы от известных функций [Myshkis, 1964: р. 392]. Следовательно, тип «уравнения, интегрируемые в квадратурах» включает в себя, в частности, уравнения с разделенными переменными: X(t)dt + Y(y)dy = 0. Равенство ∫X(t)dt + ∫Y(y)dy = const является общим решением таких уравнений [Matveyev, 1963: р. 32-65]. Но, например, $\int {(t^2+1)}^{\frac{1}{3}}dt$ и ∫sin^½(y)dy не выражаются через элементарные функции [Myshkis, 1964: р. 322]. Поэтому если X(t) = ${(t^2+1)}^{\frac{1}{3}}$ и Y(y) = sin^½(y), то, записав $\int {(t^2+1)}^{\frac{1}{3}}dt$ + ∫sin^½(y)dy = const, мы, конечно, показали, что соответствующее дифференциальное уравнение ${(t^2+1)}^{\frac{1}{3}}dt$ + sin^½(y)dy = 0 в общепринятом смысле интегрируется в квадратурах, однако, аналитического выражения y(t), дающего зависимость свойства системы (y) от времени (t), мы не нашли, и найти его невозможно, а потому в смысле Коробова-Очкова оно в квадратурах не интегрируется.

⁸ Другие англоязычные исследователи называют их ‘transformably linear models’ [Ryan, 1997, p. 72].

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

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

¹¹ Например, Excel (Add-In “Solver” – надстройка «Поиск решения» в локализованной русскоязычной версии) [Davis, 2024: p. 201-205]; Maple (команда fsolve) [Govorukhin, Tsibulin, 2001: p. 97-98], MATLAB (Optimization Toolbox: функция fsolve) [Anufriyev et al., 2005: p. 733-735, 742-746; Gartman, Klushin, 2020: p. 203-209], библиотека NAG (функция c05nbf) [Glagolev, Smagin, 2005: p. 38-51] и др.

¹² Действительно, ведь экспериментальные данные без них хорошо описались простой моделью (6), которая не содержала параметра k₄. Следовательно, информация, необходимая для определения k_4, если где-то и содержится, то именно в правых точках.

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

¹⁴ Конечно, эти параметры имеют четкий биологический смысл, но он нам сейчас не важен.

¹⁵ Действительно, в ПраЧ исходных уравнений его и так не было (поскольку уравнения – автономные), следовательно, не будет его и в ПраЧ нового уравнения, являющегося отношением ПраЧ исходных уравнений. А в левые части исходных уравнений время входило в виде дифференциала dt, но при делении одного уравнения на другое этот дифференциал сократится.

¹⁶ На всякий случай напомним, что В = 0.

¹⁷ К сожалению, в разных версиях Excel доступ к средствам регрессии осуществляется по-разному, но мы надеемся, что читатель умеет выполнять множественную линейную регрессию в той версии среды Excel, которая у него имеется.

¹⁸ К сожалению, эти авторы допустили досадную описку, могущую запутать неопытного читателя. Довольно подробно описывая статистическое моделирование методом Монте-Карло (именно из-за такого подробного, но емкого и четкого описания мы все-таки решились рекомендовать эту книгу в числе прочих), горе-авторы почему-то называют его «методом бутстрэпа». Метод с таким названием действительно существует, но это вовсе не Монте-Карло, а совершенно другой метод (по своей идеологии, пожалуй, более близкий к методу «складного ножа»).

¹⁹ Мы пишем во множественном числе, потому что на случайную величину заменяется разность в каждый момент времени, когда проводились измерения: С₁(t₂) - φ₁(t₂) = ΔХ₂, С₁(t₃) - φ₁(t₃) = ΔХ₃, …, С₁(t₁₂) ‑ φ₁(t₁₂) = ΔХ₁₂.

²⁰ Если не принять специальных мер, которые мы здесь не будем описывать.

²¹ Т.е. параметрической идентификации, аппроксимирующей функции с одной независимой переменной t вида С = φ(t) [Gartman, Klushin, 2020: p. 288]. Однако эту независимую переменную не стоит путать с оптимизирующими переменными, за счет изменения которых осуществляется процесс минимизации невязки. Разумеется, раз речь идет о решении задач многомерной параметрической идентификации, то оптимизирующих переменных может быть более одной. Выше все такие переменные (k₀, k₁, …, ) мы формально рассматривали в качестве элементов некоторого вектора k, и тогда более развернуто можно записать С = φ(k, t). Например, если изучаемый процесс зависит от одной влияющей на его состояние переменной t, причем предполагается, что зависимость эта описывается полиномом 2-го порядка, то С = k₀ + k₁·t + k₂·t². В данном случае путем варьирования элементов вектора k = {k₀, k₁, k₂} необходимо добиться соответствия между значениями С, измеренными в эксперименте при разных значениях t, и значениями вышеприведенного полинома при тех же значениях t.

²²Раздел математики, посвященный теории и методам решения задач о нахождении экстремумов функций на множествах, определяемых некоторыми ограничениями (равенствами или неравенствами), называется математическим программированием или, иначе, оптимальным программированием (такое «программирование» следует отличать от программирования на ЭВМ). Если изучаемая функция линейна и ищется на множестве, заданном линейными равенствами и неравенствами, то соответствующий раздел математического программирования называется линейным программированием [Prokhorov, 1983: p. 769]. В противном случае имеет место задача нелинейного программирования [Bazaraa, Shetty, 1979: p. 7].

²³ Переменная options на самом деле является структурой. В данном случае умение работать со структурами не требуется, можно просто следовать приведенным правилам заполнения, просмотра и использования управляющей структуры options [Anufriyev et al., 2005: p. 258]. Любознательные читатели могут подробнее освоить работу со структурами данных по многочисленной литературе (см., например, [Govorukhin, Tsibulin, 2001: p. 315-316; Anufriyev et al., 2005: p. 403-421; D’yakonov, 2008: p. 264-268]).

²⁴ Ограничивать количество вызовов функций и число итераций имеет смысл, если есть опасение, что получить решение не удастся из-за расхождения вычислительного процесса. В некоторых случаях приходится уменьшать точность, например, если вычисление исследуемой функции занимает много времени, а требуется получить только несколько первых значащих цифр ответа [Anufriyev et al., 2005: p. 260]. (Здесь же заметим, что, на наш взгляд, используемый Ануфриевым и др. термин «расхождение» применяется в отношении вычислительных процессов крайне редко; в литературе читатель гораздо чаще встретит «расходимость» или «отсутствие сходимости» вычислительного процесса).

²⁵ Для некоторых задач может оказаться полезной дополнительная возможность функции fminbnd (и, кстати, fminsearch) поддерживать работу с функциями, зависящими от параметров. Эта возможность описана, например, в [Anufriyev et al., 2005: p. 264-265].

²⁶ Предполагаем, что это значение также удовлетворяет экспериментальным данным с учетом их погрешности.

²⁷ В нем реализован модифицированный алгоритм Нелдера – Мида [Gartman, Klushin, 2020: p. 296].

About the authors

M. V. Glagolev

Lomonosov Moscow State University; Institute of Forest Science, Russian Academy of Sciences; Yugra State University

Author for correspondence.
Email: m_glagolev@mail.ru
Russian Federation, Moscow; Uspenskoe, Moscow Region; Khanty-Mansiysk

A. F. Sabrekov

Yugra State University

Email: m_glagolev@mail.ru
Russian Federation, Khanty-Mansiysk

D. V. Il’yasov

Yugra State University

Email: m_glagolev@mail.ru
Russian Federation, Khanty-Mansiysk

References

Supplementary files