Given $d$ generators $G = \{ g_1, g_2, \cdots, g_d \} \subset \mathbb{R}^2_+$. One can obtain the Zonotope via $Z = CH(\{ \sum_{g \in S} g|S \subset G \})$, here $CH$ is convex hull operation.
And the extreme points of $Z$ are also called its vertices. In the 2-dimension, we know the number of vertices of $Z$ is $2d$.
And for some reason, all the 2+ dimensional Zonotope don't have explicit forms for number of vertices.
I've run some simulation for $G' \subset \mathbb{R}^3_+$, and the number of vertices, at least for $d\leq 20$ is $d^2-d+2$. I'm wondering why there is no general form for the number of vertices based on the number of generators.
And I'm wondering does $d^2-d+2$ generalize to $d \ge 20$ cases in $\mathbb{R}^3_+$ ?
The number of vertices of a zonotope is not a function of the number of generators. As a simple example in three dimensions, consider what happens if $G$ consists of three vectors on a plane, and a fourth not on that plane. The generated zonotope will be a hexagonal prism - with 12 vertices.
However, if you chose the generators to be in general position (i.e. no three on a plane), you would get a rhombic dodecahedron with 14 vertices instead. You can imagine that if you took four generators in general position and moved one towards the plane of two others, a part of the dodecahedron would start to flatten - and, eventually, on either side of the polyhedron, three parallelogram faces would merge into a single hexagonal face, eliminating the vertex inside.
In general, the number of vertices in a zonotope generated by $g_1,\ldots,g_d$ is not only a function of the number of generators, but also their incidence structure (i.e. their (oriented) matroid structure - if that's a familiar branch of mathematics). If the generators are in general position, the number of vertices will be
$$2\sum_{k=0}^n{d - 1\choose n}$$
where $n$ is one less than the dimensions of the ambient space and $d$ is the number of generators. This is a polynomial of degree $n$ and agrees with your answer for $n=2$. However, if the generators are not in general position, there will be fewer vertices. Note that incidence structure is the only thing that matters - and it turns out that requiring the vectors to have positive components doesn't affect the possibilities at all.
You can sort of imagine generating vertices of a zonotope as follows: fix some linear function $L:\mathbb R^n\rightarrow \mathbb R$. We want to figure out on what part of the zonotope $L$ is maximized. This is easy: we just add together every generator $g_i$ on which $L(g_i)>0$, leave out those for which $L(g_i)<0$, and can choose what we want to do if $L(g_i)=0$ (in which case, we will get a whole edge or face of the zonotope maximizing $L$). This gives some sense of why we might care if some $g_i$ are coplanar - then we'd find behaviors like a linear function vanishing on three points that we would not otherwise see, and this absolutely affects the resulting zonotope.
If we merely want to count the number of vertices, this idea of generating vertices provides a useful perspective: we just need to count how many linear functions $L$ give different vertices. We can do this in a straightforwards way: consider the space $\hat {\mathbb R}^n$ of linear functionals on $\mathbb R^n$ and the planes $\hat g_1,\ldots, \hat g_d$ of linear functionals where $L(g_i)=0$ for the corresponding $g_i$. This is an arrangement of $d$ planes in an $n$ dimensional space. Each time we cross a plane, we change which vertex maximizes the considered linear functional - so the number of vertices is just the number of regions* in $\hat {\mathbb R}^n \setminus (\hat g_1 \cup \ldots \cup \hat g_d)$. One can, from this, derive the above formula for points in general position in a fairly straightforwards inductive argument - though remembering that the incidence structure of the generators can always decrease the count of vertices.
(*In matroid language, that these regions represent maximal covectors of the oriented matroid for the generators - they're the regions we get if we're allowed to specify the sign of each $L(g_i)$. Covectors are the more general concept where we can, for each $g_i$, specify whether we want $L(g_i)$ to be positive, negative, or zero)