并行优化ODE-MATLAB&Simulink-MathWorks한국万博1manbetx - 万博1manbetx,s manbetx 845,万博尤文图斯

并行优化ODE

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

当ODE解返回目标函数和非线性约束函数时，如何避免两次计算。这个例子比较模式搜索和ga运行求解器的时间和解决方案的质量。万博尤文图斯

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

步骤1。定义问题。

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

空气阻力使炮弹减速。阻力与速度的平方成正比，比例常数为0.02。重力作用于抛射体，使其以恒定的9.81 m/s向下加速²．因此，轨迹的运动方程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)以形成一个4-D轨迹向量。用增广向量表示，微分方程就变成

$\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）]；%初始值，使f（1）和f（2）正确nrm=标准（x（3:4））*.02；速度的范数乘以常数f (3) = - x(3) *全国抵抗运动;%的水平加速度f（4）=-x（4）*nrm-9.81；%的垂直加速度

设想ODE的解从离墙30米的角度出发π/ 3．

生成图形的代码

步骤3。使用patternsearch解决。

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

使用模式搜索要解决此问题，必须提供目标、约束、初始猜测和选项。

这两个文件是目标函数和约束函数。复制它们到您的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（end）]；然后求出该时刻的x坐标f=德瓦尔（索尔，零fnd，1）；f=-f；%取距离的负数以达到最大函数[c，ceq]=cannonconstraint（x）ceq=[]；x0=[x（1）；0；300*cos（x（2））；300*sin（x（2））]；sol=ode45（@cannonfoeder[0,15]，x0）；如果溶胶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 (“模式搜索”，“UseCompletePoll”,真正的);

调用模式搜索解决这个问题。

抽搐%确定解决方案的时间[xsolution，distance，eflag，outpt]=patternsearch（@cannonobjective，x0，．．．[]、[]、[]、[]、lb、ub、@cannonconstraint、opts）toc

优化已终止：网格大小小于options.MeshTolerance和约束冲突小于options.ConstraintTolerance。xsolution=-28.8123 0.6095 distance=-125.9880 eflag=1 outpt=struct带字段：function:@cannonobjective problemtype:'nonlinearconstr'pollmethod:'gpspospositivebasis2n'maxconstraint:0 searchmethod:[]迭代次数：5次函数计数：269网格大小：8.9125e-07 rngstate:[1×1结构]消息：“优化已终止：网格大小小于options.MeshTolerance↵ 约束冲突小于选项。ConstraintTolerance。'经过的时间为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）deval（sol，r，2），[sol.x（2），sol.x（end）]；t=linspace（0，zerofnd）；%等积时间x =德瓦尔(sol t 1);%插值x值y =德瓦尔(sol t 2);%插值y值绘图（xs，ys）保持在情节([0,0],[0,20],“k”）%画墙包含(“水平距离”) ylabel (的轨道高度)传奇(“轨迹”，“墙”，“位置”，“西北”)ylim（[0 70]）保持关

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

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

函数[x，f，eflag，outpt]=runcannon（x0，opts）如果纳金==1%没有提供选项选择= [];终止xLast = [];最后一位ode求解器被调用索尔= [];% ODE溶液结构fun=@objfun；%目标函数，嵌套在下面cfun=@constr；%约束函数，嵌套在下面磅= (-200;0.05);乌兰巴托=(1;π/ 2 . 05);%叫patternsearch[x，f，eflag，outpt]=模式搜索（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（end）]；%然后计算当时的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，distance，eflag，outpt]=runcannon（x0，opts）；toc

优化已终止：网格大小小于options.MeshTolerance和约束冲突小于options.ConstraintTolerance。运行时间为0.670715秒。

解算器比以前运行得更快。如果检查解决方案，您会发现输出是相同的。

第五步。并行计算。

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

parpool

正在使用“本地”配置文件启动parpool。。。已连接到并行池（工作线程数：6）。ans=ProcessPool，属性为：已连接：true NumWorkers:6群集：本地附件文件：{}AutoAddClientPath:true IdleTimeout:30分钟（剩余30分钟）SpmdEnabled:true

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

选择= optimoptions (“模式搜索”选择,“UseParallel”，对）；x0=[-30；pi/3]；tic[xsolution，distance，eflag，outpt]=runcannon（x0，opts）；toc

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

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

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

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

如中所述目标约束与非线性约束在同一函数中，使用时无法节省时间ga通过使用的嵌套函数技术模式搜索在步骤4。相反,叫ga同时设置适当的选项。

选择= optimoptions (“嘎”，“UseParallel”，对）；rng默认的%的再现性抽搐%确定解决方案的时间[xsolution，distance，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解决之道不如解决之道模式搜索答案:122米对126米。ga需要更多的时间：大约12秒，而不是2秒以下；模式搜索串行和嵌套运行所需时间更少，运行时间不到1秒ga连续运行的时间更长，一次测试大约需要30秒。