The four polynomials $p_i(x,y,z)$ should be a combination of the following terms:
\begin{align}
\left\{xyz,yx^2,y^2x,zx^2,z^2x,zy^2,z^2y,x^3,y^3,z^3\right\}
\end{align}
other terms would create the different left hand side. I have tried to use brutal force with the help of the deterministic global optimizer BARON starting with 40 unknowns coefficients by minimizing the objective of the form
\begin{align}
((a_{003}^2+b_{003}^2+d_{003}^2+c_{003}^2)-1)^2+
((c_{030}^2+b_{030}^2+d_{030}^2+a_{030}^2)-1)^2+\dots \nonumber
\end{align}
(where $a_{ijk},b_{ijk},c_{ijk},d_{ijk}$ are coefficients to polynomials $p_1,p_2,p_3,p_4$). The calculation recovered the following form of the polynomials $p_1,p_2,p_3,p_4$:
\begin{align}
p_1(x,y,z)&=a_1(yz(y+z)-y^3-z^3) + a_2x^2(y+z)-a_3x(y^2+z^2) \\
p_2(x,y,z)&=b_1yz(z-y)+b_2x^2(y-z)+b_3x(z^2-y^2) \\
p_3(x,y,z)&=c_1yz(y-z)+c_2x^2(y-z)+c_3x(z^2-y^2)+c_4(y^3-z^3) \\
p_4(x,y,z)&=d_1xyz+d_2x(y^2+z^2)-x^3+d_3(y^3+z^3-yz(y+z))
\end{align}
where
\begin{align}
a_1&=0.44778749472139178777752 \\
a_2&=1.6556891275393716966846 \\
a_3&=1.2079562282512847914973 \\
b_1&=-1.9014217828335169269138 \\
b_2&=-0.55113503404381225525555 \\
b_3&=-0.5511281347996555002311 \\
c_1&=1.1890222040768996247806 \\
c_2&=0.75843958935644073537929 \\
c_3&=1.3269533603627028384153 \\
c_4&=0.56886478882077529117822 \\
d_1&=1.8623260525794982367387 \\
d_2&=0.31014344808027582978127 \\
d_3&=0.68984000724248439873065
\end{align}
with error of the objective function $0.43$E-$014$. Hope it helps further. However, analytical expressions of these coefficients should be recovered to check if this solution is really correct.
