Efficiently compute the sensitivity of ignition delay and which algorithm is used for computing the local composition sensitivity in Cantera?

[Weiqi Ji]

Hi all,

I guess everyone whichever working on combustion chemistry development/experiment has struggled with calculating the sensitivity of ignition delay time (IDT). The usual way is the brute force approach (finite difference). Cantera provides an easy way to perturb the rate constant, thus is superior to Chemkin in performing this task.

The major problem is that the brute force approach is really time-consuming.  You have to compute the IDT thousands of times for a mechanism with thousands of reactions. And it is only for one condition.

Recently, we figured out a way for this job. The sensitivity direction of the IDT is the same as the local temperature sensitivity at the ignition point. [1] The local composition (include temperature) sensitivity to all of the reactions can be computed with one run of the sensitivity calculations, which has been implemented in Chemkin and Cantera.

Therefore, the job is reduced to compute the local sensitivity, see example: http://cantera.org/docs/sphinx/html/cython/examples/reactors_sensitivity1.html. Then the superior of easily perturb the rate constant is unnecessary and it comes to the question that which one is faster in computing the local composition sensitivity during auto-ignition, between Chemkin and Cantera.

The sensitivity is solved via solve the governing equations of the sensitivity matrix. Chemkin uses the Decoupled Direct Method (DDM). For a mechanism with hundreds of reactions, the cost is only 2~3 times of solving the governing equations of the compositions. [2]

For  Cantera, does anyone knows what algorithm it is using?

Best,

Weiqi

Ref: 

[1] Weiqi Ji, Zhuyin Ren, Chung K. Law, Evolution of Sensitivity Directions during Autoignition, Proceedings of the Combustion Instituts, 2018, doi:10.1016/j.proci.2018.07.005 

[2] T. Turányi, A.S. Tomlin, Analysis of Kinetic Reaction Mechanisms, Springer, Berlin, Heidelberg, 2014

[3] Chemkin Theory menu.

[Ray Speth]

Weiqi,

Cantera uses the Sundials CVODES solver for the sensitivity equations, using the “staggered corrector” approach described in Section 2.6.1 of the CVODES User Documentation.

I have also noticed that the solver performance degrades quite significantly when calculating sensitivities with respect to many rate constants, almost to the point where it’s no faster than the brute force approach, but it’s not something I’ve had a chance to dig into. My suspicion is that it comes from relying on CVODES to calculate the matrix of derivatives with respect to the sensitivity parameters, , which I think ends up happening by finite difference, requiring at each timestep a number of additional residual evaluations equal to the number of sensitivity parameters. CVODES provides methods for allowing user-specified routines that do this calculation more efficiently, but Cantera does not currently make use of these.

Regards,
Ray

[CatDog]

Hi Ray,

I saw in CVODES's handbook  ch 2.7, there is adjoint sensitivity analysis which is capable of evaluting ignition delay as a `free end time` problem using almost the same API (solve_adjoint) of 1D flow reactor of cantera. Why there is no such an API for 0D reactors or networks? The ignition delay can be evaluated as `g(T,y,p)=T` with proper variable end time condition such as `Temp(t=T)=T0+200K` or `d^2 Temp(t=T) / dt^2 = 0` or `d c[OH] /dt= 0`.

[Ray Speth]

Hi,

Yes, the adjoint method is certainly well-suited for ignition delay problems, since the question is the sensitivity of a scalar output to many input parameters. However, it's a bit more complicated to apply to time-dependent problems than steady-state problems. Most of this complexity can be handled by CVODES, but the user still needs to supply a lot of partial derivatives. Calculating these in the naive way, i.e. by finite difference, basically wipes out the computational advantage of the adjoint method. I had coded up a sample problem implementing this at one point, but the performance was disappointing.

Longer term, there is a goal of making Cantera compatible with automatic differentiation tools, which would make calculating these matrices of partial derivatives much more efficient. Implementing the adjoint method at that point would be one of the key applications.

Regards,

Ray

[CatDog]

Hi Ray

For time dependent problem, according to eq 2.21 to 2.22 in the CVODES handbook. For ignition delay problem, reverse integration of only one variable is needed. I am not very sure if it handles `free end time` constraint correctly. It is different from what I know. Assuming it is correct, for a very simple `end time` functional `g(T,y,p)=T` and choosing parameter as `p`=reaction multiplier. The reverse eq 2.21 should be simplified and the only derivative needed the Jacobian (df/dy) which is also needed in forward integration of state equation.

Regards,
Di

[CatDog]

Hi Ray,

I agree with you that the ignition delay problem is highly complicated and the CVODES' treatment actually do not include final time constraint. It seems 

I tried to derive a free final time adjoint. According to table 1 of this lecture note on calculus of variation. The derivative of final time $t_f$ with respect to parameter/control variable is an integral identical to fixed end time problem. But the boundary condition (a.k.a transversality condtion) is different. However, I cannot find workable boundary condition of adjoint/costate equation for free final time problem. 

Another approach to derive free end time adjoint is Long's method to convert free final time problem to fixed final time problem. The final time is treated as an additional parameter. But the literature is optimal control oriented so the results are usually heavily mixed with optimal condition.

Mathias Lemke's recent paper describes something different. But I do not understand how they derive the sensitivity of ignition delay time using a functional like a $H^1_0$ norm. It seems its correctness relies on their numerical experiment.

Best regards,

Di

On Thursday, August 30, 2018 at 11:12:58 PM UTC-4, Ray Speth wrote:

[Weiqi Ji]

Hi Di,

I am not familiar with deriving the adjoint equation. But I think the discussion has gone beyond the development of Cantera, i.e., the numerical algorithms and code development. It is now a physical problem ...

In my opinion, the problem with computing the sensitivity of the ignition delay is that the definition of the "ignition delay time" is not an algebra function of the state variables.

Mathias Lemke's recent paper proposed an expression (or you can call it correlations) for the ignition delay time. However, explaining such correlations requires further works. Our paper mentioned in the above post might help explain the correlation, but the logic not very straightforward.

Best,

Weiqi 

[Ray Speth]

Di,

The reverse integration is for a system of the same size as the original problem -- both lambda and mu are vectors the same size as y.

While you do need df/dy in the forward integration as well, the difference is that there you can use an approximation for the sake of the Newton iteration, i.e. re-use the same Jacobian for multiple time steps. In the reverse integration, I think CVODES wants to have the current value of the Jacobian at each step, which means many more Jacobian evaluations than you have in the forward integration. I'm not sure how much error would be introduced by updating this less frequently, or whether that's possible/advisable within CVODES. In addition, you need the derivatives of f with respect to the parameters to do the quadrature in Eq. 2.20 or 2.21.

Regards,

Ray

[CatDog]

Weiji, Ray,

Yes, it is a little out of scope. 

There are several definition of IDT end time constraint AFAIK: 

1. 0th order condition: T(t_final) = T_0 + 200K. This will result in a standard `free end time, fixed end value` adjoint problem.

2. 1st order condition: d C[OH](t_final) / dt = 0. maximum of C[OH]. It can be reduced to an explicit algebraic expression using the evolution equation of state vector. 

3. 2nd order condition:  d^2 T/ dt^2 = 0, inflection point of T as an indicator of thermal runaway. It can be reduced to an explicit algebraic expression too if you want.

Both 1st-order and 2nd order condition can result in an end time constraint curve in state space, which can defines an end time boundary condition for adjoint problem. I agree with Ray it is very difficult to compute the derivatives without automatic differentiation. However, I just found none of the resulted adjoint problems fit into the CVODES' adjoint formulation which do not include end time constraint's expression.

Mathias Lemke's recent paper 's expression used a slightly shifted solution as reference profile. Its effect is numerical differentiation, that's why I think it is a semi-norm in Sobolev space. I just do not know why it relates to sensitivity so well after scaling.

The adjoint solver implemented in 1D flame in cantera v2.4 is not using CVODES' backward integration, right? I see get_flame_speed_reaction_sensitivities() invokes solve_adjoint in a pyx file, which invokes solveAdjoint in a cpp file, which used evalSSJacobian to update Jacobian.

Regards.

Di.

[CatDog]

Ray,

Sorry, I made a mistake. I thought the backward integration is just one dimensional according to CVODES document. It is not well documented for dimensions of equations and conditions for vector state function. And the Meyer term and Lagrange term in the functional used the same symbol `g`, which is very confusing.

The costate/adjoint equation of IDT for different types of final time constraint is the same. The final boundary condition is the only difference. Jacobian evaluation is always the performance bottle neck if no AD is used.

Evaluation of boundary condition is difficult for inflection point of temperature because third order derivative evaluation is needed at least at final end time.

Best regards,

Di

[CatDog]

Weiqi,

Finally I found something fits IDT adjoint problem with final time constraint.

reference: case IV of free final time problem (a.k.a problem 8 in his table 5-1) in section 5.1 of Kirk's book: optimal control theory: an introduction.

 

Mathias Lemke's approach is fantastic because it represent a point of inverse function/solution in functional form. Usually, Riesz representation theorem of Hilbert space only guarantees that linear bounded functional can be represented as inner product form using elements in the same Hilbert space. It is really difficult to understand how an inverse function is expressed in that semi-norm as a functional.

Best regards,

Di