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

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

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

Нелинейное сопротивление воздуха замедляет снаряд. Сила сопротивления пропорциональна квадрату скорости с коэффициентом пропорциональности 0.02. Сила тяжести действует на снаряд, ускоряя его вниз с постоянными 9.81 m/s^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, траектория неосуществима.

Чтобы составлять этот infeasibility, установите целевую функцию для неосуществимой траектории к высокому значению. Как правило, выполнимое решение будет иметь отрицательную величину. Таким образом, значение целевой функции 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'. Запустите оптимизацию. Чтобы понять сюжет '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.