Assume that $\Omega\subset R^2$ is an open bounded set with a smooth boundary, $g:\partial\Omega\to R$ is a continuous map and $\{b_i \ | \ i=1,2,\ldots,d\}$ is a finite subset of $\Omega$.
$$\left\{\begin{align*}& \Delta \Phi = 2\pi \sum_{i=1}^d \delta_{b_i} &&\textrm{in} \ \ \Omega \\ &\frac{\partial\Phi}{\partial\nu} = g &&\textrm{on} \ \ \partial\Omega\end{align*}\right.$$
I would like to ask, if maybe someone of you know how to prove the existence of solutions of the above problem. I tried with Green function but I don't know how to use it to obtain the specific Neumann boundary conditions.
You need an additional condition for there to be solutions. Think of this in terms of heat flow.
You're asking for a steady-state temperature distribution in $\Omega$ with given heat sources and given heat flows across the boundary. The rate of change of the total amount of heat in the region is then the sum of the sources plus the flux across the boundary. If that is not $0$, there can't be a steady state.