Двухштаммная модель туберкулеза

Модель туберкулеза с двумя штаммами введена в [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. Установите эти значения как нижнюю и верхнюю границы для СН.

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;

Чтобы проверить на потенциальные числовые проблемы, проверьте свою модель предсказания, пользовательские функции и якобианцы, используя 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 из Optimization Toolbox™ используется в качестве решателя 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.