12/7/2020

Ali M. Sahlodin
Department of Chemical Engineering
AmirKabir University of Technology
1397 S.H

 BDF Solvers
 Variable order, variable step size
 Different differentiation methods
 Finite differences
 Symbolic differentiation
 Automatic differentiation

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 2

1
 12/7/2020

 Parametric ODE  t )  f(x(t ;p),p, t ) = 0

x(
x(t0 ;p)  h(p)

 Parametric DAE f(x, x,p,

 t) = 0

Initial conditions can be a function of parameters

 Parameters (p) can be

 Process constants (e.g., kinetic constants, heat
transfer coefficients)
 Optimization parameters (e.g., flow rate)

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 3

 What is sensitivity?
 Partial derivatives of functions w.r.t. parameters

 Why is it needed?
 For optimization: (gradient-based) optimization
algorithms need to know when the derivatives
become zero.
f ( p)

df *
(p )  0
dp
p*
Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 4

2
 12/7/2020

 Predict system’s behavior for different parameter

values without simulating it again.

≈F(P*+ΔP)

P* P*+ΔP
How much would the reactor yield change if the reaction volume
changed by 5%?

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 5

ODE systems  t;p)  f (t , x(t;p), p), x(t0 ;p)  h(p)

x(

Initial conditions are functions of parameters in general

Some function of interest J (p) :  (x(t f );p)

dJ
Desired: sensitivity (p) at a certain parameter
dp
value

Parameters affect J(p) through the differential equations

(this complicates sensitivity calculations).

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 6

3
 12/7/2020

 Given first-order reactions ABC, find the

kinetic constants k1 and k2.

Minimize error between the model

and the data
Optimal k1 and k2

This is a dynamic optimization

problem because we have a
dynamic model.

Sensitivities (derivatives) of the objective function w.r.t

k1 and k2 are needed.
Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 7

 Use 1st-order forward difference for each

parameter

J J ( p1 , p2 ,, p j   p j ,, pn p )  J (p)

(p) 
p j  pj

 Integrate ODEs
 t ; p)  f (t , x(t; p), p)
x(
J (p)   (x(t f );p)
x(t0 ; p)  h(p)

 How to obtain J ( p1, p1,, p j   p j ,, pn ) ? p

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 8

4
 12/7/2020

 For each parameter, integrate ODEs with p j   p j

 Keep all other parameters unchanged.
 Then, evaluate J ( p1, p2 ,, p j   p j ,, pn p )

np+1 integrations are needed to get all sensitivities

(using first-order finite differences)

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 9

 Advantages
 Easy to implement and calculate.
 Mathematical dynamic model is not needed.

pj  pj Black box dynamic model

(e.g., Process Simulators)
J ( p1,, p j   p j ,, pn p )
 Disadvantages
 Inaccuracy: finite differences
 Inefficiency: many integrations are needed
 First-order finite difference: np+1

 Second-order finite difference: 2np+1

 Difficult to find a good perturbation value

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 10

5
 12/7/2020

 Perturbation value
 Too small: accuracy affected by
rounding errors
 Too big:
 Accuracy affected by truncation Model not
errors valid here

 Perturbation may push the

p p+δp
model outside its domain of
validity  Wrong derivatives
Direction of perturbation may Use p-δp instead!
also be important: p+δp or p-δp?

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 11

 t;p)  f (t , x(t;p), p),

Dynamic system (ODE) x( x(t0 ;p)  h(p)
Given function J(p) :  (x(t f );p)
1) Construct forward sensitivity equations for each parameter
x p j (t ; p)  f x (t , x(t; p), p)x p j (t; p)  f p j (t , x(t; p), p)
x p j (t0 ; p)  h p j (p)

Linear, non-homogenous ODEs (why?)

2) Integrate the sensitivity equations until final time.

• Augment the original ODEs with all sensitivity equations

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 12

6
 12/7/2020

Augmented system of ODEs

 t ;p)  f (t , x(t;p), p), x(t0 ;p)  h(p)
x(
x p (t ; p)  f x (t , x(t; p), p)x p (t; p)
j j

 f p j (t , x(t ; p), p), x p j (t0 ; p)  h p j (p), j  1, , n p

3) Obtain sensitivities of objective function.

dJ
(p) : x (x(t f );p)T x p j   p j (x(t f );p)
dp j

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 13

 Size of augmented ODEs (np+1)nx.

 Very large ODE system when np and nx are
large
 computationally inefficient
 Partial derivatives fx fp should be obtained
from symbolic or automatic differentiation.
 Avoid finite differences where possible
 Forward sensitivities for DAEs

f(x, x,p,
 t) = 0 f x x p j  f x x p j  f p j = 0

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 14