这个例子展示了如何在代理优化中包含非线性不等式约束。这个例子求解了一个带有非线性约束的ODE。这个例子并行优化ODE演示如何使用接受非线性约束的其他求解器来解决相同的问题。
有关此示例的视频概述,请参见代理优化.
问题是改变火炮的位置和角度,让炮弹尽可能地飞出墙外。炮口初速为300米/秒。墙有20米高。如果加农炮离墙太近,就会以太陡的角度发射,射弹也不会飞得足够远。如果加农炮离墙太远,炮弹就飞得不够远。
非线性空气阻力使弹丸减速。阻力与速度的平方成正比,比例常数为0.02。重力作用于抛射体,使其向下加速,速度常数为9.81 m/s^2。因此,轨迹的运动方程x(t)
在哪里 .
初始位置x0
和初始速度xp0
是二维向量。不过,初始高度x0 (2)
是0,所以初始位置是由标量给出的x0 (1)
.初始速度的大小为300(初速),因此,只取决于初始角度,这是一个标量。对于初始角度th
,初速度是xp0 = 300 * (cos (th), sin (th))
.因此,优化问题只依赖于两个标量,使其成为一个二维问题。使用水平距离和初始角度作为决策变量。
ODE求解器要求您将模型表述为一阶系统。增加轨迹向量 它的时间导数 形成一个4-D轨迹向量。用增广向量表示,微分方程就变成
的cannonshot
文件实现了这个微分方程。
类型cannonshot
函数f = cannonshot (~ x) f = [x (3), (4); x (3); x (4)];NRM = norm(x(3:4)) * .02;速度范数乘以常数f(3) = -x(3)*nrm;%水平加速度f(4) = -x(4)*nrm - 9.81;%的垂直加速度
设想这个ODE的解从离墙30米开始,初始角度为π/ 3
.的plotcannonsolution
函数使用数值
解微分方程。
类型plotcannonsolution
函数dist = plotcannonsolution(x) %改变初始2- d点x为4-D x0 x0 = [x(1);0;300*cos(x(2));300*sin(x(2))];索尔=数值(@cannonshot [0, 15], x0);%求弹丸落地时间zerofnd = fzero(@(r)deval(sol,r,2),[sol.x(2),sol.x(end)]);t = linspace (0, zerofnd);对于plot xs = deval(sol,t,1);%插值x值ys = deval(sol,t,2);% interpolated y values plot(xs,ys) hold on plot([0,0],[0,20],'k') % Draw the wall xlabel('Horizontal distance') ylabel('Trajectory height') ylim([0 100]) legend('Trajectory',' wall ','Location','NW') dist = xs(end);title(sprintf('Distance %f',dist)) hold off
plotcannonsolution
使用fzero
求弹丸落地的时间,也就是它的高度为0。炮弹在15秒前落地,所以plotcannonsolution
使用15作为ODE解决方案的时间。
x0 =(-30;π/ 3);dist = plotcannonsolution (x0);
为了优化初始位置和角度,编写一个类似于前面绘图例程的函数。计算从任意水平位置和初始角度开始的轨迹。
包括合理的约束。水平位置不能大于0。设上界为-1。同样,水平位置不能低于-200,所以设置下限-200。初始角度必须是正的,所以将其下限设为0.05。初始角度不应超过π
/ 2;设它的上界为π
/ 2 - 0.05。
磅= (-200;0.05);乌兰巴托=(1;π/ 2 . 05);
编写一个目标函数,在给定初始位置和角度的情况下,返回到墙的结果距离的负。如果弹道以低于20的高度穿过墙,则弹道不可行;这个约束是一个非线性约束。的cannonobjcon
函数实现了目标函数的计算。为了实现非线性约束,函数调用fzero
求抛射体的x值为零的时刻。该函数解释了在fzero
函数通过检查,在时间15之后,弹丸的x值是否大于零。如果没有,则该函数跳过计算弹丸通过墙壁的时间这一步。
类型cannonobjcon
函数f = cannonobjcon(x) %改变初始2- d点x为4-D x0 x0 = [x(1);0;300*cos(x(2));300*sin(x(2))];%求解轨迹sol = ode45(@cannonshot,[0,15],x0);%求轨迹高度= 0时的时间t zerofnd = fzero(@(r)deval(sol,r,2),[1e-2,15]);dist = deval(sol,zerofnd,1);当x = 0时,炮弹穿过墙壁的高度是多少?If deval(sol,15) > 0 wallfnd = f0 (@(r)deval(sol, 15),[0,15]); / /身高=德瓦尔(sol wallfnd 2);Else height = deval(sol,15,2);f.Ineq = 20 - height;f.Fval = -dist; end
你已经算出了一个可行的初始轨迹。使用该值作为初始点。
fx0 = cannonobjcon (x0);fx0。X = x0;
surrogateopt
集surrogateopt
选择使用初始点。为了再现性,将随机数生成器设置为默认的
.使用“surrogateoptplot”
图的功能。运行优化。了解“surrogateoptplot”
情节,看到解释surrogateoptplot.
选择= optimoptions (“surrogateopt”,“InitialPoints”x0,“PlotFcn”,“surrogateoptplot”);rng默认的[xsolution、距离、exitflag、输出]= surrogateopt (@cannonobjcon、磅、乌兰巴托、选择)
surrogateopt停止,因为它超过了'options.MaxFunctionEvaluations'设置的函数求值限制。
xsolution =1×2-28.3776 - 0.6165
距离= -125.9986
exitflag = 0
输出=结构体字段:Elapsedtime: 50.7602 funccount: 200 construct: 7.6395e-04 ineq: 7.6395e-04 rngstate: [1x1 struct] message: 'surrogateopt停止了,因为它超过了由…
画出最后的轨迹。
图dist = plotcannonsolution(xsolution);
的patternsearch
解决方案并行优化ODE表示最后的距离125.9880
,几乎和这个一样surrogateopt
解决方案。