Now, I understand that if I draw the values of A_i randomly, I would not have a well-defined integer particle number so this might be a problem. But what if the A_i are obtained through some unitary transformation, such that \sum_i A_i\in \N? Is there a possible way to do this? Or do I just need to do the simulation without QN conservation.
I would guess that, most likely, your state just doesn’t conserve the symmetry, or rather that our system cannot be sure that it conserves it from the way you are inputting the numbers.
(A semi-related example of this same kind of limitation is that sometimes our users want to write the Heisenberg model as S^x_i S^x_j + S^y_i S^y_j + S^z_i S^z_j and we have to ask them to write it as \frac{1}{2} S^+_i S^-_j + \frac{1}{2} S^-_i S^+_j + S^z_i S^z_j, with the reason being that the S^+ and S^- operators individually have a well-defined “flux” whereas operators like S^x do not.)
I think the route you mentioned of making the state by unitary rotations sounds like a very good way to go. You could start from a product state that you and our system can be sure respects the symmetry, then also make the single-site operators as sums of operators that also individually have well-defined “flux”. If all the steps work, you will have your result but if along the way you find actually there is no such single-site operator for general values of A_i you will know that it’s not possible.