В этом примере показано, как включить нелинейные ограничения неравенства в суррогатную оптимизацию. В примере решается ОДУ с нелинейным ограничением. Пример Оптимизация ОДУ параллельно показывает, как решить ту же задачу с помощью других решателей, которые принимают нелинейные ограничения.
Видеообзор этого примера см. в разделе Суррогатная оптимизация.
Проблема состоит в том, чтобы изменить положение и угол пушки для стрельбы снаряда как можно дальше за стеной. Пушка имеет дульную скорость 300 м/с. Высота стены - 20 м. Если пушка находится слишком близко к стене, она стреляет под слишком крутым углом, и снаряд недостаточно далеко движется. Если пушка находится слишком далеко от стены, снаряд едет недостаточно далеко.
Нелинейное сопротивление воздуха замедляет снаряд. Сила сопротивления пропорциональна квадрату скорости с постоянной пропорциональности 0,02. Сила тяжести воздействует на снаряд, ускоряя его вниз при постоянных 9,81 м/с ^ 2. Следовательно, уравнения движения для траектории x (t) равны
0,9,81),
где )/dt.
Исходное положение x0 и начальная скорость xp0 являются 2-D векторами. Однако начальная высота x0(2) равно 0, поэтому начальное положение задается скаляром x0(1). Начальная скорость имеет величину 300 (дульная скорость) и, следовательно, зависит только от начального угла, который является скалярным. Для начального угла th, начальная скорость равна xp0 = 300*(cos(th),sin(th)). Поэтому задача оптимизации зависит только от двух скаляров, что делает ее задачей 2-D. В качестве переменных принятия решения используйте расстояние по горизонтали и начальный угол.
Решатели ОДУ требуют, чтобы модель была сформулирована как система первого порядка. Увеличьте вектор траектории (t)) с его времени
x4 (t)) ‖ x4 (t) -9,81].
cannonshot файл реализует это дифференциальное уравнение.
type cannonshotfunction 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
Визуализировать решение данного ОДУ, начиная с 30 м от стены с начальным углом pi/3. plotcannonsolution функции использует ode45 для решения дифференциального уравнения.
type plotcannonsolutionfunction 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 в качестве времени для раствора ОДУ.
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 снаряда равно нулю. Функция учитывает возможность отказа в fzero , проверяя, больше ли по истечении времени 15 значение x снаряда, чем ноль. Если нет, то функция пропускает этап нахождения времени, когда снаряд проходит стенку.
type cannonobjconfunction 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' график см. в разделе Интерпретировать суррогатеоптплот.
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.3776 0.6165
distance = -125.9986
exitflag = 0
output = struct with fields:
elapsedtime: 38.6263
funccount: 200
constrviolation: 7.6395e-04
ineq: 7.6395e-04
rngstate: [1x1 struct]
message: 'surrogateopt stopped because it exceeded the function evaluation limit set by ...'
Постройте график конечной траектории.
figure dist = plotcannonsolution(xsolution);
patternsearch решение в разделе Оптимизация ОДУ параллельно показывает конечное расстояние 125.9880, что почти то же самое, что и это surrogateopt решение.