Estimating the parameters of a system of differential equations from inaccurate observations

remark

This work is based on material licensed under the Creative Commons Attribution 4.0 (CC BY 4.0) International License. Source material: “ESTIMATION OF THE PARAMETERS OF DIFFERENTIAL EQUATIONS FROM INEXACT OBSERVATIONS” by G. Sh. Tsitsiashvili and M. A. Osipov, available at https://journals.vsu.ru/sait/article/view/10247/9468 (date of access: 07/14/2023). No changes were made to the original algorithm/method, but only its implementation was presented and an experiment was carried out.

I note that this work (code) uses a string calculator that I previously created, the operation of which was described in a previous article (https://habr.com/ru/articles/747920/). In this regard, we will not dwell on its functions in this article.

Motivation and about the article

The article mentioned above describes a method that allows one to estimate the parameters of a system of differential equations. Since this article aroused my interest, I decided to conduct a numerical experiment in the Scala language, which is closest to me, in order to see real results and get to know the algorithm in more detail.

This is what the Lorentz system looks like, consisting of three differential equations. Usually it is described in such notation, but I will rewrite it so that the notation matches those given in the article and displayed in the code.

Here yi represent variables, and bi are the constants we are going to evaluate from the observations.

The article considers such systems, and it is easy to verify that it corresponds to the Lorentz system (m = 3).

It is stated that if there are some inaccurate observations yi’ with a time step hi,

That

From the system, *10, estimates of our parameters bi are derived, which are approximately equal to bi0

The first step is data preparation.

The first thing I will do is define the variables we are setting and define the data formats.

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 – this is a necessary crutch, since I have not yet figured out how to solve arbitrary systems of equations (perhaps my future article will be devoted to this issue). Matrix solutions are suitable only for systems of linear equations.

For further work, it is also necessary to set up a string calculator so that it can work with all the operations that are in systems And reverseSystemsnamely: “+”, “-“, “*”, “/”.

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)))

}  
  

For this we use classes StringCalculator, ConstantFabric AndBinaryOperationFabric (see my article about the string calculator) and add the appropriate operations.

To add bi into equations, we represent them as constants in a string calculator.

Step 2 – generating inaccurate observations

first, let’s make a small method that will calculate the values ​​​​of functions with our variables 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)})")))
  }

This method replaces all yi in formulas for the corresponding values ​​from values and then calculates the value of the function at the point Mi(yk,1, ..., yk,m) (method calculate – part of the string calculator).

Since the Lorentz *1,2 system does not have an analytical solution, we will solve the system numerically. This is beneficial in a sense, since we need Fuzzy Observations.

As a numerical method, we will use the most elementary method – the Euler method. To understand exactly how I will use it, one formula is enough.

Let us separately generate observations for positive 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))
  })

The algorithm of this code is as follows:

1. A variable is created positivePreciseObservationswhich is initialized with a list containing the initial observation yi0 as the first element. This list will contain a sequence of precise observations of the system.

2. Then the method is applied foldLeft to the range from 1 to n. In every iteration itr represents the current iteration number.

3. Within each iteration, the following happens:

   Using the method zipWithIndex applied function fun_index to each element of the list systems (list of system functions) along with the corresponding index.

   For every pair (fun, index) a new value is calculated using the last value from acc (last exact observation) for the corresponding index index. Here fun represents a function, and acc.last(fun_index._2) gets the last value from acc for index index.

   The result of the calculation is added to acc using the operator :+. This updates the list accadding a new accurate observation of the system.

4. At the end of all iterations, a list is obtained positivePreciseObservationscontaining a sequence of “exact” system values ​​*6 after each iteration in time steps hkWhere k takes values ​​from -n before n.

To generate observations for negative values kperform the following steps:

For the reverse step, the formula is derived from Equation *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

Then concatenate the values ​​for negative and positive values k and perform the transposition of the list of observations.

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

Step 3 – generating inaccurate observations.

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)
}

This step adds an error of up to 20% of the value to each value in the “exact values” list.
Math.pow(-1, rnd.nextInt()) defines the direction of deflection (positive or negative)
rnd.nextDouble() * 0.2 determines the amount of deviation.

Step 4 – Implementing the Parameter Estimation Algorithm

Implementation of the parameter estimation algorithm:

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 = {
    ...
  }

}

We will create a trait that will be responsible for this functionality. To work with the string calculator, do not forget to add the “extendsStringCalculator” trait.

estimator – the method itself, which will return parameter estimates. He takes inaccurate observations nonPreciseObservationssystem solution *10 reverseSystems And h – step.

replaceVariableValues is a method that replaces si to the corresponding values. In our case, this will be used to replace yi And 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)})"))
  }

To calculate the value n *9 from the available data, given that the number of steps is always equal to 2 * n + 1the formula will take the following form.

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

Next, we will calculate the averages using *11.

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

We will also calculate the approximate values ​​of the functions for yi0 using *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
})

Everything here is exactly according to the formula, the only thing worth noting is that the range (-n to n) is used to go through all the steps.

The last thing left is to solve the *10 system and get approximate values bi.

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

We take the equation (as a string), replace the values Fi to the relevant Fi0and also replace yi on their average values ​​at the zero point (k = 0).

Now the class looks like this:

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)})"))
  }

}

Last step

Let’s mix a trait SystemParameterEstimator to the class where the experiment is being carried out using the keyword with SystemParameterEstimator.

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

Then we call the method estimate on the received data and display the result.

  println(estimator(nonPreciseObservations, reverseSystems, h))

Result

with parameters:

  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)

i got the result:

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

Well, if the results of the estimation do not differ much from the real values ​​of the parameters, this indicates the good performance of the algorithm and the accuracy of the estimation. This may indicate that your method is effectively estimating the parameters of a system of differential equations based on inaccurate observations.

Conclusion

I know that one example is not enough and I could do experiments on different systems with different parameters, plot frequency histograms, plot the dependence of the average deviation on the variance, and so on.

This is exactly what I’m going to do, but in the next article)

Since I have never done such work, I will ask the commentators to write on which systems it is worth doing, and what exactly is worth demonstrating in the article (frequency histograms, graphs, and so on) for a larger and more accurate study and creating a more complete picture.

PS The project can be found on GitHub

https://github.com/PicoPicoRobotWoman/estimation_of_parameters_of_differentia_equations