DIANA A short introduction Michael Krasnyk Max Planck Institute - - PowerPoint PPT Presentation

11 March 2009 ICVT, Universit¨ at Stuttgart

MAX−PLANCK−INSTITUT TECHNISCHER SYSTEME MAGDEBURG DYNAMIK KOMPLEXER

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

DIANA — A short introduction

Michael Krasnyk

Max Planck Institute for Dynamics of Complex Technical Systems, PSPD group Otto-von-Guericke-University, IFAT

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Motivation

High demand for first principle modeling of chemical processes

Complexity of processes and models Structure of implemented models

Goals of computer-based modeling

computer-aided process engineering systematical modeling approaches reusable and transparent models

Physically motivated concepts for structuring of balance based models Ponton (1991), Marquardt (1996), Gilles (1997), Mangold (2002)

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 2/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Modeling tool ProMoT

Equation-based modeling

differential-algebraic systems (differential index 1) Petri networks

• bject-oriented concepts for models

multiple inheritance aggregation

model implementation

text-based in modeling language MDL graphical modeling with GUI

Equation analysis and optimization

elimination of explicit algebraic relations structure analysis of the complete system for solvability

Symbolic differentiation in ProMoT

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 3/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Simulation tool Diana

Promot Diana Python

C++ model (define-module :class "Test" ... ) SWIG CapeDAESO CapeDAESOSolver PetriNetworks Continuation import diana import solver mod=Model("Test") ...

Data

Model description Initial values Simulation results

User

XML Make Create Run

GUI

PyGTK NumPy SciPy

Current activity Future plans

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 4/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Diana architecture

C++ Model

• encapsulates model data:
• Parameters
• State Variables
• Provides interface to
• Equations/Help-Vars
• Event Functions
• Petri network
• Assertions

Model Loader

• loads a model from a file
• calls model-specific initialization

routines

Result Files Solver Factory

• loads numerical algorithms

from file

Diana Numerical Algorithms Integrators Solvers Event Handlers Optimizers Continuators Sensitivity Analysers Parameter Analysers Numerical Libraries LAPACK BLAS SUNDIALS Umfpack Harvell ARPACK Python Scripting

• Control of simulator activity
• Model Data Access
• Callback Functions

SciPy/ Numeric

• Scientific Programming in

Python

• Matlab-like Functionality

Python Python-embedded shared libraries Swig-Wrapping

Online Plotter Gnuplot / wxPython

• Online / Offline Plots
• Selection of Variables
• Generation of standard Plots

Output Service

• collection of data during

calculation

• Access by scripts/plotting
• Output to Simulation-Log

Model Files

Promot OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 5/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Diana dependencies

Required dependencies

gcc (version 3.3.1) python (version 2.4) swig (version 1.3.31) xerces (version 2.7.0) cppunit (version 1.11.0)

External numerical libraries:

Required linear algebra libraries BLAS http://www.netlib.org/blas/ LAPACK http://www.netlib.org/lapack/ UMFPACK http://www.cise.ufl.edu/research/sparse/umfpack/ Optional linear algebra libraries ARPACK http://www.cse.scitech.ac.uk/nag/hsl/ Harwell Subroutine Library http://www.cse.scitech.ac.uk/nag/hsl/ Differential algebraic solvers IDA/Sundials http://acts.nersc.gov/sundials DASPK http://www.cs.ucsb.edu/~cse/software.html

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 6/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Diana installation

Latest sources can be found here: http://promottrac.mpi-magdeburg.mpg.de/dist/ Doxygen-generated documentation for the Diana: http://promottrac.mpi-magdeburg.mpg.de/doc/Diana/ Build and install:

./ config/bootstrap ./ configure

• -prefix=< installation

prefix > \

• -with -ufsparse=<ufsparse

installation dir > \

• -with -arpack=<arpack

lib path > \

• -with -blas=<blas

lib path > \

• -with -lapack=<lapack

lib path > \

• -with -sundials=<sundials

installation dir > \

• -with -cppunit -prefix=<cppunit

installation dir > \

• -with -xerces=<xerces

installation dir > \

• -with -mdl2diana=<mdl2diana

path > make make check make install

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 7/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Simulation models

CAPE-OPEN interface extension IDianaDAESO defines a dynamical model as a DAE system in implicit form f (t, x, ˙ x, ν) = 0, f : R × Rn × Rn × Rp → Rn ∈ C∞ extension interface gives an access to derivatives of the model ∂k+l+m f (t, x, ˙ x, ν) ∂kx ∂l ˙ x ∂mν , k, l, m 0, that are obtained with help of CAS Maxima ProMoT command mdl2diana produces C++ model code

mdl2diana <module > [-f <mdl -file >] [-d <diana -name >] [-c] [-g <generation -dir >] [-sd <triples

• f integers >]

[-sd -sing] [-sd -sens]

Diana command dianac compiles a model to shared library

dianac <model_name > [--clean] [--rebuild ] [-O <OPT >]

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 8/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Simulation model interface

Interface methods

Get/Set(All)[Variables|Derivatives|Parameters] access methods to the state vector x, derivatives vector ˙ x and parameters vector ν. Get(All)Residuals methods return the residual vector f Get[Lower|Upper]Bounds methods return user defined maximal and minimal values for the state variables Get/SetIndependentVar are access methods to the independent variable t Get(All)(Diff|Par)JacobianValues methods return values of the Jacobian matrices ∂f /∂x, ∂f /∂ ˙ x or ∂f /∂ν GetHighOrderJacobian method returns values of the higher order Jacobian ma- trix ∂(k)f /∂{x, ˙ x, ν}(k) Save/LoadState save and restore a state of the ESO instance

ESO additional parameters

SymbolicJacobian shows whether a symbolically or numerically computed deriva- tives will be returned by Get(*)JacobianValues FDPartition, FDOrder, FDEpsilon numerical differentiation parameters

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 9/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Model: Continuous Stirred Tank Reactor

cin ˙ qin c T, c

The mass balances

• f

the model read [Zeyer et al., 1999] ˙ cH2O2 = ˙ qin/V (cH2O2,in − cH2O2) − (r1 + r2 + r3) ˙ cCH3CHO = ˙ qin/V (cCH3CHO,in − cCH3CHO) + (r1 − r2) ˙ cCH3COOH = ˙ qin/V (cCH3COOH,in − cCH3COOH) + r2 ˙ ccat = ˙ qin/V (ccat,in − ccat) − (r4 − r5) Reaction rates ri, i = 1, . . . , 5 are 2 6 6 6 6 4 k1e−E1/(RT)ccat cH2O2 k2e−E2/(RT)ccat cH2O2 cCH3CHO k3e−E3/(RT)ccat cH2O2 k4e−E4/(RT)ccat √cCH3CHO k5e−E5/(RT)(cF,ges − ccat) 3 7 7 7 7 5 The energy balance is V ρcp ˙ T = ρcp ˙ qin(Tin − T) + (UA)cool(Tcool − T) + V P3

i=1 ri(−∆hR)i

Vcoolρcp ˙ Tcool = ρcp ˙ qcool(Tcool,in − Tcool) + (UA)cool(T − Tcool)

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 10/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

MDL simulation model

( define-module :class" Hafke_Reactor" : super-classes("module" ) : documentation " Ruehrkesselreaktor zur Oxidation von Ethanol zu Essigsaeure" :variables ( ("NC" : documentation "Anzahl der (C) Komponenten" : system-theoretic " structure-parameter" :value "4") ("i" : system-theoretic "index") ("Tkzu" : documentation " Zulauftemperatur Kuehlmittel [K]" : system-theoretic " real-parameter" :value "314.0") ("qknormal " : documentation " Zulaufstrom Kuehlmittel [m^3/s]" : system-theoretic " real-parameter" :value "4.44 d-05" ) ... ("c_" : documentation " Konzentration der Komponenten P E S F [mol/m^3]" :is-a "variable " : system-theoretic "state" :indices ((: index "i" :lower "1" :upper "NC")) :value "0.1" :minimum "0.0" :maximum "3000.0" : absolute-error "1.0 D-09") ("Temp" : documentation "Temperatur Reaktionsgemisch [K]" ... ) ("Tk" : documentation "Temperatur Kuehlmittel [K]" ... ) ("R_reak" : documentation " Reaktionsgeschwindigkeit der Teilreaktion " :indices ((: index "i" :lower "1" :upper "NR")) :is-a "variable " : system-theoretic "help" :value "k_reak[i] * :cond(i==1, cP_cF , i==2, cP_cF * c_[2], i==3, cP_cF , i==4, c_[4] * :sqrt(c_[3]) , c_F_ges - c_ [4])") ... ) :equations ( ("KMB_Eq" :is-a "equation " :indices ((: index "i" :lower "1" :upper "NC")) :relation ":diff <t> c_[i] == qVR * (c_zu[i] - c_[i]) + RCumsatz [i]") ("Temp_Eq " :is-a "equation" :relation ":diff <t> Temp == qVR*( Tzu - Temp) + kappa/VR/rho_cp *(Tk - Temp) + :sum <i,1,NSUM >( R_reak[i]* DHR[i])/ rho_cp") ("Tk_Eq" :is-a "equation " :relation ":diff <t> Tk == qVk *( Tkzu - Tk) + kappa/Vk/rho_cp *( Temp - Tk)" ) ) ) OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 11/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Simulation models in Python

Diana interpreter is launched with

diana [<python arguments >] [<script name >] [<arguments >]

Diana module import and initialization

import diana main=diana. GetDianaMain ()

model loaded as a shared library

mm=main. GetModelManager () model=mm. CreateModel (CAPE_CONTINUOUS , " HafkeReactor .so")

acquiring equation set object (ESO) with references to states and parameters

eso=model. GetActiveESO () evar=eso. GetStateVariables () epar=eso. GetRealParameters ()

access to ESO data

• eso. GetAllVariables ()

epar[’tkzu ’]. GetValue () evar[’tk’]. SetValue (300.0)

• eso. GetAllJacobianValues ()

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 12/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Dynamic simulation

Dynamic solvers find solutions of the Cauchy problem ϕ(t, x0, ν), such that f (t, ϕ, ˙ ϕ, ν) ≡ 0, s.t. ϕ(t0, x0, ν) = x0 Dynamic simulation is presented by the following solvers IDASolver and DASPKSolver for implicit DAE systems with differential index 1 (libraries ida.so and daspk.so) OdessaSolver for ODE systems with ∂f /∂ ˙ x = I (library odessa.so) The linear systems in the integrators are solved by the direct dense LAPACK or sparse UMFPACK linear algebra solvers For integrating a DAE initial-value problem, an important requirement is that the pair of vectors x0 and ˙ x0 are both initialized to satisfy the DAE residual f (t0, x0, ˙ x0, ν) = 0 For semi-explicit differential index-one systems, IDA provides a routine that computes consistent initial conditions from a user’s initial guess

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 13/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Dynamic solver interface

Interface methods

GetParameters method returns collections of solver parameters Solve method starts the solution of an ESO that is associated with the solver GetSolution method returns the solution vector x

Solver parameters

Start parameter specifies the start of a new simulation (True) or continuation

• f the current one (False)

CalcIC parameter controls whether consistent initial conditions are computed at the initial time (True) or not (False) T is the current value of the independent variable t T0 is the starting value of the independent variable Tend is the final value of the independent variable Intermediate parameter with the True value tells the solver to take one internal step and to return the solution at the point reached by that step, otherwise integration proceeds to the parameter value Tend without interruption VerboseLevel controls a verbosity level of the solver LASolver specifies whether the dense or sparse linear algebra solver will be used

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 14/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Dynamic simulation in Python

loading of IDA solver form SUNDIALS package

sf=main. GetSolverFactory () solver=sf. CreateSolver (CAPE_DAE , model , ’ida.so’)

• solver. Initialize ()

spar=solver. GetParameters ()

initializing of a reporting interface and linking to the solver

ri=main. CreateReportingInterface (’basic ’)

• solver. SetReportingInterface (ri)

ri.Add(spar[’T’]) ri.Add(epar[’qknormal ’])

setting of solver parameters

spar[’VerboseLevel ’]. SetValue (2) spar[’T0’]. SetValue (0) spar[’Tend ’]. SetValue (1000.0)

performing of a simulation and storing of output data to ’out.m’

epar[’qknormal ’]. SetValue (1e -4) solver.Solve ()

• ri. WriteDataMatlab (’DynamicSimulation /Example .m)

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 15/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Dynamic simulation

plotting of the dynamic simulation in Matlab

run(’DynamicSimulation /Example .m’) figure; hold on; box on; grid on; plot(data1 (:,1), data1 (: ,2)) figure; hold on; box on; grid on; plot(data1 (:,1), data1 (: ,7))

results of the tutorial case tutorial/hafke/simulation

1 2 3 4 5 6 ×104 0.5 1 1.5×10−4 t, [s] qcool, [m3/s] 1 2 3 4 5 6 ×104 300 320 340 360 380 400 420 t, [s]

• T. [K]

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 16/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Parameter continuation and nonlinear analysis

Parameter continuation is presented by the following continuation of steady-state points and analysis of local stability (library sstate.so) continuation of limit point curves and analysis of singularities of such points (library sanalyser.so) continuation of Hopf points (library hopf.so) continuation of periodic solutions and analysis of local stability based on Flo- quet multipliers (library periodic.so)

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 17/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Continuation solver interface

Interface methods

Solve method solves a nonlinear task f (x, ν) = 0 for constant parameter ν Continuate method performs a continuation with respect to parameter ν Add/RemoveFreeParameter method adds or removes λ to the nonlinear system

Continuation solver parameters

Parametrization parameter specifies the parametrization type (PseudoArclength

• r Local)

Predictor is the predictor type (Tangent or Chord) StepSize is the current step size σ(k) InitialStepSize is the initial step size σ(1) InitialDirection is the initial direction of a continuation MinStepSize is the minimal step size σmin MaxStepSize The is the maximal step size σmax MaxStepsAmount is the maximal number of steps kmax Tol relative tolerance in the argument space

Steady-state continuation solver parameters

StabilityCheck, Stability parameters specify local stability check and stability

• f the current point

ConditionCheck check for limit point SteadyStateZCE or Hopf point SteadyStateZRE conditions

Singularity analysis solver parameters

ConditionCheck check for zeros of test functions (SingularityG[x|xx|p|xp]) Gx,Gxx,Gp,Gxp values of test function derivatives ConditionEquations set of adjoint test function equations

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 18/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Parameter continuation of steady-states

create model with higher order derivatives (mdl2diana option -sd-sing) loading of continuation solver

sf=main. GetSolverFactory () conti=sf. CreateSolver (diana.CAPE_CONTI , model , "sstate.so") conti.Initialize () cpar=conti. GetParameters ()

initializing of a reporting interface and linking to the solver setting of solver parameters

cpar[’VerboseLevel ’]. SetValue (0) cpar[’MaxStepsAmount ’]. SetValue (5000) cpar[’MaxStepSize ’]. SetValue (5.0) cpar[’StabilityCheck ’]. SetValue (True) cpar[’ConditionCheck ’]. SetValue (diana. SingularityNone )

adding of a continuation parameter

• conti. AddFreeParameter ("tkzu", 250.0 ,

400.0)

performing of a continuation

conti.Continuate ()

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 19/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Parameter continuation of steady-states

Results of the tutorial case tutorial/hafke/continuation

285 290 295 300 305 310 290 300 310 320 330 340 350 360 Tcool,in [K] T [K] 290 292 294 296 298 300 30 35 40 45 50 55 60 65 Tcool,in [K] ˙ qcool [l/h] 0.5 1 1.5 2 2.5 3 3.5 300 305 310 315 320 325 330 335 ˙ qcool [l/h] T [K] 4.6 4.8 5 5.2 5.4 5.6 5.8×10−7 310 312 314 316 318 ˙ qcool [l/h] T [K]

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 20/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Optimization Problems Supported in Diana

Optimization of explicit functions Parameter Estimation Maximum likelihood approach for Gaussian white noise Optimal Experimental Design Sigma point approach

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 21/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Global Nonlinear Continuous Optimization Packages

Global Random Search: Genetic Quite flexible implementation of Genetic Algorithms with self-tuning of ran- domization parameters Statistical Global Optimization: DIRECT v2.0 Dividing rectangular global optimization method by D. R. Jones and J. Gablon- sky BBOWDA Black box global optimization method with data analysis by K. Kofler GMFL

routine bayes1 Bayesian global optimization method routine unt The global method of extrapolation type by A. Zilinskas routine lbayes The local Bayesian method by J. Mockus

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 22/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Local Nonlinear Continuous Optimization Packages

Non-gradient Methods routine NMSimplex Nelder-Mead downhill simplex method implementation by D. E. Shaw Gradient Based Methods routine L-BFGS-B Limited-memory quasi-Newton code for large-scale bound-constrained or un- constrained optimization by C. Zhu and J. Nocedal. IPOpt Package for large-scale nonlinear optimization of continuous systems, imple- ments a primal-dual interior point method, and uses line searches based on Filter methods and Hessian approximation using BFGS update.

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 23/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Capabilities of optimizers in Diana

Feature SO-library genetic direct bbowda bayes1 unt lbayes nmsimplex lbfgsb ipopt Global search X X X X X Stochastic objectives X X X X X X Requirements Need of gradient X X Support of constraints handling bound constraints X X X X X X X X direct constraints X X inequality constraints X X X X equality constraints X X

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 24/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Running Parameter Estimation

Create optimizer

sfactory = main. GetSolverFactory () solverPE = sf. CreateSolver (diana.CAPE_NLP , None , ’direct ’)

Load measured data from specified file (details [Schenkendorf et al., 2009])

md = DianaMeasuredData (model , ["c_x", "c_s"]) md.load(’./ observations .dat ’)

Specify estimated parameters

sps = [

• diana. DianaNLPRealParameterSpec ("k_s", "k_s", 1., 1., 3.),
• diana. DianaNLPRealParameterSpec ("mu_max", "mu_max",

4., 4., 6.)]

Create parameter fitting task

taskPE = diana. ParameterFittingTask (main , ’ida.so’, model , md , sps)

Assign task to optimizer and run optimization

solverPE . SetReportingInterface (report ); solverPE . SetNLPTask (taskPE ); solverPE . Initialize (); solverPE .Solve ()

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 25/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Running Optimal Experimental Design

Describe design variables

colDesignVars = [

• diana. DianaNLPRealParameterSpec ("q_0", "q_0", 0.09 , 0.07 ,

0.08) ]

Create OED task

taskOED = nlptaskfactory . CreateOEDTask (taskPE , solverPE , ... colDesignVars , ’oedsigmapoint ’)

• edpar=taskOED . GetParameters ()

# use E* optimality criterion

• edpar[" OptimalityCriterion "]. SetValue (diana. DIANA_OED_EStar )

taskOED .Initialize ()

Assign task to optimizer and run optimization

solverOED . SetReportingInterface (report ); solverOED .SetNLPTask (taskOED ); solverOED .Initialize (); solverOED .Solve ()

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 26/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Contact

Promot/Diana homepage http://www.mpi-magdeburg.mpg.de/projects/promot/ ProMoT and network theory questions Michael Mangold mangold@mpi-magdeburg.mpg.de ProMoT and GUI questions Sebastian Mirschel mirschel@mpi-magdeburg.mpg.de Optimization in Diana Sergiy Gogolenko gogolenk@mpi-magdeburg.mpg.de Diana questions Michael Krasnyk miha@mpi-magdeburg.mpg.de Thank you for your attention!

OvGU, MPI Modeling tool ProMoT / Simulation tool Diana 27/27

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

References

Schenkendorf, R., Kremling, A., and Mangold, M. (2009). Optimal experimental design with the sigma point method. IET Systems Biology, 3(1):10–23. Zeyer, K. P., Mangold, M., Obertopp, T., and Gilles, E. D. (1999). The iron(III)-catalyzed oxidation of ethanol by hydrogen peroxide: a thermoki- netic oscillator. Journal of Physical Chemistry, 103A(28):5515–5522.

OvGU, MPI 1/1