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

Трейлер

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

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

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

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

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

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

В статье рассматриваются такие системы, и нетрудно убедиться, что она соответствует системе Лоренца (m = 3).

Мы говорим, что если имеются неточные наблюдения yi’ с временным шагом hi,

Этот

Из системы *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 (см. мой пост в String Calculator) и добавьте соответствующие операции.

Добавлять 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 не имеет аналитического решения, будем решать систему численно. В некотором смысле это выгодно, поскольку нам нужны нечеткие наблюдения.

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

Отдельно сгенерируем наблюдения для 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))
  })

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

1. Создается переменная positivePreciseObservationsкоторый инициализируется списком, содержащим начальное наблюдение yi0 как первый пункт. Этот список будет содержать последовательность конкретных наблюдений за системой.

2. Затем применяется метод foldLeft в диапазоне от 1 до n. На каждой итерации itr представляет текущий номер итерации.

3. В каждой итерации происходит следующее:

   Использование метода zipWithIndex прикладная функция fun_index к каждому пункту в списке systems (список системных функций) с соответствующим индексом.

   Для каждой пары (fun, index) новое значение вычисляется с использованием последнего значения acc (последнее точное наблюдение) для соответствующего индекса index. Здесь fun представляет собой функцию, и acc.last(fun_index._2) получить последнее значение acc для индекса index.

   Результат расчета добавляется к acc с помощью оператора :+. Это обновляет список accдобавление нового точного наблюдения за системой.

4. В конце всех итераций получается список 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 до 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 соответственно

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

Заключение

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

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

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

PS Проект можно найти на GitHub