Short version;
Nonlinearstatic is complex and the incremental formulation is problematic. Moving the load level from the engineering model to the Load will probably make things a lot better.

Long version;
Almost all of the complexity in nonlinearstatic comes from indirect load control and the incremental load formulation.
In particular the incremental formulation leads to

1. Much more complex code. Load balancing, changing dirichlet b.c.s etc. are very complex. It is the reason for the "VM_RhsTotal" hack in order to be able to store the values in nodes. All is only due to the incremental formulation.

2. Forces one to keep track of the prescribed vectors as well (for reaction forces and changing b.c.), even if it turns out to never be used.

3. It means the problem can never change in size; No XFEM. No contact simulation. No remeshing or adaptivity (No error estimates).

4. VM_Incremental for use with Loads is currently broken since it doesn't respect "isImposedTimeFunction" and other corner cases.

But I understand the need to be able to do the simulations of the form;
Residual = F_0 + F_ref * lambda - F_internal = 0.
However, I don't think this actually implies that one needs to do an incremental formulation.

By introducing a load-level parameter on the Load class itself + a new vector type to assemble "ReferenceExternalForcesVector" which is used whenever the load-level is zero. After analysis has converged, the reached load-level value is set for all the Loads that were actived and had zero load level set previously.
Whether or not to use CALM or NRSolver wouldn't be directly controlled by the user, but would be automatically detected if the reference load vector is nonzero. This would solve all 4 concerns I have with nonlinearstatic.
In this formulation you could simply specify a certain load as being the reference load with a single flag. If you want the reference load to start the second step of the analysis, then just apply an "isImposedTimeFunction".

Edit: I've realized that nonlinearstatic doesn't actually seem to support changing essential boundary conditions either, though more or less all the required code is already there through the pack/unpack functions (though these are parallel mode only).