Example 4 Let real numbers $a \geqslant b \geqslant c \geqslant d>0$. Find the minimum value of the function
$$
\begin{array}{l}
f(a, b, c, d) \\
=\left(1+\frac{c}{a+b}\right)\left(1+\frac{d}{b+c}\right)\left(1+\frac{a}{c+d}\right)\left(1+\frac{b}{d+a}\right)
\end{array}
$$
Solution: It is clear that $f$ has no upper bound, because when $a=b=1, c+d \rightarrow 0$, $f \rightarrow+\infty$.
Also note that $f$ is a homogeneous function of degree zero, and when $a=b=c=d$, $f=\left(\frac{3}{2}\right)^{4}$.
Let $\sum$ denote the cyclic sum, and $\Pi$ denote the cyclic product.

We now prove: For any positive numbers $a$, $b$, $c$, $d$ satisfying the conditions, we have
$$
f(a, b, c, d) \geqslant\left(\frac{3}{2}\right)^{4},
$$

which is to prove
$$
\prod \frac{a+b+c}{3} \geqslant \prod \frac{a+b}{2} .
$$

From the given, we have $\frac{a}{d} \geqslant \frac{b}{d} \geqslant \frac{c}{d} \geqslant 1$.
Let $\frac{a}{d}=1+6 x, \frac{b}{d}=1+6 y, \frac{c}{d}=1+6 z$. Then $x \geqslant y \geqslant z \geqslant 0$.
Equation (1) becomes
$$
\begin{array}{l}
(1+2(x+y+z)) \prod(1+2(x+y)) \\
\geqslant(1+3(x+y))(1+3(y+z)) . \\
\quad(1+3 z)(1+3 x) .
\end{array}
$$

Activate a constant, and redefine 1 as $t$, and set
$$
\begin{aligned}
\varphi(t)= & (t+2(x+y+z)) \prod(t+2(x+y)), \\
\psi(t)= & (t+3(x+y))(t+3(y+z)) . \\
& (t+3 z)(t+3 x) .
\end{aligned}
$$

Then $\varphi(t)$ and $\psi(t)$ are both quartic polynomials in $t$, and $p(t)=\varphi(t)-\psi(t)$ is a quadratic polynomial in $t$. Let
$$
p(t)=A t^{2}+B t+C .
$$

To prove that equation (2) holds, we need to prove that $p(1) \geqslant 0$, so we only need to prove that $A \geqslant 0, B \geqslant 0, C \geqslant 0$.
Notice that,
$$
\begin{array}{l}
16(x+z)(x+y+z) \\
\geqslant 16(x+z)(x+2 z) \\
\geqslant 16 \times 2 \sqrt{x z} \times 2 \sqrt{2 x z} \\
=64 \sqrt{2} x z \geqslant 81 \sqrt{x z} .
\end{array}
$$

Then $c=(x+y)(y+z)$.
$$
(16(x+z)(x+y+z)-81 x z)
$$
$\geqslant 0$.
$$
\begin{array}{l}
\text { Also, } 7 \sum x^{3} \cdot 21 x y z \geqslant 20 x y z+x z^{2} \text {, } \\
x^{3} \geqslant x^{2} z \text {, } \\
13\left(x^{2} y+x y^{2}\right) \geqslant 13\left(x^{2} z+x z^{2}\right), \\
\text { so } A=\left(8\left(\sum x\right)^{2}+\sum 4(x+y)(y+z)\right)- \\
\left(\sum 9(x+y)(y+z)+9 x z\right) \\
=3 \sum x^{2}+x y+y z-8 x z \\
\geqslant 3\left(x^{2}+2 z^{2}\right)-7 x z \\
\geqslant 6 \sqrt{2} x z-7 x z \geqslant 0 \text {, } \\
B=\left(8(x+y+z) \sum(x+y)(y+z)+\right. \\
8 \Pi(x+y))-(27 \Pi(x+y)+ \\
27 x z(x+z+2 y)) \\
=8 \sum x(x+y)(x+z)+5 \prod(x+y)- \\
27 x z(x+z+2 y) \\
=8 \sum x^{3}+13\left(x^{2} y+x y^{2}+y^{2} z+y z^{2}\right)- \\
14\left(x^{2} z+x z^{2}\right)-20 x y z \\
\end{array}
$$

Thus, $p(1) \geqslant 0$.
Therefore, the minimum value of the function $f(a, b, c, d)$ is $\left(\frac{3}{2}\right)^{4}$.