Given $x,y,z \in \Bbb R^+$ satisfying
$$
x+y+z = xyz \tag{1},
$$
find the minimum value of the sum $$ S=ax+by+cz, $$ with $a,b,c \in \Bbb R^+$ are known parameters. For example, $$ \color{red}{(a,b,c) = \left(\sqrt 7, \frac{\sqrt{7}}{2},\frac{2\sqrt{7}}{7}\right)} $$
The problem looks very simple. It suffices to find a plane $ax+by+cz = S$ which is tangent to the surface $x+y+z = xyz$.
I know there exists a very clever (but unluckily not intuitive):
From $(1)$, there exists a triangle $ABC$ satisfying: $x = \tan A,y =\tan B$ and $z = \tan C$.
The sum $S$ becomes $S = a \tan A+b \tan B + c\tan C $
Applying the inequality: $\tan x \ge \tan y +(\tan y)'(x-y)$ for $(x,y) = (A,p), (B,q)$ and $(C,r)$. These values $p,q,r$ are determined such that $a(\tan p)' = b(\tan q)' = c(\tan r)'$ and $p+q+r = \pi$
The minimum of $S$ is equal to $a \tan p +b\tan q +c\tan r+...$
For $(a,b,c) = \left(\sqrt 7, \frac{\sqrt{7}}{2},\frac{2\sqrt{7}}{7}\right)$, $S$ reaches its minimum of $\color{red}{\frac{15}{2}}$ when $$ \color{red}{(x,y,z) =\left(\frac{3}{\sqrt{7}}, \frac{5}{\sqrt{7}},\sqrt{7}\right)} $$
I believe there exists a more intuitive solution without using trigonometry.
Here is my attempt (not beautiful yet)
Let denote $$P = ax+by+cz + \lambda(x+y+z - xyz)$$ From the first derivatives of $P$, we deduce $$ \begin{cases} \frac{a}{\lambda}+1 = yz \\ \frac{b}{\lambda}+1 = zx \\ \frac{c}{\lambda}+1 = xy\\ \tag{2} \end{cases} $$ From $(1),(2)$, we have a cubic equation of $\lambda$ $$\frac{1}{yz}+\frac{1}{zx}+\frac{1}{xy} =1 \implies \lambda^3 +\frac{a+b+c}{2}\lambda^2 -\frac{abc}{2} = 0 \tag{3}$$
Remark: At this step, I realize that the solution depends on the parameters $a,b,c$ (and then the equation $(3)$ can have 1 or $3$ roots). For $ \color{red}{(a,b,c) = \left(\sqrt 7, \frac{\sqrt{7}}{2},\frac{2\sqrt{7}}{7}\right)} $, we have $\lambda = -\sqrt{7}/4$, $\frac{1}{4}(-\sqrt{7}-\sqrt{14})$ or $\frac{1}{4}(-\sqrt{7}+\sqrt{14})$.
The second derivarives of $P$ is a Hessian matrix: $$\begin{pmatrix} 0 & \lambda z&\lambda y\\ \lambda z &0 & \lambda x\\ \lambda y & \lambda x &0\\ \end{pmatrix} \tag{4}$$ The equation of 3 eigenvalues of $(4)$ $$p^3 - (x^2+y^2+z^2)p-2xyz =0 \tag{5}$$
Remark: It's very strange because the equation $(5)$ can't have $3$ positive roots nor $3$ negative roots, this means that the Hessian matrix can't be definite positive (negative). And all critical points of $(3)$ aren't extremum points of $S = ax+by+cz$.
Where is my error?
Here is a solution without using trigonometry.
Note: when I wrote the conclusion of this method, I found a second method which is simpler. I don't have time to write both two methods here.
Let denote $$P = ax+by+cz + \lambda(x+y+z - xyz)$$ From the first derivatives of $P$, we deduce $$ \begin{cases} \frac{a}{\lambda}+1 = yz \\ \frac{b}{\lambda}+1 = zx \\ \frac{c}{\lambda}+1 = xy\\ \tag{2} \end{cases} $$ From $(1),(2)$, we have a cubic equation of $\lambda$ $$\frac{1}{yz}+\frac{1}{zx}+\frac{1}{xy} =1 \implies \lambda^3 +\frac{a+b+c}{2}\lambda^2 -\frac{abc}{2} = 0 \tag{3}$$
The equation $(3)$ has one positive root. We will prove that this root is unique in $\Bbb R$. Indeed, we know that $(3)$ can't have $3$ positive roots, then we suppose that there exists a negative root $\lambda'$, then $$yz = 1-\frac{a}{\lambda'}<1 \implies x+y+z = xyz <x \implies y+z<0 \text{ : contradiction}$$
So, there is only one critical point $(x,y,z,\lambda)$ satisfying $(2)$ and $(3)$. ( $(x,y,z,\lambda)$ can be solved analytically).
The second derivarives of $P$ is a bordered Hessian matrix: $$\begin{pmatrix} 0 & -\frac{a}{\lambda} &-\frac{b}{\lambda} &-\frac{c}{\lambda} \\ -\frac{a}{\lambda}& 0 & -\lambda z & -\lambda y \\ -\frac{b}{\lambda} & -\lambda z &0 & -\lambda x\\ -\frac{c}{\lambda} & -\lambda y & -\lambda x &0 \\ \end{pmatrix} \tag{4}$$
Remark: In my attempt above, I made a mistake by not including the $4$-th dimension with the variable $\lambda$. And the bordered Hessian matrix is not necessary to be definite positive in the case of minimum with constraints.
The point $(x,y,z,\lambda)$ can be proved as a minimum point by checking the sufficient condition for local extrema or here . The idea is to calculate the determinant of $2$ minors of the Hessian matrix. The $2$ determinants must be negative (same sign as $(-1)^1 = -1$).
Conclusion: $P$ reaches its minimum equal to $$P_{\text{min}} = 2\frac{x_0y_0z_0}{\lambda}$$
with $(x_0,y_0,z_0,\lambda)$ is the solution of $(2)$ and $(3)$( $(x_0,y_0,z_0,\lambda)$ can be calculated analytically).
Remark: When I was writing the conclusion, I just found a second method which seems simpler than this method. Thanks to the second method, I could have easily the minimum value $P_{\text{min}} = 2\frac{x_0y_0z_0}{\lambda}$ (with the first method above, it's not easy to get $P_{\text{min}}$). I'll write the second method later.