Anton A. Kiss, Edwin Zondervan, Richard Lakerveld, Leyla Özkan (Eds.

)
Proceedings of the 29th European Symposium on Computer Aided Process Engineering
June 16th to 19th , 2019, Eindhoven, The Netherlands. 
c 2019 Elsevier B.V. All rights reserved.
http://dx.doi.org/10.1016/B978-0-12-818634-3.50129-6

Nonlinear dynamic analysis of chemical engineering

processes described by differential-algebraic equations
systems
Ataı́de S. Andrade Netoa,* , Argimiro R. Secchia and Prı́amo A. Meloa
a Chemical
Engineering Program, COPPE, Universidade Federal do Rio de Janeiro, Brazil
ataide@peq.coppe.ufrj.br

Abstract
Chemical processes are subject to parametric variation due to, among others things, desired changes
in process conditions, instrumentation malfunction or unpredicted disturbances. The outcome of
these adversities is plural, the entire system may collapse or the output variables might not be
affected at all. Thus, the theoretical analysis of processes prone to nonlinear responses plays a key
role in design and operation of chemical plants either for safety or economic purposes. Because
of the widespread use in the last couple decades of the differential-algebraic equations (DAE) ap-
proach in chemical processes modeling, the development of computational algorithms and tools
specifically designed for the analysis of DAE systems becomes necessary, specially for high-index
cases. In this work, we present a novel and index-free approach for the nonlinear dynamic analysis
of mathematical models based on DAEs. A methodology for the direct detection of Hopf bifurca-
tion points based on an optimization procedure, in which we rewrite the necessary condition for
the Hopf bifurcation point occurrence as a constrained optimization problem was developed. In
this approach, we do not compute the generalized eigenvalues of the system, but their Hurtwitz
determinants using a new algorithm for the computation of characteristic polynomials of matrix
pencils based on the La Budde’s method. The resulting optimization problem was solved with
a hybrid global optimization technique, coupling particle swarm (PSA) and active-set (ASA) al-
gorithms. The parametric continuation of steady-state branches was conducted with the pseudo
arc-length method and the stability analysis based on the spectral theory for linear matrix pencils.
As an illustrative example we evaluate the dynamic behavior of a reactive flash drum formulated
as an index-2 DAE system. The results show that the proposed methodology can handle DAE
systems of high index in a fast and accurate way.

Keywords: Hopf bifurcation, differential-algebraic equations, stability analysis.

1. Introduction
Since the advent of calculus, the mathematical modeling of dynamic physicochemical processes
has been based on ordinary and partial differential equations, with a very few exceptions (no-
tably, the works of Lagrange on constrained mechanics dated from 1788 and Kirchhoff on electric
circuits from 1847). However, with the growing studies on the differential-algebraic equations
systems, started in the early 1970s, and also because their vast applicability this paradigm has
been changing. Several authors have already pointed out numerous advantages of working di-
rectly with DAE; but we reinforce that formulating differential-algebraic models in many practical
applications occurs naturally, specially in chemical engineering, as thermodynamic and equilib-
rium constraints, which are all imposed by algebraic equations, may be directly incorporated in
770 A. S. Andrade Neto et al.

the process modeling. Because of this feature, the DAE approach is largely exploited by the ma-
jority of dynamic simulation software used in the context of Computer Aided Process Engineering
(CAPE).
Regarding the dynamic analysis of chemical processes, we can safely say that it is a well-developed
field for systems governed by ordinary differential equations (ODE) or discrete maps, specially
the bifurcation theory which contemplates the study of dynamic systems under variation in their
evolution laws. From a historical point of view, it is possible to recognize the importance of the
bifurcation theory developed for ODE as a fundamental tool for nonlinear analysis. Unfortunately,
for systems described by DAE, the same is not true. Because it is a rather new topic compared to
ODE systems, computational algorithms and tools specifically designed for DAE have not evolved
so steadily yet. It is noticed, nevertheless, that some works in the literature have already addressed
the bifurcation analysis for DAE systems using either classical algorithms for differential equa-
tions or extensions of them for index-1 systems (Clausbruch et al., 2006), which are far much
simpler than high-index DAE.
Although it is possible to reduce some high-index DAE to index-1 systems, this can be a tedious
error-prone procedure and, in many cases, the solution of the reduced system may significantly
differ from that of the original problem due to the well-known drift-off effect. Additionally, Har-
ney et al. (2013) showed that the stability analysis of the reduced system may not represent the
stability of the higher-index system from which it is derived, often being necessary to apply some
stabilization methodologies, which has been developed for DAEs in Hessenberg form of size 2.
Knowing beforehand the precise location where nonlinear phenomena, such as bifurcation points,
occur is of great importance in real-life applications, as not only safety but also economic issues
may be avoided. Also, because the DAE approach is the current paradigm in the theoretical anal-
ysis of chemical processes, computational methods that efficiently handle DAE systems becomes
necessary. Driven by this, we present, in this work, a novel and index-free approach for the direct
computation of Hopf bifurcation points (HBPs), without the necessity of applying index reduction
algorithms or stabilization procedures.

2. Theoretical framework
2.1. Linearization of differential-algebraic equations systems
The core of bifurcation theory of nonlinear system inevitably falls back to the dynamic analysis
of linear ones. Because of that, the fundamental question one may ask is if there exist a linearized
DAE system with the same qualitative behavior around fixed points of its nonlinear counterpart.
Reich (1995) proved that the existence of such systems depends on some conditions which will be
presented in the sequence.
Consider the nonlinear index-ν DAE system in the fully implicit form
f (x(t), ẋ(t); p) = 0 (1)
with f : Rn × Rn × Rm → Rn and p ∈ Rm a vector of time-invariant parameters. A linearization of
f around a fixed point xe in the form
B(p)ẏ(t) = A(p)y(t) (2)
where A = ∇x f T (xe , 0; p), B = ∇ẋ f T (xe , 0; p) and y = x − xe , will exist only if the original DAE
in Eq. (1) is regular (or solvable) and fulfill the following conditions:
rank [∇ẋ f T (x, ẋ; p)] = constant for all t (3a)
rank [∇x f (xe , 0; p)] = n
T
(3b)
rank [∇x gTν (xe , 0; p)] = n (3c)
Nonlinear dynamic analysis of chemical engineering processes described by 771
differential-algebraic equations systems

in which gν (x, ẋ; p) is a derivative-array obtained by stacking the vector f with all the equations
that arise in an index reduction procedure. According to Campbell (1995), those conditions are
simple to evaluate and often observed in physical systems.

2.2. Stability in DAE systems

Several authors have shown that the stability of a linear DAE is closely related to the spectrum,
σ (A, B), of the matrix pair (A, B) as defined in Eq. (2) (see Reich (1995) for a detailed overview on
this subject). The spectrum of a DAE have exactly nd finite values; the other n − nd values is said
to be at the infinity. Here, nd represent the number of Dynamics Degree of Freedom (DDoF) and
n the dimension of the system. If all the finite values in σ (A, B) lie in the left complex half-plane,
then, the fixed point is stable.

2.3. Detection of Hopf bifurcation points

Hopf bifurcation points can be detected by an indirect method by tracing a branch of steady state
solutions of an ODE system, usually done by homotopy continuation methods, while monitoring
the eigenvalues of Jacobian matrix at the fixed points. When a single pair of eigenvalues crosses
the imaginary axis a trapping method, such as the bisection or the secant method, can be used
to refine the HBP. This procedure, despite being effective and easily extended to handle DAE
systems, has a high computational cost. Alternatively, HBPs can be calculated directly. The
first successful algorithm for the direct computation of Hopf bifurcation points in ODE systems
was presented by Griewank and Reddien (1983) by formulating a set of 3n + 2 equations that
algebraically determines these points. Later, Reich (1995) extended this algorithm to the DAE
systems.
We present a different approach in which instead of an augmented algebraic system, we solve a
constrained nonlinear optimization problem. The main advantage of this method is that we do not
restrict the Hopf bifurcation points search to a single parameter of the model. In fact, our algorithm
works for any number of parameters. Also, we do not compute the spectrum of the DAE system,
but its nd − 1 Hurwitz determinants, Dk .
The Hurwitz determinants can be calculated from the normalized characteristic polynomial:
P(λ ; p) = det (λ B − A)/ānd (p) = a0 (p) + a1 (p)λ + a2 (p)λ 2 + · · · + λ nd (4)
by building the Hurwitz matrix:
⎡ ⎤
a1 (p) a0 (p) 0 ··· 0
⎢a3 (p) a2 (p) a1 (p) · · · 0⎥
⎢ ⎥
⎢ 0⎥
L(p) = ⎢a5 (p) a4 (p) a3 (p) · · · ⎥ (5)
⎢ .. .. .. .. .. ⎥
⎣ . . . . .⎦
0 0 0 ··· 1
then, evaluating its first nd − 1 principal subdeterminants:
⎡ ⎤
l1,1 · · · l1,k
⎢ .. ⎥
Dk (p) = det(Lk ) = det ⎣ ... ..
. . ⎦ (6)
lk,1 ··· lk,k

where ai = āi /ānd and āi (p) are the coefficients of P̄(λ ; p) = det(λ B − A); L ∈ Rnd ×nd , Lk ∈ Rk×k ,
Dk ∈ R and k = 1, 2, 3, . . . , nd − 1.
Extending the theorem of Liu (1994) for DAE systems, a Hopf bifurcation point arises when:
a0 > 0, D1 > 0, D2 > 0, . . . , Dnd −2 > 0, Dnd −1 = 0 (7)
772 A. S. Andrade Neto et al.

This allows us to formulate the following constrained optimization problem:

min S(p) = D2nd −1 subject to:
p∈Pm , x∈Rn

a0 > 0, D1 > 0, . . . , Dnd −2 > 0 (8)

f (x, ẋ; p) = 0
where Pm ⊆ Rm is a selected parametric space and S : Pm → R. The global minimum of this
problem (S(p) = 0) is always a Hopf bifurcation point.

2.4. Characteristic polynomials of matrix pencils

A general matrix pair (A, B) can be decomposed by Householder reflections into the form A =
QHZ and B = QT Z where Q and Z are orthogonal, H upper Hessenberg and T upper triangular
n × n matrices. Then,
λ B − A = λ QT Z − QHZ ⇒ λ QT QT ZZ T − QT QHZZ T = λ T − H (9)
Solving det(λ T − H) is numerically more efficient than solving det(λ B − A), because it can be
rapidly evaluated by the following recursion formula derived from the La Budde method (Rehman
and Ipsen, 2011):
i−1  i−1 
P̄i (λ ; p) = (γi λ − αi )P̄i−1 (λ ; p) + ∑ (ti−k,i λ − hi−k,i )P̄i−k−1 (λ ; p) ∏ βm (10)
k=1 m=i−k

with P̄0 (λ ; p) = 1 and P̄1 (λ ; p) = γ1 λ − α1 ; i = 2, 3, . . . , n. The coefficients αi and βi are,

respectively, the diagonal and the subdiagonal of H; γi the diagonal of T ; ti, j and hi, j are the
(i, j)-elements of T and H; P̄n (λ ; p) = P̄(λ ; p) = det(λ T − H) = det(λ B − A).

3. Results
In order to evaluate the proposed methodology, it was applied to a benchmark model in chemical
engineering. Consider the index-2 system describing the chemical reaction A → B → C conducted
in reactive flash drum:
⎡    ⎤
−g1
⎢ x˙A − x A f + x A  + KA φ + D a1 exp ⎥
⎢ u+1
⎢      ⎥ ⎥
⎢ −g2 −g1 ⎥
⎢ x˙B − xB f + xB  + KB φ + Da2 exp − xA Da1 exp ⎥
⎢ u + 1 u + 1 ⎥
⎢     ⎥
⎢ −g −g ⎥
f (x, ẋ; p) = ⎢ u̇ − 1 + u + φ δ − Q − H D x exp 1
− H2 Da2 xB exp
2 ⎥ = 0 (11)
⎢ 1 a1 A ⎥
⎢ u+1 u+1 ⎥
⎢ ⎥
⎢ KA xA + KB xB + KC xC − 1 ⎥
⎢ ⎥
⎢ ⎥
⎣ xA + xB + xC − 1 ⎦
+φ −1
where xA , xB and xC are molar fractions, u the dimensionless temperature and φ and  the vapor
and liquid fractions; the vector of search parameters is pT = [Da1 , Da2 , g1 , g2 , δ , Q, H1 , H2 ].
The feeding compositions xA f and xB f were set to 0.8 and 0.2, respectively, and the equilibrium
constants Ki = Pisat /Pref . The saturation pressure was modeled by the Antoine equation
 
βi
Pisat = exp αi + (12)
(u + 1)Tref − γi
for some reference chemical species with the following parameters: α T : [21.3, 23.2, 25.1]; β T :
[−2428.2, −3835.2, −6022.2]; γ T = [35.4, 45.3, 25.3]. The reference temperature and pressure
were set to 298 K and 101, 325 Pa.
Nonlinear dynamic analysis of chemical engineering processes described by 773
differential-algebraic equations systems

The optimization problem in Eq. (8) for the Table 1: Average performance of optimization
model described by Eq. (11) was solved us- algorithms† .
ing a hybrid approach by coupling the parti-
cle swarm (PSA) and active-set (ASA) algo- Parameters PSA†† ASA†† Hybrid††
rithms. The PSA was applied at the start of the
optimization procedure until a predetermined all 129.0/4 −/10 12.1/0
value of the objective function was reached, Da1 , Q 53.1/0 1.9/0 13.6/0
S(pm ) = εm . Then, the ASA was called to D a1 , δ 56.6/0 1.9/0 12.9/0
refine the solution to S(p f ) = ε f , where ε f Q 51.1/0 1.8/0 11.5/0
is a predefined tolerance. This technique can † ε = 0.099; ε = 10−8 ; random initial guess;
m f
†† Time (s) / Number of failures in 10 attempts.
be more efficient and effective for non con-
vex problems with high number of optimiza-
tion variables. In Table 1, we present a comparison between the hybrid approach and the stan-
dalone algorithms in terms of number of successful attempts and spent computational time that
support this assertion. All simulations were done with Matlab on a 8GB RAM, Intel CORE-i7
laptop running Linux. With the exception of the standalone ASA, the PSA and the hybrid ap-
proaches were able to attain the global minimum in every scenario.
An appropriated selection of Pm is necessary for the success of the optimization. The vectors
LT = [0, 0, 0, 0, 0, 0, 0, 0] and U T = [10, 10, 10, 10, 100, 10, 100, 100] represent feasible values for
all the parameters, as they are dimensionless groupings resulting from the combination of physical
properties and operating conditions, therefore, they were chosen to bound the set Pm . As a result
of the optimization problem, we found that a Hopf bifurcation point occurs at:

xT = [0.2920, 0.3469, 0.0681, 0.7466, 0.3611, 0.2534]

pT = [8.6807, 9.9707, 6.7821, 3.8784, 15.5261, 1.3430, 35.9778, 99.8517]

The periodic solutions arising at this HBP are stable, as can be seen in the phase diagram and
the time series in Fig. 1. Another HBP was found at Q = 2.1741, using Q as the single search
parameter while fixing all other parameters. In this case, the arising periodic solutions were found
to be unstable.
0.348 0.37

0.3478 0.36

0.3476 0.35

0.3474 0.34
Molar fraction

0.3472 0.33 xA
xB

xB
0.347 0.32
xC
0.3468 0.31
0.3466 0.3
0.3464 Trajectory
Starting point
0.29
0.3462
0.286 0.288 0.29 0.292 0.294 0.296 0.28
0 200 400 600 800 1000
xA
Time
(a) Phase portrait xB vs. xA . (b) Composition time profile.

Figure 1: Simulation in the neighborhood of Hopf bifurcation point at Q = 1.3430, the stable case.
Total simulation time t = 1000.

The pseudo arc-length method was applied in order to make the one-parameter steady-state con-
tinuation of the model in Eq. (11). The resulting bifurcation diagram of this procedure is presented
in Fig. 2.
774 A. S. Andrade Neto et al.

0.7

0.6

0.5

xA
0.4

0.3

0.2

0.1
0 0.5 1 1.5 2 2.5 3
Q

Figure 2: Bifurcation digram for the model in Eq. (11) using Q as the continuation parameter.
Solid line: stable branch; dotted line: unstable branch; : HBP around stable periodic solutions;
: HBP around unstable periodic solutions; •: bifurcation point.

4. Conclusion
A novel methodology for the direct computation of Hopf bifurcation points in differential-algebraic
equations systems was presented. This strategy is based on the formulation of the constrained opti-
mization problem in Eq. (8) with multiple search parameters. With the proposed methodology, we
were able to find the occurrence of HBPs in a benchmark model with up to eight search parameters
in a fast and accurate manner.
Acknowledgments
This work was supported by the National Council for Scientific and Technological Development
(CNPq) grant numbers 302893/2013-0 and 152572/2016-3. This study was financed in part by the
Coordenação de Aperfeiçoamento de Pessoal de Nı́vel Superior - Brasil (CAPES) - Finance Code
001.
References
K. E. Brenan, S. L. Campbell, L. R. Petzold, 1996. Numerical Solution of Initial-Value Problems in Differential-
Algebraic Equations. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA.
S. L. Campbell, 1995. Linearization of daes along trajectories. Zeitschrift für angewandte Mathematik und Physik ZAMP
46 (1), 70–84.
B. C. Clausbruch, E. C. Biscaia, P. A. Melo, 2006. Stability analysis of differential-algebraic equations in AUTO DAE.
Vol. 21 of Computer Aided Chemical Engineering. Elsevier, pp. 297 – 302.
A. Griewank, G. Reddien, 1983. The calculation of hopf points by a direct method. IMA Journal of Numerical Analysis
3 (3), 295–303.
D. Harney, T. Mills, N. Book, 2013. Numerical evaluation of the stability of stationary points of index-2 differential-
algebraic equations: Applications to reactive flash and reactive distillation systems. Computers & Chemical Engi-
neering 49, 61 – 69.
I. Hyanek, J. Zacca, F. Teymour, W. H. Ray, 1995. Dynamics and stability of polymerization process flow sheets.
Industrial & Engineering Chemistry Research 34 (11), 3872–3877.
W. Liu, 1994. Criterion of hopf bifurcations without using eigenvalues. Journal of Mathematical Analysis and Applica-
tions 182 (1), 250 – 256.
M. Mangold, A. Kienle, E. Gilles, K. Mohl, 2000. Nonlinear computation in diva — methods and applications. Chemical
Engineering Science 55 (2), 441 – 454.
R. Rehman, I. C. F. Ipsen, 2011. La Budde’s Method for Computing Characteristic Polynomials. ArXiv e-prints.
S. Reich, 1995. On the local qualitative behavior of differential-algebraic equations. Circuits, Systems and Signal
Processing 14 (4), 427–443.
B. Simeon, 2017. On the History of Differential-Algebraic Equations. Springer International Publishing, Cham, pp.
1–39.
R. P. Soares, A. R. Secchi, 2005. Direct initialisation and solution of high-index DAE systems. Vol. 20 of Computer
Aided Chemical Engineering. Elsevier, pp. 157 – 162.