Re: Re: Solving nonlinear inequality constraints
- To: mathgroup at smc.vnet.net
- Subject: [mg91446] Re: [mg91243] Re: [mg91209] Solving nonlinear inequality constraints
- From: "Stuart Nettleton" <Stuart.Nettleton at uts.edu.au>
- Date: Fri, 22 Aug 2008 03:11:31 -0400 (EDT)
- Organization: University of Technology, Sydney
- References: <op.uf6ehqonidel0u@linux-vmzf.site> <48ADADDA.1010109@wolfram.com>
Thanks, Daniel, for your suggestions below which I very much appreciate.
Fortunately William Nordhaus has provided most of the domain specific
information we need for these equations. You are quite correct that s[0]
gets a bit out of hand. It is an artifact of me "releasing s", which I did
to demonstrate that it was more than just mu[t] that needed to be
optimised (as Nordhaus does). Probably Nordhaus doesn't correct for s[0]
because it is a beginning effect that doesn't adversely affect his solver
(conopt) and he never subsequently uses s[0] as a meaningful data point.
In our case it is quite appropriate to set an initial value for s[0] and
s[1] which I have done in the slightly amended notebook below. This
notebook also takes up your suggestion that the simple variables are
precalculated rather than left to be plumbed. This is done by defining a
list of variables and evaluating them in the first part of the processing
phase. It only takes 0.4 second and 6Mb to do either 10 or 60 periods. I
had high hopes for these changes on my first run of 10 periods because the
compile and processing time were cut by 2/3! However, memory is still the
limiting factor and the new efficiency wasn't sufficient to permit 11
periods to be evaluated in 16Gb.
Using the _NumberQ version of the notebook (below), the comparison
statistics are:
Approach Periods Compilation Total Memory
Minutes
Minutes Mb
Old plumbed 10 8.3 mins 66.8
8901
New precalculated 10 6.0 mins 21.8 9051
New precalculated 11 23.9 mins Ran out of memory
Thanks also for looking at Nordhaus' equations and book! This really
exceeds anything I could imagine in help on this group!
I am proceeding to examine your suggestion of using rational numbers for
"simple" decimals and will let you know what I find.
Nevertheless, I wonder if we are getting down to secondary issues and
missing a lurking primary issue that is causing Mm to use so much memory?
Many thanks for your continued interest!
Stuart
> I cannot really help you much. The model requires serious
> domain-specific knowledge even to begin to check for, say,
> implementation errors. I can only make a few general remarks that might
> help a bit.
>
> I strongly suspect there are in fact issues with your code
> implementation. Having run the case of 5 periods, and looked at both
> constraints/objective function and the results, I see two red flags.
>
> First is that the result has a huge value for s[0].
>
> In[268]:= (*prepare optimising function*)
> Timing[
> soln = NMinimize[Join[objectimous, endogenous], extendedvariables,
> MaxIterations -> 50];]
>
> Out[268]= {50.8672, Null}
>
> In[269]:= soln
>
> Out[269]= {-283917., {k[2] -> 403.596, k[3] -> 319.585,
> k[4] -> 278.518, k[5] -> 208.303, s[0] -> 6.66695*10^6,
> s[1] -> 0.640164, s[2] -> 0.198601, s[3] -> 0.173787,
> s[4] -> 0.106845,
> s[5] -> 0.039112, \[Mu][2] -> 0.129617, \[Mu][3] ->
> 0.113797, \[Mu][4] -> 0.0955035, \[Mu][5] -> 0.0550562}}
>
> And this remains the case when I allow more iterations, use fifferent
> methods, etc. So I do not think it is an artifact of the NMinimize
> internals.
>
> The second issue is that the objective, and some of the constraints, are
> deeply nested expressions of their variables. I realize that the
> recurrence nature of the definitions can allow this to happen. But I do
> not believe that Nordhaus' implementation can suffer from this, or else
> I believe he'd have no serious chance of handling 600 years (60 periods)
> using any software.
>
> The upshot is I would recommend careful proof reading of the system.
> I've looked a bit at p 2008 of his "The Challenge of Global Warming..."
> and confess the equations for y[t], etc. do indeed look dauntingly
> recursive. But still I suspect there is something he does to unravel the
> mess.
>
> Finally, if any of the equalities are simple e.g. linear in some
> variable(s), I'd say solve for such variables, use that to plug into
> everything else, in the hope that some of the deep nestedness might be
> removable.
>
> A last recommendation is that you use rational numbers for "simple"
> decimals like 0.0 or .05. This may help to collapse things a bit, thus
> possibly reducing the complexity.
>
> Daniel
(*AppendTo[$Echo,"stdout"];*)
(*SetOptions["stdout",PageWidth->110] ; \
*)
(* Nordhaus Brief Climate Change Policy Model May 2008 *)
(* \
Stuart Nettleton July 2008 *)
starttime = AbsoluteTime[];
periods = 10; (* projection periods *)
maxit = 5000; (* maximum \
iterations *)
(* objective function *)
(* program always minimises, so negative for \
maximisation *)
obj = {-cumu[periods]};
(* optimisation variables: topolology leaves *)
(* .. for NMinimize \
provide both upper and lower bounds *)
(* .. for FindMinimum provide \
a single start estimate *)
(* hold function prevents premature \
evaluation *)
(*opt = {{\[Mu][t],0.001,1},{k[t],100,10000}};*)
opt = {{\[Mu][t], 0.001, 1}};
(* exogenous parameters *)
pop0 = 6514;
popg = 0.35;
popa = 8600;
dela = 0.001;
d\[Sigma]1 = 0.003;
d\[Sigma]2 = 0.000;
pback = 1.17;
backrat = 2;
gback = 0.05;
\[Rho] = 0.015;
fex0 = -0.06;
fex1 = 0.30;
\[Kappa]1 = 1;
\[Kappa]2 = 1;
\[Kappa]21 = 1;
d\[Kappa] = 0;
\[Theta] = 2.8;
sr = 0.22;
(* exogenous variables *)
exogvars = {gfacpop, l, ga, a, g\[Sigma], \[Sigma], \[CapitalTheta],
eland, r, fex, \[Kappa], \[CapitalPi]};
gfacpop[t_] :=
gfacpop[t] = (Exp[popg*(t - 1)] - 1)/Exp[popg*(t - 1)];
l[t_] := l[t] = pop0*(1 - gfacpop[t]) + gfacpop[t]*popa;
ga[t_] := ga[t] = ga[0]*Exp[-dela*10*(t - 1)];
a[t_] := a[t] = a[t - 1]/(1 - ga[t - 1]);
g\[Sigma][t_] :=
g\[Sigma][t] =
g\[Sigma][0]*
Exp[-d\[Sigma]1*10*(t - 1) - d\[Sigma]2*10*(t - 1)^2];
\[Sigma][t_] := \[Sigma][t] = \[Sigma][t - 1]/(1 - g\[Sigma][t]);
\[CapitalTheta][
t_] := \[CapitalTheta][
t] = (pback*\[Sigma][t]/\[Theta])*((backrat - 1 +
Exp[-gback*(t - 1)])/backrat);
eland[t_] := eland[t] = eland[0]*(1 - 0.1)^(t - 1);
r[t_] := r[t] = 1/(1 + \[Rho])^(10*(t - 1));
fex[t_] :=
fex[t] = fex0 + If[t < 12, 0.1*(fex1 - fex0)*(t - 1), 0.36];
\[Kappa][t_] := \[Kappa][t] =
If[t >= 25, \[Kappa]21, \[Kappa]21 + (\[Kappa]2 - \[Kappa]21)*
Exp[-d\[Kappa]*(t - 2)]];
\[CapitalPi][t_] := \[CapitalPi][t] = \[Kappa][t]^(1 - \[Theta]);
(*s[t_]:= s[t] =sr;*)
(* initial values of exogenous variables *)
s[1] = s[0] = sr;
ga[0] = 0.092;
g\[Sigma][0] = -0.0730;
eland[0] = 11;
a[1] = a[0] = 0.02722;
\[Sigma][1] = \[Sigma][0] = 0.13418;
\[Kappa][1] = \[Kappa][0] = 0.25372;
(* endogenous parameters *)
\[Alpha] = 2.0;
\[Gamma] = .30;
\[Delta] = 0.1;
\[Eta] = 3.8;
t2xco2 = 3;
\[Psi]1 = 0.00000;
\[Psi]2 = 0.0028388;
\[Psi]3 = 2.00;
\[Xi]1 = 0.220;
\[Xi]2 = \[Eta]/t2xco2;
\[Xi]3 = 0.300;
\[Xi]4 = 0.050;
\[Phi]12 = 0.189288;
\[Phi]12a = 0.189288;
\[Phi]11 = 1 - \[Phi]12a;
\[Phi]23 = 0.05;
\[Phi]23a = 0.05;
\[Phi]21 = 587.473*\[Phi]12a/1143.894;
\[Phi]22 = 1 - \[Phi]21 - \[Phi]23a;
\[Phi]32 = 1143.894*\[Phi]23a/18340;
\[Phi]33 = 1 - \[Phi]32;
mat1750 = 596.4;
k0 = 137;
y0 = 61.1;
c0 = 30;
mat0 = 808.9;
mup0 = 1255;
mlo0 = 18365;
tat0 = 0.7307;
tlo0 = 0.0068;
\[Mu]0 = 0.005;
ceind0 = 0;
cumu0 = 381800;
scale1 = 194;
(* endogenous equality constraints*)
ceind[t_?NumberQ] := eind[t - 1] + ceind[t - 1];
eind[t_?NumberQ] :=
10 *\[Sigma][t] *(1 - \[Mu][t]) *ygr[t] + eland[t];
for[t_?NumberQ] := \[Eta]*(Log[(matav[t] + 0.000001)/mat1750]/
Log[2]) + fex[t];
mat[t_?NumberQ] :=
eind[t - 1] + \[Phi]11*mat[t - 1] + \[Phi]21*mup[t - 1];
matav[t_?NumberQ] := (mat[t] + mat[t + 1])/2;
mlo[t_?NumberQ] := \[Phi]23*mup[t - 1] + \[Phi]33*mlo[t - 1];
mup[t_?NumberQ] := \[Phi]12*mat[t - 1] + \[Phi]22*
mup[t - 1] + \[Phi]32*mlo[t - 1];
tat[t_?NumberQ] :=
tat[t - 1] + \[Xi]1*(for[t] - \[Xi]2*
tat[t - 1] - \[Xi]3*(tat[t - 1] - tlo[t - 1]));
tlo[t_?NumberQ] := tlo[t - 1] + \[Xi]4*(tat[t - 1] - tlo[t - 1]);
ygr[t_?NumberQ] := a[t]* k[t]^\[Gamma] *l[t]^(1 - \[Gamma]);
dam[t_?NumberQ] :=
ygr[t]*(1 - 1/(1 + \[Psi]1*tat[t] + \[Psi]2*(tat[t]^\[Psi]3)));
\[CapitalLambda][t_?NumberQ] :=
ygr[t] *\[CapitalPi][t] *\[CapitalTheta][t] *\[Mu][t]^\[Theta];
y[t_?NumberQ] :=
ygr[t]*(1 - \[CapitalPi][t]*\[CapitalTheta][
t]*\[Mu][t]^\[Theta])/(1 + \[Psi]1*
tat[t] + \[Psi]2*(tat[t]^\[Psi]3));
inv[t_?NumberQ] := (y[t] + 0.001)*s[t];
(*k[t_?NumberQ]:=10*inv[t-1]+((1-\[Delta])^10)*k[t-1];*)
ri[t_?NumberQ] := \[Gamma]*y[t]/k[t] - (1 - (1 - \[Delta])^10)/10;
c[t_?NumberQ] := y[t] - inv[t];
u[t_?NumberQ] := ((c[t]/l[t])^(1 - \[Alpha]) - 1)/(1 - \[Alpha]);
cumu[t_?NumberQ] := cumu[t - 1] + (l[t]*u[t]*r[t]*10)/scale1;
cpc[t_?NumberQ] := c[t]*1000/l[t];
pcy[t_?NumberQ] := y[t]*1000/l[t];
(* initial values of endogenous variables *)
cumu[1] = cumu[0] = cumu0;
ceind[1] = ceind[0] = ceind0;
mat[1] = mat[0] = mat0;
mlo[1] = mlo[0] = mlo0;
mup[1] = mup[0] = mup0;
tat[1] = tat[0] = tat0;
tlo[1] = tlo[0] = tlo0;
y[0] = y0;
k[1] = k[0] = k0;
c[0] = c0;
\[Mu][1] = \[Mu][0] = \[Mu]0;
(* endogenous inequality constraints*)
endog = {
k[t] <= 10*inv[t - 1] + ((1 - \[Delta])^10)*k[t - 1],
0.02*k[periods] <= inv[periods],
100 <= k[t],
20 <= c[t],
10 <= mat[t],
100 <= mup[t],
1000 <= mlo[t],
-1 <= tlo[t] <= 20,
0 <= tat[t] <= 20,
ceind[t] <= 6000,
0 <= y[t],
0 <= inv[t],
0 <= ygr[t],
0 <= eind[t],
0 <= matav[t],
0 <= \[Mu][t]
};
(* processing ... *)
(* calculate exogenous variables *)
Do[Map[exogvars[[i]], Range[periods]], {i, Range[Length[exogvars]]}];
(*prepare the objective function*)
objvars = Simplify[obj];
objectimous =
Simplify[
Flatten[Union[Map[objvars /. t -> # &, Range[periods]]] /.
x_Symbol[i_Integer /; i < 0] -> 0]];
(*prepare the endogenous constraints*)
endovars = Simplify[endog];
endogenous =
Simplify[
Flatten[Union[Map[endovars /. t -> # &, Range[periods]]] /.
x_Symbol[i_Integer /; i < 0] -> 0]];
(*prepare the independent optimising variables*)
optvars = Simplify[opt];
optimous =
Union[Cases[
Flatten[Map[optvars /. t -> # &, Range[periods]],
1], {x_Symbol[_Integer], _Integer | _Real} | \
{x_Symbol[_Integer], _Integer | _Real, _Integer | _Real} | \
{x_Symbol[_Integer], _Integer | _Real, _Integer | _Real, _Integer | \
_Real}, Infinity]];
optimousvars =
Union[Cases[optimous, x_Symbol[i_Integer] -> x[i], Infinity]];
(*include any additional optimising variables arising from the \
endogenous constraints*)
primafacievars =
Union[Cases[objectimous, x_Symbol[i_Integer] -> x[i], Infinity],
Cases[endogenous, x_Symbol[i_Integer] -> x[i], Infinity]];
extendedvariables =
Union[optimous, Complement[primafacievars, optimousvars]];
(*prepare optimising function*)
Print[Round[(AbsoluteTime[] - starttime)/60, 0.1],
" minutes compilation time, commencing optimisation ..."];
soln = NMinimize[Join[objectimous, endogenous], extendedvariables,
MaxIterations -> maxit];
(*output results*)
Print["Solution: ", soln];
Print[Round[(AbsoluteTime[] - starttime)/60, 0.1],
" minutes from start to completion"];
Print[Round[N[MaxMemoryUsed[]*10^-6], 1], " Mb memory used"];
--
UTS CRICOS Provider Code: 00099F
DISCLAIMER: This email message and any accompanying attachments may contain
confidential information. If you are not the intended recipient, do not
read, use, disseminate, distribute or copy this message or attachments. If
you have received this message in error, please notify the sender
immediately and delete this message. Any views expressed in this message
are those of the individual sender, except where the sender expressly, and
with authority, states them to be the views the University of Technology,
Sydney. Before opening any attachments, please check them for viruses and
defects.