Allow gradient descent to go beyond constraint but punish for it Disclaimer: I considered posting this on mathSE but thought maybe this is more fitting here (no pun intended).
Status Quo
I have a list of data points $\{x_k, y_k\}$ and a set of functions $f_i(x)$ with some coefficents $c_i$. I fitted the following function
$$F(x) = \left|\sum \limits_i^n c_i f_i(x)\right|^2 \tag{1}$$
to all data points using gradient descent on least squares w.r.t. to each coefficient $c_i$. The "hard" constraint is that $0 \leq c_i \leq 1$. I have implemented this with a mapping onto a sigmoid function, i.e.
$$\hat{F}(x) = \left|\sum \limits_i^n \frac{1}{1+\exp{\hat{c}_i}} f_i(x)\right|^2 \tag{2}$$
where $\hat{c}_i$ is the appropriately transformed coefficient.
My loss function looks like this:
$$\mathcal{L} = \frac{1}{n}\sum \limits_i^n \left(\frac{y_i - \hat{F}(x_i)}{\sigma_i}\right)^2 \tag{3}$$
where $\sigma_i$ is the error for each data point. This works nicely.
What I want
However, I have an additional "weak" constraint, in the sense that it may be violated, that is $ 0 \leq a_i \leq c_i \leq b_i \leq 1$ for some constants $a_i, b_i$.
What I want to achieve is that every $c_i$ "tries" to stay as much as possible between the recommended values $a_i$ and $b_i$ but is allowed to go beyond it, if it improves the fit, up to a certain point.
A naive approach would be to somehow factor this into the loss function $(3)$, say by also doing a least squares approach to the middle of the interval $[a_i, b_i]$. However, I'm not convinced of this, because the middle point of the interval is not necesarrily the best or most correct one.
My poor attempt/idea is the follwoing function I came up with
$$\zeta_i(x) = e^{-k(x-a)}+e^{k(x-b)}-2 e^{\frac{k}{2}(a_i-b_i)}$$

which still favours the center of the interval, but increases only slightly within the interval and exponentially outside of it. But I'm not sure how to figure out a good scaling factor ($k$), so that the deviation does not dominate or becomes negligible compared the the actual data points.
The unconstraint fit converges to a loss of about $\mathcal{L} \approx 9$.
I'm wondering if this has been explored and maybe an "established" solution exists, or if there's a nice way to obtain what I want.

You could try as penalty for $c_i$:
$$
p(c_i) = \max (0, \left|c_i - \frac{b_i+a_i}{2}\right|-\frac{b_i-a_i}{2}),
$$
which, for $a_i = 1/4$ and $b_i = 3/4$ would look like this:

I don't know how you do gradient descent, but those two points where $p$ is not differentiable should not pose a problem (e.g. just choose as derivative zero if the unlikely event happens that you need the derivative exactly at those points).
This loss function can of course be pulled back to $\hat c_i$. You are free to choose a stronger slope to increase the penalty.
