In this answer we analyze a generalization of OP's question. In the last section 6 we will argue heuristically that one should expect an upper bound of critical points of the form
$$c(M)~\leq~~5n-11 \qquad {\rm for} \qquad n\geq 3. $$
1) Let us identify the plane $\mathbb{R}^2\cong \mathbb{C}$ with the complex plane $z=x+\mathrm{i} y$. Let
$$Z:=\{z_1, \ldots, z_n\}\subseteq \mathbb{C}$$
be a set of $n$ different punctures in the complex plane, where $n\in\mathbb{N}$. Consider the punctured complex plane $M:=\mathbb{C} \backslash Z$ and the Riemann sphere $S^2:=\mathbb{C} \cup \{\infty\}$ with Betti numbers
$$b_0(M)~=~1, \qquad b_1(M)~=~n, \qquad b_2(M)~=~0,$$
$$b_0(S^2)~=~1, \qquad b_1(S^2)~=~0, \qquad b_2(S^2)~=~1,$$
and Euler characteristics
$$\chi(M)~=~b_0(M)-b_1(M)+b_2(M)~=~ 1-n, $$
$$\chi(S^2)~=~b_0(S^2)-b_1(S^2)+b_2(S^2)~=~ 2, $$
respectively.
2) Let $p>0$ and $k_1, \ldots, k_n >0$ be $n+1$ positive constants.
Let the potential $V:M \to]0,\infty[$ and its extension $\tilde{V}: S^2 \to [0,\infty]$ be
$$V(z)~:=~\sum_{i=1}^n \frac{k_i}{p|z-z_i|^p},$$
$$ \tilde{V}(z)~:=~\left\{ \begin{array}{rcl} V(z) &{\rm for}& z \in M, \cr +\infty &{\rm for}& z \in Z, \cr 0 &{\rm for}& z \in\{\infty\}. \end{array} \right. $$
Let $$ c_0~:=~ \#{\rm minimum~pts}, \qquad c_1~:=~ \#{\rm saddle~pts}, \qquad c_2~:=~ \#{\rm maximum~pts}, $$
So
$$c_0(S^2)~=~c_0(M)+1, \qquad c_1(S^2)~=~c_1(M), \qquad c_2(S^2)~=~c_2(M)+n, $$
because $z=\infty$ is a minimum point and $Z$ are maximum points for $\tilde{V}$.
3) Define two positive functions $E, F:M \to\mathbb{R}_{+}$ as
$$E(z)~:=~\sum_{i=1}^n\frac{ k_i(x-x_i)^2}{|z-z_i|^{p+4}}~>~0,\qquad
F(z)~:=~\sum_{i=1}^n\frac{ k_i(y-y_i)^2}{|z-z_i|^{p+4}}~>~0, \qquad z \in M. $$
The $2\times 2$ Hessian matrix $H$ for the potential $V$ is
$$H_{xx}~=~\sum_{i=1}^n k_i\frac{(p+1)(x-x_i)^2-(y-y_i)^2}{|z-z_i|^{p+4}}
~=~(p+1)E-F,\qquad z \in M, $$
$$H_{yy}~=~\sum_{i=1}^n k_i\frac{(p+1)(y-y_i)^2-(x-x_i)^2}{|z-z_i|^{p+4}}
~=~(p+1)F-E,\qquad z \in M, $$
$$H_{xy}~=~(p+2)\sum_{i=1}^n \frac{k_i(x-x_i)(y-y_i)}{|z-z_i|^{p+4}}, \qquad z \in M,$$
with positive trace
$$ {\rm tr}H ~=~H_{xx}+H_{yy}~=~p(E+F)
~=~p\sum_{i=1}^n \frac{k_i}{|z-z_i|^{p+2}}~>~0, \qquad z \in M.$$
and determinant
$$ {\rm det}H ~=~H_{xx}H_{yy}-H^2_{xy}
~=~\underbrace{(p^2+2p+2)EF}_{>0}
- \underbrace{\left((p+1)(E^2+F^2)+H^2_{xy}\right)}_{>0}.$$
A sufficient condition for negative determinant is:
$$F> (p+1)E \qquad \vee \qquad E> (p+1)F\qquad\Rightarrow \qquad {\rm det}H<0. $$
The off-diagonal element $H_{xy}$ is bounded by
$$ \frac{|H_{xy}|}{p+2} ~\leq~\sum_{i=1}^n \frac{k_i|x-x_i||y-y_i|}{|z-z_i|^{p+4}}
~\leq~\frac{E+F}{2}. $$
4) The positive trace ${\rm tr}H>0$ implies that there cannot be any maximum points in $M$,
$$c_2(M)~=~0 \qquad \Leftrightarrow \qquad c_2(S^2)~=~n .$$
A minimum point $z$ must have positive determinant ${\rm det}H(z)>0$, while a saddle point $z$ must have non-positive determinant ${\rm det}H(z)\leq 0$. Generically, one may assume that all critical points $z$ in $M$ are non-degenerate ${\rm det}H(z)\neq 0$, so that $V$ is a Morse function. This implies that all critical points in $M$ are isolated points. Morse theory yields that$^{1)}$
$$ c_2(M)-c_1(M)+c_0(M)~=~\chi(M)\qquad \Leftrightarrow \qquad
c_2(S^2)-c_1(S^2)+c_0(S^2)~=~\chi(S^2),$$
or equivalently,
$$ c_1(M)-c_0(M)~=~n-1\qquad \Leftrightarrow \qquad c_1(S^2)-c_0(S^2)~=~n-2.$$
5) Lemma: The determinant ${\rm det}H(z)<0$ is negative sufficiently close to the set $Z$.
The open set
$$ U := \{z\in M \mid {\rm det}H(z)>0\} ~=~\cup_a U_a $$
consists of a number of connected components $U_a$. Since the restriction $V|_{U}$ is concave, each connected component $U_a$ can contain at most one minimum point.
6) Up until now our analysis is rigorous. In the rest of the answer we speculate. Heuristically, there should exists a triangularization of the Riemann sphere $S^2$ with the set $Z$ as vertices, such that each face contains at most one minimum point.
Assume that this is true. Let the number of faces, edges, and vertices be $f$, $e$, and $n$, respectively. Assume $n\geq 3$. Then
$$ f-e+n ~=~\chi(S^2)~=~ 2, \qquad 2e~=~\sum_{j\geq 3} j f_j~\geq~3f, \qquad f~=~\sum_{j\geq 3} f_j, $$
implies that $\frac{f}{2}~=~\frac{3f}{2}-f~ \leq~ e-f~=~n-2$, and hence,
$$c_0(M)+1~=~c_0(S^2)~\leq~f~\leq~ 2(n-2)\qquad \Rightarrow \qquad c_0(M)~\leq~ 2n-5.$$
Then an upper bound for the number of critical points in $M$ becomes
$$c(M)~=~ c_{0}(M) + c_{1}(M) + c_{2}(M) ~=~ 2c_{0}(M)+n-1~\leq~ 5n - 11. $$
--
$^{1)}$ Naively, the Morse inequality implies that $c_{0}(M) \geq b_{0}(M)=1$, which is obviously wrong for $n=1$. This can be cured, if one takes boundary contributions in $M$ properly into account. The formulation in terms of the Riemann sphere $S^2$ doesn't have this problem, because it has no boundary.