主要内容

基于自定义求解器的非线性MPC优化结核病治疗

这个例子使用了:

打开脚本

这个例子展示了如何通过构造一个非线性MPC问题来寻找治疗双株结核病(TB)人群的最优策略。然后，非线性MPC控制器使用默认求解器和自定义求解器计算最优解。

双株结核模型

在[1]中引入了一个双株结核模型。模型描述如下:

“在缺乏有效疫苗的情况下，目前的结核病控制项目集中在化疗上。对活动性结核病(药物敏感株)患者进行抗生素治疗所需的时间和费用比那些感染了敏感结核病但尚未发病的患者要长得多。不遵守药物治疗不仅可能导致复发，而且可能导致耐药性结核病的发展——这是当今社会面临的最严重的公共卫生问题之一。药物敏感结核病病例的减少可以通过“病例保持”或“病例发现”来实现，“病例保持”指的是为确保在足够的时间内有规律地服用药物而采取的活动和技术，或通过“病例发现”来实现，“病例发现”指的是(例如，通过筛查)确定潜在感染敏感结核病的个人，这些人具有发展为该病的高风险，可能受益于预防性干预措施。这些预防性治疗将减少药物敏感结核病的发病率(每单位时间的新病例)，从而间接减少耐药结核病的发病率。”

在本例中使用的动态模型中，总宿主种群N(这是一个常数)被分为六个不同的流行病学类别。其中五个类被定义为状态变量:

  • x (1) =年代，易感人群的数量

  • x (2) =T有效治疗的个体数量

  • x (3) =L2潜伏期，感染耐药结核病，但不具传染性

  • x (4) =I1与典型的结核病有传染性

  • x (5) =I2可感染耐药结核病

第六节课，L1(潜伏期，感染典型结核病但不具传染性)计算为N - (s + t + l2 + i1 + i2)．

可以使用两个操纵变量(mv)减少耐药结核病病例:

  • U(1) -“病例发现”，花费相对廉价的精力来确定需要治疗的人。

  • U(2) -“病例保持”，维持有效治疗相对昂贵的努力。

有关动态模型的更多信息，请参见ResistantTBStateFcn.m．

非线性控制问题

结核病治疗的目标是在五年期间减少耐药结核病的潜伏(L2)和传染性(I2)个体，同时保持低成本。为了实现这一目标，可以使用一个成本函数来计算以下五年的价值。

F = L2 + I2 + 0.5*B1*u1²+ 0.5*B2*u2²

在这里,重量B1是50，权重B2为500．这些权重强调，由于其成本影响，案件查找优先于案件保存。

总人口是N＝30000．在初始条件下:

年代= 19000T= 250L2= 500I1= 1000I2= 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);
在标准成本函数中，默认情况下对一个或多个OVs应用零权重，因为mv比OVs少。

假设治疗政策只能每三个月调整一次。因此，将控制器采样时间设置为0.25年。因为您想找到五年的最优策略，所以将预测范围设置为20步(5年除以0.25)。

年= 5年;Ts = 0.25;p =年/Ts;nlobj。Ts = Ts; nlobj.PredictionHorizon = p;

对于这个规划问题，您希望使用决策变量的最大数量。为此，将控制水平设置为与预测水平相等。

nlobj。ControlHorizon = p;

中定义了预测模型ResistantTBStateFcn.m．指定此函数作为控制器状态函数。

nlobj.Model.StateFcn =“ResistantTBStateFcn”；

最好的做法是为预测模型和成本/约束函数指定解析雅可比函数。有关状态方程的雅可比矩阵计算的详细信息，请参见ResistantTBStateJacFcn.m．将此文件设置为状态雅可比函数。

nlobj.Jacobian.StateFcn =“ResistantTBStateJacFcn”；

因为所有的状态都是个体的数量，所以它们必须是非负值。的最小界0适用于所有州。

为ct = 1:nx nlobj.States(ct)。Min = 0;结束

因为在组(状态)之间有很大的总体变化，所以使用它们各自的名义值缩放状态变量。这样做提高了优化问题的数值鲁棒性。

为ct = 1:nx nlobj.States(ct)。ScaleFactor = x0(ct);结束

“寻找”和“持有”控制都有一个操作范围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;

方法检查潜在的数值问题，请验证预测模型、自定义函数和雅可比矩阵validateFcns命令。

validateFcns (nlobj x0, [0.5, 0.5])
模型。年代tateFcn 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从优化工具箱™被用作默认的NLP求解器。

lastMV = 0 (nu,1);[~，~，Info] = nlmpcmove(nlobj,x0,lastMV);
松弛变量未使用或自定义成本函数中的零加权。所有的约束都将是艰难的。

绘制和检查最优解决方案。

ResistantTBPlot(信息,Ts)
最优成本= 5194.8最终L2 = 175.9最终I2 = 862.9

最优解的成本为5195，最后一次感染耐药结核病的总人数为L2+I2＝1037．

使用自定义非线性规划求解器寻找最佳处理方法

如果您想在模拟中使用第三方NLP求解器，请编写一个接口文件来转换由定义的输入nlmpc输入由NLP求解器定义的输入，并将其指定为CustomSolverFcn在nlmpc对象。

在本例中，假设您有一个“XYZ”求解器，其用户界面与fmincon．一个ResistantTBSolver.m文件来转换定义的优化问题nlmpc对象到“XYZ”求解器所需的适当接口。检查ResistantTBSolver.m．

函数[zopt,cost,exitflag] =电阻ttbsolver (FUN,z0,A,b,Aeq,beq,lb,ub,NONLINCON)这是一个调用XYZ非线性编程的接口函数。%求解器来解决由nlmpc控制器定义的NLP问题。的此函数的输入和输出定义与fmincon相同。％这个接口函数将fmincon输入转换为所需的格式%，并将结果转换回fmincon输出。％%的输入% FUN:非线性成本。FUN接受输入z(决策变量)和%返回成本f(在z处计算的标量值)及其%梯度g(一个在z处求值的nz × 1向量)。% z0: z的初始猜测% A,b: A*z<=b% Aeq,beq: Aeq*z==beq% lb,ub: z的上下界NONLINCON:非线性约束。NONLINCON接受输入z并返回%，向量C和Ceq作为前两个输出，表示%的非线性不等式和不等式分别在其中% FUN在C(z) <= 0和Ceq(z) = 0的条件下最小化。NONLINCON也返回C (a nz × ncin)的雅可比矩阵%稀疏矩阵)和Ceq(一个nz-by-nceq稀疏矩阵)。％%输出% zopt: z的最优解%成本:最佳z处的最佳成本% exitflag:% 1满足一阶最优条件。% 0太多的函数计算或迭代。% -1输出/绘图函数停止。% -2找不到可行点。％您可以使用此示例作为模板来实现接口文件%你自己的NLP求解器。的解算器必须是M文件或MEX文件MATLAB路径。％%不要编辑上面的行。％版权所有The MathWorks, Inc.设置线性和非线性约束的尺寸信息num_lin_ineq = size(A,1);num_lin_eq = size(Aeq,1);[in,eq] = NONLINCON(z0);Num_non_ineq =长度(in);Num_non_eq =长度(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,“情商”, logicals_eq);设置决策变量边界选项。Lb = Lb;选项。Ub = Ub;设置非线性约束的RHS选项。Cl = [-inf(num_non_ineq,1);zero (num_non_eq,1)];选项。Cu = [zeros(num_non_ineq,1);zeros(num_non_eq,1)];设置线性约束的RHS选项。Rl = [-inf(num_lin_ineq,1);beq];Options.ru = [b;beq];设置一个矩阵选项。一个＝sparse([A; Aeq]);设置XYZ求解选项选项。算法= struct(“print_level”0,“max_iter”, 200,.．.“max_cpu_time”, 1000,“托尔”1.0000 e-06.．.“hessian_approximation”，内存有限的）;设置XYZ求解器使用的函数句柄Jstr = sparse(ones(num_non_ineq+num_non_eq,length(z0)));Funcs = struct(“目标”@ (x) fval(有趣,x).．.“梯度”@ (x) gval(有趣,x).．.“约束”@ (x) conval (NONLINCON x),.．.的雅可比矩阵@ (x) jacval (NONLINCON x),.．.“jacobianstructure”Jstr @ ().．.）;致电XYZ并返回成本和状态[zopt,output] = XYZsolver(z0,funcs,options);成本=乐趣(zopt);exitflag = convertStatustoExitflag(output.status);实用函数函数F = fval(fun,z)回报率非线性成本[f，~] = fun(z);函数G = gval(fun,z)返回成本梯度[~，g] = fun(z);函数C = conval(nonlcon,z)返回非线性约束[in,eq] = nonlcon(z);C = [in;eq];函数J = jacval(nonlcon,z)%返回约束雅可比矩阵作为稀疏矩阵的nc × nzJin是稀疏的，Jeq是稀疏的[~，~，Jin,Jeq] = nonlcon(z);J = [Jin Jeq]';函数exitflag = convertStatustoExitflag(状态)开关(状态)情况下0%的信息。年代tatus = 'Success';Exitflag = 1;情况下1%的信息。年代tatus = 'Solved to Acceptable Level';Exitflag = 1;情况下2%的信息。年代tatus = 'Infeasible';Exitflag = -1;情况下3.%的信息。年代tatus = 'Search Direction Becomes Too Small';Exitflag = -2;情况下4%的信息。年代tatus = 'Diverging Iterates';Exitflag = -2;情况下6%的信息。年代tatus = 'Feasible Point Found';Exitflag = 1;情况下－1%的信息。年代tatus = 'Exceeded Iterations';Exitflag = 0;情况下4%的信息。年代tatus = 'Max Time Exceeded';Exitflag = 0;否则%的信息。年代tatus = 'IPOPT Error';Exitflag = -2;结束

您可以使用此文件作为模板，实现您自己的NLP求解器的接口文件。求解器必须是MATLAB脚本或MATLAB路径上的MEX文件。

中的自定义求解器，可以将其指定为nlmpc对象。

nlobj.Optimization.CustomSolverFcn = @电阻ttbsolver;

只要“XYZ”求解器是可靠的，它的选项是正确选择的，重新运行模拟应该会产生类似的结果。

参考文献

[1]荣格，E.， S.伦哈特，冯铮。“双株结核模型治疗的最优控制”。离散和连续动力系统， B2系列，2002，第479-482页。

另请参阅

|

相关的话题