使用非线性优化结核病治疗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) +。
您可以使用两个操纵变量减少耐药结核病例(MVs):
u(1)——“发现”,相对廉价的精力来识别那些需要治疗。
u(2)——“控股”,相对昂贵的努力维持有效的治疗。
动态模型的更多信息,请参阅ResistantTBStateFcn.m。
非线性控制问题
结核病治疗的目的是减少潜在(L2)和传染性(I2)耐药结核病患者5年期间,同时保持低成本。为了实现这一目标,使用一个成本函数,总结以下值超过五年。
F = L2 + I2 u1 ^ 2 + 0.5 + 0.5 * B1 * * B2 * u2 ^ 2
在这里,重量B1是50和重量B2500年。这些权重强调偏爱发现在案件由于其成本的影响。
人口总数是N=30000年。初始条件:
年代= 19000T= 250L2= 500I1= 1000I2= 250
这让L1= 9000。
N = 30000;x0 = [76;1;2;4;1)* N / 120;
在这个例子中,找到最优使用非线性MPC控制器控制策略。创建nlmpc对象的正确数量的州、输出和输入。
nx = 5;纽约= nx;ν= 2;nlobj = nlmpc (nx、纽约、ν);
在标准成本函数、零重量默认应用到一个或多个ov因为有MVs比ov少。
假设治疗政策只能每三个月调整。因此,设置控制器样品时间0.25年。既然你想找到最优政策在5年内,将预测地平线20个步骤(5年除以0.25)。
年= 5;t = 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)。最小值= 0;结束
因为有一个庞大的人口变化的组织(州),规模状态变量使用各自的标称值。这样做提高了数值优化问题的鲁棒性。
为ct = 1: nx nlobj.States (ct)。ScaleFactor = x0 (ct);结束
“发现”和“持有”控制操作区间0.05和0.95。这些值设置为低,MVs上界。
nlobj.MV (1)。最小值= 0.05;nlobj.MV (1)。Max = 0.95;nlobj.MV (2)。最小值= 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,数据)和(X (2:,:), 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函数。在初始条件,MVs为零。默认情况下,fmincon从优化工具箱™用作默认的NLP的能手。
lastMV = 0(ν,1);[~,~,信息]= nlmpcmove (x0, nlobj lastMV);
松弛变量未使用或zero-weighted定制成本函数。所有约束将是非常困难的。
情节和检查的最佳解决方案。
ResistantTBPlot(信息,Ts)
最优成本= 5194.8最终L2最终I2 = 862.9 = 175.9
收益成本的最优解决方案5195年,个体的总数感染耐药型结核病的出现在最后的时间L2+I2=1037年。
找到最佳的治疗使用自定义非线性规划求解器进行求解
如果你想使用第三方NLP仿真解算器,写一个接口文件,该文件将定义的输入nlmpc为定义的输入你的NLP的能手,并指定的CustomSolverFcn在nlmpc对象。
在这个例子中,假设您有一个“XYZ”解算器,有一个不同的用户界面fmincon。一个ResistantTBSolver.m创建文件将定义的优化问题nlmpc对象的接口所需的解算器“XYZ”。检查ResistantTBSolver.m。
函数[zopt、成本exitflag] = ResistantTBSolver(乐趣、z0 A、b Aeq,说真的,磅,乌兰巴托,NONLINCON)%,这是一个接口函数调用XYZ非线性规划%解算器来解决一个NLP问题定义为一个nlmpc控制器。的%的输入和输出定义这个函数fmincon完全相同。%%这个接口函数转换fmincon输入所需的格式%的XYZ解算器并将结果返回给fmincon输出。%%的输入%的乐趣:非线性成本。有趣的z(决定variabls)和接受输入%返回成本f(标量值评估在z)和它的%梯度g (nz-by-1向量计算z)。% z0:初始猜测z% A, b: A * z < =% Aeq,说真的:Aeq * z = =说真的%磅,乌兰巴托:低和z的上界% NONLINCON:非线性约束。NONLINCON接受输入z和返回%的向量C测查,前两个输出,代表%的非线性不等式和等式%的乐趣是最小化C (z) < = 0和量表(z) = 0。% NONLINCON也返回C (nz-by-ncin的雅可比矩阵%稀疏矩阵)测查,(nz-by-nceq稀疏矩阵)。%%输出% zopt: z的最优解z %成本:最优成本最优% exitflag:% 1满足一阶最优性条件。% 0太多函数评估或迭代。% 1停止输出/情节功能。% 2没有可行点。%%可以使用这个例子作为一个模板来实现一个接口文件%你自己的NLP的能手。解算器必须一个M文件或文件的墨西哥人% MATLAB路径。%%不编辑线以上。%% 2018年版权MathWorks公司。% %设置维度信息的线性和非线性约束num_lin_ineq =大小(1);num_lin_eq =大小(Aeq, 1);(eq) = NONLINCON (z0);num_non_ineq =长度();num_non_eq =长度(eq);总= num_non_ineq + num_non_eq + num_lin_ineq + num_lin_eq;logicals_nlineq = false(总,1);logicals_nlineq (1: num_non_ineq) = true;logicals_nleq = false(总,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);% %组决策变量选项。磅=磅;选项。乌兰巴托=乌兰巴托;% %设置RHS非线性约束选项。cl =[无穷(num_non_ineq 1); 0 (num_non_eq 1)];选项。铜= [0 (num_non_ineq 1); 0 (num_non_eq 1)];% %设置RHS线性约束选项。rl =(负无穷(num_lin_ineq 1); beq);options.ru = [b; beq];% %设置一个矩阵选项。一个=sparse([A; Aeq]);% % XYZ求解任一用途选项。算法=结构(“print_level”0,“max_iter”,200,…“max_cpu_time”,1000,“托尔”1.0000 e-06…“hessian_approximation”,内存有限的);% %集函数处理使用XYZ解算器Jstr =稀疏的(num_non_ineq + num_non_eq,长度(z0)));func =结构(“目标”@ (x) fval(有趣,x)…“梯度”@ (x) gval(有趣,x)…“约束”@ (x) conval (NONLINCON x),…的雅可比矩阵@ (x) jacval (NONLINCON x),…“jacobianstructure”Jstr @ ()…);% %调用XYZ和返回成本和地位[zopt、输出]= XYZsolver (z0、函数、期权);成本=乐趣(zopt);exitflag = convertStatustoExitflag (output.status);% %效用函数函数f = fval(有趣,z)%返回非线性成本[f ~] =乐趣(z);函数g = gval(有趣,z)%返回成本梯度(~ g) =乐趣(z);函数c = conval (nonlcon, z)%返回非线性约束(eq) = nonlcon (z);c = (, eq);函数J = jacval (nonlcon, z)%返回约束雅可比矩阵作为nc-by-nz稀疏矩阵%金nz-by-ncin稀疏,Jeq nz-by-nceq稀疏[~,~,金,Jeq] = nonlcon (z);J =(金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脚本或墨西哥人文件的路径。
你可以插入的解算器指定的自定义解决nlmpc对象。
nlobj.Optimization。CustomSolverFcn=@ResistantTBSolver;
只要“XYZ”解决者是可靠的,其选择是正确选择,重新运行仿真将产生类似的结果。
引用
[1]荣格,E。,年代。Lenhart, and Z. Feng. "Optimal Control of Treatments in a Two-Strain Tuberculosis Model."离散和连续动力系统系列B2, 2002年,页479 - 482。
