Solve four equations system considering that the picked solution should minimize a function

I wrote equations that represents the rotation on a gimbal system for pointing a telescope with 3 DOF,
I have a system in four equations and 4 unknowns, but as you can figure out easly the system has infinite solutions ( you can point a telescope that has only 2 DOF, so with another axis you can combine three rotations in several ways ) , so I’d like to force the program to choose the solution that provides the minimum value of the Function MinimalFunction
aa bb cc and teta0 rho0 phi0 are known

<< Quaternions`
Qz = Quaternion[Cos[teta/2], 0*Sin[teta/2], 0*Sin[teta/2], 1*Sin[teta/2]];
Qx = Quaternion[Cos[rho/2], Cos[gamma]*Sin[rho/2], 0*Sin[rho/2], Sin[gamma]*Sin[rho/2]];
Qy = Quaternion[Cos[phi/2], 0*Sin[phi/2], 1*Sin[phi/2], 0*Sin[phi/2]];
QQ = Qz ** Qx ** Qy ;
quaternionx = Quaternion[0, K, 0, 0];

Finale = FullSimplify[QQ** quaternionx ** Conjugate[QQ]];
AsseXTele = {Finale[[2]], Finale[[3]], Finale[[4]]}

Espressione1 =  AsseXTele [[1]] - aa == 0;
Espressione2 =  AsseXTele [[2]] - bb == 0;
Espressione3 =  AsseXTele [[3]] - cc == 0;
MinimalFunction=  KTeta*(teta - teta0)^2 + KRho*(rho - rho0)^2 + KPhi*(phi - phi0)^2
Solve[Espressione1 && Espressione2 && Espressione3 && aa^2 + bb^2 + cc^2 == K, { teta, rho, phi,K}]

<< Quaternions`
Qz = Quaternion[Cos[teta/2], 0*Sin[teta/2], 0*Sin[teta/2], 1*Sin[teta/2]];
Qx = Quaternion[Cos[rho/2], Cos[gamma]*Sin[rho/2], 0*Sin[rho/2],Sin[gamma]*Sin[rho/2]];
Qy = Quaternion[Cos[phi/2], 0*Sin[phi/2], 1*Sin[phi/2], 0*Sin[phi/2]];
QQ = Qz ** Qx ** Qy;
quaternionx = Quaternion[0, k, 0, 0];

Finale = FullSimplify[QQ ** quaternionx ** Conjugate[QQ]];
AsseXTele = {Finale[[2]], Finale[[3]], Finale[[4]]}
KTeta = 1; teta0 = 0; KRho = 1; rho0 = 1; KPhi = 1; phi0 = 0;
Espressione1 = AsseXTele[[1]] - aa;
Espressione2 = AsseXTele[[2]] - bb;
Espressione3 = AsseXTele[[3]] - cc;
MinimalFunction = KTeta*(teta - teta0)^2 + KRho*(rho - rho0)^2 + KPhi*(phi - phi0)^2;
restrs = {Espressione1 == 0, Espressione2 == 0, Espressione2 == 0, aa^2 + bb^2 + cc^2 - k == 0}
sol = NMinimize[Join[{MinimalFunction}, restrs], {gamma, teta, rho, phi, k, aa, bb, cc}]
{gamma, teta, rho, phi, k, aa, bb, cc}/.sol[[2]]