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>
- Re: Re: Solving nonlinear inequality constraints