使用自定义求解器的非线性MPC优化结核病治疗
这个例子展示了如何通过构造一个非线性MPC问题来寻找治疗双菌株结核(TB)的最优策略。然后非线性MPC控制器使用默认求解器和自定义求解器来计算最优解。
双菌株结核模型
在[1]中引入了一个双菌株结核模型。模型描述如下:
“在没有有效的疫苗的情况下，TB的当前控制程序专注于化疗。活性TB（用药物敏感菌株）患者的抗生素治疗需要更长的时间和更高的成本，而不是那些感染敏感的TB，但尚未开发这种疾病。缺乏抑制药物治疗不仅可能导致复发，而且对抗生素抗性TB的发育 - 今天面临的社会最严重的公共卫生问题之一。案件的减少可以通过“壳体持有”来实现药物敏感的TB，这是指用于确保药物摄入量规律的活动和技术，以获得治愈的持续时间，或者通过“案例发现”是指识别（通过例如，筛选的个体潜伏期患有敏感的TB，患有高风险的疾病，可能会受益于第一代码BLOC的预防性干预描述k。这些预防治疗将减少药物敏感TB的发病率（每单位时间新的案例），因此间接降低耐药性TB的发生率。“
在该示例中使用的动态模型中，总宿主N(这是一个常量)被分为六个不同的流行病学类别。其中五个类被定义为状态变量:
x (1) =
S.，易感个体数量x（2）=
T.，有效治疗的人数x（3）=
L2感染了耐药结核病，但不具传染性x (4) =
I1是典型结核的传染性疾病x (5) =
I2，具有抗性TB的感染性
第六堂课，L1(潜伏期，感染典型结核但不具传染性N - (S + T + L2 + I1和I2) +．
您可以使用两个可操纵变量(MVs)减少耐药结核病病例:
U（1） - “案例发现”，相对廉价的努力占识别需要治疗的努力。
U（2） - “案例持有”，相对昂贵的努力维持有效治疗。
有关动态模型的更多信息，请参见ResistantTBStateFcn.m．
非线性控制问题
结核病治疗的目标是在五年内减少潜伏(L2)和传染性(I2)耐药结核病患者，同时保持低成本。为了实现这一目标，使用一个成本函数，该函数在五年内将以下值相加。
F = L2 + I2 + 0.5*B1*u1^2 + 0.5*B2*u2^2
在这里，重量B1.是50，重量b2是500.．这些权重强调，由于对成本的影响，发现案件比保留案件更优先。
总数为N=30000．在初始条件下:
S.= 19000T.= 250.L2= 500I1= 1000I2= 250.
哪个叶子L1= 9000。
N = 30000;x0 = [76;1;2;4;1) * N / 120;
在本例中，利用非线性MPC控制器找出最优控制策略。创建nlmpc对象，具有正确的状态、输出和输入数量。
nx = 5;纽约= nx;ν= 2;nlobj = nlmpc (nx、纽约、ν);
在标准成本函数中，一个或多个OVs默认为零权值，因为mv比OVs少。
假设治疗政策只能每三个月调整一次。因此，将控制器采样时间设置为0.25年。因为您想要找到5年的最优策略，所以将预测范围设置为20步(5年除以0.25)。
年= 5;t = 0.25;p =年/ Ts;nlobj。T.s = 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．将这些值设置为MV的较低和上限。
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 =“抵抗bcostfcn”;nlobj.Optimization.ReplaceStandardCost = true;
为加快仿真速度，给出了代价的解析雅可比矩阵抵制抵抗．
nlobj.jacobian.customcostfcn =.“ResistantTBCostJacFcn”;
自L1总体定义为N减去所有国家的总和，必须确保(s + t + l2 + i1 + i2) - n < 0总是满意。在nlmpc对象，则使用匿名函数将此条件指定为不等式约束。
nlobj. optimizer . customineqconfcn = @(X,U,e,data) sum(X(2:end，:)，2)-30000;
要检查潜在的数值问题，请使用validateFcns命令。
validateFcns (nlobj x0, [0.5, 0.5])
model.statefcn是好的。jacobian.statefcn还可以。没有指定输出函数。假设预测模型中的“y = x”。优化.customcostfcn是可以的。jacobian.customcostfcn还可以。优化.customineqconfcn是可以的。分析用户提供的模型，成本和约束函数完成。
计算最优控制策略，使用nlmpcmove.函数。在初始条件下，mv是零。默认情况下,fmincon使用优化工具箱™中的NLP求解器作为默认的NLP求解器。
lastMV = 0(ν,1);[~, ~,信息]= nlmpcmove (x0, nlobj lastMV);
闲置变量未使用或零加权您的定制成本函数。所有的约束都将是艰难的。
绘制并检查最优解决方案。
抵抗力（Info，TS）
最佳成本= 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．
函数[zopt，cost，ExitFlag] =抵制（乐趣，Z0，A，B，AEQ，BEQ，LB，UB，非林肯）这是一个调用XYZ非线性编程的接口函数％求解器解决NLMPC控制器定义的NLP问题。这%这个函数的输入和输出定义与fmincon相同。％％此接口功能将Fmincon输入转换为所需的格式XYZ求解器的％并将结果转换回Fmincon输出。％%的输入% FUN:非线性成本。FUN接受输入z(决策变量)和%返回成本f(在z处计算的标量值)及其％梯度g（在z时评估的NZ-BY-1载体）。% z0: z的初始猜想％a，b：a * z <= b% Aeq,说真的:Aeq * z = =说真的% lb,ub: z的下界和上界% NONLINCON:非线性约束。NONLINCON接受输入z并返回％vector c和ceq作为前两个输出，代表％非线性不等式和平等分别在哪里当C(z) <= 0和Ceq(z) = 0时，% FUN被最小化。% NONLINCON也返回C (a nz-by-ncin)的雅可比矩阵%稀疏矩阵)和Ceq (nz-by-nceq稀疏矩阵)。％%输出% zopt: z的最优解%成本:在最优z处的最优成本% exitflag:% 1满足一阶最优条件。% 0太多的函数计算或迭代。% -1由输出/绘图函数停止。% -2没有找到可行点。％%你可以使用这个例子作为模板来实现一个接口文件%你自己的NLP求解器。解算器必须是M文件或MEX文件％matlab路径。％%不要编辑上面的行。％The MathWorks, Inc.版权所有%%设置线性和非线性约束的维度信息num_lin_ineq =大小(1);num_lin_eq =大小(Aeq, 1);(eq) = NONLINCON (z0);num_non_ineq =长度();num_non_eq =长度(eq);Total = 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];%%设置A矩阵选项。一种=sparse([A; Aeq]);%%设置XYZ解算器optonsoptions.algorithm = struct（“print_level”0,“max_iter”，200，．．．“max_cpu_time”, 1000,'tol'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)% Return非线性约束(eq) = nonlcon (z);c = (, eq);函数J = jacval (nonlcon, z)％返回约束Jacobian在稀疏矩阵中为nc-by-nz％金是nz-by-ncin稀疏，jeq是nz-by-nceq稀疏[〜，〜，金，jeq] = nonlcon（z）;j = [金杰]';函数ExitFlag = ConvertStatustoExitFlag（状态）开关(状态)情况下0.%的信息。S.tatus = 'Success';exitflag = 1;情况下1％info.status ='解决了可接受的水平';exitflag = 1;情况下2％info.status ='不可行的';ExitFlag = -1;情况下3.%的信息。S.tatus = 'Search Direction Becomes Too Small';ExitFlag = -2;情况下4.％info.status ='发散迭代';ExitFlag = -2;情况下6.%的信息。S.tatus = 'Feasible Point Found';exitflag = 1;情况下－１%的信息。S.tatus = 'Exceeded Iterations';exitflag = 0;情况下－4%的信息。S.tatus = 'Max Time Exceeded';exitflag = 0;否则%的信息。S.tatus = 'IPOPT Error';ExitFlag = -2;结尾
您可以将此文件用作模板以将接口文件实现到您自己的NLP求解器。求解器必须是MATLAB路径上的MATLAB脚本或MEX文件。
您可以通过将其指定为自定义求解器来插入求解器nlmpc对象。
nlobj.Optimization.CustomSolverFcn = @ResistantTBSolver;
只要“XYZ”求解器是可靠的，并且选择其选项都被正确选择，重新运行模拟应该产生类似的结果。
参考
荣格(Jung, E. Lenhart, S. Lenhart, Z. Feng);双菌株结核模型治疗的最优控制离散和连续动力系统，B2,2002系列，第479-482号。
也可以看看
相关的话题
你也可以从以下列表中选择一个网站:
