It turns out that the integral evaluation within ParametricNDSolve is slow due to symbolic processing. I removed it from ParametricNDSolve and evaluated the integral with NIntegrate turning SymbolicProcessing->0. It was >100x faster on my machine doing it this way. I also like using non-dimensional scaling to keep values closer to 1.
I use the subscript $d$ to relate parameters with dimension to their dimensionless counterparts as shown below.
$$\begin{gathered}
L = \frac{{{L_d}}}{{{L_{\max }}}} \hfill \\
a = \frac{{{a_d}{L_{\max }}{t_{\max }}}}{{{n_0}}} \hfill \\
b = {b_d}{n_0}{t_{\max }} \hfill \\
c = {c_d}n_0^2{t_{\max }} \hfill \\
{n_0} = \frac{1}{{\sqrt {K{t_{\max }}} }} \hfill \\
\end{gathered} $$
Since the data appear fairly evenly spread out on a log-log scale, I will also do the fit on a log-log scale to achieve a better weighting. Here is how I recast the example given in the OP.
data = {{0, 0}, {3.69 10^12, 50.65}, {5.67 10^12,
134.875}, {7.05 10^12, 225.275}, {9.03 10^12,
381.65}, {1.09 10^13, 509.5}, {1.25 10^13, 595.3}, {1.45 10^13,
815.325}, {1.76 10^13, 1225}, {2.10 10^13, 1624.125}, {2.46 10^13,
2018.725}, {2.84 10^13, 2488.775}, {3.18 10^13,
2942.9}, {3.68 10^13, 3630}, {4.39 10^13, 4558.65}, {5.52 10^13,
5925.925}, {6.45 10^13, 7044.075}, {7.23 10^13,
7972.2}, {8.18 10^13, 9119.575}, {9.38 10^13, 10545}, {1.06 10^14,
11749.1}, {1.28 10^14, 13760.475}, {1.42 10^14,
15055.65}, {1.57 10^14, 16484.475}};
(* Scale data max L value *)
sdata = data;
sdata[[All, 1]] = data[[All, 1]]/Max@data[[All, 1]];
(* Data looks equally spaced on log-log scale *)
(* so let's scale it (Note: remove {0,0}) *)
lldata = N@Log@sdata[[2 ;; -1]];
(* create parameters for non-dim scaling *)
tmax = 30;
Kfac = 1.04 10^-36;
scalefac = 1;
Lmax = Max@data[[All, 1]];
n0 = 1/Sqrt[Kfac tmax];
(* Scaled equations *)
eqns = {n'[t] ==
tmax/n0 (22500 Lmax L - 10^-22 n0^2 n[t]^2 -
2.0 10^-40 n0^3 n[t]^3) , n[0] == 0};
pfun = ParametricNDSolveValue[ eqns, n, {t, 0, 1}, {L}];
f = Function[{l},
Log@NIntegrate[(pfun[Exp[l]][t])^2, {t, 0, 1},
Method -> {Automatic, "SymbolicProcessing" -> 0}], Listable];
Show[Plot[f[l], {l, -4, 0.5},
AxesLabel -> {l, K \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(30\)]\(n[l,
t]^2 \[DifferentialD]t\)\)}, ImageSize -> Large,
PlotStyle -> Red], ListPlot[lldata, PlotStyle -> Black],
Plot[Log[0.45 10^-23*Lmax^2*Exp[l]^2], {l, -4, 0.5}], LabelStyle ->
Directive[Black, Bold, Medium]]
(* Show scaled equations *)
eqns // Simplify
(* {(n^\[Prime])[t]\[Equal]591.9441354553654` \
L-0.0005370861555295748` n[t]^2-0.00019230769230769236` n[t]^3,n[0]\
\[Equal]0} *)

It looks we have recast the problem to the mimic the original state. Now, we can formulate the problem to fit $a,b,c,sf$ just using FindFit out-of-the-box.
(* Dimensionless Equations *)
eqns = {n'[t] == a L - b n[t]^2 - c n[t]^3 , n[0] == 0};
pfun = ParametricNDSolveValue[ eqns, n, {t, 0, 1}, {a, b, c, L}];
ffull = Function[{a, b, c, sf, l},
sf Log@NIntegrate[(pfun[a, b, c, Exp[l]][t])^2, {t, 0, 1},
Method -> {Automatic, "SymbolicProcessing" -> 0}], Listable];
(* rewrite in more typical model form *)
model[a_, b_, c_, sf_][l_] := ffull[a, b, c, sf, l]
{time, fit} =
Quiet@FindFit[lldata,
model[a, b, c, sf][
l], {{a, 500}, {b, 0.005} , {c, 0.0005}, {sf, 1}}, l] //
AbsoluteTiming
(* {3.18,{a->440.,b->0.00719,c->0.000705,sf->1.13} *)
Show[Plot[model[a, b, c, sf][l] /. fit, {l, -4, 0.5},
AxesLabel -> {Log@l, Log[sf \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(1\)]\(n[l,
t]^2 \[DifferentialD]t\)\)]}, ImageSize -> Large,
PlotStyle -> Red], ListPlot[lldata, PlotStyle -> Black],
Plot[Log[0.45 10^-23*Lmax^2*Exp[l]^2], {l, -4, 0.5}], LabelStyle ->
Directive[Black, Bold, Medium]]
(* Rescale Parameters Back to Dimensioned Forms *)
nd = {a, b, c, L};
scalefactors = {Lmax tmax/n0, n0 tmax, n0^2 tmax, 1/Lmax};
(nd /. fit) ~Times~(1/scalefactors)
(* {16700.,1.34*10^(-21),7.33*10^(-40),(1.57*10^(14)) L} *)

The fit looks pretty good and it took about 3 seconds on my machine.
pfun
depends on 5(!) parameters{a,b,c,k,L}
. The second argument ofFindFit
should be something likek*pfun[a, b, c, k, L]
$\endgroup$