Re: Re: Solving nonlinear inequality constraints
- To: mathgroup at smc.vnet.net
- Subject: [mg92386] Re: [mg91243] Re: [mg91209] Solving nonlinear inequality constraints
- From: "Stuart Nettleton" <Stuart.Nettleton at uts.edu.au>
- Date: Tue, 30 Sep 2008 07:35:08 -0400 (EDT)
- Organization: University of Technology, Sydney
- References: <op.uf6ehqonidel0u@linux-vmzf.site> <48ADADDA.1010109@wolfram.com>
To finish this thread for the time being, my completed topological
processor for large scale nonlinear optimisation is provided below.
The main advantages of my approach are that the equations are retained in
an abstraction layer separate from the solving algorithm, thousands of
variables can be processed with minimal memory usage (eg 3,200 in 64Mb) ,
complex roots are detected to protect the solver from failure and help it
learn, the use of constraints to control the optimisation can minimised or
avoided altogether, and avoiding constraints allows the use of Brent's
Method with no derivatives. Of course, a constraint is needed for reasons
other than model stability. My example also includes the use of a penalty
function.
IMHO a topological solver similar to this as a Mathematica function would
be very useful to many people.
My thanks Daniel and others for their assistance and to William Nordhaus
for developing the DICE model of climate change.
I hope that this model may be of help to others.
Stuart
(* Nordhaus DICE Brief Climate Change Policy Model May 2008 *)
(* Stuart \
Nettleton September 2008 *)
starttime = AbsoluteTime[];
periods = 5; (* projection periods *)
optimpenalty =
0; (* optimisation return if iteration non-real *)
<< \
Combinatorica`
(* objective function *)
(* this program always minimises, so \
negative for maximisation *)
obj = {-cumu[periods]};
(* optimisation variables: topolology start nodes ... *)
(* .. to use \
the fast & robust Brent's Method by default, *)
(* which avoids the \
use of derivatives, make sure there *)
(* are no constraints and set \
two start variables. If possible *)
(* make sure there is one on \
either side of the zero crossing. *)
(* Constrain the range with two \
additional values. *)
(* If constraints are present, FindMinimum will \
use the Interior *)
(* Point Method and only one start variable \
should be present. *)
(* This should be your best estimate of the \
solution. While *)
(* Brent's Method is much faster than Interior \
Point, both execute *)
(* much faster and use considerable less \
memory than NMinimize. *)
(* Note that if an optimisation variable is \
set here but later *)
(* is defined as an initial variable, the \
latter is used. *)
opt = {
{\[Mu][t], 0.05, 0.8},
{k[t], 340, 5000}
};
(* exogenous parameters *)
exogparams = {
pop0 -> 6514, popg -> 0.35, popa -> 8600,
ga0 -> 0.092, dela -> 0.001,
g\[Sigma]0 -> -0.0730, 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, eland0 -> 11};
(* initial exogvars *)
exoginitial = {
gfacpop[1] -> 0, gfacpop[0] -> 0,
ga[0] -> ga0, g\[Sigma][0] -> g\[Sigma]0,
a[1] -> 0.02722, a[0] -> 0.02722,
\[Sigma][1] -> 0.13418, \[Sigma][0] -> 0.13418,
eland[0] -> eland0,
fex[0] -> fex0, fex[1] -> fex1,
\[Kappa][1] -> 0.25372, \[Kappa][0] -> 0.25372
};
(* exogenous equations *)
exogeqns = {
a[t] == a[t - 1]/(1 - ga[t - 1]),
r[t] == 1/(1 + \[Rho])^(10*(t - 1)),
eland[t] == eland0*(1 - 0.1)^(t - 1),
fex[t] == fex0 + 0.1*(fex1 - fex0)*If[t < 12, (t - 1), 0.36],
\[CapitalPi][t] == \[Kappa][t]^(1 - \[Theta]),
gfacpop[t] == (Exp[popg*(t - 1)] - 1)/Exp[popg*(t - 1)],
l[t] == pop0*(1 - gfacpop[t]) + gfacpop[t]*popa,
ga[t] == ga0*Exp[-dela*10*(t - 1)],
g\[Sigma][t] ==
g\[Sigma]0*Exp[-d\[Sigma]1*10*(t - 1) - d\[Sigma]2*10*(t - 1)^2],
\[Sigma][t] == \[Sigma][t - 1]/(1 - g\[Sigma][t]),
\[CapitalTheta][
t] == (pback*\[Sigma][t]/\[Theta])*((backrat - 1 +
Exp[-gback*(t - 1)])/backrat),
\[Kappa][t] ==
If[t >= 25, \[Kappa]21, \[Kappa]21 + (\[Kappa]2 - \[Kappa]21)*
Exp[-d\[Kappa]*(t - 2)]]
};
(* endogenous parameters *)
endogparams = {
\[Alpha] -> 2.0,
\[Gamma] -> 0.30,
\[Delta] -> 0.1,
\[Psi]1 -> 0, \[Psi]2 -> 0.0028388, \[Psi]3 -> 2,
\[Xi]1 -> 0.22, \[Xi]2 -> \[Eta]/t2xco2, \[Eta] -> 3.8,
t2xco2 -> 3, \[Xi]3 -> 0.3, \[Xi]4 -> 0.05,
\[Phi]11 -> 1 - \[Phi]12a, \[Phi]12 -> 0.189288, \[Phi]12a ->
0.189288,
\[Phi]21 -> 587.473*\[Phi]12a/1143.894, \[Phi]22 ->
1 - \[Phi]21 - \[Phi]23a,
\[Phi]23 -> 0.05, \[Phi]23a -> 0.05, \[Phi]32 ->
1143.894*\[Phi]23a/18340, \[Phi]33 -> 1 - \[Phi]32,
mat1750 -> 596.4,
\[Mu]lim -> 1,
ceindlim -> 6000,
scale1 -> 194
};
(* endogenous initial *)
endoginitial = {
y[0] -> 61.1, c[0] -> 30, inv[0] -> 31.1,
ygr[1] -> 55.667, ygr[0] -> 55.667,
k[1] -> 137, k[0] -> 137,
ceind[1] -> 0, ceind[0] -> 0,
\[CapitalLambda][1] -> 0.66203, \[CapitalLambda][0] -> 0.66203,
\[CapitalOmega][1] -> 0.99849, \[CapitalOmega][0] -> 0.99849,
mat[1] -> 808.9, mat[0] -> 808.9,
mup[1] -> 1255, mup[0] -> 1255,
mlo[1] -> 18365, mlo[0] -> 18365,
tat[1] -> 0.7307, tat[0] -> 0.7307,
tlo[1] -> 0.0068, tlo[0] -> 0.0068,
\[Mu][1] -> 0.005, \[Mu][0] -> 0.005,
cumu[1] -> 381800, cumu[0] -> 381800
};
(* endogenous variables *)
(* sn modifications of Nordhaus to render \
acyclic *)
endogeqns = {
u[t] == l[t]*(c[t]/l[t])^(1 - \[Alpha])/(1 - \[Alpha]),
k[t] == 10*inv[t] + ((1 - \[Delta])^10)*k[t - 1],
y[t] == \[CapitalOmega][t]*(1 - \[CapitalLambda][t])*ygr[t],
ygr[t] == a[t]* (k[t]^\[Gamma]) *(l[t]^(1 - \[Gamma])),
\[CapitalOmega][t] ==
1/(1 + \[Psi]1*tat[t] + \[Psi]2*(tat[t]^\[Psi]3)),
\[CapitalLambda][
t] == \[CapitalPi][t] *\[CapitalTheta][t] *\[Mu][t]^\[Theta],
c[t] == y[t] - inv[t],
eind[t] == 10 *\[Sigma][t] *(1 - \[Mu][t]) *ygr[t] ,
e[t] == eind[t] + eland[t],
(*ceind[t] == eind[t]+ceind[t-1],*)
mat[t] == e[t] + \[Phi]11*mat[t - 1] + \[Phi]21*mup[t - 1],
mup[t] == \[Phi]12*mat[t - 1] + \[Phi]22*mup[t - 1] + \[Phi]32*
mlo[t - 1],
mlo[t] == \[Phi]23*mup[t - 1] + \[Phi]33*mlo[t - 1],
for[t] == \[Eta]*Log[((mat[t] + mat[t - 1])/2)/mat1750]/Log[2] +
fex[t],
tat[t] ==
tat[t - 1] + \[Xi]1*(for[t] - \[Xi]2*
tat[t - 1] - \[Xi]3*(tat[t - 1] - tlo[t - 1])),
tlo[t] == tlo[t - 1] + \[Xi]4*(tat[t - 1] - tlo[t - 1]),
cumu[t] ==
cumu[t - 1] + (u[t]*r[t]*10)/scale1 -
Clip[\[Mu][t] - 1, {0, 1}]*10000000
(*dam[t] == ygr[t]*(1-\[CapitalOmega][t]),*)
(*s[t] == inv[t]/y[
t],*)
(*ri[t] == \[Gamma]*y[t]/k[t] -(1-(1-\[Delta])^10)/
10,*)
(*cpc[t] == c[t]*1000/l[t],*)
(*pcy[t] == y[t]*1000/l[
t],*)
};
(* endogenous constraints *)
endogcons = {(*
k[t] <= 10*inv[t]+((1-\[Delta])^10)*k[t-1],
0.02*k[periods] <= inv[periods],
100 <= k[t],
20 <= c[t],
0 <= mat[t],
100 <= mup[t],
1000 <= mlo[t],
-1 <= tlo[t]<= 20,
tat[t] <= 20,
ceind[t]<= ceindlim,
0<=q[t],
0<=inv[t],
0<=ygr[t],
0<=eind[t],
0<=\[Mu][t]<= 1*)
};
(* solve topological equations *)
toponodes[eqns_] := Module[
{eqnvars, flatvars, eqnlist, mysource, mysink, edges1, edges2,
edges3, edges, vertices2, vertices, forwardgraph, networkflows,
forwardflows, revisededges, revisedgraph, toposort,
sortedequations, sortedvertices, posfirstequation,
startvertices},
eqnvars =
Map[Cases[eqns[[#]], _Symbol[_Integer], Infinity] &,
Range[Length[eqns]]];
flatvars = Union[Flatten[eqnvars]];
eqnlist = Range[Length[eqns]];
f1[a_, b_] := {a, b};
edges1 =
Map[f1[mysource, flatvars[[#]]] &, Range[Length[flatvars]]];
edges2 =
Flatten[Map[Outer[f1, eqnvars[[#]], {eqnlist[[#]]}] &, eqnlist],
2];
edges3 = Map[f1[eqnlist[[#]], mysink] &, eqnlist];
edges = Join[edges1, edges2, edges3];
vertices2 = Join[flatvars, eqnlist];
vertices = Join[{mysource}, vertices2, {mysink}] ;
forwardgraph =
MakeGraph[vertices, (MemberQ[edges, {#1, #2}]) &,
Type -> Directed, VertexLabel -> True];
If[! AcyclicQ[forwardgraph],
Print["*** ERROR: FORWARD GRAPH IS NOT ACYCLIC SO CHECK THE \
EQUATIONS ***"]];
networkflows =
NetworkFlow[forwardgraph, 1, Length[vertices], Edge];
forwardflows =
Cases[networkflows[[All, 1, All]], {x_ /; x > 1,
y_ /; y < Length[vertices]}];
forwardedges = Map[vertices[[#]] &, forwardflows];
revisededges =
Join[Complement[edges2, forwardedges],
Map[Reverse, forwardedges]];
revisedgraph =
MakeGraph[vertices2, (MemberQ[revisededges , {#1, #2}]) &,
Type -> Directed, VertexLabel -> True];
If[! AcyclicQ[revisedgraph],
Print["*** ERROR: REVISED GRAPH IS NOT ACYCLIC SO CHECK THE \
EQUATIONS ***"]; Exit[]];
(* ShowGraph[revisedgraph] *)
toposort = TopologicalSort[revisedgraph];
sortedvertices =
Cases[vertices2[[toposort]], _Symbol[_Integer], 1];
sortedequations = Cases[vertices2[[toposort]], _Integer, 1];
posfirstequation =
Apply[Plus, First[Position[vertices2[[toposort]], _Integer, 1]]];
startvertices =
vertices2[[toposort[[Range[posfirstequation - 1]]]]];
(*startvertices =vertices2[[Select[vertices,InDegree[
revisedgraph,#]==0&]]];*)
Return[{sortedequations, sortedvertices, startvertices}]
];
(* calculate exogenous variables *)
exoginitialextended =
Cases[Union[
Flatten[Map[exoginitial /. t -> # &, Range[periods]]]] //.
exogparams /. x_Symbol[i_Integer /; i < 0] -> 0 ,
Except[False | True | Null]];
exogextended =
Cases[Union[Flatten[Map[exogeqns /. t -> # &, Range[periods]]]] //.
Join[exogparams, exoginitialextended] /.
x_Symbol[i_Integer /; i < 0] -> 0 , Except[False | True | Null]];
exogtoposolver[equations_] := Module[
{eqnorder, soleqn, solvar, outputs = {}, soltest1, soltest2},
eqnorder = toponodes[equations /. Equal -> Subtract][[1]];
For[i = 1, i <= Length[eqnorder], i++,
soleqn = equations[[eqnorder[[i]]]] //. outputs;
solvar = Cases[soleqn, _Symbol[_Integer], Infinity];
If[Length[solvar] != 0,
soltest1 =
Select[Chop[
NSolve[soleqn, solvar]], (FreeQ[solvar /. #, Complex] ) &];
If[Length[soltest1] == 0,
Print[
"*** ERROR: DURING EXOGENOUS CALCULATIONS A VARIABLE HAD NO \
SOLUTION ***"]; Exit[],
soltest2 = Select[soltest1, (solvar /. #) > 0 &];
If[Length[soltest2] == 0,
outputs = Join[outputs, First[Sort[soltest1, solvar /. # &]]],
outputs = Join[outputs, Last[Sort[soltest2, solvar /. # &]]]
];
];
];
];
outputs
];
exogaugmented =
Join[exoginitialextended, exogtoposolver[exogextended]];
Print["The exogenous variables calculate as: ", Sort[exogaugmented]];
(* calculate endogenous variables *)
interimparams = Join[exogparams, exogaugmented, endogparams];
endoginitialextended = endoginitial //. interimparams ;
allparams = Join[interimparams , endoginitialextended];
endogextended =
Cases[Union[
Flatten[Map[endogeqns /. t -> # &, Range[periods]]]] //.
allparams /. x_Symbol[i_Integer /; i < 0] -> 0 ,
Except[False | True | Null]];
endogtoponodes = toponodes[endogextended /. Equal -> Subtract];
endogeqnorder = endogtoponodes[[1]];
Print["Directed acyclic graphs and topological processing have been \
completed.... optimisation commencing..."];
(*Print["Topological order of variables:" ,endogtoponodes[[2]]];*)
Print["Please note that start vertices of the endogenous equation \
tolopogy have not been automatically included as optimisation \
variables. This is for flexibility as you may wish to use a surrogate \
based on your observation of an alternative topological sort order. \
So please check the endogenous start vertices here to confirm that \
these variables (or your surrogates) have been included with \
optimisation variables at the start: ",
endogtoponodes[[3]]
];
lenendogeqnorder = Length[endogeqnorder];
endogextendedordered = endogextended[[endogeqnorder]];
(* calculate endogenous constraints *)
endogconextended =
Cases[Union[
Flatten[Map[endogcons /. t -> # &, Range[periods]]]] //.
allparams /. x_Symbol[i_Integer /; i < 0] -> 0 ,
Except[False | True | Null]];
endogconextvars =
Union[Cases[endogconextended, _Symbol[_Integer] , Infinity]];
(* calculate objective vars *)
objvar =
Cases[Union[Flatten[Map[obj /. t -> # &, Range[periods]]]] //.
allparams /. x_Symbol[i_Integer /; i < 0] -> 0 ,
Except[False | True | Null]];
(* prepare the independent optimising variables *)
optimous =
Union[Cases[
Apply[List, Map[opt /. t -> # //. allparams &, Range[periods]]],
{_Symbol[_Integer], _Integer | _Real} |
{_Symbol[_Integer], _Integer | _Real, _Integer | _Real} |
{_Symbol[_Integer], _Integer | _Real, _Integer | _Real, _Integer \
| _Real},
Infinity]];
optimousvars =
Union[Cases[optimous //. allparams, _Symbol[_Integer] , Infinity]];
(* include any additional optimising variables arising from the \
leaves and endogenous constraints *)
variables =
Union[Join[optimous,
Partition[Complement[endogconextvars, optimousvars], 1]]];
Print["The optimising variables being used are: ", variables];
finalvars = Union[Cases[variables, _Symbol[_Integer] , Infinity]];
(* commence solve *)
(* objective function ... *)
endogoptimsolver[nmvars_] := Module[
{soleqn, solvar, outputs = {}, soltest1, soltest2},
For[i = 1, i <= lenendogeqnorder, i++,
soleqn = endogextendedordered[[i]] /. outputs;
solvar = Cases[soleqn, _Symbol[_Integer], Infinity];
If[Length[solvar] != 0,
soltest1 =
Select[Chop[
NSolve[soleqn, solvar]], (FreeQ[solvar /. #, Complex] ) &];
If[Length[soltest1] == 0,
Print["*** Warning: during optimisation ", solvar,
" became complex or null in the equation ", soleqn,
" so the specified optimisation penalty of ", optimpenalty,
" was applied ***"]; Return[optimpenalty],
soltest2 = Select[soltest1, (solvar /. #) > 0 &];
If[Length[soltest2] == 0,
outputs =
Join[outputs, First[Sort[soltest1, solvar /. # &]]],
outputs = Join[outputs, Last[Sort[soltest2, solvar /. # &]]]
];
];
];
];
Apply[Plus, objvar /. outputs]
] /; VectorQ[nmvars, NumberQ];
(* optimisation phase ... *)
(* ... use FindMinimim to optimise the \
endogenous equations .. NMinimize exhausts 16Gb of memory *)
endogoptimsolution = If[(Length[endogconextvars] == 0),
FindMinimum[endogoptimsolver[finalvars], variables],
FindMinimum[{endogoptimsolver[finalvars], endogconextended},
variables]
];
Print["The solution to the endogenous optimising variables is: ",
endogoptimsolution];
(* calculate final outputs by back substitution *)
endogoutputsolver[nmvars_] := Module[
{soleqn, solvar, outputs = nmvars, soltest1, soltest2},
For[i = 1, i <= lenendogeqnorder, i++,
soleqn = endogextendedordered[[i]] /. outputs;
solvar = Cases[soleqn, _Symbol[_Integer], Infinity];
If[Length[solvar] != 0,
soltest1 =
Select[Chop[
NSolve[soleqn, solvar]], (FreeQ[solvar /. #, Complex] ) &];
If[Length[soltest1] == 0,
Print[
"*** ERROR: DURING BACKSUBSTITUTION A VARIABLE HAD NO SOLUTION \
***"]; Exit[],
soltest2 = Select[soltest1, (solvar /. #) > 0 &];
If[Length[soltest2] == 0,
outputs = Join[outputs, First[Sort[soltest1, solvar /. # &]]],
outputs = Join[outputs, Last[Sort[soltest2, solvar /. # &]]]
];
];
];
];
outputs
];
endogaugmented =
Join[endoginitialextended,
endogoutputsolver[endogoptimsolution[[2]]]];
Print["The final outputs of the endogenous equations are: " ,
Sort[endogaugmented]];
Print["Execution time: ",
Round[N[AbsoluteTime[] - starttime]/60, 0.01], " minutes using ",
Round[N[MaxMemoryUsed[]/10^6], 0.01], " Mb; with ",
Length[finalvars], " optimising variables and ",
Length[exogaugmented] + Length[endogaugmented],
" final variables in total; ",
Length[exogparams], " exogenous parameters; ",
Length[exoginitial], " exogenous initial variables; ",
Length[exogeqns], " exogenous equations; ",
Length[exogaugmented], " final exogenous variables; ",
Length[endogparams], " endogenous parameters; ",
Length[endoginitial], " endogenous initial variables; ",
Length[endogeqns], " endogenous equations; ",
Length[endogaugmented], " final endogenous variables"
];
--
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.