Модель масса-пружина-демпфер¶
В этом примере мы изучим два подхода к моделированию классической динамической модели, включающей один точечный объект с заданной массой, соединенный с опорой при помощи пружины и демпфера. Запустим модели и сравним результаты при выборе решателей с постоянным и с переменным шагом интегрирования.
Pkg.add(["Interpolations"])
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" );
Посмотрим, какие параметры интегратора прописаны в обеих моделях.
using DelimitedFiles, DataFrames
# Собрем параметры моделей
mc_dict = engee.get_param( "msd_causal" );
ma_dict = engee.get_param( "msd_acausal" );
df_causal = DataFrame("Свойство"=>collect(keys(mc_dict)), "Направленная модель"=>collect(values(mc_dict)))
df_acausal = DataFrame("Свойство"=>collect(keys(ma_dict)), "Физическая модель"=>collect(values(ma_dict)))
# Таблица для сравнения параметров моделей
df = outerjoin(df_causal, df_acausal, on="Свойство")
# Отсеим свойства модели, которые нас интересуют
params_to_display = ["StartTime", "StopTime", "SolverType", "SolverName", "FixedStep"]
df[in(params_to_display).(df."Свойство"), :]
Мы видим, что обе модели имеют одинаковые настройки интегратора.
На самом деле, алгоритм численного решения физической модели задан в блоке Solver Configuration, но верхнеуровневый интегратор модели, который принимает сигналы от физической подсистемы и возвращает измерения в виде векторов, выполняет их дискретизацию с постоянным шагом 0.01 с.
Выполнение моделей с переменным шагом¶
Выполним модели и сохраним полученные результаты в переменные rc и ra для направленной и физической моделей соответственно.
rc = engee.run( "msd_causal" );
ra = engee.run( "msd_acausal" );
Сравнение результатов выполнения:
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 LaTeXStrings
d_pos = plot( ra["p"].value - rc["p"].value, title=L"\varepsilon_{положение}", legend=false)
d_vel = plot( ra["v"].value - rc["v"].value, title=L"\varepsilon_{скорость}", legend=false )
plot( d_pos, d_vel, size=(600,200) )
Разница между результатами, очевидно, очень небольшая, порядок значения $10^{-3}$.
Посмотрим, что будет, если переключить "алгоритмическую" модель в тот же режим, в котором работает решатель физической модели.
Переключимся на решатель с постоянным шагом¶
Переключим характер работы обеих моделей. Теперь все решатели подбирают оптимальный размер шага и выдают данные с разной частотой дискретизации.
# Настроим направленную модель так, чтобы она решалась с переменным шагом интегрирования
engee.set_param!( "msd_causal", "SolverType"=>"variable-step" )
engee.set_param!( "msd_acausal", "SolverType"=>"variable-step" )
# Запустим модели
rc = engee.run( "msd_causal" );
ra = engee.run( "msd_acausal" );
Чтобы сравнить графики с разным разрешением нам потребуется интерполяция.
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)
# Вычисляем разницу между интерполяцией и расчетом физической системы
e_pos = plot( ra["p"].time, ra["p"].value - rc_p_interp_linear.(ra["p"].time), title=L"\varepsilon_{положение}", legend=false )
e_vel = plot( ra["v"].time, ra["v"].value - rc_v_interp_linear.(ra["v"].time), title=L"\varepsilon_{скорость}", legend=false )
plot( e_pos, e_vel, size=(600,200) )
На этот раз, разница между физической моделью и направленной уже гораздо меньше - расхождение имеет порядок значений $10^{-8}$.
В Engee встроен Инспектор данных, который позволит сравнить разные результаты запусков между собой.
Закрываем модели без сохранения¶
Если мы впоследствие захотим запустить этот пример еще раз с самого начала, то нужно закрыть модели без сохранения внесенных изменений.
engee.close( "msd_causal", force=true );
engee.close( "msd_acausal", force=true );
Заключение¶
В отличие от направленных моделей, физические модели, как правило, невозможно решить при помощи решателей с постоянным шагом интегрирования. Сложности в подстановке неявных начальных условий могут привести к численной нестабильности – например, к бесконечным градиентам.
Возможность исследовать модели разной природы и автоматизировано применять различные методы интегрирования позволяет нам выбрать оптимальное решение в обоих отношениях: наиболее понятный формализм моделирования и наиболее стабильный решатель.