Искусственная модель поджелудочной железы

artificialPancreasWithFISTreeControl Модель Simulink® реализует искусственную систему поджелудочной железы.

model = 'artificialPancreasWithFISTreeControl';
load_system(model)

Эта модель содержит следующие подсистемы.

Чтобы просмотреть, как уровень глюкозы крови страдающего от диабета пациента изменяется, когда тело поглощает глюкозу без корректирующего вливания инсулина, симулируйте модель с постоянным нулевым действием управления. Для этого откройте цикл управления установкой closeLoop к 0.

closeLoop = 0;
openLoopOutput = sim(model);
plotAbsorbedAndBloodGlucose(openLoopOutput)

Без корректирующего вливания инсулина терпеливый уровень глюкозы крови увеличивается и остается в гипергликемическом состоянии.

Закройте цикл управления для настройки и симуляции.

closeLoop = 1;

Создайте древовидную структуру контроллера FIS

У контроллера есть следующие входные параметры, как описано в [4].

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

Произвести оптимальную дозу инсулина на основе наблюдаемого binputs, нечеткий контроллер, описанный в [4] экспертные знания использования, чтобы создать один FIS с 75 правилами. Однако создание большой основы правила, использующей экспертные знания, является сложным процессом из-за ручной конструкции каждого нечеткого правила для всех комбинаций входных функций принадлежности (MFS).

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

Этот пример использует инкрементный подход проекта, чтобы объединить входные параметры контроллера с помощью двух объектов Mamdani FIS в инкрементной древовидной структуре. Для получения дополнительной информации о нечетких древовидных структурах смотрите Нечеткие Деревья.

Уровень глюкозы крови и его скорость изменения и способствуют больше действиям управления по сравнению с ускоряющим уровнем, который часто мал и может создать шум в выходе. Поэтому на первом уровне дерева FIS, вы предварительно вычисляете скорость введения инсулина Precalculated_Dose путем объединения эффектов уровня глюкозы крови BG_Level и его скорость изменения BG_Rate. Ускоряющий уровень BG_Accel включен во второй слой дерева FIS.

Создайте FIS (fis1) для первого уровня древовидной структуры. Входные параметры fis1 каждое использование три равномерно распределенных треугольных MFS. Выход fis1 использование пять таких MFS, названных можно следующим образом:

% Specify the maximum dose level.
maxDose = 2;

% Define membership function names for the input variables.
mfNames1 = ["L","M","H"]; % Low, Medium, High
mfNames2 = ["N","Z","P"]; % Negative, Zero, Positive

% Create first FIS.
fis1 = mamfis('Name','fis1','NumInputs',2,'NumOutputs',1, ...
    'NumInputMFs',3,'NumOutputMFs',5);

% Configure input and output variables.
fis1 = updateInput(fis1,1,'BG_Level',[80 120],mfNames1);
fis1 = updateInput(fis1,2,'BG_Rate',[-0.5 0.5],mfNames2);
fis1 = updateOutput(fis1,1,'Precalculated_Dose',[0 maxDose]);

figure
plotfis(fis1)

Создайте FIS (fis2) для второго уровня дерева. Используя fis2, вы генерируете итоговую дозу инсулина путем объединения предрасчетной дозы от первого слоя с эффектом ускоряющего уровня глюкозы крови. В этом случае входные параметры и выход также используют три и пять равномерно распределенных треугольных MFS, соответственно.

% Create second FIS.
fis2 = mamfis('Name','fis2','NumInputs',2,'NumOutputs',1, ...
    'NumInputMFs',3,'NumOutputMFs',5);

% Configure input and output variables.
fis2 = updateInput(fis2,1,'Precalculated_Dose',[0 maxDose],mfNames1);
fis2 = updateInput(fis2,2,'BG_Accel',[-0.005 0.005],mfNames2);
fis2 = updateOutput(fis2,1,'Insulin_Dose',[0 maxDose]);

figure
plotfis(fis2)

Объедините fis1 и fis2 в древовидную структуру FIS.

connection = [fis1.Name + "/" + fis1.Outputs(1).Name ...
    fis2.Name + "/" + fis2.Inputs(1).Name];
fisTInit = fistree([fis1 fis2],connection);

figure
plotfis(fisTInit)

Начальные нечеткие системы создаются с нечеткими правилами по умолчанию, которые не настраиваются, чтобы произвести оптимальные дозы инсулина.

Настройте правила контроллера

Если у вас есть древовидная структура FIS, можно оптимизировать поведение контроллера путем настройки правил и параметров MF объектов FIS компонента. Для этого можно использовать tunefis функция.

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

Чтобы изучить правила, сначала получите настраиваемые настройки от нечетких систем.

[in,out,rule] = getTunableSettings(fisTInit);

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

for rId = 1:numel(rule)
    rule(rId).Antecedent.Free = false;    
end

Создайте набор опции для настраивающего процесса.

options = tunefisOptions;

Используйте настраивающий метод генетического алгоритма по умолчанию для изучения правил. Определите максимальный номер поколений к 3 и используйте численность населения 100.

options.MethodOptions.MaxGenerations = 3;
options.MethodOptions.PopulationSize = 100;

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

В данном примере функция стоимости (costFcn.m) симулирует artificialPancreasWithFISTreeControl модель с помощью кандидата управляет основой. Используя получившуюся симуляцию выход, функция стоимости вычисляет стоимость с помощью следующих шагов.

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

% Specify the nominal and minimum glucose levels.
refLevel = 90;
minLevel = 80;

% Calculate error from the nominal value.
err = glucose - refLevel;

% Specify high error values for the glucose levels below the nominal level.
err(glucose<minLevel) = 100;

% Calculate cost as the root mean square of the error.
errSquare = err.^2;
meanSquare = mean(errSquare);
cost = sqrt(meanSquare);

Настройка является длительным процессом, таким образом, для этого примера, загрузите предварительно настроенное дерево FIS. Чтобы настроить дерево FIS самостоятельно вместо этого, установите runtunefis к true.

runtunefis = false;

% Load pretuned FIS tree data
data = load('fuzzyPancreasExampleData.mat');

minData = MinCostData;
wsVars = ["fisT" "closeLoop"];
minLevel = 80;
refLevel = 90;
if runtunefis
    rng('default')
    fisTRuleTuned = tunefis(fisTInit,rule,...
        @(fis)costFcn(fis,model,minLevel,refLevel,wsVars,minData),options);
else
    fisTRuleTuned = data.fisTRuleTuned;
    minCost = costFcn(fisTRuleTuned,model,minLevel,refLevel,wsVars)
end

minCost = 26.3335

Симулируйте модель с помощью дерева FIS с настроенными основами правила.

fisT = fisTRuleTuned;
ruleTunedOutput = sim(model);

Постройте получившиеся отрегулированные уровни глюкозы и скорость введения инсулина.

plotGlucoseAndInsulin(ruleTunedOutput,...
    'Blood Glucose and Insulin Dosage with Learned Rule Base')

С настроенной основой правила уровень глюкозы теперь обеспечен ниже 160 мг/дл, и это обосновывается близко к 90 мг/дл после третьей еды. Контроллер генерирует краткосрочную высокую дозу инсулина (режим шарика) в каждый раз еды и долгосрочную уменьшаемую дозу инсулина (основной режим) в периоды поста.

Анализируйте и измените основу правила

Чтобы визуализировать поведение настроенной основы правила, постройте поверхность управления каждой нечеткой системы в дереве FIS.

figure('Position',[300 300 600 300]);
subplot(1,2,1)
gensurf(fisTRuleTuned.FIS(1))
title('Control Surface - fis1')
subplot(1,2,2)
gensurf(fisTRuleTuned.FIS(2))
title('Control Surface - fis2')

Следующие таблицы показывают соответствующие основы правила fis1 и fis2.

showRuleBase(fisTRuleTuned.FIS(1))

Rule base of fis1:
                         BG_Rate: N               BG_Rate: Z                BG_Rate: P      
                   ______________________    _____________________    ______________________

    BG_Level: L    Precalculated_Dose: VH    Precalculated_Dose: M    Precalculated_Dose: M 
    BG_Level: M    Precalculated_Dose: VL    Precalculated_Dose: M    Precalculated_Dose: H 
    BG_Level: H    Precalculated_Dose: M     Precalculated_Dose: H    Precalculated_Dose: VL

showRuleBase(fisTRuleTuned.FIS(2))

Rule base of fis2:
                               BG_Accel: N         BG_Accel: Z         BG_Accel: P   
                             ________________    ________________    ________________

    Precalculated_Dose: L    Insulin_Dose: VH    Insulin_Dose: VL    Insulin_Dose: H 
    Precalculated_Dose: M    Insulin_Dose: VL    Insulin_Dose: L     Insulin_Dose: L 
    Precalculated_Dose: H    Insulin_Dose: L     Insulin_Dose: VL    Insulin_Dose: VH

Эти таблицы показывают, что некоторые действия управления неинтуитивны. Например:

Обновите правила путем изменения их последовательных значений.

% Update fis1 rules.
fisTRuleUpdate = fisTRuleTuned;
fisTRuleUpdate.FIS(1).Rules(1).Description = ...
    "BG_Level==L & BG_Rate==N => Precalculated_Dose=VL";
fisTRuleUpdate.FIS(1).Rules(2).Description = ...
    "BG_Level==M & BG_Rate==N => Precalculated_Dose=M";
fisTRuleUpdate.FIS(1).Rules(3).Description = ...
    "BG_Level==H & BG_Rate==N => Precalculated_Dose=H";
fisTRuleUpdate.FIS(1).Rules(9).Description = ...
    "BG_Level==H & BG_Rate==P => Precalculated_Dose=VH";

% Update fis2 rules.
fisTRuleUpdate.FIS(2).Rules(1).Description = ...
    "Precalculated_Dose==L & BG_Accel==N => Insulin_Dose=VL";
fisTRuleUpdate.FIS(2).Rules(3).Description = ...
    "Precalculated_Dose==H & BG_Accel==N => Insulin_Dose=M";
fisTRuleUpdate.FIS(2).Rules(6).Description = ...
    "Precalculated_Dose==H & BG_Accel==Z => Insulin_Dose=M";
fisTRuleUpdate.FIS(2).Rules(7).Description = ...
    "Precalculated_Dose==L & BG_Accel==P => Insulin_Dose=L";

Просмотрите модифицированные основы правила.

showRuleBase(fisTRuleUpdate.FIS(1))

Rule base of fis1:
                         BG_Rate: N               BG_Rate: Z                BG_Rate: P      
                   ______________________    _____________________    ______________________

    BG_Level: L    Precalculated_Dose: VL    Precalculated_Dose: M    Precalculated_Dose: M 
    BG_Level: M    Precalculated_Dose: M     Precalculated_Dose: M    Precalculated_Dose: H 
    BG_Level: H    Precalculated_Dose: H     Precalculated_Dose: H    Precalculated_Dose: VH

showRuleBase(fisTRuleUpdate.FIS(2))

Rule base of fis2:
                               BG_Accel: N         BG_Accel: Z         BG_Accel: P   
                             ________________    ________________    ________________

    Precalculated_Dose: L    Insulin_Dose: VL    Insulin_Dose: VL    Insulin_Dose: L 
    Precalculated_Dose: M    Insulin_Dose: VL    Insulin_Dose: L     Insulin_Dose: L 
    Precalculated_Dose: H    Insulin_Dose: M     Insulin_Dose: M     Insulin_Dose: VH

Визуализируйте получившиеся поверхности управления FIS.

figure('Position',[300 300 600 300]);
subplot(1,2,1)
gensurf(fisTRuleUpdate.FIS(1))
title('Control Surface - fis1')
subplot(1,2,2)
gensurf(fisTRuleUpdate.FIS(2))
title('Control Surface - fis2')

Поверхности управления соответствуют более интуитивному поведению контроллера.

Чтобы проверять, улучшают ли обновленные правила производительность контроллера, симулируйте модель и постройте результаты. Сравните результаты с теми из контроллера настроенными параметрами MF.

fisT = fisTRuleUpdate;
ruleUpdatedOutput = sim(model);

plotGlucoseAndInsulin([ruleTunedOutput ruleUpdatedOutput],...
    'Blood Glucose and Insulin Dosage with Updated Rule Base',...
    {'Learned rules','Updated rules'})

Контроллер с обновленными правилами снижает уровни глюкозы крови по сравнению с настроенным древовидным контроллером FIS.

Обновление правил уменьшает значение функции стоимости.

minCost = costFcn(fisTRuleUpdate,model,minLevel,refLevel,wsVars)

minCost = 23.2702

Настройте параметры функции принадлежности

Чтобы далее улучшать производительность контроллера, можно настроить параметры MF дерева FIS.

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

options.Method = 'patternsearch';
options.MethodOptions.MaxIterations = 10;

По умолчанию каждая входная переменная имеет три равномерно распределенных треугольных MFS. Например, просмотрите MFS для первого входа fis1.

figure
plotmf(fisTRuleUpdate.FIS(1),'input',1)

Сконфигурируйте настраиваемые настройки для входных переменных, таким образом, что крайний левый и самый правый peaks остается неизменным во время настройки.

for i = 1:4
    in(i).MembershipFunctions(1).Parameters.Free = [0 0 1];
    in(i).MembershipFunctions(end).Parameters.Free = [1 0 0];
end

Точно так же сконфигурируйте настраиваемые настройки для выходных переменных. Каждая выходная переменная имеет пять треугольных функций принадлежности.

for i = 1:2
    out(i).MembershipFunctions(1).Parameters.Free = [0 0 1];
    out(i).MembershipFunctions(end).Parameters.Free = [1 0 0];
end

Настройте значения параметров MF с помощью обновленных настраиваемых настроек.

if runtunefis
    figure
    reset(minData)
    rng('default')
    fisTMFTuned = tunefis(fisTRuleUpdate,[in;out],...
        @(fis)costFcn(fis,model,minLevel,refLevel,wsVars,minData),options);
else
    fisTMFTuned = data.fisTMFTuned;
    minCost = costFcn(fisTMFTuned,model,minLevel,refLevel,wsVars)
end

minCost = 22.0627

Настраивающий процесс MF уменьшает значение функции стоимости далее.

Симулируйте модель с помощью контроллера настроенными параметрами MF.

fisT = fisTMFTuned;
mfTunedOutput = sim(model);

Постройте получившиеся отрегулированные уровни глюкозы и скорость введения инсулина. Сравните результаты с теми для контроллера с настроенными основами правила.

plotGlucoseAndInsulin([ruleUpdatedOutput mfTunedOutput],...
    'Blood Glucose and Insulin Dosage with Tuned MFs',...
    {'Updated rules','Tuned MFs'})

Настроенные параметры MF улучшают производительность и уменьшают минимальную величину затрат.

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

% Close model.
close_system(model)

Ссылки

[1] Предоставьте, Пол. “Новый Подход к Диабетическому Управлению: Нечеткая логика и Технология Насоса Инсулина”. Medical Engineering & Physics 29, № 7 (сентябрь 2007): 824–27. https://doi.org/10.1016/j.medengphy.2006.08.014.

[2] Wilinska, M.E., Л.Дж. Чассин, Х.К. Шаллер, Л. Шопп, Т.Р. Пибер и Р. Ховорка. “Кинетика инсулина при диабете 1 типа: Непрерывный и Поставка Шарика Быстрого Действующего Инсулина”. Транзакции IEEE на Биоинженерии 52, № 1 (январь 2005): 3–12. https://doi.org/10.1109/TBME.2004.839639.

[3] Hovorka, римлянин, Валентина Кэнонико, Людовик Ж Шассен, Ульрих Хаютер, Массимо Масси-Бенедетти, Марко Орсини Федеричи, Томас Р Пибер, и др. “Нелинейное Прогнозирующее Управление Модели Концентрации Глюкозы в Предметах с диабетом 1 типа”. Физиологическое Измерение 25, № 4 (1 августа 2004): 905–20. https://doi.org/10.1088/0967-3334/25/4/010.

[4] Mauseth, Ричард, Юкинг Ван, Eyal Dassau, Роберт Киркэр, Дональд Мэтисон, Говард Зиссер, Лоис Jovanovič и Фрэнсис Дж. Дойл. “Предложенное Клиническое Приложение для Настройки Контроллера Нечеткой логики Искусственной Поджелудочной железы, Использующей Фактор Персонализации”. Журнал Науки и техники Диабета 4, № 4 (июль 2010): 913–22. https://doi.org/10.1177/193229681000400422.

Функции помощника

function fis = updateInput(fis,id,name,range,mfNames)
% Update FIS input with the specified parameter values.

fis.Inputs(id).Name = name;
fis.Inputs(id).Range = range;

for mfId = 1:length(mfNames)
    fis.Inputs(id).MembershipFunctions(mfId).Name = mfNames(mfId);
    params = range(1) + ...
        diff(range)*fis.Inputs(id).MembershipFunctions(mfId).Parameters;
    fis.Inputs(id).MembershipFunctions(mfId).Parameters = params;
end

end

function fis = updateOutput(fis,id,name,range)
% Update FIS output with the specified parameter values.

rangeDiff = diff(range);
fis.Outputs(id).Name = name;

% MF names - Very Low, Low, Medium, High, Very High
mfNames = [...
    "VL","L","M","H","VH"];

for mfId = 1:length(mfNames)
    fis.Outputs(id).MembershipFunctions(mfId).Name = mfNames(mfId);
    params = range(1) + ...
        rangeDiff*fis.Outputs(id).MembershipFunctions(mfId).Parameters;
    fis.Outputs(id).MembershipFunctions(mfId).Parameters = params;
end

% Extend output range values to fit the output MFs.
left = fis.Outputs(id).MembershipFunctions(1).Parameters(1);
right = fis.Outputs(id).MembershipFunctions(end).Parameters(end);
fis.Outputs(id).Range = [left right];

end