Re: Solving nonlinear inequality constraints
- To: mathgroup at smc.vnet.net
- Subject: [mg91243] Re: [mg91209] Solving nonlinear inequality constraints
- From: "Stuart Nettleton" <Stuart.Nettleton at uts.edu.au>
- Date: Mon, 11 Aug 2008 06:06:28 -0400 (EDT)
- Organization: University of Technology, Sydney
- References: <200808091147.HAA18525@smc.vnet.net>
I have taken a completely different approach by using active data for
auto-topology and using NMaximize to solve only for the leaves. This means
small numbers of periods calculate very quickly. So after setting a 25
period version running on our university high performance cluster I went
to sleep at 3am. Waking up I felt jubilant at the elegance, truth and
beauty in my new solution. However, I was dismayed to find the Mathematica
blue screen of death: "No more memory available. Mathematica kernel has shut
down. Try quitting other applications and then retry". Of course I was the
only person on the cluster node and all I had running was Mathematica. I
am not sure of the memory size on our nodes so I will check this. With Red Hat
linux it should be a Windows problem.
Would somebody be able to help me by having a brief look at this code and
see if I have made an error somewhere or if the code can be made both
beautiful and run for up to 60 periods?
Many thanks,
Stuart
AppendTo[$Echo, "stdout"];
SetOptions["stdout", PageWidth -> 110] ;
(* Nordhaus Brief Climate Change Policy Model May 2008 *)
(* Stuart \
Nettleton July 2008 *)
(* ClearSystemCache[] *)
starttime = AbsoluteTime[];
periods = 25 ; (* projection periods *)
maxit = 5000; (* maximum \
iterations *)
(* objective function *)
(* program always minimises, so negative for \
maximisation *)
(* hold function prevents premature evaluation *)
obj \
= Hold[-cumu[periods]];
objvars = Apply[List, Map[Hold, obj]];
(* 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 = Hold[{\[Mu][t],0.001,1},{k[t],100,10000}];*)
opt = Hold[{\[Mu][t], 0.001, 1}];
optvars = Apply[List, Map[Hold, opt]];
(* exogenous parameters and variables *)
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; (* savings ratio *)
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)]; ga[0] = 0.092;
a[t_] := a[t] = a[t - 1]/(1 - ga[t - 1]); a[1] = a[0] = 0.02722;
g\[Sigma][t_] :=
g\[Sigma][t] =
g\[Sigma][0]*Exp[-d\[Sigma]1*10*(t - 1) - d\[Sigma]2*10*(t - 1)^2];
g\[Sigma][0] = -0.0730;
\[Sigma][
t_] := \[Sigma][t] = \[Sigma][t - 1]/(1 - g\[Sigma][t]); \[Sigma][
1] = \[Sigma][0] = 0.13418;
\[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); eland[0] = 11;
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)]]; \[Kappa][1] = \[Kappa][0] = 0.25372;
\[CapitalPi][t_] := \[CapitalPi][t] = \[Kappa][t]^(1 - \[Theta]);
s[t_] := s[t] = sr;
(* endogenous parameters and variables *)
\[Alpha] = 2.0;
\[Gamma] = 0.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;
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_] := eind[t - 1] + ceind[t - 1];
ceind[1] = ceind[0] = ceind0;
eind[t_] := 10 *\[Sigma][t] *(1 - \[Mu][t]) *ygr[t] + eland[t];
for[t_] := \[Eta]*(Log[(matav[t] + 0.000001)/mat1750]/Log[2]) +
fex[t];
mat[t_] := eind[t - 1] + \[Phi]11*mat[t - 1] + \[Phi]21*mup[t - 1];
mat[1] = mat[0] = mat0;
matav[t_] := (mat[t] + mat[t + 1])/2;
mlo[t_] := \[Phi]23*mup[t - 1] + \[Phi]33*mlo[t - 1];
mlo[1] = mlo[0] = mlo0;
mup[t_] := \[Phi]12*mat[t - 1] + \[Phi]22*mup[t - 1] + \[Phi]32*
mlo[t - 1]; mup[1] = mup[0] = mup0;
tat[t_] :=
tat[t - 1] + \[Xi]1*(for[t] - \[Xi]2*
tat[t - 1] - \[Xi]3*(tat[t - 1] - tlo[t - 1]));
tat[1] = tat[0] = tat0;
tlo[t_] := tlo[t - 1] + \[Xi]4*(tat[t - 1] - tlo[t - 1]);
tlo[1] = tlo[0] = tlo0;
ygr[t_] := a[t]* k[t]^\[Gamma] *l[t]^(1 - \[Gamma]);
dam[t_] :=
ygr[t]*(1 - 1/(1 + \[Psi]1*tat[t] + \[Psi]2*(tat[t]^\[Psi]3)));
\[CapitalLambda][t_] :=
ygr[t] *\[CapitalPi][t] *\[CapitalTheta][t] *\[Mu][t]^\[Theta];
y[t_] :=
ygr[t]*(1 - \[CapitalPi][t]*\[CapitalTheta][
t]*\[Mu][t]^\[Theta])/(1 + \[Psi]1*
tat[t] + \[Psi]2*(tat[t]^\[Psi]3)); y[0] = y0;
inv[t_] := (y[t] + 0.001)*s[t];
k[t_] := 10*inv[t - 1] + ((1 - \[Delta])^10)*k[t - 1];
k[1] = k[0] = k0;
ri[t_] := \[Gamma]*y[t]/k[t] - (1 - (1 - \[Delta])^10)/10;
c[t_] := y[t] - inv[t]; c[0] = 30;
cpc[t_] := c[t]*1000/l[t];
pcy[t_] := y[t]*1000/l[t];
u[t_] := ((c[t]/l[t])^(1 - \[Alpha]) - 1)/(1 - \[Alpha]);
cumu[t_] := cumu[t - 1] + (l[t]*u[t]*r[t]*10)/scale1;
cumu[1] = cumu[0] = cumu0;
\[Mu][1] = \[Mu][0] = \[Mu]0;
(* endogenous inequality constraints *)
(* hold function prevents \
premature evaluation *)
endog = Hold[
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,
10 <= mat[t],
ceind[t] <= 6000,
0 <= y[t],
0 <= inv[t],
0 <= ygr[t],
0 <= eind[t],
0 <= matav[t],
0 <= \[Mu][t]
];
endovars = Apply[List, Map[Hold, endog]];
(* processing ... *)
optimous = Flatten[Map[optvars /. t -> # &, Range[1, periods]]];
endogenous = Flatten[Map[endovars /. t -> # &, Range[1, periods]]];
(* global minimisation ... *)
soln = NMinimize[
Apply[List, Map[ReleaseHold, Join[objvars, endogenous]]],
Union[Cases[Apply[List, Map[ReleaseHold, optimous]],
{x_Symbol[_Integer], _Integer | _Real, _Integer | _Real} | \
{x_Symbol[_Integer], _Integer | _Real},
Infinity]], MaxIterations -> maxit];
(* local minimisation ... *)
\
(*soln=FindMinimum[Apply[List,Map[ReleaseHold,Join[objvars,endogenous]\
]],
Union[Cases[Apply[List,Map[ReleaseHold,optimous]],
{x_Symbol[_Integer],_Integer|_Real,_Integer|_Real}|{x_Symbol[_Integer]\
,_Integer|_Real},
Infinity]],MaxIterations->maxit,Method->"InteriorPoint"];*)
(* print output to console ... *)
Print[soln];
interimtime4 = AbsoluteTime[];
(*timedata1="FindMinimum completed in: \
"<>ToString[interimtime4-interimtime3]<>" seconds or \
"<>ToString[(interimtime4-interimtime3)/3600]<>" hours or \
"<>ToString[(interimtime4-interimtime3)/60]<>" minutes";*)
timedata2 =
"Total elapsed time : " <> ToString[interimtime4 - starttime] <>
" seconds or " <> ToString[(interimtime4 - starttime)/3600] <>
" hours or " <> ToString[(interimtime4 - starttime)/60] <>
" minutes";
Print[timedata2];
(* print output to file ... *)
(*PutAppend[timedata2,outfilename];*)
outfilename =
"temp-" <> DateString[{"YearShort", "Month", "Day"}] <>
"-nmax-inpoint-it" <> ToString[maxit] <> "-p" <>
ToString[periods];
Put[soln, outfilename];
(*PutAppend[timedata1,outfilename];*)
\
(*PutAppend[timedata2,outfilename];*)
--
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.
- References:
- Solving nonlinear inequality constraints
- From: "Stuart Nettleton" <Stuart.Nettleton@eng.uts.edu.au>
- Solving nonlinear inequality constraints