Модель масса-пружина-демпфер¶
В этом примере мы изучим два подхода к моделированию классической динамической модели, включающей один точечный объект с заданной массой, соединенный с опорой при помощи пружины и демпфера. Запустим модели и сравним результаты при выборе решателей с постоянным и с переменным шагом интегрирования.
using Plots
gr()
Описание моделей¶
Динамическая система, включающая массу, пружину и демпфер, очень часто приводится для изучения природы колебательных процессов, кинематики и эффектов хаотического движения (если объектов несколько) или методов численного интегрирования.
Мы рассмотрим направленную модель msd_causal.engee, которая осуществляет численное интегрирование уравнения $m \ddot x + k x + b \dot x = 0$.
Можно наглядно убедиться, что если выразить из этой формулы уравнение для ускорения груза, то именно оно и приходит на вход интегратора $\frac{1}{s^2}$.
$$\ddot x = \frac{-k x - b \dot x}{m}$$
Вторая модель msd_acausal.engee создана из ненаправленных (акаузальных) блоков библиотеки физического моделирования Engee. В верхней части находится датчик, который измеряет перемещение массы относительно точки опоры, а также блок SolverConfiguration, который указывает глобальному решателю модели, что подсоединенные к нему сигналы и блоки будет иметь собственный интегратор с собственными настройками.
В нижней части модели находится группа элементов, напоминающая стандартные учебные иллюстрации – масса $m$ соединена с опорой при помощи пружины с жесткостью $k$ и демпфера с коэффициентом поглощения $b$ (изменяемые параметры прописаны в свойствах этих блоков).
Благодаря графической природе модели эти элементы легко расположить самым разным образом, подключить параллельно или последовательно, или многократно дублировать.
Задание параметров моделей¶
Обе модели используют идентичные переменные параметры: масса груза m, жесткость пружины k и коэффициент поглощения демпфера b.
b,k,m = 1,2,3;
Когда нужные переменные созданы, можно будет управлять моделями, меняя значения переменных в Окне переменных. При значениях b=1, k=2, m=3 мы должны наблюдать довольно медленный, инерционный процесс с заметным затуханием.
Загрузка и изучение моделей¶
При помощи однострочной команды мы либо загрузим (load), либо откроем уже загруженную модель (open). Здесь мы пользуемся синтаксисом однострочных условных операторов x = условие ? если_истинно : если_ложно.
# Однострочная команда для открытия/загрузки модели
mc = "msd_causal" in [m.name for m in engee.get_all_models()] ? engee.open( "msd_causal" ) : engee.load( "$(@__DIR__)/msd_causal.engee" );
ma = "msd_acausal" in [m.name for m in engee.get_all_models()] ? engee.open( "msd_acausal" ) : engee.load( "$(@__DIR__)/msd_acausal.engee" );
Посмотрим, какие параметры интегратора прописаны в обеих моделях.
# Отсеим свойства модели, которые нас интересуют
mc_s = split( string( engee.get_param( "msd_causal" ) ), "#Параметры интегратора:\n")[2];
ma_s = split( string( engee.get_param( "msd_acausal" ) ), "#Параметры интегратора:\n")[2];
# Приведем их в соответствие стандартному форматированию CSV файлов
mc_s = replace( strip(mc_s, [' ', '\n', ')'] ), "=>"=>";" )
ma_s = replace( strip(ma_s, [' ', '\n', ')'] ), "=>"=>";" )
# Таблица для сравнения параметров моделей
using DelimitedFiles, DataFrames, PrettyTables
pretty_table(
coalesce.(outerjoin(
DataFrame( readdlm( IOBuffer(mc_s), ';' ), ["Свойство", "Направленная модель"] ),
DataFrame( readdlm( IOBuffer(ma_s), ';' ), ["Свойство", "Физическая модель"] ),
on="Свойство" ), "-"),
backend=:html, nosubheader=true, alignment=[:r, :c, :c]
)
Мы видим, что обе модели имеют одинаковые настройки интегратора на верхнем уровне. Интегратор физической модели задан в блоке Solver Configuration, на верхнем уровне опрос модели производится с постоянным шагом 0.01 с.
Выполнение моделей с переменным шагом¶
Выполним модели и сохраним полученные результаты в переменные rc и ra для направленной и физической моделей соответственно.
rc = engee.run( mc, verbose=false );
ra = engee.run( ma, verbose=false );
Сравнение результатов выполнения:
plot( rc["p"].time, rc["p"].value, title="Положение" )
p_pos = plot!( ra["p"].time, ra["p"].value )
plot( rc["v"].time, rc["v"].value, title="Скорость" )
p_vel = plot!( ra["v"].time, ra["v"].value )
plot( p_pos, p_vel, size=(600,200) )
В общем случае, на графиках будет разное количество точек, поэтому применим интерполяцию чтобы сравнить между собой результаты выполнения разных моделей. Таким действием мы вносим дополнительную ошибку, но, хочется надеяться, несущественную.
using Interpolations
rc_p_interp_linear = linear_interpolation(rc["p"].time, rc["p"].value)
rc_v_interp_linear = linear_interpolation(rc["v"].time, rc["v"].value)
using LaTeXStrings
d_pos = plot( ra["p"].value - rc_p_interp_linear.(ra["p"].time), title=L"\varepsilon_{положение}", legend=false)
d_vel = plot( ra["v"].value - rc_v_interp_linear.(ra["v"].time), title=L"\varepsilon_{скорость}", legend=false )
plot( d_pos, d_vel, size=(600,200) )
Разница $\varepsilon$ между результатами расчета, полученными при помощи физической модели с переменным шагом интегрирования, и направленной моделью с постоянным шагом достигает значений порядка $10^{-3}$.
Переключимся на решатель с постоянным шагом¶
Теперь мы получим решение направленной модели методом с переменным шагом интегрирования и посмотрим, как изменится ошибка относительно физической модели.
# Настроим направленную модель так, чтобы она решалась с переменным шагом интегрирования
engee.set_param!( "msd_causal", "SolverType"=>"variable-step" )
# Запустим модели
rc = engee.run( mc, verbose=false );
ra = engee.run( ma, verbose=false );
Чтобы сравнить графики с разным разрешением нам потребуется интерполяция.
rc_p_interp_linear = linear_interpolation(rc["p"].time, rc["p"].value)
rc_v_interp_linear = linear_interpolation(rc["v"].time, rc["v"].value)
e_pos = plot( ra["p"].value - rc_p_interp_linear.(ra["p"].time), title=L"\varepsilon_{положение}", legend=false )
e_vel = plot( ra["v"].value - rc_v_interp_linear.(ra["v"].time), title=L"\varepsilon_{скорость}", legend=false )
plot( e_pos, e_vel, size=(600,200) )
Расхождение моделей тоже имеет колебательный характер. Но, на этот раз, разница между физической моделью и направленной уже гораздо меньше. Расхождение достигает значений порядка $10^{-12}$.
В Engee встроен Инспектор данных, который позволит сравнить разные результаты запусков между собой.
Закрываем модели без сохранения¶
Если мы впоследствие захотим запустить этот пример еще раз с самого начала, то нужно закрыть модели без сохранения внесенных изменений.
#engee.close( "msd_causal", force=true );
#engee.close( "msd_acausal", force=true );
Заключение¶
В отличие от направленных моделей, физические модели, как правило, невозможно решить при помощи решателей с постоянным шагом интегрирования. Сложности в подстановке неявных начальных условий могут привести к численной нестабильности – например, к бесконечным градиентам.
Возможность исследовать модели разной природы и автоматизировано применять различные методы интегрирования позволяет нам выбрать оптимальное решение в обоих отношениях: наиболее понятный формализм моделирования и наиболее стабильный решатель.