Skip to Main Content
Skip Nav Destination

Physical models of reacting systems aim to generalise experimental observations and mechanistic process understanding solely through the application of chemical kinetics principles and without the incorporation of any data-driven approaches. At this chapter’s core are the fundamental concepts of the mass-action law and pseudo-steady-state hypothesis, the application of which is exemplified through the derivation of both chemical and biochemical kinetic models. In virtue of its indispensable role during model development, deterministic optimisation theory is also introduced to illustrate the formulation and solution of constrained optimisation problems, with particular focus on nonlinear parameter estimation problems. This lays the theoretical foundation for model construction methodologies explored in latter chapters.

The construction of an accurate and reliable kinetic model is the foremost task in (bio)chemical reactor design. From an analytical standpoint, the model and its kinetic parameters provide insight on the process by describing the concentration and temperature dependences that control the rate of reaction, potentially allowing us to identify mechanistic information such as inhibiting effects and rate-limiting steps in the underlying reaction network. From a practical application standpoint, the kinetic model provides a mathematical representation of the reaction dynamics, thus linking the reactor performance (e.g. conversion or selectivity) with the operating conditions and feed compositions; this, in turn, makes it an essential tool for reactor optimisation and control.

The complexity and degree of detail of a kinetic model obviously depend on its desired application. In this chapter, we exemplify the physical model construction procedure (from both empirical and mechanistic perspectives) for simple chemical and biochemical reactions. This will serve as an introduction for more complex model development methodologies explored in other chapters.

Empirical kinetic models refer to rate expressions that have been developed through experimental observations in the absence of deep mechanistic understanding of the reaction process. The most common example of one such model is the power law approximation, whereby the rate expression is dependent on the concentrations of reactants and products raised to some apparent reaction order coefficients, and the rate temperature dependence is expressed in terms of the Arrhenius equation. For example, for some trivial reaction A → B:
(1.1)
where r is the rate of reaction (mol m−3 s−1), Eapp is the apparent activation energy of the reaction (J mol−1), Rg is the universal gas constant (8.314 J mol−1 K−1), T is temperature (K), CA and CB are the concentrations (mol m−3) of species A and B, respectively, and α and β are their apparent reaction orders. The pre-exponential factor k0 is temperature independent and its units depend on the reaction orders.

It must be noted that empirical models are not necessarily void of physical interpretation. Indeed, if the reaction in consideration were an elementary reaction, the power law expression in eqn (1.1) would simply reduce to a mass-action law model where the reaction order of the reactant corresponds to its stoichiometric coefficient and the apparent activation energy corresponds to the real activation energy of the reaction. More often than not, however, power laws are limiting approximations of the true, complex underlying kinetics that govern the reaction. For instance, in transport-limited heterogeneous catalytic reactions the apparent characteristics used to approximate the rate equation can be shown to depend not only on the intrinsic chemical kinetics of the reaction but also on the transport parameters of the system.1  Indeed, the apparent kinetic characteristics of a reaction may vary significantly across different operating ranges (i.e. concentrations, temperatures, and flow regimes, among others), meaning that empirical kinetic models typically exhibit poor extrapolability beyond the conditions on which they have been trained. Moreover, their model structures and parameters do not capture the mechanistic details of the reaction, and so their applicability for tasks such as model-based design of experiments, catalyst development or reactor design is limited. Nonetheless, empirical rate laws are extensively used in industrial process monitoring and control owing to their simplicity, low computational cost, and ease of implementation alongside heat transfer and catalyst deactivation models.2 

Although the main focus of this chapter is the development of more rigorous kinetic models, the dynamic parameter estimation techniques described in Section 1.3 are equally applicable to empirical models.

Mechanistic kinetic models aim to capture the concentration and temperature dependences of the rate by considering the sequence of elementary steps that constitute the overall reaction. In this section, we introduce the fundamental concepts of the mass-action law and the pseudo-steady-state hypothesis. Rate expressions for some common reaction mechanisms are also developed, thus setting the basis for the construction of microkinetic models for complex heterogeneous reactions in Chapter 8.

Reaction mechanisms are composed of several elementary reactions. An elementary reaction refers to a chemical transformation that occurs in one step where the reactants directly form the product, passing through a single transition state without the involvement of any intermediates. The rate expressions for such elementary reactions are generalised by the well-known mass-action law, which states that the rate is proportional to the concentration of the reactants raised to their stoichiometric coefficients. For example, for an irreversible elementary reaction as shown in eqn (1.2), the rate of reaction is given by eqn (1.3):
(1.2)
(1.3)
The rate constant k is a temperature dependent parameter which reflects the fact that the reaction rate increases exponentially with temperature. This is typically assumed to be an Arrhenius dependence in the form:
(1.4)
where A is the temperature-independent pre-exponential factor (same units as k), Ea is the activation energy of the elementary reaction (J mol−1), Rg is the universal gas constant (8.314 J mol−1 K−1), and T is temperature (K).
For homogeneous reactions, the rate r (mol m−3 s−1) is defined as the number of chemical transformations taking place per unit time per unit volume of the reaction mixture; note that the molar rate of disappearance or formation of any of the species involved can then be obtained by multiplying the rate of reaction by their corresponding stoichiometric coefficient (by convention, these are negative for the reactants and positive for the products):
(1.5)
(1.6)
These rate expressions are sufficient to simulate the reactant and product concentration–time profiles in simple batch systems. For more complex reactor configurations, we must couple the kinetic rate laws alongside material balances to obtain a reactor model. For example, for the previous elementary reaction in an ideal steady-state plug flow unit, we would obtain the rate of change of species concentration along the space-time of the reactor (assuming constant density of the reaction mixture):
(1.7)
where τ(s) is the space-time, defined as the amount of time required for the reaction mixture to pass through a certain length of the reactor’s spatial dimension. Namely, at the reactor entrance τ = 0 and at the reactor exit τ = Vr/qv where Vr is the reactor volume and qv is the feed volumetric flowrate.
For heterogeneous reactions, it is more appropriate to express the rate of elementary steps in terms of the fractional coverages of the surface species involved. The coverage of a surface species simply refers to its concentration normalised against the total concentration of (both occupied and unoccupied) active sites in the catalyst. Consider, for instance, the following irreversible, heterogeneous elementary reaction:
(1.8)
Here Ai and Bi refer to the fluid-phase reactants and products, while Xj and Yj are the surface reactants and products, respectively. Then, for the rate of the irreversible elementary reaction, the mass-action law can be formulated in terms of the gas-phase concentrations and the coverage of surface species:
(1.9)
where
(1.10)
where θXj is the coverage of species Xj corresponding to a surface concentration SXj (mol g−1) and SθT is the total number of surface species per gram of catalyst (mol g−1). For heterogeneous catalytic reactions, the rate rw (mol g−1 s−1) is typically defined as the number of chemical transformations taking place per unit time per unit mass of catalyst. Thus, if this elementary reaction were carried out in an ideal packed-bed tubular unit (and for constant reaction mixture density), we would obtain the following equation for the rate of disappearance of the fluid-phase species:
(1.11)
Here the spatial dimension of interest is W (g), the catalyst mass over which the reaction mixture is contacted. That is, at the reactor entrance W = 0, whereas at the exit W = Vr · ρb where ρb is the bulk density of the catalyst (g m−3).

Now that we have elucidated the application of the mass-action law to the development of rate and (ideal) reactor equations for single elementary reactions, let us consider some reaction mechanisms involving multiple elementary steps. Of particular interest are the series and parallel reaction mechanisms shown in Figure 1.1; albeit these simple schemes are rarely found in practice, they are illustrative examples given that complex reaction mechanisms combine both parallel and series elementary steps.

Figure 1.1

Series and parallel reaction mechanisms. ki refers to the reaction constant for each elementary step i.

Figure 1.1

Series and parallel reaction mechanisms. ki refers to the reaction constant for each elementary step i.

Close modal

The mass-action law can once again be applied to develop rate equations for the elementary steps in the mechanisms, after which the net rate of consumption (or formation) of each species can be worked out from a simple material balance. These are summarised in Table 1.1.

Table 1.1

Kinetic models for simple mechanisms.

Reaction scheme Rate equations Species molar balance
Series  r1 = k1CA   dCAdt=r1  
r2 = k2CB   dCBdt=r1r2  
  dCCdt=r2  
Parallel  r1 = k1CA   dCAdt=(r1+r2)  
r2 = k2CA   dCBdt=r1  
  dCCdt=r2  
Series–parallel  r1 = k1CA   dCAdt=r1  
r2 = k2CB   dCBdt=r1(r2+r3)  
r3 = k3CB   dCCdt=r2  
  dCDdt=r3  
Reaction scheme Rate equations Species molar balance
Series  r1 = k1CA   dCAdt=r1  
r2 = k2CB   dCBdt=r1r2  
  dCCdt=r2  
Parallel  r1 = k1CA   dCAdt=(r1+r2)  
r2 = k2CA   dCBdt=r1  
  dCCdt=r2  
Series–parallel  r1 = k1CA   dCAdt=r1  
r2 = k2CB   dCBdt=r1(r2+r3)  
r3 = k3CB   dCCdt=r2  
  dCDdt=r3  
It should be noted that in many instances the concentration of some of the intermediate species is not observable. This is particularly true for heterogeneous catalytic reactions, where the fractional coverage of active site surface species cannot be measured in practical applications. In such cases, the pseudo-steady-state hypothesis (PSSH) is applied; this allows for the elimination of non-observable states from the kinetic model, albeit at the cost of incorporating the assumption that the net rate of formation of intermediates is negligible. Consider, for instance, the series–parallel reaction scheme in Figure 1.1. If the concentration of reaction intermediate B is not measurable, then we apply the PSSH such that
(1.12)
This, in turn, allows for the concentration of the intermediate CB to be reformulated in terms of the measurable state CA and the rate constants:
(1.13)
This result can then be substituted into the expression for the net rate of product formation to obtain a kinetic model that only depends on the concentration of observable species:
(1.14)
(1.15)
Interestingly, we notice that the numerators of eqn (1.14) and (1.15) correspond to the product of the rate constants of elementary steps involved in the overall reactions A → C and A → D, respectively. On the other hand, the denominator of the expressions is the sum of the rate constants of the elementary reactions B → C and B → D, thus reflecting the mutual inhibiting effect (or kinetic resistance) that arises due to the “competition” between these two elementary steps. For more complex mechanisms involving multiple intermediates and reversible elementary reactions, the application of the PSSH yields kinetic model structures akin to the well-known Langmuir kinetic model; these can be similarly interpreted as the ratio between a reaction driving force and a kinetic resistance. This is further discussed in Chapter 8.

Our discussion so far has only concerned intrinsic kinetics whereby the reaction chemistry is the rate-limiting factor. In reality, both homogeneous and heterogeneous reactions encompass other physical phenomena (such as heat and mass transport) that may kinetically control the reaction process. For example, solid-catalysed reactions may be internal mass transfer limited whereby the diffusion of species within the catalyst pores is rate-limiting, or external mass transfer limited whereby interfacial transport between the catalyst surface and the bulk fluid is rate-limiting. In either case, there is a mismatch between the observed reaction rate and the intrinsic rate that would be expected in the absence of mass-transfer limitations.

To illustrate this, consider a first-order catalytic reaction that is external mass transfer limited. For simplicity, we assume that the rate of diffusion from the bulk fluid to the catalyst surface is given by the linear driving force (LDF) model:3 
(1.16)
where rd is the molar rate of diffusion per unit mass of catalyst (mol g−1 s−1), ke (m s−1) is the LDF mass transfer coefficient, as (m2 g−1) is the specific surface area of the catalyst, Cb (mol m−3) is the concentration of the reactant at the fluid bulk, and Cs (mol m−3) is the concentration of the fluid in equilibrium with the catalyst surface.
At steady state, the rate of diffusion of the reactant must be equal to its rate of disappearance due to chemical transformations at the catalyst. Therefore, for a first-order elementary reaction we have
(1.17)
where kw (m3 g−1 s−1) is the intrinsic rate constant for the first-order elementary reaction.
Rearranging eqn (1.17), we can obtain an expression for the concentration of the reactant at the catalyst surface in terms of its concentration at the fluid bulk:
(1.18)
This result can then be substituted into the intrinsic kinetics rate equation to yield the reaction rate observed in the external mass transfer limited case:
(1.19)
Now, we introduce the concept of the external effectiveness factor. This is simply defined as the ratio between the external mass transfer limited reaction rate in eqn (1.19) and the intrinsic rate that would correspond to bulk conditions (i.e. in the absence of transport limitations Cs = Cb):
(1.20)
The term kwkeas is known as the second Damköhler number, a dimensionless ratio between the intrinsic chemical reaction and external mass transport timescales. Namely, we notice that if the intrinsic reaction is significantly slower than external mass transport, then kw ≪ keas and ηex ≈ 1, such that the observed reaction rate is kinetically controlled. Conversely, if external mass transport is slower than the surface chemistry kinetics, we then have kw ≫ keas and ηex ≪ 1, such that the observed rate is dominated by the diffusion to the catalyst. Similar conclusions can be obtained for other heterogeneous systems such as gas–liquid reactions.

In the field of bioprocess kinetics, most modeling methodologies can be classified into structured and unstructured approaches. Structured models are detailed mathematical representations of the metabolism of microorganisms, including information such as microbial structure and physiology to describe cell cultures whose morphology and composition are strongly time-dependent.4  These are particularly useful for metabolic engineering and cellular regulatory process studies.5  On the other hand, unstructured approaches assume the cell culture to be a homogeneous biomass and the dynamics of cell growth, substrate consumption and extracellular product formation are all treated from a macroscopic perspective. These approaches therefore provide a compromise between the mathematical complexity of the model and the incorporation of fundamental microbial kinetics knowledge, thus serving as a powerful tool for the control and optimisation of bioprocesses.6 

In this section of the chapter, we will discuss some common unstructured kinetic models and show how these can be mechanistically derived by approximating the cellular metabolism as an enzymatic reaction mechanism. Given the semi-empirical nature of unstructured models, the reader should note that many theoretical derivations have been proposed in the past by incorporating different assumptions related to kinetic, thermodynamic and/or substrate transport characteristics of the biochemical system.7  As such, some of the kinetic parameters that arise in unstructured models do not have a universal physical meaning but rather depend on the assumptions made during model construction.

The natural starting point for the construction of biochemical kinetic models is to consider biomass growth and substrate consumption, as in the absence of product inhibition effects these should be a function of substrate and biomass states only. The growth of a microorganism on a single substrate can, for instance, be approximated as a simple enzymatic reaction in the form:
(1.21)
where S is the substrate promoting growth and Ysx is its corresponding mass stoichiometric coefficient. XA and XS refer to active and saturated biomass states, respectively, which represent the cellular metabolism through which the substrate is assimilated for growth. Needless to say, this approximation is not free from criticism as we are assuming that the metabolism of biomass can be classified into two extremes, namely, saturated and unsaturated biomass. Furthermore, the careful reader may have noticed that this is analogous to an irreversible enzymatic mechanism where the product of the reaction is the participating enzyme itself. This is clearly a generalisation, as in reality microbial growth is controlled by a large number of rate-limiting metabolic pathways rather than two irreversible elementary steps. Following the same procedure shown in Section 1.1.2, we can apply the mass-action law to develop the rate expressions in eqn (1.22) and (1.23) for the biomass species; as bioprocess state measurements are in units of mass per volume, these have been expressed accordingly:
(1.22)
(1.23)
where XA, XS and S are the concentrations (g L−1) of active biomass, saturated biomass and substrate, respectively. k1 (L g−1 h−1) is the reaction constant for substrate assimilation, while k2 (h−1) is the product desorption rate constant. It is obvious that the active and saturated biomass states cannot be observed separately, rather the total concentration of viable biomass X (g L−1) is measured:
(1.24)
(1.25)
Together with eqn (1.24), the pseudo-steady-state hypothesis can be applied to the biomass–substrate intermediate (setting eqn (1.23) to zero) to obtain expressions for XS and XA that are only in terms of the observable states X and S:
(1.26)
The expression for the saturated biomass intermediate can be substituted into the rate equation for total biomass (eqn (1.25)) to yield the desired biomass growth kinetic equation:
(1.27)
Eqn (1.27) has a model structure equivalent to the well-known Monod model. The physical meaning of the kinetic parameters can now be evaluated from a mechanistic perspective:
  • The rate constant k2 (h−1) is analogous to the maximum specific biomass growth rate. We notice that, for cultures where the substrate is present in excess (i.e.Sk2k1), the biomass growth rate equation simply reduces to dXdt=k2X; this of course corresponds to a situation where all of the biomass is saturated with the substrate (X = XS).

  • The term k2k1(gL1) is typically known as the Monod constant, and it is representative of the affinity of the substrate to the cellular metabolism. We observe that if the assimilation of the substrate is much faster than the desorption of growth products (i.e. k1 ≫ k2), we once again obtain the maximum specific growth rate. This constant is also referred to as the half-saturation constant, given that when S=k2k1 the growth rate is half the maximum (XS=12X).

A similar procedure can be carried out to obtain the typical substrate consumption model:
(1.28)
The term Ysx has units of (gs gx−1) and it is commonly known as the substrate yield coefficient. Simply put, this is a mass-stoichiometric ratio that links the rate of substrate utilisation to the rate of biomass growth. In many instances, this model structure is sufficient to describe substrate consumption of growing microbial cultures, although in some cases further consideration must be given to maintenance-related substrate utilisation.

Owing to physiological dissimilarities between different microorganisms, it is ultimately not possible to derive a “master” unstructured model that captures the growth and substrate utilisation dynamics of any microbial population. In some cases, the same microbial culture might even exhibit different dynamic behaviours over a range of conditions due to the influence of metabolic regulations and other intracellular mechanisms.8  Table 1.2 summarises some typical unstructured models in the literature, many of which can also be interpreted as kinetic models for enzymatic reactions.

Table 1.2

Unstructured models for biomass growth and substrate utilisation.

Description Model structure Comments and underlying assumptions
Monod model with cellular decay  dXdt=μmaxSS+KsXμdX   The rate of endogenous biomass decay is assumed to be first order with respect to biomass concentration. 
Multiplicative Monod model  dXdt=μmaxS1S1+K1S2S2+K2X   Both substrates are required and hence the growth is co-limited. 
Additive Monod model  dXdt=(μmax1S1S1+K1+μmax2S2S2+K2)X   Substrates are substitutable and contribute individually to growth. 
Additive Monod model with competitive substrate inhibition  dXdt=μmax1S1S1+K1+ki,1S2X+μmax2S2S2+K2+ki,2S1X   Can be derived by assuming the substrates compete for the same enzyme intermediate.9  Notice the similarity with Langmuirian model structures. 
Contois model  dXdt=μmaxSS+KcXX   Mechanistic derivations considering diffusional barriers and cell flocculation have been proposed.10   
Haldane/Aiba/Andrew model for substrate inhibition  dXdt=μmaxSS+Ks+S2KiX   High substrate concentrations can inhibit growth; for instance, this structure has been used to model photoinhibition in photosynthetic organisms.11   
Tesser model  dXdt=μmaxSnSn+KsX   Analogous to a Monod model with exponential concentration dependences. Can be derived by considering cooperative binding mechanisms.12   
Substrate utilisation with maintenance  dSdt=YsxμmaxSS+KsXmsX   The rate of maintenance-related substrate utilisation is assumed proportional to biomass concentration. 
Description Model structure Comments and underlying assumptions
Monod model with cellular decay  dXdt=μmaxSS+KsXμdX   The rate of endogenous biomass decay is assumed to be first order with respect to biomass concentration. 
Multiplicative Monod model  dXdt=μmaxS1S1+K1S2S2+K2X   Both substrates are required and hence the growth is co-limited. 
Additive Monod model  dXdt=(μmax1S1S1+K1+μmax2S2S2+K2)X   Substrates are substitutable and contribute individually to growth. 
Additive Monod model with competitive substrate inhibition  dXdt=μmax1S1S1+K1+ki,1S2X+μmax2S2S2+K2+ki,2S1X   Can be derived by assuming the substrates compete for the same enzyme intermediate.9  Notice the similarity with Langmuirian model structures. 
Contois model  dXdt=μmaxSS+KcXX   Mechanistic derivations considering diffusional barriers and cell flocculation have been proposed.10   
Haldane/Aiba/Andrew model for substrate inhibition  dXdt=μmaxSS+Ks+S2KiX   High substrate concentrations can inhibit growth; for instance, this structure has been used to model photoinhibition in photosynthetic organisms.11   
Tesser model  dXdt=μmaxSnSn+KsX   Analogous to a Monod model with exponential concentration dependences. Can be derived by considering cooperative binding mechanisms.12   
Substrate utilisation with maintenance  dSdt=YsxμmaxSS+KsXmsX   The rate of maintenance-related substrate utilisation is assumed proportional to biomass concentration. 
With reference to product formation kinetics, the Leudeking–Piret equation13  was first developed to model lactic acid accumulation in bacterial cultivations and since then it has been extensively used to describe product synthesis in various fermentation systems. The premise behind this model is that product accumulation is a result of growth-related and growth-independent metabolite synthesis, such that
(1.29)
where P is the product concentration (g L−1), α is the growth-associated product yield coefficient (gp gx−1), and β is the growth-independent specific product formation rate (gp gx−1 h−1). Several modifications of the classical Leudeking–Piret model have been proposed, for instance, to simulate product synthesis in xanthan gum fermentation,14  carotenoid accumulation in yeast cultivations,15  or waste-product formation in mammalian cell cultures,16  among many others.

Optimisation plays a vital role in both development and application of mathematical process models. This section aims to provide an introduction to the formulation and solution of optimisation problems, two recurring themes in this book due to their role in model parameter estimation, model-based design of experiments and reactor optimisation.

Constrained optimisation problems have three main components: (i) the first is an objective function to be minimised or maximised with respect to its optimisation variables. In the context of reaction engineering, this could for example be the minimisation of by-product yield or maximisation of conversion by finding the optimum operating conditions or feed concentrations in model-based design of experiments. (ii) The second is equality constraints. These are unit operation or process models (such as mass–energy balances and rate equations) that capture the input–output nature of the system being optimised. (iii) The third is inequality constraints, which bound the variables of the objective function to a feasible region. These may, for instance, constrain the optimisation variables to only take physically meaningful values (e.g. bounding temperature or concentrations to non-negative values). By convention, optimisation is formulated as a minimisation problem, as shown in eqn (1.30)(1.32). We note however that any maximisation problem can be equivalently transformed to a minimisation problem by taking the negative of the objective function.
(1.30)
(1.31)
(1.32)
where x ∈ n are the optimisation variables, f : n →  is the objective function to be minimised, h : n → m are the equality constraints, and g : n → r are the inequality constraints. Following the Karush–Kuhn–Tucker (KKT) relaxation approach, the constrained problem can be reformulated into a Lagrangian function which explicitly penalises the violation of the constraints:
(1.33)
where λ ∈ m and μ ∈ r are the Lagrange multipliers for the equality and inequality constraints, respectively. The stationarity KKT condition states that, if a local minimum x* exists for the constrained problem, then there exist multipliers λ* and μ* that minimise the Lagrangian function such that
(1.34)
Additionally, the following conditions which are related to the constrained nature of the problem must also be satisfied at the optimum:
(1.35)
(1.36)
(1.37)
The equations in eqn (1.35) are known as the primal feasibility conditions and simply correspond to the constraints of the original problem. The dual feasibility condition (eqn (1.36)) together with the complementary slackness condition (eqn (1.37)) indicates whether an inequality constraint is active at the optimum; namely, we notice that for an active inequality constraint gi(x*) = 0 and the corresponding dual variable μi*>0.

For most constrained optimisation problems in reaction engineering, the KKT conditions in eqn (1.33)(1.37) cannot be solved analytically. In such cases, optimisation algorithms must be used, whereby the original optimisation problem is iteratively solved through a sequence of relaxed subproblems, in an attempt to find a solution that satisfies the KKT conditions. To demonstrate this, interior-point optimisation algorithms are discussed in the next subsection.

A wide variety of approaches can be adopted to solve nonlinear programming problems. Here we will deal with deterministic (or gradient-based) optimisation methods, as they guarantee the optimality of the solution and converge satisfactorily for optimisation problems with a relatively low number of decision variables.17  This is of course a very general statement and it should be noted that stochastic methods have also found applications in the chemical engineering field, particularly for optimisation problems involving non-differentiable process models and highly stiff objective functions.18 

Within the class of deterministic approaches, interior-point methods have been successfully applied in reaction engineering to solve nonlinear optimisation problems such as bioprocess kinetic parameter estimation15  and model predictive control of jacketed continuous stirred tank reactors.19  To illustrate their application, consider the following constrained optimisation problem:
(1.38)
(1.39)
(1.40)
where x are the optimisation variables, f : n →  is the objective function to be minimised, and h : n → m are the equality constraints (with n ≥ m). Interior point methods first reformulate eqn (1.38)(1.40) as a barrier problem where the inequality constraints have been relaxed:20 
(1.41)
(1.42)
where β is the barrier parameter. The purpose of the barrier term is quite intuitive; as we approach the boundaries of the feasible set (i.e. as xi ≈ 0 in this example), the term approaches infinity, and so the expanded objective function is continuous only within the feasible region enclosed by the inequality constraints. Given that the gradient of the barrier term “pushes” the search direction away from the constraint boundaries, this also implies that descent methods will traverse the interior of the feasible region (hence the name “interior-point”).
The programming problem with relaxed inequality constraints can be solved approximately for a given value of β. To demonstrate this, let us consider the Karush–Kuhn–Tucker conditions for the reformulated barrier problem; we have a stationarity condition for the first derivatives of the Lagrangian, and a primal feasibility condition that enforces the equality constraints:
(1.43)
(1.44)
where L(x,λ) is the Lagrangian of the reformulated optimisation problem, λ ∈ m are the Lagrange multipliers of the equality constraints. We introduce the dual variable z such that zi=βxi to obtain a set of modified (or perturbed) KKT conditions:
(1.45)
(1.46)
(1.47)
where X and Z are the diagonal matrices of the vectors x and z, and e is a vector of ones [1, 1, 1,…] so that the resulting vector dimensions of eqn (1.47) are appropriate. It can be verified that, for β = 0, the perturbed KKT conditions are equivalent to the KKT conditions for the original constrained problem in eqn (1.38)(1.40). Consequently, we can converge to a local solution for the original problem by solving a sequence of barrier problems whereby the barrier parameter β decreases (approaching 0) after each iteration.21  Interior-point algorithms numerically solve the set of perturbed KKT conditions for a fixed value of β, after which the procedure is repeated with an updated barrier parameter and initial guess; for instance, the system of nonlinear equations (eqn (1.45)(1.47)) can be linearly approximated such that a single step of Newton’s method yields
(1.48)
(1.49)
(1.50)
(1.51)
where k denotes the iterate of the Newton sub-algorithm, j denotes the iterate in the overall algorithm, α is a damping factor that scales the computed search directions d k x , d k λ , d k z for each iteration. Naturally, the approximated solution of the jth barrier problem is used to initialise the numerical approximation sub-algorithm for the (j + 1)th problem.

There are a number of considerations that fall beyond the scope of this discussion. Namely, the careful reader will have already noticed that we have not elucidated how the barrier parameter β is updated after each iteration, or what convergence criteria are used to terminate both the numerical approximation sub-algorithm and the overall optimisation algorithm. Moreover, state-of-the-art interior point optimisation methods such as open source IPOPT incorporate further improvements, like line-search filter methods and second-order corrections, to provide more robust convergence. For a detailed breakdown of interior-point methods, we would guide the reader to ref. 22.

To verify that a proposed kinetic model is indeed an accurate representation of the reacting system in question, the model parameters must be fitted against some experimental data. Parameter estimation can be formulated as an optimisation problem whereby the sum of squared residuals is the objective function to be minimised (eqn (1.52)), subject to nonlinear process constraints (eqn (1.53)) and parameter bound constraints (eqn (1.54)):
(1.52)
(1.53)
(1.54)
where N is the number of datapoints, xi,E and xi,M refer to experimental measurements and model estimations of the state variables, respectively, and f(x,k) is the kinetic model describing the rate of change of states. klb and kub are the lower and upper bounds of the optimisation variables k, which can be set on the basis of physical understanding (e.g. rate constants can only be non-negative) or prior literature estimates if available. The weighting matrix Λ is used to scale the residuals, such that the optimal solution is not biased by the relative magnitudes of the states.

More often than not, we do not know a priori the ground truth that comprises the underlying reaction mechanism and its associated rate constants. Moreover, multiple model structures can provide similar fitting accuracy despite incorporating vastly different mechanistic assumptions in their derivations. It is then essential to combine the parameter estimation framework alongside statistical model discrimination methods and model-based design of experiments in order to identify the most plausible model structure. This is further discussed in Chapter 4.

In some instances, modified parameter estimation problems must be formulated through explicit regularisation techniques in order to avoid overfitting. This involves the incorporation of an additional term to the objective function which penalises the violation of a “prior belief” we may have on the system. For instance, consider the case where two kinetic data sets have been gathered at similar (but not equal) temperatures; it is obvious that the fitted kinetic parameters will assume different values across the data sets owing to their temperature dependence, but we can introduce a regularisation term to the objective function to penalise large discrepancies between them (as the temperature range is narrow):
(1.55)
where Nk is the number of kinetic parameters in the model and the subscripts 1 and 2 are used to distinguish between the parameters for the first and second data sets, respectively. The penalty weight ω must be carefully tuned such that the model accuracy is not compromised.

Such modified parameter estimation problems are adopted in Chapter 3 to enforce smoothness of time-varying kinetic parameters by preventing them from changing drastically over short time intervals, or in Chapter 8 where a regularisation term is used to penalise the number of non-zero reaction constants in a complex microkinetic model (thus aiding in model structure identification).

To solve the parameter estimation problem (eqn (1.52)(1.54)) formulated in the previous section, the differential process constraints in eqn (1.53) must be first converted into algebraic equality constraints that correspond to the time-profile of the state variables. Generally, the analytical solution of the differential constraints is not known, and so numerical discretisation methods must be used. The choice of discretisation algorithm is not trivial as generally there is a trade-off between the numerical accuracy of the approximated derivative and the computational cost of the algorithm.23  In this section, we will discuss orthogonal collocation over finite elements; this discretisation method is well-established and allows for accurate approximation of ODEs with relatively few finite elements, therefore being very favourable for optimisation problems with multiple state variables.24 

To illustrate the application of orthogonal collocation on finite elements, consider the initial value problem defined by the differential process constraints and the initial conditions of the measured states:
(1.56)
Here x is the unknown function of the states with respect to time that we want to approximate. The discretisation scheme divides the continuous time domain t into finite intervals [ti,ti+1]; over each of these finite intervals, the time-profile of the states is interpolated across collocation points j = 0, 1,…, K using a Lagrange interpolation polynomial of degree K:25 
(1.57)
where
(1.58)
Here x i K (t) is the Lagrange interpolation polynomial used to approximate the states over the finite element i, x i,j refers to the value of the states at collocation point j of finite element i, and Lj(τ) is the Lagrangian basis function defined in terms of a dimensionless time variable τ.
So far we only have a polynomial description of the time dependence of the state variable x(t). A necessary condition is that, at the collocation points (τi,j, x i,j ), the derivative of this polynomial approximation is equivalent to the true derivative of the state variable. Thus, we have
(1.59)
Substituting the definition of the Lagrange interpolation polynomial into eqn (1.59) and rearranging yields the following expression for the discretised profile of the state variables over a finite element i:
(1.60)
Noting that, for a chosen number of collocation points and finite elements, the term dLjdτ|τj is simply a constant that can be pre-computed for each collocation point. Finally, we must enforce continuity of our discretised profiles across the different finite elements over which we are collocating. That is, the initial value of the states at a finite element must be equal to the value of the states at the last collocation point of the previous finite element:
(1.61)
This continuity constraint together with eqn (1.60) allows for the discretisation of the differential expression in eqn (1.56) to yield algebraic equations for the states at the collocation points xi,j which are only in terms of the initial conditions x0, the optimisation variables k and the finite discretisation step size (ti+1 − ti). Thus, the parameter estimation problem previously formulated in eqn (1.52)(1.54) can be rewritten as
(1.62)
(1.63)
(1.64)
where h(x,k) are the algebraic equality constraints that result from the discretisation of the differential process model, thus transforming the parameter estimation problem into a non-linear programming problem. This, in turn, allows us to apply the interior-point methodology outlined in Section 1.2.2 to estimate the kinetic parameter values k* that minimise the constrained least squares function.

The principles of physical model construction for chemical and biochemical reaction processes have been described in this chapter. First, the development of mechanistic kinetic models for simple chemical reaction networks was elucidated through the application of the mass-action law and the pseudo-steady-state hypothesis. We then discussed how some typical biochemical kinetic models can be similarly derived by approximating the cellular metabolism to enzymatic reaction schemes. Finally, the fundamentals of deterministic optimisation were introduced with particular focus on the formulation and solution of constrained optimisation problems.

1
Rajadhyaksha
 
R. A.
Doraiswamy
 
L. K.
Falsification of kinetic parameters by transport limitations and its role in discerning the controlling regime
Catal. Rev.
1976
, vol. 
13
 
1
(pg. 
209
-
258
)
2
Motagamwala
 
A. H.
Dumesic
 
J. A.
Microkinetic modeling: A tool for rational catalyst design
Chem. Rev.
2021
, vol. 
121
 
1
(pg. 
1049
-
1076
)
3
Sircar
 
S.
Hufton
 
J. R.
Why does the linear driving force model for adsorption kinetics work?
Adsorption
2000
, vol. 
6
 (pg. 
137
-
147
)
4
Montesinos
 
J. L.
Lafuente
 
J.
Gordillo
 
M. A.
Valero
 
F.
Solà
 
C.
Charbonnier
 
S.
Cheruy
 
A.
Structured modeling and state estimation in a fermentation process: Lipase production by Candida rugosa
Biotechnol. Bioeng.
1995
, vol. 
48
 
12
(pg. 
573
-
584
)
5
Ramkrishna
 
D.
Song
 
H.-S.
Dynamic models of metabolism: Review of the cybernetic approach
AIChE J.
2012
, vol. 
58
 
4
(pg. 
986
-
997
)
6
Birol
 
G.
Kirdar
 
B.
Onsan
 
Z. I.
A simple structured model for biomass and extracellular enzyme production with recombinant saccharomyces cerevisiae ypb-g
J. Ind. Microbiol. Biotechnol.
2002
, vol. 
29
 
9
(pg. 
111
-
116
)
7
Liu
 
Y.
Overview of some theoretical approaches for derivation of the Monod equation
Appl. Microbiol. Biotechnol.
2007
, vol. 
73
 
1
(pg. 
1241
-
1250
)
8
Adnan Jouned
 
M.
Kager
 
J.
Herwig
 
C.
Barz
 
T.
Event driven modeling for the accurate identification of metabolic switches in fed-batch culture of s. cerevisiae
Biochem. Eng. J.
2022
, vol. 
180
 
3
pg. 
108345
 
9
Yoon
 
H.
Klinzing
 
G.
Blanch
 
H. W.
Competition for mixed substrates by microbial populations
Biotechnol. Bioeng.
1977
, vol. 
19
 
8
(pg. 
1193
-
1210
)
10
Wang
 
Z.-W.
Li
 
Y.
A theoretical derivation of the contois equation for kinetic modeling of the microbial degradation of insoluble substrates
Biochem. Eng. J.
2014
, vol. 
82
 
1
(pg. 
134
-
138
)
11
Zhang
 
D.
Dechatiwongse
 
P.
del Rio-Chanona
 
E. A.
Maitland
 
G. C.
Hellgardt
 
K.
Vassiliadis
 
V. S.
Modelling of light and temperature influences on cyanobacterial growth and biohydrogen production
Algal Res.
2015
, vol. 
9
 
5
(pg. 
263
-
274
)
12
Stefan
 
M. I.
Le Novère
 
N.
Cooperative Binding
PLoS Comput. Biol.
2013
, vol. 
9
 
6
pg. 
e1003106
 
13
Luedeking
 
R.
Piret
 
E. L.
A kinetic study of the lactic acid fermentation. batch process at controlled ph
J. Biochem. Microbiol. Technol. Eng.
1959
, vol. 
1
 
12
(pg. 
393
-
412
)
14
Weiss
 
R. M.
Ollis
 
D. F.
Extracellular microbial polysaccharides. I. Substrate, biomass, and product kinetic equations for batch xanthan gum fermentation
Biotechnol. Bioeng.
1980
, vol. 
22
 
4
(pg. 
859
-
873
)
15
Ramon
 
F. V.
Zhu
 
X.
Savage
 
T. R.
Petsagkourakis
 
P.
Jing
 
K.
Zhang
 
D.
Kinetic and hybrid modeling for yeast astaxanthin production under uncertainty
Biotechnol. Bioeng.
2021
, vol. 
118
 
12
(pg. 
4854
-
4866
)
16
Zeng
 
A.-P.
A kinetic model for product formation of microbial and mammalian cells
Biotechnol. Bioeng.
1995
, vol. 
46
 
5
(pg. 
314
-
324
)
17
Liberti
 
L.
Kucherenko
 
S.
Comparison of deterministic and stochastic approaches to global optimization
Int. Trans. Oper. Res.
2005
, vol. 
12
 
5
(pg. 
263
-
285
)
18
J. R.
Banga
and
W. D.
Seider
, Global optimization of chemical processes using stochastic algorithms, in
State of the Art in Global Optimization
,
Springer
,
1996
, pp.
563
583
.
19
Albuquerque
 
J.
Gopal
 
V.
Staus
 
G.
Biegler
 
L. T.
Ydstie
 
B. E.
Interior point sqp strategies for large-scale, structured process optimization problems
Comput. Chem. Eng.
1999
, vol. 
23
 
5
(pg. 
543
-
554
)
20
Byrd
 
R. H.
Gilbert
 
J. C.
Nocedal
 
J.
A trust region method based on interior point techniques for nonlinear programming
Math. Program.
2000
, vol. 
89
 
11
(pg. 
149
-
185
)
21
Biegler
 
L. T.
An overview of simultaneous strategies for dynamic optimization
Chem. Eng. Process.
2007
, vol. 
46
 
11
(pg. 
1043
-
1053
)
22
Wächter
 
A.
Biegler
 
L. T.
On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming
Math. Program.
2006
, vol. 
106
 
3
(pg. 
25
-
57
)
23
Wu
 
H.
Xue
 
H.
Kumar
 
A.
Numerical discretization-based estimation methods for ordinary differential equation models via penalized spline smoothing with applications in biomedical research
Biometrics
2012
, vol. 
68
 
6
(pg. 
344
-
352
)
24
Carey
 
G. F.
Finlayson
 
B. A.
Orthogonal collocation on finite elements
Chem. Eng. Sci.
1975
, vol. 
30
 
5
(pg. 
587
-
596
)
25
Nicholson
 
B.
Siirola
 
J. D.
Watson
 
J.-P.
Zavala
 
V. M.
Biegler
 
L. T.
Pyomo.dae: a modeling and automatic discretization framework for optimization with differential and algebraic equations
Math. Program. Comput.
2018
, vol. 
10
 
6
(pg. 
187
-
223
)
Close Modal

or Create an Account

Close Modal
Close Modal