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

Шаг 1. Определите задачу.

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

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

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

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

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

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

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

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. Решить с помощью pattersearch.

Задача состоит в том, чтобы найти начальное положение 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. Не вызывайте дорогую стандартную подпрограмму дважды.

Как целевая, так и нелинейная функция ограничения вызывают решатель ОДУ, чтобы вычислить их значения. Используйте метод в Objective и нелинейных ограничениях в той же Функции, чтобы избежать вызова решателя дважды. The 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 с за один тестовый запуск.