使用自定义求解器的非线性MPC优化肺结核治疗
这个例子展示了如何通过构造一个非线性MPC问题来找到治疗两株结核病（TB）人群的最优策略。然后，非线性MPC控制器使用默认解算器和自定义解算器来计算最优解。
双毒株结核模型
在[1]中引入了一个双菌株结核模型。模型描述如下:
在缺乏有效疫苗的情况下，目前的结核病控制计划集中于化疗。对活动性结核病（具有药物敏感性菌株）的抗生素治疗与那些感染了敏感结核病但尚未发病的患者相比，患者需要更长的时间和更高的成本。缺乏药物治疗的依从性不仅可能导致复发，还可能导致耐药性结核病的发展——这是当今社会面临的最严重的公共卫生问题之一药物敏感性结核病病例的减少可以通过“病例保持”（case holding）来实现，这是指用于确保在足够的时间内有规律地服用药物以实现治愈的活动和技术，或者通过“病例发现”（case finding）来实现，这是指鉴定（例如通过筛查）潜在感染敏感型结核病的高危人群，他们可能受益于第一代码块的预防性干预说明。这些预防性治疗将降低药物敏感型结核病的发病率（每单位时间的新病例），从而间接降低耐药型结核病的发病率。”
在本例中使用的动态模型中，宿主总数N(这是一个常量)被分为六个不同的流行病学类别。其中五个类被定义为状态变量:
x（1）=
s，易感个体的数量x（2）=
T，有效治疗的人数x（3）=
L2，潜伏，感染耐药结核，但不具传染性x (4) =
I1，感染典型结核病x (5) =
I2，感染耐药结核病
第六节课，L1（潜伏，感染典型结核，但不具有传染性）计算为N - (S + T + L2 + I1和I2) +。
你可以使用两个操纵变量(MVs)来减少耐药结核病病例:
u（1）-“病例发现”，花费相对低廉的努力来确定需要治疗的患者。
u（2）-“病例保留”，维持有效治疗的相对昂贵的努力。
有关动态模型的更多信息，请参见电阻TbStatefCn.m。
非线性控制问题
结核病治疗的目标是在五年期间减少潜伏(L2)和传染性(I2)耐药结核病患者，同时保持低成本。为了实现这一目标，使用一个成本函数，它在五年内合计以下价值。
F=L2+I2+0.5*B1*u1^2+0.5*B2*u2^2
这里，重量地下一层是50，权重B2是500. 由于成本影响，这些权重强调案例发现优先于案例持有。
总人口为N=30000A.t the initial condition:
s= 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=Years/Ts；nlobj.Ts=Ts；nlobj.PredictionHorizon=p；
对于这个规划问题，您希望使用最大数量的决策变量。为此，将控制视界设置为与预测视界相等。
nlobj.ControlHorizon=p；
预测模型定义在电阻TbStatefCn.m。将此函数指定为控制器状态函数。
nlobj.Model.StateFcn =“ResistantTBStateFcn”；
最好为预测模型和成本/约束函数指定分析雅可比函数。有关状态方程雅可比计算的详细信息，请参阅电阻TbStateJacfcn.m. 将此文件设置为状态雅可比函数。
nlobj.Jacobian.StateFcn =“ResistantTBStateJacFcn”；
因为所有的状态都是个体的数目，它们必须是非负值。指定的最小界限0适用于所有国家。
对于ct=1:nx非直瞄状态（ct）.Min=0；终止
因为在组(状态)之间存在很大的总体变异，使用它们各自的名义值对状态变量进行缩放。这样做可以提高优化问题的数值鲁棒性。
对于ct=1:nx nlobj.States（ct）.ScaleFactor=x0（ct）；终止
“查找”和“保持”控件的工作范围都在0.05和0.95。将这些值设置为mv的上下限。
nlobj.MV(1)。最小值= 0.05;nlobj.MV(1)。Max = 0.95;nlobj.MV(2)。最小值= 0.05;nlobj.MV(2)。Max = 0.95;
使结核病人口和治疗成本最小化的成本函数定义为电阻TbCostFCn.m.由于此规划问题不需要参考跟踪或干扰抑制，因此使用此成本函数替换标准成本。
nlobj.Optimization.CustomCostFcn =“电阻TbCostFcn”；nlobj.Optimization.ReplaceStandardCost = true;
此外，为了加快仿真速度，中提供了成本的分析雅可比矩阵电阻TbCostJacfcn。
nlobj.Jacobian.CustomCostFcn=“电阻TbCostJacfcn”；
自从L1人口的定义是N减去所有状态的和，你必须确保（S+T+L2+I1+I2）-N<0总是满意的。在nlmpc对象，使用匿名函数将此条件指定为不等式约束。
nlobj.Optimization.CustomIneqConFcn = @(X,U,e,data) sum(X(2:end，:)，2)-30000;
要检查潜在的数值问题，请使用验证EFCNS命令
validateFcns (nlobj x0, [0.5, 0.5])
Model.StateFcn正常。Jacobian.StateFcn正常。未指定输出函数。假设预测模型中的“y=x”。Optimization.CustomCostFcn正常。Jacobian.CustomCostFcn正常。Optimization.CustomIneqConFcn正常。用户提供的模型、成本和约束函数的分析完成。
要计算最优控制策略，请使用nlmpcmove函数。在初始条件下，mv为零。默认情况下,fmincon从优化工具箱™ 用作默认NLP解算器。
lastMV = 0(ν,1);[~, ~,信息]= nlmpcmove (x0, nlobj lastMV);
闲置变量未使用或自定义成本函数中的零加权。所有的限制都是困难的。
绘制并检查最佳解决方案。
电阻TbPlot（信息，Ts）
最优成本=5194.8最终L2=175.9最终I2=862.9
最优解决方案的成本为5195，最终感染耐药结核病的总人数为L2+I2=1037。
使用自定义非线性规划求解器找到最佳处理方法
如果您想在模拟中使用第三方NLP求解器，请编写一个接口文件来转换由nlmpc输入NLP求解器定义的输入，并将其指定为CustomSolverFcn在nlmpc对象
在本例中，假设您有一个“XYZ”解算器，该解算器的用户界面与fminconA.电阻TbSolver.m创建文件以转换由定义的优化问题nlmpc对象转换为“XYZ”求解器所需的适当接口。检查电阻TbSolver.m。
作用[zopt，cost，exitflag]=电阻TbSolver（乐趣，z0，A，b，Aeq，beq，lb，ub，非林肯）%这是一个调用XYZ非线性编程的接口函数%用于解决nlmpc控制器定义的NLP问题的解算器该函数的输入和输出定义与fmincon相同。%%此接口功能将fmincon输入转换为所需格式%通过XYZ解算器，并将结果转换回fmincon输出。%%投入% FUN:非线性成本。FUN接受输入z(决策变量)和%返回成本f（在z处计算的标量值）及其%梯度g（在z处计算的nz-x-1矢量）。%z0：z的初始猜测%A，b:A*z<=b%Aeq，beq:Aeq*z==beq%lb，ub:z的上下限% NONLINCON:非线性约束。NONLINCON接受输入z并返回%向量C和Ceq作为前两个输出，表示%非线性不等式和等式，其中当C(z) <= 0和Ceq(z) = 0时，% FUN最小。% NONLINCON也返回C的雅可比矩阵(一个nz-by-ncin%稀疏矩阵）和Ceq（nceq稀疏矩阵）。%%输出%zopt：z的最优解%成本:在最优z下的最优成本%exitflag：%1满足一阶最优性条件。%0太多的函数求值或迭代。%-1由输出/绘图功能停止。%-2未找到可行点。%%您可以使用这个示例作为模板来实现接口文件%您自己的NLP解算器。解算器必须是计算机上的M文件或MEX文件% MATLAB路径。%%不编辑上面的行。%%版权所有2018 MathWorks，Inc。%%设置线性和非线性约束的尺寸信息num_lin_ineq=大小（A，1）；num_lin_eq=大小（Aeq，1）；[in，eq]=NONLINCON（z0）；num_non_ineq=长度（英寸）；num_non_eq=长度（eq）；总计=num_non_ineq+num_non_eq+num_lin_ineq+num_lin_eq；logicals_nlineq=假（总计1）；逻辑线性（1:num非线性）=真；logicals_nleq=假（总计1）；逻辑逻辑（num_non_ineq+（1:num_non_eq））=true；logicals_ineq=假（总计1）；逻辑逻辑（num_non_ineq+num_non_eq+（1:num_lin_ineq））=true；logicals_eq=false（总计1）；逻辑均衡（num_non_ineq+num_non_eq+num_lin_ineq+（1:num_lin_eq））=true；选项=结构(“nlineq”，逻辑线，“nleq”logicals_nleq,．..“ineq”，logicals_ineq，“情商”，逻辑（eq）；%%设置决策变量界限options.lb=lb；options.ub=ub；%%设置非线性约束的RHS选项。cl =[无穷(num_non_ineq 1); 0 (num_non_eq 1)];选项。铜= [0 (num_non_ineq 1); 0 (num_non_eq 1)];%%设置线性约束的RHSoptions.rl=[-inf（num_lin_ineq，1）；beq]；options.ru=[b；beq]；%%集合A矩阵选项A=稀疏（[A；Aeq]）；%%设置XYZ求解器选项options.algorithm=struct(“打印级别”，0，“max_iter”,200,．..“最大cpu时间”, 1000,“托尔”，1.0000e-06，．..“hessian_approximation”,内存有限的);%%设置XYZ求解器使用的函数句柄Jstr =稀疏的(num_non_ineq + num_non_eq,长度(z0)));func =结构(“目标”，@（x）fval（FUN，x），．..“梯度”，@（x）gval（FUN，x），．..“约束”@ (x) conval (NONLINCON x),．..“雅可比”，@（x）jacval（NONLINCON，x），．..“雅可比结构”Jstr @ ()．..);%%呼叫XYZ并返回成本和状态[zopt，output]=XYZsolver（z0，funcs，options）；成本=乐趣（zopt）；exitflag=convertStatustoExitflag（output.status）；%%效用函数作用f=fval（乐趣，z）%收益非线性成本[f，~]=fun（z）；作用g = gval(有趣,z)%收益成本梯度(~ g) =乐趣(z);作用c = conval (nonlcon, z)%返回非线性约束(eq) = nonlcon (z);c = (, eq);作用J = jacval (nonlcon, z)在稀疏矩阵中返回约束雅可比矩阵为nc-by-nz%Jin是nz的ncin稀疏，Jeq是nz的nceq稀疏[~, ~,金,Jeq] = nonlcon (z);J = ';作用exitflag=convertStatustoExitflag（状态）转换(状态)情况下0%的信息。status = 'Success';exitflag = 1;情况下1.%信息状态='已解决至可接受水平'；exitflag = 1;情况下2.%信息状态='不可行'；exitflag=-1；情况下3.%的信息。status = 'Search Direction Becomes Too Small';exitflag=-2；情况下4.%info.Status='发散迭代'；exitflag=-2；情况下6.%的信息。status = 'Feasible Point Found';exitflag = 1;情况下-1%info.Status='超过迭代次数'；exitflag = 0;情况下-4%信息状态='超过最大时间'；exitflag = 0;否则%的信息。status = 'IPOPT Error';exitflag=-2；终止
您可以使用此文件作为模板，将接口文件实现到您自己的NLP解算器。解算器必须是MATLAB路径上的MATLAB脚本或MEX文件。
控件中的自定义求解器可以插入该求解器nlmpc对象
nlobj.Optimization.CustomSolverFcn = @ResistantTBSolver;
只要“XYZ”解算器可靠且其选项选择正确，重新运行模拟应该会产生类似的结果。
参考文献
荣格，E.， S. Lenhart, Z. Feng。“双菌株结核模型治疗的最佳控制”。离散和连续动力系统，B2系列，2002年，第479-482页。
另见
相关话题
您还可以从以下列表中选择网站：
