Suppose you have $n$ physical objects that you want to determine the mass and charge of. You do not have access to any reference object with known mass and charge (that also includes things like the Earth which we typically assume to be neutral: we would like to determine its charge empirically, not assuming it as an additional axiom). The whole point of this set up is that we start from a mathematical model of electrostatic and gravitational interactions, and we want to discover a set of experiments that allow us to fit all the masses and charges of elementary particles (or larger objects such as the Earth). The mathematical model is described next.
Assumptions
We assume that there is a well-defined system of units for length and time.
For simplicity, we will assume that all objects can be treated as point masses and point charges for the purposes of calculating gravitational and Coulomb interactions. This allows us to represent the $i$-th object by a tuple $(m_i, M_i, Q_i)$ corresponding to its inertial mass $m_i$, gravitational mass $M_i$ and electrical charge $Q_i$ (for now we do not assume that $m_i = M_i$, as we hope this can be discovered from the experiments).
We assume that the gravitational and electrostatic interactions obey Newton's and Coulomb's laws, respectively, and that the values of $k$ and $G$ are fixed (knowing $k$ and $G$ amounts to a choice of units for charge and gravitational mass). We will neglect the electromagnetic interactions due to magnetic fields, radiation phenomena and other relativistic effects for simplicity (or you can view this as the observation that the corrections are on the order of $v/c \ll 1$, so that they lie within experimental error given the instruments available to us).
We also assume that Newton's laws are applicable, and that we work in an inertial frame. In that case, all the kinematic observations are summarized by the following system of $n$ ODEs:
$$m_i\ddot{\mathbf{x}}_i = \sum\limits_{j \neq i} \frac{kQ_iQ_j - GM_iM_j}{|\mathbf{x}_i-\mathbf{x}_j|^3}(\mathbf{x}_i-\mathbf{x}_j)$$
Finally, we assume that $M_i \gt 0$ and $m_i \gt 0$ for all $i$.
Define the interaction constants $\alpha_{ij} = \frac{1}{m_i}(kQ_iQ_j - GM_iM_j)$.
Our system of ODEs now becomes:
$$\ddot{\mathbf{x}}_i = \sum\limits_{j \neq i} \alpha_{ij}\frac{\mathbf{x}_i-\mathbf{x}_j}{|\mathbf{x}_i-\mathbf{x}_j|^3}$$
As you can see, the constants $\alpha_{ij}$ determine completely the kinematics of the system. Using kinematic experiments, we can measure the values of all $\alpha_{ij}$ for $1 \leq i \neq j \leq n$.
Question
Is there a set of kinematic experiments satisfying the assumptions above that allows us to recover the masses and charges of each object uniquely (up to experimental error)? If yes, what is the minimal size of $n$ for which such experiments exist? If not, what additional model corrections, assumptions and measurement capabilities were historically present for the experimenters that defined the system of units?
Comments
A positive answer should not simply mention or describe known historical experiments. The answer should demonstrate (or provide reference to a such a demonstration) that the experiment allows unique determination of all masses and charges participating (up to the resolution of the instruments), without making additional assumptions other than the ones stated above. It would be nice, but not necessary, to also mention whether this demonstration of uniqueness was considered by the experimenters, historically.
A negative answer should prove that for any set of values $\alpha_{ij}$, either there are no solutions, or there is more than one solution (even accounting for the limited resolution), and cite which assumptions were made historically to resolve this ambiguity. It would also be nice to mention whether the experimenters had this problem in mind when designing their experiments.
This is not as obvious as it might look like. Remember that all the experiments can determine are the values of the $\alpha_{ij}$. Finding the $m_i, M_i$ and $Q_i$ from the $\alpha_{ij}$ involves solving a constrained system of $n^2 - n$ polynomial equations in $3n$ unknowns (where some of the equations might not actually be independent). The experimenters do not know the value of any of these parameters prior to the measurements (otherwise, this is simply pushing the problem further down). It is likely a necessary but not sufficient condition that you need at least as many equations as unknowns, so you would need $n^2 - n \geq 3n$, or at least $n \geq 4$. This involves solving a system of at least $12$ polynomial equations in $12$ unknowns! The experimenters could have been lucky and found one simple solution to the system of equations (which would be amazing on its own), but how can we be certain that the solution is unique without doing some nasty computation?
My intuition is that this simple model is a bit too simple, in the sense that charge and mass occur in very symmetric ways in the equations, so isolating charge from mass is tricky (either impossible, or algebraically difficult). With added corrections like magnetic forces, the charges show up separately in different terms that don't have the same inverse-square law form, so it is more likely that some measurable parameters depend only on the $Q_i$. So I am suspecting that these corrections were likely necessary for the experimenters to have unambiguous results. But these corrections are on the order of $v/c$ so the precision of the instruments would have to be quite good. Also note that magnetic interactions are difficult to model, because they rely on both the Lorentz force $\mathbf{F} = q\mathbf{v}\times\mathbf{B}$ and the Liénard-Wiechert formulas for $\mathbf{B}$ in terms of the source particles' trajectories.
