Is there a robust algorithm for solving material balances? One of the simplest cases of a material balance consists of a steady state balance without reactions. The general balance equations
$$\frac{dm}{dt}=Inputs-Outputs+Generation-Consumption$$
Reduce to simply $Inputs-Outputs=0$
What you have, is a balance equation for each chemical species present, a sum of each mass or mole fractions, some specified data and process specifications, which usually come as linear combinations of the data. For example, the following problem comes from Reklaitis' book Introduction to material and energy balances

The resulting set of equations is:
$$w_{H_2O}^1F^1+w_{H_2O}^2F^2-w_{H_2O}^3F^3-w_{H_2O}^4F^4=0$$
$$w_{NaOH}^1F^1+w_{NaOH}^2F^2-w_{NaOH}^3F^3-w_{NaOH}^4F^4=0$$
$$w_{NaAlO_2}^1F^1+w_{NaAlO_2}^2F^2-w_{NaAlO_2}^3F^3-w_{NaAlO_2}^4F^4=0$$
$$w_{solids}^1F^1+w_{solids}^2F^2-w_{solids}^3F^3-w_{solids}^4F^4=0$$
Where the $w_i^j$ are the mass fractions of component i in stream j. The sum of the mass fractions in each stream equals 1.
$$w_{H_2O}^j+w_{NaOH}^j+w_{NaAlO_2}^j+w_{solids}^j=1$$
Together these equations are equivalent to the total mass balance:
$$F^1+F^2-F^3-F^4=0$$
With the data given in the problem, the system becomes
$$0.98F^1+0.63F^2-w_{H_2O}^3F^3-0.95F^4=0$$
$$0.02F^1+0.11F^2-w_{NaOH}^3F^3-w_{NaOH}^4F^4=0$$
$$0.16F^2+w_{NaAlO_2}^3F^3-w_{NaAlO_2}^4F^4=0$$
$$0.1F^2-0.2F^3=0$$
With $F^1=1000$ and the two specifications $w_{NaOH}^3=0.8w_{NaOH}^4$ and $w_{NaAlO_2}^3=0.8w_{NaAlO_2}^4$
These problems are usually solved with algebraic substitution, but it can be quite tedious and requires some experience to find the best sequence. These books were written before linear algebra was a bigger part of chemical engineering, so I was wondering if there was a more methodic way of solving these systems of equations.
 A: First you have the mass balance equations $WF=0$ where $W$ is a matrix whose $i$th column contains the mole fractions in the $i$th stream. Here I have absorbed the minus signs into $F$ itself, so that the a component of $F$ in an input stream has the opposite sign of one in an output stream. In this system, each row corresponds to a particular substance, so there are $n_{substance}$ equations in $n_{stream}$ unknowns, namely the flow rates.
Next you have the normalization equations $(1,1,\dots,1) W=(1,1,\dots,1)$. After taking the transpose on both sides you have a system where each row corresponds to a particular stream, so there are $n_{stream}$ equations in $n_{stream} n_{substance}$ unknowns, namely the mole fractions.
Let's say for now that $n_{stream}=n_{substance}=n$.
In total that is only $2n$ equations in $n^2+n$ unknowns. You need at least $n^2-n$ more equations to have a hope of pinning down unique weights and flow rates.
As a digression, this accounting is actually not perfect at the moment. In actuality, $F$ cannot be completely determined from $WF=0$, because then it would have to be that $F=0$. However, for now with $W$ square, when $F$ is only determined up to $k$ constants from $WF=0$, we pick up $k$ restrictions on $W$ itself through the requirement that $W$ has $0$ as an eigenvalue with multiplicity $k$. So $WF=0$ still counts as $n_{substance}$ equations overall.
Getting back to the point, in your example problem, you are still short by $4 \cdot 3=12$ equations. Where do these come from? They come from the mole fraction and flow rate relationships that you imposed externally. In total you imposed 14 of those, but two of them were redundant with the normalization equations for the first two streams (since you specified all four fractions of the input streams), so you gained 12 independent equations as expected.
Now what's the structure of solving this whole thing? Your specification is more or less canonical: you have $n$ mass conservation equations, $n$ normalization equations, $n^2-n-1$ mole fraction requirements, and finally $1$ flow rate equation to set the scale of the whole process. The first and last groups, together, pin down $F$, once $W$ is determined. The second and third groups, together with the requirement that $W$ is singular that is inherited from the mass conservation, pin down $W$.
To actually solve the whole thing now, you first symbolically solve the $n^2-1$ linear equations from normalization and externally imposed mole fractions, and then you substitute that into $\operatorname{det}(W)=0$.   This is a 1-variable, potentially $n$th degree polynomial equation, so you can solve it in principle, although it is not obvious that there is a unique solution or even a unique physically meaningful solution (meaning a solution with all nonnegative entries).
Then you turn around and solve $WF=0$ with the additional constraint that sets the scale of the system, and you're done.
If $n_{stream}$ and $n_{substance}$ aren't equal then I'm not entirely sure what happens. Somehow you get a different situation depending on whether $WF=0$ is generically overdetermined ($n_{substance}>n_{stream}$) or generically undetermined ($n_{substance}<n_{stream}$). In the first case you pick up even more requirements on $W$ for it to have a nontrivial kernel; in the second case I guess you need more flow rate equations.
This dimension reduction is not necessarily beneficial to the numerical analysis; you could very well be better off just attacking the original nonlinear system directly. Indeed with modern computers this is probably the easier approach.
That said, in your particular case most of your equations for the mole fractions just straight up specify what they are, so you actually just have 4 equations in 5 unknowns to solve symbolically, followed by the singularity equation, and then finally you have the 5 equations in 4 unknowns (one of which is redundant) to get the flow rates.
