Signal decomposition: getting phase with known amplitude and frequency of main components
Trying to estimate the phase directly from the following equation requires nonlinear optimization.
$$
y(t) = \sum_{i=0}^N{\alpha_i\cos(2\pi f_i t + \phi_i)}
$$
Since the amplitudes and frequencies are known you can transform your problem to
$$
y(t) = \sum_{i=0}^N{\cos(\phi_i) \alpha_i \cos(2\pi f_i t)} - \sum_{i=0}^N{ \sin(\phi_i) \alpha_i  \sin(2\pi f_i t)}
$$
If you have $M$ samples for $t=t_0,... t_{M-1}$ and define
 $$p_i=\alpha_i \cos(\phi_i)$$ $$q_i=- \alpha_i \sin(\phi_i)$$ as your new parameters. Then you have a classic least square problem with $2N$ parameters.
$$\bf{Y} = A X + \text{noise} $$
where
$${\bf{Y}}=[y(t_0),....,y(t_{M-1})]^T$$

$A$ is a $M\times (2N)$ matrix where 
$$A(j,i)=\alpha_i \cos(2\pi f_i t_{j}) ~~~\text{for}~~~ i=0,..,N-1$$
$$A(j,i)=\alpha_i \sin(2\pi f_i t_{j}) ~~~\text{for}~~~ i=N,..,2N-1$$
and
$${\bf{X}}=[p_0,....,p_{N-1}, q_0,....,q_{N-1}]^T$$
The unknown vector $X$ (and thus all $p_i$ and $q_i$) can be estimated by normal equation (see for https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)):

$$X=(A^{T}A)^{-1}A^{T}Y$$

After estimating $p_i$ and $q_i$ you can estimate $\phi_i$, with an ambiguity of $\pm 2\pi k$ by taking 4 quadrant inverse $\tan$:
$$\phi_i=\text{atan2}(-q_i,p_i)$$

Also the amplitude $a_i$ can be estimated:
$$a_i =\sqrt{p^2_i+q^2_i}$$