Building up from @Dan's comment, the problem seems friendly enough that you only seem to need the equations of motion at time $0$ to find a geometry that follows the desired rules. In reality, things turned out to be trickier.
First, we start by setting the centroid as the origin of coordinates. Then we place charge $1$ at $\mathbf{r_1}\left( t = 0\right) = (1, 0)$, charge $2$ at $\mathbf{r_2}\left( t = 0\right) = (\cos(\phi),\sin(\phi))$ and charge $3$ at $\mathbf{r_3}\left( t = 0\right) = - \mathbf{r_1}\left( t = 0\right) - \mathbf{r_2}\left( t = 0\right) $.
For the condition established in the problem to be satisfied, the net force acting on charge $1$, $\mathbf{f_1}$, must remain horizontal at all times. In fact, all net forces must remain parallel to the corresponding position vectors: $\mathbf{f_i} \parallel \mathbf{r_i}$. In particular this is true at $t=0$. The equation $0 = (0, 1) \cdot \mathbf{f_1}\left( t = 0\right)$ happens to be enough to give you a value for $\phi$, which, in this setup, fully determines the geometry of the system.
After replacing this value and letting the system evolve, it turns out that this is not a general enough ansatz to have a triangle that genuinely remains symmetrical at all times. The initial position of charge $2$ should be less constrained that in my previous ansatz, so we go instead for
\begin{eqnarray}
\mathbf{r_1}\left( t = 0\right) &=& (1, 0) \, ,
\\
\mathbf{r_2}\left( t = 0\right) &=& l(\cos(\phi),\sin(\phi)) \, ,
\\
\mathbf{r_3}\left( t = 0\right) &=& - \mathbf{r_1}\left( t = 0\right) - \mathbf{r_2}\left( t = 0\right) \, ,
\end{eqnarray}
with $l$ some positive number. We exhausted the degree of freedom associated with scaling invariance by setting the norm of $\mathbf{r_1}\left( t = 0\right)$ equal to $1$, and we exhausted rotational invariance by fixing its orientation along the $x$ axis. So setting the norm of $\mathbf{r_2}\left( t = 0\right)$ to $1$ was a mistake in that first ansatz.
With the more general ansatz we are forced to use an additional equation, so that we can solve for both $\phi$ and $l$. This can come, for instance, from $\mathbf{f_2}\left( t = 0\right)$ having no component orthogonal to $\mathbf{r_2}\left( t = 0\right)$. These equations are not particularly nice, but Mathematica is able to provide a solution:
\begin{eqnarray}
\phi &\approx& 1.802119183 \, \, \mathrm{rad} \approx 103.2538234°
\\
&\mathrm{and}&
\\
l &\approx& 1.381620315 \, .
\end{eqnarray}
Notice that $\phi$ is not an angle of the triangle, but the solution allows you to find the largest triangle angle which happens to be
\begin{equation}
\arccos \left[ \frac{(\mathbf{r_2} - \mathbf{r_1}) \cdot (\mathbf{r_3} - \mathbf{r_1})}{|\mathbf{r_2} - \mathbf{r_1}| |\mathbf{r_3} - \mathbf{r_1}|} \right] \approx 1.470038873 \, \, \mathrm{rad} \approx 84.22702317° \, .
\end{equation}