Model calibration with phase space data

There are several ways to do this. I am going to show a more-or-less straightforward way. Some adjustments probably have to be done when using real data.

General procedure

1. Using the data find an interpolation (or use a fit with appropriate basis) of y[x].

2. Define a function ParDist that measures the deviation of the fit. That function takes as arguments:

   • the found parametric functions,

   • a pair of parameters,

   • an interpolation function for y[x],

   • a vector of time points (to be used to compute the measure).

3. Select a set of time points for ParDist.

4. Minimize the measure function ParDist over an appropriate domain for the parameters.

5. Compute new data points with the found parameters over the selected time points.

6. Compare the new data with the initial data.

7. Repeat the steps 1-6 using different minimization methods, norms in ParDist, and sets of time points for ParDist.

It might happen that the interpolation in step 1 is not an adequate approach, and some fitting procedure has to be used. For simplicity I am using interpolation below.

Code

Interpolation function for y[x]:

xyFunc = Interpolation[data[[All, 2 ;; 3]]]

The measure function to be minimized, ParDist:

Clear[ParDist]
ParDist[{x_ParametricFunction, y_ParametricFunction}, {d_?NumberQ, v0_?NumberQ}, xyFunc_, ts : {_?NumberQ ..}] :=
  Block[{points, xMin, xMax},
   points = Map[{x[d, v0][#], y[d, v0][#]} &, ts];
   {xMin, xMax} = xyFunc["Domain"][[1]];
   points = Select[points, xMin <= #[[1]] <= xMax &];
   If[Length[points] == 0, 10^10.,
    Norm[Map[(xyFunc[#[[1]]] - #[[2]])/Abs[xyFunc[#[[1]]]] &, points]]
   ]
  ];

(The selection of the points to be within the interpolation domain of xyFunc is important.)

Select a set of time points (does not have to be a regular grid):

ts = Range[0.5, 3, 0.13];

Minimize ParDist over an appropriate domain for the parameters:

nsol = 
 NMinimize[{ParDist[{x, y} /. solPN, {d, v0}, xyFunc, ts], 0 <= d <= 20, 5 <= v0 <= 25}, {d, v0}]

(* {0.654588, {d -> 9.2491, v0 -> 19.9544}} *)

Compute new data points with the found parameters:

Block[{d, v0},
  {d, v0} = {d, v0} /. nsol[[2]];
  fitData = Table[{t, x[d, v0][t] /. solPN, y[d, v0][t] /. solPN}, {t, ts}]
];

Compare the plots of the fitted data and initial data:

ListPlot[{data[[All, 2 ;; 3]], fitData[[All, 2 ;; 3]]}, 
 PlotRange -> All, PlotTheme -> "Detailed", 
 PlotLegends -> {"initial", "fitted"}]