Описание задачи

Задача состоит в том, чтобы изменить положение и угол пушки, чтобы выпустить снаряд как можно дальше за стенкой. Пушка имеет дульную скорость 300 м/с. Стенка составляет 20 м. Если пушка слишком близко к стенке, она стреляет под слишком крутым углом, и снаряд не проходит достаточно далеко. Если пушка слишком далеко от стенки, снаряд не проходит достаточно далеко.

Нелинейное сопротивление воздуха замедляет снаряд. Сила сопротивления пропорциональна квадрату скорости с константой пропорциональности 0,02. Гравитация воздействует на снаряд, ускоряя его вниз с константой 9,81 м/с ^ 2. Поэтому уравнения движения для траектории x (t)

где $$v(t)=dx(t)/dt$$.

Начальное положение x0 и начальная скорость xp0 являются 2-D векторами. Однако начальная высота x0(2) равен 0, поэтому начальное положение задается скаляром x0(1). Начальная скорость имеет величину 300 (дульная скорость) и, следовательно, зависит только от начального угла, который является скаляром. Для начального угла th, начальная скорость xp0 = 300*(cos(th),sin(th)). Поэтому задача оптимизации зависит только от двух скаляров, что делает ее 2-D задачей. Используйте горизонтальное расстояние и начальный угол в качестве переменных принятия решений.

Формулировка модели ОДУ

Решатели ОДУ требуют от вас сформулировать модель как систему первого порядка. Увеличение вектора траектории $$(x_1(t),x_2(t))$$ с его производной по времени $$(x_1\prime (t),x_2\prime (t))$$ для формирования 4-D вектора траектории. В терминах этого дополненного вектора дифференциальное уравнение становится

The 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

Визуализируйте решение этой ОДУ начиная с 30 м от стенки с начальным углом pi/3. The 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 в качестве количества времени для решения ОДУ.

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, траектория недопустима; это ограничение является нелинейным ограничением. The cannonobjcon функция реализует вычисление целевой функции. Чтобы реализовать нелинейное ограничение, функция вызывает fzero чтобы найти время, когда значение X снаряда равняется нулю. Функция учитывает возможность отказа в fzero функция путем проверки, является ли после времени 15 значение X снаряда больше нуля. Если нет, то функция пропускает шаг нахождения времени, когда снаряд проходит стенку.

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.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);

The patternsearch решение в Optimize a ODE in Parallel показывает окончательное расстояние 125.9880, что почти так же, как и это surrogateopt решение.