[Edit: I found this method a rather pleasing application of analytic geometry, so I rewrote the explanation hopefully to do it justice. In fact, it can be applied in any dimension, and I've updated the code to be general]
Here's my way:
faces[simplex_] := Partition[simplex, Length@simplex - 1, 1, 1]; (* outward-oriented unit normals to each of faces[a,b,c,d] *) normals[simplex_] := With[{df = Transpose[Transpose @ Rest[#] - First[#]]}, Normalize[(-#.Last@df) # &[Cross @@ (df[[;; Length @ simplex - 2]])]]] & /@ Partition[simplex, Length @ simplex, 1, 1]; insphere[simplex_?(Subtract @@ Dimensions[#] == 1 &)] := Module[{x, center, radius, tangentpoints}, With[{normals = normals[simplex]}, center = Array[x, Length @ simplex - 1]; tangentpoints = center + radius # & /@ normals; insphere[center, radius, tangentpoints] /. First @ Solve[ MapThread[ Function[{f, n}, (center + radius n - First[f]) . n == 0], {faces[simplex], normals}], Flatten[{center, radius}]] ]] (* some interfaces to the insphere data *) insphere[center_, radius_, tangentpoints_]["Sphere"] /; Length[center] == 2 := Circle[center, radius]; insphere[center_, radius_, tangentpoints_]["Sphere"] /; Length[center] == 3 := Sphere[center, radius]; insphere[center_, radius_, tangentpoints_]["Sphere"] := {center, radius}; insphere[center_, radius_, tangentpoints_]["center"] := center; insphere[center_, radius_, tangentpoints_]["radius"] := radius; insphere[center_, radius_, tangentpoints_]["tangentpoints"] := tangentpoints;
The idea is that the condition that defines the insphere is that the perpendiculars dropped from the center to the faces are all equal. This leads to a system of linear equations that is easy for Solve to deal with.
For a given face $f$, let $\bf n$ be its unit normal vector pointing out of the tetrahedron. Let $(x, y, z)$ be the unknown center of the insphere and $r$ be its radius. Then we seek a condition that $(x, y, z) + r\,{\bf n}$ lies on the face $f$. Let ${\bf v}_1=(a_1, b_1, c_1)$ and ${\bf v}_2=(a_2, b_2, c_2)$ be the displacement vectors of two edges of $f$, and let ${\bf P}_0$ be one of the vertices of $f$. Then the unit outward normal $\bf n$ will be the normalization of $\pm({\bf v}_1 \times {\bf v}_2) $. The standard point-normal form for the equation of a plane containing a point ${\bf P}_0$ and perpendicular to a vector $\bf n$ is, for an arbitrary point ${\bf P}$ on the plane, $$({\bf P}-{\bf P}_0) \cdot {\bf n} = 0\,.$$ With our variables, in which ${\bf P}$ equals center + radius n and ${\bf P}_0$ equals First[f], this is equivalent to
(center + radius n - First[f]) . n == 0
To get the orientation of the normal $\bf n$ right, we take the three vertices of a face plus the fourth vertex of the tetrahedron. Subtract the first face vertex from the other three. The first two vectors represent the edges of the face coming out of the first vertex. The last vector points toward the interior side of the face, toward the fourth vertex, and its dot product with the outward normal will be negative. Scaling the cross product of the edge vectors by the negative of dot product of the cross product and the last vector, -#.Last@df, ensures the normal points outward. Normalizing the result yields our unit outward normal $\bf n$.
The function insphere[tet] returns a database in the form
insphere[center, radius, tangentpoints]
which may be accessed through functions like those at the end of the code above. Alternatively, to get a return value like the one specified by the OP, replace it in the code for insphere by
{center, radius, tangentpoints} /. First @ Solve[...]
Example
a = {0, 0, 0}; b = {1, 0, 0}; c = {2, 1, 0}; d = {1, 3/2, 1}; midpoints[tet_] := Mean /@ faces[tet]; tet0 = {a, b, c, d}; in0 = insphere[tet0]; Graphics3D[{ Opacity[0.4], Polygon[faces[tet0]], MapThread[Arrow[{##}] &, Accumulate @ {midpoints[tet0], normals[tet0]}], Red, in0["Sphere"], Black, Point @ in0["tangentpoints"], White, Opacity[1], Polygon[faces @ in0["tangentpoints"]] }]
Additional examples
2D
SeedRandom[1]; With[{in0 = insphere[#]}, Graphics[{Point @ in0["tangentpoints"], Opacity[0.4], Line @ faces @ #, in0["Sphere"]}]] &@ RandomReal[1, {3, 2}]
4D
SeedRandom[3]; proj0 = N @ Orthogonalize[IdentityMatrix[4] ~Prepend~ {1, 1, 1, -2}][[2 ;; 4]]; proj[pts_?VectorQ] := proj0.pts; proj[pts_?MatrixQ] := Transpose[proj0.Transpose @ pts]; sim = RandomReal[1, {5, 4}]; in0 = insphere[#] &@ sim; in0["Sphere"] Graphics3D[{Point @ proj @ in0["tangentpoints"], Opacity[0.5], Sphere[proj @ in0["center"], in0["radius"]], Opacity[0.25], MapIndexed[{Hue[First @ #2 / 4], Polygon[faces @ proj @ #1]} &, faces @ sim] }]
{{0.335297, 0.4023, 0.634872, 0.392616}, 0.059207}
Graphics3D@{Opacity[0.5],Polygon[{a,b,c,d}~Subsets~{3}]}$\endgroup$