Hello @samuele
Have you tried debugging the code to see how the matrix and vector of the Gauss-Newton problem look ?
The actual errors of the measurements should not lead to a singular problem since they only impact the residuals computation.
The errors in input (the standard deviations you give for your measurements) do intervene in the computation of the Jacobian (derivatives matrix) but only as weights, so they should not lead to a singular problem neither.
The “singular problem” comes from the fact that your system cannot be decomposed with the QR decomposition in the Gauss-Newton algorithm.
It means that the pseudo-inverse of the Jacobian J of your problem cannot be computed.
There’s a threshold in the Hipparchus Gauss-Newton class to check if the problem is singular but it is fixed to 1e-11, so you cannot play around with that and see the effect.
The Jacobian J is such as: J_{i,j}= \frac{\partial m_i}{\partial y_{0j}}
Where m_i is the ith measurement and y_{0j} the jth component of the initial orbit y_0. For example if y_0 is in Cartesian coordinates: y_{01}=X_0, y_{02}=Y_0... and so on.
This matrix is decomposed as: J_{i,j}= \frac{\partial m_i}{\partial y_{0j}} = \sum_{k=1}^{k=6}{\frac{\partial m_i}{\partial y_{k}}.\frac{\partial y_k}{\partial y_{0j}}}
Where y_k is the state vector evaluated at measurement’s date.
In your case, y_k is states[0] in the theoreticalEvaluation method.
And the \frac{\partial m_i}{\partial y_{k}} are set in your derivative array.
The \frac{\partial y_k}{\partial y_{0j}} are the components of the so-called state transition matrix that is computed through integration along with the orbit.
In your case \frac{\partial m_i}{\partial y_{k}} seems to be a constant: \frac{\partial m_i}{\partial y_{k}} = [n_x, n_y, n_z, 0, 0, 0] = [\vec{n}, \vec{0}].
So I’m afraid that if the measurements are too close in time and the orbit model vary slowly (which seems to be the case since the period is 7 days), then the components of the state transition matrix \frac{\partial y_k}{\partial y_{0j}} do not change much between two measurements, and, since the \frac{\partial m_i}{\partial y_{k}} are constant, the two lines of the Jacobian are almost identical, which is a classic cause of singular problem (non pseudo-inversible matrix) exception.
That being said, I may be wrong, it’s just a feeling since I don’t know the actual data.
Could you explain us a bit more what is behind your vector \vec{n} and why it can be considered constant in your case ?
Can you try using measurements that are further apart (say with 90° anomaly separation from each other) to see what you get ?
I suggest you try to use the automatic differentiation feature of Orekit to compute the derivatives of your measurement. It may help having different derivatives for different measurements. And could also show that there may actually be a component on the velocity that wasn’t intuitively easy to apprehend.
Last question, I don’t see any light travel-time correction in your code, is that normal ? I’m guessing it must already be corrected somewhere upstream otherwise you wouldn’t manage to get some orbit determinations to work 
Hope this helps,
Maxime