主要内容

并行优化ODE

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

当ODE解返回目标函数和非线性约束函数时,如何避免两次计算。这个例子比较PatternSearch.GA.运行求解器的时间和解决方案的质量。金宝搏官方网站

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

步骤1。定义问题。

问题是改变火炮的位置和角度,让炮弹尽可能地飞出墙外。炮口初速为300米/秒。墙有20米高。如果加农炮离墙太近,就必须以太陡的角度发射,射弹就不能飞得足够远。如果加农炮离墙太远,炮弹也不能飞得足够远。

空气阻力使炮弹减速。阻力与速度的平方成正比,比例常数为0.02。重力作用于抛射体,使其以恒定的9.81 m/s向下加速2.因此,轨迹的运动方程xt)

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(枪口速度),因此仅取决于初始角度,标量。初始角度thxp0300 * (cos (th), sin (th)).因此,优化问题只依赖于两个标量,因此是一个二维问题。使用水平距离和角度作为决策变量。

步骤2。建立ODE模型。

ODE求解器要求您将您的模型作为一阶系统制定。增强轨迹矢量(x1t),x2t))及时衍生物(x1t),x2t)以形成一个4-D轨迹向量。用增广向量表示,微分方程就变成

d d t x t 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

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

函数f = cannonfodder(t,x)f = [x(3); x(4); x(3); 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米的角度出发π/ 3

生成图形的代码

步骤3。使用patternsearch解决。

问题是找到初始位置x0 (1)和初始角度x0 (2)为了最大限度地提高离墙的距离,炮弹落地。因为这是一个最大化的问题,最小化距离的负(见最大化与最小化)。

使用PatternSearch.要解决此问题,您必须提供目标,约束,初始猜测和选项。

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

函数f = cannonobjective x0 (x) = (x (1); 0; 300 * cos (x (2)); 300 * sin (x (2)));索尔=数值(@cannonfodder [0, 15], x0);%求y_2(t) = 0时的时间tzerofnd = fzero(@(r)deval(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)));索尔=数值(@cannonfodder [0, 15], x0);如果sol.y(1,结束)<= 0投射物永远达不到墙C = 20 - sol.y(2,结束);别的当射弹十字架x = 0时发现%zerofnd = fzero (@ (r)德瓦尔(溶胶,r, 1), [sol.x (2), sol.x(结束)]);%然后在那里找到高度,并从20中减去C = 20 - deval(sol,zerofnd,2);结尾

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

初始位置x0 (1)不能高于0,并且它是徒劳的,它低于-200。(它应该在-20附近,因为没有空气阻力,最长的轨迹将以-20开始以一定角度开始π/ 4。)类似地,初始角度x0 (2)不能低于0,也不能高于0π/ 2.设置边界稍微远离这些初始值:

磅= (-200;0.05);乌兰巴托=(1;π/ 2 . 05);x0 =(-30,π/ 3);%初始猜测

设定UseCompletePoll选项真正的.这提供了更高质量的解决方案,并实现了与并行处理的直接比较,因为并行计算需要此设置。

选择= optimoptions ('patternsearch'“UseCompletePoll”,真正的);

调用PatternSearch.解决问题。

抽搐%时间解决方案[xsolution、距离、eflag, outpt] = patternsearch (x0, @cannonobjective...[],[],[],[],LB,UB,@ Cannonstraint,OPTS)TOC
优化终止:少于options.meshtolerance和约束违规的网格尺寸小于选项.ConstraintTolerance。xsolution = -28.8123 0.6095距离= -125.9880 eflag = 1 outpt = struct with字段:@cannonobjective问题rype:'nonlinearconstr'pollmethod:'gpspositivebasis2n'maxconstraint:0 searchmethod:[]迭代:5 funccount:269网格化:8.9125e-07 rngstate:[1×1结构]消息:'优化终止:网格尺寸小于options.meshtolance‖和约束违规小于选项。“constrainttolerance。经过时间为0.792152秒。

以0.6095弧度的角度从距离壁面约29米的地方开始发射,结果产生了最远的距离,约126米。报告的距离是负的,因为目标函数是到墙的距离的负。

可视化的解决方案。

x0 = [xsolution(1); 0; 300 * cos(xsolution(2)); 300 * sin(xsolution(2))];索尔=数值(@cannonfodder [0, 15], x0);%找到射弹土地的时间zerofnd = fzero(@(r)deval(sol,r,2),[sol.x(2),sol.x(端)]);t = linspace(0,zerofnd);情节的%相同时间x =德瓦尔(sol t 1);%插值x值y =德瓦尔(sol t 2);%内插y值绘图(xs,ys)持有情节([0,0],[0,20],'K'%绘制墙壁包含('水平距离') ylabel (的轨道高度) 传奇(“轨迹”“墙”“位置”“西北”)ylim([070])保持离开

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

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

函数[x,f,eflag,outpt] = runcannon(x0,opts)如果输入参数个数= = 1%没有选择选择= [];结尾xLast = [];最后一位ode求解器被调用索尔= [];% ODE溶液结构有趣= @objfun;%目标函数,嵌套在下面cfun = @constr;%约束函数,嵌套在下面磅= (-200;0.05);乌兰巴托=(1;π/ 2 . 05);%叫patternsearch[x,f,Eflag,Outpt] = Patternsearch(Fun,X0,[],[],[],[],LB,UB,CFUN,OPTS);函数y = objfun(x)如果〜isequal(x,xlast)检查是否需要计算x0 = [x (1), 0; 300 * cos (x (2)); 300 * sin (x (2)));索尔=数值(@cannonfodder [0, 15], x0);xLast = x;结尾%现在计算目标函数首先找出弹丸什么时候击中地面zerofnd = fzero(@(r)deval(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)));索尔=数值(@cannonfodder [0, 15], x0);xLast = x;结尾%现在计算约束函数首先求弹丸什么时候越过x = 0zerofnd = fzero(@(r)deval(sol,r,1),[sol.x(1),sol.x(end)]);%然后在那里找到高度,并从20中减去C = 20 - deval(sol,zerofnd,2);结尾结尾

重新初始化问题并确定调用的时间runcannon

x0 = [-30; pi / 3];tic [xsolution,距离,eflag,outpt] = runcannon(x0,opts);TOC.
优化终止:少于options.meshtolerance和约束违规的网格尺寸小于选项.ConstraintTolerance。经过时间为0.670715秒。

求解器比以前更快。如果检查解决方案,则会看到输出相同。

步骤5.并行计算。

尝试通过并行计算更多时间。首先打开一个平行的工人池。

parpool
使用“本地”配置文件启动Parpool ...连接到并行池(工人数量:6)。ANS = ProcessPool具有属性:Connected:True NumWorkers:6群集:本地连接文件:{} autoaddclientPath:true idledimeout:30分钟(剩余30分钟)spmded:true

设置选项以使用并行计算,并重新运行求解器。

选择= optimoptions ('patternsearch'选择,“UseParallel”,真正的);x0 = [-30; pi / 3];tic [xsolution,距离,eflag,outpt] = runcannon(x0,opts);TOC.
优化终止:少于options.meshtolerance和约束违规的网格尺寸小于选项.ConstraintTolerance。运行时间为1.894515秒。

在这种情况下,并行计算速度较慢。如果检查解决方案,则会看到输出相同。

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

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

如上所述目标约束与非线性约束在同一函数中,使用时无法节省时间GA.通过使用的嵌套功能技术PatternSearch.在步骤4。相反,叫GA.同时设置适当的选项。

选择= optimoptions ('Ga'“UseParallel”,真正的);RNG.默认的%的再现性抽搐%时间解决方案[xsolution,距离,eflag,outpt] = ga(@ cannonobjective,2,...[],[],[],[], 磅,乌兰巴托,toc @cannonconstraint、期权)
优化终止:适应度值的平均变化小于选项。函数容忍度和约束违背小于options. constraintolerance。xsolution = -37.6217 0.4926 distance = -122.2181 eflag = 1 outpt = struct with fields: problemtype: 'nonlinearconstr' rngstate: [1×1 struct] generations: 4 funccount: 9874 message: '优化终止:健康值的平均变化小于选项。函数容忍↵和约束违背小于options. constraintolerance。' maxconstraint: 0 hybridflag: [] Elapsed time is 12.529131秒。

GA.解决之道不如解决之道PatternSearch.答案:122米对126米。GA.需要更多的时间:大约12岁,而不是2秒;PatternSearch.串行和嵌套的时间更少,不到1秒。跑步GA.连续运行的时间更长,一次测试大约需要30秒。

相关的例子

更多关于