Methods of computational optimization for automated insulin therapy control

Full Text

ВВЕДЕНИЕ

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

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

Такой подход является обоснованным, поскольку естественные механизмы управления концентрацией глюкозы в крови даже у пациентов с такими тяжелыми нарушениями, как сахарный диабет 1-го типа, обладают значительной устойчивостью [8]. В то же время данные механизмы обеспечивают большое отклонение установившегося значения от целевого, а применение описанных выше методов внешнего управления хотя и позволяет обеспечивать приемлемое устоявшееся значение, является далеким от оптимального [12].

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

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

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

Одним из наиболее перспективных методов автоматизации управления является комбинация управления по обратной связи с непрерывного датчика глюкозы и управления с прогнозирующими моделями [14]. В то время как в большинстве разрабатываемых систем используется краткосрочное прогнозирование для анализа трендов в динамике концентрации глюкозы в крови и предсказания потенциальных опасных ситуаций [20], мы предложили применять прогнозирующее управления для расчета оптимальных управляющих воздействий [15], а также для коррекции показаний монитора глюкозы, подверженных дрейфу.

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

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

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

1. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ МАТЕМАТИЧЕСКОЙ МОДЕЛИ

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

Разработанная математическая модель содержит 23 ОДУ. Более половины из них являются нелинейными. Наиболее частым случаем нелинейности является наличие членов вида ax/(b + x), соответствующих уравнениям ферментативной кинетики. Наиболее сложное из таких уравнений, описывающее динамику количества запасенной глюкозы в виде гликогена, имеет вид:

$\begin{array}{l}\frac{\text{d}{x}_1\left(t\right)}{\text{d}t}=\frac{k_1{x}_2\left(t\right)x_3\left(t\right)}{k_2+x_3\left(t\right)}-\\ -\frac{k_3{x}_1\left(t\right)x_4\left(t\right)}{k_4+k_5{x}_1\left(t\right)+k_6{x}_5\left(t\right)+k_7{x}_1\left(t\right)x_5\left(t\right)},\end{array}$

где x – переменные;

k – константы.

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

Как правило при работе с моделями систем и объектов управления производится линеаризация моделей, однако линеаризация уравнений разработанной модели изменяет качество моделируемой динамики, что приводит к неадекватной оценке прогнозируемой эффективности того или иного управляющего воздействия на объект управления. Таким образом для решения уравнений модели нельзя применять аналитические и линейные методы решения систем ОДУ, и требуется численное решение.

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

Два уравнения модели содержат вынуждающие функции, описывающие прием пищи и введение инсулина, и имеют вид:

$\frac{\text{d}x\left(t\right)}{\text{d}t}=-kx\left(t\right)+f\left(t\right).$

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

Поскольку для обозначенной задачи нет общепринятого высокоэффективного метода решения, предложено исследовать современные неявные методы и алгоритмы решения ОДУ. В рамках исследования рассматривались обратный метод Эйлера, неявные методы Рунге–Кутты, в т.ч. метод Радау [9], методы Адамса–Мултона, неявные методы Гира (метод обратного дифференцирования, BDF), метод Милна, алгоритм LSODA [10], комбинирующий методы Адамса и BDF.

Критериями эффективности являлась вычислительная устойчивость в представленной задаче, достигаемая точность решения и скорость решения.

Для исследования задан типовой сценарий функционирования объекта управления, включающий внешнее возмущение в виде приема пищи с массой углеводов 90 г с прямоугольным профилем длительностью 20 мин, а также введение двух прямоугольных управляющих воздействий в виде введения инсулина инсулиновой помпой с прямоугольным профилем. Каждое воздействие характеризуется временем начала, длительностью и общей дозой инсулина. Производился перебор по равномерной 6-мерной сетке порядка 10³ комбинаций данных параметров и для каждой комбинации вычислялось значение функции эффективности управления. Эффективность определялась по критериям [18] максимального значения целевой величины выше заданного верхнего порогового значения (соответствует степени гипергликемии), минимального значения целевой величины ниже заданного верхнего порогового значения (соответствует степени гипогликемии), интегрального превышения профилем целевой величины верхнего порогового значения (соответствует образованию гликированного гемоглобина).

Предварительно оптимальные параметры вычислялись в ходе перебора 10⁶ параметров методом равномерного поиска всеми методами по отдельности с допустимой ошибкой 10^–7, результаты для которого совпали для всех рассматриваемых методов. После чего оценивалось время вычисления, отклонение от оптимального значения и гладкость полученной гиперповерхности для каждого метода при заданной допустимой ошибке. Наилучшие результаты показали методы LSODA, Радау и метод BDF. Тем не менее метод BDF не обеспечил гладкость решения [4], что не позволяет применять его для дальнейшего решения задач оптимизации, в то же время метод LSODA, частично использующий метод BDF обеспечивает гладкость и за счет использования метода Адамса на участках, далеких от переключений, обеспечивает лучшую скорость вычисления. Таким образом, для решения дальнейших вычислительных задач применяется метод LSODA с допустимой ошибкой 0,015%, обеспечивающий погрешность вычисления менее 0,01%.

2. ПАРАМЕТРИЧЕСКАЯ ИДЕНТИФИКАЦИЯ

Другой вычислительной задачей является параметрическая идентификация математической модели. Модель содержит 51 параметр. Часть параметров модели представляет собой измеримые величины, к которым относятся масса тела человека, базальные уровни некоторых веществ, входящие в модель как пороговые значения. Некоторые параметры не обладают специфичностью для отдельного человека. Так низкой вариабельность и изменчивостью от человека к человеку характеризуются константы Михаэлиса, входящие в ферментативные уравнения. В то же время большая часть параметров, в первую очередь эмпирических, не может быть прямо измерена, либо не имеет прямого физического смысла. Значения таких параметров необходимо найти численно.

Для решения этой задачи необходимо найти набор параметров, обеспечивающий наилучшую сходимость к некоторому набору экспериментальных данных. Учитывая, что количество подбираемых параметров в наборе составляет 43, данная задача имеет высокую размерность, что требует применения методов оптимизации вычислений. В качестве задачи оптимизации предложено минимизировать функционал:

$\mathrm{\Phi }\left(k_j,\mathrm{...},k_n\right)=\underset{T}{\int }\left|g_p\left(k_j\mathrm{,...},k_n,t\right)-{\tilde{g}}_p\left(t\right)\right|\text{\hspace{0.17em}d}t,$ (1)

где k_i, … , k_n – оптимизируемые параметры;

g̃_p – экспериментальные данные о концентрации глюкозы в плазме;

Т – временной промежуток, на котором представлены экспериментальные данные.

В качестве источника экспериментальных данных использовался симулятор пациента с сахарным диабетом первого типа, разработанный в университетах Вирджинии и Падуи UVA/Padova T1DMS [7]. Симулятор содержит 300 «виртуальных пациентов», разделенных на взрослых подростков и детей, характеризующихся набором скрытых и открытых параметров, а также встроенную математическую модель, на основе которой по заданному сценарию можно строить профили концентрации глюкозы в крови и других измеримых величин. Дальнейшие эксперименты проводились для 10 случайных «взрослых» пациентов.

С использованием симулятора для каждого «виртуального пациента» осуществлялась генерация профилей глюкозы при нескольких сценариях, включающих один прием пищи разной амплитуды и управляющие воздействия инсулином различных амплитуд, вводимым двумя импульсами. Измеримые параметры модели определялись из документации соответствующего «виртуального пациента» симулятора. Часть параметров, относящихся в первую очередь к уравнениям ферментативной кинетики, взята из литературных источников.

Начальное предположение значений большинства параметров выбиралось на основе типовых значений, представленных в литературе. Поскольку разработанная модель в значительной степени основана на модели К. Далла Ман с соавт. [6], типовые значения 23 коэффициентов содержатся в статьях авторов. Еще 9 параметров, к которым относятся скорости деградации и секреции глюкагона, максимальные скорости производства и убыли глюкозы при гликогеногенезе, глюконеогенезе и гликогенолизе, а также параметры секреции глюкагона, были определены из иных литературных данных. Начальные предположение значений остальных 11 коэффициентов, главным образом константы скорости активации и деактивации ферментов, подбирались в ходе численных экспериментов для обеспечения соответствия моделируемой динамики частных процессов литературным данным.

Численное получение функции g_p(k_i, … , k_n, t) для заданного сценария требует решения задачи Коши, для чего требуется массив начальных условий модели. Ввиду того, что большинство переменных модели являются недоступными для экспериментального измерения, а некоторые не имеют аналогов в других уже идентифицированных моделях, необходимо для каждой итерации определять корректные начальные условия.

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

${\overline{x}}_k=\underset{t\to \infty }{\lim }{x}_k\left(t\right).$

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

Таким образом, одним из способов вычисления начальных условий для нахождения значения функционала (1) для заданного набора условий является предварительное вычисление значения функции через некий длительный интервал времени T_b, при котором конечное значение можно считать устоявшимся, а переходной процесс завершенным по какому-либо уровню, например, 0,05.

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

Тогда алгоритм нахождения начальных условий для оценки m-го набора параметров модели для заданного n-го сценария состоит из следующих шагов:

• задание m-го набора проверяемых значений параметров модели;
• задание n-го сценария;
• решение модифицированной системы ОДУ на интервале T_b с заданным грубыми оценками начальных условий;
• принятие в качестве устоявшегося x̄_k = x_k(T_b).

Экспериментальные исследования показали, что T_b порядка 2000 мин достаточно для определения начальных условий, поскольку, за это время переходной процесс при верных значениях искомых коэффициентов закончится.

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

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

• подсчет параметров сценариев;
• подсчет массивов экспериментальных данных, сгенерированных в симуляторе UVA/Padova T1DMS, соединение их в единый массив;
• вычисление начальных условий всех переменных;
• задание начальноого предположения;
• решение системы ОДУ для каждого сценария и соединение результатов решений для концентрации глюкозы в единый массив;
• оценивание точности набора параметров путем оценки значения функционала;
• повторение, пока точность не достигнет целевого значения в соответствии с алгоритмом оптимизации или время вычисления не превысит 240 минут;
• регистрация значения минимума функционала и время его вычисления.

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

Сравнивались такие алгоритмы и метода глобальной оптимизации как метод полного перебора, метод «прыжков по бассейну» [21], метод дифференциальной эволюции, алгоритм разделения прямоугольников DIRECT [13], алгоритм Бройдена–Флетчера–Гольдфарба–Шанно (BFGS) и др. Ввиду высокой размерности задачи неэффективные методы, такие как полный перебор, не позволяют за разумное время приблизиться к целевым результатам. Методы безусловной оптимизации, такие как DIRECT и BFGS, также не позволяют решить представленную многомерную задачу за приемлемое время. Наибольшую скорость решения задачи показал метод «прыжков по бассейну», основанный на методе Монте-Карло, использующий локальную минимизацию в точках «прыжков» и таким образом сильно зависящий от эффективности локального минимизатора.

Среди локальных минимизаторов рассматривались, метод сопряженных градиентов, метод Пауэлла, усеченный метод Ньютона (TNC), метод Нелдера–Мида и др. Метод сопряженных градиентов показал неустойчивую сходимость к ошибочным точкам для представленной в задаче функции, метод Пауэлла показал чрезвычайно низкую скорость нахождения локального оптимума, метод Ньютона в свою очередь показывал в ряде случае ошибку, близкую к бесконечной. Наилучшую скорость при высокой точности показал метод Нелдера–Мида.

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

3. ОПТИМИЗАЦИЯ УПРАВЛЯЮЩИХ ВОЗДЕЙСТВИЙ

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

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

При этом звено непрерывного управления может осуществлять либо программное управление, либо управление на основе ПИД-регулятора [19]. В последнем случае оптимизация может осуществляться стандартными методами оптимизации ПИД-управления и в настоящей работе не рассматривается. Управляемыми параметрами для описанных болюсных воздействий являются амплитуда и время первого воздействия, амплитуда, время начала и длительность второго воздействия.

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

 

Рис. 1. Симуляция типовой динамики концентрации глюкозы в крови при трехкратном приеме пищи

Fig. 1. Simulation of typical dynamics of blood glucose concentration during a single meal

 

Задача оптимизации управления концентрацией глюкозы в крови при наличии внешних возмущений при известном состоянии системы y₀, известных параметрах внешних возмущений С, заключается в определении оптимального набора параметров P, обеспечивающего минимизацию функционала Ф(y₀, C, P).

Множество значений каждого параметра можно ограничить. Так амплитуды первого и второго воздействия целесообразно представить общей дозой инсулина и долей инсулина, приходящейся на первый импульс. Общая доза инсулина ограничивается лечащим врачом, а доля распределения принадлежит отрезку [0, 1]. Время начала воздействий целесообразно отсчитывать от момента фактического начала приема пищи. Первый импульс опасно вводить ранее, чем за 30 минут до начала приема пищи, таким образом множество значений составляет [–30, 0] минут. Второе воздействие вводится как можно раньше после окончания приема пищи, а его длительность нецелесообразно выбирать более 60 минут, поскольку пик отклика объекта управления за это время уже будет близок к наступлению. В связи с ограниченностью множеств значений параметров, задача является задачей условной оптимизации.

В ходе решения задачи оптимизации исследована эффективность различных методов оптимизации. В ходе исследований экспериментально установлена унимодальность минимизируемого функционала, в связи с чем возможно применение методов локальной многомерной оптимизации. Среди таких методов выбраны как одни из наиболее эффективных методов градиентный спуск, метод сопряженных направлений (метод Пауэлла) [16] и метод условной оптимизации с помощью линейной аппроксимации (constrained optimization by linear approximations, COBYLA) [3; 17]. Сравнение проводилось путем численного решения уравнений модели для каждого набора оптимизируемых параметров. Ввиду неизменности параметров модели начальные условия, соответствующие устоявшемуся состоянию, определялись один раз для всего исследования. Через 60 мин после начала моделирования действовало внешнее возмущение (прием пищи) прямоугольной формы длительностью 20 мин, соответствующее 90 г углеводов. Рассчитывались оптимальные параметры управляющего воздействия. В момент начала управляющего воздействия, непрерывное звено управления отключалось на 4 ч. Оценивалось отклонение от предварительно рассчитанного методом равномерного поиска оптимума, время нахождения оптимума, отклонение значений оптимальных параметров.

Наиболее высокую эффективность показал метод COBYLA, обеспечивающий нахождение близких к оптимальным значений за время менее 10 с на персональном компьютере средней мощности. Достижение отклонения менее 0,01% по главному параметру, общему количеству инсулина, занимает порядка 15 с.

4. ОБСУЖДЕНИЕ И ЗАКЛЮЧЕНИЕ

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

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

Как видно из рис. 2, типичный суточный профиль глюкозы пересекает пороговые значения, вблизи которых явные методы решения дифференциальных уравнений дают взрыв погрешности, всего порядка 10 раз. Отсюда следует, что большая часть вычислений проходит в нежесткой области. Этим обосновывается эффективность применения алгоритма LSODA, осуществляющего переключения между одним из наиболее быстрых явных методов Адамса и относительно медленным неявным методом обратного дифференцирования.

 

Рис. 2. Динамика концентрации глюкозы в крови с порогами выработки гликогена и почечной экскреции

Fig. 2. Dynamics of blood glucose concentration with thresholds of glycogen production and renal excretion

 

Алгоритм LSODA позволяет производить расчет модели с ошибкой менее 0,01% на 720 мин за время порядка 0,2 с на процессоре Intel Core i5-4570 при реализации алгоритма в среде PyCharm. Другие испытанные методы для решения аналогичной задачи требуют более 0,5 с.

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

Другой важной проблемой работы с комплексными прогнозирующими моделями является большое количество неизмеримых параметров, значений которых нужно определить для каждого пациента. Параметрическая идентификация должна быть основана на экспериментальных данных о динамике целевой величины с максимально полной информацией о внешних возмущениях и управляющих воздействиях. В целях безопасности сбор таких данных должен производиться пассивно без применения «необученной» системы управления. Предложенный вычислительный метод параметрической идентификации, основанный на глобальном алгоритме «прыжков по бассейну» с локальным минимизатором Нелдера–Мида позволяет решать задачу на указанных выше программно-технических средства за время от 4 до 10 ч при оптимизации по 24 часам входных экспериментальных данных.

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

Третьей проблемой, рассматриваемой в работе, является разработка методов многомерной оптимизации управляющего воздействия. Метод условной оптимизации с помощью линейной аппроксимации COBYLA обеспечивает время оптимизации приемлемое для расчета оптимального воздействия без распараллеливания на невысоких вычислительных мощностях, что принципиально позволяет производить вычисления на смартфоне пользователя без использования серверных мощностей.

Ускорение оптимизации открывает новые возможности, в частности исследования оптимальных временных распределений управляющего воздействия в заданной области в противовес стандартной двухволновой схеме введения болюса. Для поиска первого приближения целесообразно разбить время допустимого введения инсулина (от –30 до +60 минут от момента начала приема пищи) на интервалы и провести поиск оптимального распределения заданной дозы по данным временным отрезкам с учетом предельного объема инсулина, который может быть введен до подтверждения приема пищи. Такой поиск целесообразно провести при возмущениях, происходящих в устоявшемся состоянии, в точке пика отклика на инсулин и в точке второго пика отклика на прием пищи. Если получившееся оптимальное распределение может быть с достаточной точностью аппроксимировано некоторой аналитической кривой, то на второй стадии можно провести оптимизацию с более высокой точностью параметров полученной функции.

About the authors

Kirill V. Pozhar

National Research University of Electronic Technology (MIET)

Author for correspondence.
Email: pozhar@bms.zone
ORCID iD: 0000-0001-9879-0220
SPIN-code: 6609-8070

Cand. Sci. (Eng.), Associate Professor; associate professor, Institute of Biomedical Systems

Russian Federation, Zelenograd, Moscow

Dmitry A. Chuprakov

National Research University of Electronic Technology (MIET)

Email: 89120209984d@gmail.com
ORCID iD: 0009-0002-9384-2049
SPIN-code: 2753-4276

lab assistant, Laboratory of Systems of Artificial Biomedical Regulation, Institute of Biomedical Systems

Russian Federation, Zelenograd, Moscow

References

Supplementary files