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

Solving nonlinear inequality constraints

• To: mathgroup at smc.vnet.net
• Subject: [mg91209] Solving nonlinear inequality constraints
• From: "Stuart Nettleton" <Stuart.Nettleton at eng.uts.edu.au>
• Date: Sat, 9 Aug 2008 07:47:31 -0400 (EDT)

Some advice on my approach here would be much appreciated as I am a Mathematica 
newbie and have not been able to find much on this in the archive. Over  
the last 6 weeks I have been learning Mathematica to work on solving a set of  
linear and nonlinear equalities in discrete time for up to 60 periods. The  
nonlinear equality constraints are relatively straightforward using  
FindRoot. Solve also works for a small number of periods. However, the  
nonlinear inequalities have proven to be quite time consuming. After  
trying may approaches only FindMinimum with Method-> InteriorPoint  
reliably works and I am really impressed with it. With small numbers of  
periods I can get NMaximize and FindMaximum (which doesn't seem to have  
InteriorPoint) to occasionally do the job.

I have success but so far not very efficiently. The processing time is  
0.56 hours for 10 periods, 4.6 hours for 25 periods and something over 24  
hours for 60 periods (which is still running on our university high  
performance cluster of 3Ghz intels as I write). In contrast, I have used  
the student edition of GAMS with a few minutes for 10 periods using the  
CONOPT solver. It says there are 1381 variables, 1263 constraints, 3767  
Jacobean. One author (Nordhaus) claims just 4.5 minutes for 60 periods and  
at the same time doing a fair bit more work in the program. Submitting the  
problem to the NEOS server at http://neos.mcs.anl.gov/neos/ solves 60  
periods in 2.4 minutes. However, these non-Mathematica approaches are not so  
desirable because of the old monolithic style of program and batch  
approaches which lack flexibility and direct connection for the further  
postprocessing of data.

Some advice on whether I have gone about coding the right way would be  
very much appreciated. Also, would there be a better approach to this sort  
of problem in whole or in part? For example, could there be a better  
approach in somehow using the equations as functions rather than  
equalities and inequalities?

Lastly, the advantage in GAMS  seems to be preparing acyclic networks as  
described at http://www.mat.univie.ac.at/~neum/ms/dag.pdf  and, going back  
a few years, Numerica  
(http://mitpress.mit.edu/book-home.tcl?isbn=0262720272 ) A "directed  
acyclic graph" would eliminate all the intermediate variables to map the  
variables to just the unknown leaves of the acyclic tree (in a way similar  
to a spreadsheet topology). This would reduce the number variables to be  
optimised and presumably cut the processing time quite considerably. I  
have tried Reduce but it doesn't seem to cope with the scale of the  
problem, even allowing infinite recursion. Does anyone know of a module I  
could try to map equations in Mathematica?

Many thanks for your consideration of these issues! Stuart

Here is my Notebook, which is at present arranged for unattended command  
line processing of the associated package file:

AppendTo[$Echo,"stdout"];
SetOptions["stdout",PageWidth->110] ; (* most people don't use 80 columns \
anymore; set the default to 110 *)
ZeroTol = 10^-10 ;
(*PrecisionTol = $MachinePrecision ;*)
PrecisionTol = $MachinePrecision ;
(*AccuracyGoalTol = PrecisionTol*.5;
PrecisionGoalTol = AccuracyGoalTol*.5;*)

(* Nordhaus Brief Climate Change Policy Model May 2008 *)
(* Stuart \
Nettleton July 2008 *)
(*ClearSystemCache[]*)

starttime=AbsoluteTime[];
periods = 5 ; (* projection periods *)
maxit = 5000; (* maximum iterations *)
\
objvars = {-cumu[periods]}; (* Objective function *)

(* Exogenous variables *)
exogvars={
   gfacpop[t] == (Exp[popg*(t-1)]-1)/Exp[popg*(t-1)],
   l[t]      == pop0*(1-gfacpop[t])+gfacpop[t]*popa,
   ga[t]     == ga[0]*Exp[-dela*10*(t-1)],
   a[t]      == If[t==1,a[0],a[t-1]/(1-ga[t-1])],
   g\[Sigma][t]     == \
g\[Sigma][0]*Exp[-d\[Sigma]1*10*(t-1)-d\[Sigma]2*10*(t-1)^2],
   \[Sigma][t]      == If[t==1,\[Sigma][0],\[Sigma][t-1]/(1-g\[Sigma][t])],
   \[CapitalTheta][t]      == \
(pback*\[Sigma][t]/\[Theta])*((backrat-1+Exp[-gback*(t-1)])/backrat),
   eland[t]  == eland[0]*(1-0.1)^(t-1),
   r[t]      == 1/(1+\[Rho])^(10*(t-1)),
   fex[t]    == fex0 + If[t<12,0.1*(fex1-fex0)*(t-1),0.36],
   \[Kappa][t]      == If[t==1,\[Kappa][0],If[t>= 25, \[Kappa]21,  
\[Kappa]21 + \
(\[Kappa]2- \[Kappa]21)*Exp[-d\[Kappa]*(t-2)]]],
   \[CapitalPi][t]      == \[Kappa][t]^(1-\[Theta]),
   s[t]      == sr
};

(* Constraints, note that all variables should start with lower *)
(* case to \
avoid Mathematica's built-in functions *)
endovars={
   ceind[t+1] == eind[t]+ceind[t],
   k[t+1] <= 10*inv[t]+((1-\[Delta])^10)*k[t],
   0.02*k[periods] <= inv[periods],
   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+1] == e[t] +\[Phi]11*mat[t] +\[Phi]21*mup[t],
   matav[t] == (mat[t]+mat[t+1])/2,
   mlo[t+1] == \[Phi]23*mup[t]+\[Phi]33*mlo[t],
   mup[t+1] == \[Phi]12*mat[t]+\[Phi]22*mup[t] +\[Phi]32*mlo[t],
   tat[t+1] ==  
tat[t]+\[Xi]1*(for[t+1]-\[Xi]2*tat[t]-\[Xi]3*(tat[t]-tlo[t])),
   tlo[t+1] == tlo[t]+\[Xi]4*(tat[t]-tlo[t]),
   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)),
   s[t] == inv[t]/(0.001+y[t]),
   ri[t] == \[Gamma]*y[t]/k[t] -(1-(1-\[Delta])^10)/10,
   c[t] == y[t]-inv[t],
   (*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+1] == cumu[t]+(l[t+1]*u[t+1]*r[t+1]*10)/scale1,
   100 <= k[t],
   20 <= c[t],
   10 <= mat[t],
   100 <= mup[t],
   1000 <= mlo[t],
   tlo[t] <= 20,
   -1 <= tlo[t],
   tat[t] <= 20,
   \[Mu][t] <= \[Mu]lim,
   ceind[t] <= ceindlim
};

   (*\[CapitalLambda][t] == \[CapitalPi][t] *\[CapitalTheta][t] *\[Mu][t]^\
\[Theta],
   \[CapitalOmega][t] == 1/(1+\[Psi]1*tat[t]+\[Psi]2*(tat[t]^\[Psi]3)),
   y[t] == \[CapitalOmega][t]*(1-\[CapitalLambda][t])*ygr[t]*)

(* positive variable conditions *)
\
posvars={q[t],inv[t],ygr[t],eind[t],matav[t],\[Mu][t]};

(* parameter values *)
parametervals={
     \[Alpha]-> 2.0,
     \[Gamma]-> 0.30,
     \[Delta]-> 0.1,
     \[Rho]-> 0.015,
     \[Eta]-> 3.8,
     \[Theta]-> 2.8,
     \[Psi]1-> 0.00000,
     \[Psi]2-> 0.0028388,
     \[Psi]3-> 2.00,

     \[Xi]1-> 0.220,
     (*\[Xi]2-> 3.8/3,*)
     \[Xi]2-> \[Eta]/t2xco2,
     t2xco2 ->3,
     \[Xi]3-> 0.300,
     \[Xi]4-> 0.050,

     pop0-> 6514,
     popg-> 0.35,
     popa-> 8600,
     dela-> 0.001,
     d\[Sigma]1-> 0.003,
     d\[Sigma]2-> 0.000,

     \[Phi]11-> 1-\[Phi]12a,
     (*\[Phi]11-> 1-0.18928,*)
     (* \[Phi]11->0.810712, *)
     \[Phi]12a-> \
0.189288,
     \[Phi]12-> 0.189288,
     (* \[Phi]21->0.097213, *)
     (* \[Phi]22->0.852787, *)
     \[Phi]23-> \
0.05,
     \[Phi]23a-> 0.05,
     (* \[Phi]32->0.003119, *)
     (* \[Phi]33->0.996881, *)
     \[Phi]21-> \
587.473*\[Phi]12a/1143.894,
     (*\[Phi]21-> 587.473*0.189288/1143.894,*)
     \[Phi]22-> 1-\[Phi]21 - \
\[Phi]23a,
     (*\[Phi]22-> 1-(587.473*0.189288/1143.894) - 0.05,*)
     \[Phi]32-> \
1143.894*\[Phi]23a/18340,
     (*\[Phi]32-> 1143.894*0.05/18340,*)
     \[Phi]33-> 1-\[Phi]32,
     (*\[Phi]33-> 1-(1143.894*0.05/18340),*)

     fex0-> -0.06,
     fex1-> 0.30,
     mat1750-> 596.4,
     pback-> 1.17,
     backrat-> 2,
     gback-> 0.05,
     \[Kappa]1 ->1,
     \[Kappa]2 ->1,
     \[Kappa]21 ->1,
     d\[Kappa] ->0,
     \[Mu]lim-> 1,
     ceindlim-> 6000,
     scale1-> 194,
     sr-> 0.22 (* savings ratio *)
};

(* initial values *)
initialvalues={
     a[0]-> 0.02722,
     k[0]-> 137,
     y[0]-> 61.1,
     c[0]-> 30,
     eland[0]-> 11,
     ga[0]-> 0.092,
     \[Sigma][0]-> 0.13418,
     g\[Sigma][0]-> -0.0730,
     mat[0]-> 808.9,
     mup[0]-> 1255,
     mlo[0]-> 18365,
     tat[0]-> 0.7307,
     tlo[0]-> 0.0068,
     \[Kappa][0]-> 0.25372,
     \[Mu][0]-> 0.005,
     ceind[0]->0,
     cumu[0]->381800,
     ceind[1] -> ceind[0],
     k[1] -> k[0],
     mat[1]-> mat[0],
     mup[1] -> mup[0],
     mlo[1] -> mlo[0],
     tat[1] -> tat[0],
     tlo[1] -> tlo[0],
     \[Mu][1] -> \[Mu][0],
     cumu[1]->cumu[0]
};

(* Data conversions based on perturbationAIM by Gary Anderson and Eric  
Swanson \
http://www.ericswanson.org *)

VarNum1[equations_]:= VarNum1[equations] = Length[VarNames1[equations]] ;
VarNum2[equations_]:= VarNum2[equations] = Length[VarNames2[equations]] ;
VarNum3[equations_]:= VarNum3[equations] = Length[VarNames3[equations]] ;
VarNames1[equations_]:= VarNames1[equations] = \
Union[Cases[equations,x_Symbol[t+i_Integer]-> x,Infinity]] ;
VarNames2[equations_]:= VarNames2[equations] = \
Union[Cases[equations,x_Symbol[i_Integer]-> x[i],Infinity]] ;
VarNames3[equations_]:= VarNames3[equations] = \
Complement[Complement[Union[Cases[equations,x_Symbol[i_Integer]-> x[i],
                         Infinity]],variables],VarNames2[initialvalues]] ;
Vars1[equations_]:= Vars1[equations] = \
Union[Cases[equations,x_Symbol[i_Integer] -> x[i],Infinity]] ;
Vars2[equations_]:= Vars2[equations] = \
Union[Cases[equations,x_Symbol[i_Integer] -> x[i]>=0,Infinity]] ;
MaxLag[equations_]:= MaxLag[equations] = \
-Min[Cases[equations,x_Symbol[t+i_Integer]->i,Infinity]] ;
MaxLead[equations_]:= MaxLead[equations] = \
Max[Cases[equations,x_Symbol[t+i_Integer]->i,Infinity]] ;
maxlead=Max[0,MaxLead[endovars]];
SNSS[periods_,equations_,initialvalues_List,  
parameters_List,opts___Rule]:= \
SNSS[periods,equations,initialvalues, parameters,ZeroTol,PrecisionTol] =
  Module[{ssguess,sseqns,sseqns11,sseqns22,sseqns33,sseqns44,sseqns55},
   sseqns = equations ; (* /.Equal->Subtract ; *)
   (* Print["I am sseqns \
",sseqns]; *)
   sseqns11 = Map[sseqns /.t -> # &, Range[1, periods + \
maxlead]];
   sseqns22 = Union[Flatten[sseqns11]];
   (* Print["I am sseqns22 ",sseqns22]; *)
   sseqns33 = sseqns22 //.parameters \
; (* substitute parameters so starting conditions are entered *)
   (* \
Print["I am sseqns33 ",sseqns33]; *)
   sseqns44 = sseqns33 //.initialvalues ; \

   (* Print["I am sseqns44 ",sseqns44] ; *)
   sseqns55 = sseqns44 \
/.x_Symbol[i_Integer /;i <= 0] -> 0 ;
   (* Print["I am sseqns55 ",sseqns55] ; *)
   sseqns55
];

exogvarsext=SNSS[periods, exogvars, initialvalues, \
parametervals]/.x_Real:>SetPrecision[x,PrecisionTol];
(* exogenous = \
Block[{$RecursionLimit=Infinity},Flatten[Solve[exogvarsext,Vars1[exogvarsext]]\
]]; *)
exoguess=Transpose[{Vars1[exogvarsext],Table[.2,{i,Length[exogvarsext]}\
]}];
(*exogenous = FindRoot[exogvarsext,exoguess, MaxIterations -> 5000, \
WorkingPrecision->PrecisionTol, AccuracyGoal->AccuracyGoalTol, \
PrecisionGoal->PrecisionGoalTol];*)
exogenous = \
FindRoot[exogvarsext,exoguess, MaxIterations -> maxit, \
WorkingPrecision->PrecisionTol];
interimtime1=AbsoluteTime[];
(*Print["exogenous completed in: ",interimtime1-starttime];*)
(* substitute \
values into equations: *)
objective = objvars \
//.Join[parametervals,exogenous]; (*/.x_Real:>SetPrecision[x,PrecisionTol]  
;*)
\
endogenous = SNSS[periods, endovars, initialvalues, \
Join[parametervals,exogenous]];
interimtime2=AbsoluteTime[];
(*Print["endogenous completed in: ",interimtime2-interimtime1];*)
\
(*Print["elapsed time: ",interimtime2-starttime];*)
\
variables=Vars1[endogenous];
Print["number of variables: ",Length[variables]];
positivevars = Vars2[Intersection[SNSS[periods,posvars,{},{}],variables]];
interimtime3=AbsoluteTime[];
(*Print["positivevars completed in: ",interimtime3-interimtime2];*)
\
(*Print["elapsed time: ",interimtime3-starttime];*)

finalvars=Join[objective,endogenous,positivevars];
soln=Block[{$RecursionLimit=Infinity},FindMinimum[finalvars,variables,\
MaxIterations->maxit,Method->"InteriorPoint"]];
interimtime4=AbsoluteTime[];
outfilename="temp-"<>DateString[{"YearShort","Month","Day"}]<>"-nmax-inpoint-\
it"<>ToString[maxit]<>"-p"<>ToString[periods];
Put[soln,outfilename];
interimtime4=AbsoluteTime[];
timedata1="FindMinimum completed in: \
"<>ToString[interimtime4-interimtime3]<>" seconds or \
"<>ToString[(interimtime4-interimtime3)/3600]<>" hours or \
"<>ToString[(interimtime4-interimtime3)/60]<>" minutes";
PutAppend[timedata1,outfilename];
timedata2="Total elapsed time      :  
"<>ToString[interimtime4-starttime]<>" \
seconds or "<>ToString[(interimtime4-starttime)/3600]<>" hours or \
"<>ToString[(interimtime4-starttime)/60]<>" minutes";
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.

• Follow-Ups:
  • Re: Re: Solving nonlinear inequality constraints
    • From: "Stuart Nettleton" <Stuart.Nettleton@uts.edu.au>
  • Re: Re: Solving nonlinear inequality
    • From: danl@wolfram.com
  • Re: Solving nonlinear inequality constraints
    • From: "Stuart Nettleton" <Stuart.Nettleton@uts.edu.au>
  • Re: Solving nonlinear inequality constraints
    • From: "Stuart Nettleton" <Stuart.Nettleton@uts.edu.au>
  • Re: Solving nonlinear inequality constraints
    • From: "Stuart Nettleton" <Stuart.Nettleton@uts.edu.au>