ClimateMARGO / ClimateMARGO.jl

Julia implementation of MARGO, an idealized climate-economic modelling framework for Optimizing trade-offs between emissions Mitigation, Adaptation, carbon dioxide Removal, and solar Geoengineering.
https://margo.plutojl.org/
MIT License
67 stars 13 forks source link

Coherent & concise JuMPification of MARGO? #18

Open hdrake opened 4 years ago

hdrake commented 4 years ago

As mentioned in https://github.com/hdrake/ClimateMARGO.jl/issues/14, it would be nice to reduce the JuMP implementation of MARGO to a few single-line statements like:

for i in N; @constraint{ T(M,R,G,A) < Tstar }; end
@NLobjective( model_optimizer, Min, net_present_cost(M,R,G,A) );

which reuse the diagnostic functions such as T which fully describe the model (see https://github.com/hdrake/ClimateMARGO.jl/tree/updating-structure) but are complicated functions of both the four controls {M,R,G,A} as well as a number of prescribed inputs functions and parameters. In the current implementation https://github.com/hdrake/ClimateMARGO.jl/commit/af16f490a3fa55754cca4738bead7c1877a68064, these optimization constraints are hard-coded from scratch, span many lines, and are repeated several times over.

hdrake commented 4 years ago

The first thing to do is set up some very simple tests of JuMP and see how far we can push the use of functions to simplify the optimization before we break JuMP.

hdrake commented 4 years ago

Looks like some of this is possible with the user-defined functions API in JuMP. I'm a bit nervous about whether this will extend to the cumsum implementation of the Green's functions kernels in the energy balance model, since it's not just a combination of elementary functions...

fonsp commented 3 years ago

Hey Henri! With only a small modifcation, we can use diagnostic functions like control cost, damages, and final temperature inside JuMP. These are Vector -> scalar functions, and there’s a trick for that: https://jump.dev/JuMP.jl/v0.21.1/nlp/#User-defined-functions-with-vector-inputs-1

Using this trick, we can minimize control costs subject to temp[2200] <= T_max (i.e. overshoot allowed), all using MARGO's diagnostic functions.

But we won’t be able to use Vector -> Vector functions, which is needed for the global temperature constraint. You need to write out the cumsum behaviour like we currently do, but I think it can still be made more concise by partially reusing the diagnostic code.

I wrote a more detailed notebook about it here, which is just the mitigation optimization. To run the notebook, you need to checkout the forward-diffable branch of ClimateMARGO, so I also made a PDF.

Have a look at "Conclusions" on page 6 🎈 jump test margo.jl ⚡ Pluto.jl ⚡.pdf

https://github.com/fonsp/disorganised-mess/blob/master/jump%20test%20margo.jl

hdrake commented 3 years ago

Awesome progress @fonsp! I'll take a look at it in the next couple of days!

hdrake commented 3 years ago

I feel like there should be some trick to generalize the clunky cumsum formulation we have now, even if it has to be done separately for individual Vector -> Vector operators...

fonsp commented 3 years ago

To add to this: JuMP uses forward mode automatic differentiation (dual numbers). This means that it can automatically differentiate our diagnostic functions 🎉, and also user-defined functions 🎊, even though they are a 'black box' from the perspective of ClimateMARGO.jl. As long as they are not 'too complex' (not sure what that means yet), and as long as they accept the more general Real type instead of the Float64 type, it will work.

If we ever have a user whose custom functions are not auto-diffable, JuMP does not include finite differencing by default, but it seems like we can add that functionality later, either by supplying a finite difference method as the 'derivative', or by switching packages.

fonsp commented 3 years ago

In the end, I was able to also use Vector->Vector functions with an additional wrapper function, but the performance impact is very large (~ 10x - 500x slower). I tried lots of things, but I can't figure out exactly why it is this much slower...

If anyone is interested, I can document what I tried. Code is here, here and here.