Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2008

[Date Index] [Thread Index] [Author Index]

Search the Archive

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.


  • Prev by Date: Re: More Inquiries
  • Next by Date: Re: Printing and Exporting graphics images from Mathematica 6 to PDF?
  • Previous by thread: Re: Solving nonlinear inequality constraints
  • Next by thread: Re: Re: Solving nonlinear inequality