CC BY
129
УДК 004.9
atrofimovalist. ru
Национальный исследовательский ядерный университет "МИФИ” 2Московский государственный технический университет им. Н.Э. Баумана
Введение
Сахарный диабет - самое распространенное эндокринное заболевание во всем мире. По данным Всемирной организации здравоохранения диабет занимает четвертое место среди причин преждевременной смерти и согласно прогнозам в следующие 10 лет количество смертей от диабета возрастет более чем на 50%, если не будут приняты неотложные меры.
На сегодняшний день известно множество способов лечения сахарного диабета. Лечение инсулином (инсулинотерапия) является жизненно необходимым лицам с диабетом 1-го типа и может применяться в ряде случаев для лиц с диабетом 2-го типа. Целью инсулинотерапии является контроль уровня глюкозы в крови пациента и удержание его в физиологичном интервале путём введения доз инсулина.
Оптимальные дозы вводимого инсулина зависят от многих факторов и индивидуальны для каждого пациента. Подбор этих доз является сложной задачей, с которой могут справиться далеко не все пациенты. Для решения указанной задачи созданы системы непрерывного измерения уровня глюкозы - Continuous Glucose Monitoring Systems (CGM-системы), а также системы непрерывного подкожного введения инсулина (инсулиновые помпы). На основе этих систем разработаны и интенсивно разрабатываются системы автоматического управления уровнем глюкозы в крови пациента.
Необходимым этапом при построении систем управления уровнем глюкозы является построение модели объекта управления - механизма переработки глюкозы клетками организма под действием инсулина. Несмотря на чрезвычайную сложность этого механизма, на сегодняшний день существует ряд упрощённых аналитических моделей, описывающих динамику инсулина и глюкозы в крови пациента (модели Бергмана, Стариса, Никита, Энгельборга и др.). Как правило, эти модели содержат множество эмпирически оцениваемых параметров, индивидуальных для каждого пациента. Детальный обзор существующих моделей уровня глюкозы приведён в [1].
Высокая сложность моделируемого объекта, индивидуальность его характеристик и отсутствие точных математических моделей делает эффективным использование для его моделирования аппарата искусственных нейронных сетей [2,3]. Известно множество архитектур нейронных сетей, применяемых для моделирования динамических объектов, а также моделей на их основе. В [1] приводится описание методик моделирования уровня глюкозы с помощью сетей прямого распространения, рекуррентных нейронных сетей, нечетких нейронных сетей.
В настоящей работе рассматривается подход к моделированию уровня глюкозы с использованием самообучающихся нейронных сетей (SOM, Self-Organizing Maps). Применение самообучающихся сетей для моделирования динамических объектов является относительно новым и малоизученным направлением теории нейронных сетей, о чём свидетельствует небольшое число публикаций в этой области. Однако способность самообучающихся нейронных сетей адаптивно выделять кластеры (паттерны) в пространстве входных переменных позволяет рассматривать их как эффективный инструмент обработки данных, в том числе, динамических.
Статья состоит из трёх разделов. В разделе 1 приводится постановка задачи идентификации динамического объекта, в разделе 2 описывается класс используемых
математических моделей для идентификации и алгоритмы их расчёта, в разделе 3 приводятся результаты идентификации уровня глюкозы с применением самообучающихся нейронных сетей.
1. Постановка задачи
Основная проблема, с которой сталкивается исследователь при решении задачи нейросетевой идентификации, состоит в выборе архитектуры нейронной сети, моделирующей объект. Однозначного решения этой проблемы нет. В каждом конкретном случае выбор класса математических моделей для идентификации определяется интуицией исследователя и результатами сравнительного анализа точности различных моделей.
Известно, что в качестве модели динамического объекта может быть выбрана статическая нейросеть [2]. При этом в сеть некоторым образом должна быть привнесена память. Обычно для этой цели используются линии задержек входного сигнала сети и/или выходного сигнала объекта (TDL, Time-Delayed Line). Обоснованность такого подхода доказывается в рамках теории вложения [4].
Пусть u(t) и y(t), t = 1, T, - наблюдаемые управляющее воздействие и выход объекта соответственно, T - число дискретных отсчётов. Тогда согласно теории вложения выход y(t) модели в момент времени t может быть представлен в следующем виде:
y(t) = F (y(t - ^..^ y(t - ny X u(t X u (t - O^.. u (t - nu)), (1)
где nu, ny - число звеньев в линиях задержки по управлению и по выходу соответственно.
В рассматриваемой схеме идентификации, называемой последовательнопараллельной (рис. 1а), в каждый момент времени на вход модели подаются известные измеренные значения входного и выходного сигналов объекта. Известна другая схема идентификации - параллельная (рис. 16), где в каждый момент времени на вход модели подаются измеренные значения входного сигнала и выходы модели в предыдущие моменты времени:
y(t) = F (y(t - ^..^ y(t - ny X u(tX u(t - O^.. u(t - nu)) . (1')
В обоих случаях задача идентификации сводится к задаче аппроксимации функции F, имеющей для реальных динамических систем, встречающихся в биологии, экономике и
других областях науки и техники, в частности, для системы переработки глюкозы
клетками организма, сложный и априорно неизвестный нелинейный вид.
Рис. 1. Схемы идентификации динамических объектов: а) последовательнопараллельная, б) параллельная.
В настоящей работе рассматривается применение последовательно-параллельной схемы для идентификации уровня глюкозы в крови больных сахарным диабетом 1-го типа.
2. Метод идентификации динамических систем в классе самообучающихся нейронных сетей
Все нейросетевые архитектуры можно разделить на две группы: обучаемые с учителем и самообучающиеся. Задача сетей, обучаемых с учителем (многослойный персептрон, RBF-сети и др.), состоит в установлении соответствия между входными и выходными переменными и аппроксимации тем самым неизвестной зависимости по имеющимся обучающим примерам. Задача самообучающихся сетей, в отличие от сетей, обучаемых с учителей, состоит в выделении схожих в соответствии с некоторым критерием близости групп данных среди имеющихся наблюдений. Популярность самообучающихся сетей связана с их вычислительной простотой и удобством интерпретации результатов самообучения. Как показано в [5], сети SOM демонстрируют результаты не хуже, а в некоторых случаях даже лучше, чем MLP и RBF-сети.
Известно, что основное назначение самообучающихся сетей состоит в решении задач кластеризации данных. В [5,6] предложена методика использования самообучающихся сетей для решения задачи аппроксимации функции F в уравнении (1) посредством кластеризации, называемая VQTAM (Vector-Quantized Temporal Associative Memory). Методика VQTAM для идентификации динамических систем является альтернативой сетям, обучаемым с учителем.
1. Методика VQTAM
Рассмотрим подход, который позволяет сети SOM решать задачу аппроксимации функции F в уравнении (1). Особенностью подхода является то, что в качестве входного
вектора сети Кохонена используется вектор x(t) = ( xm (t), xout (t)) , где
x'n (t) = (y(t - ^..^ y(t - ny X u(t X u(t - ^..^u (t - nu if, xout (t) = y (t),
т.е. желаемое значение функции F подаётся на вход сети. Таким образом, обучающая выборка для сети состоит из набора примеров (xm(t), xout(t)), t = t0,T, где
t0 = max(n„ , ny ) + 1.
Обозначим вектор синаптических коэффициентов /-го нейрона слоя Кохонена через wl(Г) = (н>г"(Г),Н’°и‘(Г)) , / = 1,N, где wгjn(Г) = (Н’г’П(),...,м>™п0)) - вектор
синаптических коэффициентов от входов хт({) размерности п = пу + пи +1, w0iut(^) -
синаптический коэффициент от входа х°ш(£), N - число нейронов в слое. В процессе самообучения для определения нейрона-победителя используются только значения хш(^ и
wT ^):
/ * ^щт Ц х1п ^) - w^n ^ )| 1}, (2)
/=1,N^11 И/
где У - некоторая метрика, как привило, евклидова.
Настройке подвергаются как значения wг" ^), так и w°ut ^), в соответствии с классическим алгоритмом для самообучающихся карт [2]:
Дwi ($) = а^^(/*, /, t) (х(t) - wi ^)), / = 1, N, (3)
где а(0 - параметр скорости самообучения, *, /, t) - функция соседства нейронов на карте, / * - номер нейрона-победителя. В качестве функции соседства ^/*, /, t) может быть выбрана, например, гауссова функция
^і *, і, t) = ехр
( II || 2 Л
Г (t) - г. (tі
2о2 (і)
(4)
где г (і) и Г .(і) - векторы координат нейронов с номерами і(і) и і*(і) на
самообучающейся карте соответственно, параметр <г(і) > 0 определяет радиус функции
соседства на самообучающейся карте на шаге і. В качестве метрики ||-|| на
самообучающейся карте могут быть выбраны евклидова метрика, метрика максимального координатного смещения, суммарного координатного смещения или другие.
После прогона всех примеров обучающей выборки с подстройкой синаптических коэффициентов после каждого примера завершается эпоха самообучения сети. Процесс самообучения, как правило, предполагает проведение нескольких эпох. С началом новой эпохи осуществляется повторный прогон всех примеров обучающей выборки и подстройка синаптических коэффициентов в соответствии с формулой (3). В качестве критерия останова процедуры самообучения может быть предложен останов по достижении заданного числа эпох или заданного значения параметра скорости самообучения.
После того, как сеть обучена, она может быть использована для оценки выхода объекта у(і) в момент времени і:
У(і) = (іX і =1,т, (5)
где индекс і нейрона-победителя определяется по формуле (2).
Таким образом, задача идентификации динамических систем с помощью самообучающихся нейронных сетей состоит в установлении ассоциативной связи между векторами хіп(і) и х°ш(ї) без явного расчёта ошибок идентификации в отличие от использования сетей, обучаемых с учителем, предполагающих аппроксимацию функциональной зависимости Г в уравнении (1) путём минимизации ошибки идентификации.
Схема применения самообучающихся сетей для идентификации динамических объектов приведена на рис. 2, а архитектура используемой сети Кохонена - на рис. 3.
Рис. 2. Схема идентификации на основе самообучающихся моделей.
Рис. 3. Архитектура сети Кохонена, используемая для идентификации (подход VQTAM).
К сложностям подхода относится проблема определения числа нейронов в слое Кохонена. В связи с тем, что подход VQTAM предполагает векторную квантизацию выходных значений, выход у(^) модели может принимать лишь дискретные значения, число которых определяется числом нейронов сети N. Таким образом, для обеспечения точности моделирования может потребоваться довольно большое число нейронов, что нарушает "компактность" модели и увеличивает объём вычислений.
В настоящей работе предложен метод линейной аппроксимации дискретных значений выходов сети Кохонена, обеспечивающий непрерывность выхода самообучающейся модели, что позволяет выбирать меньшее число нейронов сети для обеспечения требуемой точности.
2. Линейная аппроксимация выходов самообучающейся сети.
Алгоритм применяется после завершения самообучения сети Кохонена в соответствии с методикой VQTAM, рассмотренной в предыдущем разделе.
Пусть на вход сети в момент времени I подаётся вектор х1П(() размерности п. Обозначим /1(0,...,4(0 - индексы k ближайших к вектору хш(0 нейронов с точки зрения используемой в (2) метрики, k > п +1. Тогда для векторов wk(0,...,wh (^) синаптических коэффициентов этих нейронов может быть построена аппроксимирующая их
гиперплоскость в ("+1)-мерном пространстве. Выход сети w0Ut (t) определим из условия
принадлежности вектора (xm(t), wout (t)) этой гиперплоскости: wout (t) = a„ (t) + fjaJ (t) xj (t) = (1, x;" (t )T) a(t).
(6)
Расчёт коэффициентов а = ( а0,..., ап) этой гиперплоскости может быть
осуществлён по методу наименьших квадратов [7]: а(1) = ( RT (I) R(I) )* Ег (I) р(1),
(7)
где
R(t) =
Г1 w“(t) ... w£„ (t)^
1 w;;,i(t) ... w""(t)
(8)
vi w;",i(t) ... w;",„(o
P(t) = ( W;0Ut (t),..., WOU (t ))' , а оператор + означает взятие псевдообратной матрицы. Отметим, что при к = " +1 матрица R является квадратной матрицей размерности ("+1)х("+1). В этом случае рассматриваемая гиперплоскость интерполирует центры кластеров ;1(t),_,;k(t), а формула (6) преобразуется в формулу
a(t) = R_1(t) P(t). (9)
Иллюстрация работы алгоритма для двумерного случая (" = 1) при к = 2 приведена на рис. 4.
Рис. 4. Иллюстрация работы линейного аппроксиматора выходов сети Кохонена. Точечной линией обозначена аппроксимируемая функция, пунктирной - выход сети Кохонена, сплошной - результат линейной аппроксимации (к = 2). Знаками "о" обозначены центры кластеров.
Из рисунка видно, что метод позволяет устранить дискретность выходов сети и повысить точность аппроксимации.
3. Результаты экспериментальных исследований
Методика VQTAM была применена для идентификации уровня глюкозы в крови больных сахарным диабетом 1-го типа. В результате четырёхдневного наблюдения над пациентом были получены данные об его уровне глюкозы, уровне вводимого инсулина и уровне потреблённых углеводов, регистрируемые с интервалом 5 мин. в течение каждых суток. Общее число временных отсчетов Т = 1040. Перед построением нейросетевых моделей имеющиеся данные были разделены на обучающую (первые 80% данных, 831 отсчёт) и тестовую (последние 20% данных, 209 отсчётов) выборки.
Вектор управляющих воздействий состоит из двух компонентов и(1) = (м1(!), и2(1)),
где и1(1) и и2(1) - уровень вводимого инсулина и уровень потреблённых углеводов в момент времени I соответственно, наблюдаемый вектор выходов у(1) - уровень глюкозы в крови в момент времени I, I = 1, Т. Графики исходных данных приведены на рис. 5.
1 5 1 О
5
О
I
в
6
4
2
О
6
4
2
О
О 200 400 воо 800 10ОО 1200
Glucose
Insulin
Carbohydrates
Рис. 5. Исходные данные для идентификации: вверху - уровень глюкозы, в центре - уровень вводимого инсулина, внизу - уровень потреблённых углеводов. Для построения нейросетевой модели использованы первые 831 отсчёт, оставшиеся 209 отсчётов использованы для тестирования. По оси абсцисс отложены номера временных отсчётов.
Для оценки качества идентификации уровня глюкозы введены следующие критерии качества.
1. RMSE (Root Mean Squared Error) - среднеквадратическая ошибка:
RMSE =
V
1 T 2
Z( y(t) - y(t))
T -1
1 l0 t =0
2. MAXE (Maximum Absolute Error) - максимальная абсолютная ошибка:
MAXE = max Iy(t) - y(t)|.
t=to ,T
Были построены нейросетевые модели с различными архитектурными параметрами - "u, "y и N и исследована их точность. Значения "u варьировались в пределах [0; 15], "у - в пределах [1; 15] и число нейронов сети N - в пределах [1; 100]. Нейроны были организованы в самообучающуюся карту, подстройка синаптических коэффициентов осуществлялась в соответствии с формулой (3). Зависимость параметра скорости самообучения a(t) от времени выбрана экспоненциальной:
а(0) = а0,
а(!) = уа(!-1), I = 1,2,...
Самообучение большинства рассматриваемых ниже сетей проводилось с параметрами а0 = 0.7, у = 0.974. При указанных параметрах самообучение сети практически останавливается после 100-й эпохи (а(100) = 0.05). В результате экспериментов установлено, что указанное число эпох достаточно для качественной кластеризации при различных значениях пи, пу и N из указанного диапазона изменения. На рис. 6 приводятся графики зависимости ошибок RMSE и МАХЕ на обучающей выборке для 100 эпох (рис. 6а) и 150 эпох (рис. 6б) при различных параметрах скорости самообучения. Из графиков видно, что в обоих случаях (а и б) ошибка устанавливается на некотором уровне ещё до завершения процедуры самообучения. Топология карты выбрана одномерной (рис. 7а), в качестве функции соседства нейронов на карте использована гауссова функция вида (4) с параметром сг(!) = 1, метрика на карте определена как модуль разности между номерами нейронов:
Г 0) - ) = \К!) - 1
1 1 к
и Ый; 4 і і II її
Ґ ' • І 1 1 ■ 1 * 1 и ■ У' І'і
мі г-і1-11
20
40
60
80
100
а)
б)
Рис. 6. Графики зависимости ошибок RMSE (сплошная линия) и МАХЕ (пунктирная линия) на обучающей выборке от номера эпохи самообучения: а) у = 0.974, число эпох 100, а(100) « 0.05, б) у = 0.983, число эпох 150, а(150) « 0.05. Параметры модели: пи = 3, пу = 3, N = 40, начальная скорость самообучения а0 = 0.7.
Проведена серия экспериментов, посвящённых исследованию влияния топологии самообучающейся карты при фиксированном числе нейронов N на точность идентификации. В таблице приведены значения ошибок RMSE и МАХЕ на обучающей и тестовой выборках для различных топологий двумерной карты (рис. 6б). В качестве метрики на карте использовано евклидово расстояние на плоскости. Функция соседства -гауссова вида (4).
о—о—о—о—о—о
о—о—о—о—о—о
с>- -о- -о- -о- -о- -о
о- -о- -о- -о- -о- -о
с>- -о- -о- -о- -о- -о
о—о—о—о—о—о
а) б)
Рис. 7. Организация нейронов в одномерную (а) и двумерную (б) карты Кохонена.
Ошибки RMSE и МАХЕ на обучающей и тестовой выборках для различных топологий самообучающейся карты. Параметры модели Пи = 3, Пу = 3, N = 40.
Размер карты 1*40 2*20 4*10 5*8
RMSE обуч. 0,647 0,633 0,629 0,644
МАХЕ обуч. 2,580 2,373 2,579 2,844
RMSE тест. 0,917 0,917 0,853 0,884
МАХЕ тест. 3,913 3,913 3,385 3,713
Из таблицы видно, что топология карты практически не оказывает влияния на точность идентификации. Подобная картина наблюдается и при других значениях параметров Пи, Пу и N из указанного выше диапазона изменения. В связи с этим в дальнейших экспериментах рассматривается лишь одномерная карта Кохонена.
Основная проблема при построении моделей динамических систем вида (1) по экспериментальным данным состоит в выборе оптимального числа звеньев Пи и Пу в линиях задержек по управлению и по выходу. На рис. 8 приведены графики зависимостей ошибки RMSE на тестовой выборке от параметров Пи и Пу при фиксированном числе нейронов N.
а) б)
Рис. 8. Графики зависимости ошибки RMSE на тестовой выборке от параметров п и Пу: а) N = 50, б) N = 70. Параметр функции соседства нейронов на карте ) = 1,
параметры скорости самообучения у = 0.974, число эпох 100.
Из графиков видно, что с ростом п и Пу точность идентификации падает, оставаясь удовлетворительной при Пи < 5 и Пу < 5, что связано с ухудшением обобщающих способностей самообучающейся сети в пространстве признаков высокой размерности. Наименьшая ошибка идентификации наблюдается при п = 0 и Пу = 1, т.е. при минимальном числе аргументов функции F в выражении (1). Соответствующие этим параметрам результаты идентификации уровня глюкозы на тестовой выборке приведены на рис. 9.
а) б)
Рис. 9. Результаты тестирования нейросетевой модели. Пунктирной линией показаны наблюдаемые значения уровня глюкозы, сплошной - результат работы нейросетевой модели на тестовой выборке. Параметры модели: а) пи = 0 , пу = 1, N = 50, б) пи = 0, пу = 1, N = 70. По оси абсцисс отложены номера временных отсчётов.
На рис. 10 показаны графики зависимости ошибки RMSE на тестовой выборке от числа нейронов N при фиксированных значениях пи и пу.
а) б)
Рис. 10. Графики зависимости ошибки RMSE на тестовой выборке от числа нейронов N. Параметры модели: а) пи = 0, пу = 1, б) пи = 3 , пу = 3.
Из графиков видно, что с ростом числа нейронов примерно до N = 40 ошибки идентификации на тестовой выборке падают, при этом дальнейшее увеличение числа нейронов практически не оказывает влияния на точность. Это обстоятельство говорит о принципиальной невозможности добиться лучшей точности аппроксимации функции F в формуле (1) за счёт более детального квантования выходов самообучающейся сети. Отметим, что с ростом числа нейронов в сети ошибка её работы на обучающих данных лишь уменьшается.
Применение линейного аппроксиматора результатов работы самообучающейся сети в соответствии со схемой, изложенной в разд. 2.2, позволило улучшить точность идентификации. На рис. 11 приведены результаты идентификации уровня глюкозы, полученные с использованием линейного аппроксиматора. В проведённых экспериментах ошибка работы самообучающейся сети на тестовой выборке составила RMSEтест = 0,521, с
использованием линейного аппроксиматора при К = 4 - RMSEтест = 0,411 (рис. 11а), при К = 10 - RMSEmест = 0,426 (рис. 11 б).
а) б)
Рис. 11. Результаты тестирования нейросетевой модели с последующей линейной аппроксимацией. Точечной линией показаны наблюдаемые значения уровня глюкозы, пунктирной - результат работы нейросетевой модели на тестовой выборке, сплошной -результат последующей линейной аппроксимации: а) К = 4, б) К = 10. Параметры модели: пи = 0, пу = 1, N = 50. По оси абсцисс отложены номера временных отсчётов.
Построение самообучающихся нейросетевых моделей уровня глюкозы и проведение экспериментальных исследований выполнялось с помощью разработанного специализированного программного обеспечения IdeпtSOM для идентификации динамических систем в классе самообучающихся нейронных сетей. Программная реализация осуществлена на языке программирования С++ в среде QDevelop с использованием библиотеки Qt и принципов объектно-ориентированного программирования. Визуальное представление результатов осуществлялось средствами системы МАПАБ. На рис. 12 приведены некоторые экранные снимки работы программы.
Рис. 12. Экранные снимки программной системы IdeпtSOM.
Заключение
Основной методический результат, представленный в работе, состоит в определении метода идентификации динамических систем в классе самообучающихся нейронных сетей. С использованием подхода VQTAM, позволяющего применять статические карты Кохонена для решения задач идентификации, была построена
самообучающаяся нейросетевая модель уровня глюкозы в крови больных сахарным диабетом 1 -го типа.
Исследована зависимость точности идентификации от порядка модели и числа нейронов. Установлено, что самообучающуюся модель целесообразно применять лишь при небольших порядках, с ростом порядка модели обобщающие способности нейронной сети ухудшаются. В результате проведённых исследований установлено, что построенная самообучающаяся модель уровня глюкозы по точности не уступает моделям, основанным на многослойных нейронных сетях и сетях радиально-базисных функций.
Преимуществом самообучающихся моделей является их вычислительная простота. В отличие от большинства других нейросетевых моделей и методов их обучения, метод самообучения не требует использования алгоритмов оптимизации, предполагающих трудоёмкий расчёт градиента целевой функции.
Основным недостатком самообучающихся моделей является дискретность их выходных значений. Для решения этой проблемы в работе предложен метод линейной аппроксимации выходов модели, позволивший повысить точность идентификации.
Работа выполнена при поддержке гранта РФФИ №10-08-00816-а.
Литература
1. Гоменюк С.М., Емельянов А.О., Карпенко А.П., Чернецов С.А. Методы
прогнозирования оптимальных доз инсулина для больных сахарным диабетом I типа. Обзор // Наука и образование. №4, 2009 [Электр. журн.
http://technomag.edu.ru/doc/119663. html].
2. Осовский С. Нейронные сети для обработки информации. - М.: Финансы и статистика, 2002.
3. Хайкин С. Нейронные сети: полный курс. - М.: Издательский дом "Вильямс",
2006.
4. Льюнг, Л. Идентификация систем. Теория для пользователя / Л. Льюнг. - М.: Гл. ред. физ.-мат. лит., 1991.
5. Barreto G., Araujo A. Identification and control of dynamical systems using the selforganizing map. IEEE Transactions on Neural Networks. V.15, No.5. P.1244-1259. 2004.
6. Souza L., Barreto G. System identification using local linear models based on the selforganizing map. Learning and Nonlinear Models. V.4. No.2. P.112-123. 2008.
7. Айвазян С. А., Мхитарян B.C. Теория вероятностей и прикладная статистика. М.: ЮНИТИ-ДАНА, 2001.
