Модель туберкулеза с двумя штаммами

Двухтеневая модель туберкулеза представлена в [1]. Модель описывается следующим образом:

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

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

Шестой класс, L1 (латентный, инфицированный типичным туберкулезом, но не инфекционный) рассчитывается как N - (S+T+L2+I1+I2).

Можно уменьшить число случаев устойчивого ТБ с помощью двух управляемых переменных (MV):

Дополнительные сведения о динамической модели см. в разделе ResistantTBStateFcn.m.

Проблема нелинейного управления

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

F = L2 + I2 + 0.5*B1*u1^2 + 0.5*B2*u2^2

Вот, вес B1 является 50, а вес B2 равен 500. Эти веса подчеркивают предпочтение поиска дела, а не хранения дела из-за его влияния на стоимость.

Общая численность населения составляет N = 30000. При исходном условии:

S = 19000 T = 250 L2 = 500 I1 = 1000 I2 = 250

которая уходит L1 = 9000.

N = 30000;
x0 = [76; 1; 2; 4; 1]*N/120;

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

nx = 5;
ny = nx;
nu = 2;
nlobj = nlmpc(nx,ny,nu);

In standard cost function, zero weights are applied by default to one or more OVs because there are fewer MVs than OVs.

Предположим, что политика лечения может корректироваться только каждые три месяца. Поэтому установите время выборки контроллера равным 0,25 года. Так как вы хотите найти оптимальную политику в течение пяти лет, установите горизонт прогнозирования на 20 шагов (5 лет, разделенных на 0,25).

Years = 5;
Ts = 0.25;
p = Years/Ts;
nlobj.Ts = Ts;
nlobj.PredictionHorizon = p;

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

nlobj.ControlHorizon = p;

Модель прогнозирования определена в ResistantTBStateFcn.m. Укажите эту функцию в качестве функции состояния контроллера.

nlobj.Model.StateFcn = 'ResistantTBStateFcn';

Рекомендуется указывать аналитические якобинские функции для модели прогнозирования и функции затрат/ограничений. Дополнительные сведения о вычислении якобиана для уравнений состояний см. в разделе ResistantTBStateJacFcn.m. Установите этот файл как функцию якобиана состояния.

nlobj.Jacobian.StateFcn = 'ResistantTBStateJacFcn';

Поскольку все состояния - это количество индивидуумов, они должны быть неотрицательными значениями. Укажите минимальную границу 0 для всех штатов.

for ct = 1:nx
    nlobj.States(ct).Min = 0;
end

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

for ct = 1:nx
    nlobj.States(ct).ScaleFactor = x0(ct);
end

Как «поиск», так и «удержание» органов управления имеют рабочий диапазон между 0.05 и 0.95. Установите эти значения в качестве нижней и верхней границ для MV.

nlobj.MV(1).Min = 0.05;
nlobj.MV(1).Max = 0.95;
nlobj.MV(2).Min = 0.05;
nlobj.MV(2).Max = 0.95;

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

nlobj.Optimization.CustomCostFcn = "ResistantTBCostFcn";
nlobj.Optimization.ReplaceStandardCost = true;

Кроме того, для ускорения моделирования аналитический якобиан стоимости предусмотрен в ResistantTBCostJacFcn.

nlobj.Jacobian.CustomCostFcn = "ResistantTBCostJacFcn";

С момента L1 население определяется как N минус сумма всех состояний, вы должны убедиться, что (S+T+L2+I1+I2) - N < 0 всегда удовлетворен. В nlmpc укажите это условие как ограничение неравенства с помощью анонимной функции.

nlobj.Optimization.CustomIneqConFcn = @(X,U,e,data) sum(X(2:end,:),2)-30000;

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

validateFcns(nlobj,x0,[0.5;0.5])

Model.StateFcn is OK.
Jacobian.StateFcn is OK.
No output function specified. Assuming "y = x" in the prediction model.
Optimization.CustomCostFcn is OK.
Jacobian.CustomCostFcn is OK.
Optimization.CustomIneqConFcn is OK.
Analysis of user-provided model, cost, and constraint functions complete.

Для вычисления оптимальной политики управления используйте nlmpcmove функция. В исходном состоянии MV равны нулю. По умолчанию fmincon из Toolbox™ Оптимизация (Optimization) используется в качестве решателя NLP по умолчанию.

lastMV = zeros(nu,1);
[~,~,Info] = nlmpcmove(nlobj,x0,lastMV);

Slack variable unused or zero-weighted in your custom cost function. All constraints will be hard.

Постройте график и изучите оптимальное решение.

ResistantTBPlot(Info,Ts)

Optimal cost =   5194.8
Final     L2 =    175.9
Final     I2 =    862.9

Оптимальное решение дает стоимость 5195и общее число людей, инфицированных резистентным туберкулезом в последнее время, составляет L2 + I2 = 1037.

Поиск оптимальной обработки с помощью пользовательского решателя нелинейного программирования

Если вы хотите использовать сторонний решатель NLP в моделировании, напишите файл интерфейса, который преобразует входные данные, определенные nlmpc во входы, определенные решателем NLP, и укажите его как CustomSolverFcn в nlmpc объект.

В этом примере предположим, что у вас есть решатель «XYZ», который имеет интерфейс пользователя, отличный от fmincon. A ResistantTBSolver.m создается файл для преобразования задачи оптимизации, определенной nlmpc объект к соответствующему интерфейсу, требуемому решателем «XYZ». Просмотрите ResistantTBSolver.m.

function [zopt,cost,exitflag] = ResistantTBSolver(FUN,z0,A,b,Aeq,beq,lb,ub,NONLINCON)
% This is an interface function that calls a XYZ nonlinear programming
% solver to solve an NLP problem defined by an nlmpc controller. The
% input and output definitions of this function are identical to fmincon.
%
% This interface function converts fmincon inputs to the format required
% by the XYZ solver and converts the results back to the fmincon outputs.
%
% Inputs
%       FUN: nonlinear cost. FUN accepts input z (decision variabls) and
%            returns cost f (a scalar value evaluated at z) and its
%            gradient g (a nz-by-1 vector evaluated at z).  
%        z0: initial guess of z
%       A,b: A*z<=b
%   Aeq,beq: Aeq*z==beq
%     lb,ub: lower and upper bounds of z
% NONLINCON: nonlinear constraints. NONLINCON accepts input z and returns
%            the vectors C and Ceq as the first two outputs, representing
%            the nonlinear inequalities and equalities respectively where
%            FUN is minimized subject to C(z) <= 0 and Ceq(z) = 0.
%            NONLINCON also returns Jacobian matrices of C (a nz-by-ncin
%            sparse matrix) and Ceq (a nz-by-nceq sparse matrix).
%
% Outputs
%     zopt:  optimal solution of z
%     cost:  optimal cost at optimal z
% exitflag:
%            1  First order optimality conditions satisfied.
%            0  Too many function evaluations or iterations.
%           -1  Stopped by output/plot function.
%           -2  No feasible point found.
%
% You can use this example as a template to implement an interface file to
% your own NLP solver. The solver must be an M file or a MEX file on the
% MATLAB path.
%
% DO NOT EDIT LINES ABOVE.
%
% Copyright 2018 The MathWorks, Inc.

%% Set dimensional information of linear and nonlinear constraints
num_lin_ineq = size(A,1);
num_lin_eq = size(Aeq,1);
[in,eq] = NONLINCON(z0);
num_non_ineq = length(in);
num_non_eq = length(eq);
total = num_non_ineq + num_non_eq + num_lin_ineq + num_lin_eq;
logicals_nlineq = false(total,1);
logicals_nlineq(1:num_non_ineq) = true;
logicals_nleq = false(total,1);
logicals_nleq(num_non_ineq+(1:num_non_eq)) = true;
logicals_ineq = false(total,1);
logicals_ineq(num_non_ineq+num_non_eq+(1:num_lin_ineq)) = true;
logicals_eq = false(total,1);
logicals_eq(num_non_ineq+num_non_eq+num_lin_ineq+(1:num_lin_eq)) = true;
options = struct('nlineq',logicals_nlineq,'nleq',logicals_nleq,...
                 'ineq',logicals_ineq,'eq',logicals_eq);
%% Set decision variable bounds
options.lb = lb;
options.ub = ub;
%% Set RHS of nonlinear constraints
options.cl = [-inf(num_non_ineq,1);zeros(num_non_eq,1)];
options.cu = [zeros(num_non_ineq,1);zeros(num_non_eq,1)];
%% Set RHS of linear constraints
options.rl = [-inf(num_lin_ineq,1);beq];
options.ru = [b;beq];
%% Set A matrix
options.A = sparse([A; Aeq]);
%% Set XYZ solver optons
options.algorithm = struct('print_level',0,'max_iter',200,...
                           'max_cpu_time',1000,'tol',1.0000e-06,...
                           'hessian_approximation','limited-memory');
%% Set function handles used by XYZ solver
Jstr = sparse(ones(num_non_ineq+num_non_eq,length(z0)));
funcs = struct('objective',@(x) fval(FUN,x),...
              'gradient',@(x) gval(FUN,x),...
              'constraints',@(x) conval(NONLINCON,x),...
              'jacobian',@(x) jacval(NONLINCON,x),...
              'jacobianstructure',@() Jstr...
              );
%% Call XYZ and return cost and status
[zopt,output] = XYZsolver(z0,funcs,options);
cost = FUN(zopt);
exitflag = convertStatustoExitflag(output.status);

%% Utility functions
function f = fval(fun,z)
% Return nonlinear cost
[f,~] = fun(z);

function g = gval(fun,z)
% Return cost gradient
[~,g] = fun(z);

function c = conval(nonlcon,z)
% Return nonlinear constraints
[in,eq] = nonlcon(z);
c = [in;eq];

function J = jacval(nonlcon,z)
% Return constraints Jacobian as nc-by-nz in sparse matrix
% Jin is nz-by-ncin sparse, Jeq is nz-by-nceq sparse
[~,~,Jin,Jeq] = nonlcon(z); 
J = [Jin Jeq]';

function exitflag = convertStatustoExitflag(status)
switch(status)
    case 0
        %info.Status = 'Success';
        exitflag = 1;
    case 1
        %info.Status = 'Solved to Acceptable Level';
        exitflag = 1;
    case 2
        %info.Status = 'Infeasible';
        exitflag = -1;
    case 3
        %info.Status = 'Search Direction Becomes Too Small';
        exitflag = -2;
    case 4
        %info.Status = 'Diverging Iterates';
        exitflag = -2;
    case 6
        %info.Status = 'Feasible Point Found';
        exitflag = 1;
    case -1
        %info.Status = 'Exceeded Iterations';
        exitflag = 0;
    case -4 
        %info.Status = 'Max Time Exceeded';
        exitflag = 0;
    otherwise        
        %info.Status = 'IPOPT Error';
        exitflag = -2;
end

Этот файл можно использовать в качестве шаблона для реализации файла интерфейса в собственном решателе NLP. Решатель должен быть скриптом MATLAB или файлом MEX по пути MATLAB.

Можно подключить решатель, указав его как пользовательский решатель в nlmpc объект.

nlobj.Optimization.CustomSolverFcn = @ResistantTBSolver;

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

Ссылки

[1] Юнг, Э., С. Ленхарт и З. Фэн. «Оптимальный контроль лечения в модели туберкулеза с двумя штаммами». Дискретные и непрерывные динамические системы, серия B2, 2002, стр. 479-482.