Документация

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

Проблема состоит в том, чтобы сменить положение и угол орудия, чтобы запустить снаряд в максимально возможной степени за стеной. Орудие имеет скорость морды 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];

Вы уже вычислили одну выполнимую начальную траекторию. Включайте то значение в структуру, которая дает и исходное положение и угол, и задает отрицание получившегося расстояния. Дайте отрицание расстояния, потому что вы хотите максимизировать расстояние, что означает, минимизируют отрицание расстояния.

x0 = struct('X',[-30,pi/3],'Fval',-dist);

Запишите целевую функцию, которая возвращает отрицание получившегося расстояния от стены, учитывая исходное положение и угол. Будьте осторожны при использовании fzero функция, потому что, если вы запускаете ОДУ на очень отрицательном расстоянии или очень мелком углу, траектория никогда не может пересекать стену. Кроме того, если траектория пересекает стену на высоте меньше чем 20, траектория неосуществима.

Чтобы составлять эту недопустимость, установите целевую функцию для неосуществимой траектории к высокому значению. Как правило, выполнимое решение будет иметь отрицательную величину. Таким образом, значение целевой функции 0 представляет плохое решение.

cannonobjconstraint функционируйте реализует вычисление целевой функции, учитывая нелинейное ограничение путем обнуления целевой функции при неосуществимых точках. Функция составляет возможность отказа в fzero функция при помощи try/catch операторы.

type cannonobjconstraint

function f = cannonobjconstraint(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
    try
        zerofnd = fzero(@(r)deval(sol,r,2),[sol.x(2),sol.x(end)]);
    catch
        zerofnd =  0;
    end
    % Find the horizontal position at that time
    f = deval(sol,zerofnd,1);
    % What is the height when the projectile crosses the wall at x = 0?
    try
        wallfnd = fzero(@(r)deval(sol,r,1),[sol.x(1),sol.x(end)]);
    catch
        wallfnd = 0;
    end
    height = deval(sol,wallfnd,2);
    if height < 20
        f = 0; % Objective for hitting wall
    end
    % Take negative of distance for maximization
    f = -f; 
end

Решите оптимизацию Используя surrogateopt

Установите surrogateopt опции, чтобы использовать начальную точку. Для воспроизводимости, набор генератор случайных чисел, чтобы принять значение по умолчанию. Используйте 'surrogateoptplot' функция plot. Запустите оптимизацию. Изучать 'surrogateoptplot' постройте, смотрите, Интерпретируют surrogateoptplot.

opts = optimoptions('surrogateopt','InitialPoints',x0,'PlotFcn','surrogateoptplot');
rng default
[xsolution,distance,exitflag,output] = surrogateopt(@cannonobjconstraint,lb,ub,opts)

Surrogateopt stopped because it exceeded the function evaluation limit set by 
'options.MaxFunctionEvaluations'.

xsolution = 1×2

  -26.6884    0.6455

distance = -125.8115

exitflag = 0

output = struct with fields:
       rngstate: [1×1 struct]
      funccount: 200
    elapsedtime: 29.0432
        message: 'Surrogateopt stopped because it exceeded the function evaluation limit set by ↵'options.MaxFunctionEvaluations'.'

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

figure
dist = plotcannonsolution(xsolution);

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