Solution Multiplicity of Inversion Problems in Distributed Systems

1. Solution Multiplicity of Inversion Problems in Distributed Systems AIChE 2006 Annual Meeting San Francisco, CA November 12 – 17, 2006 Session 10D04: Applied Mathematics in Bioengineering Paper 541e Kedar Kulkarni, Jeonghwa Moon, Libin Zhang Andreas A. Linninger 11/17/2006 Laboratory for Product and Process Design, Department of Chemical Engineering, University of Illinois, Chicago, IL 60607, U.S.A.

2. Motivation: Why Inversion? Data available: Biological Reactor: (at specific locations) Characteristics:- Partial information available- The system cannot be studied as a whole To determine: Reaction: • Kinetic constants: • Rate of the reaction (k1) • Other constants • Thermal constants: • Thermal conductivity (k) • Porosity (ε) • Diffusivity (D) Microbial growth: Heat transfer: What are the values of k1, μmax, Ks and S0 ?

3. Why Multiplicity of Solutions? • Tang et al (2005) suggest that multiplicity in inversion solution is due to errors in experimental data caused by instruments and experimental procedures- Local optimizers get stuck at the first available minima- Using different initial guesses randomly is NOT an efficient way to explore all the optima • Contour plot of the residual error in the space of the kinetic parameters (reaction constants) k1 and k11 for the H2/O2 combustion reactions- Used 9 concentration data-points to solve for these parameters using the Physically-bounded Gauss-Newton (PGN) method- Shows physically meaningful multiple minima

4. Outline: • Review of methodologies used • Global terrain methods (Lucia and Feng, 2002) for Transport and kinetic inversion problems • Evolutionary methods (Niche genetic algorithms) • Case studies • Multiplicity of solutions in catalytic pellet reactor • Multiplicity of inversion solution in a simple bio-reactor • Discovery of metabolic properties in the human brain • Conclusions and Future work

5. Global Terrain Method • Basic concept of Global Terrain Method (Lucia and Feng,2002) • A method to find all physically meaningful solutions and singular points for a given (non) linear system of equations (F=0) • Based on intelligent movement along the valleys and ridges of the least-squares function of the system (FTF) • The task : tracing out lines that ‘connect’ the stationary points of FTF. • Mathematical background • Valleys and ridges in the terrain of FTF could be represented as the solutions (V) to: V = opt gTg such that FTF = L, for all L єL F: a vector function, g = 2JTF, J: Jacobian matrix, L: the level-set of all contours

6. Global Terrain Method • Applying KKT conditions to the this optimization problem we get the following Eigen value problem Hi : The Hessian for the i th function • Thus solutions or stationary points are obtained as solutions to an eigen-value problem where the Eigen values are identical to the KKT multipliers • Initial movement • It can be calculated from M or H using Lanzcos or some other eigenvalue-eigenvector technique (Sridhar and Lucia, 2001) • Direction • Downhill: Eigendirection of negative Eigenvalue • Uphill: Eigendirection of positive Eigenvalue

7. Equations Feasible region Starting point (1.1, 2.0) Global Terrain Method (example) 3D space of case 1

8. Use of Global Terrain Methods for inversion Mathematical Problem Formulation for Inversion Given: Temperature/Concentration Field (from experimental data e.g. MRI) Find: s.t.

9. Use of Global Terrain Methods for inversion (formulation) The set of equations to be solved: - i th model predicted state variable • i th experimental datum for the state variable The Jacobian and the Hessians:

10. Use of Finite Volume Method (FVM) to obtain Jacobian and the Hessians Consider the 2D transport problem with diffusion and reaction Transport Model Integrating over a control volume N n e w P E W s S Discretization Model

11. Formulation to obtain Jacobian and the Hessians After discretization, steady-state transport problem can be expressed as follows: Implicitly differentiating this set of equations, we find Thus, First term Second term

12. Formulation to obtain Jacobian and the Hessians Solve the following linear algebraic system to obtain the jacobians: Implicitly differentiating this set of equations again, we obtain: (A) (B) The linear algebraic system (A) could be solved for the diagonal elements of the Hessians and (B) could be used to obtain the off-diagonal elements

13. Evolutionary methods Simple Genetic algorithms (GA) • GA is an optimization and search technique based on the principals of genetics and natural selection. • It works on a population of possible solutions, while other heuristic methods use a single solution in their iterations. • GAs are probabilistic (stochastic), not deterministic. • They seek the global optimum Define cost function, cost, variables Select GA parameters Generate initial population Find cost for each chromosome Select mates Mating Mutation Convergence Check Solution Flow chart of a continuous GA

14. True error surface Starting point  p1 p2 Use of reduced space hybrid evolutionary methods (rs-HEM) for inversion Gradient based methods rs - HEM Evolutionary methods Reduced space gradient methods + Evolutionary methods Advantages:- Globally optimal - Can handle non- differentiabilities- Can be parallelized- Can handle large scale problemsDisadvantages:- Slow- Cannot handle constraints Advantages:- FastDisadvantages:- Locally optimal- Sensitivity to initial guess- Cannot handle non-differentiabilities Advantages:- Fast- Reduced space approach- Globally optimal - Can handle non- differentiabilities- Can handle constraints- Can be parallelized- Can handle large scale problems

15. Evolutionary mehods Generate initial guess for all variables, θ Generate Lagrangian and apply optimality conditions Generate initial population of the model parameters Define cost function, cost, variables Select GA parameters Generate search direction (Δθ) using gradient based NLP methods Solve for the state variables using the Finite volume method - Gradient based solvers (Newton’s method) Generate initial population Generate correct stepsize by line search Find cost for each chromosome Find least squares error (LSQ) Select mates Select the parameters with least LSQ Update: θnew = θnew +αΔθ Mating Mating, Mutation No Mutation Convergence? No Convergence? Convergence Check Yes Yes Output the optimal parameter set Output the optimal parameter set Use of reduced space hybrid evolutionary methods (rs-HEM) for inversion Gradient based methods rs - HEM Output the optimal parameter set

16. Niche Evolutionary Methods • Genetic algorithms for problem with multiple extrema (multi-modal) • Use the concept of fitness sharing • Restricts the number of individuals within a given niche by “sharing” their fitness, so as to allocate individuals to niches in proportion to the niche fitness niche radius No fitness sharing Fitness sharing

17. Niche Evolutionary Methods (example) • Spots 4 maxima and 1 saddle point

18. Case study I Multiplicity of Solutions in Catalytic Pellet Reactor

19. Cooling Outlet Multiscale Model B A Tubular Reactor Cooling inlet Packed Catalytical Pellet Bed Catalyst Pellet Micro Pores of Catalyst Catalytic Pellet Reactor s.t. Darcy’s law Mass and energy balance Pellet model Optimal Parameters: D= 7.32  10-5 (m2/s) ; k =0.82 10-4 (1/s); E=5268.62 (J/mol)

20. Multiplicity of inversion solution Reaction constant, k Diffusivity, D Contours of the least squares error function in the parameter space • Different initial guesses yield different minima • Multiple solutions obtained when the experimental data has error • The solution obtained will be the mimimum that is nearest to the initial guess

21. Case study II Simple Bio-reactor

22. The transport problem Species B Species A BC(A)2 BC(B)2 BC(A)3 BC(A)4 BC(B)4 BC(B)3 BC(A)1 BC(B)1 BC(B)1: BC(A)1: BC(B)2: BC(A)2: BC(B)3: BC(A)3: BC(B)4: BC(A)4:

23. Multiplicity of inversion solution Concentration profile along the length of the reactor • S1 is a minimum • S2 is a saddle point Simultaneous diffusion and reaction S1 x S2 x

24. Case study III Discovery of unknown transport/metabolic properties in the human brain

25. Discovery of unknown transport and metabolic/kinetic properties Problem formulation: s. t. Convection-Diffusion equation Transport of the bulk fluid Concentration field obtained from imaging data Model predicted concentration field of the drug Covariance matrix of the measurements To find:

26. NH 2 | CH – C – CO H OH 2 2 | CH 3 OH Discovering metabolic properties for F-Dopa Brain Tissue Blood Circulation System [18F] Methyl F-DOPA [18F] Methyl F-DOPA Methylation Methylation INJECTION [18F]FDA [18F]FDOPA [18F]FDOPA L-Dopa Methyl Dopa L-Dopa Dopamine To determine: Blood-tissue clearance K1 Tissue-blood clearance k2 The diffusivities DD , DM s. t. F-dopa clearance Methyl F- dopa clearance

27. 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 Clinical concentration field of F-dopa Unstructured computational grid s. t. Fdopa clearance Methyl - Fdopa clearance Discovering metabolic properties for F-Dopa (cont’d) Putamen Predicted F-Dopa concentration field with optimal parameters ( Gjedde, A. et. al. Neurobiology, 88, 2721-2725, 1991 ) Best-fit Parameters:

28. Conclusions • Solution multiplicity in distributed systems can be efficiently handled using global terrain methods and niche evolutionary methods • Experimental error does not necessarily lead to multiplicity • Future Work • Consideration of dynamic systems • Higher number of parameters

29. Thank you!