非线性约束代理优化

打开实时脚本

这个例子展示了如何在代理优化中包含非线性不等式约束。该示例求解带有非线性约束的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 ））$ 形成一个四维轨迹向量。关于增广向量，微分方程就变成了

$\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)];% initial，得到f(1)和f(2)正确的NRM = norm(x(3:4)) * .02;速度的%范数乘以常数f(3) = -x(3)*nrm;%水平加速度f(4) = -x(4)*nrm - 9.81;垂直加速度%

可视化的解决这个ODE开始30米从墙的初始角度π/ 3．的plotcannonsolution函数使用数值来解微分方程。

类型plotcannonsolution

将初始二维点x改为四维点x0 x0 = [x(1);0;300*cos(x(2));300*sin(x(2))];Sol = ode45(@cannonshot，[0,15]，x0);找到抛射物降落的时间zerofnd = fzero(@(r)deval(sol,r,2)，[sol.x(2)，sol.x(end)]);T = linspace(0,zerofnd);= deval(sol,t,1);%插值x值ys = deval(sol,t,2);%插值y值plot(xs,ys)按住plot([0,0]，[0,20]，'k') %绘制墙壁xlabel('水平距离')ylabel('弹道高度')ylim([0 100]) legend('弹道'，'墙壁'，'位置'，'NW') dist = xs(end);title(sprintf('距离%f'，dist))暂停

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

X0 = [-30;pi/3];Dist = plotcannonsolution(x0);

图中包含一个轴对象。标题为Distance 76.750599的axis对象包含2个类型为line的对象。这些物体代表轨迹、墙。

准备优化

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

包括合理的约束。水平位置不能大于0。设置上限为-1。同样，水平位置不能低于-200，因此设置下限为-200。初始角度必须为正，因此将其下限设置为0.05。初始角度不应超过π/ 2;设它的上界为π/2 - 0.05。

Lb = [-200;0.05];Ub = [-1;pi/2-.05];

写一个目标函数，在给定初始位置和角度的情况下，返回与墙的最终距离的负数。如果弹道在小于20的高度穿过墙壁，则弹道不可行的;这个约束是一个非线性约束。的cannonobjcon函数实现了目标函数的计算。为了实现非线性约束，函数调用fzero求出抛射物的x值为零的时间。该函数说明了失败的可能性fzero通过检查时间15后，弹丸的x值是否大于零来函数。如果不是，则函数跳过寻找抛射物通过墙壁的时间这一步。

类型cannonobjcon

function f = cannonobjcon(x) %将初始二维点x改为四维点x0 x0 = [x(1);0;300*cos(x(2));300*sin(x(2))];%求解轨迹sol = ode45(@cannonshot，[0,15]，x0);%找到时间t时，弹道高度= 0 zerofnd = f0 (@(r)deval(sol,r,2)，[1e-2,15]);找到当时的水平位置dist = deval(sol,zerofnd,1);x = 0时，弹丸穿过墙壁时的高度是多少?如果deval(sol,15,1) > 0 wallfnd = f0 (@(r)deval(sol,r,1)，[0,15]);高度= deval(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．

Opts = optimoptions(“surrogateopt”，“InitialPoints”x0,“PlotFcn”，“surrogateoptplot”）;rng默认的[xsolution,distance,exitflag,output] = surrogateopt(@cannonobjcon,lb,ub,opts)

图优化图函数包含一个轴对象。标题为Best: -125.999现任:-125.913现任:-123.14的坐标轴对象包含9个line类型的对象。这些对象表示最佳、现任、初始样本、随机样本(Infeas)、随机样本、自适应样本(Infeas)、自适应样本、现任(Infeas)、代理重置。

surrogateopt停止，因为它超过了'options.MaxFunctionEvaluations'设置的函数计算限制。

xsolution =1×2-28.4012 - 0.6161

Distance = -125.9987

Exitflag = 0

输出=带字段的结构:Elapsedtime: 58.5070 funccount: 200 constrviolation: 8.0630e-04 ineq: 8.0630e-04 rngstate: [1x1 struct] message: 'surrogateopt停止，因为它超过了函数计算限制…'

画出最后的轨迹。

图dist = plotcannonsolution(xsolution);

图中包含一个轴对象。标题为Distance 125.998707的axis对象包含2个类型为line的对象。这些物体代表轨迹、墙。

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

另请参阅

surrogateopt