Fitting experimental data with differential equations : Optimization

[madoko]

Hello !!!

I'm a physics student trying to solve an experimental problem in fluid dynamics and here is the issue I'm having.
The data output of my experiment is a 2D trajectory ([X,Y] array). I also have a theoretical model in the form of 3 coupled differential equations, solved using Runge Kutta 4, which also gives me a 2D trajectory ([x,y] array). This model depends mainly on 3 constants (a,G,B) of unknown values.
Here is an equivalent python script of my problem:

import matplotlib.pyplot as plt
import numpy as np

XY = np.loadtxt("experimental_data.txt")
plt.plot(XY[:,0],XY[:,0],label="experimental data")

x,y = Solver_diff_eqt(a,B,G) #Solution of my model
plt.plot(x,y,label="model")

plt.legend()
plt.show()

Obviously my code is much more complicated, but this is the basic idea.
So finally my question. I'm looking for a tool or an algorithm which could optimize the values of (a,G,B) by trying to minimize the difference between my XY graph and my xy graph.

(I guess that a better way to say it is : "I've 2 distributions, one of which depends on 3 constants. I'm looking for a way to calculate these 3 constants so that the 2 distributions are the same.", but it's harder to understand)

To give you some context on what I already tried.
I tried doing it by hand with sliders but it's way too complicated.
I tried to use scipy.optimize.curve_fit with an intelligent hack, but I'm obviously not intelligent enough .

Now you know everything ... Can you please help me ?
Thanks

[micseydel]

I think your question is clear, but I'm not familiar with existing tools which do it. I think sympy is a good bet, though I'm not 100% sure you'll find what you're looking for there.

[Gribouillis]

You could define a functional that measures the distance between the two curves such as

X, Y = ...

def distance(a, B, G):
    x, y = Solver(a, B, G)
    return np.linalg.norm(y - Y)

Then minimize this function of 3 variables with scipy.optimize.minimize(distance, **args).

[madoko]

Thank you both for answering

You could define a functional that measures the distance between the two curves such as

X, Y = ...

def distance(a, B, G):
    x, y = Solver(a, B, G)
    return np.linalg.norm(y - Y)

Then minimize this function of 3 variables with scipy.optimize.minimize(distance, **args).

This looks promising Gribouillis ! It's going to take me some time implement because optimize.minimize has a lot of options and subtleties. I will post the solution to this problem if I'm successful.

[scidam]

Hm..., this is quite interesting problem. Traditional approach to solving such problems consist in
building iterative process converging to exact (or almost exact) solution (if this is possible).
The approach proposed by Gribouillis is correct in general, but being applied in practice, it likely lead
to very expensive and long computations. This is because Scipy's minimize function will
likely try to compute gradient of the residual function norm(y-Y) that, in turn, depends on solution of
the system of ODE (ordinary differential equations). Therefore, getting the gradient estimation will require a lot of computations.

Another approach assumes the following steps:

1) Problem statement. Let we have (three ODE's as stated above) a system of ODEs and observations:

dx/dt = F(x, y, p, a, B, G)
dy/dt = G(x, y, p, a, B, G)
dp/dt = H(x, y, p, a, B, G)
O(t) = (x(t), y(t)).T # O - observation vector, .T means transposed,i.e. O is column vector

Since a, B, G are constants the system above could be expanded as follows:

dx/dt = F(x, y, p, a, B, G)
dy/dt = G(x, y, p, a, B, G)
dp/dt = H(x, y, p, a, B, G)
da/dt = 0
dB/dt = 0
dG/dt = 0
O(t) = (x(t), y(t)).T

This is a classical dynamical system, which input parameters, i.e. initial conditions of the ODE, are needed to be found. These parameters are: x0, y0, p0, a0, B0, G0. The last triple a0, B0, G0 doesn’t change in time (they are constants).
So, we can rewrite the problem in vector-matrix notation as follows:

dX/dt = Z(X)
O(t) = M(X) + e(t)

where X = (x, y, p, a, B, G).T; Z -- nonlinear function that returns a vector of size 6 and depends on X;
O is a vector of available observations. It depends on time, and might be known only at specific time moments, e.g. t_k, k=1,...,N. In your case M is quite trivial, it returns just x and y. i.e. M = (X[0], X[1]).T = (x, y).T; e(t) -- error of measurements vector (we don't need to know it, we will use it
to state the optimization problem in least square sense).

In general, Z(X) could be of non-linear and/or non-stationary (i.e. depends on time, Z = Z(X,t)) form. I hope that your ODEs is stationary one.

So, if we have measurements (O(t_k)) on the trajectory, we need to find initial conditions X0=X(t0)
that leads to the trajectory that, in turn, satisfies (in least square sense)
to all of our observations. This means that the trajectory
minimizes the sum(e(t_k)^2).

sum e(t_k)^2 = sum (O(t_k) - M(X(X0, t_k))^2 => min by X0

This is non-linear (in general) optimization problem, but it could be solved
in iterative way as a series of linear ones.

Suppose that we have a good approximation of the solution, we denote it as R0. We can
rewrite the system in variational form (full derivation of this I'm leaving out):

dP/dt = A(R0) * P(t)
dO_k = O(t_k) - O(R0, t_k) = H * P(t_k) + e(t_k)

where P(t) = X(new_R0, t) - X(R0, t); X(new_R0, t) is a better estimation of the solution that we are trying to find; If we find P, we can get the better estimation from the current one, e.g. X(new_R0, t) = X(R0, t) + P or new_dR0 = R0 + P(t_0).
O(t_k) -- our observations, we already know them; so, dO(t_k) is known vector of
observation residuals; A(R0) is Jacobian matrix of Z(X), A(R0) = dZ/dX (when X=R0) (it doesn't
depend on time if the system of ODEs is stationary); H is Jacobian matrix of the measurement function;
in your case it would be very simple, H = [[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0]].T, since you're observing raw coordinates x and y, that are first two elements of X vector.

The latter system is linear and it can be solved using matrix exponential. So, we have

P(t_k) = EXP(A(R0)*(t_k-t0)) * P(t_0),

and, finally, substituting the latter into equation of measurements we get an
overdetermined system of linear equations:

dO_k = H * EXP(A(R0)*(t_k-t0)) * P(t_0) + e(t_k), k = 1, ..., N (assume that we have N measurements)

This is a system of linear algebraic equations (SLAE). It can be solved in least square sense, e.g.
using Moore-Penrose transfrom/inverse (np.linalg.pinv); Once we get its solution P(t_0) we can improve
our initial conditions: new_R0 = R0 + P(t_0) and repeat this process with new new_R0.

Depending on choice of R0 this process could converge to the solution we are trying to find.

This approach allows to get good estimation of initial conditions that leads to observed trajectory.
Also, if you try to estimate not initial X(t0) but current X(t) state of the system, then you probably invent the Kalman filter...

Hope that helps...

You can read about linear/nonlinear dynamical systems, e.g. in classical work "Topics in mathematical system theory" by Kalman, Falb, Arbib.

[madoko]

First scidam, thank you so much for your answer. It's so complete and precise !

I will just add a few precisions:

• Yes my ODE system is indeed stationary.
  
• My initial conditions are known.
  

I'm not sure I understand what you mean by:

Suppose that we have a good approximation of the solution, we denote it as R0.

The R0 here confuses me, especially in the case of this equation:

dO_k = O(t_k) - O(R0, t_k) = H * P(t_k) + e(t_k)

What you mean is that R0 = [a_approx,B_approx,G_approx], refining the approximation at each iteration ?

Again thanks very much !

[scidam]

Because I didn't know that initial conditions are known
I considered extended system, so R0 = [x0_approx, y0_appros, p0_approx, a_approx, B_approx, G_approx];

Therefore, dO_k -- residual between known measurement at time t_k and measurement at t_k obtained for trajectory started at (R0, t0).

From now, we know that (x0, y0, p0) are exactly known
(I wrote x0, y0, p0 because we have 3 coupled differential equations as you denoted that in the first post),
so our iterative process should be slightly changed:

There are several ways we can adopt our iterative process:

1) don't change x0, y0, p0 when iterating; apply changes to a, B, G only;
2) If you have measurement at t0, you can consider weighted least squares and assign very big weight to this measurement. Since measurements are raw coordinates
big weight for t0-measurement will account that x0 and y0 are very important measurements and shouldn't be changed over iterations;

Finally, everything is depending on the system of ODEs you are trying to solve. Is your system really hard to solve analytically?
You are likely needing to construct problem-specific algorithm to get estimations for a, B, G.

Have you some a priory information regarding a, B, G values...?
If you have fast numerical procedure for solving the ODE, you can try, e.g. simulated annealing approach.
At last, did you try the way proposed by Gribouillis?