并行优化ODE

本示例展示了如何优化ODE的参数。

它还展示了当ODE解返回目标函数和非线性约束函数时，如何避免计算两次。这个例子比较patternsearch而且遗传算法在运行求解器的时间和解的质量方面。万博尤文图斯

您需要并行计算工具箱™许可证才能使用并行计算。

步骤1。定义问题。

问题是改变大炮的位置和角度，使炮弹尽可能地发射到墙外。加农炮初速300米/秒。这面墙有20米高。如果大炮离墙太近，它就必须以太陡的角度发射，而炮弹的飞行距离就不够远。如果大炮离墙太远，弹丸也不能飞得足够远。

空气阻力使弹丸变慢。阻力与速度的平方成正比，比例常数为0.02。重力作用于弹丸，使其以恒定的9.81米/秒向下加速²．因此，轨迹的运动方程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))．因此，优化问题只依赖于两个标量，因此它是一个二维问题。以水平距离和角度作为决策变量。

步骤2。建立ODE模型。

ODE求解器要求您将模型表述为一阶系统。增大轨迹向量(x₁（t)，x₂（t))及其时间导数(x＇₁（t)，x＇₂（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 =炮灰(t,x) F = [x(3);x(4);x(3);x(4)];% Initial，使f(1)和f(2)正确NRM = norm(x(3:4)) * .02;%速度模乘以常数F (3) = -x(3)*nrm;水平加速度%F (4) = -x(4)*nrm - 9.81;%垂直加速度

可视化的解的ODE开始30米从墙壁的角度π/ 3．

生成图形的代码

步骤3。使用patternsearch解决。

问题是求初始位置x0 (1)初始角x0 (2)为了使弹丸与墙的距离最大化。因为这是一个最大化问题，所以最小化距离的负值(参见最大化与最小化）.

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

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

函数f = cannonobjective x0 (x) = (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 = deval(sol,zerofnd,1);F = -f;%取距离的负数为最大值函数[c,ceq] = cannonconstraint(x) ceq = [];X0 = [x(1);0;300*cos(x(2));300*sin(x(2))];Sol = ode45(@cannonfodder，[0,15]，x0);如果Sol.y (1,end) <= 0投掷物从未击中墙壁C = 20 - sol.y(2,end);其他的找到当弹丸穿过x = 0时zerofnd = fzero (@ (r)德瓦尔(溶胶,r, 1), [sol.x (2), sol.x(结束)]);然后求出这里的高，再减去20C = 20 - deval(sol,zerofnd,2);结束

注意，目标函数和约束函数设置了它们的输入变量x0到ODE求解器的4-D初始点。如果弹丸击中墙壁，ODE求解器不会停止。相反，约束函数只是变成正的，表明一个不可行的初始值。

初始位置x0 (1)不能大于0，小于-200是徒劳的。(它应该在-20附近，因为如果没有空气阻力，最长的轨迹将以一个角度从-20开始π/ 4)。同样，初始角x0 (2)不能小于0，也不能大于0π/ 2．设置边界稍微远离这些初始值:

Lb = [-200;0.05];Ub = [-1;pi/2-.05];X0 = [-30,pi/3];%初始猜测

设置UseCompletePoll选项真正的．这提供了一个更高质量的解决方案，并支持与并行处理进行直接比较，因为并行计算需要这种设置。

Opts = optimoptions(“patternsearch”，“UseCompletePoll”,真正的);

调用patternsearch解决问题。

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

优化终止:网格尺寸小于选项。MeshTolerance和约束违反小于options. constraintolerance。xsolution = -28.8123 0.6095 distance = -125.9880 eflag = 1 outpt = struct with fields: function: @cannonobjective problemtype: '非线性constr' pollmethod: 'gpspositivebasis2n' maxconstraint: 0 searchmethod: [] iterations: 5 funccount: 269 meshsize: 8.9125 -07 rngstate: [1×1 struct] message: '优化终止:网格尺寸小于选项。MeshTolerance(和约束违反小于options. constraintolerance。)运行时间为0.792152秒。

以0.6095弧度的角度从距离墙壁约29米的地方发射，可以获得最远的距离，约126米。报告的距离是负数，因为目标函数是到墙的距离的负数。

可视化解决方案。

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);%等于plot的时间Xs = deval(sol,t,1);%插值的x值Ys = deval(sol,t,2);%插值的y值情节(x, y)在情节([0,0],[0,20],“k”）%画墙包含(的水平距离) ylabel (的轨道高度)传说(“轨迹”，“墙”，“位置”，“西北”) ylim([0 70])保持从

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

目标函数和非线性约束函数都调用ODE求解器来计算它们的值。在同一函数中的客观约束与非线性约束避免两次调用求解器。的runcannonFile实现了这项技术。将此文件复制到MATLAB路径下的文件夹中。

函数[x,f,eflag,outpt] = runcannon(x0,opts)如果Nargin == 1%未提供选项Opts = [];结束xLast = [];调用了ode求解器的最后一个位置Sol = [];ODE解结构Fun = @objfun;目标函数，嵌套在下面Cfun = @constr;%约束函数，嵌套在下面Lb = [-200;0.05];Ub = [-1;pi/2-.05];调用patternsearch[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(@cannonfodder，[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(@cannonfodder，[0,15]，x0);xLast = x;结束现在计算约束函数首先找到当弹丸穿过x = 0时zerofnd = fzero (@ (r)德瓦尔(溶胶,r, 1), [sol.x (1) sol.x(结束)]);然后求出这里的高，再减去20C = 20 - deval(sol,zerofnd,2);结束结束

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

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

优化终止:网格尺寸小于选项。MeshTolerance和约束违反小于options. constraintolerance。运行时间为0.670715秒。

求解器比以前跑得快了。如果检查解决方案，您会发现输出是相同的。

第5步。并行计算。

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

parpool

使用“本地”配置文件启动parpool…连接到并行池(number of workers: 6). ans = ProcessPool with properties: Connected: true NumWorkers: 6 Cluster: local AttachedFiles: {} AutoAddClientPath: true IdleTimeout: 30 minutes(剩余30分钟)SpmdEnabled: true

将选项设置为使用并行计算，并重新运行求解器。

Opts = optimoptions(“patternsearch”选择,“UseParallel”,真正的);X0 = [-30;pi/3];Tic [xsolution,distance,eflag,outpt] = runcannon(x0,opts);toc

优化终止:网格尺寸小于选项。MeshTolerance和约束违反小于options. constraintolerance。运行时间为1.894515秒。

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

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

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

如在同一函数中的客观约束与非线性约束，使用时不能节省时间遗传算法所采用的嵌套函数技术patternsearch步骤4。相反,叫遗传算法同时通过设置适当的选项。

选项= optimoptions(“遗传算法”，“UseParallel”,真正的);rng默认的%用于再现性抽搐%溶液时间[xsolution,distance,eflag,outpt] = ga(@cannonobjective,2，.．.[],[],[],[], 磅,乌兰巴托,toc @cannonconstraint、期权)

优化终止:适应度值的平均变化小于选项。FunctionTolerance和constraint违例要小于options. constraintolerance。xsolution = -37.6217 0.4926 distance = -122.2181 eflag = 1 outpt = struct with fields: problemtype: '非线性constr' rngstate: [1×1 struct] generations: 4 funccount: 9874 message: '优化终止:健身值的平均变化小于选项。FunctionTolerance(和约束违反小于options. constraintolerance)。' maxconstraint: 0 hybridflag:[]运行时间为12.529131秒。

的遗传算法解决办法不如解决办法patternsearch解决方案:122米对126米。遗传算法需要更多时间:大约12秒vs不到2秒;patternsearch串行和嵌套的时间更短，不到1秒。运行遗传算法串行运行需要更长的时间，一次测试运行大约30秒。