Something to investigate: if I enter breaks=1, I get a result that says:
Optimal 1-segment partition:
Call:
breakpoints.formula(formula = formula(response ~ trend + harmon),
h = 46, breaks = 1, data = bfastpp(bfts))
Breakpoints at observation number:
NA
Corresponding to breakdates:
NA
That is weird, because breaks = 1 means a 2-segment partition, not a 1-segment partition. breaks=2 works as expected.