I think I have found the answer :). TheLagrange multipliers technique increases the condition number of the stiffness matrix. This is caused by the different physical units of the Lagrange multipliers (force) and displacements (length). The condition number of stiffness is decreased by normalizing the Lagrange multipliers with the maximum diagonal entries of stiffness matrix. For more info see "Jorg Unger, Stefan Eckardt, Multiscale modelling of concrete, 2011, page 371".

Using Lagrange multipliers, you generally impose an exact verification of your boundary conditions constraints (primal method). With a penalty method, you relax your boundary condition by assuming an error on the condition, scaled by the penalty coefficient (dual method).

It is quite classical that the application of Largange multipliers can cause lock-in in your model (overconstrained interfaces with incompatible meshes). Besides, the resolution algorithms with Largange multipliers are specific and can be numerically sensitive (fixed point methods).

On the other hand the application of a penalization allows using standard modified Newton Raphson algorithms whose convergence patterns are better known in the literature.

Be aware that if your penalization coefficient is too small, the accepted error on the verification of your constraint will be too large and your results may become inaccurate. It is classically necessary to check the solution convergence as function of the penalty coefficient (or to have a calibration method) to acknowledge your results accuracy.

I'm using Newton-Raphson solver for my coupled nonlinear system using mesh-less method. Unfortunately in mesh-less methods the only way to control the boundaries are Lagrange multipliers and beta penalty method and I'm having this non-optimal converging rate problem now. I don't know if it's considered to be OK. I don't care if it doesn't converge optimally but I'm afraid it is converging to a non-solution field.