In[]:=
cons1=allcons/.{r->s,x->x[1],y->y[1],z->z[1]};
In[]:=
cons2=allcons/.{r->s/2,x->x[2],y->y[2],z->z[2]};
In[]:=
dist[pt1_,pt2_]:=(pt1-pt2).(pt1-pt2)
In[]:=
cntrVars[i_]:={x[i],y[i],z[i]}
In[]:=
consDJ[i_,j_,s_,r1_,r2_]:=
2
(s*r1+s*r2)
<=dist[cntrVars[i],cntrVars[j]]
In[]:=
consdj=consDJ1,2,s,1,
1
2

Out[]=
9
2
s
4
≤
2
(x[1]-x[2])
+
2
(y[1]-y[2])
+
2
(z[1]-z[2])
In[]:=
allCons=Join[cons1,cons2,{consdj}];
In[]:=
sln=NMaximize[{s,allCons},{s,x[1],y[1],z[1],x[2],y[2],z[2]}]
Out[]=
{0.204124,{s0.204124,x[1]2.72643×
-14
10
,y[1]7.82214×
-15
10
,z[1]-1.8116×
-8
10
,x[2]-0.144338,y[2]-0.25,z[2]-0.102062}}
In[]:=
Show[​​Graphics3D@{Opacity[0.3],Tetrahedron[]},​​Graphics3D@Ball[{x[1],y[1],z[1]},s]/.sln[[2]],​​Graphics3D@Ball[{x[2],y[2],z[2]},s/2]/.sln[[2]]​​,ImageSize->Small,Axes->True]
Out[]=