According to Petr-Douglas-Neumann theorem (Wikipedia article) after all n-2 steps (n-number of sides/vertices in the polygon) are complete, the result is a regular polygon whose centroid coincides with the centroid of the original polygon.
Trying to check with GeoGebra the theorem for a quadrilateral, applying angle $\pi$ (k=2, midpoints of sides) and then angle $\frac{\pi}{2}$ ($k=1$, base angles $\frac{\pi}{4}$).
As the result I really get a square, but its centroid (i.e. centre) doesn't match the centroid of the original quadrilateral. To find the centroid of the original quadrilateral I use method suggested at another math.stackexchange thread.
BTW, when I connect the opposite midpoints of the original quadrilateral the intersection does coincide with the centre of constructed square. However, as far as I understand it, the point obtained by connecting opposite midpoints isn't in general a centroid.
At the attached image: $A_1 A_2 A_3 A_4$ — original quadrilateral, $B_1 B_2 B_3 B_4$ — after first step (applied angle $\pi$), $C_1 C_2 C_3 C_4$ — the result, $N$ — centroid of the original quadrilateral, $O$ — centroid of the resulting square.
GeoGebra file: https://drive.google.com/file/d/19BOvzUcFPmxl_d6zcKE6Bi-7yxmpItGC/view?usp=sharing
In the context of Petr-Douglas-Neumann (PDN), "centroid" means "vertex-centroid" (ie, the average of the vertices). To see why the preservation of the vertex-centroid is "obvious" (and the least-interesting thing about PDN), let's consider a $n$-gon $p_0p_1p_2\cdots p_{n-1}$ with vertices $p_k$ in the complex plane.
Each step of the PDN process effectively applies a common operation to rotate-and-scale each side $p_kp_{k+1}$ about vertex $p_k$; we then shift focus to consider the polygon whose vertices are the images of the $p_{k+1}$. Each such image is given by $$p_k + t (p_{k+1}-p_k) = (1-t) p_{k}+ tp_{k+1} \tag1$$ for some complex $t := s e^{i\theta}$, where $s$ is the (real) scale factor and $\theta$ the (real) rotation angle. (The exact values of, or relations between, $s$ and $\theta$ don't matter for this discussion.) We can represent this effect across the entire polygon by using matrices. Defining $$P := \begin{bmatrix}p_0 & p_1 & p_2 & \cdots & p_{n-1}\end{bmatrix} \qquad\text{and}\qquad T_t := \begin{bmatrix} 1-t & 0 & 0 & \cdots & t \\ t & 1-t & 0 & \cdots & 0 \\ 0 & t & 1-t & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1-t \end{bmatrix} \tag2$$ we have that the entries of $PT_t$ are the vertices of the resulting polygon; likewise, the entries of $PT_t T_u$, $PT_t T_u T_v$, etc, are the vertices of polygons from further iterations of the process.
Now, the vertex-centroid of a polygon is the average of its vertices; ie, $$k := \frac1n(p_0+p_1+p_2+\cdots+p_{n-1}) = \frac1n P \,\mathbf{1}_n \tag3$$ where $\mathbf{1}_n$ is the "all-$1$s" column vector. But note: $T_t\,\mathbf{1}_n$ is just $\mathbf{1}_n$ again, so the vertex-centroid of the polygon obtained from applying the $T_t$ operation to $P$ is $$\frac1n(P\,T_t)\,\mathbf{1}_n = \frac1n P\,(T_t\,\mathbf{1}_n) = \frac1n P\,\mathbf{1}_n = k \tag4$$ That is to say: the vertex-centroid of the original polygon $P$ is preserved under the $T_t$ operation; and clearly, it's preserved under every $T$-operation applied thereafter. The thing never moves! $\square$
Of course, centroid-preservation is a worthwhile property to mention, but its significance pales in comparison to the fact that the PDN process ultimately results in a (standard, convex) regular $n$-gon from any arbitrary $n$-gon, thus generalizing Napoleon's theorem for triangles. Also interesting: (1) the order of the steps doesn't matter; and (2) if we replace a step with one based on that "final" $n$-gon, then the new "final" figure corresponds to the step we replaced! (I'm being a bit loosey-goosey with the description here. In any case, showing these aspects requires getting into the details of the various $s$- and $\theta$-values, so I'll leave that as an exercise to the reader.)