I just want to add to the discussion (without confusing the original post) that all the different routes will work if there is only one differential equation:
data = NDSolveValue[{
x''[t] - k1*(1 - x[t]^2)*x'[t] + k2*x[t] == 0,
x[0] == 2, x'[0] == 0} /. {k1 -> 1., k2 -> 1.},
Table[{t, x[t] + RandomReal[{-.3, .3}]}, {t, 0, 10, .2}], {t, 10}];
dataT = data\[Transpose];
ti = dataT[[1, All]];
di = dataT[[2, All]];
pfun = ParametricNDSolveValue[{
x''[t] - k1*(1 - x[t]^2)*x'[t] + k2*x[t] == 0,
x[0] == 2, x'[0] == 0},
x, {t, 0, 10}, {k1, k2}];
FindFit finds the best-fit parameters sucessfully:
fit = FindFit[data, pfun[k1, k2][t], {{k1, 2}, {k2, 0}}, t]
Out[1]= {k1 -> 1.09028, k2 -> 1.02729}
FindMinimum found the identical result:
FindMinimum[Total[(di - Map[pfun[k1, k2], ti])^2], {k1, k2}]
Out[2]= {1.41041, {k1 -> 1.09028, k2 -> 1.02729}}
And the Module approach produced the same result:
objfun[k_] := Module[{},
fun = NDSolveValue[{
x''[t] - k[[1]]*(1 - x[t]^2)*x'[t] + k[[2]]*x[t] == 0,
x[0] == 2, x'[0] == 0},
x, {t, 0, 10}];
Total[(di - Map[fun, ti])^2]
]
FindMinimum[objfun[{k1, k2}], {{k1, 1.0}, {k2, 0.0}}]
Out[3]= {1.41041, {k1 -> 1.09028, k2 -> 1.02729}}