非线性约束的替代优化
此示例说明如何在替代优化中包含非线性不等式约束。该示例解决了具有非线性约束的 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 模型
ODE 求解器要求您将模型制定为一阶系统。用轨迹向量 的时间导数 来增强该轨迹向量,以形成四维轨迹向量。就这个增广向量而言,微分方程变成
cannonshot
文件实现了这个微分方程。
type cannonshot
function f = cannonshot(~,x) f = [x(3);x(4);x(3);x(4)]; % initial, gets f(1) and f(2) correct nrm = norm(x(3:4)) * .02; % norm of the velocity times constant f(3) = -x(3)*nrm; % horizontal acceleration f(4) = -x(4)*nrm - 9.81; % vertical acceleration
可视化该 ODE 的解,从距墙壁 30 米处开始,初始角度为 pi/3
。plotcannonsolution
函数使用 ode45
来求解微分方程。
type plotcannonsolution
function dist = plotcannonsolution(x) % Change initial 2-D point x to 4-D x0 x0 = [x(1);0;300*cos(x(2));300*sin(x(2))]; sol = ode45(@cannonshot,[0,15],x0); % Find the time when the projectile lands zerofnd = fzero(@(r)deval(sol,r,2),[sol.x(2),sol.x(end)]); t = linspace(0,zerofnd); % equal times for plot xs = deval(sol,t,1); % interpolated x values 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;pi/3]; dist = plotcannonsolution(x0);
准备优化
为了优化初始位置和角度,编写一个类似于前面绘图例程的函数。从任意水平位置和初始角度开始计算轨迹。
包括合理的边界约束。水平位置不能大于 0。设置上界为 -1。同样,水平位置不能低于 -200,因此设置下界为 -200。初始角度必须为正数,因此将其下界设置为 0.05。初始角度不应超过 pi
/2;将其上界设置为 pi
/2 – 0.05。
lb = [-200;0.05]; ub = [-1;pi/2-.05];
编写一个目标函数,在给定初始位置和角度的情况下,返回与墙壁距离的负数。如果轨迹在小于 20 的高度处穿过墙壁,则轨迹不可行;此约束是非线性约束。cannonobjcon
函数实现目标函数计算。为了实现非线性约束,函数调用 fzero
来查找射弹的 x 值为零的时间。该函数通过检查时间 15 之后射弹的 x 值是否大于零来解释 fzero
函数失败的可能性。如果不是,那么该函数将跳过查找射弹穿过墙壁的时间的步骤。
type cannonobjcon
function f = cannonobjcon(x) % Change initial 2-D point x to 4-D x0 x0 = [x(1);0;300*cos(x(2));300*sin(x(2))]; % Solve for trajectory sol = ode45(@cannonshot,[0,15],x0); % Find time t when trajectory height = 0 zerofnd = fzero(@(r)deval(sol,r,2),[1e-2,15]); % Find the horizontal position at that time dist = deval(sol,zerofnd,1); % What is the height when the projectile crosses the wall at x = 0? if deval(sol,15,1) > 0 wallfnd = fzero(@(r)deval(sol,r,1),[0,15]); height = deval(sol,wallfnd,2); else height = deval(sol,15,2); end f.Ineq = 20 - height; % height must be above 20 % Take negative of distance for maximization f.Fval = -dist; end
您已经计算出一条可行的初始轨迹。使用该值作为初始点。
fx0 = cannonobjcon(x0); fx0.X = x0;
使用 surrogateopt
解决优化问题
设置 surrogateopt
选项以使用初始点。为了实现可再现性,请将随机数生成器设置为 default
。使用 'surrogateoptplot'
绘图函数。运行优化。要了解 'surrogateoptplot'
图,请参阅 解释 surrogateoptplot。
opts = optimoptions('surrogateopt','InitialPoints',x0,'PlotFcn','surrogateoptplot'); rng default [xsolution,distance,exitflag,output] = surrogateopt(@cannonobjcon,lb,ub,opts)
surrogateopt stopped because it exceeded the function evaluation limit set by 'options.MaxFunctionEvaluations'.
xsolution = 1×2
-28.4081 0.6160
distance = -125.9988
exitflag = 0
output = struct with fields:
elapsedtime: 24.7447
funccount: 200
constrviolation: 8.4347e-04
ineq: 8.4347e-04
rngstate: [1x1 struct]
message: 'surrogateopt stopped because it exceeded the function evaluation limit set by ...'
绘制最终轨迹。
figure dist = plotcannonsolution(xsolution);
并行优化 ODE 中的 patternsearch
解显示最终距离为 125.9880
,与这个 surrogateopt
解几乎相同。