Exact two layer solution

[hdrake]

I tried to implement the exact solution to the two-layer problem (sum of a fast-timescale exponential and a slow-timescale exponential) rather than the one-exponential deep-layer model we have now.

For some reason, introducing a second timescale makes the problem unsolvable by Ipopt via JuMP. The algorithm solves successfully for time horizons less than about 2150, but then I get errors of "search_direction_too_small" or "iterations_diverging". I'm going to try out NLopt and see if it's any better.

[hdrake]

NLopt unfortunately no longer supported in JuMP because they didn't switch to MathOptInterface, which is required as of JuMP v0.19.0

It may be possible to speed up MARGO (and deal with divergence issues) by being smart about warm starts, e.g. by setting the variables

set_start_value(M[1], 0.0)

[hdrake]

Another unfortunate problem with the exact two-layer solution. Forward finite differences results in huge errors if the timestep is on the order of– or longer– than one-third of the fast timescale (~4 years for CMIP5 values). In practice, I found I need to keep the timestep below about 0.5 years to achieve good results. This could be improved somewhat with a better time-stepping scheme but with the 10–20 year timesteps we want for the web-app, I'm not sure we will be able to get around this.

We might be able to write discrete approximate solutions to the integral kernels... Not sure if it would be worth it.