These are 3 circles radius h1-3, of different centers (shifted in z which is immaterial). A solution could be be any intersections, but no equation yet creates coupling/dependency between these three equations of 6 unknowns. Bye the way Mathematica could graph the first 3 circles and even give the points using it's library of new geometry helper functions)
"wrong symmetry of placement", no go
In[384]:= x1 = 1; y1 = 2; z1 = 4;
x2 = 1; y2 = 2; z2 = 4;
x3 = 1; y3 = 2; z3 = 4;
will be True
In[387]:= x1 = 1; y1 = 1; z1 = 1;
x2 = 3; y2 = 3; z2 = 3;
x3 = 7; y3 = 7; z3 = 7;
In[390]:= exp1 = (xd1 - x1)^2 + (yd1 - y1)^2 + z1^2 == h1^2
Out[390]= 1 + (-1 + xd1)^2 + (-1 + yd1)^2 == h1^2
In[391]:= exp2 = (xd2 - x2)^2 + (yd2 - y2)^2 + z2^2 == h2^2
Out[391]= 9 + (-3 + xd2)^2 + (-3 + yd2)^2 == h2^2
In[392]:= exp3 = (xd3 - x3)^2 + (yd3 - y3)^2 + z3^2 == h3^2
Out[392]= 49 + (-7 + xd3)^2 + (-7 + yd3)^2 == h3^2
That these are all true "says something". That intersection of the spheres is not the solution.
But symmetry in placement is a rule, likely there is some nifty geometric interpretation for each
In[393]:= exp4 = ((x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2)^(1/
2) + ((x2 - x3)^2 + (y2 - y3)^2 + (z2 - z3)^2)^(1/
2) == ((x1 - x3)^2 + (y1 - y3)^2 + (z1 - z3)^2)^(1/2)
Out[393]= True
In[394]:= exp5 = ((z2 - z1)/(((x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2)^(1/2))) == ((z3 -
z2)/(((x2 - x3)^2 + (y2 - y3)^2 + (z2 - z3)^2)^(1/2)))
Out[394]= True
In[395]:= exp6 = ((z3 - z1)/(((x1 - x3)^2 + (y1 - y3)^2 + (z1 - z3)^2)^(1/2))) == ((z3 -
z2)/(((x2 - x3)^2 + (y2 - y3)^2 + (z2 - z3)^2)^(1/2)))
Out[395]= True
But it also says something else: that even when the last 3 are True first three "are still not coupled", they do not depend on each other, in any way. (it was my purpose so far to bench check if the equations are coupled and make sense, btw). So really, you could say, you have 3 equations in 12 unknowns, or 3 eq. in 6 unknowns if your statement 6 are variable must hold, ...
I hope that was not too rudimentary and helpful. I do not have a formula to compute "any system of non-linear matrix". I can say you might try reducing the spheres to points (get that solution, then ask about how to parlay that into one of circles and spheres "extra padding" rules). There is another approach which is to convert your spheres into planes by mapping, solve that, then transform them back into spheres or circles, which would be a length discussion.
In[400]:= sol = NSolve[{exp1, exp2, exp3, exp4, exp5, exp6}, {xd1, yd1, xd2, yd2, xd3,
yd3}, Reals];
During evaluation of In[400]:= NSolve::svars: Equations may not give solutions for all "solve" variables.
yd1 -> ConditionalExpression[
1. - 1. Sqrt[-2. + h1^2 + 2. xd1 - 1. xd1^2], Or[
And[d`h1 > 1., d`h2 > 3., d`h3 > 7.,
Inequality[
1. - 1. (-1. + d`h1^2)^Rational[1, 2], Less, d`xd1, Less,
1. + (-1. + d`h1^2)^Rational[1, 2]],
<SNIP (a very long list of conditions), and large number of suggest solutions.>