Matlab ode45 извлечение параметров

Я экспериментирую с ode45 в Matlab. Я научился передавать параметры в функцию ode, но у меня все еще есть вопрос. Предположим, что я хочу вычислить траекторию (профиль скорости) автомобиля, и у меня есть функция, например. getAcceleration, что дает мне ускорение автомобиля, а также правильную передачу: [acceleration, gear] = getAcceleration(speed,modelStructure), где modelStructure представляет модель автомобиля.

Функция оды будет:

function [dy] = car(t,y,modelStructure)

dy           = zeros(2,1);
dy(1)        = y(2);
[dy(2),gear] = getAcceleration(y(1),modelStructure);

Затем я вызываю интегратор Ode45 таким образом:

tInit = 0;
tEnd  = 5,
[t,y] = ode45(@car,[tInit tEnd], [speedInitial,accelerationInitial],options,modelStructure);

Проблема в том, как мне получить вектор, хранящий шестерни? Должен ли я иметь что-то вроде [t,y,gear]=ode45(....) или gear должен быть внутри вектора y?

Я работал над своим кодом и, используя функцию событий, теперь могу получать изменения «передач» автомобиля (как события). Теперь у меня есть новая проблема, связанная с тем же кодом. Представьте, что когда я оцениваю вектор 'dy', я могу получить дополнительное значение Z, которое позволяет мне значительно ускорить вызов вычисления ускорения (getAcceleration):

function [dy] = car(t,y,modelStructure)

dy           = zeros(2,1);
dy(1)        = y(2);
[dy(2),Z(t)] = getAcceleration(y(1),modelStructure,Z(t-1));

и предположим, что я также могу вычислить Z в начальных условиях. Проблема в том, что я не могу вычислить производную Z.

Есть ли способ передать значение Z, бросив степпинг без его интеграции?

Спасибо, парни.


parameters matlab ode
person Matlab User    schedule 26.11.2012    source источник

Ответы (1)


arrow_upward
5
arrow_downward

Прежде всего: почему начальными значениями дифференциального уравнения являются начальная скорость (speedInitial) и начальное ускорение (accelerationInitial)? Это означает, что дифференциальное уравнение car будет вычислять ускорение (y(1)) и рывок (y(2)), производную от ускорения по времени, в каждый момент времени t. Это кажется неверным... Я бы сказал, что начальными значениями должны быть начальное положение (positionInitial) и начальная скорость (speedInitial). Но я не знаю вашу модель, могу ошибаться.

Теперь, получая gear в решении напрямую: вы не можете, не взломав ode45. Это также логично; вы также не можете всегда получать dy напрямую, не так ли? Просто ode45 устроен не так.

Здесь я вижу два выхода:

Глобальная переменная

ВНИМАНИЕ: не используйте этот метод. Это только для того, чтобы показать, что большинство людей сделали бы с первой попытки.

Вы можете сохранить gear в глобальной переменной. Вероятно, это наименьшее количество кода, но и наименее удобный результат:

global ts gear ii

ii    = 1;
tInit = 0;
tEnd  = 5,
[t,y] = ode45(...
    @(t,y) car(t,y,modelStructure), ...
    [tInit tEnd], ...
    [speedInitial, accelerationInitial], options);

...

function [dy] = car(t,y,modelStructure)
global ts gear ii

dy    = zeros(2,1);
dy(1) = y(2);
[dy(2),gear(ii)] = getAcceleration(y(1),modelStructure);

ts(ii) = t;
ii = ii + 1;

Но из-за природы ode45 вы получите массив значений времени ts и связанного с ним gear, который содержит промежуточные точки и/или точки, которые были отклонены ode45. Итак, вам придется отфильтровать их позже:

ts( ~ismember(ts, t) ) = [];

Повторю еще раз: это НЕ метод, который я бы рекомендовал. Используйте глобальные переменные только при тестировании или выполнении каких-то быстрых и грязных вещей, но всегда очень быстро переключайтесь на другие решения. Кроме того, глобальные переменные увеличиваются на каждой (под)итерации ode45, что является неприемлемым снижением производительности.

Лучше использовать следующий метод:

Звонок после решения

Это также не слишком сложно для вашего случая, и я бы порекомендовал вам пойти. Сначала измените дифференциальное уравнение, как показано ниже, и решите как обычно:

tInit = 0;
tEnd  = 5,
[t,y] = ode45(...
    @(t,y) car(t,y,modelStructure), ...
    [tInit tEnd], ...
    [speedInitial, accelerationInitial], options);

...

function [dy, gear] = car(t,y,modelStructure)    

dy    = [0;0];
dy(1) = y(2);
[dy(2),gear] = getAcceleration(y(1),modelStructure);

а затем после завершения ode45 сделайте следующее:

gear = zeros(size(t));
for ii = 1:numel(t)
    [~, gear(ii)] = car(t(ii), y(ii,:).', modelStructure); 
end

Это даст вам все передачи, которые машина могла бы иметь в разы t.

Единственный недостаток, который я здесь вижу, заключается в том, что у вас будет намного больше вычислений функций car, чем ode45 будет использовать сам по себе. Но это реальная проблема только в том случае, если каждое вычисление car занимает порядка секунд или дольше, что, я подозреваю, не так в вашей настройке.

person Rody Oldenhuis    schedule 27.11.2012
comment
Привет, Роди! Спасибо за отличный ответ! Прежде всего, вы совершенно правы насчет начального условия: они должны быть speedInitial и positionInitial. Я очень удивлен из-за этого ограничения, а что, если бы у меня было событие, связанное с механизмом? Я имею в виду, что было бы возможно и осуществимо иметь своего рода «управление передачей» во время интеграции, как я мог бы это сделать? Я использую неправильный метод интеграции? Еще раз спасибо. - person Matlab User; 27.11.2012
comment
@MatlabUser: Ну, о ode45 можно многое узнать. Например, то, что вы предлагаете, возможно с использованием так называемых функций событий; см. документацию. Это, вероятно, займет некоторое время, чтобы освоить, но оно того стоит; это очень мощная концепция, с которой я еще не сталкивался ни в одном другом инструменте на любом другом языке. - person Rody Oldenhuis; 27.11.2012
comment
@MatlabUser: Кроме того, как я уже сказал, это не совсем ограничение. Просто то, что вы делаете, уже выходит далеко за рамки любого обычного заурядного ОДУ, для чего предназначен ode45. - person Rody Oldenhuis; 27.11.2012
comment
хорошо, Роди! Я буду следовать вашему совету и исследовать использование событий. Спасибо. - person Matlab User; 27.11.2012