Prove $\sum_{\text{cyc}}^{}\sqrt[3] {1+2ac} \le 3\sqrt [3] {3}$.

TLDR

A standard computer-assisted (but rigorous) proof is to use Lagrange multipliers method together with interval arithmetic libraries such as IntervalRoots.jl.

---

We are optimizing within a compact set in $\mathbb R^3$ as shown below

So there exist maximum points, either in the interior, or on the boundaries.

We can use Lagrange method for the interior. Let $$ f(a, b, c)=\sqrt[3]{2 a b+1}+\sqrt[3]{2 a c+1}+\sqrt[3]{2 b c+1}-3 \sqrt[3]{3} $$ and $$ g(a,b,c,l) = f(a,b,c)+l (a b c+a+b+c-4). $$ Then we just to find the critical points of $g$, i.e., solve $\nabla g = 0$, i.e., \begin{align} \frac{2 b}{3 (2 a b+1)^{2/3}}+\frac{2 c}{3 (2 a c+1)^{2/3}}+l (b c+1)&=0, \\ \frac{2 a}{3 (2 a b+1)^{2/3}}+l (a c+1)+\frac{2 c}{3 (2 b c+1)^{2/3}}&=0, \\ l (a b+1)+\frac{2 a}{3 (2 a c+1)^{2/3}}+\frac{2 b}{3 (2 b c+1)^{2/3}}&=0, \\ a b c+a+b+c-4&=0. \end{align} Look at this for a while and you will see that one solution is $$ a=b=c=1, l = -\frac{2}{3\ 3^{2/3}} $$ and this should give us the maximum of $f(1,1,1)=0$. To rule out other solutions, we can use rigorous numerics libararies like IntervalRoots.jl.

It's not difficult to see that the solution $(a,b,c,l)$ can only be within $[0,4]^3 \times [-55,0]$. The following Julia code finds all such solutions rigorously.

using IntervalArithmetic, IntervalRootFinding, ForwardDiff

f((a, b, c)) = -3*3^(1/3) + (1 + 2*a*b)^(1/3) + (1 + 2*a*c)^(1/3) + (1 + 2*b*c)^(1/3)
g((a, b, c, l))=f((a, b, c))+l*(a + b + c + a*b*c - 4)
∇g = ∇(g)
box = IntervalBox(0..4,3)×(-55..0)
rts = roots(∇g, box, Krawczyk, 1e-5)
println(rts)

And the result is

Root{IntervalBox{4,Float64}}[Root([0.999999, 1.00001] × [0.999999, 1.00001] × [0.999999, 1.00001] × [-0.3205, -0.320499], :unique)]

To see why only chekcing $l \in [-55,0]$ suffices, note that $$ l = -\frac{2 \left(c (2 a b+1)^{2/3}+b (2 a c+1)^{2/3}\right)}{3 (2 a b+1)^{2/3} (2 a c+1)^{2/3} (b c+1)} $$ Using $a, b, c \ge 0$ at the bottom and $a, b, c \le 4$ at the top shows that $l > -55$.

This in fact proves (not just verifies) our conjecture that there is only one solution of $\nabla g=0$ (if there is no bug in the library).

However, to be sure that maximum points do not appear on the boundary, we still need to check, for instance $a=0$. This reduces to find maximum of $$ \sqrt[3]{2 (4-c) c+1}-3 \sqrt[3]{3}+2, $$ which is $$ 2-3 \sqrt[3]{3}+3^{2/3} < 0 $$ when $c = 2$.

I deleted my previous solution. I give another solution.

WLOG, assume that $c = \min(a,b,c)$.

Since $x\mapsto \sqrt[3]{x}$ is concave on $(0, \infty)$, we have \begin{align} \sqrt[3]{1+2ac} + \sqrt[3]{1+2ba} + \sqrt[3]{1+2cb} &\le 2\sqrt[3]{\frac{1+2ac + 1+2cb}{2}} + \sqrt[3]{1+2ba}\\ &= 2\sqrt[3]{1+ac + cb} + \sqrt[3]{1+2ba}. \end{align} It suffices to prove that $$2\sqrt[3]{1+ac + cb} + \sqrt[3]{1+2ba} \le 3\sqrt[3]{3}.$$

We split into two cases:

1) $ba \le 1$: We have $ac + cb \le 2ba\le 2$ and thus $2\sqrt[3]{1+ac + cb} + \sqrt[3]{1+2ba} \le 3\sqrt[3]{3}$.

2) $ba > 1$: From $a+b+c+abc = 4$, we have $c = \frac{4-a-b}{ab+1}$. Also, we have $a+b \ge 2\sqrt{ab} > 2$. Thus, we have $$ac + cb = \frac{(4-a-b)(a+b)}{ab + 1} = \frac{4 - (a+b - 2)^2}{ab+1} \le \frac{4 - (2\sqrt{ab} - 2)^2}{ab+1}.$$ Thus, it suffices to prove that $$2\sqrt[3]{1+ \frac{4 - (2\sqrt{ab} - 2)^2}{ab+1}} + \sqrt[3]{1+2ba} \le 3\sqrt[3]{3}.$$ Let $x = \sqrt{ba}$. Then $1 < x \le 2$. It suffices to prove that for $1< x\le 2$, $$2\sqrt[3]{1+ \frac{4 - (2x - 2)^2}{x^2+1}} + \sqrt[3]{1+2x^2} \le 3\sqrt[3]{3}$$ or $$2\sqrt[3]{1 + \frac{2(3x-1)(1-x)}{3(x^2+1)}} + \sqrt[3]{\frac{1+2x^2}{3}} \le 3.$$ Note that $\sqrt[3]{1+u} \le 1 + \frac{u}{3}$ for all $u > -1$ since $(1+\frac{u}{3})^3 - (1+u) = \frac{1}{27}u^2(u+9)$. Thus, it suffices to prove that $$2\left(1 + \frac{1}{3} \frac{2(3x-1)(1-x)}{3(x^2+1)}\right) + \sqrt[3]{\frac{1+2x^2}{3}} \le 3 \tag1$$ or $$\sqrt[3]{\frac{1+2x^2}{3}} \le \frac{21x^2-16x+13}{9(x^2+1)}$$ or $$\frac{1+2x^2}{3} \le \left(\frac{21x^2-16x+13}{9(x^2+1)}\right)^3$$ or $$\frac{2(x-1)^2(-243x^6-486x^5+3051x^4-3996x^3+4527x^2-2102x+977)}{729(x^2+1)^3}\ge 0.$$ It is easy to prove that $-243x^6-486x^5+3051x^4-3996x^3+4527x^2-2102x+977 > 0$ for $1 < x\le 2$.
EDIT: Indeed, By letting $x = 1 + v$ for $0 < v \le 1$, we have \begin{align} &-243x^6-486x^5+3051x^4-3996x^3+4527x^2-2102x+977\\ =\ & -243v^6-1944v^5-3024v^4-1512v^3+2340v^2+3280v+1728\\ \ge \ & -243v^2-1944v^2-3024v^2-1512v^2+2340v^2+3280v+1728\\ =\ & -4383v^2+3280v+1728\\ \ge \ & -4383v+3280v+1728\\ =\ & -1103v+1728\\ >\ & 0. \end{align} We are done.