This answer is trying to answer the question with a full solution, but moreover the questions behind the question, which are:
- $\bbox[yellow]{(0)}$ General discussion, which is the right framework?
- $\bbox[yellow]{(1)}$ Which is the minimal value of the function $f:[0,\infty)^3\to\Bbb R$,
$$
f(a,b,c)=\sum_{\text{cyclic}}\sqrt{5a+5b+8ab}\ ,
$$
constrained to the relation / on the subset of all $(a,b,c)$ with
$$
1 = ab+bc+ca\ ?
$$
- $\bbox[yellow]{(2)}$ Which is the complexity involved in the computation?
- $\bbox[yellow]{(3)}$ Which methods are available to investigate humanly this minimal value?
$\bbox[yellow]{(0)}$
In a first part i will use CAS (computer algebra software, sage) in order to have a first firm answer with intermediate check points that may serve as a further orientation. The complexity will be transparent from the computer data. It remains to draw conclusions for the human point of view, and go into the one or other way to humanly attack the problem.
First of all, for the given kind of problem we have already a method, it uses Lagrange multiplicators, and universities do want its applications (unfortunately only in cheaper situations), and sites like MSE insist to avoid it (arguably because it is not elementary enough, pretending that an olimpiad-like solution is only wanted. Yes, if there is such a solution it is welcome. But notice that very often a problem like ours here will never survive into the final list of an NMO / IMO, say. Except for special counties, i've got such problems in the live stream. And there are reasons for this elimination.)
In my eyes, there is no reason for rejecting such an efficient analysis tool in problems of extremum in situations like ours here.
So let us combine the human and the computer force to have a first verdict.
$\bbox[yellow]{(1)}$
We associate the function
$F:(0,\infty)^3\times\Bbb R$ defined by
$$
F(a,b,c,m)=f(a,b,c) - m(ab+bc+ca-1)\ .
$$
Its $\inf$ for points on the domain
$$
D=\{\ (a,b,c)\ :\ a,b,c>0\ ,\ ab+bc+ca=1\ \}
$$
are
- either at the boundary of $D$, i.e. points with $abc=0$, $ab+bc+ca=1$,
- or points of $(0,\infty)^3\times \Bbb R$ where $F'=(F'_a,F'_b,F'_c,F'_m)$ vanishes. (The condition $F'_m=0$ makes such a point belong to $D$.)
$\bbox[lightyellow]{(1a)}$
Let us examine the boundary points first. Using the symmetry, we may and do assume $c=0$. The constraint on $(a,b,c)$ implies $ab=1$. Then $f(a,b,c)/\sqrt 5$ is
$$
\sqrt{a+b+sab} + \sqrt{a}+\sqrt{b}\ ,
$$
with $\bbox[lightgreen]{\ s=8/5\ }$ for short, and we minimize this function.
Apply again Lagrange multipliers.Let $G$ be the function
$$(a,b,m)\to F(a,b,0,m)/\sqrt 5
=G(a,b,m)
=\sqrt{a+b+s} + \sqrt {a}+\sqrt {b}-m(ab-1)\ ,
$$
and work with the condition $ab=1$. Boundary (with $ab=1$) examination. For $a\to 0$ or $b\to 0$ we have $b\to\infty$ respectively $a\to\infty$, so $F(a,b,0,m)\to\infty$, no infimum at the boundary. For the local extremal points we ask for $G'_a=G'_b=G'_m=0$.
The condition $G'_m=0$ is $ab=1$, i.e. $a=1/b$. We obtain:
$$
\begin{aligned}
\frac 1{2\sqrt{a+b+s}} +\frac 1{2\sqrt a}-mb
&=0
=\frac 1{2\sqrt{a+b+s}} +\frac 1{2\sqrt b}-ma
\ ,\text{ i.e.}
\\
\frac a{\sqrt{a+b+s}} + \sqrt a
&=2m
=\frac b{\sqrt{a+b+s}} +\sqrt b\ .
\end{aligned}
$$
We ignore $m$, use the equation resulted from above by subtracting the right part from the left side, and can factor $\sqrt a-\sqrt b$. The remained factor is $\frac {\sqrt a+\sqrt b}{\sqrt{a+b+s}} + 1>0$, so there is only one point that comes in discussion. It must correspond to the one absolute minimum (which exists, since the boundary discussion shows that $G\to \infty$ at boundary). So $f(a,b,0)$ is minimal for $a=b=1$, we obtain the value $3\sqrt 2+2\sqrt 5$ for this point. This is a proof for the minimum value on the boundary so far.
$\bbox[lightyellow]{(1b)}$
We will compare this value with local extremal (minimal) values.
When computing derivatives it is useful to use $f(a,b,c)/\sqrt 8$,
which is the cyclic sum of the expression $\sqrt{ab+s(a+b)}$, where this
time $\bbox[lightgreen]{\ s=5/8\ }$ for short.
Then $f$ and $f/\sqrt 8$ are minimal at the same points. We use the latter version, for typing reasons. Instead of $F$ we have its version
We set and compute:
$$
\begin{aligned}
F(a,b,c;m)&:=
\underbrace{\sqrt{bc+s(b+c)}}_{\text{Notation: }x=x(a,b,c)} +
\underbrace{\sqrt{ca+s(c+a)}}_{\text{Notation: }y=y(a,b,c)} +
\underbrace{\sqrt{ab+s(a+b)}}_{\text{Notation: }z=z(a,b,c)}
-m(bc+ca+ab-1)\ ,\\
&\qquad\text{ Then the system $0=F'_a=F'_b=F'_c=F'_m$ has the equations:}
\\
0&=\frac 1{2x}\cdot 0 +\frac 1{2y}(c+s) +\frac 1{2z}(b+s) - m(b+c)\ ,\\
0&=\frac 1{2x}(c+s) +\frac 1{2y}\cdot 0 +\frac 1{2z}(a+s) - m(c+a)\ ,\\
0&=\frac 1{2x}(b+s) +\frac 1{2y}(a+s) +\frac 1{2z}\cdot 0 - m(a+b)\ ,\\
1 &=bc+ca+ab\ .
\end{aligned}
$$
The first three equations may be seen as a linear system in $1/(2x)$,
$1/(2y)$, $1/(2z)$,
$$
\begin{bmatrix}
0 & c+s & b+s\\
c+s & 0 & a+s\\
b+s & a+s & 0
\end{bmatrix}
\begin{bmatrix}
1/(2x)\\
1/(2y)\\
1/(2z)
\end{bmatrix}
=
m
\begin{bmatrix}
b+c\\
c+a\\
a+b
\end{bmatrix}
\ ,
$$
the involved $3\times 3$-matrix has determinant $2(a+s)(b+s)(c+s)$,
and its solution isolates the corresponding values with this common denominator, so we obtain four equations in $a,b,c;m$ of the shape:
$$
\begin{aligned}
\frac1{2\sqrt{bc+s(b+c)}}
=\frac 1{2x}
&=\frac{m\color{gray}{(a+s)}(bc+as)}{\color{gray}{(a+s)}(b+s)(c+s)}\ ,
\\
\frac1{2\sqrt{ca+s(c+a)}}
=\frac 1{2y}
&=\frac{m\color{gray}{(b+s)}(ca+bs)}{(a+s)\color{gray}{(b+s)}(c+s)}\ ,
\\
\frac1{2\sqrt{ab+s(a+b)}}
=\frac 1{2z}
&=\frac{m\color{gray}{(c+s)}(ab+cs)}{(a+s)(b+s)\color{gray}{(c+s)}}\ ,\\
1 &= bc+ca+ab\ .
\end{aligned}
$$
We can eliminate $m$ now, there are three expressions equal to the
square of $(a+s)(b+s)(c+s)/(2m)$, well, squaring may introduce solutions (corresponding to $\pm m$, but humanly we do so), and obtain:
$$
(bc+s(b+c))(a+s)^2(bc+as)^2 =
(ca+s(c+a))(b+s)^2(ca+bs)^2 =
(ab+s(a+b))(c+s)^2(ab+cs)^2 \ ,
$$
this leads to two equations of degree six in $a,b,c$.
Subtracting any two of the three parts of the double equality
introduces one of the factors $(a-b)$, $(b-c)$, $(c-a)$,
but the remained factor has huge degree for the human hand.
It is time to face the reality, let us ask the machine which are the
real positive solutions of the system. This is a giant leap for the man, but a small one for the mankind, so let us keep it small:
R.<a,b,c,X,Y,Z,m> = PolynomialRing(QQ)
# we are using X, Y, Z instead of 1/x, 1/y, 1/z
S.<t> = PolynomialRing(QQ)
s = 5/8
J = R.ideal([
X/2*0 + Y/2*(c + s) + Z/2*(b + s) - m*(b + c) ,
X/2*(c + s) + Y/2*0 + Z/2*(a + s) - m*(c + a) ,
X/2*(b + s) + Y/2*(a + s) + Z/2*0 - m*(a + b) ,
-1 + X^2*(b*c + s*(b + c)) ,
-1 + Y^2*(c*a + s*(c + a)) ,
-1 + Z^2*(a*b + s*(a + b)) ,
-1 + b*c + c*a + a*b , ])
pa = J.elimination_ideal([X,Y,Z, b,c, m]).gens()[0]
print(f'After elimination of all other variables, a satisfies the equation:')
print(f'0 = {pa.subs(a=t).factor()}')
And the prints deliver the complexity of the problem:
After elimination of all other variables, a satisfies the equation:
0 = (3221225472000)
* (t^2 - 1/3)
* (t^3 - 21/8*t^2 + t + 45/32)
* (t^3 + 11/8*t^2 + t - 5/32)
* (t^5 + 7/80*t^4 + 501/512*t^3 + 421/1280*t^2 - 4021/25600*t + 25481/2048000)
* (t^5 + 5/8*t^4 - 75/128*t^3 - 25/128*t - 125/512)
And each of the involved non-constant factor has a positive root in $(0,\infty)$,
we have expected the root $a=1/\sqrt 3$, from $t^2-1/3=0$, but also other factors come with positive roots:
sage: for f, mul in pa.subs(a=t).factor():
....: print(f'{f}\n\tvanishes in:\n\t{[r for r in f.roots(ring=AA, multiplicities=False) if 0 < r < 1]}')
....:
t^2 - 1/3
vanishes in:
[0.5773502691896258?]
t^3 - 21/8*t^2 + t + 45/32
vanishes in:
[]
t^3 + 11/8*t^2 + t - 5/32
vanishes in:
[0.1305786882139770?]
t^5 + 7/80*t^4 + 501/512*t^3 + 421/1280*t^2 - 4021/25600*t + 25481/2048000
vanishes in:
[0.12193648467888863?, 0.1632193275966720?]
t^5 + 5/8*t^4 - 75/128*t^3 - 25/128*t - 125/512
vanishes in:
[0.8500133943883169?]
A "manual" consideration has to overcome with thus difficulties.
$\bbox[yellow]{(2)}$ The involved complexity is thus too much for the pen and paper approach. Even with a clever handling of the above system in $a,b,c;m$, we still have to consider (algebraically) complicated points. One speculation would be that some positivity argument may rule out some of the many possible $a$-values above. Well, let us see which are all the solutions of the given system. We are again using computer support, alternatively, we use the system in $a,b,c$ that we have finally obtained, and restrict to positive solutions.
For each solution we also compute the value of $f$ in it. Code:
s = 5/8
R.<a,b,c> = PolynomialRing(QQ)
def E(a, b, c):
return (b*c + s*(b + c)) * (a + s)^2 * (b*c + a*s)^2
J = R.ideal([
a*b + b*c + c*a - 1
, E(a,b,c) - E(b,c,a), E(a,b,c) - E(c,a,b)])
for point in J.variety(ring=AA):
A, B, C = point.values()
val = ( sqrt(8*A*B + 5*A + 5*B) +
sqrt(8*B*C + 5*B + 5*C) +
sqrt(8*C*A + 5*C + 5*A) )
if A >= B and B >= C and C > 0:
print("The following point (a, b, c) is an ordered solution:")
print(f"\ta = {A} with minimal polynomial = {A.minpoly()}")
print(f"\tb = {B} with minimal polynomial = {B.minpoly()}")
print(f"\tc = {C} with minimal polynomial = {C.minpoly()}")
print("The value of the initial function in this point is:")
print(f"\tf(a, b, c) = {val} with minimal polynomial {val.minpoly()}\n")
This gives:
The following point (a, b, c) is an ordered solution:
a = 0.5773502691896258? with minimal polynomial = x^2 - 1/3
b = 0.5773502691896258? with minimal polynomial = x^2 - 1/3
c = 0.5773502691896258? with minimal polynomial = x^2 - 1/3
The value of the initial function in this point is:
f(a, b, c) = 8.71559087079392? with minimal polynomial x^4 - 48*x^2 - 2124
The following point (a, b, c) is an ordered solution:
a = 0.8500133943883169? with minimal polynomial = x^5 + 5/8*x^4 - 75/128*x^3 - 25/128*x - 125/512
b = 0.8500133943883169? with minimal polynomial = x^5 + 5/8*x^4 - 75/128*x^3 - 25/128*x - 125/512
c = 0.1632193275966720? with minimal polynomial = x^5 + 7/80*x^4 + 501/512*x^3 + 421/1280*x^2 - 4021/25600*x + 25481/2048000
The value of the initial function in this point is:
f(a, b, c) = 8.74927159655177? with minimal polynomial x^10 - 163/4*x^8 - 568387/256*x^6 - 613979/32*x^4 - 6078993/4*x^2 - 4578338
Both values (possible local extremal values) are greater than the boundary minimum $f(0,1,1)=f(1,0,1)=f(1,1,0)$:
sage: QQbar(sqrt(18) + 2*sqrt(5))
8.71477664211887?
So the absolute minimum is the boundary minimum $3\sqrt 2+2\sqrt 5$.
$\bbox[yellow]{(3)}$
Can we do better with human eyes and hands?
One idea would be to use convexity arguments. But unfortunately near the boundary, especially near the minimal point, we have in all directions an increasing function, and "somewhere in the middle" in $(m,m,m)$ with $m=1/\sqrt3$ we also have a local minimum. So there will be portions (on the line / on any path) from $(0,1,1)$ to $(m,m,m)$ of concavity. It is hard to make them explicit.
Maybe a picture of the situation may help to visualize the problem.


This is one and the same spherical plot picture shown from different angles.
Our function $f$ is red, but to make the picture focus on the critical part the plot was done for $\max(10, f)$, this red plot has full opacity.
A spherical plot of a (positive) constant function is a part of the sphere with the radius given by the constant.
For orientation we have in the same pictures such spherical constant plots for the following values:
- ten, the gray truncation of $f$,
8.714776642118
, which corresponds numerically to (almost) $\min f$ taken at the boundary, it is green, and we allow some transparency to see through it,
- and finally
8.716
, the cyan piece of the sphere, which is opaque.
The intransparent red and cyan surfaces show where we expect the locally minimal value. In that portion of the plot, red, cyan, and green are almost indistinguishable as surfaces, we need colors to observe the difference.
The picture was computed as follows.
def E(a, b): return sqrt(8*a*b + 5*a + 5*b)
def f(a, b, c): return E(b, c) + E(c, a) + E(a, b)
def myplot(u, v):
A, B, C = cos(u)*sin(v), sin(u)*sin(v), cos(v)
den = sqrt(B*C + C*A + A*B)
a, b, c = A/den, B/den, C/den # thus bc + ca + ab = 1
return min(10, f(a, b, c))
def ten(u, v): return 10
def mymin(u, v): return 8.714776642118
def mymid(u, v): return 8.716
eps = 0.000000000000001
R = (0 + eps, pi/2 - eps)
g1 = spherical_plot3d(myplot, R, R, color='red', opacity=1)
g2 = spherical_plot3d( ten, R, R, color='lightblue', opacity=0.2)
g3 = spherical_plot3d( mymin, R, R, color='green', opacity=0.5)
g4 = spherical_plot3d( mymid, R, R, color='cyan', opacity=1)
(g1 + g2 + g3 + g4).show(aspect_ratio=1)
The above lines were presenting a full computer aided solution to the problem.
Is there any human way to go? I have to doubt convexity arguments, since there cannot be any convexity - at least near the absolute minimal values at the boundary. Also, going algebraically is also a hard elimination task. Critical points of algebraic degree five (as algebraic numbers, numbers in $\bar{\Bbb Q}$) do occur, and we have to say something about these numbers, yes, the true local minimum is "a" simple algebraic number $1/\sqrt3$ (repeated three times on the components), but we still have to compute it and compare its value with the other critical points of degree five, or at least compute the Hessian / the convexity in these degree five critical points. So there is no "easy" way to go for a human (that fits the reasonable framework), i am for legitimating computer assistance in such cases. (Well, i would have done it first also without legitimation.)