• KEYWORDS: Physical Systems

Systems of first-order differential equations

Nonlinear and linear models
• We consider the most general case of dynamical system that arise in the modelling of
physical systems. The majority of the models takes some differentail equation form since
it is build on the notion of states time evolution (e.g. the tank level is a state that changes
in the real world industrial plant). In particular, the setting takes the following form:
dx
=F ( x ( t ) , u ( t ) ) , x ( 0 )=x 0
dt
where x ( 0 )=x 0 is initial condition and u ( t ) is some input to the system or forcing function.
Hence, the x (t ) -states start from the x 0 state at time t=0 and evolve further in time. The above
equation is also equiped with the some expression that usually represents output- y ( t ) of the
system so that one considers the following:

y ( t ) =G ( x ( t ) , u ( t ) )

• Finally, the large number of the system takes the form as it is given by above
equations.

• Therefore, the systems description of underlying dynamics considers states x (t ) ,

inputs u ( t ), outputs y ( t ) and initial conditions in order to completly describe
systems dynamics.

• In general, the consideration of the forward time evolution of the nonlinear system
is still open research topic.

Linear models
• Linear models in general arise by the process of linearization of nonlinear models
around operating point of interest. Hence, the linear models is obtained by
considering the following proecedure:

• Obtain the steady state of interest, this is obtained by considering the solution of the
nonlinear algebraic equation
dx
=0=F ( x s s , us s ) → ( x s s , us s )
dt
• Apply the lienarization at the steady state of interest x s s , us s and express the state as
value of the steady state augmented with the perturbation of the state x p (e.g.
x s s + x p ( t ) =x ( t ), u s s +u p ( t ) =u (t ), y ( t ) = y s s + y p ( t ) ), which yields:

d xp ∂ F
=
dt ∂x x
❑ | ss ,us s x p+
∂F
∂u ) x s s , us s
up
 y p (t)=
∂G
❑
∂x x | ss ,us s x p+
∂G
∂u ) xs s ,us s
up

• Finally, the above expression provides well known linear model representation in
terms of matrics (A,B,C,D) or in the equation form
ẋ p ( t ) ¿ A x p ( t ) + B u p ( t )
y p ( t ) ¿ C x p ( t ) + Du p ( t )

with initial condition x p ( 0 )=x 0.

• Taylor expansion of F ( x ,u ) around the steady state:

[ )
f 1 ( x1 , ⋯ , x n ,u 1 , ⋯ , un )
⋯
F ( x ,u )=
⋯
f n ( x 1 , ⋯ , x n , u1 , ⋯ , un )

[ )[ [ )[
∂f1 ∂f1 ∂f1 ∂f1

[
⋯ ⋯

) ) )
f 1 ( x 1 s s , ⋯ , x n s s , u1 s , ⋯ ,u n s s ) ∂ x1 ∂ x n x 1 ( t ) − x1 s s ∂u1 ∂u n u1 ( t ) −u1 s s
F ( x ,u )= ⋯ + ⋯ ⋯ ⋯ ⋯ + ⋯ ⋯ ⋯ ⋯ + H . O. T .
f n ( x 1 s s ,⋯ , x n s s ,u 1 s s ,⋯ , un s s ) ∂ f n ∂ f n xn ( t ) − xn s s ∂fn ∂ f n un ( t ) −un s s
⋯ ⋯
∂ x1 ∂ xn ∂u1 ∂u n

and the same applys to G ( x ( t ) ,u ( t ) ), so that by using Taylor expansion one can transform
nonlinear to linear models.

Linear models characteristics

• Linear models have important characteristics which we will examine and
charaterize. For example:

• Linear models without the input force function (e.g. u p =0) have unique
equilibrium (or steady state point at zero) providng that matrix A is full rank

• Linear models reveal the fundamental characteristics of the systems such as

stability, controllability and observability (we will discuss in addition controllability
and observability)

• Linear dynamics and state evolution are completely characterized by the

eigenvalues of the matrix A
– Dynamics is charactererized by the distribution of the eigenvalues
• eigenvalues of A are distinct and negative λ i< 0, stable system
• eigenvalues of A are purely complex number (marginaly stable,
λ i=± j ω) and other eigenvalues of A negative
 • eigenvalues of A are at least one positive λ i> 0, unstable
• The linear system dynamics can be described in the following simple computations.

import numpy as np
import math
from numpy import linalg as LA

A = np.array([[0, -1, 0], [1, 0, 0],[0,0,-2]])

print(A)

w, v = LA.eig(A)

print(np.round(w))
print(v)

t0=0
x0=np.array([[1],[1],[1]])

t, h = np.linspace(t0, 15, 1000, retstep=True) # Note the optional

argument to get the stepsize.
x = x0
t[0] = t0
x[[0,1,2],0]=1
print(x )

# Implementation of Euler's method

for n in range(0, len(t) - 1):
x0=x0 + h*np.matmul(A,x0)
x=np.concatenate((x,x0),axis=1)

print(t)
print(x.transpose())

%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(t, x.transpose())
plt.xlabel('t')
plt.ylabel('x')
#plt.xlim([x.min(), x.max()])

[[ 0 -1 0]
[ 1 0 0]
[ 0 0 -2]]
[ 0.+1.j 0.-1.j -2.+0.j]
[[0.70710678+0.j 0.70710678-0.j 0.
+0.j ]
[0. -0.70710678j 0. +0.70710678j 0.
+0.j ]
[0. +0.j 0. -0.j 1.
+0.j ]]
[[1]
[1]
[1]]
[ 0. 0.01501502 0.03003003 0.04504505 0.06006006
0.07507508
0.09009009 0.10510511 0.12012012 0.13513514 0.15015015
0.16516517
0.18018018 0.1951952 0.21021021 0.22522523 0.24024024
0.25525526
0.27027027 0.28528529 0.3003003 0.31531532 0.33033033
0.34534535
0.36036036 0.37537538 0.39039039 0.40540541 0.42042042
0.43543544
0.45045045 0.46546547 0.48048048 0.4954955 0.51051051
0.52552553
0.54054054 0.55555556 0.57057057 0.58558559 0.6006006
0.61561562
0.63063063 0.64564565 0.66066066 0.67567568 0.69069069
0.70570571
0.72072072 0.73573574 0.75075075 0.76576577 0.78078078
0.7957958
0.81081081 0.82582583 0.84084084 0.85585586 0.87087087
0.88588589
0.9009009 0.91591592 0.93093093 0.94594595 0.96096096
0.97597598
0.99099099 1.00600601 1.02102102 1.03603604 1.05105105
1.06606607
1.08108108 1.0960961 1.11111111 1.12612613 1.14114114
1.15615616
1.17117117 1.18618619 1.2012012 1.21621622 1.23123123
1.24624625
1.26126126 1.27627628 1.29129129 1.30630631 1.32132132
1.33633634
1.35135135 1.36636637 1.38138138 1.3963964 1.41141141
1.42642643
1.44144144 1.45645646 1.47147147 1.48648649 1.5015015
1.51651652
1.53153153 1.54654655 1.56156156 1.57657658 1.59159159
1.60660661
1.62162162 1.63663664 1.65165165 1.66666667 1.68168168
1.6966967
1.71171171 1.72672673 1.74174174 1.75675676 1.77177177
1.78678679
1.8018018 1.81681682 1.83183183 1.84684685 1.86186186
1.87687688
1.89189189 1.90690691 1.92192192 1.93693694 1.95195195
1.96696697
1.98198198 1.996997 2.01201201 2.02702703 2.04204204
2.05705706
 2.07207207 2.08708709 2.1021021 2.11711712 2.13213213
2.14714715
2.16216216 2.17717718 2.19219219 2.20720721 2.22222222
2.23723724
2.25225225 2.26726727 2.28228228 2.2972973 2.31231231
2.32732733
2.34234234 2.35735736 2.37237237 2.38738739 2.4024024
2.41741742
2.43243243 2.44744745 2.46246246 2.47747748 2.49249249
2.50750751
2.52252252 2.53753754 2.55255255 2.56756757 2.58258258
2.5975976
2.61261261 2.62762763 2.64264264 2.65765766 2.67267267
2.68768769
2.7027027 2.71771772 2.73273273 2.74774775 2.76276276
2.77777778
2.79279279 2.80780781 2.82282282 2.83783784 2.85285285
2.86786787
2.88288288 2.8978979 2.91291291 2.92792793 2.94294294
2.95795796
2.97297297 2.98798799 3.003003 3.01801802 3.03303303
3.04804805
3.06306306 3.07807808 3.09309309 3.10810811 3.12312312
3.13813814
3.15315315 3.16816817 3.18318318 3.1981982 3.21321321
3.22822823
3.24324324 3.25825826 3.27327327 3.28828829 3.3033033
3.31831832
3.33333333 3.34834835 3.36336336 3.37837838 3.39339339
3.40840841
3.42342342 3.43843844 3.45345345 3.46846847 3.48348348
3.4984985
3.51351351 3.52852853 3.54354354 3.55855856 3.57357357
3.58858859
3.6036036 3.61861862 3.63363363 3.64864865 3.66366366
3.67867868
3.69369369 3.70870871 3.72372372 3.73873874 3.75375375
3.76876877
3.78378378 3.7987988 3.81381381 3.82882883 3.84384384
3.85885886
3.87387387 3.88888889 3.9039039 3.91891892 3.93393393
3.94894895
3.96396396 3.97897898 3.99399399 4.00900901 4.02402402
4.03903904
4.05405405 4.06906907 4.08408408 4.0990991 4.11411411
4.12912913
4.14414414 4.15915916 4.17417417 4.18918919 4.2042042
4.21921922
4.23423423 4.24924925 4.26426426 4.27927928 4.29429429
4.30930931
4.32432432 4.33933934 4.35435435 4.36936937 4.38438438
4.3993994
4.41441441 4.42942943 4.44444444 4.45945946 4.47447447
4.48948949
4.5045045 4.51951952 4.53453453 4.54954955 4.56456456
4.57957958
4.59459459 4.60960961 4.62462462 4.63963964 4.65465465
4.66966967
4.68468468 4.6996997 4.71471471 4.72972973 4.74474474
4.75975976
4.77477477 4.78978979 4.8048048 4.81981982 4.83483483
4.84984985
4.86486486 4.87987988 4.89489489 4.90990991 4.92492492
4.93993994
4.95495495 4.96996997 4.98498498 5. 5.01501502
5.03003003
5.04504505 5.06006006 5.07507508 5.09009009 5.10510511
5.12012012
5.13513514 5.15015015 5.16516517 5.18018018 5.1951952
5.21021021
5.22522523 5.24024024 5.25525526 5.27027027 5.28528529
5.3003003
5.31531532 5.33033033 5.34534535 5.36036036 5.37537538
5.39039039
5.40540541 5.42042042 5.43543544 5.45045045 5.46546547
5.48048048
5.4954955 5.51051051 5.52552553 5.54054054 5.55555556
5.57057057
5.58558559 5.6006006 5.61561562 5.63063063 5.64564565
5.66066066
5.67567568 5.69069069 5.70570571 5.72072072 5.73573574
5.75075075
5.76576577 5.78078078 5.7957958 5.81081081 5.82582583
5.84084084
5.85585586 5.87087087 5.88588589 5.9009009 5.91591592
5.93093093
5.94594595 5.96096096 5.97597598 5.99099099 6.00600601
6.02102102
6.03603604 6.05105105 6.06606607 6.08108108 6.0960961
6.11111111
6.12612613 6.14114114 6.15615616 6.17117117 6.18618619
6.2012012
6.21621622 6.23123123 6.24624625 6.26126126 6.27627628
6.29129129
6.30630631 6.32132132 6.33633634 6.35135135 6.36636637
6.38138138
6.3963964 6.41141141 6.42642643 6.44144144 6.45645646
6.47147147
 6.48648649 6.5015015 6.51651652 6.53153153 6.54654655
6.56156156
6.57657658 6.59159159 6.60660661 6.62162162 6.63663664
6.65165165
6.66666667 6.68168168 6.6966967 6.71171171 6.72672673
6.74174174
6.75675676 6.77177177 6.78678679 6.8018018 6.81681682
6.83183183
6.84684685 6.86186186 6.87687688 6.89189189 6.90690691
6.92192192
6.93693694 6.95195195 6.96696697 6.98198198 6.996997
7.01201201
7.02702703 7.04204204 7.05705706 7.07207207 7.08708709
7.1021021
7.11711712 7.13213213 7.14714715 7.16216216 7.17717718
7.19219219
7.20720721 7.22222222 7.23723724 7.25225225 7.26726727
7.28228228
7.2972973 7.31231231 7.32732733 7.34234234 7.35735736
7.37237237
7.38738739 7.4024024 7.41741742 7.43243243 7.44744745
7.46246246
7.47747748 7.49249249 7.50750751 7.52252252 7.53753754
7.55255255
7.56756757 7.58258258 7.5975976 7.61261261 7.62762763
7.64264264
7.65765766 7.67267267 7.68768769 7.7027027 7.71771772
7.73273273
7.74774775 7.76276276 7.77777778 7.79279279 7.80780781
7.82282282
7.83783784 7.85285285 7.86786787 7.88288288 7.8978979
7.91291291
7.92792793 7.94294294 7.95795796 7.97297297 7.98798799 8.003003
8.01801802 8.03303303 8.04804805 8.06306306 8.07807808
8.09309309
8.10810811 8.12312312 8.13813814 8.15315315 8.16816817
8.18318318
8.1981982 8.21321321 8.22822823 8.24324324 8.25825826
8.27327327
8.28828829 8.3033033 8.31831832 8.33333333 8.34834835
8.36336336
8.37837838 8.39339339 8.40840841 8.42342342 8.43843844
8.45345345
8.46846847 8.48348348 8.4984985 8.51351351 8.52852853
8.54354354
8.55855856 8.57357357 8.58858859 8.6036036 8.61861862
8.63363363
8.64864865 8.66366366 8.67867868 8.69369369 8.70870871
8.72372372
 8.73873874 8.75375375 8.76876877 8.78378378 8.7987988
8.81381381
8.82882883 8.84384384 8.85885886 8.87387387 8.88888889
8.9039039
8.91891892 8.93393393 8.94894895 8.96396396 8.97897898
8.99399399
9.00900901 9.02402402 9.03903904 9.05405405 9.06906907
9.08408408
9.0990991 9.11411411 9.12912913 9.14414414 9.15915916
9.17417417
9.18918919 9.2042042 9.21921922 9.23423423 9.24924925
9.26426426
9.27927928 9.29429429 9.30930931 9.32432432 9.33933934
9.35435435
9.36936937 9.38438438 9.3993994 9.41441441 9.42942943
9.44444444
9.45945946 9.47447447 9.48948949 9.5045045 9.51951952
9.53453453
9.54954955 9.56456456 9.57957958 9.59459459 9.60960961
9.62462462
9.63963964 9.65465465 9.66966967 9.68468468 9.6996997
9.71471471
9.72972973 9.74474474 9.75975976 9.77477477 9.78978979
9.8048048
9.81981982 9.83483483 9.84984985 9.86486486 9.87987988
9.89489489
9.90990991 9.92492492 9.93993994 9.95495495 9.96996997
9.98498498
10. 10.01501502 10.03003003 10.04504505 10.06006006
10.07507508
10.09009009 10.10510511 10.12012012 10.13513514 10.15015015
10.16516517
10.18018018 10.1951952 10.21021021 10.22522523 10.24024024
10.25525526
10.27027027 10.28528529 10.3003003 10.31531532 10.33033033
10.34534535
10.36036036 10.37537538 10.39039039 10.40540541 10.42042042
10.43543544
10.45045045 10.46546547 10.48048048 10.4954955 10.51051051
10.52552553
10.54054054 10.55555556 10.57057057 10.58558559 10.6006006
10.61561562
10.63063063 10.64564565 10.66066066 10.67567568 10.69069069
10.70570571
10.72072072 10.73573574 10.75075075 10.76576577 10.78078078
10.7957958
10.81081081 10.82582583 10.84084084 10.85585586 10.87087087
10.88588589
10.9009009 10.91591592 10.93093093 10.94594595 10.96096096
10.97597598
10.99099099 11.00600601 11.02102102 11.03603604 11.05105105
11.06606607
11.08108108 11.0960961 11.11111111 11.12612613 11.14114114
11.15615616
11.17117117 11.18618619 11.2012012 11.21621622 11.23123123
11.24624625
11.26126126 11.27627628 11.29129129 11.30630631 11.32132132
11.33633634
11.35135135 11.36636637 11.38138138 11.3963964 11.41141141
11.42642643
11.44144144 11.45645646 11.47147147 11.48648649 11.5015015
11.51651652
11.53153153 11.54654655 11.56156156 11.57657658 11.59159159
11.60660661
11.62162162 11.63663664 11.65165165 11.66666667 11.68168168
11.6966967
11.71171171 11.72672673 11.74174174 11.75675676 11.77177177
11.78678679
11.8018018 11.81681682 11.83183183 11.84684685 11.86186186
11.87687688
11.89189189 11.90690691 11.92192192 11.93693694 11.95195195
11.96696697
11.98198198 11.996997 12.01201201 12.02702703 12.04204204
12.05705706
12.07207207 12.08708709 12.1021021 12.11711712 12.13213213
12.14714715
12.16216216 12.17717718 12.19219219 12.20720721 12.22222222
12.23723724
12.25225225 12.26726727 12.28228228 12.2972973 12.31231231
12.32732733
12.34234234 12.35735736 12.37237237 12.38738739 12.4024024
12.41741742
12.43243243 12.44744745 12.46246246 12.47747748 12.49249249
12.50750751
12.52252252 12.53753754 12.55255255 12.56756757 12.58258258
12.5975976
12.61261261 12.62762763 12.64264264 12.65765766 12.67267267
12.68768769
12.7027027 12.71771772 12.73273273 12.74774775 12.76276276
12.77777778
12.79279279 12.80780781 12.82282282 12.83783784 12.85285285
12.86786787
12.88288288 12.8978979 12.91291291 12.92792793 12.94294294
12.95795796
12.97297297 12.98798799 13.003003 13.01801802 13.03303303
13.04804805
13.06306306 13.07807808 13.09309309 13.10810811 13.12312312
13.13813814
 13.15315315 13.16816817 13.18318318 13.1981982 13.21321321
13.22822823
13.24324324 13.25825826 13.27327327 13.28828829 13.3033033
13.31831832
13.33333333 13.34834835 13.36336336 13.37837838 13.39339339
13.40840841
13.42342342 13.43843844 13.45345345 13.46846847 13.48348348
13.4984985
13.51351351 13.52852853 13.54354354 13.55855856 13.57357357
13.58858859
13.6036036 13.61861862 13.63363363 13.64864865 13.66366366
13.67867868
13.69369369 13.70870871 13.72372372 13.73873874 13.75375375
13.76876877
13.78378378 13.7987988 13.81381381 13.82882883 13.84384384
13.85885886
13.87387387 13.88888889 13.9039039 13.91891892 13.93393393
13.94894895
13.96396396 13.97897898 13.99399399 14.00900901 14.02402402
14.03903904
14.05405405 14.06906907 14.08408408 14.0990991 14.11411411
14.12912913
14.14414414 14.15915916 14.17417417 14.18918919 14.2042042
14.21921922
14.23423423 14.24924925 14.26426426 14.27927928 14.29429429
14.30930931
14.32432432 14.33933934 14.35435435 14.36936937 14.38438438
14.3993994
14.41441441 14.42942943 14.44444444 14.45945946 14.47447447
14.48948949
14.5045045 14.51951952 14.53453453 14.54954955 14.56456456
14.57957958
14.59459459 14.60960961 14.62462462 14.63963964 14.65465465
14.66966967
14.68468468 14.6996997 14.71471471 14.72972973 14.74474474
14.75975976
14.77477477 14.78978979 14.8048048 14.81981982 14.83483483
14.84984985
14.86486486 14.87987988 14.89489489 14.90990991 14.92492492
14.93993994
14.95495495 14.96996997 14.98498498 15. ]
[[ 1.00000000e+00 1.00000000e+00 1.00000000e+00]
[ 9.84984985e-01 1.01501502e+00 9.69969970e-01]
[ 9.69744519e-01 1.02980458e+00 9.40841743e-01]
...
[-1.58071422e+00 -7.32077611e-02 6.28078545e-14]
[-1.57961500e+00 -9.69422089e-02 6.09217327e-14]
[-1.57815942e+00 -1.20660152e-01 5.90922513e-14]]

Text(0, 0.5, 'x')

 • The salient characteristic of linear system is oscilatory behaviour if the eigenvalues are
complex. Along the same line, we can infere that scalar systems can have only
exponential in time behavior (stable or unstable), while the second order system (matrix
A is second and higher order) can have oscillatory behavior. One of the special cases is
the situation when eigenvalues are purely complex (that is λ 1, 2=± j ω), which
characterizes pure oscillatory behavior.

Linear systems - Reactor stability

• Model based analysis is beneficial since the model parameters with well defined
physical meaning can take a role in the analysis of the dynamical system. Namely,
one can consider the following example of reaction engineering system with the
following reactor system dynamics. We consider the simple chemical reaction taking
place in the CSTR:
A→ pr o duc ts
and we are aware that this reaction is characterized by the heat release, hence the
material and energy balance can be written. It will be demonstrated that the
coupling of material and energy balances for the CSTR may give a rise to some
surprisingly complex and interesting behavior. We consider, the following balances:
d Ca C a f − Ca
=¿ − k Ca
dt θ
dT U T −T △ H R
=¿ ( T a− T )+ f − kC a
dt C ps θ Cp s
where the simple first order kinetics is taking place r a =k C a with the
heating/cooling coil at temperature T a. Hence, the following holds:
d Ca
=f 1 ( C a ,T )
dt
dT
=f 2 ( C a ,T )
dt
With reactor being at steady state, we are interested in finding out if small
perturbations in temperature or concentrations will induce the change to the new
steady state conditions or the system will be no sensitive to perturbations. We
linearize equations around the steady state of interest which is obtained by
calculating the following pair ( C a s , T s ):

0=f 1 ( C a s ,T s )
0=f 2 ( C a s ,T s )

so that the following steady states are obtained (with assumption that adiabatic
reactor is considered U =0) by solving the system of equations:
Caf
C a s=
1+k ( T s ) θ

Ca kθ
and (one can also look at the steady state conversions x a=1 − , or x a= )
Caf 1+k θ

T s=T f −θ
[ △H
Cps
k (T s )
Caf
1+ k ( T s ) θ )
− E (1 /T − 1/ T )
where k ( T s ) =k m e s
. One needs to solve for T s values. This can be done with
m

Newton-Raphson method. Hence, by the Taylor series expansion we obtain:

∂f1 ∂f1
f 1 ( C a , T ) =f 1 ( C a s ,T s ) + ( C a− C as)+ ( T − T s ) +h i g h e r − o r d e r t e r m s
∂ Ca ∂T
∂f2 ∂f
f 2 ( C a ,T ) =f 2 ( C a s ,T s ) + ( C a − C a s ) + 2 ( T −T s ) +h i g h e r − o r d e r t e r m s
∂ Ca ∂T

and we define the perturbation variables as C a p=C a − C a s and T p=T − T s. Therefore

the following approximate linear equations for the perturbation variables are
generated
d Ca p ∂f1 ∂f1
¿ Ca p+ T
dt ∂C a ∂T p
dTp ∂f2 ∂f2
¿ Ca p+ T
dt ∂C a ∂T p
 In the compact notation, with x p=[ C a p T p ) the above system becomes:

d xp
=J x p
dt
where J -Jacobian matrix is given as follows:

[ )
1 E
− −k − k Ca s
2
θ Ts
J=
△ HR U 1 △ HR E
− k − − − k Cas 2
C ps C ps θ C ps Ts

We consider the following parameters (parameters are taken from example 6.3.1.
steady-state multiplicity, Chemical Reactor Analysis and Design Fundamentals,
James Rawlings, John Ekerdt, Nob Hill Publishing, 2002. https://epdf.pub/chemical-
reactor-analysis-and-design-fundamentals.html):

Parameters Value Units

Tf 298 K
Tm 298 K
Cp 4.0 K j/k g K
Caf 2.0 k m o l/m
3

km 0.001 m in
−1

E 8.0×10 3 K
ρ 10
3
k g/m
3

△HR −3.0 × 10
5
kj/kmol
U 0

Let us look at the distribution of the eigenvalues versus the steady state temperature.

import numpy as np
import math
from numpy import linalg as LA
from numpy.linalg import matrix_rank
from scipy.optimize import fsolve

k_m=0.001#min^{-1}
Tm=298. #K
Tf=298.
E=8.0e3 #K
Caf=2.00
H=-3.e5
Cp=4.
rho_f=1.e3
Cps=rho_f*Cp
U=0

theta =np.array([1.79 ,15 ,20 ,25 ,30.9, 25 , 15 , 10, 2.0,

1.79, 10, 15 ,28])
T_guess =np.array([298.2,300.5,305,307,311 , 324, 338, 360,380,
421 , 440, 446.5,447. ])

def k(T,k_m):
E=8.0e3
Tm=298.
y=k_m*math.exp(-E*(1./T-1./Tm))
return y
T_s=np.zeros(len(theta))
k_T=np.zeros(len(theta))
Cas=np.zeros(len(theta))
J11=np.zeros(len(theta))
J12=np.zeros(len(theta))
J22=np.zeros(len(theta))
J21=np.zeros(len(theta))
lambda_1=np.zeros(len(theta))
lambda_2=np.zeros(len(theta))

##############################################
### Solving for steady state
##############################################
def f(T,theta):
k_m=0.001#min^{-1}
Tm=298. #K
Tf=298.
E=8.0e3 #K
Caf=2.00
H=-3.e5
Cp=4.
rho_f=1.e3
Cps=rho_f*Cp
k=k_m*math.exp(-E*(1./T-1./Tm))
f= Cps*(Tf-T)/theta-k*H*Caf/(1.+k*theta)
return f

#############################################
for i in range (0,len(theta)):
T_s[i]=fsolve(f,T_guess[i],theta[i])
k_T[i]=k(T_s[i],k_m)
Cas[i]=Caf/(1.+k_T[i]*theta[i])
J11[i]=(-1./theta[i])-k_T[i]
J12[i]=-k_T[i]*Cas[i]*E/pow(T_s[i],2.)
J21[i]=-H*k_T[i]/Cps
J22[i]=-(1./theta[i])-H*k_T[i]*Cas[i]/Cps*E/pow(T_s[i],2.)
A=np.matrix([[J11[i],J12[i]],[J21[i],J22[i]]])
 w, v = LA.eig(A)
#print w
lambda_1[i]=w[0]
lambda_2[i]=w[1]

%matplotlib inline
import matplotlib.pyplot as plt
plt.figure()
plt.xlabel('$T_{ss}$')
plt.ylabel('$\lambda_{1,2}$')
plt.ylim(-1,1)
plt.axhline(0, color='green', lw=2, alpha=0.5)
plt.plot(T_s,lambda_1,'-bo',T_s,lambda_2,'-bo')
plt.figure()
plt.ylabel('$T_{ss}$')
plt.xlabel('$\Theta$')
plt.plot(theta,T_s,'-r+')
plt.text(x = 1.97, y = 298,s = 'A',fontsize = 14, weight='bold',
backgroundcolor = 'white')
plt.text(x = 30, y = 320,s = 'B',fontsize = 14, weight='bold',
backgroundcolor = 'white')
plt.text(x = 2, y = 435,s = 'C',fontsize = 14, weight='bold',
backgroundcolor = 'white')
plt.text(x = 30, y = 440,s = 'D',fontsize = 14, weight='bold',
backgroundcolor = 'white')

Text(30, 440, 'D')

• Above figure shows the multiplicity of steady states when it comes to variation of
parameter θ . Namely, as it can be seen the unique value for θ induces three different
values for the T s s . Also, one can find that for different values of θ corresponding pairs of
steady state temperatures and concentrations induce matrix J which has the positive
eigenvalues. In other words, the corresponding steady states which belong to the AB-
branch are the stable ones, while the BC-branch steady states are unstable ones, and
finally the CD-branch is again the stable one. The points where the steady-state curve
turns are known as ignition (point B) and extinction (point A) points.
– Another argument to access the stability from the engineering point of view
is to consider the mass and energy balance and to substitute the expression
of the mass balance for the steady state in the steady state expression for the
energy balance, so that:

k C ps
0=−
1+ k θ
Caf △ HR+
θ
( T f −T )=Q g e n+ Qr e m

where the first term represent the heat generation rate Q g en while the second
term represents the enthalpy difference between the inflow and outflow
streams which is denoted as removal rate Qr e m. One can notice that the
steady state represents condition when these two rates are equal in
magnitude. Hence, we consider to plot these two functions as temperature
varies with the understanding that the first is the nonlinear function of
temperature by dependence on the rate constant, and the second term is
C
linear function with slope p s .
θ
import numpy as np

T=np.linspace(295.,450.,30)
theta=15
Q_gen=np.zeros(len(T))
Q_rem=np.zeros(len(T))

for i in range (0,len(T)):

Q_gen[i]=-k(T[i],k_m)/(1+k(T[i],k_m)*theta)*Caf*H
Q_rem[i]=Cps/theta*(T[i]-Tf)

%matplotlib inline
import matplotlib.pyplot as plt
plt.figure()
plt.xlabel('T(K)')
plt.ylabel('$Q_{rem},Q_{gen}$')
plt.plot(T,Q_gen,'b-',T,Q_rem,'r-')
plt.text(x = 320., y = 100,s = 'A - stable steady state',fontsize =
14, weight='bold', backgroundcolor = 'white')
plt.text(x = 350, y =9500,s = 'B - unstable steady state',fontsize =
14, weight='bold', backgroundcolor = 'white')
plt.text(x = 430, y =30000,s = 'C - stable steady state ',fontsize =
14, weight='bold', backgroundcolor = 'white')

Text(430, 30000, 'C - stable steady state ')

• As it can be seen from the figure, at the point A the slope of the heat removal rate is
higher than the slope of heat generation terms and therefore if the reactor temperature
is perturbed slightly by the feed it will induce the reactor to cool due to higher removal
rate and the temperature will move back to the stable point by system resisting to the
 temperature perturbation. The similar behavior can be observed in the case of C stable
steady state point where the temperature in the reactor will also resist to change to the
small perturbation by cooling off due to the higher rate removal. On the contrary, the
situation for the steady state point B reveals that the heat generation term has a higher
slope than the heat removal and any small perturbation of the temperature on the
branch where heat generation is higher than removal will be amplified and reactor
temperature will either escape from the point B up to the higher conversion and
temperature steady state or the perturbations can bring it to the lower stable steady
state. Finally, any small temperature perturbation around unstable steady state B will
lead to the transition to either upper or lower steady state.

Limit Cycles and sustained oscillations

• The dynamical behavior of the CSTR with the seemingly simple irreversible kinetics of
type A → p r o d u c t can be more complex than the analysis of stable and unstable steady
states. Namely, the operating parameters of the system might be of the type to support
the sustained oscillations or so called limit cycle behavior. We will consider the following
parameters:
Parameters Value Units
Tf 298 K
Tm 298 K
Cp 4.0 K j/k g K
Caf 2.0 k m o l/m
3

km 0.004 m in
−1

E 1.5×10 4 K
ρ 10
3
k g/m
3

△HR −2.2 ×10

5
kj/kmol
U 0
• It can be notices that these parameters are slightly different in the values of activation
energy E and constant k m.
import numpy as np
import math
from scipy.integrate import odeint

tspan=np.linspace(0,100,35000)
x=np.zeros(len(tspan))
T=np.zeros(len(tspan))
def k(T,k_m):
E=1.5e4
Tm=298.
y=k_m*math.exp(-E*(1./T-1./Tm))
return y
##############################################
### Integration subroutine
##############################################
def fun(t,X):
theta=35
Ca=X[0]
T =X[1]
k_m=0.004#min^{-1}
Tf=298.
Caf=2.00
H=-2.2e5
Cp=4.
rho_f=1.e3
Cps=rho_f*Cp
dCdt=(Caf-Ca)/theta-k(T,k_m)*Ca
dTdt=(Tf -T )/theta-H/Cps*k(T,k_m)*Ca
y=np.array([dCdt,dTdt])
return y
#############################################
init=[2.00,298.]

def ode45_step(fun, x, t, dt):

"""
One step of 4th Order Runge-Kutta method
"""
k = dt
k1 = k * fun(t, x)
k2 = k * fun(t + 0.5*k, x + 0.5*k1 )
k3 = k * fun(t + 0.5*k, x + 0.5*k2)
k4 = k * fun(t + dt, x + k3)
return x + 1/6. * (k1 + 2*k2 + 2*k3 + k4)

def ode45(fun, t, x0):

"""
4th Order Runge-Kutta method
"""
n = len(t)
x = np.zeros((n, len(x0)))
x[0] = x0
for i in range(n-1):
dt = t[i+1] - t[i]
x[i+1] = ode45_step(fun, x[i], t[i], dt)
return x

Xres=ode45(fun,tspan,init)

Caf=2.0
x=1.-Xres[:,0]/Caf
T=Xres[:,1]
#print(T)
#print(Xres[:,0])
%matplotlib inline
import matplotlib.pyplot as plt
fig=plt.figure()
ax1=fig.add_subplot(211)
ax1.plot(tspan[0:10000],T[0:10000],'b-')
ax1.set_ylabel('$T$(K)')
ax1.set_xlabel('time')

ax2=fig.add_subplot(212)
ax2.plot(tspan[0:10000],x[0:10000],'r-')
ax2.set_ylabel('$x$')
ax2.set_xlabel('$time$')

Text(0.5, 0, '$time$')

• In physical terms the oscillating behavior that we observe can be explained as follows:
Since the reaction is exothermic, the heat is liberated as the reaction proceeds. For the
proper combination of system parameters, this causes an increase in the reactor
temperature which increases the reaction rate constant, k, and therefore the rate at
which reactant is converted. This autocatalytic phenomenon proceeds at an accelerating
rate which becomes much faster than the rate at which reactant is supplied to the
system. Therefore, at some point the effective reactant concentration in the system
drops to a very small value and the reaction rate term, k V C a , becomes very small. At
the same time the rate of heat generation caused by the reaction must drop to a low
value. Then heat is removed from the system by convective flow and the cooling coil,
while the effective reactant concentration is building up. Although the decreasing
temperature tends to reduce the magnitude of the reaction rate term, the increasing
concentration increases its value. Then at some point the concentration becomes large
enough so that the reaction rate term becomes appreciable and the system again follows
 an autocatalytic path. Thus the composition and temperature in the reactor undergo
continuous oscillations, and the reactor outputs will be periodic even though the inputs
are maintained constant.

Limit Cycle example

• In order to demonstrate the limit cycle behavior, we considered the following
example given in Chemical Engineering Science publication (Unsteady state process
operation, J. M. DOUGLAS and D. W. T. RIPPIN, Chemical Engineering Science, 1966,
Vol. 21, pp. 305-315). We consider the slightly different system which account for
the mass and energy balance, so the dynamics is written as:
d z1 q
dt
¿
v
(1 − z1 )− k z1
d z2 q U a K qc ( z2 − zc )
¿ ( Zf − Z2) − + k z1
dt v V C p ( 1+ K q c )

T Cp
where z 1= A/ A f , ( A f is entry feed concentration of species A), z 2= ,
(− △ H ) Af
TcC p Tf Cp ρc
Z c= , zf = , K=2C c .
(− △ H ) Af (− △ H ) A f Ua

import numpy as np
import math
from scipy.integrate import odeint

V=1000.0 #
A_f=0.0065 #
Tf=350
Tc=408.6
E=28000 #K
R=1.987
del_H=27000.
rho=1.0
C=1.0
Ua=50
K=0.3
k0=math.exp(29.6298)
q=10
q_c=20

tspan=np.linspace(0,1000,1000)

z1=np.zeros(len(tspan))
z2=np.zeros(len(tspan))

def k_T(T,E):
R=1.987
y=math.exp(-E/(R*T)+29.6298)
 return y
##############################################
### Integration subroutine
##############################################
def fun(X):
z1=X[0]
z2=X[1]
V=1000.0 #
A_f=0.0065 #
Tf=350.
Tc=408.6
E=28000. #K
R=1.987
del_H=27000.
rho=1.0
C=1.0
Ua=50.
K=0.3
#K=2.*C*rho/Ua
#k0=math.exp(29.6298)
q=10.
q_c=5.
Zf=Tf*C*rho/del_H/A_f
Zc=Tc*C*rho/del_H/A_f
T=z2*del_H*A_f/(C*rho)
dz1dt=( 1.-z1)*q/V-k_T(T,E)*z1
dz2dt=(Zf-z2)*q/V-Ua*K*q_c*(z2-Zc)/V/C/rho/(1.+K*q_c)+k_T(T,E)*z1
y=np.array([dz1dt,dz2dt])
return y
#############################################

def ode45_step(fun, x, t, dt):

"""
One step of 4th Order Runge-Kutta method
"""
k = dt
k1 = k * fun(x)
k2 = k * fun(x+0.5*k1)
k3 = k * fun( x + 0.5*k2)
k4 = k * fun( x + k3)
return x + 1/6. * (k1 + 2*k2 + 2*k3 + k4)

def ode45(fun, t, x0):

"""
4th Order Runge-Kutta method
"""
n = len(t)
x = np.zeros((n, len(x0)))
x[0] = x0
 for i in range(n-1):
dt = t[i+1] - t[i]
x[i+1] = ode45_step(fun, x[i], t[i], dt)
return x

init=np.array([0.1,2.6])
Xres=ode45(fun,tspan,init)

init=np.array([0.1,1.4])
Xres1=ode45(fun,tspan,init)

x=Xres[:,0]
T=Xres[:,1]

#print(T)
#print(x)
%matplotlib inline
import matplotlib.pyplot as plt
fig=plt.figure()
ax1=fig.add_subplot(211)
ax1.plot(tspan,x,'b-')
ax1.set_ylabel('$x$')
ax1.set_xlabel('time')

ax2=ax1.twinx()
ax2.plot(tspan,T,'r-')
ax2.set_ylabel('$T$(K)')

ax3=fig.add_subplot(212)
ax3.plot(x,T,'b-',x[0],T[0],'rx',Xres1[:,0],Xres1[:,1],'b-',Xres1[0,0]
,Xres1[0,1],'rx')
ax3.set_ylabel('$z_1$')
ax3.set_xlabel('$z_2$')

Text(0.5, 0, '$z_2$')
Limit Cycle Example
• Let us consider a simple kinetic model of self-oscillation in glycolysis. We are lookin gin
the seminal paper of E.E. Selkov, Self-Oscillations in Glycolysis, European J. Biochem. 4
(1968) 79-86. The paper describes a simple kinetic model of an open monosubstrate
enzyme reaction with substrate inhibition and product activation. The model yields a
kinetic model of phosphofrucktokinase reaction, and represents a momosubstrate and
monoproduct reaction, with enzyme that is inhibited by the substrate and activated by
the product.
– Simple kinetic model of enzymatic catalysis with the product activation of the
enzyme is given as
γ γ
S1 + E S 2 ↔ S1 E S 2
γ γ
S1 E S 2 → E S 2 + S 2
γ
γ S 2+ E↔ E S2

where the substrate S1 (ATP) supplied by a certain source is irreversibly

converted to the product S2. The product is removed by an irreversible sink.
The free enzyme E (phosphofructokinase) is interactive by itself but becomes
active combining with γ product molecules to form the complex E Sγ2. We
introduce a notation that s1=[ S 1) , s2= [ S 2) , e=[ E ) and x 1=[ E S2 ), x 2=[ S1 E S 2 ).
γ γ

Therefore, according to the mathematical the law of mass action and the law
of mass conservation the model is:
 d s1
=v 1 − k +1 s s1 x 1+ k −1 x2
dt
d s2 γ
=k +2 x 2 − k +3 s 2 +k − 3 x 1 − k 2 s 2
dt
d x1
=− k +1 s1 x 1+ ( k −1 +k +2 ) x 2 +k +3 s 2 e − k −3 x 1
dt
d x2
=k +1 s 1 x 1 − ( k −1 +k +2 ) x 2
dt
de
=− k +3 s 2 e +k − 3 x 1
dt
– After some order reduction and system approximations the above system of
reactions becomes:
dx γ
=1 − x y
dt
dy
=α y ( x y −1 )
γ −1
dt
where at γ >1 system represents a mathematical model of a product-
activated and substrate-inhibited reaction (Remark: The above system is a
generalization of Lotka system and coincides with it for γ =1).

import numpy as np
import math
from scipy.integrate import odeint

alpha=1.1 #
gamma=2 #

tspan=np.linspace(0,20,30)

x=np.zeros(len(tspan))
y=np.zeros(len(tspan))

##############################################
### Integration subroutine
##############################################
def fun(X,t):
x=X[0]
y=X[1]
alpha=1.1
gamma=2
dxdt=1.-x*pow(y,gamma)
dydt=alpha*y*(x*pow(y,gamma-1)-1)
 y=np.array([dxdt,dydt])
return y
#############################################

init=np.array([0.1,2.6])
Xres=odeint(fun,init,tspan)

%matplotlib inline
import matplotlib.pyplot as plt
fig=plt.figure()
ax1=fig.add_subplot(211)
ax1.plot(tspan,Xres[:,0],'b-')
ax1.set_ylabel('x')
ax1.set_xlabel('time')

ax2=ax1.twinx()
ax2.plot(tspan,Xres[:,1],'r-')
ax2.set_ylabel('y')

ax3=fig.add_subplot(212)
ax3.plot(Xres[:,0],Xres[:,1],'b-',Xres[0,0],Xres[0,1],'rx')
ax3.set_ylabel('$x$')
ax3.set_xlabel('$y$')

Text(0.5, 0, '$y$')
Final Note on Nonlinear and Linear Systems
• Linear systems (e.g. models of type A , B ,C , D ) analysis is well known and established
engineering material.
– We will realize that the liner system properties (stability, controllability and
observability) can be easily established and quantified.
– The link between local linear analysis and full nonlinear dynamical description is
established.
– We consider different types of nonlinear behaviour among each the limit cycle is
the most transparent and emerges from the simple model.

Summary
• We learned how to starting from the general nonlinear equations obtain linear
equations.
– We learned characteristics of linear models and we explore numerical solutions
of the linear models