I have $m$ oscillators, where $m$ is unknown, with periods $\vec p = p_1, p_2, \ldots, p_{m}$. Each of the oscillators $j$ has associated with it a vector of sine coefficients $\vec A_j$ and angle coefficients $\vec \phi_j$ of unknown length.
The oscillators combine at time $t$ to make
$$F(t, \vec p, \vec A, \vec \phi) = \sum_j^{m} \sum_{k}^{|A_j|} A_{jk} \sin \left(2\pi k \left[ \frac{t}{p_j} \bmod 1 \right] + \phi_{jk} \right)$$
Now, I unevenly sample $n$ times with noise from this function and to get a matrix with rows $(t_i, F(t_i), \epsilon_i)$ where $\epsilon_i$ is the error on the observation.
My question: how do I find the smallest $\hat p, \hat A,$ and $\hat \phi$ such that the approximation is minimal and further that $$\left (\hat F(t_i, \hat p, \hat A, \hat \phi) - F(t_i, \vec p, \vec A, \vec \phi)\right)^2 \leq \epsilon_i \; \text{for all } i = 1, \ldots, n?$$
My understanding is that when $m = 1$, people try tackling this problem by first solving for $p$ e.g. using the Lomb-Scargle periodogram, then seek $|\vec A|$ e.g. using auto-correlation (Baart 1982), and finally solve for $\vec A$ and $\vec \phi$ using least-squares. But I don't think this approach generalizes to multiple periods.
In the multiperiodic case, I think the problem is one of multi-objective constrained optimization, since we have three objectives (sparsest $\hat p$, sparsest $\hat A$, and minimal loss) under $n$ constraints (loss is smaller than the sampling error for all samples). Is that the case? Would some MCOP toolbox be able to find $\hat p, \hat A$ and $\hat \phi$ for this problem? Are there some known solutions to this problem?