luis logo

Control Theory for Reactor Automation


Luis López Reviriego
11-12-2025

1. Introduction

The Laplace Transform is a mathematical tool that allows linear differential equations to be converted into simpler algebraic problems, with the aim of computing the solutions of the differential equations from the solution of the corresponding algebraic problem (Cánovas, 2025).

Let f{f} be an integrable function, that is, the Riemann integral of f{f} exists over every compact interval. The Laplace Transform of the function is defined as:

L{f(t)}=F(s)=0estf(t)dt\mathcal{L}\{f(t)\} = F(s) = \int_0^\infty {e^{-st}}f(t)dt

Fortunately, the integral is a tabulated integral, that is, it has been solved for most commonly used functions, with the results collected in tables, which significantly simplifies the process.

In summary, in order to solve the linear equations that model the dynamical system, one can:

  1. Transform the differential equation into the Laplace domain.
  2. Solve the resulting algebraic expression in the Laplace domain, where it is simpler to handle.
  3. Use the table of Laplace Transforms to perform the inverse transformation and return the equation to the time domain.

The first step would be to define the differential equations that model the system. According to Seborg (2010), there are empirical, theoretical, and semi-empirical models. Empirical models are obtained using experimental data from the system and are easier to develop, but they are not robust over a wide range of operating conditions. On the other hand, purely theoretical models are highly complex. Semi-empirical models are the most widely used in industry.

In conclusion, a semi-empirical model would be preferred, as it combines the simplicity of implementation of empirical models with the robustness of theoretical models under different operating conditions.

The Laplace Transform is used to obtain closed-form solutions to linear differential equations. In general, chemical processes are highly nonlinear and unstable. It is often very difficult, or even impossible, to solve these equations directly. Even when the system parameters are known, powerful computational tools are required to obtain solutions. To simplify the task of description and modeling, the system is frequently linearized and matrix algebra is applied to solve the linearized equations (Woolf, 2025).

For example, the reaction rate of first-order, gas producing chemical reactions (AB+gas{A \rightarrow B + gas}) that occur in certain types of reactors depends on a kinetic constant that is highly temperature-dependent. This dependence is given by the Arrhenius equation (Seborg, 2010):

r=k0eEATr = k_{0}e^{-\frac{E_A}{T}}

Where k0{k_{0}} is the Arrhenius constant, EA{E_A}​ is the activation energy, and T{T} is the temperature.

This relationship is nonlinear. Furthermore, the reaction can be exothermic or endothermic (releasing or consuming heat, respectively), directly affecting the temperature of the tank. On the other hand, the pressure may refer to the pressure of a gas in a gaseous phase above the liquid in the reactor, which is regulated by a gas valve acting as an actuator. In this case, if there is a large pressure difference between the inside and outside of the reactor, the gas flow will depend on the pressure gradient across the valve and will be subject to compressibility effects; that is, it will not be linear.

All of these conditions make the mathematical models of chemical processes highly nonlinear. For this reason, in order to apply the Laplace Transform in the analysis and design of the control system, these models must first be linearized around an operating or equilibrium point.

Equilibrium points can be stable or unstable. If the system starts near equilibrium, it may move toward the equilibrium or away from it (North Carolina State University, Department of Mathematics, 2017).

2. Linearization Around the Operating Point

According to Woolf (2025), the first step consists of performing a Taylor series expansion:

f(x)=n=0fn(a)n!(xa)nf(x) = \sum_{n=0}^\infty \frac{f^{n}(a)}{n!}(x-a)^n f(x)=f(a)+f(a)(xa)+f2(a)2!(xa)2++fn(a)n!(xa)nf(x) = f(a) + f'(a)(x-a) + \frac{f^{2}(a)}{2!}(x-a)^2+\dots +\frac{f^n(a)}{n!}(x-a)^n

Where a{a} is the operating point around which the linearization will be performed, and (xa){(x-a)} is the deviation of the point x{x} from its equilibrium value. Since most processes are intended to operate near a steady state, this point corresponds to the steady-state or equilibrium point. Only the first two terms of the expansion (zeroth and first order) are used, which ensures a linear polynomial. As the zeroth-order term, f(a){f(a)}, corresponds to the function at the equilibrium point, its value is 0 (f(a)=0{f(a) = 0}).

For the practical case, assuming a CSTR (Continuous Stirred-Tank Reactor) with constant volume, two internal variables (pressure and temperature), two actuators (uP{u_P}​ and uT{u_T}, for pressure and temperature, respectively), and ambient temperature as a disturbance (dT{d_T}​), it is possible to express its differential equations in the following form, using Taylor series expansions, up to first order:

{dPdt=fPPΔP+fPTΔT+fPuPΔuP+fPuTΔuT+fPdTΔTambdTdt=fTPΔP+fTTΔT+fTuPΔuP+fTuTΔuT+fTdTΔTamb \left\{\begin{aligned} \frac{dP}{dt} = \frac{\partial f_P}{\partial P}\Delta P + \frac{\partial f_P}{\partial T}\Delta T + \frac{\partial f_P}{\partial u_P}\Delta u_P + \frac{\partial f_P}{\partial u_T}\Delta u_T +\frac{\partial f_P}{\partial d_T}\Delta T_{amb} \\ \frac{dT}{dt} = \frac{\partial f_T}{\partial P}\Delta P + \frac{\partial f_T}{\partial T}\Delta T + \frac{\partial f_T}{\partial u_P}\Delta u_P + \frac{\partial f_T}{\partial u_T}\Delta u_T + \frac{\partial f_T}{\partial d_T}\Delta T_{amb} \end{aligned}\right.

Where:

ΔP=PPsetΔT=TTsetΔuP=uPuP,setΔuT=uTuT,set\begin{aligned} &\Delta P = P - P_{set} \\ &\Delta T = T - T_{set} \\ &\Delta u_P = u_P - u_{P, set} \\ &\Delta u_T = u_T - u_{T, set} \end{aligned}

All variables with the “set” subscript correspond to values that place the system at equilibrium (which we fix at the setpoint).

When working with multiple variables and, consequently, multiple differential equations, Jacobians in matrix form are used. These matrices contain the constants required to describe the linearity of the system of differential equations.

J(x1,...,xn)=[y1x1y1xnynx1ynxn]J(x_1,...,x_n) = \begin{bmatrix} \frac{\partial y_1}{\partial x_1} & \dots & \frac{\partial y_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial y_n}{\partial x_1} & \dots & \frac{\partial y_n}{\partial x_n} \\ \end{bmatrix}

Thus, the linearized system, for two variables (pressure and temperature), can be expressed as:

[P˙T˙]=[fPPfPTfTPfTT]J[PPsetTTset]\begin{bmatrix} \dot{P} \\ \dot{T} \end{bmatrix} = \underbrace{ \begin{bmatrix} \frac{\partial f_P}{\partial P} & \frac{\partial f_P}{\partial T} \\ \frac{\partial f_T}{\partial P} & \frac{\partial f_T}{\partial T} \end{bmatrix} }_{J} \begin{bmatrix} P-P_{set} \\ T - T_{set} \end{bmatrix}

This corresponds to a system without actuators or disturbances, that is, a system that oscillates around the operating point on its own and may be either stable or unstable. Actuators allow the system to be maintained near the operating point.

When actuators and disturbances are included, their effects are added to the linear model. Assuming a pressure actuator (uP{u_P}), a temperature actuator (uT{u_T}), and ambient temperature as a disturbance (dT{d_T}):

[P˙T˙]=[fPPfPTfTPfTT][PPsetTTset]+[fPuPfPuTfTuPfTuT][upuP,setuTuT,set]+[fPdTfTdT]ΔTamb\begin{bmatrix} \dot{P} \\ \dot{T} \end{bmatrix} = \begin{bmatrix} \frac{\partial f_P}{\partial P} & \frac{\partial f_P}{\partial T} \\ \frac{\partial f_T}{\partial P} & \frac{\partial f_T}{\partial T} \end{bmatrix}\begin{bmatrix} P-P_{set} \\ T - T_{set} \end{bmatrix} + \begin{bmatrix} \frac{\partial f_P}{\partial u_P} & \frac{\partial f_P}{\partial u_T} \\ \frac{\partial f_T}{\partial u_P} & \frac{\partial f_T}{\partial u_T} \end{bmatrix}\begin{bmatrix} u_p - u_{P,set}\\ u_T - u_{T,set} \end{bmatrix} +\begin{bmatrix} \frac{\partial f_P}{\partial d_T} \\ \frac{\partial f_T}{\partial d_T} \end{bmatrix} \Delta T_{amb}

In short:

x˙=Aδx+Bδu+EδTamb\dot{x} = A\delta x + B\delta u + E\delta T_{amb}

3. Laplace Transform and Transfer Functions

With the linearized model around the operating point, we proceed to transform it into the Laplace domain using the Laplace Transform. The Laplace form of a time derivative is:

L{f(t)}=sF(s)f(0)\mathcal{L}\{f'(t)\} = sF(s) - f(0)

In this case, the function at time 0 corresponds to the equilibrium, so f(0)=0{f(0) = 0} (we start at equilibrium with no initial deviations).

L{Δx˙}=sΔX(s)\mathcal{L}\{\Delta\dot{x}\} = s\Delta X(s)

Applied to the full state equation, and using the variable-separation property of the Laplace Transform:

L{Δx˙}=AL{Δx}+BL{Δu}+EL{Δd}\mathcal{L}\{\Delta\dot{x}\} = A\mathcal{L}\{\Delta x\}+B\mathcal{L}\{\Delta u\} + E\mathcal{L}\{\Delta d\}

Which results in:

sΔX(s)=AΔX(s)+BΔU(s)+EΔD(s)s\Delta X(s) = A\Delta X(s) + B\Delta U(s) + E\Delta D(s)

Once the linearized model has been transformed, it is possible to obtain the transfer functions:

ΔX(s)=(sIA)1(BΔU(s)+EΔD(s))=Gu(s)ΔU(s)+Gd(s)ΔD(s)\begin{aligned} \Delta X(s)&=(sI-A)^{-1}(B\Delta U(s) + E \Delta D(s)) \\\\ &= G_u(s)\Delta U(s) + G_d(s)\Delta D(s) \end{aligned}

The control transfer function, which relates deviations of the actuators to deviations in the internal states of the reactor, is:

Gu(s)=(sIA)1BG_u(s) = (sI-A)^{-1}B

Where:

  • (sIA)1{(sI-A)^{-1}} captures the internal dynamics of the system.
  • B{B} includes the effect of the actuators.

Similarly, the disturbance transfer function is:

Gd(s)=(sIA)1EG_d(s) = (sI-A)^{-1}E

Where:

  • (sIA)1{(sI-A)^{-1}} captures the internal dynamics of the system.
  • E{E} includes the effect of disturbances.

Each element of the control transfer function matrix gives the relationship between a deviation in an actuator and a deviation in a reactor internal state variable in a simple algebraic form. For example, the relationship between deviations in the temperature actuator and deviations in reactor temperature is:

GT,uT=ΔT(s)ΔuT(s)G_{T,u_T} = \frac{\Delta T(s)}{\Delta u_T(s)}

This corresponds to open-loop transfer functions, since feedback on the inputs has not been considered. Considering a two-block closed-loop system, the following functions in the Laplace domain are defined:

  • A(s){A(s)}: reference input to the controller (setpoint values).
  • U(s){U(s)}: actuator actions.
  • P(s){P(s)}: PID controller.
  • D(s){D(s)}: disturbances.
  • V(s){V(s)}: controller output.
  • Gu{G_u}(s): actuator dynamics (open-loop model).
  • Gd{G_d}(s): dynamics of disturbances affecting the internal states.
  • X(s){X(s)}: feedback signal.
  • H(s){H(s)}: feedback block, modeling the effect of sensors on the output readings (signal gain or attenuation). If sensors do not alter the signal, its value is 1.
  • R(s){R(s)}: error signal, obtained as the difference between the reference A(s){A(s)} and the feedback signal X(s){X(s)}.

The following relations are obtained:

R(s)=A(s)X(s)X(s)=H(s)V(s)U(s)=P(s)R(s)V(s)=Gu(s)U(s)+Gd(s)D(s)\begin{aligned} &R(s) = A(s) - X(s) \\\\ &X(s) = H(s)V(s) \\\\ &U(s) = P(s)R(s) \\\\ &V(s) = G_u(s)U(s) + G_d(s)D(s) \end{aligned}

We want to obtain the final expression for the outputs in our model. Substituting into the equation for V{V}:

V(s)=Gu(s)[P(s)[A(s)H(s)V(s)]]+Gd(s)D(s)V(s)=Gu(s)P(s)A(s)Gu(s)P(s)H(s)V(s)+Gd(s)D(s)[1+Gu(s)P(s)H(s)]V(s)=Gu(s)P(s)A(s)+Gd(s)D(s)V(s)=Gu(s)P(s)1+Gu(s)P(s)H(s)A(s)+Gd(s)1+Gu(s)P(s)H(s)D(s)\begin{aligned} &V(s) = G_u(s)[P(s)[A(s) - H(s)V(s)]] + G_d(s)D(s) \\\\ &\Rightarrow V(s) = G_u(s)P(s)A(s)-G_u(s)P(s)H(s)V(s) + G_d(s)D(s) \\\\ &\Rightarrow [1+G_u(s)P(s)H(s)]V(s) = G_u(s)P(s)A(s) + Gd(s)D(s) \\\\ &\boxed{V(s) = \frac{G_u(s)P(s)}{1+G_u(s)P(s)H(s)}A(s)+ \frac{G_d(s)}{1+G_u(s)P(s)H(s)}D(s)} \end{aligned}

The effect of D(s){D(s)} on V(s){V(s)} is attenuated by the feedback term (1+Gu(s)P(s)H(s){1+G_u(s)P(s)H(s)}), reflecting the closed-loop system’s ability to compensate for external disturbances.

The overall transfer function for the closed-loop model, including disturbances, is:

M(s)=V(s)A(s)=Gu(s)P(s)1+Gu(s)P(s)H(s)Inputs’ effect+Gd(s)A(s)[1+Gu(s)P(s)H(s)]D(s)Disturbances’ effect\boxed{M(s) = \frac{V(s)}{A(s)} = \underbrace{\frac{G_u(s)P(s)}{1+G_u(s)P(s)H(s)}}_{\text{Inputs' effect}} + \underbrace{\frac{G_d(s)}{A(s)[1+G_u(s)P(s)H(s)]}D(s)}_{\text{Disturbances' effect}}}

This equation describes the total transfer function of the closed-loop system, combining the effect of the reference input A(s){A(s)} and the effect of disturbances D(s){D(s)}, both attenuated by the closed-loop denominator 1+Gu(s)P(s)H(s){1+G_u(s)P(s)H(s)}, which represents the feedback action. It constitutes the closed-loop model of the linearized system, including the effects of external disturbances.

4. PID Controller

With the linearized differential model, it is possible to analyze the system stability as well as its dynamic behavior around the operating point (setpoint), considering small deviations from equilibrium (Woolf, 2025).

By using the Laplace Transform and obtaining the transfer functions, PID controllers can be designed employing frequency-domain control techniques. This approach allows improving the accuracy, robustness, and stability of the control system in industrial applications (Caman Cruz & Sánchez Márquez, 2024).

PID controllers are used in more than 90% of automated control systems in industrial processes. Nowadays, PID controllers are implemented in PLC (Programmable Logic Controller) and DCS (Distributed Control System) software (Hägglund & Guzmán, 2024).

The mathematical formula of the PID controller is as follows:

u(t)=Ke(t)+1Ti0te(τ)dτ+Tdde(t)dtu(t) = Ke(t) + \frac{1}{T_i} \int_0^t{e(\tau)d\tau} + T_d \frac{de(t)}{dt}

In Laplace's Domain:

C(s)=Kp+Kis+KdsC(s) = K_p+\frac{K_i}{s}+K_ds

Where:

  • Kp{K_p} is the proportional gain.
  • Ki{K_i} is the integral gain.
  • Kd{K_d} is the derivative gain.

To determine the optimal values of the PID parameters, tuning techniques such as Ziegler-Nichols are employed, which derive the parameters from the study of the system’s response to a step input (Caman Cruz & Sánchez Márquez, 2024).

When designing the controller, it is also important to consider the effect of actuators on the output signals. For instance, a servomotor acts as an integrator, so it may be necessary to differentiate the output to compensate for the integration introduced by the actuator. This compensation must be implemented directly in the controller’s code (Bhandari & Csurcsia, 2022).

It is equally important to account for sensor measurement noise. In industrial applications, noise typically occurs at high frequencies, while the useful information is found at lower frequencies. Therefore, it is advisable to incorporate a low-pass filter to reduce the noise reaching the controller. Moreover, the derivative component of the PID controller is sensitive to this high-frequency noise, as well as to abrupt changes in the setpoint, which cause rapid variations in the error and significantly increase its derivative.

Regarding the integral component of the controller, measures should be taken to prevent actuator saturation. The integrator accumulates error over time, so if the actuator reaches its physical limits (e.g., a fully open valve), the error will continue to accumulate without affecting the controlled variable. When a subsequent change in direction is required, the controller must compensate for this accumulation, generating a delay known as “integrator windup.” An effective strategy to mitigate this phenomenon is to detect actuator saturation and temporarily limit the integral action until the system returns to an acceptable operating range.

5. References