2D напрягите модель туберкулеза

Модель туберкулеза 2D деформации введена в [1]. Модель описана можно следующим образом:

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

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

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

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

Для получения дополнительной информации о динамической модели смотрите ResistantTBStateFcn.m.

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

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

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

Здесь, вес B1 50, и весом B2 является 500. Эти веса подчеркивают настройку выявлению заболевания по случаю, содержащему из-за его удара стоимости.

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

S = 19 000 T = 250 L2 = 500 I1 = 1 000 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. Установите эти значения как нижние и верхние границы для MVS.

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 функция. При начальном условии MVS является нулем. По умолчанию, 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. 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] Юнг, E., С. Ленхарт и Цз. Фэн. "Оптимальное управление над Обработками в Модели Туберкулеза 2D Деформации". Дискретные и Непрерывные Динамические Системы, Серия B2, 2002, стр 479-482.