- Email: [email protected]
Sensitivity analysis of oscillatory systems
Sensitivity oscillatory
analysis sys terns
of
Mark A. Kramer* Department of Chemical 08544, USA
Engineering,
Princeton
University,
Princeton,
New Jersey
Princeton
University,
Princeton,
New Jersey
Herschel Rabitz Department
of Chemistry,
08544,
USA
Joseph M. Calo Division of Engineering, Brown University, Providence, (Received March 1983; revised August 1983}
Rhode
Island 02912,
USA
Methods of sensitivity analysis are extended to find the parametric dependencies of systems of ordinary differential equations which exhibit limit cycle oscillations. The quantitative relations between the system parameters and the observable period, amplitude, phase and cycle shape are developed. These formulae, presented for both the first and second order, are applicable to systems of arbitrary size and complexity. The techniques are used here to develop correlations for period and amplitude in a non-isothermal oscillating stirred-tank chemical reactor, and to find and optimize a beneficial periodic operating strategy for a lumpedparameter catalytic reacting system. Key words: mathematical oscillatory systems
Sensitivity analysis, the systematic study of the effects of parameter values on the predictions of mathematical models, is a useful tool for model validation and evaluation, as well as for quantifying the effects of parametric uncertainty and variability. The numerical link between parameters and predictions in systems of ordinary differential equations (ODES) is provided by the linear sensitivity coefficients’ aui(t)/aoll, where yi is the ith dependent variable, CX~ is the ith system parameter, and t (time) is the variable of integration. These coefficients may be obtained by several numerical techniques ‘-’ if the model is not analytically solvable. Higher order sensitivities, such as a2yi(t)/$ a&, are also calculable, and provide information on nonlinearity and parametric interactions. The interpretation and application of sensitivity information is a rapidly-developing area.6 However, the available analysis is not directly applicable to oscillatory systems. *Present address: chusetts Institute USA.
328
Department of Chemical Engineering, Massaof Technology, Cambridge, Massachusetts 02139,
Appl. Math. Modelling,
1984,
Vol.
8, October
model,
sensitivity
analysis,
limit
cycles,
Furthermore, quantities of obvious interest for oscillatory trajectories, such as frequency and amplitude, are not directly addressed by the sensitivity coefficients. These complications will be shown to arise from the parametric dependencies of the velocity of the trajectory around the limit cycle. A method of analysis will be presented which isolates the parametric sensitivity of the cycle phase, thereby rendering explicit the parametric effects on period, amplitude, shape, size, and dynamics of the oscillation, and the associated transient. The procedures developed in this work are applicable only to trajectories which evolve to limit cycles, that is, to systems in which the equilibrium orbits are independent of the initial conditions. This type of behaviour arises in many problems in biology, catalytic kinetics, and reactor dynamics. For non-limit cycle oscillations, several of the formulae presented here are valid - in particular, equations (lo)-(1 5). However, it is impossible to define the sensitivity of the phase of the oscillation, as is done below, which is crucial to subsequent developments. 0307-904X/84/05328-13/$03.00 Q 1984 Butterworth & Co. (Publishers)
Ltd
Sensitivity
Parametric studies of limit cycles have usually focused on determination of parametric effects ‘in-the-large’, that is, determination of the ranges of parameter values for which a particular system may have oscillatory solutions. In contrast, relatively little effort has been directed at ‘finegrained’ analysis, in which the effect of parameter perturbations on specific equilibrium orbits is investigated. This information can be useful in fitting a model to data, in uncertainty analysis of modelling predictions, in experimental design, in identification of governing or underlying processes, or in optimization of periodically controlled processes. These applications, excepting the last one, are discussed in detail in Demiralp and Rabitz.’ An example of periodic optimization is included in this paper. The techniques presented here are applicable to any system described by ordinary differential equations. Previous efforts to describe limit cycle behaviour relied on assumptions of near-linearity and near-harmonic oscillations,7-g or used singular perturbation techniques which become inaccurate for large values of the perturbation variable.1o~‘1 Numerical comparisons between these methods and sensitivity analysis will be presented.
Sensitivity equations The system under consideration is a set of coupled ODES, which may be written in the following general form:
(1) Here, y is the dependent variable vector (y E R”), a is the parameter vector (ar E Rm), and f is the variable of integration. The parameter vector may include the initial conditions y(O), although y(0) usually does not enter explicitly into f. The latter vector is generally a nonlinear function of y and the non-initial value parameters. This equation can be differentiated with respect to the jth parameter CY~ to yield’:
d a~ dt
aai
= J(t)
2 (t) + E(r) J
(2)
J
where J is the n x IZ Jacobian matrix (J&t) = afi(t)&& and ay(t)/aol, is an n-vector of linear sensitivities. The initial condition for ay/aol, in equation (2) is the zero vector, unless aj is the initial condition for yi (i.e. CX~ f yi(O), in which case ayi(O)/aai = 1). Equations (1) and (2) can be integrated simultaneously, with an appropriate differential equation solver, to yield the linear sensitivities. Alternatively, an n x n Green’s function matrix, K(t, t’), may be constructed by solving3: ; K(t, t’) - J(t) K(t, t’) = 0
K(t’, t’) = I, I 2 t’
analysis of oscillatory
The relevant equation
systems: M. A. Kramer et al.
for the adjoint is simply:
K+(t’, t) + K+(t’, t) J(i) = 0
5
K+(t, t) = I, t’ G t
(5)
Equations (3)-(S) are exact, and yield the same ‘raw’ sensitivity values as equation (2). These ‘raw’ coefficients will be manipulated to yield information pertinent in the limit cycle case.
Sensitivity of phase to initial state The system (1) will be assumed to have a limit cycle solution, i.e. the steady-state orbit is independent of the initial conditions. We will first consider sensitivities of the type [avi(t)/ayi(0)]., which represent the change inYj(t) induced by a differential charge in yi(O), all other yk(O) and all parameters held fixed. The system state at t = 0 may lie anywhere within the range of attraction of the limit cycle. Since the steady-state orbit is unaffected by such a perturbation, only the position along the orbit (phase) can be influenced by yi(O), after an initial transient period. This phase offset can be quantified as the coefficient at/ayi(0), where at is the ultimate differential delay or advance along the limit cycle resulting from the initial differential perturbation ayi(0), measured in time units. This is shown schematically in Figure 1. The formula for calculating this sensitivity is easily derived by taking the total differential ofyi(or, t’, y(O)), then dividing by dyi(O) with dvi = 0 and at constant cr and yk(O), k #j. The imposition that dyi = 0 follows from the assumption that the limit cycle is ultimately achieved. Solving for the phase sensitivity, we obtain: Qj(0)
at
= -=-r!h+n[$)/f~)
(6)
auj (0)
The limit appearing in this formula is necessitated by the fact that the trajectory never actually reaches the limit cycle, but only approaches it asymptotically. Numerically, the number of revolutions required to converge to an accurate value of Q is governed by the stiffness of the system.” Equation (6) iS independent of the choice of i, since the
STEADY - STATE
#ORBIT
(3) DIFFERENTIAL
The linear sensitivities may then be expressed in integral, rather than differential form: r $
(t) = K(t, 0) 2 J
(0) + j- K(t, t’) a: J
0
(t’) dt’
(4)
J
Since t’ G t, it is actually more practical to solve for the adjoint Green’s function, K+(t’, t), subsequently applying the equality K+(t’, t) = K(t, t’) for use in equation (4).7
Figure 1 Schematic representation in two dimensions of the sensitivity of the ultimate cycle phase to initial conditions
Appl.
Math.
Modelling,
1984,
Vol. 8, October
329
Sensitivity
analysis of oscillatory
systems:
M. A. Kramer
et al.
cycle phase is a property of the entire system, rather than single variables. One further property to be observed here is that @i(t)/ayi(0) must pass through zero whenever dyj(t)/dt is zero, since the phase sensitivity must be finite. Special numerical attention (perhaps application of L’Hospital’s rule) may be needed at these points. Sensitivity of the ultimate cycle phase to perturbations of the dependent variables (y) is of particular importance when the initial system state is numerically ‘on’ the steadystate orbit. These on-cycle phase sensitivities are closely related to the Lyaponov stability of the cycle, and are, in fact, linearly related to the stability eigenvectors.13 These values will be used subsequently to calculate the sensitivity of the system phase to non-initial condition parameters. For this application, initial condition phase sensitivities are required at each point the parametric phase sensitivities are to be calculated. If we were to use equation (3) to obtain the initial condition phase sensitivities, as suggested by equation (6) (since Kii(f, 0) = avi(t)/vj(O)), a new solution of the matrix differential equation (3) would be necessary, using each point of interest in turn as the zero time in equation (3), each time the parametric phase sensitivities were to be calculated. Instead, it is preferable to use equation (5) - the adjoint Green’s function - which need only be solved once to obtain the entire set of on-cycle phase sensitivities. The solution of (5) proceeds backwards along the trajectory beginning at an arbitrary point t’ on the steady-state cycle, and produces Kt(t, t’), where t’> t. The phase sensitivity for a perturbation at t is then calculated from:
Qdt>= ,!e_ - K&(t,t’)
I
dvi@‘I dt
(7)
Note that only one row of Kt is required, so equation (5) can be solved as a vector rather than a matrix, thereby saving considerable computational effort. As a matter of practicality, the limit appearing in equation (7) needs only to be observed to the extent that a problem-dependent transient period is passed. This transient has been found, in several stiff systems, to be somewhat shorter than one period. Specific times may be checked to ensure they lie far enough from t’ for equation (7) to be applicable. The following relationship must hold for the calculated on-cycle phase sensitivities: f
Qk(t)
y
= -
(8)
1
k=l
The relationship is derived in the following way. Take the total differential of yi(t’) = yi(t’, y(t)), and divide by dt with t’ fixed (0~is also fixed for the purposes of this derivation). Thus:
The phase sensitivities can be used to identify the directions in which the system may be perturbed without disturbing the phase. A perturbation dy applied at time t leaves the phase of the steady state orbit unchanged if and only if: n
c
&c(t) dYk= 0
k=l
We will define @o(t)to be the collection of all vectors y’(t) evolving to the same orbital phase as y(t) at long times. Any differential dy satisfying the previous summation lies in 4(t). A perturbation not lying in g(t) changes the system phase by the amount dt, where y(t) + dy lies in (((t + dt). Mapping dy back on to G(t) removes phase shift effects and leaves a differential dy’ which satisfies the previous summation. This mapping process will be used subsequently to isolate the pathindependent parts of parametric sensitivities. The dynamics of response to differential perturbations in the dependent variables are contained in the transient period of K. To extract these dynamics, the phase shift component Qj(O)(dy(t)/dt), is subtracted from the ith column of K(t, 0). This leaves the component of the perturbation along the constant-phase locus g(t). The norm of this remainder is the distance to the nominal trajectory along G(t). The damping of the original perturbation may be examined by preparing a plot of this distance against time. The function must asymptotically approach zero for stable limit cycles. A word should be said at this point about sustained oscillations which are not limit cycles. The initial condition sensitivities for these oscillations do not demonstrate the simple behaviour of equation (6). Instead, the sensitivities contain information both on the phase sensitivities and on the sensitivity of the shape or location of the steady-state cycle. Unfortunately, these two parts cannot be separated unambiguously (although equation (14) can be used to determine the period sensitivity), since the current definition of initial condition phase sensitivities requires that the nominal and differentially-perturbed trajectories converge to the same path, so one trajectory is strictly leading or following the other. If the two trajectories follow different paths, there is no obvious way of judging which is ‘ahead’ and which ‘behind’. Work is continuing on this problem, since it would be useful to generalize the analysis presented here to non-limit cycle oscillations, and in fact, to altogether non-oscillatory trajectories. In the latter case, it is hoped that conventional sensitivities could be divided into two parts, representing the sensitivity of the velocity of the trajectory and the sensitivity of the shape of the trajectory, which would be analogous to the current treatment of the limit cycle sensitivities.
Sensitivity d@) -= dt
’ aYi (t’) dYk(t) c--
k=
1 ayk(t>
dt
Dividing both sides by dyj(t’)/dt, and applying equation (7), we obtain equation (8). If equation (8) is violated, either t is too close to t’, or the numerical accuracy of the adjoint kernel, and hence Qk, is poor. To get Qk(t) for t’- t less than the transient time, equation (5) must be integrated around the entire cycle (backwards) to time t - f. Then &(t) = &(t - r), since t and t - T refer to the same position on the cycle.
330
Appl.
Math.
Modelling,
1984,
Vol. 8, October
to non-initial
condition
parameters
Unlike initial conditions within the attractor region, parameters appearing explicitly in equation (1) can affect the period, amplitude, shape, and location of the limit cycle. This leads to additional difficulties in the interpretation of the raw sensitivity coefficients, since the sensitivity of these features, as well as the phase-shift information, are together contained in the raw coefficients. In order to deconvolute this information, one must separate the parts of the raw sensitivities which result from a parametric phase dependence, and those which refer only to the shape and position of the trajectory in phase space. The former will be
Sensitivity
referred to as the path-dependent part of the raw sensitivity, and the latter the position-dependent (or pathindependent) parts. As in the case of initial condition sensitivities, the contribution to the raw sensitivity due to the phase lead or lag can be written as the product of a differential phase shift, (at/sol,),, and the local slope,13a so the total sensitivity may be written:
analysis
of oscillatory
systems:
M. A. Kramer
et al.
In equation (13b), (dyi/daj), is the derivative of equation (10) taken with r (but not y(O))-fixed. Equation (11) may therefore be written:
t
a7 dJJi
7
aai dt
(14aJ
which, because (dyk(O)/dQj), = dYk(O)/dOlj, is equivalent to:
(9) where (ayi/&xj)@ refers to the path-independent portion of the tbtal sensitivity. The formula for calculating (at,/a~~)G will be developed shortly, and the meaning of the subscript $ will be discussed, but first it will be demonstrated that division of the raw sensitivities, as indicated in equation (9), is a natural consequence of partial differentiation of y with respect to the system parameters. If the system (1) has a periodic solution, yi(t) can be represented as a Fourier series: m
yi(t)
=
2nd Qi, cos -
C
2nnt + b’,
7
n=O
(10)
sin 7
where 7 is the period of the oscillation. The coefficients af, bi, and 7 are implicit functions of the system parameters. Differentiating equation (10) with respect o+, the linear sensitivity coefficient can be expressed as: dyi(t) dai
ar -
2nt
1 (mzI, sin ?Z aiinzo 7
=72 -
aat Csol-I
+n:o
2nW cos -
7
+
abi
2nN
-
-=
which, in general, is not zero. The sensitivity as previously defined is the partial derivative of yi with respect to CY~ at constant y(O), i.e. (ayi(t)/&q),(,), which is zero at t = 0. TO correct equation (1 l), consider Yi = yi [CY,t, yo(~l)] , which is the functional form implicit in equation (10). Taking the total differential, and dividing by daj, with dt = 0 and all other parameters fixed, we obtain:
dori
am
n
dYi(t)
CPP
k=l
aYk(o)
d&O)
(12)
d&i
The final quantity here cpmes from equation (1 l), evaluated at t = 0. Before applying this expression, note that: = sio
[-nab
sin (z)
2
Qk(t’)
7)
$
k=l
and that: [$cosfT)+$sin(y)] (13b)
(16) I
7 (13a)
(gi=,Z,
+ a0ij
Equation (15) is valid for all t, if equation (10) is satisfied. The summation term in equation (15) damps to zero for large t, and when this term can be ignored, only the sensitivities ayJ?Ior, are required to calculate the sensitivity of the period. Edelson14 used this approach, in analysis of the ‘Oregonator’, a prototype oscillator. The sensitivity of the phase to non-initial condition parameters will now be derived. To begin, the role of y and t as dependent and independent variables, respectively, must be reversed. Along any phase-space trajectory, a point may be described either by its coordinates y, or the time r at which the point is reached. There is, in fact, a one-to-one correspondence between y and t along a given path because phase-space trajectories never cross or coalesce, provided that the initial conditions lie off the limit cycle. Therefore, along a given path, time may be written as a function of position (i.e. t = t(y(a)). The parametric sensitivity of this function may be probed. Taking the total differential and then the partial differential with respect to aj:
+ nb’, cos (El] 7
avi(f
(15)
=
z
-
aCYi r 1 (11)
Note that the total differential symbol is used in this expression; indeed, this dyi(t)/daj is not the sensitivity coefficient one calculates using equations (2) or (4). This can be seen by evaluating equation (11) at t = 0, which gives: m 2-i dyj(0)/daj = c ““” n=O aai
@i(t) -=_-sol,
Equation (14b) was also given by Tomovic’ and !Larter,13 and is also known as a result of multiple scale analysis in perturbation theory. There is a general similarity in the appearance of equations (14b) and (9), however, (ay/acu), and (ay/acu)G are not equivalent, and (at/&x), # t/@-/aa). It will be shown subsequently that equation (14b) is deficient because only the overall phase shift behaviour is taken into account. Nonetheless, it can be seen from this derivation that application of the chain rule in partial differentiation of y quite naturally gives rise to the division of the raw sensitivity into two parts. Equation (14a), incidentally, suggests a method of isolating the period sensitivity. Since (dyi(t)/daj)T is periodic (see equation (13b)), the period sensitivity can be calculated from the sensitivities for adjoining periods”:
a7 __hi(t) sol, i\ aai
- nb; cos ““‘) 7
-sin
(14b)
The sensitivity (at(t’)/&+)+ is a measure of the phase lead or lag, in time units, at time t’, resulting from the differential kj. The subscript C$is included because this phase sehsitivity gives the phase shift in reference to the constarit-phase locus @(t’). It will be recalled that @(t’) is the set of all phase-space points, which, considered as initial conditions, ultimately give rise to the same phase on the
Appl.
Math.
Modelling,
1984,
Vol. 8, October
331
Sensitivity analysis of oscillatory systems: M. A. Kramer et al.
mentary sensitivity coefficients. If the initial conditions are numerically ‘on’ the nominal steady-state orbit, the following formula is valid:
+aCY y2
Figure 2 Phase shift due to perturbation in the non-initial condition parameter 0~. The phase difference indicated is at = t’ - r”
steady-state orbit as the nominal solution (evolving from y = y(t’) at t = t’). Considering a nominal and a parametrically perturbed trajectory, if the perturbed trajectory crosses tit’) at t = t”, then the phase shift given by equation (16) is (at/da), = t’- t” (see Figure 2). An alternative derivation of equation (16) is as follows. The parametric phase sensitivity at t = t’ can be found by (a), solving the model using the parameter value oj + doj, from time zero to t = t’, and (b), returning the parameter to its nominal value, and solving the model until the trajectory closely approaches the steady-state cycle at t = t”‘. The phase sensitivity is then given by comparison with the nominal solution of the model at this time, using equation (6). After step (a), the difference between the nominal and perturbed trajectories at t = t” is dy = (ay(t’)/aaj) daj. In step (b), dy is considered a non-parametric perturbation from the nominal traje’ctory, and the ultimate phase shift arising from dy is calculated from:
This is surprising in view of the fact that, even when initial coditions are on-cycle, there is a transient in ay/aaj, prior to establishment of regular behaviour, as given by equation (15). There is, however, a structure to the components of ay/acri during this transient, which makes at/iJol, regular from time zero (i.e. equation (17a) valid for t’ -+ O+). A geometric interpretation of equation (16) can be used to explain this effect. Equation (16) maps the perturbed trajectory back to the nominal trajectory along a line of constant phase $, and then measures the phase shift. Therefore, the same phase sensitivity will be calculated whether the perturbed trajectory lies between the nominal limit cycle and the perturbed limit cycle (corresponding to the transient in ay/aoj), or the perturbed trajectory lies on the perturbed steady-state cycle (corresponding to regular behaviour of ay/&rj, indicated by complete damping of the summation term in equation (15)). Thus, the transient in the linear sensitivities has no effect on the phase sensitivities. We may now return to the comparison of equations (9) and (14b). It should be clear that the period sensitivity is only one aspect of the total parametric phase sensitivity. The period sensitivity is the overall period-to-period phase shift, while the parametric phase sensitivity (at/aaj) encompasses not only this overall behaviour, but also the local variations in the phase lead or lag. These local variations, being functions of position rather than path, can be described by a periodic function (at/sol,),:
(18) if the trajectory is effectively on the limit cycle. Separation of the total time scale sensitivity into these two components allows calculation of at/sol, at any time once a single cycle has been analysed. The total sensitivity is then used in equation (9) to determine the path-independent portion of the raw parametrjc sensitivity. It can be shown that the path-independent sensitivities (ay(t)/aor,), do indeed lie on the constant-phase locus tit), justifying the subscript 4 in equation (16). Beginning with equations (9) and (16) and applying equation (8):
-f
ti?k(6
(y)@
= ’
(19)
k=l
dt = f
Q,&‘) d2’k
k=l
which follows from the definition of Q. Substituting for dyk in this expression, and dividing by doj, gives equation (16). A formula for the sensitivity of the period of the oscillation can now be derived. This sensitivity is simply the accumulated phase lead or lag after completion of one revolution around the steady state orbit, SO that:
(17a) where t’ is any time longer than the transient
332
Appl.
Math.
Modelling,
1984,
Vol.
in the ele-
8, October
These sensitivities therefore conform to the definition of g(t). The dynamics of response to step changes in the parameters, when the trajectory is on cycle, is contained in the transient period of (ay/acu)Q The approach to the perturbed equilibrium orbit is given by the difference (ay(t + 7)/h), - (ay(t)/aa),. This will be demonstrated in the example below. The sensitivity of amplitude to parameters is now easily determined. Since at an extreme in yi, dyi/dt = 0, the path-independent sensitivities are equal to the raw coefficients at these points. The sensitivity of the amplitude to first order is therefore simply the difference of the raw coefficients at the extrema.
Sensitivity
analysis of oscillatory
systems:
M. A. Kramer
et al.
Formulae for the second order Second order sensitivity information permits analysis of parametric interaction, and/or increased accuracy in predictions of period, amplitude, or cycle shape. These coefficients may be obtained by direct solution of the governing equations (equations (21) and (24)) or by using a packaged sensitivity solver such as AIM,4P5which uses a modification of the Green’s function method.3 The AIM algorithm is particularly useful for second order calculations, since the integral tranformation of the Green’s function method reduces the computational cost of the second order coefficients to the same order as first order coefficients. Direct solution of equations (21) and (24) will, in general, be much more expensive. Second order sensitivity coefficients are treated in a similar fashion to the first order coefficients. Differentiation of equation (6) with respect to a second initial condition, yk(0), yields the following expression for second or phase sensitivities:
+;
(t,
0)
=
J(t)
kjk(t,
+
L(t)
Kj(t,
= aJ,,/ay,,
’
foli~j(‘)
sjk(o) = 0
aJ&)/aqW
(fajak@))i = a2fi(t)/aaj(t) a&k(t)
0)
Kk(t,
0)
and
and Sj(t) = ay(t)/aaj Equation (23) may be used to obtain an explicit expression for the second order period sensitivity. As in the first order case, this is done by comparing the sensitivities obtained on two successive periods. Noting that:
-=
a27
a&j aark
a27
_= aaj ask One may also obtain second order phase shifts directly for all points around the orbit by solving a single differential equation. This is similar to the adjoint procedure (equations (5) and (7)) used to obtain the full set of on-cycle first order phase sensitivities. The relevant equation, from differentiation of equation (5) and application of equation (7), is:
qjZ(t)Jdt)
a2t(t’ + 7) -- a?(i) a&jask
affj aak
where t’ is any time referring to a numerically position, we obtain:
It will be observed that, due to symmetry:
‘y =,ilQdt)L,,(t) -
+ J,(t) Sk(t)
(24)
(Jaj(Ohp =
(21)
where L,,
= J(t) S;k@
+ J&k(f) sj(t> + L(t) sj(t) Sk(t)
ia%)QjW) 0)
f$k(t)
;
aYk@)
for any species i. The second order initial condition sensitivities (kjk(t, O))i = a2yi(t)/dyi(0) ayk(O) are defined by: *jk 7
The derivation of this expression is given in the Appendix. Again, the subscript @refers to the path-independent portion of the total sensitivity. As in the first order case, it can be shown that the (ay(t)/hx, acq)@lie on N(t), where (t(t) now includes the curvature due to second order effects. The r’aw second order barametric derivatives Sik(t) = ay(t)/acq sol, are defined by the equation3:
where
a2Yi(t) aYj@)
(23)
(22)
where qjk(t) = a2t/@j(t) r&(t), and q(0) = 0. The solution of this equation proceeds backwards around the limit cycle. The starting point may be anywhere numerically ‘on’ the steady-state orbit, and, as in the first order case, there is a transient period associated with equation (22). This transient corresponds to the time required for the damping of the perturbations ayj(t) and ayk(t). For non-initial condition parameters, a decomposition formula may be obtained from differentiation of equation (9) yielding:
‘dyi -r a2yi(t’) F( dt 1 aaj hk
a2yi(t’ ‘f T) a&j a&
+ a7 a7 d2Yi _ 37 d affj aak, dt2
---
‘on’-cycle
Mt’)
&vj dt i ask i
-
(25)
where i is arbitrary. The derivation of this expression is also given in the Appendix. To evaluate equation (23), an expression for (a2t/ &q a&& must be developed. Differentiating equation (163, we obtain:
(26) Unfortunately, an additional differential equation must be solved to obtain’aQ/aak, since from equation (7) it can be seen that the explicit expression for a&(t)/&, Will involve the term:
APPI.
Math.
Modelling,
1984,
Vol. 8, October
333
Sensitivity analysis of oscillatory systems: M. A. Kramer et al.
where i is arbitrary and t’ is a time sufficiently removed from t such that the transient in equation (7) has run its course. Because f appears as a reference time in this term, it cannot be calculated from any of the raw sensitivities. The differential equation directly controlling aQl/acu, can be derived from partial differentiation of equation (5) with respect to 0~~:
+
aQm(t) ~ J,nl(t)
(27)
aLwk
Solution of this equation, like that of equations (5) and (22), proceeds in the reverse direction around the limit cycle, and provides the full set of on-cycle derivatives. Equation (26) can then be used to calculate the second order phase sensitivity to non-initial parameters. As in the first order case, this phase sensitivity can be compared on successive periods, in a manner analogous to equation (18), to calculate the second order period sensitivities. If the initial conditions are numerically ‘on’ the limit cycle, the second order period sensitivities are simply -(at(T)/ acvi aaj)e To determine the second order amplitude sensitivities, the effect of the parameters on the times of occurrence of the extrema must first be calculated. Since dyi(ol, r,)/dt = 0, where t, is the time of occurrence of the extrema, we obtain:
(28) by differentiation with, respect to oj, and solution for at,/&.xj. We wish now to find the component of a’yi/ sol, aa, in the direction of change of the extrema. This is given by (see Appendix): _
a2Yi(fe))
ateat,d2Yi(te>
aaj ih!,
a@j acUk dt2
(29)
The amplitude sensitivity is then the difference between the second order sensitivities in the direction of change of the extrema: -=
a2ai
a2Yi(fmin)
a’Yi(LJ
a@j aa,
aaj ibk t--
at,,, ihi
%,,
h!j a&!k d2Yi(fmax)
aaj
af,,,i, - afmin -aaj aa,
dt2 d’Yi(trnin) dt2
(30)
The procedure, then, for a full second order analysis is as follows. First, the raw coefficients are obtained by solving equation (24) or an equivalent method. No more than two cycles need be traversed to provide enough information for all the analyses. The second order period sensitivities may be obtained immediately via equation (25). Solution of equation (27) is then required to provide the necessary data for calculation of the time scale sensitivity by equation (26). The full decomposition is then executed using equation (23). The path-independent sensitivities thus isolated are used for estimating the parametric effects on the shape of the equilibrium orbit. Second order
334
Appl.
Math.
Modelling,
1984,
Vol. 8, October
phase sensitivities may be obtained by solving equation (22). Second order amplitude sensitivities are calculated from equations (28) and (30).
Examples Prediction of amplitude, period, and dynamic behaviour of an oscillating CSTR
To demonstrate some of the applications of the preceding techniques, a simple perfectly-mixed, lumped parameter chemical reactor model will be used. Although this model has been thoroughly studied in terms of relations between parameters and stability, the relations between specific limit cycle trajectories and parameters have not been systematically analysed. In this regard, sensitivity analysis offers some new insights. The specific system under consideration is a well-mixed continuous non-adiabatic stirred tank reactor (CSTR) in which a single first order, irreversible exothermic reaction occurs. This system is represented by the following two coupled differential equations,15 representative of the dimensionless transient mass and energy balances for the reactor: j1 = -Y1
j,=b
+ Da(l -yJ
-py2
Y2
exp 1+ fY2
+ Da(1 -ye)
(31) Y2
exp 1+ EY2
where Da = us exp [---u&p + e~~)]/( 1 - ~3. In this equation, b is the dimensionless enthalpy of reaction, E is the dimensionless activation energy, U, is the dimensionless concentration at the system steady state for a stable system, or at the limit cycle focus for unstable systems, and p is the dimensionless heat transfer parameter. The nominal parameter values b = 20, E = 0.010375, us = 0.7 and p = 0.417 give rise to an unstable condition, indicated by det J = 14.4, tr J = 1.8, at the point at whichjl =i2 = 0. (The conditions det J > 0 and tr J < 0 are necessary and sufficient for the point to be locally asymptotically stable in this particular system).16 The trajectory for the nominal parameter values is shown in Figure 3. An arbitrary point on the cycle y = (0.8851, 1.678) was chosen as the reference point for the sensitivity analysis (i.e. t = 0 at this point). Beginning at the reference point, the adjoint kernel equation (5) was solved by coupling it with equation (3 1). The resulting values of K 7 were divided by dy(O)/dt according to equation (7) to produce the initial condition phase sensitivities (see Figure 4). The definition of the constant-phase locus 4(t) (zQk dyk = 0) was then used to determine the approximate course of these loci in the immediate vicinity of the orbit. The further requirement that all @((t)meet at the unstable focus of the limit cycle, which results from the fact that a randomly directed perturbation away from this point can lead to any on-cycle phase, allows the position of the constant-phase loci to be estimated for the entire cycle interior. The results are shown in Figure 3. The raw coefficients were then obtained using equation (2). The sensitivities with respect to the parameter p are shown in Figure 5. The strong peaks in these values result from the ‘saw-tooth’ (relaxation oscillation) character of the individual trajectories. The non-periodicity predicted by equation (9) is easily observed. If one naively uses these values in a first order Taylor series expansion (y(/* + AP) =
Sensitivity
1
6.0
analysis of oscillatory
systems:
M. A. Kramer
et al.
60C
7.0
6.0
YI figure 3
Steady state orbit for CSTR. @(t) are shown every l/10 period
Lines of constant phase
I
I 3.0 I
I
0.0
2.0
1.0
3.0
t II
ei-
Figure 5 Raw first order sensitivity over two periods (7 = 1.488)
I
/
I
o*“*
-1.0
0.0
0.5
w
1.0
1.5
t Figure 4
Phase shift coefficients
atlay
for steady state orbit
y(p) + @y/all) Ap) to predict trajectories for perturbed parameter values, absurd results can be obtained, even for very small Ag. The physical requirement that yl< 1 .O is violated (Figure 6).
coefficients
with respect to I.C,
To obtain the path-independent parts of the raw coefficients, the parametric phase sensitivity is determined using equation (16). The sensitivity of the period can be read directly from the plot of these values (Figure 7). The pathindependent sensitivities (Figure 8) can then be found via equation (9). A plot of y(p + A/J) similar to Figure 6, but using the adjusted coefficients, reveals that not only the shrinkage of the limit cycle, but also the transient trajectory, is accurately predicted (Figure 9). The only deviation occurs at the peak of the trajectory, where slopes and sensitivities change quite rapidly. The second order correction should improve the prediction in this area. To demonstrate the utility of sensitivity analysis in predicting amplitudes and periods, first and second order estimates were prepared using the amplitude and period sensitivities in Taylor series expansions for various paramt.er
Appl.
Math. Modelling,
1984, Vol. 8, October
335
Sensitivity analysis of oscillatory systems: M. A. Kramer et al.
8.0
’
7.0 ’
6.0..
5.0.
Y2 4.0 .-
3.0 .. 4
2.0..
‘.O/
1.0 .’
0.5
0.6
0.7
0.8
0.9
1.0
YI Figure 6 Trajectory prediction using raw sensitivities with Ap = 0.05 in a first order Taylor series. Actual trajectory indicated by dashed line
1
00
1.0
3.t
2.0
t Figure 8
-2.0
Path-independent
values moderately removed from the nominal values (Tables I and 2). These predictions compare favourably with those of linearized theory17 and singular perturbation theory.” Frequency predictions from sensitivity analysis are, in fact, quite as good as thosk provided by a regression equation,15 derived by fitting data from 500 separate kinetic simulations. Amplitudes presented in TabZe 2 are modified values as defined by:
*.
-3x)
1 -4.0
‘-
-5.0
..
a’=-
2 [ cl+
I I I I I I I I*
L 0.0
1.0
2.0
3.f
t
Appl.
Math.
Modelling,
1984,
Y2 min
Y2 max
v2
mad
-
C1 + WY2mid
1
Synthesis of a beneficial periodic operating strategy for catalytic kinetics
The following reactions are assumed to occur in a wellmixed isothermal catalytic reactor: (i) A+S+AS
Figure 7 Sensitivity of phase to parameter period growth is indicated by dashed line
336
sensitivities
p. Linear period-to-
Vol. 8, October
(ii) B + S + BS
Sensitivity analysis of oscillatory systems: M. A. Kramer et al.
(iii) AS + AS 3 CS + S k, (iv) BS + CS -+ P + 2s
8.0
where S are catalytic surface sites and P is the desired product. The adsorption and desorption of A, B, and P are assumed to be very fast compared with the rate of reaction (iii). Given the reaction rate constants, adsorption equilibrium constants K*and KB, and the total molar flow rate per unit reactor volume, F, it is desired to maximize to yield of P, by varying only the feed composition. Using KA= KB = k4 = F = 1, and k3 = 0.1, the optimum fraction of A in the feed is Af= 0.812, giving P = 0.01158, in the steady state. This state may be calculated directly by coupling equations (1) and (2) with CX~ = Af, setting the entire LHS and aP/aA, to zero, and solving the resultant set of eight algebraic equations in the eight unknowns, A, B, CS, P, aA/i3Af, aB/aA,, XS/all,, and Af Sensitivity analysis was then undertaken to determine if the optimum steady state process is indeed the optimum process. Previous work approached this question by calculating the second variation about the optimal steady state by asymptotic methods in the limits of high or low frequency forcing,” or by the so-called ‘pi-criterion’ test of Bittanti et a1.19 Alternatively, perturbation approaches have been used (e.g. Skeirek and Devera),” and these methods most closely resemble the sensitivity treatment, particularly in their ability to handle any cycling frequency. Horn and LinzO showed that the first variations in the periodic variables around an optimum steady state process are zero, and therefore as a first step in synthesizing a periodic control, the second order sensitivities a2P/a(u2 must be evaluated. In the current process, the controllable parameters are the fractions of A and B in the feed stream, and the cycling frequency w = 21r/r. Periodic flows of the type Af(t) = Af( 1 + fA sin wt) and Bf(t) = Bf( 1 + en sin(ot + 6)) were
7.0
6.0
5.0
Y2 4.0
3.0
2.0
1.0
0.0
0.5
0.6
0.7
0.8
0.9
1.0
YI figure 9 Actual and first-order predicted trajectories with Ap = 0.05 using path-independent sensitivities. Actual trajectory indicated by dashed line
Table 1
Exact and estimated
periods for a CSTR
in cases considered by Hugo and Wirgesls
Parameters
0.8 0.8 0.6 0.8 0.75 0.6 0.55 aSensitivity
Table 2
Period, r
E
B
IJ
Exact
Linear theoryI
Singular perturbation”
First ordera sensitivity
Second ordera sensitivity
0 0 0 0.02075 0.02075 0.02075 0.02075
10 10 10 30 30 30 30
0.280 0.2509 0.261 0.3118 0.310 0.3038 0.2716
2.611 3.278 6.698 1.267 1.453 1.989 2.625
2.566 2.967 11.20 1.823 3.050 Imaginary Imaginary
_ 0.48 0.56 0.76 1.02
2.818 3.308 4.029 0.7045 0.9576 1.731 2.496
3.225 4.032 5.064 1 .a21 2.004 2.285 2.618
results based on Taylor
Exact and estimated
series with nominal parameter
modified
amplitudes
for CSTR
values us = 0.7, E = 0.010375,
in cases considered Modified
Parameters
amplitude,
b = 20, p = 0.284
by Hugo and Wirges” a’
us
E
B
fi
Exact
Empirical approximation12
Singular perturbation”
First order sensitivitya
Second order sensitivitya
0.8 0.6 0.8 0.7 0.65
0 0 0.02075 0.02075 0.02075
10 IO 30 30 30
0.297 0.261 0.573 0.467 0.445
0.275 2.21 1.19 3.82 4.45
1.22 1.79 1.44 3.84 4.33
2.00 3.20 2.83 4.10 4.61
2.16 3.64 2.77 5.71 4.62
0.89 2.23 1.37 4.39 4.02
‘Sensitivity
results based on Taylor
series with nominal parameter
values us = 0.7, E = 0.010375,
Appl.
Math.
b = 20, p = 0.417
Modelling,
1984,
Vol.
8, October
337
Sensitivity
analysis of oscillatory
systems: M. A. Kramer et al.
examined, where d;f+ Bf = 1. Taking two derivatives of P, the average of P over one period r. one obtains:
a2P -= ad
Qa7ap 7 Gas aff2 +,
1 i
(t + 7) + 1 (P(t + 7) - P) ; r
ar 2dP G --(If7)--1
2
ap a7
r
aa aa
(32)
Evaluated at eA = fu = 0, only the first two terms of this expression are non-zero, and only the first is non-zero if (Yf w. The full expression, obtained by Liebnitz’s rule for differentiation of integrals, is actually the time average of the second order path-independent sensitivities, which reinforces the concept that these quantities are the true indicator of the parametric effects on the shape of the equilibrium orbit. However, with the aforementioned simplifications and because ar/acr is either -r/w or zero, depending on whether or not (Y= o, direct use of the elementary sensitivities is straightforward. The values of:
calculated by the AIM method4” are listed in Table 3 at various w. For this problem, both a’P/&iand a’P/aei were vanishingly small, although this is not the case generally. The AIM program, incidentally, provided an extremely efficient method of determining these sensitivities, for two reasons. First, since the sensitivities were evaluated around a steady state, the Jacobian was constant, and thus it only had to be diagonalized once. Diagonalization is usually a major part of the computational effort in AIM. Secondly, various values of w and 6 were handled in a single simulation, by considering each numerical value to be a distinct parameter. Using AIM, the marginal cost of the sensitivity to a new parameter is very small, and in fact the computational cost of all the data in Table 3 was only about twice that of a single kinetic simulation. Positive values of a2P/ &A a+, were found for 6 = 7~,Af = 0.812, in a range of frequencies from about 0.2 to 50, with a maximum of +0.0007 at w = 0.50. Extrapolation using a second order Taylor series indicates an improvement of about 6% over the steady optimum value with this frequency, and eA = es= 1. A simulation was performed to determine whether a periodic state existed at these parameters values,‘l and the true extent of yield augmentation. A limit cycle of period r = 2n/c; was indeed present, and a 5.7% increase in P was
Tab/e 3 Second order sensitivitie? state of example
around the optimal
steady
w
a*iaeA
1o-4 10-S 10-Z 10-l 1 10 IO2 103 lo4
-0.007 -0.007 -0.007 -0.002 + 0.005 +0.001 -0.00004 -0 -0
aEvaluated
338
Appl.
at e,q = EB = 0, S = rr, &
Math.
Modelling,
= 0.812
1984,
Vol.
8, October
acg
calculated. The question now arises as to whether this periodic mode is the most beneficial. The easiest way to answer this question is to evaluate the normalized sensitivity of P to the various parameters, on this limit cycle. This calculation gave a In P/a hr eA = t0.8, a In P/ ~ln~~~+0.05,~~ln~/~ll~~~<~O~4,~~~n~/~ln~~< 10d4, and I a In P/a In 6;,1 < 1Oe4. The first two sensitivities suggest further increases in eA and en, however eA = eu = 1 are the maximum permissible values. The final three sensitivities are approximately zero, indicating local maxima. Therefore, this is the optimum strategy for sinusoidal feed variations. A physical interpretation of this result is that the maximum segregation of A and B provided by 6 = 7~promotes a high concentration of AS, which stimulates the slow reaction (iii). during the first half of the period, whereas during the second half period, B is introduced in quantity to react with CS to form the product. This amounts to serially running reactions (iii) and (iv) while providing maximum surface coverage of the respective reactants by removing the competition between A and B for surface sites. This suggests that square wave feed profiles, providing further segregation of A and B, may be even more advantageous than sinusoids, and indeed, an increase to P = 0.01289 was obtained using bang-bang (on-off) control of A and B, using the same Af, Bf, and w as previously. Further simulations indicated that this was the optimum control of its type. The advantages of the sensitivity approach in synthesizing periodic operating strategies are two-fold. First, a high degree of automation is introduced, with the potential to scan through a large number of possible control schemes at relatively low cost. without restriction on the complexity of the system. Second, and most importantly, one is not restricted to the immediate neighbourhood of the optimum steady state in making sensitivity estimates. The first estimate of the value of a periodic scheme, provided by the second derivatives around the steady state, may be followed by further simulations and sensitivity calculations on the limit cycle generated by the predicted best combination of parameters. Thus, an iterative scheme may be constructed to truly optimize the periodic control. In addition, it may be possible to use the methods presented earlier in this work to further the understanding of the physicochemical processes controlling the forced oscillation, and use this information to construct waveforms and strategies better suited to the process in question.
Conclusion The techniques of sensitivity analysis have been extended to treat systems which exhibit limit cycle behaviour. A method was developed to separate the raw sensitivity coefficients into path-dependent and path-independent portions. The path-dependent data was quantified in terms of the parametric dependencies of the phase of the system, and could be used to determine the sensitivity of the cycle period, and the sensitivity of the time of occurrence of features such as extrema. The path-independent sensitivities were shown to contain information on perturbative dynamics and the shape of the equilibrium orbit. Each of these analyses were developed to the second order. A simple example involving the limit cycle behaviour of a well-mixed chemical reactor was presented in which these methods were used to predict parametric sensitivities of limit cycle phase, dynamics, shape, amplitude, and period.
Sensitivity
Additionally, it was shown that the sensitivity approach can be used to design beneficial periodic operating strategies. Another method suggested as a possible approach to these problems, known as feature sensitivity analysis,22 deserves mention at this point. In this procedure, one finds a functional ‘fit’ for each oscillatory trajectory, in terms of a set of new parameters p, viz., y = g(p), where each of these parameters has a specific significance in terms of the shape or dimension of a particular ‘feature’ of the trajectory y. The partial derivaties @/acr quantitatively link the parameters of the original system to the features of the output y. These quantities are determined by solving a set of linear algebraic equations involving the sensitivities ay/ap and the conventional elementary sensitivities. The utility of this approach seems to be limited by the ability to find an appropriate function g which not only accurately fits the oscillations, but also has a clear physicality in terms of its parametric dependencies. Unless there is an a priori reason for selecting a particular functional form for g, finding a suitable function may require significant effort, and the current approach is probably preferable. The sensitivity analysis approach provides modellers with a method of probing the detailed relationship between the various aspects of limit cycle behaviour and the system parameters, in a quantitative fashion. Such information can be of significant value in both theoretical and experimental studies involved in limit cycle characterization, or in determination of the phenomena responsible for oscillatory behaviour. Furthermore, sensitivity analysis provides a natural method of developing correlations for limit cycle period and amplitude, with an economy of effort. Significantly, the method is not limited by the degree of nonlinearity in the model, and does not depend on any assumption of near-harmonic or small frequency oscillations. The sensitivity approach shows great promise as synthesis method for periodic operation strategies, and work continues on this problem.
analysis
of oscillatory
systems:
M. A. Kramer
et al.
Fife, P. C.J. Phys. Chem. 1981,85, 2861 Larter, R., Rabitz, H. and Kramer, M. A, ‘Sensitivity analysis of limit cycles with applications to the Brusselator’,J. Chem. Phys. 1984,80,4120 13b Larter, R. ‘Sensitivity analysis of autonomous oscillators: separation of secular terms and determination of structural stabiJity’,J. Phys. Chem. 1983, 87, 3114 14 Edelson. D. and Thomas. V. M. J. Phvs. Chem. 1981.85. 1555 ’ 15 Hugo, P. and Wirges, H. P. 5th Int. Symp. Chem. Reaction Eng., Houston, 1978, p. 498 16 Uppal, M., Ray, W. H. and Poore, M. B., Chem. Eng. Sci. 12 13a
1914,29,961 17
18 19 20 21
22
Wicke, E. ‘Grundlgen der Chemische Prozessregelung’ (Eds Oppelt and Wicke), Oldenbourge Verlag, Miinchen, 1964, p. 54 Bailey, J. E. ‘Periodic optimization’ (Ed. A. Marzollo), Springer-Verlag/CISM, Wien, 1972 Bittanti, S., Gronza, C. and Guardabassi, G. IEEE Trans.
Auto. Control 1973. AC-18. 33 Horn, F. J. M. and iin, R. d. Ind. Eng. Chem. Des. Dev. 1967, 6, 21 Periodic states are not guaranteed, and theoretical proof of existence is problematic. See Sincic, D. and Bailey, J. E., Chem. Eng. Sci. 1977, 32,281 Skumanich, M. and Rabitz, H., Commun. J. Molecular Sci.
1982,2,79
Appendix Derivation
of selected
second order expressions
The detailed derivations of several of the expressions used in the paper are presented below. Equation
(23). We begin with equation
and differentiate
(9), written as:
with respect to &k at constant
r$:
Acknowledgement
(AlI
The authors are grateful to the Office of Naval Research for support of this work.
References Tomovic, R. and Vukobratovic, M. ‘General sensitivity theory’, American Elsevier, New York, 1972 Dickinson, R. P. and Gelinas, R. J. J Comp. Phys. 1976, 21, 123 Hwang, J. T., Dougherty, E. P., Rabitz, S. and Rabitz, H. J.
Chem. Phys. 1978,69,5 Kramer,
180
M. A., Calo, J. M. and Rabitz,
8 9
Kramer, M. A., Calo, J. M., Rabitz, H. and Kee, R. J. ‘AIM: the analytically integrated Magnus method for linear and second order sensitivity coefficients’, Sandia National Laboratories report 1982, SAND82-8231 Demiralp, M. and Rabitz, H.J. Chem. Phys. 1981,75, 1810 Gilles, E. D. ‘Gundlagen der Chemische Prozessrogelung’ (Eds Oppelt and Wicke), Oldenbourg Verlag, Miinchen, 1964, p. 89 Beck, J. Chem. Eng. Set. 1977, 32, 159 Douglas, J. M. and Caitondo, N. Y. Jnd. Chem. Fund. 1967,6,
265 10 11
--a
aYi
i affk a+
) =“‘+(_EJLJ%) G
aai ask
I
H. Appl. Math. Model-
ling 1981, 5, 432
6 7
To simplify this expression, it must be realized that yi is a dummy variable in equation (9) and that equation (9) is the general definition of partial differentiation at constant $. Applying equation (9) to the first term on the RHS of equation (Al), we obtain:
Chang, H. C. and Calo, J. M. Chem. Eng. Cammun. 1979,3, 431 Skeirek. R. D. and Devera, A. L. Chem. Eng. Sci. 1982, 37, 1015
The second term may be similarly evaluated:
Substituting
into equation
(Al), we obtain equation (23).
Equation (25). We begin with equation (23), evaluated at t + T. Several terms arise which will need to be simplified, namely :
Appl.
Math.
Modelling,
1984,
Vol. 8, October
339
Sensitivity analysis of oscillatory systems: M. A. Kramer et al.
which, after simplification,
and
i($f+?)) =-g$&Z)
g
(t + 7) i
gives:
&
k
(t) k
I
=$(2 (t) j-g dz(t) In the following, all terms are evaluated at time t, unless otherwise indicated. From equation (23), we have:
ar
&r d2yi
dyi
+aol_cdt2---I k
a2yi (t + 7) aai ask
a27
dt affj aok
Solving for a2r/aoj aok we obtain equation
(25).
Equation (29). Begin with an expression analogous to equation (9), with ‘e’ appearing in place of 4, indicating decomposition in the direction of change of the extrema:
avi _ ayi~at,dYi
0aai e
Differentiate
Subtracting a*Yi aaj
(t
equation +
r)
&Yj at
(-42)
with respect to ok:
(23) evaluated at time t, we obtain:
a33 aai ask
ask
aai
+dyi
(‘)
a2te
dt sol, aok The last term in this expression is zero because dyi(tJ dt = 0. The second-to-last term is also zero, because along e (the locus of all phase-space points through which the extreme passes as a function of Q) dYi/dt is by definition zero. The sole remaining term can be simplified using equation (A2): (A3)
dyi
a27
dt
aaiask
Appl.
Math.
--p
340
Modelling,
Applying equation (28), it can easily be seen that equation (A3) is the same as equation (29).
1984,
Vol. 8, October
analysis sys terns
of
Mark A. Kramer* Department of Chemical 08544, USA
Engineering,
Princeton
University,
Princeton,
New Jersey
Princeton
University,
Princeton,
New Jersey
Herschel Rabitz Department
of Chemistry,
08544,
USA
Joseph M. Calo Division of Engineering, Brown University, Providence, (Received March 1983; revised August 1983}
Rhode
Island 02912,
USA
Methods of sensitivity analysis are extended to find the parametric dependencies of systems of ordinary differential equations which exhibit limit cycle oscillations. The quantitative relations between the system parameters and the observable period, amplitude, phase and cycle shape are developed. These formulae, presented for both the first and second order, are applicable to systems of arbitrary size and complexity. The techniques are used here to develop correlations for period and amplitude in a non-isothermal oscillating stirred-tank chemical reactor, and to find and optimize a beneficial periodic operating strategy for a lumpedparameter catalytic reacting system. Key words: mathematical oscillatory systems
Sensitivity analysis, the systematic study of the effects of parameter values on the predictions of mathematical models, is a useful tool for model validation and evaluation, as well as for quantifying the effects of parametric uncertainty and variability. The numerical link between parameters and predictions in systems of ordinary differential equations (ODES) is provided by the linear sensitivity coefficients’ aui(t)/aoll, where yi is the ith dependent variable, CX~ is the ith system parameter, and t (time) is the variable of integration. These coefficients may be obtained by several numerical techniques ‘-’ if the model is not analytically solvable. Higher order sensitivities, such as a2yi(t)/$ a&, are also calculable, and provide information on nonlinearity and parametric interactions. The interpretation and application of sensitivity information is a rapidly-developing area.6 However, the available analysis is not directly applicable to oscillatory systems. *Present address: chusetts Institute USA.
328
Department of Chemical Engineering, Massaof Technology, Cambridge, Massachusetts 02139,
Appl. Math. Modelling,
1984,
Vol.
8, October
model,
sensitivity
analysis,
limit
cycles,
Furthermore, quantities of obvious interest for oscillatory trajectories, such as frequency and amplitude, are not directly addressed by the sensitivity coefficients. These complications will be shown to arise from the parametric dependencies of the velocity of the trajectory around the limit cycle. A method of analysis will be presented which isolates the parametric sensitivity of the cycle phase, thereby rendering explicit the parametric effects on period, amplitude, shape, size, and dynamics of the oscillation, and the associated transient. The procedures developed in this work are applicable only to trajectories which evolve to limit cycles, that is, to systems in which the equilibrium orbits are independent of the initial conditions. This type of behaviour arises in many problems in biology, catalytic kinetics, and reactor dynamics. For non-limit cycle oscillations, several of the formulae presented here are valid - in particular, equations (lo)-(1 5). However, it is impossible to define the sensitivity of the phase of the oscillation, as is done below, which is crucial to subsequent developments. 0307-904X/84/05328-13/$03.00 Q 1984 Butterworth & Co. (Publishers)
Ltd
Sensitivity
Parametric studies of limit cycles have usually focused on determination of parametric effects ‘in-the-large’, that is, determination of the ranges of parameter values for which a particular system may have oscillatory solutions. In contrast, relatively little effort has been directed at ‘finegrained’ analysis, in which the effect of parameter perturbations on specific equilibrium orbits is investigated. This information can be useful in fitting a model to data, in uncertainty analysis of modelling predictions, in experimental design, in identification of governing or underlying processes, or in optimization of periodically controlled processes. These applications, excepting the last one, are discussed in detail in Demiralp and Rabitz.’ An example of periodic optimization is included in this paper. The techniques presented here are applicable to any system described by ordinary differential equations. Previous efforts to describe limit cycle behaviour relied on assumptions of near-linearity and near-harmonic oscillations,7-g or used singular perturbation techniques which become inaccurate for large values of the perturbation variable.1o~‘1 Numerical comparisons between these methods and sensitivity analysis will be presented.
Sensitivity equations The system under consideration is a set of coupled ODES, which may be written in the following general form:
(1) Here, y is the dependent variable vector (y E R”), a is the parameter vector (ar E Rm), and f is the variable of integration. The parameter vector may include the initial conditions y(O), although y(0) usually does not enter explicitly into f. The latter vector is generally a nonlinear function of y and the non-initial value parameters. This equation can be differentiated with respect to the jth parameter CY~ to yield’:
d a~ dt
aai
= J(t)
2 (t) + E(r) J
(2)
J
where J is the n x IZ Jacobian matrix (J&t) = afi(t)&& and ay(t)/aol, is an n-vector of linear sensitivities. The initial condition for ay/aol, in equation (2) is the zero vector, unless aj is the initial condition for yi (i.e. CX~ f yi(O), in which case ayi(O)/aai = 1). Equations (1) and (2) can be integrated simultaneously, with an appropriate differential equation solver, to yield the linear sensitivities. Alternatively, an n x n Green’s function matrix, K(t, t’), may be constructed by solving3: ; K(t, t’) - J(t) K(t, t’) = 0
K(t’, t’) = I, I 2 t’
analysis of oscillatory
The relevant equation
systems: M. A. Kramer et al.
for the adjoint is simply:
K+(t’, t) + K+(t’, t) J(i) = 0
5
K+(t, t) = I, t’ G t
(5)
Equations (3)-(S) are exact, and yield the same ‘raw’ sensitivity values as equation (2). These ‘raw’ coefficients will be manipulated to yield information pertinent in the limit cycle case.
Sensitivity of phase to initial state The system (1) will be assumed to have a limit cycle solution, i.e. the steady-state orbit is independent of the initial conditions. We will first consider sensitivities of the type [avi(t)/ayi(0)]., which represent the change inYj(t) induced by a differential charge in yi(O), all other yk(O) and all parameters held fixed. The system state at t = 0 may lie anywhere within the range of attraction of the limit cycle. Since the steady-state orbit is unaffected by such a perturbation, only the position along the orbit (phase) can be influenced by yi(O), after an initial transient period. This phase offset can be quantified as the coefficient at/ayi(0), where at is the ultimate differential delay or advance along the limit cycle resulting from the initial differential perturbation ayi(0), measured in time units. This is shown schematically in Figure 1. The formula for calculating this sensitivity is easily derived by taking the total differential ofyi(or, t’, y(O)), then dividing by dyi(O) with dvi = 0 and at constant cr and yk(O), k #j. The imposition that dyi = 0 follows from the assumption that the limit cycle is ultimately achieved. Solving for the phase sensitivity, we obtain: Qj(0)
at
= -=-r!h+n[$)/f~)
(6)
auj (0)
The limit appearing in this formula is necessitated by the fact that the trajectory never actually reaches the limit cycle, but only approaches it asymptotically. Numerically, the number of revolutions required to converge to an accurate value of Q is governed by the stiffness of the system.” Equation (6) iS independent of the choice of i, since the
STEADY - STATE
#ORBIT
(3) DIFFERENTIAL
The linear sensitivities may then be expressed in integral, rather than differential form: r $
(t) = K(t, 0) 2 J
(0) + j- K(t, t’) a: J
0
(t’) dt’
(4)
J
Since t’ G t, it is actually more practical to solve for the adjoint Green’s function, K+(t’, t), subsequently applying the equality K+(t’, t) = K(t, t’) for use in equation (4).7
Figure 1 Schematic representation in two dimensions of the sensitivity of the ultimate cycle phase to initial conditions
Appl.
Math.
Modelling,
1984,
Vol. 8, October
329
Sensitivity
analysis of oscillatory
systems:
M. A. Kramer
et al.
cycle phase is a property of the entire system, rather than single variables. One further property to be observed here is that @i(t)/ayi(0) must pass through zero whenever dyj(t)/dt is zero, since the phase sensitivity must be finite. Special numerical attention (perhaps application of L’Hospital’s rule) may be needed at these points. Sensitivity of the ultimate cycle phase to perturbations of the dependent variables (y) is of particular importance when the initial system state is numerically ‘on’ the steadystate orbit. These on-cycle phase sensitivities are closely related to the Lyaponov stability of the cycle, and are, in fact, linearly related to the stability eigenvectors.13 These values will be used subsequently to calculate the sensitivity of the system phase to non-initial condition parameters. For this application, initial condition phase sensitivities are required at each point the parametric phase sensitivities are to be calculated. If we were to use equation (3) to obtain the initial condition phase sensitivities, as suggested by equation (6) (since Kii(f, 0) = avi(t)/vj(O)), a new solution of the matrix differential equation (3) would be necessary, using each point of interest in turn as the zero time in equation (3), each time the parametric phase sensitivities were to be calculated. Instead, it is preferable to use equation (5) - the adjoint Green’s function - which need only be solved once to obtain the entire set of on-cycle phase sensitivities. The solution of (5) proceeds backwards along the trajectory beginning at an arbitrary point t’ on the steady-state cycle, and produces Kt(t, t’), where t’> t. The phase sensitivity for a perturbation at t is then calculated from:
Qdt>= ,!e_ - K&(t,t’)
I
dvi@‘I dt
(7)
Note that only one row of Kt is required, so equation (5) can be solved as a vector rather than a matrix, thereby saving considerable computational effort. As a matter of practicality, the limit appearing in equation (7) needs only to be observed to the extent that a problem-dependent transient period is passed. This transient has been found, in several stiff systems, to be somewhat shorter than one period. Specific times may be checked to ensure they lie far enough from t’ for equation (7) to be applicable. The following relationship must hold for the calculated on-cycle phase sensitivities: f
Qk(t)
y
= -
(8)
1
k=l
The relationship is derived in the following way. Take the total differential of yi(t’) = yi(t’, y(t)), and divide by dt with t’ fixed (0~is also fixed for the purposes of this derivation). Thus:
The phase sensitivities can be used to identify the directions in which the system may be perturbed without disturbing the phase. A perturbation dy applied at time t leaves the phase of the steady state orbit unchanged if and only if: n
c
&c(t) dYk= 0
k=l
We will define @o(t)to be the collection of all vectors y’(t) evolving to the same orbital phase as y(t) at long times. Any differential dy satisfying the previous summation lies in 4(t). A perturbation not lying in g(t) changes the system phase by the amount dt, where y(t) + dy lies in (((t + dt). Mapping dy back on to G(t) removes phase shift effects and leaves a differential dy’ which satisfies the previous summation. This mapping process will be used subsequently to isolate the pathindependent parts of parametric sensitivities. The dynamics of response to differential perturbations in the dependent variables are contained in the transient period of K. To extract these dynamics, the phase shift component Qj(O)(dy(t)/dt), is subtracted from the ith column of K(t, 0). This leaves the component of the perturbation along the constant-phase locus g(t). The norm of this remainder is the distance to the nominal trajectory along G(t). The damping of the original perturbation may be examined by preparing a plot of this distance against time. The function must asymptotically approach zero for stable limit cycles. A word should be said at this point about sustained oscillations which are not limit cycles. The initial condition sensitivities for these oscillations do not demonstrate the simple behaviour of equation (6). Instead, the sensitivities contain information both on the phase sensitivities and on the sensitivity of the shape or location of the steady-state cycle. Unfortunately, these two parts cannot be separated unambiguously (although equation (14) can be used to determine the period sensitivity), since the current definition of initial condition phase sensitivities requires that the nominal and differentially-perturbed trajectories converge to the same path, so one trajectory is strictly leading or following the other. If the two trajectories follow different paths, there is no obvious way of judging which is ‘ahead’ and which ‘behind’. Work is continuing on this problem, since it would be useful to generalize the analysis presented here to non-limit cycle oscillations, and in fact, to altogether non-oscillatory trajectories. In the latter case, it is hoped that conventional sensitivities could be divided into two parts, representing the sensitivity of the velocity of the trajectory and the sensitivity of the shape of the trajectory, which would be analogous to the current treatment of the limit cycle sensitivities.
Sensitivity d@) -= dt
’ aYi (t’) dYk(t) c--
k=
1 ayk(t>
dt
Dividing both sides by dyj(t’)/dt, and applying equation (7), we obtain equation (8). If equation (8) is violated, either t is too close to t’, or the numerical accuracy of the adjoint kernel, and hence Qk, is poor. To get Qk(t) for t’- t less than the transient time, equation (5) must be integrated around the entire cycle (backwards) to time t - f. Then &(t) = &(t - r), since t and t - T refer to the same position on the cycle.
330
Appl.
Math.
Modelling,
1984,
Vol. 8, October
to non-initial
condition
parameters
Unlike initial conditions within the attractor region, parameters appearing explicitly in equation (1) can affect the period, amplitude, shape, and location of the limit cycle. This leads to additional difficulties in the interpretation of the raw sensitivity coefficients, since the sensitivity of these features, as well as the phase-shift information, are together contained in the raw coefficients. In order to deconvolute this information, one must separate the parts of the raw sensitivities which result from a parametric phase dependence, and those which refer only to the shape and position of the trajectory in phase space. The former will be
Sensitivity
referred to as the path-dependent part of the raw sensitivity, and the latter the position-dependent (or pathindependent) parts. As in the case of initial condition sensitivities, the contribution to the raw sensitivity due to the phase lead or lag can be written as the product of a differential phase shift, (at/sol,),, and the local slope,13a so the total sensitivity may be written:
analysis
of oscillatory
systems:
M. A. Kramer
et al.
In equation (13b), (dyi/daj), is the derivative of equation (10) taken with r (but not y(O))-fixed. Equation (11) may therefore be written:
t
a7 dJJi
7
aai dt
(14aJ
which, because (dyk(O)/dQj), = dYk(O)/dOlj, is equivalent to:
(9) where (ayi/&xj)@ refers to the path-independent portion of the tbtal sensitivity. The formula for calculating (at,/a~~)G will be developed shortly, and the meaning of the subscript $ will be discussed, but first it will be demonstrated that division of the raw sensitivities, as indicated in equation (9), is a natural consequence of partial differentiation of y with respect to the system parameters. If the system (1) has a periodic solution, yi(t) can be represented as a Fourier series: m
yi(t)
=
2nd Qi, cos -
C
2nnt + b’,
7
n=O
(10)
sin 7
where 7 is the period of the oscillation. The coefficients af, bi, and 7 are implicit functions of the system parameters. Differentiating equation (10) with respect o+, the linear sensitivity coefficient can be expressed as: dyi(t) dai
ar -
2nt
1 (mzI, sin ?Z aiinzo 7
=72 -
aat Csol-I
+n:o
2nW cos -
7
+
abi
2nN
-
-=
which, in general, is not zero. The sensitivity as previously defined is the partial derivative of yi with respect to CY~ at constant y(O), i.e. (ayi(t)/&q),(,), which is zero at t = 0. TO correct equation (1 l), consider Yi = yi [CY,t, yo(~l)] , which is the functional form implicit in equation (10). Taking the total differential, and dividing by daj, with dt = 0 and all other parameters fixed, we obtain:
dori
am
n
dYi(t)
CPP
k=l
aYk(o)
d&O)
(12)
d&i
The final quantity here cpmes from equation (1 l), evaluated at t = 0. Before applying this expression, note that: = sio
[-nab
sin (z)
2
Qk(t’)
7)
$
k=l
and that: [$cosfT)+$sin(y)] (13b)
(16) I
7 (13a)
(gi=,Z,
+ a0ij
Equation (15) is valid for all t, if equation (10) is satisfied. The summation term in equation (15) damps to zero for large t, and when this term can be ignored, only the sensitivities ayJ?Ior, are required to calculate the sensitivity of the period. Edelson14 used this approach, in analysis of the ‘Oregonator’, a prototype oscillator. The sensitivity of the phase to non-initial condition parameters will now be derived. To begin, the role of y and t as dependent and independent variables, respectively, must be reversed. Along any phase-space trajectory, a point may be described either by its coordinates y, or the time r at which the point is reached. There is, in fact, a one-to-one correspondence between y and t along a given path because phase-space trajectories never cross or coalesce, provided that the initial conditions lie off the limit cycle. Therefore, along a given path, time may be written as a function of position (i.e. t = t(y(a)). The parametric sensitivity of this function may be probed. Taking the total differential and then the partial differential with respect to aj:
+ nb’, cos (El] 7
avi(f
(15)
=
z
-
aCYi r 1 (11)
Note that the total differential symbol is used in this expression; indeed, this dyi(t)/daj is not the sensitivity coefficient one calculates using equations (2) or (4). This can be seen by evaluating equation (11) at t = 0, which gives: m 2-i dyj(0)/daj = c ““” n=O aai
@i(t) -=_-sol,
Equation (14b) was also given by Tomovic’ and !Larter,13 and is also known as a result of multiple scale analysis in perturbation theory. There is a general similarity in the appearance of equations (14b) and (9), however, (ay/acu), and (ay/acu)G are not equivalent, and (at/&x), # t/@-/aa). It will be shown subsequently that equation (14b) is deficient because only the overall phase shift behaviour is taken into account. Nonetheless, it can be seen from this derivation that application of the chain rule in partial differentiation of y quite naturally gives rise to the division of the raw sensitivity into two parts. Equation (14a), incidentally, suggests a method of isolating the period sensitivity. Since (dyi(t)/daj)T is periodic (see equation (13b)), the period sensitivity can be calculated from the sensitivities for adjoining periods”:
a7 __hi(t) sol, i\ aai
- nb; cos ““‘) 7
-sin
(14b)
The sensitivity (at(t’)/&+)+ is a measure of the phase lead or lag, in time units, at time t’, resulting from the differential kj. The subscript C$is included because this phase sehsitivity gives the phase shift in reference to the constarit-phase locus @(t’). It will be recalled that @(t’) is the set of all phase-space points, which, considered as initial conditions, ultimately give rise to the same phase on the
Appl.
Math.
Modelling,
1984,
Vol. 8, October
331
Sensitivity analysis of oscillatory systems: M. A. Kramer et al.
mentary sensitivity coefficients. If the initial conditions are numerically ‘on’ the nominal steady-state orbit, the following formula is valid:
+aCY y2
Figure 2 Phase shift due to perturbation in the non-initial condition parameter 0~. The phase difference indicated is at = t’ - r”
steady-state orbit as the nominal solution (evolving from y = y(t’) at t = t’). Considering a nominal and a parametrically perturbed trajectory, if the perturbed trajectory crosses tit’) at t = t”, then the phase shift given by equation (16) is (at/da), = t’- t” (see Figure 2). An alternative derivation of equation (16) is as follows. The parametric phase sensitivity at t = t’ can be found by (a), solving the model using the parameter value oj + doj, from time zero to t = t’, and (b), returning the parameter to its nominal value, and solving the model until the trajectory closely approaches the steady-state cycle at t = t”‘. The phase sensitivity is then given by comparison with the nominal solution of the model at this time, using equation (6). After step (a), the difference between the nominal and perturbed trajectories at t = t” is dy = (ay(t’)/aaj) daj. In step (b), dy is considered a non-parametric perturbation from the nominal traje’ctory, and the ultimate phase shift arising from dy is calculated from:
This is surprising in view of the fact that, even when initial coditions are on-cycle, there is a transient in ay/aaj, prior to establishment of regular behaviour, as given by equation (15). There is, however, a structure to the components of ay/acri during this transient, which makes at/iJol, regular from time zero (i.e. equation (17a) valid for t’ -+ O+). A geometric interpretation of equation (16) can be used to explain this effect. Equation (16) maps the perturbed trajectory back to the nominal trajectory along a line of constant phase $, and then measures the phase shift. Therefore, the same phase sensitivity will be calculated whether the perturbed trajectory lies between the nominal limit cycle and the perturbed limit cycle (corresponding to the transient in ay/aoj), or the perturbed trajectory lies on the perturbed steady-state cycle (corresponding to regular behaviour of ay/&rj, indicated by complete damping of the summation term in equation (15)). Thus, the transient in the linear sensitivities has no effect on the phase sensitivities. We may now return to the comparison of equations (9) and (14b). It should be clear that the period sensitivity is only one aspect of the total parametric phase sensitivity. The period sensitivity is the overall period-to-period phase shift, while the parametric phase sensitivity (at/aaj) encompasses not only this overall behaviour, but also the local variations in the phase lead or lag. These local variations, being functions of position rather than path, can be described by a periodic function (at/sol,),:
(18) if the trajectory is effectively on the limit cycle. Separation of the total time scale sensitivity into these two components allows calculation of at/sol, at any time once a single cycle has been analysed. The total sensitivity is then used in equation (9) to determine the path-independent portion of the raw parametrjc sensitivity. It can be shown that the path-independent sensitivities (ay(t)/aor,), do indeed lie on the constant-phase locus tit), justifying the subscript 4 in equation (16). Beginning with equations (9) and (16) and applying equation (8):
-f
ti?k(6
(y)@
= ’
(19)
k=l
dt = f
Q,&‘) d2’k
k=l
which follows from the definition of Q. Substituting for dyk in this expression, and dividing by doj, gives equation (16). A formula for the sensitivity of the period of the oscillation can now be derived. This sensitivity is simply the accumulated phase lead or lag after completion of one revolution around the steady state orbit, SO that:
(17a) where t’ is any time longer than the transient
332
Appl.
Math.
Modelling,
1984,
Vol.
in the ele-
8, October
These sensitivities therefore conform to the definition of g(t). The dynamics of response to step changes in the parameters, when the trajectory is on cycle, is contained in the transient period of (ay/acu)Q The approach to the perturbed equilibrium orbit is given by the difference (ay(t + 7)/h), - (ay(t)/aa),. This will be demonstrated in the example below. The sensitivity of amplitude to parameters is now easily determined. Since at an extreme in yi, dyi/dt = 0, the path-independent sensitivities are equal to the raw coefficients at these points. The sensitivity of the amplitude to first order is therefore simply the difference of the raw coefficients at the extrema.
Sensitivity
analysis of oscillatory
systems:
M. A. Kramer
et al.
Formulae for the second order Second order sensitivity information permits analysis of parametric interaction, and/or increased accuracy in predictions of period, amplitude, or cycle shape. These coefficients may be obtained by direct solution of the governing equations (equations (21) and (24)) or by using a packaged sensitivity solver such as AIM,4P5which uses a modification of the Green’s function method.3 The AIM algorithm is particularly useful for second order calculations, since the integral tranformation of the Green’s function method reduces the computational cost of the second order coefficients to the same order as first order coefficients. Direct solution of equations (21) and (24) will, in general, be much more expensive. Second order sensitivity coefficients are treated in a similar fashion to the first order coefficients. Differentiation of equation (6) with respect to a second initial condition, yk(0), yields the following expression for second or phase sensitivities:
+;
(t,
0)
=
J(t)
kjk(t,
+
L(t)
Kj(t,
= aJ,,/ay,,
’
foli~j(‘)
sjk(o) = 0
aJ&)/aqW
(fajak@))i = a2fi(t)/aaj(t) a&k(t)
0)
Kk(t,
0)
and
and Sj(t) = ay(t)/aaj Equation (23) may be used to obtain an explicit expression for the second order period sensitivity. As in the first order case, this is done by comparing the sensitivities obtained on two successive periods. Noting that:
-=
a27
a&j aark
a27
_= aaj ask One may also obtain second order phase shifts directly for all points around the orbit by solving a single differential equation. This is similar to the adjoint procedure (equations (5) and (7)) used to obtain the full set of on-cycle first order phase sensitivities. The relevant equation, from differentiation of equation (5) and application of equation (7), is:
qjZ(t)Jdt)
a2t(t’ + 7) -- a?(i) a&jask
affj aak
where t’ is any time referring to a numerically position, we obtain:
It will be observed that, due to symmetry:
‘y =,ilQdt)L,,(t) -
+ J,(t) Sk(t)
(24)
(Jaj(Ohp =
(21)
where L,,
= J(t) S;k@
+ J&k(f) sj(t> + L(t) sj(t) Sk(t)
ia%)QjW) 0)
f$k(t)
;
aYk@)
for any species i. The second order initial condition sensitivities (kjk(t, O))i = a2yi(t)/dyi(0) ayk(O) are defined by: *jk 7
The derivation of this expression is given in the Appendix. Again, the subscript @refers to the path-independent portion of the total sensitivity. As in the first order case, it can be shown that the (ay(t)/hx, acq)@lie on N(t), where (t(t) now includes the curvature due to second order effects. The r’aw second order barametric derivatives Sik(t) = ay(t)/acq sol, are defined by the equation3:
where
a2Yi(t) aYj@)
(23)
(22)
where qjk(t) = a2t/@j(t) r&(t), and q(0) = 0. The solution of this equation proceeds backwards around the limit cycle. The starting point may be anywhere numerically ‘on’ the steady-state orbit, and, as in the first order case, there is a transient period associated with equation (22). This transient corresponds to the time required for the damping of the perturbations ayj(t) and ayk(t). For non-initial condition parameters, a decomposition formula may be obtained from differentiation of equation (9) yielding:
‘dyi -r a2yi(t’) F( dt 1 aaj hk
a2yi(t’ ‘f T) a&j a&
+ a7 a7 d2Yi _ 37 d affj aak, dt2
---
‘on’-cycle
Mt’)
&vj dt i ask i
-
(25)
where i is arbitrary. The derivation of this expression is also given in the Appendix. To evaluate equation (23), an expression for (a2t/ &q a&& must be developed. Differentiating equation (163, we obtain:
(26) Unfortunately, an additional differential equation must be solved to obtain’aQ/aak, since from equation (7) it can be seen that the explicit expression for a&(t)/&, Will involve the term:
APPI.
Math.
Modelling,
1984,
Vol. 8, October
333
Sensitivity analysis of oscillatory systems: M. A. Kramer et al.
where i is arbitrary and t’ is a time sufficiently removed from t such that the transient in equation (7) has run its course. Because f appears as a reference time in this term, it cannot be calculated from any of the raw sensitivities. The differential equation directly controlling aQl/acu, can be derived from partial differentiation of equation (5) with respect to 0~~:
+
aQm(t) ~ J,nl(t)
(27)
aLwk
Solution of this equation, like that of equations (5) and (22), proceeds in the reverse direction around the limit cycle, and provides the full set of on-cycle derivatives. Equation (26) can then be used to calculate the second order phase sensitivity to non-initial parameters. As in the first order case, this phase sensitivity can be compared on successive periods, in a manner analogous to equation (18), to calculate the second order period sensitivities. If the initial conditions are numerically ‘on’ the limit cycle, the second order period sensitivities are simply -(at(T)/ acvi aaj)e To determine the second order amplitude sensitivities, the effect of the parameters on the times of occurrence of the extrema must first be calculated. Since dyi(ol, r,)/dt = 0, where t, is the time of occurrence of the extrema, we obtain:
(28) by differentiation with, respect to oj, and solution for at,/&.xj. We wish now to find the component of a’yi/ sol, aa, in the direction of change of the extrema. This is given by (see Appendix): _
a2Yi(fe))
ateat,d2Yi(te>
aaj ih!,
a@j acUk dt2
(29)
The amplitude sensitivity is then the difference between the second order sensitivities in the direction of change of the extrema: -=
a2ai
a2Yi(fmin)
a’Yi(LJ
a@j aa,
aaj ibk t--
at,,, ihi
%,,
h!j a&!k d2Yi(fmax)
aaj
af,,,i, - afmin -aaj aa,
dt2 d’Yi(trnin) dt2
(30)
The procedure, then, for a full second order analysis is as follows. First, the raw coefficients are obtained by solving equation (24) or an equivalent method. No more than two cycles need be traversed to provide enough information for all the analyses. The second order period sensitivities may be obtained immediately via equation (25). Solution of equation (27) is then required to provide the necessary data for calculation of the time scale sensitivity by equation (26). The full decomposition is then executed using equation (23). The path-independent sensitivities thus isolated are used for estimating the parametric effects on the shape of the equilibrium orbit. Second order
334
Appl.
Math.
Modelling,
1984,
Vol. 8, October
phase sensitivities may be obtained by solving equation (22). Second order amplitude sensitivities are calculated from equations (28) and (30).
Examples Prediction of amplitude, period, and dynamic behaviour of an oscillating CSTR
To demonstrate some of the applications of the preceding techniques, a simple perfectly-mixed, lumped parameter chemical reactor model will be used. Although this model has been thoroughly studied in terms of relations between parameters and stability, the relations between specific limit cycle trajectories and parameters have not been systematically analysed. In this regard, sensitivity analysis offers some new insights. The specific system under consideration is a well-mixed continuous non-adiabatic stirred tank reactor (CSTR) in which a single first order, irreversible exothermic reaction occurs. This system is represented by the following two coupled differential equations,15 representative of the dimensionless transient mass and energy balances for the reactor: j1 = -Y1
j,=b
+ Da(l -yJ
-py2
Y2
exp 1+ fY2
+ Da(1 -ye)
(31) Y2
exp 1+ EY2
where Da = us exp [---u&p + e~~)]/( 1 - ~3. In this equation, b is the dimensionless enthalpy of reaction, E is the dimensionless activation energy, U, is the dimensionless concentration at the system steady state for a stable system, or at the limit cycle focus for unstable systems, and p is the dimensionless heat transfer parameter. The nominal parameter values b = 20, E = 0.010375, us = 0.7 and p = 0.417 give rise to an unstable condition, indicated by det J = 14.4, tr J = 1.8, at the point at whichjl =i2 = 0. (The conditions det J > 0 and tr J < 0 are necessary and sufficient for the point to be locally asymptotically stable in this particular system).16 The trajectory for the nominal parameter values is shown in Figure 3. An arbitrary point on the cycle y = (0.8851, 1.678) was chosen as the reference point for the sensitivity analysis (i.e. t = 0 at this point). Beginning at the reference point, the adjoint kernel equation (5) was solved by coupling it with equation (3 1). The resulting values of K 7 were divided by dy(O)/dt according to equation (7) to produce the initial condition phase sensitivities (see Figure 4). The definition of the constant-phase locus 4(t) (zQk dyk = 0) was then used to determine the approximate course of these loci in the immediate vicinity of the orbit. The further requirement that all @((t)meet at the unstable focus of the limit cycle, which results from the fact that a randomly directed perturbation away from this point can lead to any on-cycle phase, allows the position of the constant-phase loci to be estimated for the entire cycle interior. The results are shown in Figure 3. The raw coefficients were then obtained using equation (2). The sensitivities with respect to the parameter p are shown in Figure 5. The strong peaks in these values result from the ‘saw-tooth’ (relaxation oscillation) character of the individual trajectories. The non-periodicity predicted by equation (9) is easily observed. If one naively uses these values in a first order Taylor series expansion (y(/* + AP) =
Sensitivity
1
6.0
analysis of oscillatory
systems:
M. A. Kramer
et al.
60C
7.0
6.0
YI figure 3
Steady state orbit for CSTR. @(t) are shown every l/10 period
Lines of constant phase
I
I 3.0 I
I
0.0
2.0
1.0
3.0
t II
ei-
Figure 5 Raw first order sensitivity over two periods (7 = 1.488)
I
/
I
o*“*
-1.0
0.0
0.5
w
1.0
1.5
t Figure 4
Phase shift coefficients
atlay
for steady state orbit
y(p) + @y/all) Ap) to predict trajectories for perturbed parameter values, absurd results can be obtained, even for very small Ag. The physical requirement that yl< 1 .O is violated (Figure 6).
coefficients
with respect to I.C,
To obtain the path-independent parts of the raw coefficients, the parametric phase sensitivity is determined using equation (16). The sensitivity of the period can be read directly from the plot of these values (Figure 7). The pathindependent sensitivities (Figure 8) can then be found via equation (9). A plot of y(p + A/J) similar to Figure 6, but using the adjusted coefficients, reveals that not only the shrinkage of the limit cycle, but also the transient trajectory, is accurately predicted (Figure 9). The only deviation occurs at the peak of the trajectory, where slopes and sensitivities change quite rapidly. The second order correction should improve the prediction in this area. To demonstrate the utility of sensitivity analysis in predicting amplitudes and periods, first and second order estimates were prepared using the amplitude and period sensitivities in Taylor series expansions for various paramt.er
Appl.
Math. Modelling,
1984, Vol. 8, October
335
Sensitivity analysis of oscillatory systems: M. A. Kramer et al.
8.0
’
7.0 ’
6.0..
5.0.
Y2 4.0 .-
3.0 .. 4
2.0..
‘.O/
1.0 .’
0.5
0.6
0.7
0.8
0.9
1.0
YI Figure 6 Trajectory prediction using raw sensitivities with Ap = 0.05 in a first order Taylor series. Actual trajectory indicated by dashed line
1
00
1.0
3.t
2.0
t Figure 8
-2.0
Path-independent
values moderately removed from the nominal values (Tables I and 2). These predictions compare favourably with those of linearized theory17 and singular perturbation theory.” Frequency predictions from sensitivity analysis are, in fact, quite as good as thosk provided by a regression equation,15 derived by fitting data from 500 separate kinetic simulations. Amplitudes presented in TabZe 2 are modified values as defined by:
*.
-3x)
1 -4.0
‘-
-5.0
..
a’=-
2 [ cl+
I I I I I I I I*
L 0.0
1.0
2.0
3.f
t
Appl.
Math.
Modelling,
1984,
Y2 min
Y2 max
v2
mad
-
C1 + WY2mid
1
Synthesis of a beneficial periodic operating strategy for catalytic kinetics
The following reactions are assumed to occur in a wellmixed isothermal catalytic reactor: (i) A+S+AS
Figure 7 Sensitivity of phase to parameter period growth is indicated by dashed line
336
sensitivities
p. Linear period-to-
Vol. 8, October
(ii) B + S + BS
Sensitivity analysis of oscillatory systems: M. A. Kramer et al.
(iii) AS + AS 3 CS + S k, (iv) BS + CS -+ P + 2s
8.0
where S are catalytic surface sites and P is the desired product. The adsorption and desorption of A, B, and P are assumed to be very fast compared with the rate of reaction (iii). Given the reaction rate constants, adsorption equilibrium constants K*and KB, and the total molar flow rate per unit reactor volume, F, it is desired to maximize to yield of P, by varying only the feed composition. Using KA= KB = k4 = F = 1, and k3 = 0.1, the optimum fraction of A in the feed is Af= 0.812, giving P = 0.01158, in the steady state. This state may be calculated directly by coupling equations (1) and (2) with CX~ = Af, setting the entire LHS and aP/aA, to zero, and solving the resultant set of eight algebraic equations in the eight unknowns, A, B, CS, P, aA/i3Af, aB/aA,, XS/all,, and Af Sensitivity analysis was then undertaken to determine if the optimum steady state process is indeed the optimum process. Previous work approached this question by calculating the second variation about the optimal steady state by asymptotic methods in the limits of high or low frequency forcing,” or by the so-called ‘pi-criterion’ test of Bittanti et a1.19 Alternatively, perturbation approaches have been used (e.g. Skeirek and Devera),” and these methods most closely resemble the sensitivity treatment, particularly in their ability to handle any cycling frequency. Horn and LinzO showed that the first variations in the periodic variables around an optimum steady state process are zero, and therefore as a first step in synthesizing a periodic control, the second order sensitivities a2P/a(u2 must be evaluated. In the current process, the controllable parameters are the fractions of A and B in the feed stream, and the cycling frequency w = 21r/r. Periodic flows of the type Af(t) = Af( 1 + fA sin wt) and Bf(t) = Bf( 1 + en sin(ot + 6)) were
7.0
6.0
5.0
Y2 4.0
3.0
2.0
1.0
0.0
0.5
0.6
0.7
0.8
0.9
1.0
YI figure 9 Actual and first-order predicted trajectories with Ap = 0.05 using path-independent sensitivities. Actual trajectory indicated by dashed line
Table 1
Exact and estimated
periods for a CSTR
in cases considered by Hugo and Wirgesls
Parameters
0.8 0.8 0.6 0.8 0.75 0.6 0.55 aSensitivity
Table 2
Period, r
E
B
IJ
Exact
Linear theoryI
Singular perturbation”
First ordera sensitivity
Second ordera sensitivity
0 0 0 0.02075 0.02075 0.02075 0.02075
10 10 10 30 30 30 30
0.280 0.2509 0.261 0.3118 0.310 0.3038 0.2716
2.611 3.278 6.698 1.267 1.453 1.989 2.625
2.566 2.967 11.20 1.823 3.050 Imaginary Imaginary
_ 0.48 0.56 0.76 1.02
2.818 3.308 4.029 0.7045 0.9576 1.731 2.496
3.225 4.032 5.064 1 .a21 2.004 2.285 2.618
results based on Taylor
Exact and estimated
series with nominal parameter
modified
amplitudes
for CSTR
values us = 0.7, E = 0.010375,
in cases considered Modified
Parameters
amplitude,
b = 20, p = 0.284
by Hugo and Wirges” a’
us
E
B
fi
Exact
Empirical approximation12
Singular perturbation”
First order sensitivitya
Second order sensitivitya
0.8 0.6 0.8 0.7 0.65
0 0 0.02075 0.02075 0.02075
10 IO 30 30 30
0.297 0.261 0.573 0.467 0.445
0.275 2.21 1.19 3.82 4.45
1.22 1.79 1.44 3.84 4.33
2.00 3.20 2.83 4.10 4.61
2.16 3.64 2.77 5.71 4.62
0.89 2.23 1.37 4.39 4.02
‘Sensitivity
results based on Taylor
series with nominal parameter
values us = 0.7, E = 0.010375,
Appl.
Math.
b = 20, p = 0.417
Modelling,
1984,
Vol.
8, October
337
Sensitivity
analysis of oscillatory
systems: M. A. Kramer et al.
examined, where d;f+ Bf = 1. Taking two derivatives of P, the average of P over one period r. one obtains:
a2P -= ad
Qa7ap 7 Gas aff2 +,
1 i
(t + 7) + 1 (P(t + 7) - P) ; r
ar 2dP G --(If7)--1
2
ap a7
r
aa aa
(32)
Evaluated at eA = fu = 0, only the first two terms of this expression are non-zero, and only the first is non-zero if (Yf w. The full expression, obtained by Liebnitz’s rule for differentiation of integrals, is actually the time average of the second order path-independent sensitivities, which reinforces the concept that these quantities are the true indicator of the parametric effects on the shape of the equilibrium orbit. However, with the aforementioned simplifications and because ar/acr is either -r/w or zero, depending on whether or not (Y= o, direct use of the elementary sensitivities is straightforward. The values of:
calculated by the AIM method4” are listed in Table 3 at various w. For this problem, both a’P/&iand a’P/aei were vanishingly small, although this is not the case generally. The AIM program, incidentally, provided an extremely efficient method of determining these sensitivities, for two reasons. First, since the sensitivities were evaluated around a steady state, the Jacobian was constant, and thus it only had to be diagonalized once. Diagonalization is usually a major part of the computational effort in AIM. Secondly, various values of w and 6 were handled in a single simulation, by considering each numerical value to be a distinct parameter. Using AIM, the marginal cost of the sensitivity to a new parameter is very small, and in fact the computational cost of all the data in Table 3 was only about twice that of a single kinetic simulation. Positive values of a2P/ &A a+, were found for 6 = 7~,Af = 0.812, in a range of frequencies from about 0.2 to 50, with a maximum of +0.0007 at w = 0.50. Extrapolation using a second order Taylor series indicates an improvement of about 6% over the steady optimum value with this frequency, and eA = es= 1. A simulation was performed to determine whether a periodic state existed at these parameters values,‘l and the true extent of yield augmentation. A limit cycle of period r = 2n/c; was indeed present, and a 5.7% increase in P was
Tab/e 3 Second order sensitivitie? state of example
around the optimal
steady
w
a*iaeA
1o-4 10-S 10-Z 10-l 1 10 IO2 103 lo4
-0.007 -0.007 -0.007 -0.002 + 0.005 +0.001 -0.00004 -0 -0
aEvaluated
338
Appl.
at e,q = EB = 0, S = rr, &
Math.
Modelling,
= 0.812
1984,
Vol.
8, October
acg
calculated. The question now arises as to whether this periodic mode is the most beneficial. The easiest way to answer this question is to evaluate the normalized sensitivity of P to the various parameters, on this limit cycle. This calculation gave a In P/a hr eA = t0.8, a In P/ ~ln~~~+0.05,~~ln~/~ll~~~<~O~4,~~~n~/~ln~~< 10d4, and I a In P/a In 6;,1 < 1Oe4. The first two sensitivities suggest further increases in eA and en, however eA = eu = 1 are the maximum permissible values. The final three sensitivities are approximately zero, indicating local maxima. Therefore, this is the optimum strategy for sinusoidal feed variations. A physical interpretation of this result is that the maximum segregation of A and B provided by 6 = 7~promotes a high concentration of AS, which stimulates the slow reaction (iii). during the first half of the period, whereas during the second half period, B is introduced in quantity to react with CS to form the product. This amounts to serially running reactions (iii) and (iv) while providing maximum surface coverage of the respective reactants by removing the competition between A and B for surface sites. This suggests that square wave feed profiles, providing further segregation of A and B, may be even more advantageous than sinusoids, and indeed, an increase to P = 0.01289 was obtained using bang-bang (on-off) control of A and B, using the same Af, Bf, and w as previously. Further simulations indicated that this was the optimum control of its type. The advantages of the sensitivity approach in synthesizing periodic operating strategies are two-fold. First, a high degree of automation is introduced, with the potential to scan through a large number of possible control schemes at relatively low cost. without restriction on the complexity of the system. Second, and most importantly, one is not restricted to the immediate neighbourhood of the optimum steady state in making sensitivity estimates. The first estimate of the value of a periodic scheme, provided by the second derivatives around the steady state, may be followed by further simulations and sensitivity calculations on the limit cycle generated by the predicted best combination of parameters. Thus, an iterative scheme may be constructed to truly optimize the periodic control. In addition, it may be possible to use the methods presented earlier in this work to further the understanding of the physicochemical processes controlling the forced oscillation, and use this information to construct waveforms and strategies better suited to the process in question.
Conclusion The techniques of sensitivity analysis have been extended to treat systems which exhibit limit cycle behaviour. A method was developed to separate the raw sensitivity coefficients into path-dependent and path-independent portions. The path-dependent data was quantified in terms of the parametric dependencies of the phase of the system, and could be used to determine the sensitivity of the cycle period, and the sensitivity of the time of occurrence of features such as extrema. The path-independent sensitivities were shown to contain information on perturbative dynamics and the shape of the equilibrium orbit. Each of these analyses were developed to the second order. A simple example involving the limit cycle behaviour of a well-mixed chemical reactor was presented in which these methods were used to predict parametric sensitivities of limit cycle phase, dynamics, shape, amplitude, and period.
Sensitivity
Additionally, it was shown that the sensitivity approach can be used to design beneficial periodic operating strategies. Another method suggested as a possible approach to these problems, known as feature sensitivity analysis,22 deserves mention at this point. In this procedure, one finds a functional ‘fit’ for each oscillatory trajectory, in terms of a set of new parameters p, viz., y = g(p), where each of these parameters has a specific significance in terms of the shape or dimension of a particular ‘feature’ of the trajectory y. The partial derivaties @/acr quantitatively link the parameters of the original system to the features of the output y. These quantities are determined by solving a set of linear algebraic equations involving the sensitivities ay/ap and the conventional elementary sensitivities. The utility of this approach seems to be limited by the ability to find an appropriate function g which not only accurately fits the oscillations, but also has a clear physicality in terms of its parametric dependencies. Unless there is an a priori reason for selecting a particular functional form for g, finding a suitable function may require significant effort, and the current approach is probably preferable. The sensitivity analysis approach provides modellers with a method of probing the detailed relationship between the various aspects of limit cycle behaviour and the system parameters, in a quantitative fashion. Such information can be of significant value in both theoretical and experimental studies involved in limit cycle characterization, or in determination of the phenomena responsible for oscillatory behaviour. Furthermore, sensitivity analysis provides a natural method of developing correlations for limit cycle period and amplitude, with an economy of effort. Significantly, the method is not limited by the degree of nonlinearity in the model, and does not depend on any assumption of near-harmonic or small frequency oscillations. The sensitivity approach shows great promise as synthesis method for periodic operation strategies, and work continues on this problem.
analysis
of oscillatory
systems:
M. A. Kramer
et al.
Fife, P. C.J. Phys. Chem. 1981,85, 2861 Larter, R., Rabitz, H. and Kramer, M. A, ‘Sensitivity analysis of limit cycles with applications to the Brusselator’,J. Chem. Phys. 1984,80,4120 13b Larter, R. ‘Sensitivity analysis of autonomous oscillators: separation of secular terms and determination of structural stabiJity’,J. Phys. Chem. 1983, 87, 3114 14 Edelson. D. and Thomas. V. M. J. Phvs. Chem. 1981.85. 1555 ’ 15 Hugo, P. and Wirges, H. P. 5th Int. Symp. Chem. Reaction Eng., Houston, 1978, p. 498 16 Uppal, M., Ray, W. H. and Poore, M. B., Chem. Eng. Sci. 12 13a
1914,29,961 17
18 19 20 21
22
Wicke, E. ‘Grundlgen der Chemische Prozessregelung’ (Eds Oppelt and Wicke), Oldenbourge Verlag, Miinchen, 1964, p. 54 Bailey, J. E. ‘Periodic optimization’ (Ed. A. Marzollo), Springer-Verlag/CISM, Wien, 1972 Bittanti, S., Gronza, C. and Guardabassi, G. IEEE Trans.
Auto. Control 1973. AC-18. 33 Horn, F. J. M. and iin, R. d. Ind. Eng. Chem. Des. Dev. 1967, 6, 21 Periodic states are not guaranteed, and theoretical proof of existence is problematic. See Sincic, D. and Bailey, J. E., Chem. Eng. Sci. 1977, 32,281 Skumanich, M. and Rabitz, H., Commun. J. Molecular Sci.
1982,2,79
Appendix Derivation
of selected
second order expressions
The detailed derivations of several of the expressions used in the paper are presented below. Equation
(23). We begin with equation
and differentiate
(9), written as:
with respect to &k at constant
r$:
Acknowledgement
(AlI
The authors are grateful to the Office of Naval Research for support of this work.
References Tomovic, R. and Vukobratovic, M. ‘General sensitivity theory’, American Elsevier, New York, 1972 Dickinson, R. P. and Gelinas, R. J. J Comp. Phys. 1976, 21, 123 Hwang, J. T., Dougherty, E. P., Rabitz, S. and Rabitz, H. J.
Chem. Phys. 1978,69,5 Kramer,
180
M. A., Calo, J. M. and Rabitz,
8 9
Kramer, M. A., Calo, J. M., Rabitz, H. and Kee, R. J. ‘AIM: the analytically integrated Magnus method for linear and second order sensitivity coefficients’, Sandia National Laboratories report 1982, SAND82-8231 Demiralp, M. and Rabitz, H.J. Chem. Phys. 1981,75, 1810 Gilles, E. D. ‘Gundlagen der Chemische Prozessrogelung’ (Eds Oppelt and Wicke), Oldenbourg Verlag, Miinchen, 1964, p. 89 Beck, J. Chem. Eng. Set. 1977, 32, 159 Douglas, J. M. and Caitondo, N. Y. Jnd. Chem. Fund. 1967,6,
265 10 11
--a
aYi
i affk a+
) =“‘+(_EJLJ%) G
aai ask
I
H. Appl. Math. Model-
ling 1981, 5, 432
6 7
To simplify this expression, it must be realized that yi is a dummy variable in equation (9) and that equation (9) is the general definition of partial differentiation at constant $. Applying equation (9) to the first term on the RHS of equation (Al), we obtain:
Chang, H. C. and Calo, J. M. Chem. Eng. Cammun. 1979,3, 431 Skeirek. R. D. and Devera, A. L. Chem. Eng. Sci. 1982, 37, 1015
The second term may be similarly evaluated:
Substituting
into equation
(Al), we obtain equation (23).
Equation (25). We begin with equation (23), evaluated at t + T. Several terms arise which will need to be simplified, namely :
Appl.
Math.
Modelling,
1984,
Vol. 8, October
339
Sensitivity analysis of oscillatory systems: M. A. Kramer et al.
which, after simplification,
and
i($f+?)) =-g$&Z)
g
(t + 7) i
gives:
&
k
(t) k
I
=$(2 (t) j-g dz(t) In the following, all terms are evaluated at time t, unless otherwise indicated. From equation (23), we have:
ar
&r d2yi
dyi
+aol_cdt2---I k
a2yi (t + 7) aai ask
a27
dt affj aok
Solving for a2r/aoj aok we obtain equation
(25).
Equation (29). Begin with an expression analogous to equation (9), with ‘e’ appearing in place of 4, indicating decomposition in the direction of change of the extrema:
avi _ ayi~at,dYi
0aai e
Differentiate
Subtracting a*Yi aaj
(t
equation +
r)
&Yj at
(-42)
with respect to ok:
(23) evaluated at time t, we obtain:
a33 aai ask
ask
aai
+dyi
(‘)
a2te
dt sol, aok The last term in this expression is zero because dyi(tJ dt = 0. The second-to-last term is also zero, because along e (the locus of all phase-space points through which the extreme passes as a function of Q) dYi/dt is by definition zero. The sole remaining term can be simplified using equation (A2): (A3)
dyi
a27
dt
aaiask
Appl.
Math.
--p
340
Modelling,
Applying equation (28), it can easily be seen that equation (A3) is the same as equation (29).
1984,
Vol. 8, October