非线性约束下的代理优化

这个例子展示了如何在代理优化中包含非线性不等式约束。该示例解决了一个带有非线性约束的ODE。这个例子并行优化ode演示如何使用接受非线性约束的其他求解器解决相同的问题。

有关此示例的视频概述，请参见代理优化．

问题描述

问题是要改变大炮的位置和角度，使炮弹尽可能飞出墙外。炮口初速300米/秒。这面墙有20米高。如果大炮离墙太近，它发射的角度就太陡，炮弹就飞不了多远。如果大炮离墙太远，炮弹飞得不够远。

非线性空气阻力使弹丸减速。阻力与速度的平方成正比，比例常数为0.02。重力作用于弹丸，使其以恒定的9.81 m/s^2向下加速。因此，轨迹的运动方程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)．初速大小为300(初速)，因此，只取决于初角，这是一个标量。对于初始角th初速度为xp0 = 300 * (cos (th), sin (th))．因此，优化问题只依赖于两个标量，使其成为二维问题。用水平距离和初始角度作为决策变量。

制定ODE模型

ODE求解器要求您将模型表述为一阶系统。增加轨迹矢量 $（ x_{1} （ t ）， x_{2} （ t ））$ 它的时间导数 $（ x_{1} ” （ t ）， x_{2} ” （ t ））$ 形成4-D轨迹向量。根据增广向量，微分方程变成

$\frac{d}{dt} x （ t ）＝［ \begin{array}{c} x_{3.} （ t ） \\ x_{4} （ t ） \\ - 0 ． 02 ‖ （ x_{3.} （ t ）， x_{4} （ t ）） ‖ x_{3.} （ t ） \\ - 0 ． 02 ‖ （ x_{3.} （ t ）， x_{4} （ t ）） ‖ x_{4} （ t ） - 9 ． 81 \end{array} ］．$

的cannonshotFile实现了这个微分方程。

类型cannonshot

函数f = cannonshot (~ x) f = [x (3), (4); x (3); x (4)];%初始化，得到f(1)和f(2)的正确NRM =范数(x(3:4)) * .02;速度的%范数乘以常数f(3) = -x(3)*nrm;%水平加速度f(4) = -x(4)*nrm - 9.81;%的垂直加速度

想象这个ODE的解从距离墙30米开始，初始角度为π/ 3．的plotcannonsolution函数使用数值来解微分方程。

类型plotcannonsolution

将初始二维点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);%插值y值plot(xs,ys) hold on plot([0,0]，[0,20]，'k') %绘制墙壁xlabel('水平距离')ylabel('轨迹高度')ylim([0 100]) legend('轨迹'，'墙壁'，'位置'，'NW') dist = xs(end);标题(sprintf('距离%f'，dist))暂停

plotcannonsolution使用fzero求出抛射物落地的时间，即它的高度为0。弹丸在15秒前落地，所以plotcannonsolution使用15作为ODE解的时间量。

x0 =(-30;π/ 3);dist = plotcannonsolution (x0);

图中包含一个axes对象。标题为Distance 76.750599的axes对象包含两个类型为line的对象。这些物体代表了轨迹，墙壁。

准备优化

为了优化初始位置和角度，编写一个类似于前面的绘图例程的函数。计算从任意水平位置和初始角度开始的轨迹。

包括合理的约束条件。水平位置不能大于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) %改变初始二维点x到4-D x0 x0 = [x(1);0;300*cos(x(2));300*sin(x(2))];%求解弹道sol = ode45(@cannonshot，[0,15]，x0);%寻找时间t时，弹道高度= 0 zerofnd = fzero(@(r)deval(sol,r,2)，[1e-2,15]);%找出当时的水平位置dist = deval(sol,zerofnd,1);当x = 0时，弹丸穿过壁面的高度是多少?If deval(sol,15,1) > 0 wallfnd = fzero(@(r)deval(sol,r,1)，[0,15]);身高=德瓦尔(sol wallfnd 2);Else height = deval(sol,15,2);end f.e 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. maxfunctionevaluments '设置的函数求值限制。

xsolution =1×2-28.4012 - 0.6161

距离= -125.9987

exitflag = 0

输出=结构体字段:Elapsedtime: 53.0425 funccount: 200 constrviolation: 8.0630e-04 ineq: 8.0630e-04 rngstate: [1x1 struct] message: '代理选择停止，因为它超过了函数计算限制设置…'

画出最终轨迹。

Figure dist = plotcannonsolution(xsolution);

图中包含一个axes对象。标题为Distance 125.998707的axes对象包含两个类型为line的对象。这些物体代表了轨迹，墙壁。

的patternsearch解决方案并行优化ode表示最后的距离125.9880，和这个几乎一样surrogateopt解决方案。

另请参阅

surrogateopt