I'm trying to approximate the first-order loss function $$\mathbb{E}[max(d-y,0)]=\sigma \cdot (\phi(\frac{d-\mu}{\sigma})- \frac{d-\mu}{\sigma} \cdot (1- \Phi (\frac{d-\mu}{\sigma})),$$ where $d\in \mathbb{N}$ and $ y \sim \mathcal{N}(\mu,\sigma^2)$.
I need the approximation to solve an optimization problem. There is a discussion about lower & upper bounds in https://www.sciencedirect.com/science/article/pii/S0096300314000563 but I don't understand it.
Is there an easy way to approximate it with python, or in general to get a good and not to complex approximation of a non-linear function?