Оптимизируйте ОДУ параллельно

В этом примере показано, как оптимизировать параметры ОДУ.

Это также показывает, как постараться не вычислять объективную и нелинейную ограничительную функцию дважды, когда решение для ОДУ возвращает обоих. Пример сравнивает patternsearch и ga в терминах времени, чтобы запустить решатель и качество решений.

Вам нужна лицензия Parallel Computing Toolbox™, чтобы использовать параллельные вычисления.

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

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

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

d2x(t)dt2=−0.02‖v(t)‖v(t)−(0,9.81),

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

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

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

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

ddtx(t)=[x3(t)x4(t)−.02‖(x3(t),x4(t))‖x3(t)−.02‖(x3(t),x4(t))‖x4(t)−9.81].

Запишите дифференциальное уравнение как файл функции и сохраните его на вашем пути 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.

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

Шаг 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 =
         function: @cannonobjective
      problemtype: 'nonlinearconstr'
       pollmethod: 'gpspositivebasis2n'
    maxconstraint: 0
     searchmethod: []
       iterations: 5
        funccount: 269
         meshsize: 8.9125e-07
         rngstate: [1x1 struct]
          message: 'Optimization terminated: mesh size less than options.MeshTolerance…'

Elapsed time is 3.174088 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. Постарайтесь не вызывать дорогую стандартную подпрограмму дважды.

И объективный и нелинейный ограничительный вызов функции решатель ОДУ, чтобы вычислить их значения. Используйте метод в Объективных и Нелинейных Ограничениях в Той же Функции (Optimization Toolbox), чтобы не вызывать решатель дважды. 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
Elapsed time is 2.610590 seconds.

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

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

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

parpool
Starting parpool using the 'local' profile ... connected to 4 workers.

ans = 

 Pool with properties: 

            Connected: true
           NumWorkers: 4
              Cluster: local
        AttachedFiles: {}
          IdleTimeout: 30 minute(s) (30 minutes remaining)
          SpmdEnabled: true

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

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

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

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

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

runcannonga файл вызывает генетический алгоритм и избегает двойной оценки решателя ОДУ. Это напоминает runcannon, но вызовы ga вместо patternsearch, и также проверки, достигает ли траектория стены. Скопируйте этот файл в папку на вашем пути MATLAB.

function [x,f,eflag,outpt] = runcannonga(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 ga
[x,f,eflag,outpt] = ga(fun,2,[],[],[],[],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
        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
    end

end

Вызовите runcannonga параллельно.

opts = optimoptions('ga','UseParallel',true);
rng default % for reproducibility
tic
[xsolution,distance,eflag,outpt] = runcannonga(opts)
toc
Optimization terminated: average change in the fitness value less than 
options.FunctionTolerance and constraint violation is less than options.ConstraintTolerance.

xsolution =
  -17.9172    0.8417

distance =
 -116.6263

eflag =
     1

outpt = 
      problemtype: 'nonlinearconstr'
         rngstate: [1x1 struct]
      generations: 5
        funccount: 20212
          message: [1x140 char]
    maxconstraint: 0

Elapsed time is 119.630284 seconds.

ga решение не так хорошо как patternsearch решение: 117 м по сравнению с 126 м. ga занял намного больше времени: приблизительно 120 с по сравнению с под 5 с.

Связанные примеры

Больше о

Документация Global Optimization Toolbox
Поддержка
Для просмотра документации необходимо авторизоваться на сайте