PicoPicoRobotWoman 15 июл 2023 в 02:10

Оценка параметров системы дифференциальных уравнений по неточным наблюдениям

Сложный

7 мин

2.4K

Программирование*Алгоритмы*Scala*Математика*Функциональное программирование*

Ремарка

Эта работа основана на материале, распространяемом под лицензией Creative Commons Attribution 4.0 (CC BY 4.0) International License. Исходный материал: "ОЦЕНКА ПАРАМЕТРОВ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ПО НЕТОЧНЫМ НАБЛЮДЕНИЯМ" от Г. Ш. Цициашвили и М. А. Осипова, доступно по адресу https://journals.vsu.ru/sait/article/view/10247/9468 (дата обращения: 14.07.2023). В оригинальный алгоритм/метод не быловнесено изменений, а лишь представлена его реализация и проведен эксперимент.

Отмечу, что в данной работе (коде) используется ранее созданный мною строковый калькулятор, функционирование которого было описано в предыдущей статье (https://habr.com/ru/articles/747920/). В связи с этим, не будем останавливаться на его функциях в данной статье.

Мотивация и о статье

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

$\left\{ \begin{aligned} \frac{dx}{dt} &= \sigma (y - x) \\ \frac{dy}{dt} &= x (\rho - z) - y \\ \frac{dz}{dt} &= xy - \beta z \\ \end{aligned} \right. (*1)$

Вот так выглядит система Лоренца, состоящая из трех дифференциальных уравнений. Обычно она описывается именно в таких обозначениях, но я перепишу ее так, чтобы обозначения совпадали с теми, которые приведены в статье и отображены в коде.

$\left\{ \begin{aligned} \frac{dy_1}{dt} &= b_1 (y_2 - y_1) \\ \frac{dy_2}{dt} &= y_1 (b_2 - y_3) - y_2 \\ \frac{dy_3}{dt} &= y_1y_2 - b_3 y_3 \\ \end{aligned} \right. (*2)$

Здесь yi представляют переменные, а bi - константы, которые мы собираемся оценить по наблюдениям.

$\frac{dy_i}{dt} = F_i (y_1 ,\ldots , y_m , b_1 ,\ldots , b_m ), \quad i = \overline{1, m} \ \ (*3)$

Статья рассматривает подобные системы, и несложно убедиться, что системе Лоренца она соответствует (m = 3).

Утверждается, что если имеются некоторые неточные наблюдения yi' с временным шагом hi,

${y_i}^{'}(ht) = y_i(ht) + \epsilon(hi) \ \ \ (*4) \\ где \\ {y_i}^{'}(ht) - неточное\ наблюдение\ в\ момент\ ht\ (*5)\\ {y_i}^{}(ht) - точное \ значение\ значение\ переменой\ в\ момент\ ht\ (*6) \\ \epsilon(ht) - отклонение\ от\ точного\ значения\ (*7) \\ h - \ шаг\ h = n^{-\alpha} \ где\ \alpha \in (1, 1.5) (*8) \\ t = \overline{-n, n}\ , n - количество\ наблюдений (*9)$

то

$F_i (\overline{y_1^0}, \ldots, \overline{y_m^0}, \beta_1^0, \ldots, \beta_m^0) = \overline{F_i^0} (*10)\\где \\ \overline{y_i^0} = \frac{\sum_{i=-n}^{n} y_i}{2 * n + 1} - проще\ говоря\ среднее (*11)\\ \overline{F_i^0} = \frac{\sum_{t=-n}^{n} (y_i(ht)*(ht))}{\sum_{\substack{t=-n \\ }}^{n} (ht)^2 } (*12) \\ собственно\ \beta_i^0 - это\ и\ есть\ наша\ оценка\ параметра (*13)$

Из системы, *10, выводятся оценки наших параметров bi, которые примерно равны bi0

Шаг первый - подготовка данных.

Первое, что я сделаю, это определюсь с переменными, которые мы задаем, и определюсь с форматами данных.

object LorenzAttractor extends App {
  
  val n = 1000 // *9
  val alpha = 1.25 // *8
  val h = Math.pow(n, -alpha) // *8

  val bi: List[Double] = List(2, 3, 1) // фиксированные параметры bi
  val yi0: List[Double] = List(1, 2, 1) // начальне значения во временом наге равном 0
  val systems: List[String] = List("b1*(y2-y1)", "y1*(b2-y3)-y2", "y1*y2-b3*y3") // сисиема уравнений *2
  val reverseSystems: List[String] = List("F1/(y2-y1)", "(F2+y2)/y1+y3", "(y1*y2 - F3)/y3") // решение системы *10

}

reverseSystem - это необходимый костыль, так как я еще не придумала, как решать произвольные системы уравнений (возможно, этому вопросу будет посвящена моя будущая статья). Матричные решения подходят только для систем линейных уравнений.

Для дальнейшей работы также необходимо настроить строковый калькулятор, чтобы он мог работать со всеми операциями, которые есть в systems и reverseSystems, а именно: "+", "-", "*", "/".

object LorenzAttractor extends App with StringCalculator with ConstantFabric with BinaryOperationFabric {
    
  addBinaryOperation(createBinary("+", 1, (left, right) => left + right, 0))
  addBinaryOperation(createBinary("-", 1, (left, right) => left - right, 0))
  addBinaryOperation(createBinary("*", 2, (left, right) => left * right, 1))
  addBinaryOperation(createBinary("/", 2, (left, right) => left / right, 1))
  
  ...
  
  bi.zipWithIndex.foreach(b_index => addConstant(createConstant(s"b${b_index._2 + 1}", b_index._1)))

}

Для этого мы используем классы StringCalculator, ConstantFabric иBinaryOperationFabric (смотрите мою статью про строковый калькулятор) и добавляем соответствующие операции.

Чтобы добавить bi в уравнения, мы представляем их как константы в строковом калькуляторе.

Шаг 2 - генерация неточных наблюдений

для начала сделаем небольшой метод, который будет вычеслять значения функций с нашими переменными yk,i

  private def calkFun(fun: String, values: List[Double]): Double = {
    calculate(values.indices.foldLeft(fun)((acc, i) => acc.replace(s"y${i + 1}", s"(${values(i)})")))
  }

Данный метод осуществляет замену всех yi в формулах на соответствующие значения из values и затем вычисляет значение функции в точке Mi(yk,1, ..., yk,m) (метод calculate - часть строкового калькулятора).

Так как система Лоренца *1,2 не имеет аналитического решения, мы будем решать систему численными методами. Это в некотором смысле выгодно, поскольку нам требуются Неточные наблюдения.

В качестве численного метода мы воспользуемся самым элементарным методом - методом Эйлера. Чтобы понять, как именно я им воспользуюсь, достаточно одной формулы.

$y_{(k+1),i} = y_{k,i} + h f_i(y_{k,1}, ..., y_{k,m}) (*14) \\$

Отдельно сгенерируем наблюдения для положительных k.

  val positivePreciseObservations: List[List[Double]] = (1 to n).foldLeft(List(yi0))((acc, itr) => {
    acc :+ systems.zipWithIndex.map(fun_index => acc.last(fun_index._2) + h * calkFun(fun_index._1, acc.last))
  })

Алгоритм работы данного кода следующий:

Создается переменная positivePreciseObservations, которая инициализируется списком, содержащим начальное наблюдение yi0 в качестве первого элемента. Этот список будет содержать последовательность точных наблюдений системы.
Затем применяется метод foldLeft к диапазону от 1 до n. В каждой итерации itr представляет текущий номер итерации.
Внутри каждой итерации происходит следующее:
С помощью метода zipWithIndex применяется функция fun_index к каждому элементу списка systems (список функций системы) вместе с соответствующим индексом.
Для каждой пары (fun, index) выполняется вычисление нового значения, используя последнее значение из acc (последнее точное наблюдение) для соответствующего индекса index. Здесь fun представляет функцию, а acc.last(fun_index._2) получает последнее значение из acc для индекса index.
Результат вычисления добавляется в acc с помощью оператора :+. Это обновляет список acc, добавляя новое точное наблюдение системы.
По завершении всех итераций получается список positivePreciseObservations, содержащий последовательность "точных" значений системы *6 после каждой итерации во временных шагах hk, где k принимает значения от -n до n.

Для генерации наблюдений для отрицательных значений k, выполним следующие шаги:

Для обратного шага формула выводится из уравнения *14.

  val negativePreciseObservations: List[List[Double]] = (1 to n).foldRight(List(yi0))((itr, acc) => {
    systems.zipWithIndex.map(fun_index => acc.head(fun_index._2) - h * calkFun(fun_index._1, acc.head)) :: acc
  }).init

Затем объединим значения для отрицательных и положительных значений k и выполним транспонирование списка наблюдений.

val preciseObservations: List[List[Double]] = (negativePreciseObservations ::: positivePreciseObservations).transpose

Шаг 3 - генерация неточных наблюдений.

val rnd = new Random()

val nonPreciseObservations: List[List[Double]] = preciseObservations.map { nums =>
  nums.map(num => (1 + Math.pow(-1, rnd.nextInt()) * rnd.nextDouble() * 0.2) * num)
}

В этом шаге к каждому значению из списка "точных значений" добавляется ошибка до 20% от значения.
Math.pow(-1, rnd.nextInt()) определяет направление отклонения (положительное или отрицательное)
rnd.nextDouble() * 0.2 определяет величину отклонения.

Шаг 4 - реализация алгоритма оценки параметров

Реализация алгоритма оценки параметров:

trait SystemParameterEstimator extends StringCalculator {

  def estimator(nonPreciseObservations: List[List[Double]], reverseSystems: List[String], h: Double): List[Double] = {

    ...

  }

  private def replaceVariableValues(fun: String,symbol: String , values: List[Double]): String = {
    ...
  }

}

Мы создадим трейт, который будет отвечать за эту функциональность. Для работы со строковым калькулятором не забываем подмешать трейт "extendsStringCalculator".

estimator - собственно метод, который будет возвращать оценки параметров. Он принимает неточные наблюдения nonPreciseObservations, решение системы *10 reverseSystems и h - шаг.

replaceVariableValues - метод, который заменяет si на соответствующие значения. В нашем случае, это будет использоваться для замены yi и Fi.

  private def replaceVariableValues(fun: String,symbol: String , values: List[Double]): String = {
    values.indices.foldLeft(fun)((acc, i) => acc.replace(s"$symbol${i + 1}", s"(${values(i)})"))
  }

Для вычисления значения n *9 из доступных данных, учитывая, что количество шагов всегда равно 2 * n + 1, формула примет следующий вид.

val n = (nonPreciseObservations.head.size - 1 ) / 2

Далее мы рассчитаем средние значения, используя *11.

val yi0 = nonPreciseObservations.map(yi => yi.sum / yi.size)

Также мы вычислим приближенные значения функций при yi0 с использованием *12.

val Fi0 = nonPreciseObservations.map(yi => {
  (-n to n).zip(yi).map(k_y => k_y._1 * k_y._2 * h).sum / (-n to n).map(k => Math.pow(k * h, 2)).sum
})

Тут все точно по формуле, единственное, что стоит отметить, что для прохождения всех шагов используется диапазон (-n to n).

Последнее, что остается, - это решить систему *10 и получить приближенные значения bi.

reverseSystems.zipWithIndex.map( sys_index => {
  val expression1: String = replaceVariableValues(sys_index._1, "F", Fi0)
  val expression2: String = replaceVariableValues(expression1, "y", yi0)
  calculate(expression2)
})

Берем уравнение (в виде строки), заменяем значения Fi на соответствующие Fi0, а также подменяем yi на их средние значения в нулевой точке (k = 0).

Теперь класс выглядит следующим образом:

trait SystemParameterEstimator extends StringCalculator {

  def estimator(nonPreciseObservations: List[List[Double]], reverseSystems: List[String], h: Double): List[Double] = {

    val n = (nonPreciseObservations.head.size - 1 ) / 2

    val yi0 = nonPreciseObservations.map(yi => yi.sum / yi.size)

    val Fi0 = nonPreciseObservations.map(yi => {
      (-n to n).zip(yi).map(k_y => k_y._1 * k_y._2 * h).sum / (-n to n).map(k => Math.pow(k * h, 2)).sum
    })

    reverseSystems.zipWithIndex.map( sys_index => {
      val expression1: String = replaceVariableValues(sys_index._1, "F", Fi0)
      val expression2: String = replaceVariableValues(expression1, "y", yi0)
      calculate(expression2)
    })

  }

  private def replaceVariableValues(fun: String,symbol: String , values: List[Double]): String = {
    values.indices.foldLeft(fun)((acc, i) => acc.replace(s"$symbol${i + 1}", s"(${values(i)})"))
  }

}

Последний шаг

Давайте подмешаем трейт SystemParameterEstimator в класс, где проводится эксперимент, с помощью ключевого слова with SystemParameterEstimator.

object LorenzAttractor extends App with StringCalculator with ConstantFabric with BinaryOperationFabric with SystemParameterEstimator {

Затем вызовем метод estimate на полученных данных и выведем результат.

  println(estimator(nonPreciseObservations, reverseSystems, h))

Результат

при параметрах:

  val n = 1000
  val alpha = 1.25
  val h = Math.pow(n, -alpha)

  val bi: List[Double] = List(2, 3, 1)
  val yi0: List[Double] = List(1, 2, 1)

я получила результат:

List(1.9609705152167594, 3.1061318517029735, 0.9548191490247522) - b1, b2, b3 соответственно

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

Заключение

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

Этим собственно я и собираюсь заняться, но уже в следующей статье)

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

P.S. с проектом можно ознакомится на GitHub

https://github.com/PicoPicoRobotWoman/estimation_of_parameters_of_differentia_equations

Теги:

Хабы: