Fitting data using ParametricNDSolve

Question

I have a system of differential equations that are solved using ParametricNDSolve, leaving two variables undefined (m1 and m2) so that they can be fitted to data later. However, Nonlinearmodelfit does not seem to like the function I give it. Specifically, it says it's not a list of real numbers. My guess is that I am misunderstanding what the output of ParametricNDSolve actually is. Before, m1 and m2 had actual values attached to them and I simply used ndsolvevalue, got some interpolating functions for xbh1, xbh2, ybh1, ybh2, used them in my fitting function and got a solid fit. Any thoughts?

lc = {{0.005599829, -0.024633333333333396}, {0.16799487, \
-0.01863333333333339}, {0.33598974, -0.014633333333333387}, \
{0.515184268, -0.00663333333333338}, {0.705578454, 
    0.005366666666666575}, {0.890372811, 
    0.014366666666666583}, {1.058367681, 
    0.017366666666666586}, {1.259961525, 
    0.019366666666666588}, {1.45595554, 
    0.02336666666666659}, {1.6519495549999998, 
    0.019366666666666588}, {1.8199444249999999, 
    0.009366666666666579}, {2.010338611, 
    0.003366666666666629}, {2.228731942, -0.013633333333333386}, \
{2.419126128, -0.01763333333333339}, {2.603920485, \
-0.022633333333333394}, {2.794314671, -0.026633333333333398}, \
{2.979109028, -0.015633333333333388}, {3.1639033849999993, \
-0.013633333333333386}, {3.3542975709999996, -0.0006333333333333746}, \
{3.516692612, 0.005366666666666575}, {3.723886285, 
    0.014366666666666583}, {3.9254801289999994, 
    0.02236666666666659}, {4.082275341, 
    0.02136666666666659}, {4.278269356, 
    0.02036666666666659}, {4.463063713, 
    0.013366666666666582}, {4.631058583, 
    0.005366666666666575}, {4.810253111, -0.0006333333333333746}, \
{5.011846955, -0.007633333333333381}, {5.196641312, \
-0.013633333333333386}, {5.420634472, -0.01763333333333339}}; (*data*)

error = {0.0024, 0.0028, 0.0028, 0.0027, 0.0023, 0.0025, 0.0023, 
   0.002, 0.0022, 0.0022, 0.0032, 0.0022, 0.0025, 0.002, 0.0022, 
   0.0041, 0.0021, 0.0019, 0.002, 0.0021, 0.0029, 0.0026, 0.0023, 
   0.0022, 0.0022, 0.0028, 0.0029, 0.002, 0.0024, 0.0019}; (*error*)



G = 6.6743*10^-11;
    c = 2.99792458*10^8;
    Periods = 5.599829;
    P = Periods*86400;
\[Omega] = 5.3511795
    ecc = 0.018;
    a = Power[(P^2*G*(m1 + m2))/(4*\[Pi]^2), (3)^-1]
    \[Mu] = (m1 m2)/(m1 + m2);
    aph = a (1 + ecc);
    x0 = (m2*aph)/(m1 + m2)
    y0 = 0;
    X0 = -((m1*aph)/(m1 + m2))
    Y0 = 0;
    Vel1 = ((G \[Mu] )/a  m2/m1 (1 - ecc)/(1 + ecc) )^(
     1/2)(*initial velocity star 1*)
    Vel2 = -((G \[Mu] )/a  m2/m1 (1 - ecc)/(1 + ecc) )^(1/2) m1/
      m2(*initial velocity star 2*)
    
    
    system1 = {
       x1''[t] == -((G m1 m2  (x1[t] - x2[t]))/(
         m1 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2))),
       y1''[t] == -((G m1 m2  (y1[t] - y2[t]))/(
         m1 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2))),
       x2''[t] == (G m1 m2  (x1[t] - x2[t]))/(
        m2 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2)),
       y2''[t] == (G m1 m2  (y1[t] - y2[t]))/(
        m2 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2))};
    
    initials = {
       x1[0] == x0,
       x1'[0] == 0,
       y1[0] == 0,
       y1'[0] == Vel1,
       x2[0] == X0,
       x2'[0] == 0,
       y2[0] == 0,
       y2'[0] == Vel2};
    
    {xbh1, ybh1, xbh2, ybh2} = 
     ParametricNDSolve[{system1, initials}, {x1[t], y1[t], x2[t], 
       y2[t]}, {t, 0, 6.0*10^7}, {m1, m2}]
    bh1 = {(xbh1[t*86400]/(1.5*10^11)), (ybh1[
         t*86400]/(1.5*10^11))};(*larger star position vectors*)
    bh2 = {(xbh2[t*86400]/(1.5*10^11)), (ybh2[
         t*86400]/(1.5*10^11))};(*smaller star position vectors*)

\[CapitalGamma] = 0.4;
g = 0.4;
\[Alpha] = (0.15 (15 + \[CapitalGamma]) (1 + g))/(3 - \[CapitalGamma])

sev = \[Alpha] m1/
  m2 (rstar/\[Sqrt]((xbh2[t*86400] - 
         xbh1[t*86400])^2 + (ybh2[t*86400] - ybh1[t*86400])^2))^3 Cos[
   2 ((2 \[Pi])/P*(t*86400) + \[Pi] + \[Omega] + f)] Sin[\[Pi]/
    2 - \[Phi]]^2; (*function to be fitted*)

myfit = NonlinearModelFit[lc, 
   Re[0 - 2.5 Log10[((1 + (sev)))]], {{m1, 20 (1.988*10^30)}, {m2, 
     40 (1.988*10^30)}, {rstar, 22 (6.957*10^8)}, {f, 0.7}, {\[Phi], 
     1}}, {t}, Weights -> 1/error^2, ConfidenceLevel -> 0.99, 
   WorkingPrecision -> 30];