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.