Documentation: penalty for metric constraint violation

FLAML offers `metric_constraints`, which are supposed to guide the optimizer towards narrowing (as much as possible) the search space to a region specified by the user. [Here](https://github.com/microsoft/FLAML/blob/c48babd02fa9b3252cf49100599fb99c70794591/flaml/tune/searcher/blendsearch.py#L418) appears to be the main place in the code where metric constraint handling is implemented.

As per my understanding of the code, this is the basic idea of the algorithm:
- the optimizer is driven to stay away from unfeasible regions by applying a very high penalty on the metric that is used as an optimization objective
- it's desirable for the penalty to be proportional to the distance to the feasible region, so that even outside of the bounding box there's still a gradient for the optimizer to learn to follow, so that the optimization process leads towards within the user-specified bounds
- because the constrained metrics can have different magnitudes, and they are not known beforehand, the penalty for each metric should be updated based on observed metric values until it proves to provide enough of an incentive for the optimizer to avoid the unfeasible region

To this end, I think the steps taken look something like this (as much as I could understand the code):
- in the beginning, until the first combo is found that satisfies the constraints, a huge penalty (`D x 10^10`) is applied on the objective value
- as per the above formula, the penalty is also proportional to `D` which represents by how much the user-specified thresholds were missed (`D` is the difference between the metric value reported as part of the results of a trial and the constraint the user specified for that metric)
- once a feasible config is found (i.e., one that doesn't violate any constraints), the penalty coefficient is reduced from `10^10` to `1`. This means that the score of the next unfeasible combo after some feasible ones have been found will only be penalized by `D1` instead of` D1 x 10^10`
- every time a constraint is violated, the penalty is increased by `D`, so the score of the second unfeasible combo will be penalized by `(1 + D1) x D2`, the score of the third one, by `(1 + D1 + D2) x D3`, the score of the fourth one by `(1 + D1 + D2 + D3) x D4` and so on. The penalty is increased until `10^10` is reached again

Now, the above procedure at least looks interesting and seems to follow sound desiderata, but has anyone actually verified empirically that it works? Is there some theoretical proof that it converges to the feasible regions? Is there a risk for the penalties to explode and get too large? How does this affect the affinity of the optimizer towards exploring configurations on the edge of the feasibility (which might well be the best performing ones)?

In other words, is there a scientific paper in which this procedure is described into further detail? Or at least some place in the docs where it is addressed?