В этом примере показано, как оптимизировать параметры ОДУ.
Это также показывает, как постараться не вычислять объективную и нелинейную ограничительную функцию дважды, когда решение для ОДУ возвращает обоих. Пример сравнивает patternsearch
и ga
в терминах времени, чтобы запустить решатель и качество решений.
Вам нужна лицензия Parallel Computing Toolbox™, чтобы использовать параллельные вычисления.
Проблема состоит в том, чтобы сменить положение и угол орудия, чтобы запустить снаряд в максимально возможной степени за стенкой. Орудие имеет скорость морды 300 м/с. Стенка 20 м высотой. Если орудие слишком близко к стенке, оно должно выстрелить в слишком крутой угол, и снаряд не перемещается достаточно далеко. Если орудие слишком далеко от стенки, снаряд не перемещается достаточно далеко также.
Сопротивление воздуха замедляет снаряд. Сила сопротивления пропорциональна квадрату скорости с коэффициентом пропорциональности 0.02. Сила тяжести действует на снаряд, ускоряя его вниз с постоянными 9.81 m/s2. Поэтому уравнения движения для траектории x (t)
где
Исходное положение x0
и начальная скорость xp0
2D векторы. Однако начальная высота x0(2)
0, таким образом, исходное положение зависит только от скалярного x0(1)
. И начальная скорость xp0
имеет величину 300 (скорость морды), поэтому зависит только от начального угла, скаляра. Для начального угла th
, xp0
= 300*(cos(th),sin(th))
. Поэтому задача оптимизации зависит только от двух скаляров, таким образом, это - 2D проблема. Используйте горизонтальное расстояние и угол как переменные решения.
Решатели ОДУ требуют, чтобы вы сформулировали свою модель как систему первого порядка. Увеличьте вектор траектории (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(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
И объективный и нелинейный ограничительный вызов функции решатель ОДУ, чтобы вычислить их значения. Используйте метод в Объективных и Нелинейных Ограничениях в Той же Функции (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.
Решатель запустился быстрее, чем прежде. Если вы исследуете решение, вы видите, что выход идентичен.
Попытайтесь сэкономить больше времени путем вычисления параллельно. Начните путем открытия параллельного пула рабочих.
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.
В этом случае параллельные вычисления были медленнее. Если вы исследуете решение, вы видите, что выход идентичен.
Можно также попытаться решить задачу с помощью генетического алгоритма. Однако генетический алгоритм обычно медленнее и менее надежен.
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 с.