Solution of Extreme Transcendental Differential Equations

Stuart Nettleton

University of Technology Sydney
August 2013

Abstract

Extreme transcendental differential equations are found in many applications including geophysical climate change models. Solution of these systems in continuous time has only been feasible with the recent development of Chebyshev solvers such as Mathematica 9’s NDSolve function. This paper presents the challenges and means of solving the widely used DICE 2007 integrated assessment model in continuous time. Application of the solution technique in a mobile policy tool is discussed.

An example

Transcendental Differential Equation System

equations={k'[t]==-0.1` k[t]-269.3998529790092` (E^(-9.2` (E^(-0.001` t)))) (((8600-2086 (E^(-0.035` t))))^(0.7`)) If[Tat[t]≤0,1,(1/(1+0 Tat[t]+0.0028388` ((Tat[t])^2)))] ((k[t])^(0.3`)) (-1+0.0024598146458005322` (E^(2.433333333333333` (E^(-0.003` t)))) (1+(E^(-0.005` t))) ((μ[t])^(2.8`)))-c[t],Tat'[t]==0.22` (-36.7073066928394`+5.482241155378061` Log[Mat[t]]+(({{\[Piecewise], ({{-0.06`+0.0036` t, t≤100}, {0.3`, True}})}})))-0.06441748633333333` Tat[t]+0.011002153` Tlo[t],Tlo'[t]==0.0048` Tat[t]-0.0048` Tlo[t],Mat'[t]==1.1` (E^(-0.01` t))-0.0190837` Mat[t]+0.009800871837862596` Mup[t]-3.171771574386283` (E^(2.433333333333333` (E^(-0.003` t))-9.2` (E^(-0.001` t)))) (((8600-2086 (E^(-0.035` t))))^(0.7`)) ((k[t])^(0.3`)) (-1.`+μ[t]),Mup'[t]==0.0190837` Mat[t]+0.0003369934177753544` Mlo[t]-0.015203871837862596` Mup[t],Mlo'[t]==-0.0003369934177753544` Mlo[t]+0.005403` Mup[t],
μ'[t]==(1/10) μconst (1-μ[t]) μ[t]};
NDSolve[Flatten[{equations/.{c[t]->38.1,μconst->1.2},k[0]==136.7,Tat[0]==0.7307,Tlo[0]==0.0068,Mat[0]==808.9,Mup[0]==1255.0,Mlo[0]==18365.0,μ[0]==0.005}],{k'[t],k[t],Tat[t],Tlo[t],Mat[t],Mup[t],Mlo[t],μ[t]},{t,0,1400}]

Plot[Evaluate[%[[1,All,2]]],{t,0,1400}]

Transcendental Differential Equations

Transcendental equations contain transcendental or non-algebraic functions such as (e^x), Log[x] and Cos[x]. The defining characteristic of transcendental functions is that the roots are not algebraically independent, which means that the roots cannot be expressed as the solution to a polynomial equation whose coefficients are polynomials with rational coefficients.

DICE 2007 example system

It is often difficult to appreciate the sensitivity of complex social systems to changing constraints without simulating key interrelationships in the system. The DICE 2007 integrated assessment model was developed to understand the interrelationships between climate change, the social cost of carbon and efficient carbon abatement trajectories [1], [2], [3]. It has become a classic climate change policy simulation tool for evaluating the social and geophysical effects of global warming. An indicator of the success of models such as DICE 2007 in policy formation is the embedding of results and recommendations within the national climate change policies of many countries.

The DICE 2007 model is an optimization problem defined by a system of transcendental differential equations. Figure 1 summarizes a continuous 2007 DICE formulation utilizing parameters fitted through dynamic programming [4], [5].

Maximize
utility per capita=381000 + (1/194)(∫_0^(1400years))(e^(-0.015 t))(Subscript[c, t])(1-((Subscript[l, t])/(Subscript[c, t]))) dt
Where:
(Subscript[c, t])≥0 is global consumption at time t
(Subscript[l, t])=6514(e^(-0.035t))+8600(1-(e^(-0.035t))) is global population at time t
(Subscript[k, t])≥0 is global capital stock at time t
(Subscript[Tat, t]), (Subscript[Tlo, t]) is global mean surface and lower ocean temperature rise at time t
(Subscript[Mat, t]), (Subscript[Mlo, t]), (Subscript[Mlo, t]) are the masses of carbon in the atmosphere, upper & lower oceans at time t
0≤(Subscript[μ, t])≤1 is amelioration & abatement proportion at time t

Subject to:
(a) The equations of motion:
(Subscript[k, t]^('))= -(Subscript[c, t])-0.1(Subscript[k, t]^('))-(269.4(e^(-9.2(e^(-0.001t))))(Subscript[k, t]^0.3)(Subscript[l, t]^0.7)(0.00245981(e^(2.43333(e^(-0.003t))))(1+(e^(-0.005t)))(Subscript[μ, t]^2.8)-1))/If[(Subscript[Tat, t])≤0, 1, 1+0.0028388(Subscript[Tat, t]^2)]
(Subscript[Tat, t]^('))=0.22(If[t≤100,0.0036t-0.06,0.03]+5.48224Log[(Subscript[Mat, t])]-36.7073)-0.0644175(Subscript[Tat, t])+0.0110022(Subscript[Tlo, t])
(Subscript[Tlo, t]^('))=0.0048((Subscript[Tat, t])-(Subscript[Tlo, t]))
(Subscript[Mat, t]^('))=-3.17177(e^(2.43333(e^(-0.003t))-9.2(e^(-0.001t))))(Subscript[k, t]^0.3)(Subscript[l, t]^0.7)((Subscript[μ, t])-1)-0.0190837(Subscript[Mat, t])+0.00980087(Subscript[Mup, t])+1.1(e^(-0.01t))
(Subscript[Mup, t]^('))=0.0190837(Subscript[Mat, t])+0.000336993(Subscript[Mlo, t])-0.0152039(Subscript[Mup, t])
(Subscript[Mlo, t]^('))=0.005403(Subscript[Mup, t])-0.000336993(Subscript[Mlo, t])
(b) Atmospheric temperature rise constraint 
(Subscript[Tat, t])≤a maximum temperature rise level for all t (say 2°, which is the internationally agreed "safe level") 

Figure 1: Continuous Formulation of the DICE 2007 Model

Materially adding to the challenge of extreme transcendentality is the outer intertemporal joint optimization across the two unknown functions (Subscript[c, t]) and (Subscript[μ, t]) together with the continuous constraint (Subscript[Tat, t])≤ (Subscript[T, max]). Such extreme transcendentality, non-linearity and the presence of a large and possibly infinite number of local maxima between the two unknown functions to be jointly optimized has led to this apparently concise system of equations becoming renowned as extraordinarily difficult to solve.

DICE 2007 Solution

Traditional and enhanced-traditional approaches

The DICE 2007 model has always been solved by discretizing into decade time periods and introducing a helper Savings Ratio constraint to eliminate consumption (Subscript[c, t]) as an unknown function. Usually this is calculated in GAMS with the CONOPT solver [1], [2].

Recently, decade time periods have been reduced to discrete single year periods using Dynamic Programming, Approximate Dynamic Programming and Optimal Control Theory [4], [6], [7], [8].

New approach: Chebyshev solver

Mathematica 9’s NDSolve numerical differential algebraic equation function is able to solve the DICE 2007 model using in-built Chebyshev polynomial approximation [9], [10]. Furthermore, NDSolve solution of transcendental differential equations is quick and suitable for placing within the global optimizer NMaximize.

Resolving the two completely unknown equations (Subscript[c, t]) and (Subscript[μ, t])

Notwithstanding the extraordinarily powerful combination of NMaximize and NDSolve, the overwhelming interaction between two the simultaneously unknown functions (Subscript[c, t]) and (Subscript[μ, t]) is a major barrier to rapid solution.

Rather than following the traditional approach of making an arbitrary assumption to eliminate (Subscript[c, t]), it is more realistic to consider a likely form of the unknown function (Subscript[μ, t]).

Amelioration and abatement function (Subscript[μ, t]) as generic sigmoid

Policy makers may choose any number of forms for (Subscript[μ, t]), including "business as usual" leaving the form of (Subscript[μ, t]) totally to market forces (in the sense that as discomforting levels of damage begin to occur then consumers, insurers, producers and local authorities will implement amelioration, abatement and adaption measure to reduce the exposure to damage).

From the perspective of evaluating future choices, policy makers readily understand notions such as Pareto’s 80/20 rule, where 20% of effort brings 80% of the effect, while achieving the remaining 20% of effect takes 80% of the effort. This reflects the idea of low-hanging fruit, where some amelioration and abatement efforts will be low cost and readily achieved while others measure will become progressively more painful and expensive. The generic sigmoid is a well understood mathematical expression of this concept.

In the context of (Subscript[μ, t]), a generic sigmoid allows amelioration and abatement to rise very slowly through the current level (0.5% assumed by DICE 2007) to saturate at 100% after a number of decades. Such a sigmoid may be described by the analytical equation of motion: (Subscript[μ, t])' = const (1-(Subscript[μ, t])) (Subscript[μ, t]) . This equation of motion greatly simplifies the DICE 2007 optimization problem by reducing the unknown (Subscript[μ, t]) function to just one unknown scalar quantity, which is the acceleration constant that determines whether the sigmoid rises quickly or over an extended time period.

Consumption function (Subscript[c, t]) as Chebyshev polynomial

Resolving the form of (Subscript[μ, t]), albeit still an unknown function, allows the other unknown function (Subscript[c, t]) to be determined as a Chebyshev polynomial defined though a eleven  judiciously chosen handles (e.g.  (Subscript[c, 25]), (Subscript[c, 50]), (Subscript[c, 75]), (Subscript[c, 100]), (Subscript[c, 150]), (Subscript[c, 200]), (Subscript[c, 350]), (Subscript[c, 550]), (Subscript[c, 800]), (Subscript[c, 1100]), (Subscript[c, 1400])). Including with the single (Subscript[μ, t]) acceleration constant there are just twelve optimization parameters.

Memoized Optimizing Constraint Solver

Indicative code for the Memoized Optimizing Constraint Solver is provided below. The solver ensures that all constraints are correctly calculated and returned to the global optimizer in the required mode.

An important feature of the Memoized Optimizing Constraint Solver is that at each iteration the objective equation calls the Chebyshev solver to solve the System of Transcendental Differential Equations. This solution is retained within the solver’s DownValues array, where it is subsequently drawn upon by each of the constraints.  There is no limit to the number of constraints except for the need to carefully prepare the logic. In order to reduce the use of memory, the solver’s DownValues array is periodically reset.

Constraints to continuous functions provide challenges in optimization. The Memoized Optimizing Constraint Solver addresses this with in-built utilities to find maximum and minimum values across continuous functions as well as a value at a particular time period.

The latter is quite important in conditioning the optimization for a stable solution that avoids a pervasive issue that often plagues dynamic models, which is the degenerate solution of maximizing consumption through cannibalizing capital. In the continuous formulation it is possible to directly address the cause by simply requiring that capital is greater than zero at the end of the simulation. This avoids the need for an interim stabilizing constraint such DICE’s savings ratio or a DuPont Ratio-like production to capital ratio that can become less appropriate in the presence of extended improvements to Total Factor Productivity, for example with simulations over many centuries [23, 24] .

statevars={k[t],Tat[t],Tlo[t],Mat[t],Mup[t],Mlo[t],μ[t]};
solver[params_]:=Check[NDSolve[Flatten[{equations/.Thread[parsarray->params],Thread[(statevars/.t->0)==Last[xvec3][[2]]]}],Join[{k'[t]},statevars],{t,0,testperiods}, MaxSteps->(10^8),InterpolationOrder->All]//Quiet,Print["**** solver failed ****"];{}];
cons0[params_]:=cons0[params]=Module[{optsolver0,optsolver,findmax,findmin,conssol,findopt,finddiff,findval},
optsolver0=solver[params];
optsolver=If[optsolver0=={},Return[{-1,-1,1}(10^30)],optsolver0[[1]]];
findmax[var_,imtestyear_]:=Check[If[MatchQ[var/.optsolver/.{t->imtestyear},_Complex],
 Abs[Im[var/.optsolver/.t->imtestyear]](10^3),FindMaxValue[var/.optsolver,{t,imtestyear,30,testperiods}]]//Quiet,(10^30)];
findmin[var_,imtestyear_]:=Check[If[MatchQ[var/.optsolver/.{t->imtestyear},_Complex],
 -Abs[Im[var/.optsolver/.t->imtestyear]](10^3),FindMinValue[var/.optsolver,{t,imtestyear,40,testperiods}]]//Quiet,-(10^30)];
finddiff[var_,imtestyear_]:=Check[If[MatchQ[var/.optsolver/.t->imtestyear-var/.optsolver/.(t->imtestyear-1),_Complex],
 -Abs[Im[var/.optsolver/.t->imtestyear-var/.optsolver/.(t->imtestyear-1)]],
var/.optsolver/.t->imtestyear-var/.optsolver/.(t->imtestyear-1)]//Quiet,-(10^30)];
findval[var_,imtestyear_]:=Check[If[MatchQ[var/.optsolver/.t->imtestyear,_Complex],
 -Abs[Im[var/.optsolver/.t->imtestyear]],var/.optsolver/.t->imtestyear]//Quiet,-(10^30)];
findopt=((Subscript[cumutil, 0])/.sbsd)+NIntegrate[(Exp[-ρ t](((c[t]/(Subscript[l, t]))^(1-γ))-1)(Subscript[l, t])/((Subscript[scale, 1])(1-γ)))/.sbsd/.Thread[parsarray->params]/.optsolver,{t,0,testperiods},WorkingPrecision->24];
conssol={findopt,findval[k[t],testperiods],findmax[Tat[t],80]}
]//Quiet
cons1[params:{_?NumberQ..}]:=cons0[params][[1]]
cons2[params:{_?NumberQ..}]:=cons0[params][[2]]
cons3[params:{_?NumberQ..}]:=cons0[params][[3]]
starttime=AbsoluteTime[];
stepresult={};stepresult>>"teststep";
Print["Starting Memory Used: "<>ToString[N[MaxMemoryUsed[]/10^6]]<>" Mb"];
sol=Block[{ii=1,jj=1,kk=1,temp=""},
NMaximize[Join[{cons1[optpars],0≤cons2[optpars],0.995 tatmax≤cons3[optpars],cons3[optpars]≤1.005 tatmax},Reverse[Thread[Rest[Most[optpars]]>Most[Most[optpars]]]],{copt1>0,10≥μopt,μopt>0}],Map[{#[[1]],Min[#[[2]],#[[3]]],Max[#[[2]],#[[3]]]}&,Transpose[{optpars,pars*0.995,1.005pars}]],MaxIterations->3500,WorkingPrecision->24,
StepMonitor:>(stepresult={ii,jj,kk,Thread[{Join[{res},Map[ToExpression["c"<>ToString[#]]&,Range[Length[handels]]],{μ}],
Join[{cons1[optpars]},optpars]}]};
stepresult>>>"teststep";jj++;kk=1),EvaluationMonitor:>(ii++;kk++;
If[Mod[ii,100]==0,(*If[Mod[ii,1000]==0,evalresult={ii,jj,kk,Thread[{vars,varsopt}]};
evalresult>>>"m138pNum02eval"];*)DownValues[cons0]=Last[DownValues[cons0]];
NotebookDelete[temp];
temp=PrintTemporary["Step: "<>ToString[jj]<>"; Sub "<>ToString[kk]<>"; Cum "<>ToString[ii]<>"; Curr Mem: "<>ToString[N[MemoryInUse[]/10^6]]<>"; Peak: "<>ToString[N[MaxMemoryUsed[]/10^6]]<>" Mb"<>"; Hours: "<>ToString[N[(AbsoluteTime[]-starttime)/3600]]]])]]//Timing

Results

(Legended[, (LineLegend["Tat 2°C", "Tat no constraint", "Tat 3°C", "Tat 4°C", "Tat 5°C", "Tat 6°C", "μ 2°C", "μ no constraint", "μ 3°C", "μ 4°C", "μ 5°C", "μ 6°C"])])

(Legended[, (PointLegend["Consumption c/100", "Capital k/100", "Abatement Proportion μ", "Temperature Tat"])])

Discussion

A continuous solution to the DICE 2007 model is important for main reasons. The first is that the DICE 2007 model is widely used as a discrete decade period model, which is rather inconvenient for policy makers.  For example, the impact of policy over the next two to three decades cannot be properly assessed with just three calculated points where each represents an average across a decade.

More important is the assessment of industry policy in each country using a dynamic global multi-regional model that settles markets in the presence of the competition from other economies that are interconnected through trade. Such dynamic models are usually computable general equilibrium (CGE) formulations that dynamically optimize the evolution of economies from the current economic position in the presence of geophysical constraints, such as DICE 2007 geophysical constraints, and other constraints such as labor and land use.

The starting point, or current global economic position, is generally drawn from  Input-Output data provided by the Global Trade Analysis Project [11], EORA [12], EXIOBASE [13], WIOD [14], GRAM [15] or IDE-JETRO AIIOT [16].

Mobile Application

Policy making involves both modelling and a social research. It is important to put policy modelling directly into the hands of policy makers so they can test their own scenarios fast, in order to speed-up debate across informed viewpoints .

Bringing the extraordinary power of dynamic CGE models directly to policy makers (and other users such as corporate strategic planners) is a challenging task.  Firstly, CGE techniques are inherently complex and there is an enormous cross-disciplinary mathematical, economic and geophysical task in achieving meaningful outputs. Secondly, computational demands are equally challenging and technology has been a major limitation.

All of these factors have changed in the last few of years. It is now feasible to provide policy makers with the immediacy of a mobile application using back end cloud computation. While highly qualified researchers might be generally disintermediated from an in-line role in the policy making process, the experts remain of great importance in an off-line role, engaged in tasks such as geographic centering and the pre-computation of major scenarios so global optimizers are set in the region of final solutions.

The emergence of a new form of CGE for constrained non-linear systems that simultaneously settles thousands of markets in both price and volume using shadow prices (i.e. marginal utilities) has greatly enhancing the ability to develop a mobile application. Although this form of CGE is quite new, it is a continuing theme from John von Neumann [17] to Paul Samuelson [18], Wassily Leontief [19], Michael Farrell [20], Anne Carter [21] and Thijs ten Raa [22].  Nettleton demonstrated this theme in a discrete global multi-regional CGE model [23, 24, 25 & 26].

The major achievement of solving DICE 2007 in continuous time demonstrated in this presentation now provides the ability to develop a continuous global multi-regional CGE model for a mobile application incorporating fast solution techniques and high scaling. It is hoped that the Mathematica Cloud will provide a facility to back end such a mobile application for next phase of this development.

Conclusion

This presentation has demonstrated the solution of extreme transcendental differential equations using a modern Chebyshev solver, Chebyshev function fitting and a highly developed memoized optimizing constraint solver. These techniques have been applied solve the traditional DICE 2007 climate change model in continuous time. The achievement of a continuous solution to the DICE 2007 model is an important step in developing a mobile application with high scaling and fast solution.

Acknowledgement

The author would like to acknowledge the encouragement and support of Daniel Lichtblau, Kernel Technology group at Wolfram Research, Inc. in developing the DICE 2007 solution as well as Daniel’s earlier assistance in memoizing the constraint solver.

References

[1] W. D. Nordhaus, New base DICE-2006 as of November 16, 2006 (D_111606_alt_base.gms). 2006 [Online]. Available:  http://www.econ.yale.edu/~nordhaus/homepage/D_111606_alt_base.GMS. [Accessed: 04-Mar-2013]
[2] W. D. Nordhaus, “Notes on how to run the DICE model,” Oct. 2007 [Online]. Available: http://www.econ.yale.edu/~nordhaus/homepage/DICE2007.htm. [Accessed: 26-Jun-2008]
[3] W. D. Nordhaus, A Question of balance: weighing the options on global warming policies. Yale University Press, 2008.
[4] Y. Cai, K. L. Judd, and T. S. Lontzek, “Open science is necessary,” Nat. Clim. Change, vol. 2, no. 5, pp. 299–299, 2012.
[5] Y. Cai, K. L. Judd, and T. S. Lontzek, CJL -- one year. National Bureau of Economic Research, 2012 [Online]. Available: https://sites.google.com/site/openscienceletter/results/cjl----one-year/gams-code. [Accessed: 07-Mar-2013]
[6] M. Webster, N. Santen, and P. Parpas, “An approximate dynamic programming framework for modeling global climate policy under decision-dependent uncertainty,” Comput. Manag. Sci., pp. 1–24, 2012.
[7] A. Haurie, “Integrated assessment modeling for global climate change: an infinite horizon optimization viewpoint,” Environ. Model. Assess., vol. 8, no. 3, pp. 117–132, 2003.
[8] A. Smirnov, “Attainability Analysis of the DICE Model,” International Institute for Applied Sytems Analysis, Laxenburg, Austria, Interim Report IR-05-049, Sep. 2005 [Online]. Available: http://www.iiasa.ac.at/publication/more_IR-05-049.php. [Accessed: 05-Jul-2013]
[9] Wolfram Inc., “NDSolve—Wolfram Mathematica 9 Documentation,” 2013. [Online]. Available: http://reference.wolfram.com/mathematica/ref/NDSolve.html. [Accessed: 07-Jun-2013]
[10] Adduci, J. “Numerical Optimization in Mathematica “ Wolfram Research Presentation, 27 March 2013
[11] Purdue University Department of Agricultural Resources, “GTAP 8.1 Database,” 2013.  [Online]. Available: https://www.gtap.agecon.purdue.edu/databases/v8/default.asp
[12] M. Lenzen, K. Kanemoto, D. Moran, and A. Geschke, “Mapping the Structure of the World Economy,” Environ. Sci. Technol., vol. 46, no. 15, pp. 8374–8381, Aug. 2012 [Online]. Available: http://dx.doi.org/10.1021/es300171x. [Accessed: 07-Aug-2013]
[13] A. Tukker, “EXIOPOL: towards a global Environmentally Extended Input-Output Table,” presented at the GTAP 2008 Conference, Helsinki, Finland, 2008, vol. Ref 2366 [Online]. Available: https://www.gtap.agecon.purdue.edu/resources/res_display.asp?RecordID=2702. [Accessed: 14-Nov-2008]
[14] E. Dietzenbacher, B. Los, R. Stehrer, M. Timmer, and G. de Vries, “The Construction of World Input–Output Tables in the WIOD Project,” Econ. Syst. Res., vol. 25, no. 1, pp. 71–98, 2013 [Online]. Available: http://www.tandfonline.com/doi/abs/10.1080/09535314.2012.761180. [Accessed: 14-Aug-2013]
[15] M. Bruckner, S. Giljum, C. Lutz, and K. S. Wiebe, “Materials embodied in international trade–Global material extraction and consumption between 1995 and 2005,” Glob. Environ. Change, vol. 22, no. 3, pp. 568–576, 2012 [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0959378012000350. [Accessed: 14-Aug-2013]
[16] B. Meng, Y. Zhang, and S. Inomata, “Compilation and Applications of IDE-JETRO’s International Input–Output Tables,” Econ. Syst. Res., vol. 25, no. 1, pp. 122–142, 2013 [Online]. Available: http://www.tandfonline.com/doi/abs/10.1080/09535314.2012.761597. [Accessed: 14-Aug-2013]
[17]	J. Von Neumann, “A model of general economic equilibrium,” Readings Welf. Econ., vol. 13, no. 1945, pp. 1–9, 1938. 
[18]	S. Dorfman, P. A. Samuelson, and R. M. Solow, Linear Programming and Economics Analysis. New York: McGraw-Hill, 1958. 
[19]	W. W. Leontief, “Input-Output analysis and economic structure: Studies in the structure of the American economy: Theoretical and Empirical Explorations in Input-Output Analysis,” Am. Econ. Rev., vol. 45, no. 4, pp. 626–636, Sep. 1955. 
[20]	M. J. Farrell, “The measurement of productive efficiency,” J. R. Stat. Soc. Ser. A Stat. Soc., vol. 120, no. 3, pp. 253–82, 1957. 
[21] T. ten Raa and W. J. Baumol, “Anne Carter and Input-Output: Technology, Trade, and Pollution,” OEconomia, vol. 2011, no. 01, pp. 61–73, 2011. 
[22]	T. ten Raa, The Economics of Input Output Analysis. New York: Cambridge UniversityPress, 2005 [Online]. Available: www.cambridge.org/9780521841795
[23]	S. J. Nettleton, “Benchmarking climate change strategies under constrained resource usage.,” UTS theses submitted as part of the Australasian Digital Theses Program (ADT), University of Technology, Sydney, Sydney, Australia, 2010 [Online]. Available: http://utsescholarship.lib.uts.edu.au/iresearch/scholarly-works/handle/2100/1012. [Accessed: 30-Mar-2010]
[24]	S. J. Nettleton, “The Service Science of Climate Change Policy Analysis: applying the Spatial Climate Economic Policy Tool for Regional Equilibria,” presented at the International Input Output Association (June 2010), Sydney University, Australia, 2010 [Online]. Available: http://viXra.org/abs/1210.0144
[25]	Wolfram, Inc., Dynamics of Change: Modeling the Economic Effects of Global Warming with Mathematica. Champaign, IL.: Wolfram, Inc., 2011 [Online]. Available: http://www.wolfram.com/mathematica/customer-stories/modeling-the-economic-effects-of-global-warming-with-mathematica.html. [Accessed: 28-Jan-2012]
[26]	S. J. Nettleton, “Computable General Equilibrium (CGE) of Multiregion Input-Output Model,” Wolfram Demonstrations Project, 2011 [Online]. Available: http://demonstrations.wolfram.com/ComputableGeneralEquilibriumCGEOfMultiregionInputOutputModel/. [Accessed: 23-Jan-2012]

Attachment

The DICE 2007 problem requires a careful configuration of the Optimization loop as a memoizing constraint solver. Indicative code for the transcendental differential equations and the optimising constraint solver is provided below:

Transcendental Differential Equation System

equations={k'[t]==-0.1` k[t]-269.3998529790092` (E^(-9.2` (E^(-0.001` t)))) (((8600-2086 (E^(-0.035` t))))^(0.7`)) If[Tat[t]≤0,1,(1/(1+0 Tat[t]+0.0028388` ((Tat[t])^2)))] ((k[t])^(0.3`)) (-1+0.0024598146458005322` (E^(2.433333333333333` (E^(-0.003` t)))) (1+(E^(-0.005` t))) ((μ[t])^(2.8`)))-c[t],Tat'[t]==0.22` (-36.7073066928394`+5.482241155378061` Log[Mat[t]]+(({{\[Piecewise], ({{-0.06`+0.0036` t, t≤100}, {0.3`, True}})}})))-0.06441748633333333` Tat[t]+0.011002153` Tlo[t],Tlo'[t]==0.0048` Tat[t]-0.0048` Tlo[t],Mat'[t]==1.1` (E^(-0.01` t))-0.0190837` Mat[t]+0.009800871837862596` Mup[t]-3.171771574386283` (E^(2.433333333333333` (E^(-0.003` t))-9.2` (E^(-0.001` t)))) (((8600-2086 (E^(-0.035` t))))^(0.7`)) ((k[t])^(0.3`)) (-1.`+μ[t]),Mup'[t]==0.0190837` Mat[t]+0.0003369934177753544` Mlo[t]-0.015203871837862596` Mup[t],Mlo'[t]==-0.0003369934177753544` Mlo[t]+0.005403` Mup[t],
μ'[t]==(1/10) μconst (1-μ[t]) μ[t]};
NDSolve[Flatten[{equations/.{c[t]->38.1,μconst->1.2},k[0]==136.7,Tat[0]==0.7307,Tlo[0]==0.0068,Mat[0]==808.9,Mup[0]==1255.0,Mlo[0]==18365.0,μ[0]==0.005}],{k'[t],k[t],Tat[t],Tlo[t],Mat[t],Mup[t],Mlo[t],μ[t]},{t,0,1400}]

Plot[Evaluate[%[[1,All,2]]],{t,0,1400}]

Memoized Optimizing Constraint Solver

statevars={k[t],Tat[t],Tlo[t],Mat[t],Mup[t],Mlo[t],μ[t]};
solver[params_]:=Check[NDSolve[Flatten[{equations/.Thread[parsarray->params],Thread[(statevars/.t->0)==Last[xvec3][[2]]]}],Join[{k'[t]},statevars],{t,0,testperiods}, MaxSteps->(10^8),InterpolationOrder->All]//Quiet,Print["**** solver failed ****"];{}];
cons0[params_]:=cons0[params]=Module[{optsolver0,optsolver,findmax,findmin,conssol,findopt,finddiff,findval},
optsolver0=solver[params];
optsolver=If[optsolver0=={},Return[{-1,-1,1}(10^30)],optsolver0[[1]]];
findmax[var_,imtestyear_]:=Check[If[MatchQ[var/.optsolver/.{t->imtestyear},_Complex],
 Abs[Im[var/.optsolver/.t->imtestyear]](10^3),FindMaxValue[var/.optsolver,{t,imtestyear,30,testperiods}]]//Quiet,(10^30)];
findmin[var_,imtestyear_]:=Check[If[MatchQ[var/.optsolver/.{t->imtestyear},_Complex],
 -Abs[Im[var/.optsolver/.t->imtestyear]](10^3),FindMinValue[var/.optsolver,{t,imtestyear,40,testperiods}]]//Quiet,-(10^30)];
finddiff[var_,imtestyear_]:=Check[If[MatchQ[var/.optsolver/.t->imtestyear-var/.optsolver/.(t->imtestyear-1),_Complex],
 -Abs[Im[var/.optsolver/.t->imtestyear-var/.optsolver/.(t->imtestyear-1)]],
var/.optsolver/.t->imtestyear-var/.optsolver/.(t->imtestyear-1)]//Quiet,-(10^30)];
findval[var_,imtestyear_]:=Check[If[MatchQ[var/.optsolver/.t->imtestyear,_Complex],
 -Abs[Im[var/.optsolver/.t->imtestyear]],var/.optsolver/.t->imtestyear]//Quiet,-(10^30)];
findopt=((Subscript[cumutil, 0])/.sbsd)+NIntegrate[(Exp[-ρ t](((c[t]/(Subscript[l, t]))^(1-γ))-1)(Subscript[l, t])/((Subscript[scale, 1])(1-γ)))/.sbsd/.Thread[parsarray->params]/.optsolver,{t,0,testperiods},WorkingPrecision->24];
conssol={findopt,findval[k[t],testperiods],findmax[Tat[t],80]}
]//Quiet
cons1[params:{_?NumberQ..}]:=cons0[params][[1]]
cons2[params:{_?NumberQ..}]:=cons0[params][[2]]
cons3[params:{_?NumberQ..}]:=cons0[params][[3]]
starttime=AbsoluteTime[];
stepresult={};stepresult>>"teststep";
Print["Starting Memory Used: "<>ToString[N[MaxMemoryUsed[]/10^6]]<>" Mb"];
sol=Block[{ii=1,jj=1,kk=1,temp=""},
NMaximize[Join[{cons1[optpars],0≤cons2[optpars],0.995 tatmax≤cons3[optpars],cons3[optpars]≤1.005 tatmax},Reverse[Thread[Rest[Most[optpars]]>Most[Most[optpars]]]],{copt1>0,10≥μopt,μopt>0}],Map[{#[[1]],Min[#[[2]],#[[3]]],Max[#[[2]],#[[3]]]}&,Transpose[{optpars,pars*0.995,1.005pars}]],MaxIterations->3500,WorkingPrecision->24,
StepMonitor:>(stepresult={ii,jj,kk,Thread[{Join[{res},Map[ToExpression["c"<>ToString[#]]&,Range[Length[handels]]],{μ}],
Join[{cons1[optpars]},optpars]}]};
stepresult>>>"teststep";jj++;kk=1),EvaluationMonitor:>(ii++;kk++;
If[Mod[ii,100]==0,(*If[Mod[ii,1000]==0,evalresult={ii,jj,kk,Thread[{vars,varsopt}]};
evalresult>>>"m138pNum02eval"];*)DownValues[cons0]=Last[DownValues[cons0]];
NotebookDelete[temp];
temp=PrintTemporary["Step: "<>ToString[jj]<>"; Sub "<>ToString[kk]<>"; Cum "<>ToString[ii]<>"; Curr Mem: "<>ToString[N[MemoryInUse[]/10^6]]<>"; Peak: "<>ToString[N[MaxMemoryUsed[]/10^6]]<>" Mb"<>"; Hours: "<>ToString[N[(AbsoluteTime[]-starttime)/3600]]]])]]//Timing