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

Шаг 1. Определение проблемы.

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

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

где $$v(t)\text{}=dx(t)/dt.$$

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

Шаг 2. Сформулируйте модель ODE.

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

Запишите дифференциальное уравнение как файл функции и сохраните его на вашем MATLAB^® path.

function f = cannonfodder(t,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.

 Код для генерации фигуры

x0 = [-30;0;300*cos(pi/3);300*sin(pi/3)];
sol = ode45(@cannonfodder,[0,10],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')
legend('Trajectory','Wall','Location','NW')
ylim([0 120])
hold off

Шаг 3. Решите использование patternsearch.

Проблема состоит в том, чтобы найти исходное положение x0(1) и начальный угол x0(2) максимизировать расстояние от стены земли снаряда. Поскольку это - проблема максимизации, минимизируйте отрицание расстояния (см. Максимизацию по сравнению с Минимизацией).

Использовать patternsearch чтобы решить эту задачу, необходимо обеспечить цель, ограничение, исходное предположение и опции.

Эти два файла являются ограничительными функциями и целью. Скопируйте их в папку на вашем пути MATLAB.

function f = cannonobjective(x)
x0 = [x(1);0;300*cos(x(2));300*sin(x(2))];

sol = ode45(@cannonfodder,[0,15],x0);

% Find the time t when y_2(t) = 0
zerofnd = fzero(@(r)deval(sol,r,2),[sol.x(2),sol.x(end)]);
% Then find the x-position at that time
f = deval(sol,zerofnd,1);

f = -f; % take negative of distance for maximization



function [c,ceq] = cannonconstraint(x)

ceq = [];
x0 = [x(1);0;300*cos(x(2));300*sin(x(2))];

sol = ode45(@cannonfodder,[0,15],x0);

if sol.y(1,end) <= 0 % Projectile never reaches wall
    c = 20 - sol.y(2,end);
else
    % Find when the projectile crosses x = 0
    zerofnd = fzero(@(r)deval(sol,r,1),[sol.x(2),sol.x(end)]);
    % Then find the height there, and subtract from 20
    c = 20 - deval(sol,zerofnd,2);
end

Заметьте, что цель и ограничительные функции устанавливают их входную переменную x0 к 4-D начальной точке для решателя ОДУ. Решатель ОДУ не останавливается, если снаряд врезается в стену. Вместо этого ограничительная функция просто становится положительной, указывая на неосуществимое начальное значение.

Исходное положение x0(1) не может быть выше 0, и бесполезно иметь его быть ниже –200. (Это должно быть близко –20, потому что без сопротивления воздуха самая долгая траектория запустилась бы в –20 под углом pi/4.) Точно так же начальный угол x0(2) не может быть ниже 0 и не может быть выше pi/2. Установите границы немного далеко от этих начальных значений:

lb = [-200;0.05];
ub = [-1;pi/2-.05];
x0 = [-30,pi/3]; % Initial guess

Установите UseCompletePoll опция к true. Это дает решение более высокого качества и включает прямое сравнение с параллельной обработкой, потому что вычисление параллельно требует этой установки.

opts = optimoptions('patternsearch','UseCompletePoll',true);

Вызвать patternsearch решать задачу.

tic % Time the solution
[xsolution,distance,eflag,outpt] = patternsearch(@cannonobjective,x0,...
    [],[],[],[],lb,ub,@cannonconstraint,opts)
toc

Optimization terminated: mesh size less than options.MeshTolerance
 and constraint violation is less than options.ConstraintTolerance.

xsolution =

  -28.8123    0.6095


distance =

 -125.9880


eflag =

     1


outpt = 

  struct with fields:

         function: @cannonobjective
      problemtype: 'nonlinearconstr'
       pollmethod: 'gpspositivebasis2n'
    maxconstraint: 0
     searchmethod: []
       iterations: 5
        funccount: 269
         meshsize: 8.9125e-07
         rngstate: [1×1 struct]
          message: 'Optimization terminated: mesh size less than options.MeshTolerance↵ and constraint violation is less than options.ConstraintTolerance.'

Elapsed time is 0.792152 seconds.

Запуск снаряда приблизительно 29 м от стены под углом 0,6095 результата радиана в самом дальнем расстоянии, приблизительно 126 м. Расстояние, о котором сообщают, отрицательно, потому что целевая функция является отрицанием расстояния до стены.

Визуализируйте решение.

x0 = [xsolution(1);0;300*cos(xsolution(2));300*sin(xsolution(2))];

sol = ode45(@cannonfodder,[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')
legend('Trajectory','Wall','Location','NW')
ylim([0 70])
hold off

Шаг 4. Постарайтесь не вызывать дорогую стандартную подпрограмму дважды.

И объективный и нелинейный ограничительный вызов функции решатель ОДУ, чтобы вычислить их значения. Используйте метод в Объективных и Нелинейных Ограничениях в Той же Функции, чтобы не вызывать решатель дважды. runcannon файл реализует этот метод. Скопируйте этот файл в папку на вашем пути MATLAB.

function [x,f,eflag,outpt] = runcannon(x0,opts)

if nargin == 1 % No options supplied
    opts = [];
end

xLast = []; % Last place ode solver was called
sol = []; % ODE solution structure

fun = @objfun; % The objective function, nested below
cfun = @constr; % The constraint function, nested below

lb = [-200;0.05];
ub = [-1;pi/2-.05];

% Call patternsearch
[x,f,eflag,outpt] = patternsearch(fun,x0,[],[],[],[],lb,ub,cfun,opts);

    function y = objfun(x)
        if ~isequal(x,xLast) % Check if computation is necessary
            x0 = [x(1);0;300*cos(x(2));300*sin(x(2))];
            sol = ode45(@cannonfodder,[0,15],x0);
            xLast = x;
        end
        % Now compute objective function
        % First find when the projectile hits the ground
        zerofnd = fzero(@(r)deval(sol,r,2),[sol.x(2),sol.x(end)]);
        % Then compute the x-position at that time
        y = deval(sol,zerofnd,1);
        y = -y; % take negative of distance
    end

    function [c,ceq] = constr(x)
       ceq = [];
        if ~isequal(x,xLast) % Check if computation is necessary
            x0 = [x(1);0;300*cos(x(2));300*sin(x(2))];
            sol = ode45(@cannonfodder,[0,15],x0);
            xLast = x;
        end
        % Now compute constraint functions
        % First find when the projectile crosses x = 0
        zerofnd = fzero(@(r)deval(sol,r,1),[sol.x(1),sol.x(end)]);
        % Then find the height there, and subtract from 20
        c = 20 - deval(sol,zerofnd,2);
    end

end

Повторно инициализируйте проблему и время вызов runcannon.

x0 = [-30;pi/3];
tic
[xsolution,distance,eflag,outpt] = runcannon(x0,opts);
toc

Optimization terminated: mesh size less than options.MeshTolerance
 and constraint violation is less than options.ConstraintTolerance.
Elapsed time is 0.670715 seconds.

Решатель запустился быстрее, чем прежде. Если вы исследуете решение, вы видите, что выход идентичен.

Шаг 5. Вычислите параллельно.

Попытайтесь сэкономить больше времени путем вычисления параллельно. Начните путем открытия параллельного пула рабочих.

parpool

Starting parpool using the 'local' profile ...
Connected to the parallel pool (number of workers: 6).

ans = 

 ProcessPool with properties: 

            Connected: true
           NumWorkers: 6
              Cluster: local
        AttachedFiles: {}
    AutoAddClientPath: true
          IdleTimeout: 30 minutes (30 minutes remaining)
          SpmdEnabled: true

Установите опции использовать параллельные вычисления и повторно запускать решатель.

opts = optimoptions('patternsearch',opts,'UseParallel',true);
x0 = [-30;pi/3];
tic
[xsolution,distance,eflag,outpt] = runcannon(x0,opts);
toc

Optimization terminated: mesh size less than options.MeshTolerance
 and constraint violation is less than options.ConstraintTolerance.
Elapsed time is 1.894515 seconds.

В этом случае параллельные вычисления были медленнее. Если вы исследуете решение, вы видите, что выход идентичен.

Шаг 6. Сравните с генетическим алгоритмом.

Можно также попытаться решить задачу с помощью генетического алгоритма. Однако генетический алгоритм обычно медленнее и менее надежен.

Как объяснено в Объективных и Нелинейных Ограничениях в Той же Функции, вы не можете сэкономить время при использовании ga методом вложенной функции, используемым patternsearch на Шаге 4. Вместо этого вызовите ga параллельно путем установки подходящих вариантов.

options = optimoptions('ga','UseParallel',true);
rng default % For reproducibility
tic % Time the solution
[xsolution,distance,eflag,outpt] = ga(@cannonobjective,2,...
    [],[],[],[],lb,ub,@cannonconstraint,options)
toc

Optimization terminated: average change in the fitness value less than options.FunctionTolerance
 and constraint violation is less than options.ConstraintTolerance.

xsolution =

  -37.6217    0.4926


distance =

 -122.2181


eflag =

     1


outpt = 

  struct with fields:

      problemtype: 'nonlinearconstr'
         rngstate: [1×1 struct]
      generations: 4
        funccount: 9874
          message: 'Optimization terminated: average change in the fitness value less than options.FunctionTolerance↵ and constraint violation is less than options.ConstraintTolerance.'
    maxconstraint: 0
       hybridflag: []

Elapsed time is 12.529131 seconds.

ga решение не так хорошо как patternsearch решение: 122 м по сравнению с 126 м. ga занимает больше времени: приблизительно 12 с по сравнению с под 2 с; patternsearch занимает еще меньше времени в сериале и вложенный, меньше чем 1 с. Выполнение ga последовательно берет еще дольше, приблизительно 30 с в одном тестовом прогоне.