In short the two models presented are identical for the particular dataset.
"Intuition is sometimes tricky on fitting procedures." in the answer by @mikuszefsky is right on the money. So is the comment by @MarcoB about increasing the number of iterations.
The two fits are found just fine by NonlinearModelFit:
data = With[{a = 43.45582489316203`,
f = 94.92003941300389`, ϕ = 431.155471523826`},
SeedRandom[1234];
Table[{t, a Cos[2 π f t + ϕ] + RandomReal[{-0.1, 0.1}]},
{t, 13.439999656460714`, 13.479799655455281`, 0.0002}]];
m1 = a Cos[2 π f t] + b Sin[2 π f t];
m2 = a Cos[2 π f t + ϕ];
nlm1 = NonlinearModelFit[data, m1, {{a, -30}, {b, 30}, {f, 95}}, t, MaxIterations -> 1000]
nlm2 = NonlinearModelFit[data, m2, {{a, 40}, {f, 100}, {ϕ, 3.9}}, t]
Both models are essentially identical but one cannot tell that from looking at the ParameterConfidenceIntervalTable.
The estimated variances are the same:
nlm1["EstimatedVariance"]
(* 0.0037219 *)
nlm2["EstimatedVariance"]
(* 0.0037219 *)
The $AIC_c$ values are identical:
nlm1["AICc"]
(* -543.167 *)
nlm2["AICc"]
(* -543.167 *)
The predictions are identical:
ListPlot[{Transpose[{nlm1["PredictedResponse"],
nlm2["PredictedResponse"]}], {{-45, -45}, {45, 45}}}, Joined -> {False, True}]
For m1 the numerical instabilities can be also seen by examining the parameter correlation matrix:
nlm1["CorrelationMatrix"]//MatrixForm
When the values in the off-diagonal are this close to 1 or -1, that suggests that there will either be numerical instability in the estimates and/or over-parameterization of the model for the particular data set.
The correlation matrix for m2 is also not without issues:
nlm2["CorrelationMatrix"]//MatrixForm
Addition:
Given the values of the parameters for one model, one can find the parameters for the other model that give identical predictions. Consider the following slight deviation from the original definitions:
m1 = a Cos[2 π f t] + b Sin[2 π f t];
m2 = c Cos[2 π f t + ϕ];
(Only m2 is changed: a is replaced with c.) Then
{a, b} = c {Cos[ϕ], -Sin[ϕ]}
and



