Оценка параметров ДУ в Python

Одной из распространенных задач, возникающих при исследовании различных объектов — построение математической модели. Нередко математическая модель представляется в виде системы дифференциальных уравнений, однако непосредственное измерение всех, входящих в них параметров, как правило, невозможно по различным причинам. В таком случае, одним из подходов является проведение идентификационных экспериментов и оценка параметров ДУ путем решения оптимизационной задачи.

В статье рассмотрен простой способ оценки параметров системы ДУ в форме Коши на языке Python.

Постановка задачи

Разработка простого модуля Python для решения задачи оценки параметров.
Исходные данные для задачи:

  • Представление модели в форме Коши, первые N величин — наблюдаемы (порядок уравнений не имеет значения с точки зрения математики, но упорядочивание сильно упростит жизнь при разработке);
  • Имеются данные эксперимента, на основе которого будет производится оценка параметров/

Исходные данные

Исходными данными, передаваемыми при создании экземпляра класса оценки параметров являются:

  1. Функция, реализующая вычисление ДУ.
  2. Массив экспериментальных данных (Два массива: время и значения)

Таким образом, получаем следующий класс с конструктором:

Тут просто сохраняются значения параметров и функция в поля класса.

Оценка параметров

Процедуру оценки параметров можно сделать на основе функции scipy.optimize.leastsq. Эта функция минимизирует сумму квадратов полученных величин и может принимать два интересных для начального понимания аргумента: func и x0.

Функция func должна принимать один аргумент (который м.б. вектором), являющийся параметрами, а x0 — начальные значения параметров.

Метод оценки параметров можно реализовать следующим образом:

Метод получает начальные значения начальных условий ДУ и начальные значения параметров системы. Функцию вычисления вектора ошибок между решением системы ДУ с текущими значениями параметров и экспериментальными значениями необходимо реализовывать самостоятельно, потому, следующим шагом будет ее написание:

Функция получает вектор параметров системы (включая начальные условия ДУ).

Предложенная функция содержит вызов еще одной функции my_ls_func, выделенной для разделения функционала. Выглядит она следующим образом и вычисляет непосредственно решение системы ДУ при заданных параметрах:

Таким образом, процедура оценки параметров работает следующим образом:

  1. Вычисляется решение системы ДУ при текущих (или начальных) значения параметров.
  2. Рассчитывается вектор ошибки между решением и экспериментальными данными.
  3. Производится корректировка значений параметров.
  4. Переход на шаг 1, если не достигнуто максимальное число итераций или минимальное рассогласование.

Визуализация

Для визуализации напишем следующую функцию:

где calcODE — еще одна служебная функция для расчета системы ДУ:

Тестирование и пример использования

В качестве примера использования попробуем оценить параметры системы ДУ, описывающие маятник с трением.

\[ {d \theta \over dt } = \omega \\ {d \omega \over dt } = -b * \omega — c* \sin(\theta) \]

где \(\theta \) — угол, \(\omega\) — скорость.

Этот объект уже использовался в качестве иллюстрации к статье о фазовых портретах и частотному методу синтеза регуляторов.

Для получения «экспериментальных» данных воспользуемся решением системы ДУ с известными «истинными» параметрами:

Объявленная тут система будет использоваться и для оценки параметров, но уже без известных значений:

Получившиеся результаты являются достаточно точно оценкой:

Real parameter: b = 0.3, c = 5.0
Real initial condition: 1 0
Estimated parameter: b = 0.299993833612628, c = 5.000066486504925
Estimated initial condition: [ 9.99976783e-01 1.58491793e-04]

Ссылка на исходный код (GitLab)

Поделиться
  • 10
    Поделились

Добавить комментарий