How to calculate vertex normals on a mesh with non-planar polygons
If you're interested in vertex normals specifically, there's an easy answer even for non-planar polygons that avoids the question of defining what the exact surface is: for each vertex, calculate the normal of the plane formed by the two edges entering and leaving that vertex.
More formally, given vertices $\mathbf{v}_1, \mathbf{v}_2, \ldots \mathbf{v}_n$ with counterclockwise winding, define the normal at the $i$th vertex as:
$$
\mathbf{n}_i = (\mathbf{v}_i - \mathbf{v}_{i-1}) \times (\mathbf{v}_{i+1} - \mathbf{v}_i)
$$
(where the indices wrap around).
You can then proceed to accumulate normals calculated this way from all the faces that share a given vertex, as usual for smooth shading.
As has been discussed, there are multiple ways to define a smooth surface corresponding to a non-planar polygon, but any reasonable way of defining such a surface must converge to the normals as defined here near the vertices, or else the surface can not both be smooth (locally flat) and meet the straight-line edges between the vertices.
A caveat to this approach, though, is that it won't define a normal for collinear vertices (where the entering and exiting edges are parallel), since the cross product goes to zero there. If this is a problem, it might work to patch up such vertices by interpolating normals to them from the surrounding non-collinear vertices.