NonlinearModelFit increase accuracy for 2 dimensional model

Question

I am trying to fit a 2 dimensional function. For that, I am using the function NonlinearModelFit. The model is

model = (a1[1, 0] - 3 x a1[2, 1] - 3/2 (x^2 + y^2) a1[3, 0] + 
15 x^2 a1[3, 2] - 3 y (5 y a1[3, 2] + b1[2, 1]) + 
30 x y b1[3, 2])^2 + 
1/4 (2 y (a1[2, 0] + 6 a1[2, 2] - 3 x (a1[3, 1] + 30 a1[3, 3])) + 
2 b1[1, 1] - 3 x (4 b1[2, 2] + x (b1[3, 1] - 30 b1[3, 3])) - 
9 y^2 (b1[3, 1] + 10 b1[3, 3]))^2 + (a1[1, 1] - 
9/2 x^2 (a1[3, 1] - 10 a1[3, 3]) - 
3/2 y (y a1[3, 1] + 30 y a1[3, 3] + 4 b1[2, 2]) + 
x (a1[2, 0] - 6 a1[2, 2] - 3 y b1[3, 1] + 90 y b1[3, 3]))^2

where a1[i,j] and b1[i,j] are the variables.

I am getting 10000 points for the fit from the function I want my model to be approximated to (just my free choice) in the form data={{x1,y1,f[x1,y1]},{x2,y2,f[x2,y2]},...}.

For the initial conditions,

InitialConditions = 
Table[{Variables[model /. {x -> 1, y -> 1}][[i]], 1*^-6}, {i, 1, 
Length[Variables[model /. {x -> 1, y -> 1}]]}]

And finally,

data = Import["https://pastebin.com/raw/wMLHJX9m", "Package"];
NonlinearModelFit[
data, {model, 
k - 1*^-6 < Abs[b1[1, 1]] < k + 1*^-6, -1*^-6 < Abs[a1[1, 0]] < 
1*^-6, -1*^-6 < Abs[a1[1, 1]] < 1*^-6}, InitialConditions, {x, y},
MaxIterations -> 10000, WorkingPrecision -> 9, 
Method -> {NMinimize, 
Method -> {"DifferentialEvolution", "ScalingFactor" -> 0.9, 
  "CrossProbability" -> 0.1, 
  "PostProcess" -> {FindMinimum, Method -> "QuasiNewton"}}}]

There is no reason for the choice of parameters, method and working precision that I am using, I have been trying different things and follow different examples.

The problem is that I am not getting good results. The behavior of the model function for the best fit parameters is not similar to the expected one at all. Also, each time I run the function, it takes around 20 min, so it is very difficult to try different parameters and method to compare.

I think that the problem might be that the function I want to approximate has the form $f(x,y)=f_0+h(x,y)$ with $h(x,y)$ of the order of $10^{-6}$ and $f_0 \approx 1$

EDIT:

This is a DensityPlot of my model with the fit parameters minus the exact function. As you can see, the functions don't even have the same behavior. I would also like to decrease this difference.

Answer 1

I think that if you subtract the minimum $y$ value and then multiply by 10^9 and remove all of the parameter restrictions, you'll get a good fit.

model = (a1[1, 0] - 3 x a1[2, 1] - 3/2 (x^2 + y^2) a1[3, 0] + 
      15 x^2 a1[3, 2] - 3 y (5 y a1[3, 2] + b1[2, 1]) + 
      30 x y b1[3, 2])^2 + 
   1/4 (2 y (a1[2, 0] + 6 a1[2, 2] - 3 x (a1[3, 1] + 30 a1[3, 3])) + 
       2 b1[1, 1] - 3 x (4 b1[2, 2] + x (b1[3, 1] - 30 b1[3, 3])) - 
       9 y^2 (b1[3, 1] + 10 b1[3, 3]))^2 + (a1[1, 1] - 
      9/2 x^2 (a1[3, 1] - 10 a1[3, 3]) - 
      3/2 y (y a1[3, 1] + 30 y a1[3, 3] + 4 b1[2, 2]) + 
      x (a1[2, 0] - 6 a1[2, 2] - 3 y b1[3, 1] + 90 y b1[3, 3]))^2;

InitialConditions = 
  Table[{Variables[model /. {x -> 1, y -> 1}][[i]], 1*^-6}, {i, 1, 
    Length[Variables[model /. {x -> 1, y -> 1}]]}];

data = Import["https://pastebin.com/raw/wMLHJX9m", "Package"];
data2 = data;
data2[[All, 3]] = ((data[[All, 3]] - Min[data[[All, 3]]]))*10^9;
nlm = NonlinearModelFit[data2, model, InitialConditions, {x, y}];
nlm["BestFitParameters"]

(* {a1[1, 0] -> -0.219553, a1[2, 1] -> -16.0184, a1[3, 0] -> 72.2813, 
 a1[3, 2] -> -56.9915, b1[2, 1] -> 1.85967, b1[3, 2] -> 3.48335, 
 a1[1, 1] -> 0.428318, a1[2, 0] -> 32.548, a1[2, 2] -> 7.37364, 
 a1[3, 1] -> -51.3845, a1[3, 3] -> -3.57784, b1[2, 2] -> 1.918, 
 b1[3, 1] -> 183.339, b1[3, 3] -> 17.3889, b1[1, 1] -> -0.361997} *)

Show[ListPointPlot3D[data2, PlotStyle -> PointSize[0.02], 
  BoxRatios -> {1, 1, 1}],
 ListPointPlot3D[
  Transpose[{data[[All, 1]], data[[All, 2]], 
    nlm["PredictedResponse"]}], PlotStyle -> {Red, PointSize[0.02]}]]