Infinite grid of planets with Newtonian gravity

In order to say that a configuration is in an equilibrium state, you need to say that the net force acting on each planet is zero. However, in order to define net force (if you're a mathematician), you need to choose and justify a rule for regularizing a sum that is not absolutely convergent. The rule that you seem to have chosen takes a limit over concentric balls centered at the planet in question. This method gives a convenient answer of zero net force.

You might ask why this is preferable to some other regularization method. For example, if you chose a limit over increasingly large ellipsoids of equal eccentricity, for which one focus is the planet in question, and the other focus moved away in a certain direction in the plane, you'd get an attractive force in that direction. This is a valid summation method, in the sense that the contributions from all of the planets are eventually considered when computing the net force on a planet.

One advantage of the spherical regularization over the ellipsoidal method is that it is invariant under rigid motions, and in fact any regularization method that is invariant under reflections or rotations will give you an answer of zero. This makes the equilibrium an aesthetically natural answer, but that doesn't necessarily mean it gives the correct answer. In particular, I think the setup is sufficiently separated from realistic physical situations that a preference for symmetry doesn't justify a regularization, so I'd say that this question does not have a well-defined answer (but physicists might disagree).

This is a very subtle limit of general relativity in which the evaluation of indefinite forms becomes important.

First, note that because of the uniform matter density, the gravitational potential won't be well-behaved. It would have to satisfy $$\Delta \Phi \sim \rho_{\rm planets} $$ but because the right hand side is uniform (at distances longer than the grid spacing), you would need something like $$\Delta = C\vec x\cdot \vec x$$ That's bad because the gravitational potential clearly goes to infinity at infinity and there's only one place at which the gravitational potential is stationary. All other planets would feel force. However, you may also ignore the gravitational potential and consider its gradient only. Equivalently, you may add a Newtonian "cosmological constant term" by $$\Delta \Phi \sim \rho_{\rm planets} -\rho_0 $$ Here, the cosmological constant term $\rho_0$ may be chosen to cancel the average value of $\rho_{\rm planets}$. This extra addition of $\rho_0$ wouldn't affect the accelerations and forces because it's completely uniform.

So if you insist on the standard formulae for the gravitational potential and you ban any kind of cosmological constant, this setup can't be in equilibrium. However, if you allow a constant vacuum energy to be added and to subtract the average density coming from the planetary grid, you may get a stationary state similar to a non-relativistic version of the Einstein static Universe that Einstein wanted to believe in general relativity (which was proved wrong when the expansion of the Universe was observed by Hubble).

Just like the Einstein static universe, this uniform configuration of matter is unstable. A small deviation will start to destroy the grid and make the planets clump, converting a part of their potential energy to kinetic energy and chaotic motion.

By the way, if you didn't want to use the gravitational potential at all, you could just compute the accelerations. Your reasoning - ending up with an "equilibrium" answer - is based on symmetry. The acceleration from the other planets is proportional to the integral $$ \int \frac{\vec r}{r^3} r^2 dr $$ written in the spherical coordinates. By symmetry, you could argue that it vanishes as the integrand is an odd function. However, that's only true around an origin that has to be chosen differently for each planet-probe. If you shift the center of the spherical coordinates, you will inevitably get a different result. This problem is mathematically isomorphic to the "anomalies" that arise from linear divergences in quantum field theory. In some sense, the theory is ill-defined with your planetary grid.

Well, the main question about a 2D grid in a 3+1 dimensional space-time has already been answered by Lubos Motl and Scott Carnahan with the main conclusion that both the potential sum and the force sum are not absolutely convergent, and therefore not mathematically well-posed.

Concerning the bonus question about a grid of dimension $d$, clearly a 3D grid would only make things worse, which leaves us with a 1D grid, say, a single row of planets along the $x$-axis at position $x\in a\mathbb{Z}$, where $a$ is the grid size. The corresponding potential sum is logarithmically divergent, but the force sum is actually absolutely convergent, so here the calculation makes rigorous sense. If we call the unit-vector in the $x$-direction for $\hat{x}$, then the gravitational force on the planet $P$ at $\vec{r}=\vec{0}$ is

$$\vec{F}~=~\frac{Gm^2}{a^2}\hat{x}\sum_{n\in\mathbb{Z}\backslash\{0\}}\frac{n}{|n|^3}~=~\vec{0},$$

which is zero by $n\leftrightarrow -n$ symmetry, as expected.