The lengths of the broken stick pieces $(Y_1, Y_2, Y_3, Y_4)$ are spacings of ordered uniform. The joint distribution is equivalent to
$$ \left(\frac {X_1} {\sum_{i=1}^4 X_i},
\frac {X_2} {\sum_{i=1}^4 X_i},
\frac {X_3} {\sum_{i=1}^4 X_i},
\frac {X_4} {\sum_{i=1}^4 X_i}\right)$$
where $X_i$ are iid exponential random variables. If we ordered the above spacings, the distribution of $(Y_{(1)}, Y_{(2)}, Y_{(3)}, Y_{(4)})$ is equivalent to
$$ \left(\frac {X_1/4} {\sum_{i=1}^4 X_i},
\frac {X_1/4 + X_2/3} {\sum_{i=1}^4 X_i},
\frac {X_1/4 + X_2/3 + X_3/2 } {\sum_{i=1}^4 X_i},
\frac {X_1/4 + X_2/3 + X_3/2 + X_4} {\sum_{i=1}^4 X_i}\right)$$
Any $3$ of the $4$ pieces cannot form a triangle if and only if
$$ \begin{cases}
Y_{(1)} + Y_{(2)} < Y_{(3)} \\
Y_{(1)} + Y_{(2)} < Y_{(4)} \\
Y_{(1)} + Y_{(3)} < Y_{(4)} \\
Y_{(2)} + Y_{(3)} < Y_{(4)}
\end{cases}$$
Note that the second inequality is implied by the first, and the third inequality is implied by the fourth. So the probability of no triangle being formed is
$$ \begin{align}
&\Pr\{ Y_{(1)} + Y_{(2)} < Y_{(3)}, Y_{(2)} + Y_{(3)} < Y_{(4)}\} \\
=& \Pr\Bigg\{\frac {X_1} {4} + \frac {X_1} {4} + \frac {X_2} {3} <
\frac {X_1} {4} + \frac {X_2} {3} + \frac {X_3} {2}, \\
& \frac {X_1} {4} + \frac {X_2} {3} + \frac {X_1} {4} + \frac {X_2} {3}
+ \frac {X_3} {2} <
\frac {X_1} {4} + \frac {X_2} {3} + \frac {X_3} {2} + X_4\Bigg\} \\
=& \Pr\{X_1 < 2X_3, 3X_1 + 4X_2 < 12X_4\} \\
=& \int_0^{\infty} \Pr\{2X_3 > x\}\Pr\{12X_4 - 4X_2 > 3x\}e^{-x}dx
\end{align}$$
Note that
$$ \begin{align}
&\Pr\{12X_4 - 4X_2 > 3x\} \\
=& \int_0^{\infty}\Pr\{12X_4 - 4u > 3x\}e^{-u}du \\
=& \int_0^{\infty} e^{-(3x+4u)/12} e^{-u}du \\
=& e^{-x/4}\int_0^{\infty} e^{-4u/3}du \\
=& \frac {3} {4}e^{-x/4}
\end{align}
$$
So the integral become
$$ \int_0^{\infty} e^{-x/2}\frac {3} {4} e^{-x/4}e^{-x}dx
= \frac {3} {4} \int_0^{\infty} e^{-7x/4} dx = \frac {3} {4} \times \frac {4} {7} = \frac {3} {7} $$
This is not a very elegant way but at least it is doable. Looking forward to someone to post a better solution.
The probability is exactly
$$\frac{5^2\cdot 29^2}{2^2\cdot 7\cdot 11\cdot 13\cdot 19}=\frac{21025}{76076}\approx 0.2763683684736316$$
We're trying to compute the volume of a rather complicated region in 5D space, but one which is ultimately bounded by facets with rational coordinates, so we know the answer will have to be rational in the end. We just want to find a way of cutting up this big 5D polytope inside the unit cube into pieces we can integrate.
The first simplification we'll make is to assume that the five cuts are ordered smallest to largest, without loss of generality, i.e. $0\le x_1\le x_2\ldots\le x_5\le 1$. This means we'll need an extra factor of $120$ at the end, but we can now talk about the lengths of the six segments as $x_1, x_2-x_1, \ldots, x_5-x_4, 1-x_5$.
In this framing, every possible way to arrange the segments into triangles corresponds to a set of six linear inequalities, one for each side of each triangle asserting that that side is at least as large as the sum of the other two sides on its triangle.
So for any given set of ways we can arrange the segments into two triangles, we'll have a convex region described by some half-spaces corresponding to the breakpoints in the large segment such that we can make the pieces into triangles in all of those ways (e.g. we might ask about the splits such that segments (1,2,3) and (4,5,6) make triangles, as well as segments (1,4,5) and (2,3,6)).
Ultimately, we want to know the volume of the region for which at least one of these $\frac{{}_6C_3}2=10$ triangle pairs are possible, so we'll use some inclusion-exclusion counting:
$$\text{Vol}(\text{at least one region}) = \sum_{\emptyset\neq S\subset \{1,2,\ldots,10\}}(-1)^{|S|+1}\text{Vol}\left(\bigcap_{i\in S} \text{triangle pair $i$ is possible}\right)$$
For each of these $1023$ regions* we want to compute, we just need to compute the volume of a particular region in 5-space bounded by a set of rationally-parametrized half spaces. This can be done numerically in Python using scipy.spatial.HalfspaceIntersection and scipy.spatial.ConvexHull. To convert a floating-point volume $V$ to an exact value, I looked at the coordinates of the intersections of the region, found a scaling factor $F$ that made those intersections lie at integer points, and then searched for small rational values $q$ such that $qFV$ was very close to an integer. In all cases there was a rational value within a factor of $10^{-14}$ of the floating-point approximation whose denominator had no prime factors over $19$. Then I just used exact integer arithmetic on the resulting fractions to obtain the final answer. (As supporting evidence for these conversions to rational numbers, one of these regions had a volume of $2491157/41409688320$, which is much more complicated than the final sum but which one wouldn't expect to cancel out if it were a random coincidence.)
Given how comparatively nice the final answer ends up being, I suspect there is some deeper reason that things work out so well; the numerator and denominator are both suspiciously nice compared to the volumes of the intermediate regions in this sum. I have no idea what approaches might shed light on the reasons behind this, though.
*Because of the symmetry between the six segments we break the stick into, many of these regions are congruent - we only end up needing to consider 18 different possibilities.
Best Answer
Your solution does not fulfill the random protocol described in the problem. The second break is not at a random point in $[0,1]$, but at a random point in the longer piece of the first break.
Choose the first breaking point at $x\in\,\bigl]{1\over2},1\bigr]$ and the second at $xy$ with $y\in\,\bigl]{1\over2},1\bigr]$. The variable $(x,y)$ is then uniformly distributed in the square $Q:=\bigl]{1\over2},1\bigr]^2$ of area ${1\over4}$. In this way we obtain three pieces of length $$\ell_1=xy,\quad \ell_2=x-xy=x(1-y),\quad \ell_3=1-x\ .$$ The $\ell_i$ can be the sides of a triangle iff all of them are $<{1\over2}$. For $\ell_2$ and $\ell_3$ this is automatically the case, and we are left with the condition $\ell_1=xy<{1\over2}$. Let $S$ be the corresponding part of $Q$. The required probability then is $$p={{\rm area}(S)\over{\rm area}(Q)}=4\int_{1/2}^1\left({1\over2x}-{1\over2}\right)dx=\ldots=2\log2-1\approx0.3863\ .$$