Описание проблемы

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

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

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

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

Сформулируйте модель ОДУ

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

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. 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, траектория неосуществима; это ограничение является нелинейным ограничением. 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' функция plot. Запустите оптимизацию. Изучать '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.4071    0.6160

distance = -125.9990

exitflag = 0

output = struct with fields:
        elapsedtime: 28.9019
          funccount: 200
    constrviolation: 9.8552e-04
               ineq: 9.8552e-04
           rngstate: [1×1 struct]
            message: 'Surrogateopt stopped because it exceeded the function evaluation limit set by ↵'options.MaxFunctionEvaluations'.'

Постройте итоговую траекторию.

figure
dist = plotcannonsolution(xsolution);

patternsearch решение в Оптимизирует ОДУ, параллельно показывает итоговое расстояние 125.9880, который является почти тем же самым как этим surrogateopt решение.