I now know how useful analytic geometry can be in bashing geometry problems involving side lengths. Does anybody have any tips on how to fix coordinates to keep the solution from becoming too tedious?
The solution doesn't have to be so tedious. With the right choice of coordinates it becomes simple.
When dealing
Choose the origin as a point where a lot of lines intersect, or as a center of a circle, or as a center of mass. Choose the axis so that a lot of lines are parallel to them.
For the specific problem, probably best to choose $P$ as the origin. You might choose $G$, but it has a complex relationship between $A,B,C$, so you'd have to put conditions on those coordinates.
If you wait to substitute for $G=(g_1,g_2)$, you get:
$$\begin{align}GA^2+GB^2+GC^2+3GP^2 = &(g_1-x_1)^2+(g_2-x_2)^2+\\ &(g_1-y_1)^2+(g_2-y_2)^2+\\ &(g_1-z_1)^2+(g_2-z_2)^2+\\ &3g_1^2 + 3g_2^2\\ =&6g_1^2+(x_1^2+y_1^2+z_1^2)-2g_1(x_1+y_1+z_1) +\\ &6g_2^2+(x_2^2+y_2^2+z_2^2)-2g_2(x_2+y_2+z_2) \end{align}$$
Only here at the end do you use that $x_1+y_1+z_1=3g_1$ and $x_2+y_2+z_2=3g_2$. Then you see that $6g_i^2-2g_i(x_i+y_i+g_i)=0$.
This actually reveals something interesting, by the way - any equality involving the linear combination of squares of lengths and only "linear combinations" of a fixed set of otherwise independent points, is true for all points if any only if it is true on each coordinate - that is, if you assume the points are all on a line. It generalizes to a few other types of expression, as well - mostly, via the law of cosine.
These particular problems work without coordinates, using vector notation. If you define, for $x=(x_1,x_2)$ and $y=(y_1,y_2)$ that $x\cdot y = x_1y_1+x_2y_2$, then the distance between $x$ and $y$, squared, is $(x-y)\cdot(x-y)=x\cdot x -2x\cdot y +y\cdot y$. This "dot product" is nice so you get good properties, like $(a+b)\cdot c = a\cdot c + b\cdot c$, and $a\cdot b = b\cdot a$.
Then if $P=(0,0)$ and $a,b,c$ are any points, and $g=\frac{a+b+c}{3}$, then your first question is:
$$\begin{align} |AG|^2+|BG|^2+|CG^2|+3|PG|^2=\,& (g-a)\cdot (g-a) + \\ &(g-b)\cdot (g-b) +\\ &(g-c)\cdot (g-c) + \\ &3g\cdot g\\ =\,&6g\cdot g -2g\cdot (a+b+c) + a\cdot a + b\cdot b + c\cdot c\\ =\,&6g\cdot g - 6g\cdot \frac{a+b+c}{3} + a\cdot a + b\cdot b + c\cdot c\\ =\,&6g\cdot g - 6g\cdot g +a\cdot a + b\cdot b + c\cdot c\\ =\,&|PA|^2+|PB|^2+|PC|^2 \end{align}$$
By the way, this equality is actually deeply related to statistics, and can be generalized to any number of points in any number of dimensions:
$$|A_1G|^2+\cdots+|A_nG|^2 + n|PG|^2 =|A_1P|^2+\cdots +|A_nP|^2$$
where $G$ is the center of gravity of the $n$ points $A_1,A_2,\dots,A_n$.
This indicates that the mean of a bunch of values "minimizes" a certain function. If the right side is $f(P)$, then we see it as:
$$f(P)=|GP|^2+f(G)$$ and $P=G$ is the point where $f$ is minimized. The value of $\frac{1}{n}f(G)$ is called the "variance" of $A_1,\dots,A_n$.