How to estimate position from noisy angles?

Each of the $\ i=1,2,...,n\ $ bearing measurements can be expressed as, $$ \alpha_i = \angle\delta_i - \alpha_0 $$ where $\alpha_0$ is the unknown angle from $x_1$ to $\hat{x}_1$ and $$ \delta_i := P_i-P_0\\ \angle\delta_i = \tan^{-1}({\delta_i}_2\ /\ {\delta_i}_1) $$

If we define, $$ q := \begin{bmatrix} P_0 \\ \alpha_0 \end{bmatrix}\ \ \ \ \ \ z := \begin{bmatrix} \alpha_1 \\ \alpha_2 \\ \vdots \\ \alpha_n \end{bmatrix} $$ and assume our bearing sensor is subject to additive noise then we can write, $$ z = h(q) + \epsilon $$ as our complete measurement model, where $\epsilon$ is a random variable with some known characteristics, like perhaps zero mean and covariance $R$.

It would also be nice / common to have some "prior" probabilistic information about our state $q$. That is, even if we don't know it exactly, we do have a hunch about what it could be, described by a probability distribution $q \sim \rho(q)$.

Unfortunately, I don't see a way to reformulate $h(q)$ to be linear in $q$, so I can't just provide you with the closed-form optimal solution: i.e. the $q$ that maximizes the conditional probability distribution $\rho(q|z)$. As far as I can tell, we can only approximate this solution using the general approaches of nonlinear estimation theory.

I implemented an extended Kalman filter as a first attempt. This involves a linearization of $h(q)$, so I computed its Jacobian by hand, $$ \frac{dh_i}{dq} = \begin{bmatrix} \frac{{\delta_i}_2}{||\delta_i||^2_2} & \frac{{-\delta_i}_1}{||\delta_i||^2_2} & -1\end{bmatrix} $$

Another thing to note is that since we are dealing with angles, our measurements (and part of the state) are really members of $\mathbb{SO}2$ not $\mathbb{R}$. Thus we must be careful to work in the Lie algebra, which for this simply means "unwrapping" any angle differences onto $[-\pi, \pi]$.

To fully realize your problem, I have assumed that the additive bearing sensor noise is Gaussian with zero mean and a standard deviation of $5$ degrees. I arbitrarily placed $7$ landmarks. The animation below shows the EKF converging. I start with an incorrect state guess and a large covariance to reflect my prior uncertainty. You can also see that the bearing measurements (yellow rays) don't quite hit their respective landmarks (green dots). Obviously the "process model" in this case is stationary (the true state isn't moving). I linearize about the first guess, compute the best linear estimator (i.e. EKF corrector step), then relinearize about the updated mean and "reset" my covariance back to the prior before attempting another linearization & update. To be clear, only $7$ measurements total were taken (one for each landmark), and all $7$ are incorporated at the same time - the "iterations" are just relinearizations. The code is here. If you have any questions feel free to ask.

Here is one more visualization this time giving the bearing sensor rediculous noise - a standard deviation of $30$ degrees (you can see how way-off the yellow rays are). The filter still does a good job, with its posterior covariance reasonably reflecting the direction along which we still have uncertainty. Note that heading always converges well because the sensor model is actually linear in that state.