$f(x)$ is monotonically increasing in $[0,1]$, $0 \le f \le 1$ and $\int_0^1 (f(x) - x) \mathrm{d}x = 0$. Prove that $\int_0^1|f(x)-x|\mathrm{d}x \le \frac{1}{2}$.
It's easy if $f(x) \ge x$ in $[0,1]$. And even in $[a,b]$ we have $\int_a^b |f(x)-x|\mathrm{d}x \le \frac{(b-a)^2}{2}$. But the zero points of $f(x) - x$ may be infinitely many. This is where difficulty exists.
Here we present a bit different, calculus-themed approach. In this answer, we will assume that $f : [0, 1] \to [0, 1]$ is monotone-increasing. We also write $I(f) = \int_{0}^{1} |f(x) - x| \, \mathrm{d}x$ for brevity.
Step 1 - Proof under extra assumptions. Assume further that $f$ is piecewise-smooth, $f(0) = 0$, and $f(1) = 1$. Then by the formula $\int |x|\,\mathrm{d}x=\frac{1}{2}x|x|+\mathsf{C}$, we have
$$ \int_{0}^{1} |f(x) - x|(f'(x)-1) \, \mathrm{d}x = \left[ \frac{1}{2}|f(x)-x|(f(x)-x) \right]_{0}^{1} = 0. $$
In particular,
$$I(f) = \frac{1}{2}\int_{0}^{1} |f(x) - x|(f'(x)+1) \, \mathrm{d}x.$$
Now pick $\alpha \in [0, 1]$ so that $f(\alpha) + \alpha = 1$. (This is possible since $x \mapsto f(x)+x$ increases from $0$ to $2$. Then by triangle inequality,
\begin{align*} \int_{0}^{\alpha} |f(x) - x|(f'(x)+1) \, \mathrm{d}x \leq \int_{0}^{\alpha} (f(x)+x)(f'(x)+1) \, \mathrm{d}x = \frac{1}{2}. \end{align*}
Similarly, by writing $|f(x)-x| = |(1-f(x))-(1-x)| \leq (1-f(x)) + (1-x)$, we get
\begin{align*} \int_{\alpha}^{1} |f(x) - x|(f'(x)+1) \, \mathrm{d}x \leq \int_{\alpha}^{1} (2-f(x)-x)(f'(x)+1) \, \mathrm{d}x = \frac{1}{2}. \end{align*}
Therefore $\int_{0}^{1} |f(x)-x| (f'(x)+1) \, \mathrm{d}x \leq 1$, which in turn implies $I(f) \leq \frac{1}{2}$ as required.
Remark. Let $\gamma(t) = (f(t)+t, f(t)-t)$. Then $\int_{0}^{1} |f(t)-t|(f'(t)+1)\,\mathrm{d}t = \int_{\gamma} |y|\,\mathrm{d}x$ computes the area between the path $\gamma$ and the horizontal axis. Note that $\gamma$ is essentially the $-45^\circ$-rotation of the graph $y = f(x)$ up to scaling.
Then the above bounds immediately follow from the fact that the graph of $\gamma$ defines a function on $[0, 2]$ which is squeezed between lines $y = \pm x$ and $y = \pm (2-x)$.
Step 2 - General case. For the general case, let $f_n$ be the linear interpolation of the points
$$(0, 0), \quad (\tfrac{1}{n},f(\tfrac{1}{n})), \quad \cdots, \quad (\tfrac{n-1}{n}, f(\tfrac{n-1}{n})), \quad (1, 1).$$
Then by monotonicity,
\begin{align*} |I(f_n) - I(f)| &\leq \int_{0}^{1} |f_n(x) - f(x)| \, \mathrm{d}x = \sum_{k=1}^{n} \int_{\frac{k-1}{n}}^{\frac{k}{n}} |f_n(x) - f(x)| \, \mathrm{d}x \\ &\leq \frac{1}{n}\left( [f(\tfrac{1}{n})-0] + \sum_{k=2}^{n-1} [f(\tfrac{k}{n}) - f(\tfrac{k-1}{n})] + [1-f(\tfrac{n-1}{n})] \right) \\ &= \frac{1}{n}, \end{align*}
hence $I(f_n) \to I(f)$ as $n\to\infty$ and the desired inequality $I(f) \leq \frac{1}{2}$ follows from the previous step.
In case you have some familiarity with measure theory:
Take an increasing sequence of simple functions $f_k \rightarrow f$ converging pointwise. By the monotone convergence theorem we have $\lim_{k \rightarrow \infty} \int_0^1 f_k(x) dx = \int_0^1 f(x) dx$. So it suffices to show the result for simple functions. The functions $f_k(x) - x$ will have only finitely many zeroes, so you know how to do this.
Here is a sketch of a different argument. Let $K$ be the set of monotone increasing functions $f:[0,1]\to[0,1]$ that are right continuous and have limits on the left (the càdlàg functions). One may as well restrict the $f$ of the problem statement to be in $K$, as any $f$ can be modified at its (countably few) discontinuities to make it càdlàg and while preserving all the integrals. Now, $K$ is convex, and its extreme points are the cumulative distribution functions of the point masses at points $a\in[0,1],$ together with the zero function. In particular, the simple functions constant on $[0,a)$ and on $[a,1]$. (This claim is, in effect, that the point masses and the zero measure are the extreme points of the sub-probability measures on $[0,1]$.) Now, by a theorem of Dubins (see also), the extreme points of $K\cap L$, where $L$ is the set of $f$ for which $\int_0^1 f = 1/2$ are all convex combinations of at most two extreme points of $K$. Which is to say, simple functions with at most 2 discontinuities. Since $f\mapsto\int_0^1|f(x)-x|dx$ is continuous and convex, the maximum is attained at an extreme point of $K\cap L$. As the comments indicate, $\int_0^1|f(x)-x|dx$ is equal to the area of the union of finitely many triangles, the sum of whose heights is $\le 1$ and the sum of whose lengths is $\le 1$ and hence have total area $\le 1/2$.
