NMaximize scrambles constraints

There have been lots of questions around regarding the issue that NMaximize/NMinimize may give results - or at least evaluate the objective function - although constraints are violated. The explanation for this was that constraints correspond to a "badness" quantifier and are not rigidly imposed.

However, the following issue puzzles me completely. I stumbled upon it when I tried to implement constraints regarding the eigenvalues of a Hermitian matrix and the constraint function complained about (numerically significant) imaginary values. Here's an extremely condensed example:

dummy[x__?NumericQ] := (Print["called objective"]; Abort[];);
dummyConstr[x_List?(AllTrue[#, NumericQ, 2]&)] := (Print[x]; Abort[];);
NMaximize[{dummy[a, br, bi, c],
dummyConstr[{{a, br + I bi}, {br - I bi, c}}] >= 0.}, {a, br, bi, c}]


So here, I don't actually care about the optimization, I just want to find out about the values the constraint function is called with. The output is (probably depending on a random seed, but anyways):

> {{0.652468, 0.682813 + 0.63307 I}, {0.935202 - 0.566352 I, 0.976188}}
> $Aborted  The constraint function is evaluated. However, its parameters have a structure that is Hermitian, whatever the (according to the documentation, real-valued) parameters a, br, bi and c may be. But clearly, the matrix that dummyConstr is called with is not! Solve[{{a, br + I bi}, {br - I bi, c}} == (* the above matrix *)] > {{a -> 0.652468 + 0. I, bi -> 0.599711 + 0.126195 I, br -> 0.809008 + 0.0333593 I, c -> 0.976188 + 0. I}}  Apparently, NMaximize not only allows violations to the constraints, but it also evaluates them for complex parameters. While this is undocumented, I could still live with it (but I wanted to mention it here separately, since I did not find a notice of this before). And if I remember correctly, there are optimization algorithms that require constraint functions to be analytic, so why not evaluate them with complex numbers? However, now let me change the matrix so that it is of unit trace. NMaximize[{dummy[a, br, bi, c], dummyConstr[{{a, br + I bi}, {br - I bi, 1 - a}}] >= 0.}, {a, br, bi}] > {{0.652468, 0.682813 + 0.63307 I}, {0.935202 - 0.566352 I, 0.0238122}} >$Aborted
Solve[{{a, br + I bi}, {br - I bi, c}} == (* the above matrix *)]
> {{a -> 0.652468 + 0. I, bi -> 0.599711 + 0.126195 I, br -> 0.809008 + 0.0333593 I, c -> 0.0238122 + 0. I}}


Note that here when solving, I did not say 1 - a to show that the absence of a solution is not a numerical problem. As you can easily see, based on the value for a, c should be about 0.35, but it is 0.02 instead. And I cannot understand this at all. The only way to produce such an output is if during the evaluation of the constraint function, the value for the parameter a changes. How can this be the case? And even if it is a bug, it must be a very strange one.

Tested with Mathematica 11.3 and 12.