Potential Problems in Convergence Criteria for Steady State

[jqbond]

Dear @ChrisRackauckas et. al.,

I am following up on a potential issue I ran into with convergence criteria on the DynamicSS solver and/or TerminateSteadyState callback. I spent a week or so trying to track down the quirk, and my tentative conclusion is that, for my particular system, I was running into situations where the TerminateSteadyState callback would terminate integration prematurely.

The issue I have is probably a bit of a niche one, but it may come up from time to time, so I wrote up a benchmarking problem to illustrate what I think is happening. The high level version is that I'm working on a parameter estimation problem that requires me to solve for the steady state value of multiple state variables. In general, it is hard-to-impossible to provide good guesses as to the true solution, so I opt to solve the temporal problem and integrate to steady state. The solution then gets compared with laboratory measurements, and we use a least squares routine to regress parameter values.

The differential equations that describe this system tend to be an extremely stiff system with state variables and derivative values spanning many orders of magnitude. I believe what is happening with DynamicSS() is that, at times, the system will meet the steady state convergence criteria without truly being at steady state. I spoke with Chris about this, and, in theory, with good values of relative tolerance specific to each state variable, you could probably handle this using the default test criteria in TerminateSteadyState. However, in my case, the magnitude of the state variable and its derivative are highly sensitive to the parameter values, and the parameter values are unknown a priori, which makes it difficult to provide rational values for relative tolerances at the outset.

I created what is probably a fairly contrived problem (I stripped out the parameter estimation part and reduced the number of state variables to two), but I think it highlights the issue reasonably well. The two state variables approach steady state values that are 8 orders of magnitude apart on time scales that are probably 15 orders of magnitude apart.

I encounter this type of physical system in describing the accumulation of species adsorbed on surfaces, where the rates and equilibrium values can vary by many orders of magnitude. Superficially, one might conclude that a steady state value of 1e-8 or 1e-23 are both still effectively zero compared to the other state variable with a steady state value very near 1; however, in calculating rates, these values are often multiplied by proportionality constants that may be on the scale of 1e10 or 1e20, so even small numbers are important in calculating the physically observable response of the system, i.e., this makes a big difference in the parameter estimation problem.

In any event, I'm linking the benchmarking code below. I'm still not 100% sure these issues aren't operator error, but I've done my level best to suss out the details, and I provide a set of parameters where integrating the problem to t --> infinity gives a different set of values for the state variables than using the DynamicSS solver.

https://github.com/jqbond/Test-Cases/blob/master/SSPROBLEM.jl

I played around with tracking the mean value and standard deviation of the state variables as an alternate criteria for termination. It works fine, but, as Chris pointed out to me, is a bit more expensive. I have not played around much with alternative criteria, but am happy to do so if others can confirm this as a potential issue with TerminateSteadyState.

Sincerely,

Jesse

[ChrisRackauckas]

We'll keep this in mind. That tracking of the running mean could be a new callback and algorithm choice.

[jqbond]

Thanks @ChrisRackauckas! The steady state solvers are a wonderful feature in Julia, so I'll look forward to seeing more about them in the future. Please let me know if you need me to kick the tires at any point.

[ChrisRackauckas]

Any contribution would be great since I don't see myself having the time for this any time soon, but wanted to at least have this around to consider having a future summer student look at psuedo-transient methods and more robust measures of steady state handing.

[jqbond]

Sounds good - I'll probably continue playing with it since I have to deal with this problem in research so often. If I find a termination condition that I like and is not slower than just integrating to t --> inf, I'll update that test problem with a different callback and ping you. I will leave making it efficient to the professionals.

Good luck starting your semester!

[ChrisRackauckas]

Did you ever make a code for this?

[jqbond]

@ChrisRackauckas let me review it again and get back to you in the next few weeks. I don't think I did anything with it after initially creating this test case. I ended up chasing down a few different issues, and I don't think I ever went back to the steady state solver.

[ChrisRackauckas]

Solved in https://github.com/SciML/SteadyStateDiffEq.jl/pull/50