并行优化ODE

这个例子展示了如何优化ODE的参数。

当ODE解返回目标函数和非线性约束函数时，如何避免两次计算。这个例子比较patternsearch和遗传算法在运行解算器的时间和解决方案的质量方面。金宝搏官方网站

你需要一个并行计算工具箱™ 使用并行计算的许可证。

步骤1.定义问题。

问题是改变加农炮的位置和角度，以使射弹尽可能远离墙。加农炮的初速为300 m/s。墙高20 m。如果加农炮离墙太近，则必须以太陡的角度发射，射弹的射程不够远。如果加农炮离墙太远，则t他也走得不够远。

空气阻力使弹丸减速。阻力与速度的平方成正比，比例常数为0.02。重力作用于弹丸，使其向下加速，速度常数为9.81 m/s^2.．因此，轨迹的运动方程x(T)

$\frac{D^{2.} x (T)}{D T^{2.}} = - 0.02 ‖ v (T) ‖ v (T) - (0, 9.81),$

哪里 $v (T) = D x (T) / D T ．$

初始位置x0初速xp0是二维向量。但是，初始高度x0（2）为0，因此初始位置仅取决于标量x0（1）. 和初始速度xp0大小为300(枪口初速)，所以只取决于初始角度，一个标量。对于初始角度th,xp0=300*（cos（th）、sin（th））因此，优化问题只依赖于两个标量，因此它是一个二维问题。使用水平距离和角度作为决策变量。

第二步。制定ODE模型。

ODE求解器要求您将模型表述为一阶系统。增加轨迹矢量(x_1.(T)，x_2.(T)的时间导数(x'_1.(T)，x'_2.(T))形成一个四维轨迹向量。根据这个增广向量，微分方程变成

$\frac{D}{D T} x (T) = [\begin{matrix} x_{3.} (T) \\ x_{4.} (T) \\ - .02 ‖ (x_{3.} (T), x_{4.} (T)) ‖ x_{3.} (T) \\ - .02 ‖ (x_{3.} (T), x_{4.} (T)) ‖ x_{4.} (T) - 9.81 \end{matrix}] ．$

将微分方程写成函数文件，并保存在MATLAB中^®路径。

函数F = [x(3);x(4);x(4)];%初始值，使f（1）和f（2）正确NRM = norm(x(3:4)) * .02;%速度乘以常数的范数f (3) = - x(3) *全国抵抗运动;%水平加速度f（4）=-x（4）*nrm-9.81；%的垂直加速度

想象ODE的解决方案，从墙30 m处开始，角度为pi/3．

生成图形的代码

步骤3。使用patternsearch解决。

问题是找到初始位置x0（1）和初始角x0（2）为了最大限度地增加射弹与墙壁的距离。因为这是一个最大化问题，所以最小化距离的负数（请参见最大化和最小化)．

使用patternsearch要解决这个问题，您必须提供目标、约束、初始猜测和选项。

这两个文件是目标函数和约束函数。将它们复制到MATLAB路径上的文件夹中。

函数f=cannonobjective（x）x0=[x（1）；0；300*cos（x（2））；300*sin（x（2））]；sol=ode45（@cannonfodder[0,15]，x0）；%求y_2（t）=0时的时间tzerofnd = fzero (@ (r)德瓦尔(sol r 2), [sol.x (2), sol.x(结束)]);%然后找到当时的x位置f =德瓦尔(sol zerofnd 1);f = - f;%取距离的负数表示最大化函数[c,ceq] = cannonconstraint(x) ceq = [];x0 = [x (1), 0; 300 * cos (x (2)); 300 * sin (x (2)));索尔=数值(@cannonfodder [0, 15], x0);如果sol.y (1) < = 0%射弹永远达不到墙c=20-溶胶y（2，结束）；其他的求弹丸穿过x = 0的时间zerofnd=fzero（@（r）deval（sol，r，1），[sol.x（2），sol.x（end）]；然后求出这里的高度，用20减去c=20-德瓦尔（索尔，零FND，2）；结束

请注意，目标函数和约束函数设置了它们的输入变量x0到ODE解算器的4-D初始点。如果投射物撞击墙壁，ODE解算器不会停止。相反，约束函数只是变为正，表示初始值不可行。

初始位置x0（1）不能高于0，低于-200是徒劳的。（它应该接近-20，因为在没有空气阻力的情况下，最长的弹道将以-20的角度开始。）pi/4．)同样的，初始角度x0（2）不能低于0，也不能高于pi/2.将边界设置为稍微远离这些初始值：

lb=[-200；0.05]；ub=[-1；pi/2-.05]；x0 =(-30,π/ 3);%初始猜测

设置使用CompletePoll选择符合事实的．这提供了更高质量的解决方案，并支持与并行处理的直接比较，因为并行计算需要此设置。

opts=options(“patternsearch”,“UseCompletePoll”,真正的);

呼叫patternsearch来解决这个问题。

抽搐%时间溶液[xsolution，distance，eflag，outpt]=patternsearch（@cannonobjective，x0，...[],[],[],[], 磅,乌兰巴托,@cannonconstraint toc选择)

优化终止:网格尺寸小于选项。MeshTolerance和constraint violation小于options. constraintolerance。xsolution = -28.8123 0.6095 distance = -125.9880 eflag = 1 outpt = function: @cannonobjective problemtype: 'nonlinearconstr' pollmethod: ' gppositivebasis2n ' maxconstraint: 0 searchmethod: [] iterations: 5 funccount: 269 meshsize: 8.9125e-07 rngstate: [1x1 struct] message: '优化终止:meshsize小于选项。运行时间为3.174088秒。

以0.6095弧度的角度从距离墙壁约29 m的位置开始发射，会产生最远的距离，约126 m。报告的距离为负值，因为目标函数为距离墙壁的负值。

将解决方案可视化。

x0=[xsolution（1）；0；300*cos（xsolution（2））；300*sin（xsolution（2））]；sol=ode45（@cannonfodder[0,15]，x0）；%找出射弹落地的时间zerofnd = fzero (@ (r)德瓦尔(sol r 2), [sol.x (2), sol.x(结束)]);t = linspace (0, zerofnd);%等积时间x =德瓦尔(sol t 1);%插值x值ys=德瓦尔（索尔，t，2）；%插值y值情节(x, y)在…上绘图（[0,0]，[0,20]，“k”)%画墙包含(的水平距离)伊拉贝尔(“弹道高度”)传奇(“轨迹”,“墙”,“位置”,“西北”([0 70])保持从

步骤4.避免调用昂贵的子例程两次。

目标约束函数和非线性约束函数都调用ODE解算器来计算它们的值。在中使用该技术同一函数中的目标约束和非线性约束（优化工具箱）以避免调用求解器两次runcannon文件实现了这种技术。复制此文件到您的MATLAB路径上的一个文件夹。

函数[x，f，eflag，outpt]=runcannon（x0，opts）如果输入参数个数= = 1%没有提供选项opts=[]；结束xLast=[]；%上次调用ode解算器的位置溶胶=[]；%常微分方程解结构有趣= @objfun;%目标函数，嵌套在下面cfun = @constr;%约束函数，嵌套在下面lb=[-200；0.05]；ub=[-1；pi/2-.05]；%调用模式搜索[x, f, eflag, outpt] = patternsearch (x0有趣, ,[],[],[],[], 磅,乌兰巴托,cfun选择);函数y=objfun（x）如果~ isequal (x, xLast)%检查是否需要计算x0=[x（1）；0；300*cos（x（2））；300*sin（x（2））]；sol=ode45（@炮灰，[0,15]，x0）；xLast=x；结束%现在计算目标函数%第一次发现时，弹丸击中地面zerofnd = fzero (@ (r)德瓦尔(sol r 2), [sol.x (2), sol.x(结束)]);%然后计算当时的x位置y=deval（sol，zerofnd，1）；y=-y；取距离的负数结束函数[c，ceq]=constr（x）ceq=[]；如果~ isequal (x, xLast)%检查是否需要计算x0=[x（1）；0；300*cos（x（2））；300*sin（x（2））]；sol=ode45（@炮灰，[0,15]，x0）；xLast=x；结束%现在计算约束函数首先求弹丸什么时候越过x = 0zerofnd=fzero（@（r）deval（sol，r，1），[sol.x（1），sol.x（end）]；然后求出这里的高度，用20减去c=20-德瓦尔（索尔，零FND，2）；结束结束

重新初始化问题并将调用计时到runcannon．

x0 =(-30;π/ 3);Tic [xsolution,distance,eflag,outpt] = runcannon(x0,opts);toc

运行时间为2.610590秒。

解算器比以前跑得快。如果您检查这个解决方案，您会看到输出是相同的。

第5步。并行计算。

尝试通过并行计算来节省更多时间。首先打开一个并行的工作人员池。

帕尔普

正在使用“本地”配置文件启动parpool。。。连接到4名工人。ans=具有以下属性的池：已连接：true NumWorkers:4群集：本地附件文件：{}IdleTimeout:30分钟（剩余30分钟）SpmdEnabled:true

设置选项以使用并行计算，然后重新运行解算器。

opts=options(“patternsearch”，选择，“使用并行”,真正的);x0 =(-30;π/ 3);Tic [xsolution,distance,eflag,outpt] = runcannon(x0,opts);toc

运行时间为3.917971秒。

在这种情况下，并行计算速度较慢。如果检查解决方案，您会发现输出是相同的。

步骤6.与遗传算法进行比较。

你也可以尝试用遗传算法来解决这个问题。然而，遗传算法通常速度较慢，可靠性较低。

这个朗卡纳加文件调用遗传算法，避免ODE求解器的双重求值。它就像runcannon，但是遗传算法而不是patternsearch，并检查轨迹是否到达壁面。复制此文件到您的MATLAB路径上的一个文件夹。

函数[x，f，eflag，outpt]=runcannonga（opts）如果输入参数个数= = 1%没有提供选项opts=[]；结束xLast=[]；%上次调用ode解算器的位置溶胶=[]；%常微分方程解结构有趣= @objfun;%目标函数，嵌套在下面cfun = @constr;%约束函数，嵌套在下面lb=[-200；0.05]；ub=[-1；pi/2-.05]；%叫ga[x，f，eflag，outpt]=ga（fun，2，[]，[]，[]，[]，[]，[]，[]，[]，[]，[]，lb，ub，cfun，opts）；函数y=objfun（x）如果~ isequal (x, xLast)%检查是否需要计算x0=[x（1）；0；300*cos（x（2））；300*sin（x（2））]；sol=ode45（@炮灰，[0,15]，x0）；xLast=x；结束%现在计算目标函数%第一次发现时，弹丸击中地面zerofnd = fzero (@ (r)德瓦尔(sol r 2), [sol.x (2), sol.x(结束)]);%然后计算当时的x位置y=deval（sol，zerofnd，1）；y=-y；取距离的负数结束函数[c，ceq]=constr（x）ceq=[]；如果~ isequal (x, xLast)%检查是否需要计算x0=[x（1）；0；300*cos（x（2））；300*sin（x（2））]；sol=ode45（@炮灰，[0,15]，x0）；xLast=x；结束%现在计算约束函数如果sol.y (1) < = 0%射弹永远达不到墙c=20-溶胶y（2，结束）；其他的求弹丸穿过x = 0的时间zerofnd=fzero（@（r）deval（sol，r，1），[sol.x（2），sol.x（end）]；然后求出这里的高度，用20减去c=20-德瓦尔（索尔，零FND，2）；结束结束结束

呼叫朗卡纳加同时。

opts=optimoptions（'ga'，'UseParallel'，true）；rng再现性的默认百分比[X解、距离、eflag、输出]=runcannonga（opts）toc

优化终止:适应度值的平均变化小于选项。函数容忍度和约束违背小于options. constraintolerance。xsolution = -17.9172 0.8417 distance = -116.6263 eflag = 1 outpt = problemtype: 'nonlinearconstr' rngstate: [1x1 struct] generations: 5 funccount: 20212 message: [1x140 char] maxconstraint: 0运行时间是119.630284秒。

这个遗传算法解决方案不如解决方案好patternsearch答案:117米对126米。遗传算法需要更多的时间：大约120秒，而不是5秒以下。