Triangulate the Voronoi cell then write the integral as a sum over the triangles:
$$\int_{\Omega}\|P - Pi\|\,dP = \sum_{k=1}^{N}\int_{\Delta_k}\|P-P_i\|\,dP.$$
Write the integration over the triangle in barycebtric coordinates. Let the Jacobian of the transformation for triangle $k$ be $J_k$, and $|det(J_k)| = 2|Area_{\Delta_k}|$. Then $\int_{\Delta_k}f(P)\,dP = |det(J_k)|\int_0^1\int_0^{1-\beta}f(P(\alpha, \beta))\,d\beta\,d\alpha$, where $\alpha, \beta$ are two of the barycentric coordinates.
Now we need only evaluate (let $\vec{v}_i$ be the vertices of the current triangle): $$\int_0^1\int_0^{1-\beta}\|\vec{v}_1 + \alpha (\vec{v}_2 - \vec{v}_1) + \beta (\vec{v}_3 - \vec{v}_1) - P_i\|^2\,d\alpha\,d\beta$$
Expand this and integrate the polynomials.