Multi-component pilot-plant mass balances can be solved in Python by expressing the independent equations as a matrix system, (A\mathbf{x}=\mathbf{b}), and solving for unknown stream flow rates, concentrations, or reaction extents. For a fully specified linear model, NumPy’s numpy.linalg.solve() provides the solution; underdetermined, overdetermined, or nonlinear systems require different methods.

The key is not the Python command itself, but building a physically correct and independent equation set first. Define clear system boundaries, perform a degree-of-freedom check, assemble (A) and (\mathbf{b}) with consistent units, and then select a solver appropriate to the model.

Convert the Pilot Plant into Equations

Define system boundaries

For an entire pilot plant, the overall boundary includes external feeds and final products but excludes internal streams between units.

For a reactor, separator, mixer, or recycle junction, define a separate subsystem boundary and include only streams crossing that boundary. This prevents internal streams from being counted incorrectly in an overall balance.

Write component balances

For each nonreactive component at steady state:

[ \sum \dot{n}{i,\text{in}}-\sum \dot{n}{i,\text{out}}=0 ]

where (\dot{n}_i) is the molar flow rate of component (i).

For a reacting component, include generation or consumption:

[ \sum \dot{n}{i,\text{in}}-\sum \dot{n}{i,\text{out}} +\sum_j \nu_{ij}\xi_j=0 ]

Here, (\nu_{ij}) is the stoichiometric coefficient of species (i) in reaction (j), and (\xi_j) is the extent of reaction.

Choose independent variables

The unknown vector might contain:

Stream flow rates
Component flow rates
Stream concentrations or compositions
Recycle rates
Reaction extents
Unknown utility flows, such as steam-header distributions

For reacting systems, independent reaction extents can reduce the number of unknowns. If (M) independent reactions are identified, species amounts can be written as:

[ N_i=N_{i,0}+\sum_{j=1}^{M}\nu_{ij}\xi_j ]

This avoids introducing redundant reaction variables.

Form the Matrix System in Python

Interpret (A\mathbf{x}=\mathbf{b})

The matrix equation has three parts:

(A): coefficients of the unknown variables
(\mathbf{x}): unknown flow rates, concentrations, or extents
(\mathbf{b}): known feeds, specifications, and constants

For example, suppose two streams are mixed and the output contains 60 units per hour of component A. Let:

[ \mathbf{x}= \begin{bmatrix} F_1\ F_2\ F_3 \end{bmatrix} ]

where (F_1) and (F_2) are feed flow rates and (F_3) is the total outlet flow rate. The equations are:

[ F_1+F_2=100 ]

[ 0.8F_1+0.2F_2=60 ]

[ F_3-F_1-F_2=0 ]

These become:

[ \begin{bmatrix} 1 & 1 & 0\ 0.8 & 0.2 & 0\ -1 & -1 & 1 \end{bmatrix} \begin{bmatrix} F_1\ F_2\ F_3 \end{bmatrix}

\begin{bmatrix} 100\ 60\ 0 \end{bmatrix} ]

Implement the model with NumPy

import numpy as np

A = np.array([
    [1.0, 1.0, 0.0],    # Total feed balance
    [0.8, 0.2, 0.0],    # Component A balance
    [-1.0, -1.0, 1.0]   # Outlet flow specification
])

b = np.array([100.0, 60.0, 0.0])

# Check the coefficient matrix
rank_A = np.linalg.matrix_rank(A)

if rank_A < A.shape[1]:
    raise ValueError("The system is underdetermined or contains dependent variables.")

x = np.linalg.solve(A, b)

F1, F2, F3 = x

print(f"F1 = {F1:.3f}")
print(f"F2 = {F2:.3f}")
print(f"F3 = {F3:.3f}")

The solution is approximately:

F1 = 66.667
F2 = 33.333
F3 = 100.000

The correct NumPy function is np.linalg.matrix_rank(). A function named numpy.linalg.rank() is not the standard NumPy interface.

Perform a Proper Degree-of-Freedom Check

Count unknowns and independent equations

A system has zero degrees of freedom when the number of independent equations matches the number of independent unknowns.

However, simply comparing the number of rows and columns is insufficient. Equations may be duplicated, redundant, or contradictory.

For a square system, check:

n_equations, n_unknowns = A.shape
rank_A = np.linalg.matrix_rank(A)

print("Equations:", n_equations)
print("Unknowns:", n_unknowns)
print("Rank:", rank_A)

For a unique solution, the coefficient matrix must have full column rank. In a square system, this means:

rank_A == A.shape[0] == A.shape[1]

Distinguish underdetermined and overdetermined systems

If there are fewer independent equations than unknowns, the system is underdetermined. Additional specifications or measurements are required.

If there are more equations than unknowns, the system is overdetermined. This commonly occurs when pilot-plant measurements contain redundancy or noise.

Use np.linalg.lstsq() for a least-squares estimate:

x, residuals, rank_A, singular_values = np.linalg.lstsq(
    A, b, rcond=None
)

This estimates the variables that best fit the measured balances rather than forcing every measurement to agree exactly.

Check consistency, not only rank

A system can have full rank and still contain a unit or data error. After solving, calculate the residual:

residual = A @ x - b

print("Residual vector:", residual)
print("Maximum absolute residual:", np.max(np.abs(residual)))

For an exact model, the residual should be close to zero within numerical tolerance.

Extend the Method to Pilot-Plant Networks

Model multiple unit operations

A flowsheet containing mixers, reactors, separators, and recycle loops can be represented by combining the equations from each unit into one global matrix.

For example:

Write component balances around each unit.
Add equations connecting outlet streams to downstream inlet streams.
Add separator recovery or split-fraction specifications.
Add reaction stoichiometry or extent equations.
Add known feed and product specifications.
Assemble all equations into (A\mathbf{x}=\mathbf{b}).

Internal streams remain in the global unknown vector, even though they cancel from an overall plant balance.

Include recycle and purge streams

Recycle systems often produce simultaneous equations because a downstream stream affects an upstream unit through the recycle loop.

A simple linear recycle relation may look like:

[ R=\alpha(F+R) ]

which can be rearranged into:

[ (1-\alpha)R=\alpha F ]

The resulting equation can be placed directly into the matrix. Component-specific recycle fractions produce additional rows for each component.

A purge is modeled similarly by adding the purge flow or purge fraction as an unknown or specification.

Include steam and utility networks

Steam headers at different pressures can be modeled using flow-balance equations at headers and process units.

For each header:

[ \text{steam in}+\text{steam generated}

\text{steam to users}+\text{steam letdown}+\text{steam out} ]

The unknown vector may contain flows through high-, intermediate-, and low-pressure headers, process reboilers, condensate returns, and letdown valves.

The same matrix approach applies, provided pressure-dependent energy relationships are either specified separately or incorporated into a suitable model.

Handle Nonlinear Pilot-Plant Models

Recognize when `np.linalg.solve()` is not sufficient

The direct matrix solver applies only when the equations are linear in the unknowns.

Nonlinear behavior occurs with:

Equilibrium relationships
Nonlinear reaction kinetics
Temperature-dependent properties
Vapor–liquid equilibrium
Pressure-drop correlations
Flow-dependent heat-transfer coefficients
Product terms such as (F x_i) when both variables are unknown

These systems cannot be represented by a constant matrix (A) without linearization.

Use a nonlinear solver

For nonlinear algebraic balances, use SciPy:

import numpy as np
from scipy.optimize import root

def balances(x):
    F1, F2, T = x

    return np.array([
        F1 + F2 - 100.0,
        0.8 * F1 + 0.2 * F2 - 60.0,
        T - (300.0 + 0.1 * F1**2)
    ])

initial_guess = np.array([50.0, 50.0, 300.0])

solution = root(balances, initial_guess)

if not solution.success:
    raise RuntimeError(solution.message)

print(solution.x)
print("Residuals:", balances(solution.x))

The initial guess, variable scaling, and physical bounds can strongly affect nonlinear solver performance.

Validate the Results Against the Physical Process

Check conservation

After solving, independently recompute total and component balances.

For each component:

[ \text{inlet flow} - \text{outlet flow} ]

should equal zero for a nonreactive component, or the expected net reaction contribution for a reactive component.

Check physical feasibility

A mathematically valid solution may still be physically impossible. Check that:

Flow rates are nonnegative.
Mole fractions lie between 0 and 1.
Mole fractions sum to one.
Reaction extents do not consume more reactant than available.
Recovery and split fractions lie within their permitted ranges.
Pressure and temperature values are within equipment limits.

if np.any(x < 0):
    print("Warning: Negative flow rate detected.")

Reconcile experimental measurements

Pilot-plant data rarely close perfectly because of measurement noise, sampling errors, leaks, holdup, and unmeasured losses.

For redundant measurements, formulate a weighted least-squares problem. Weights can reflect measurement uncertainty so that more reliable instruments influence the reconciliation more strongly.

The reconciled result should be compared with raw measurements rather than treated as a replacement for them.

Understanding the Trade-offs

Direct solving is fast but assumes a correct model

np.linalg.solve() is efficient and precise for a well-conditioned, square, linear system.

It cannot correct missing specifications, incorrect stream signs, inconsistent units, or incorrectly selected balances.

Rank checks do not replace engineering judgment

A full-rank matrix confirms algebraic independence, not physical correctness.

A model can be full rank while using the wrong composition, an incorrect stoichiometric coefficient, or a duplicated physical relationship expressed in a misleading way.

Poor conditioning can amplify measurement errors

Even a solvable matrix may be ill-conditioned. Small errors in measured flows or compositions can then produce large changes in calculated variables.

Inspect the condition number:

condition_number = np.linalg.cond(A)
print("Condition number:", condition_number)

A very large condition number indicates that scaling, better measurements, reformulation, or regularization may be needed.

Unit inconsistency silently corrupts results

Do not mix kg/h, kmol/h, mol/s, lb/h, or volumetric units without explicit conversion.

Use one consistent unit system for each balance, and apply molecular weights when converting between mass and molar flow.

How to Apply This to Your Pilot Plant

A reliable workflow is:

If your primary focus is a fully specified linear flowsheet: Define independent balances, assemble A and b, verify rank and residuals, and solve with np.linalg.solve().
If your primary focus is experimental data reconciliation: Include redundant measurements and use np.linalg.lstsq() with appropriate weighting and residual analysis.
If your primary focus is reaction systems: Identify independent reactions and use reaction extents with stoichiometric coefficients to avoid redundant variables.
If your primary focus is recycle, equilibrium, or kinetic behavior: Use the matrix method for the linear portions and scipy.optimize for nonlinear equations.
If your primary focus is process validation: Check component conservation, nonnegative flows, composition limits, units, and agreement with measured pilot-plant data.

With disciplined equation setup and physical validation, Python turns complex multi-component pilot-plant balances into transparent, testable, and reusable process models.

Summary Table:

Modeling Situation	Recommended Python Method	Essential Checks
Fully specified linear balances	`numpy.linalg.solve()`	Confirm full rank, consistent units, and near-zero residuals
Underdetermined system	Add independent measurements or specifications	Verify degrees of freedom
Overdetermined or noisy measurements	`numpy.linalg.lstsq()`	Review residuals and apply measurement weighting
Recycle and purge networks	Global matrix model	Include internal streams and recycle relationships
Reactions	Stoichiometric coefficients and independent reaction extents	Check reactant availability and reaction consistency
Equilibrium, kinetics, or flow-dependent properties	`scipy.optimize.root()` or another nonlinear solver	Use good initial guesses, scaling, and physical bounds
Experimental pilot-plant data	Weighted least-squares reconciliation	Compare reconciled values with raw measurements
Any process model	Independent conservation and feasibility checks	Check nonnegative flows, compositions, units, and equipment limits

Need a reliable way to model your pilot plant? Contact LABPARK to discuss Educational and Vocational Unit Operations Pilot Plants for chemical engineering, bioprocess and biotech, and environmental and water treatment applications. LABPARK helps universities, research institutes, and enterprises build practical training and research platforms with measurable material balances, configurable process units, recycle and utility networks, and data suitable for Python-based modeling, validation, and process optimization. Tell us your target process, teaching or research objectives, and required scale so we can recommend a fit-for-purpose pilot-plant solution.