这个例子展示了如何优化ODE的参数。
当ODE解返回目标函数和非线性约束函数时,如何避免两次计算。这个例子比较patternsearch
和遗传算法
运行求解器的时间和解决方案的质量。金宝搏官方网站
您需要一个并行计算工具箱™许可才能使用并行计算。
问题是改变火炮的位置和角度,让炮弹尽可能地飞出墙外。炮口初速为300米/秒。墙有20米高。如果加农炮离墙太近,就必须以太陡的角度发射,射弹就不能飞得足够远。如果加农炮离墙太远,炮弹也不能飞得足够远。
空气阻力使炮弹减速。阻力与速度的平方成正比,比例常数为0.02。重力作用于抛射体,使其以恒定的9.81 m/s向下加速2.因此,轨迹的运动方程x(t)
在哪里
初始位置x0
和初始速度xp0
是二维向量。不过,初始高度x0 (2)
是0,所以初始位置只取决于标量x0 (1)
.初始速度xp0
大小为300(枪口初速),所以只取决于初始角度,一个标量。对于初始角度th
,xp0
=300 * (cos (th), sin (th))
.因此,优化问题只依赖于两个标量,因此是一个二维问题。使用水平距离和角度作为决策变量。
ODE求解器要求您将模型表述为一阶系统。增加轨迹矢量(x1(t),x2(t)的时间导数(x'1(t),x'2(t)以形成一个4-D轨迹向量。用增广向量表示,微分方程就变成
将微分方程写成函数文件,并保存在MATLAB中®路径。
函数F = [x(3);x(4);x(4)];% Initial,得到f(1)和f(2)正确NRM = norm(x(3:4)) * .02;速度的范数乘以常数f (3) = - x(3) *全国抵抗运动;%的水平加速度F (4) = -x(4)*nrm - 9.81;%的垂直加速度
设想ODE的解从离墙30米的角度出发π/ 3
.
问题是找到初始位置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)德瓦尔(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 - 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...[],[],[],[], 磅,乌兰巴托,@cannonconstraint toc选择)
优化终止:网格尺寸小于选项。MeshTolerance和constraint violation小于options. constraintolerance。xsolution = -28.8123 0.6095 distance = -125.9880 eflag = 1 outpt = struct with fields: function: @cannonobjective problemtype: 'nonlinearconstr' pollmethod: ' gppositivebasis2n ' maxconstraint: 0 searchmethod: [] iterations: 5 funccount: 269 meshsize: 8.9125e-07 rngstate: [1×1 struct] message: '优化终止:meshsize小于选项。MeshTolerance↵和约束违背小于options. constraintolerance。运行时间为0.792152秒。
以0.6095弧度的角度从距离壁面约29米的地方开始发射,结果产生了最远的距离,约126米。报告的距离是负的,因为目标函数是到墙的距离的负。
可视化的解决方案。
x0 = [xsolution (1), 0; 300 * cos (xsolution(2)); 300 *罪(xsolution (2)));索尔=数值(@cannonfodder [0, 15], x0);找出投射物落地的时间zerofnd = fzero (@ (r)德瓦尔(sol r 2), [sol.x (2), sol.x(结束)]);t = linspace (0, zerofnd);%等于plot的次数x =德瓦尔(sol t 1);%插值x值y =德瓦尔(sol t 2);%插值y值情节(x, y)在情节([0,0],[0,20],“k”)%画墙包含(的水平距离) ylabel (的轨道高度)传说(“轨迹”,“墙”,“位置”,“西北”([0 70])保持从
目标约束函数和非线性约束函数都调用ODE求解器来计算它们的值。在目标约束与非线性约束在同一函数中避免调用求解器两次。的runcannon
文件实现了这种技术。复制此文件到您的MATLAB路径上的一个文件夹。
函数[x, f, eflag, outpt] = runcannon (x0,选择)如果输入参数个数= = 1%没有提供选项选择= [];结束xLast = [];最后一位ode求解器被调用索尔= [];% ODE溶液结构有趣= @objfun;%目标函数,嵌套在下面cfun = @constr;%约束函数,嵌套在下面磅= (-200;0.05);乌兰巴托=(1;π/ 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)));索尔=数值(@cannonfodder [0, 15], x0);xLast = x;结束%现在计算目标函数首先找出弹丸什么时候击中地面zerofnd = fzero (@ (r)德瓦尔(sol r 2), [sol.x (2), sol.x(结束)]);%然后计算当时的x位置y =德瓦尔(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)德瓦尔(溶胶,r, 1), [sol.x (1) sol.x(结束)]);然后求出这里的高度,用20减去C = 20 - deval(sol,zerofnd,2);结束结束
重新初始化问题并确定调用的时间runcannon
.
x0 =(-30;π/ 3);Tic [xsolution,distance,eflag,outpt] = runcannon(x0,opts);toc
优化终止:网格尺寸小于选项。MeshTolerance和constraint violation小于options. constraintolerance。运行时间为0.670715秒。
解算器比以前跑得快。如果您检查这个解决方案,您会看到输出是相同的。
尝试通过并行计算来节省更多时间。首先打开一个并行的工作人员池。
parpool
使用“本地”配置文件启动parpool…connection to the parallel pool (number of workers: 6). ans = ProcessPool with properties: Connected: true NumWorkers: 6 Cluster: local attachdfiles: {} AutoAddClientPath: true IdleTimeout: 30 minutes (30 minutes remaining) SpmdEnabled: true
设置选项以使用并行计算,并重新运行求解器。
选择= optimoptions (“patternsearch”选择,“UseParallel”,真正的);x0 =(-30;π/ 3);Tic [xsolution,distance,eflag,outpt] = runcannon(x0,opts);toc
优化终止:网格尺寸小于选项。MeshTolerance和constraint violation小于options. constraintolerance。运行时间为1.894515秒。
在这种情况下,并行计算速度较慢。如果您检查这个解决方案,您会看到输出是相同的。
你也可以尝试使用遗传算法来解决这个问题。然而,遗传算法通常速度较慢,可靠性较差。
在解释目标约束与非线性约束在同一函数中,使用时无法节省时间遗传算法
的嵌套函数技术patternsearch
在步骤4。相反,叫遗传算法
同时设置适当的选项。
选择= optimoptions (“遗传算法”,“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秒。
的遗传算法
解决之道不如解决之道patternsearch
答案:122米对126米。遗传算法
需要更多的时间:大约12秒而不是2秒;patternsearch
在串行和嵌套中花费的时间更少,不到1秒。运行遗传算法
连续运行的时间更长,一次测试大约需要30秒。