State-Space Modeling for Process Control Using Mass and Energy Balances.

The purpose of this article is to show how to derive a rigorous state-space model directly from mass and energy balances, then linearize it into A, B, C, D matrices that can be used for simulation, estimation, and control design.

1. What a state-space model means in process engineering.

A state-space model is a compact way to write process dynamics as first-order ordinary differential equations coupled with algebraic output equations.

The general continuous-time nonlinear form is shown below.

dx/dt = f(x, u, p, t). y = g(x, u, p, t).

The vector x contains states that store “process memory,” such as holdup, composition, and temperature, and the vector u contains manipulated inputs, such as flow rates, inlet conditions, and heat duties.

For most chemical processes, the functions f and g come directly from conservation laws combined with constitutive relations, such as reaction kinetics, phase equilibrium, and heat transfer models.

2. Choosing states, inputs, and outputs without ambiguity.

2.1 State selection rules that stay consistent with balances.

A practical state choice follows the energy and mass storage terms that appear in the balances.

Typical state candidates are total holdup, component moles, internal energy or temperature, and in distributed systems, spatially discretized cell inventories and temperatures.

Note : A good state is a quantity whose time derivative appears naturally in a conservation equation, because this avoids hidden algebraic loops and reduces modeling errors.

2.2 Inputs and disturbances are separated for control use.

Manipulated inputs are variables that actuators can change, such as coolant flow, heater power, valve positions, and feed flow rates.

Disturbances are variables that influence the process but are not directly manipulated, such as feed composition changes, ambient temperature changes, and supply pressure variations.

Both can be included in u, or disturbances can be placed in a separate vector d, depending on the controller structure.

2.3 Outputs are what sensors actually measure.

Outputs y should reflect measured variables used for monitoring and feedback control, such as reactor temperature, product composition, and tank level.

If an output is derived from states using an algebraic relation, it should be placed in g(x, u, p, t) so the measurement equation remains explicit.

3. A repeatable workflow for deriving dx/dt from balances.

The derivation is easiest when done in a fixed order that prevents missing terms.

Step. Action. Deliverable.
1. Define control volume, assumptions, and sign conventions. Well-posed modeling scope and consistent directions.
2. Write total mass and component mass balances. Holdup and composition state equations.
3. Write energy balance using enthalpy or internal energy form. Temperature or energy state equation.
4. Add constitutive relations and kinetics. Closed-form f(x, u, p, t).
5. Define outputs from sensors and lab analyzers. Measurement function y = g(x, u, p, t).
6. Compute steady state, then linearize for A, B, C, D. Linear model for controller design.

4. Worked derivation for a nonisothermal CSTR with one reaction.

This example is widely used in process control because it contains coupled mass and energy storage, nonlinear reaction kinetics, and a heat removal mechanism.

4.1 Process description and assumptions.

Consider a perfectly mixed continuous stirred-tank reactor with a single irreversible reaction A → B in liquid phase.

The reactor has constant liquid density ρ and constant heat capacity Cp, and its volume V can be constant or variable depending on the level dynamics included.

A cooling jacket removes heat with a lumped heat transfer term UA(T − Tc), where T is reactor temperature and Tc is coolant temperature.

The feed stream has volumetric flow rate F, inlet concentration CA,in, and inlet temperature Tin.

4.2 Define states, inputs, and parameters.

One common state choice is x = [CA, T]ᵀ for constant volume operation.

A common input choice is u = [F, CA,in, Tin, Tc]ᵀ, where some elements may be treated as disturbances depending on the control structure.

Parameters p can include V, ρ, Cp, UA, reaction pre-exponential k0, activation energy E, gas constant R, heat of reaction ΔH, and any stoichiometric factors.

Symbol. Meaning. Typical unit.
CA.Concentration of A in reactor.mol·m⁻³.
F.Volumetric flow rate.m³·s⁻¹.
V.Reactor liquid volume.m³.
T.Reactor temperature.K.
Tin.Feed temperature.K.
Tc.Coolant temperature.K.
UA.Overall heat transfer coefficient times area.W·K⁻¹.
k(T).Reaction rate constant at temperature T.s⁻¹ for first-order.
ΔH.Heat of reaction for A → B.J·mol⁻¹.

4.3 Mass balance for species A.

The general component balance is accumulation = in − out − reaction consumption.

For a well-mixed constant-volume CSTR, the accumulation term is V dCA/dt.

The inlet term is F CA,in and the outlet term is F CA.

For a first-order reaction rate rA = k(T) CA, consumption is V rA.

The resulting differential equation is shown below.

dCA/dt = (F/V)(CA,in - CA) - k(T) CA.

4.4 Energy balance for reactor temperature.

A common sensible-enthalpy form for a constant-density liquid is accumulation = enthalpy in − enthalpy out + heat of reaction − heat removed.

The accumulation term is ρ Cp V dT/dt.

The flow enthalpy contribution becomes ρ Cp F (Tin − T) under constant Cp and negligible kinetic and potential energy changes.

The heat generated by reaction is (−ΔH) V rA, where −ΔH is positive for exothermic reactions if ΔH is defined negative for exothermicity.

The jacket heat removal is UA (T − Tc), with the sign chosen so that if T is greater than Tc, heat is removed from the reactor.

The resulting temperature equation is shown below.

dT/dt = (F/V)(Tin - T) + ( -ΔH / (ρ Cp) ) k(T) CA - (UA / (ρ Cp V)) (T - Tc).
Note : Always verify the sign convention for ΔH and the heat transfer term, because a sign error can flip stability predictions and mislead controller tuning.

4.5 Close the model with kinetics.

A widely used temperature dependence is the Arrhenius expression.

k(T) = k0 * exp( -E / (R T) ).

Substituting k(T) into the mass and energy balances yields a closed nonlinear state-space model dx/dt = f(x, u, p).

5. Writing the model explicitly as dx/dt = f(x, u).

With x = [CA, T]ᵀ and u = [F, CA,in, Tin, Tc]ᵀ, the state derivative can be written as a vector function.

dx/dt = [ f1(CA, T, F, CA,in) f2(CA, T, F, Tin, Tc) ].

Where f1 and f2 are the right-hand sides of the two balance equations.

Outputs can be chosen as y = [CA, T]ᵀ when both concentration and temperature are measured, or as y = [T] if only temperature feedback is used.

y = g(x, u) = [ CA T ].

6. Deriving the linear state-space matrices A, B, C, D by linearization.

Many control design methods, including LQR, Kalman filtering, and MPC linear prediction models, use a linear time-invariant approximation around a steady operating point.

Let (x̄, ū) be a steady state satisfying f(x̄, ū) = 0.

Define deviations δx = x − x̄ and δu = u − ū.

The first-order Taylor approximation yields the standard linear form.

d(δx)/dt = A δx + B δu. δy = C δx + D δu.

6.1 Jacobian definitions of A, B, C, D.

The matrices are Jacobians evaluated at the operating point.

A = ∂f/∂x |_(x̄,ū). B = ∂f/∂u |_(x̄,ū). C = ∂g/∂x |_(x̄,ū). D = ∂g/∂u |_(x̄,ū).

For the CSTR example, A is a 2×2 matrix where each entry is a partial derivative of f1 or f2 with respect to CA or T at (x̄, ū).

Because k(T) depends on T, derivatives involve both direct terms and chain-rule terms through k(T).

6.2 Derivatives that appear repeatedly in nonisothermal reactors.

For Arrhenius kinetics, the temperature derivative of k(T) is often needed.

dk/dT = k(T) * (E / (R T^2)).

This expression follows directly from differentiating k0 exp(−E/(R T)) with respect to T.

Note : Keep T in Kelvin for Arrhenius expressions, because using Celsius breaks the functional form and invalidates derivatives used in A and B.

6.3 Output matrices for common measurement choices.

If y = [CA, T]ᵀ, then C is the identity matrix and D is a zero matrix of compatible size.

If y = [T], then C = [0, 1] and D remains zero unless the measurement equation explicitly depends on inputs.

7. Extending the approach to variable volume and level dynamics.

If the reactor volume V is not constant, total mass balance introduces an additional state, often liquid height h or volume V.

A typical total mass balance for constant density is dV/dt = Fin − Fout, which becomes a state equation.

Component balances then use V dCA/dt + CA dV/dt = Fin CA,in − Fout CA − reaction, which ensures conservation under changing holdup.

Energy balance similarly includes the effect of changing holdup through the accumulation term ρ Cp V dT/dt plus any additional terms depending on the chosen energy variable.

Note : When volume changes, writing balances directly in extensive form, such as total moles of A and total sensible energy, reduces algebraic mistakes and improves numerical robustness.

8. Discretizing the state-space model for digital control.

Digital controllers and data historians operate in discrete time with sample time Ts.

For a linear continuous-time model, an exact discretization under zero-order hold is defined by matrix exponentials.

x[k+1] = Ad x[k] + Bd u[k]. y[k] = Cd x[k] + Dd u[k].

Where Ad = exp(A Ts) and Bd = ∫₀^{Ts} exp(A τ) dτ B, and Cd = C, Dd = D.

When exact matrix exponentials are not convenient, numerical methods such as forward Euler, backward Euler, or Runge–Kutta can discretize the original nonlinear balances directly.

Note : For stiff reacting systems, explicit Euler can become unstable at practical sample times, so higher-order or implicit integration is often required for reliable simulation and estimation.

9. Practical validation checks that prevent silent modeling errors.

9.1 Unit consistency check.

Every term in each differential equation must share the same unit as the left-hand side.

For dCA/dt, each term must have units of mol·m⁻³·s⁻¹.

For dT/dt, each term must have units of K·s⁻¹.

9.2 Steady-state sanity check.

At steady state, setting dCA/dt = 0 and dT/dt = 0 should produce physically plausible values.

For example, CA should remain nonnegative, and predicted heat removal should balance predicted heat generation plus enthalpy flow effects.

9.3 Energy sign check under limiting cases.

If reaction is turned off by setting k0 = 0, the temperature equation should reduce to a stable mixing and heat-transfer dynamic that approaches a weighted influence of Tin and Tc.

If UA = 0, the only heat exchange should be through the inlet and outlet enthalpy flow, plus reaction heat if present.

10. Implementation templates for simulation and control workflows.

10.1 Balance-based nonlinear simulator skeleton.

The following skeleton shows the structure used in most simulation tools, regardless of programming language.

Given parameters p and initial state x0. For each time step: 1) Read inputs u(t). 2) Compute kinetics k(T). 3) Evaluate f(x, u, p). 4) Integrate ẋ = f to obtain next x. 5) Compute outputs y = g(x, u, p).

10.2 Linearization workflow skeleton for A, B, C, D.

For models defined by balances, linearization is a deterministic sequence.

1) Choose operating point (x̄, ū). 2) Solve f(x̄, ū) = 0 for steady state. 3) Compute Jacobians A = ∂f/∂x and B = ∂f/∂u at (x̄, ū). 4) Define outputs and compute C = ∂g/∂x and D = ∂g/∂u. 5) Validate linear model by comparing small-signal responses to nonlinear simulation.

10.3 Example symbolic differentiation pattern for Jacobians.

When f contains products such as k(T) CA, the chain rule must be applied.

∂(k(T) CA)/∂CA = k(T). ∂(k(T) CA)/∂T = (dk/dT) CA.

11. Common pitfalls and how to avoid them in state-space derivation.

Pitfall. What goes wrong. Prevention.
Mixing intensive and extensive forms inconsistently. Hidden missing terms when holdup varies. Write balances in extensive form first, then convert to intensive variables.
Wrong ΔH sign or heat-transfer direction. Stability and temperature trajectories become unphysical. Test limiting cases with reaction off and with UA off.
Using Celsius in Arrhenius expressions. Rate constants and derivatives become incorrect. Use Kelvin and document the unit convention in parameters.
Ignoring actuator constraints in input definition. Controller design uses infeasible moves. Include saturation limits and rate limits in simulation and MPC formulation.
Linearizing too far from the true operating point. Linear model fails for realistic disturbances. Compute multiple linear models or use gain scheduling when needed.

FAQ

How to choose states when multiple components and reactions exist.

Choose states that represent conserved inventories, such as moles of each independent component or a reduced set based on reaction invariants, plus at least one energy state such as temperature or internal energy.

For N components with R independent reactions, the number of independent composition states is often N−R when expressed in terms of reaction extents and invariants, but a full component-state model is also valid if it remains numerically stable and well-posed.

How to include heat exchanger or jacket dynamics in the state-space model.

Add an additional energy balance for the jacket fluid or wall temperature as another state, and replace the algebraic coolant temperature Tc with a dynamic state such as Tj.

This turns UA(T − Tc) into UA(T − Tj) and adds a differential equation for Tj based on coolant flow and jacket holdup.

How to derive A and B when the model is large.

Use a Jacobian-based approach where each entry is a partial derivative, and compute them analytically when possible or by automatic differentiation when available.

For very large models, sparse Jacobian structures often exist because each balance only depends on a local set of states, which can be exploited for speed and numerical conditioning.

How to handle algebraic constraints and still obtain a state-space form.

If algebraic equations exist, the model is a differential-algebraic system rather than a pure ODE state-space model.

A common approach is index reduction or substitution to eliminate algebraic variables, or using DAE-capable simulators and linearization tools that produce generalized state-space representations.

How to validate that the derived state-space model is correct.

Verify unit consistency, confirm steady-state solutions satisfy the original balances, and compare transient responses against plant data or a trusted high-fidelity simulator for several input steps.

For linear models, validate by applying small step changes around the operating point and matching the linear response to the nonlinear simulation response.

How to connect the derived model to common control methods.

The nonlinear balance model is used for digital twins and nonlinear MPC, and the linearized A, B, C, D model is used for LQR, Kalman filtering, observer design, and linear MPC.

When operation spans multiple regimes, multiple linearizations can support gain scheduling or piecewise-linear MPC.