How do we design a neural network such that the $L_1$ norm of the outputs is less than or equal to 1?
Penalty (barrier function) is perfectly valid and simplest method for simplex type constraint (L1 norm is simplex constraint on absolute values). Any type of barrier function may work, logarithmic, reciprocal or quadratic. All of them supported by any major framework(pytorch, tensorflow), just add them to loss function. You would need some hyperparameter tuning for the scale factor of penalty.

There is more efficient, though more complex way to do it. Instead of putting constraint you can automatically output value wich satisfy simplex constraint:

Assume that L1 norm constraint is $\left \|v\right \|_1 \leq 1$, $v \in \mathbb{R}^n$


put $sigmoid(v_i)$ activation on output to norm elements to [-1, 1]
add slack (fake) variable element $v_{n+1} = 1 - \sum_{1}^{n}  v_i $
project new $v{}'\in \mathbb{R}^{n+1}$, $v{}'_i = |v_i|,1\leq  i  \leq n+1$ onto unit simplex with standard algorithm (also here)


Backpropagation of last step may require differentiable sorting, which is missing in most of frameworks, you may have to look for open sourced implementation, for example extract it from here or use some automatic differentiation package. Both require some substantial code reading/debugging. However in my experience assuming constant $\Delta$ also works in many cases, in that case  differentiable sorting is not needed. Intuition behind constant $\Delta$  is that  $\Delta$  could be chosen such way that there is some interval on which it's value doesn't affect sorting order.