Given f(x,y,z) = x^2+y^2-z^2-1, and p = (x,y,z), the distance between point p_0=(1,3,1) and p in f(x,y,z)=0 is given by d = norm(p-p_0). Here instead of d we will consider d^2=(p-p_0)*(p-p_0) observing that the minimum of d is also the minimum of d^2.
The problem restated is:
Find the point p^+ = (x^+,y^+,z^+) such that
min d^2(x,y,z)=d^2(x^+,y^+,z^+) and obeying f(x^+,y^+,z^+)=0
We know that d^2=(x-x_0)^2+(y-y_0)^2+(z-z_0)^2.
First we formulate the Lagrangian which is:
L(x,y,z,lambda)=f(x,y,z)+lambda d^2(x,y,z)
L(x,y,z,lambda) = x^2+y^2-z^2-1+lambda((x-x_0)^2+(y-y_0)^2+(z-z_0)^2)
Second we calculate the partial derivatives of L regarding (x,y,z,lambda) so
L_x=2 (-1 + x) + 2 x lambda
L_y=2 (-3 + y) + 2 y lambda
L_z=2 (-1 + z) - 2 z lambda
L_{lambda} =-1 + x^2 + y^2 - z^2
Third we equate L_x=0,L_y=0,L_z=0,L_{lambda}=0
solving for (x^+,y^+,z^+) finding
{(x_1^+= -0.322242), (y_1^+ = -0.966725), (z_1^+ = 0.195953), (lambda_1^+= -4.10326):}
and
{(x_2^+= 0.678696), (y_2^+ = 2.03609), (z_2^+ = 1.89902), (lambda_2^+= 0.473413):}
Fourth we qualify the found points. Calculating the Hessian matrix
of d^2(x,y,z) obtaining
H(x,y,z) = ((f_{x x},f_{xy},f_{xz}),(f_{yx},f_{yy},f_{yz}),(f_{zx},f_{zy},f_{zz})) = ((2,0,0),(0,2,0),(0,0,2)),
with positive eigenvalues, qualifying the point as a local minimum.
Fifth we calculate d_1=sqrt(d^2(x_1^+,y_1^+,z_1^+)) = 4.2579
and d_2=sqrt(d^2(x_2^+,y_2^+,z_2^+)) = 1.35669
as local minimal distances.