Re: Re: Solving nonlinear inequality constraints
- To: mathgroup at smc.vnet.net
- Subject: [mg91405] Re: [mg91243] Re: [mg91209] Solving nonlinear inequality constraints
- From: "Stuart Nettleton" <Stuart.Nettleton at uts.edu.au>
- Date: Wed, 20 Aug 2008 06:22:17 -0400 (EDT)
- Organization: University of Technology, Sydney
Daniel's suggestions have led me to refine my notebook (below), which has
reduced memory usage slightly and thereby allowed a marginally increased
projection period from 9 to 10 periods:
1. Introducing a function for mu[t_] - this has not been achieved. Perhaps
I am missing something very major but I cannot see how this can be done!
As soon as mu[1], mu[2] are generated in variables, they are immediately
evaluated by the arbitrary function (such as mu[t_]:=0.5). The effect of
providing functions for all the variables is to evaluate the objective and
constraints. Therefore, NMinimize merely has a numeric value for the
objective function, a set of true and false for the constraints, and no
variables for optimising. The notebook below shows that not only mu[t] but
also variables like k[t] and s[t] need to be "freed" for optimization (as
in the standard Nordhaus model), which I can only see means that variables
cannot be given their own functions.
2. Removing Hold - the need for Hold is removed by using _?NumberQ in
functions to prevent the symbolic evaluation that would otherwise occur.
The first notebook below has this feature. The second notebook continues
to use Holds, which are eliminated before the NMinimize function. I have
some preference for the Hold approach because answers with small numbers
of periods can be different than when _?NumberQ is used and I tend to
think symbolic answers are more accurate.
3. Using the Hold approach (ie. without using _?NumberQ in functions)
provides the following performance and memory usage on an intel 64bit/16Gb:
Projection Compilation Total Memory
Periods Time Minutes Mb
5 1.0 sec 1.2 234
6 2.4 secs 2.4 535
7 8.9 secs 24.8** 1234
8 0.5 mins 12.8 2572
9 2.3 mins 40.9 4943
10 9.4 mins 66.3 8908
11 33.5 mins* runs out of memory
* 10.5Gb used in compilation phase
** looks odd but retested and this was the same result, so perhaps its due
to the shape of the curve
Usually at 6 periods and upwards there were one or two complaints that
"NMinimize::incst: NMinimize was unable to generate any initial points \
satisfying the inequality constraints" but the function kept running to a
final solution.
4. As previously observed, the result using _?NumberQ in functions is
surprisingly similar to that with Hold. There is a shorter compilation
phase of 1 minute offset by longer optimisation phase of 1.5 minutes,
resulting in a slightly longer calculation period. The comparative
statistics for 10 periods are:
Approach Compilation Total Memory
Minutes Minutes Mb
Hold 9.4 mins 66.3 8908
_?NumberQ 8.3 mins 66.8 8901
Conclusion:
It is some success that the number of projection periods has increased
from 8 to 10.
However, I still need to ramp-up to 60 periods and would like to achive
this in less than an hour if possible.
The continued interest of Daniel and participants on this list is warmly
appreciated!
Thanks,
Stuart
**** Notebook using _?NumberQ ****
(*AppendTo[$Echo,"stdout"];*)
(*SetOptions["stdout",PageWidth->110] ; \
*)
(* Nordhaus Brief Climate Change Policy Model May 2008 *)
(* \
Stuart Nettleton July 2008 *)
starttime = AbsoluteTime[];
periods = 5; (* projection periods *)
maxit = 5000; (* maximum \
iterations *)
(* objective function *)
(* program always minimises, so negative for \
maximisation *)
obj = {-cumu[periods]};
(* 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 = {{\[Mu][t],0.001,1},{k[t],100,10000}};*)
opt = {{\[Mu][t], 0.001, 1}};
(* exogenous parameters *)
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;
(* exogenous variables *)
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)];
a[t_] := a[t] = a[t - 1]/(1 - ga[t - 1]);
g\[Sigma][t_] :=
g\[Sigma][t] =
g\[Sigma][0]*
Exp[-d\[Sigma]1*10*(t - 1) - d\[Sigma]2*10*(t - 1)^2];
\[Sigma][t_] := \[Sigma][t] = \[Sigma][t - 1]/(1 - g\[Sigma][t]);
\[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);
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)]];
\[CapitalPi][t_] := \[CapitalPi][t] = \[Kappa][t]^(1 - \[Theta]);
(*s[t_]:= s[t] =sr;*)
(* initial values of exogenous variables *)
ga[0] = 0.092;
a[1] = a[0] = 0.02722;
g\[Sigma][0] = -0.0730;
\[Sigma][1] = \[Sigma][0] = 0.13418;
eland[0] = 11;
\[Kappa][1] = \[Kappa][0] = 0.25372;
(* endogenous parameters *)
\[Alpha] = 2.0;
\[Gamma] = .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;
c0 = 30;
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_?NumberQ] := eind[t - 1] + ceind[t - 1];
eind[t_?NumberQ] :=
10 *\[Sigma][t] *(1 - \[Mu][t]) *ygr[t] + eland[t];
for[t_?NumberQ] := \[Eta]*(Log[(matav[t] + 0.000001)/mat1750]/
Log[2]) + fex[t];
mat[t_?NumberQ] :=
eind[t - 1] + \[Phi]11*mat[t - 1] + \[Phi]21*mup[t - 1];
matav[t_?NumberQ] := (mat[t] + mat[t + 1])/2;
mlo[t_?NumberQ] := \[Phi]23*mup[t - 1] + \[Phi]33*mlo[t - 1];
mup[t_?NumberQ] := \[Phi]12*mat[t - 1] + \[Phi]22*
mup[t - 1] + \[Phi]32*mlo[t - 1];
tat[t_?NumberQ] :=
tat[t - 1] + \[Xi]1*(for[t] - \[Xi]2*
tat[t - 1] - \[Xi]3*(tat[t - 1] - tlo[t - 1]));
tlo[t_?NumberQ] := tlo[t - 1] + \[Xi]4*(tat[t - 1] - tlo[t - 1]);
ygr[t_?NumberQ] := a[t]* k[t]^\[Gamma] *l[t]^(1 - \[Gamma]);
dam[t_?NumberQ] :=
ygr[t]*(1 - 1/(1 + \[Psi]1*tat[t] + \[Psi]2*(tat[t]^\[Psi]3)));
\[CapitalLambda][t_?NumberQ] :=
ygr[t] *\[CapitalPi][t] *\[CapitalTheta][t] *\[Mu][t]^\[Theta];
y[t_?NumberQ] :=
ygr[t]*(1 - \[CapitalPi][t]*\[CapitalTheta][
t]*\[Mu][t]^\[Theta])/(1 + \[Psi]1*
tat[t] + \[Psi]2*(tat[t]^\[Psi]3));
inv[t_?NumberQ] := (y[t] + 0.001)*s[t];
(*k[t_?NumberQ]:=10*inv[t-1]+((1-\[Delta])^10)*k[t-1];*)
ri[t_?NumberQ] := \[Gamma]*y[t]/k[t] - (1 - (1 - \[Delta])^10)/10;
c[t_?NumberQ] := y[t] - inv[t];
u[t_?NumberQ] := ((c[t]/l[t])^(1 - \[Alpha]) - 1)/(1 - \[Alpha]);
cumu[t_?NumberQ] := cumu[t - 1] + (l[t]*u[t]*r[t]*10)/scale1;
cpc[t_?NumberQ] := c[t]*1000/l[t];
pcy[t_?NumberQ] := y[t]*1000/l[t];
(* initial values of endogenous variables *)
cumu[1] = cumu[0] = cumu0;
ceind[1] = ceind[0] = ceind0;
mat[1] = mat[0] = mat0;
mlo[1] = mlo[0] = mlo0;
mup[1] = mup[0] = mup0;
tat[1] = tat[0] = tat0;
tlo[1] = tlo[0] = tlo0;
y[0] = y0;
k[1] = k[0] = k0;
c[0] = c0;
\[Mu][1] = \[Mu][0] = \[Mu]0;
(* endogenous inequality constraints*)
endog = {
k[t] <= 10*inv[t - 1] + ((1 - \[Delta])^10)*k[t - 1],
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,
ceind[t] <= 6000,
0 <= y[t],
0 <= inv[t],
0 <= ygr[t],
0 <= eind[t],
0 <= matav[t],
0 <= \[Mu][t]
};
(* processing ... *)
(*prepare the objective function*)
objvars = Simplify[obj];
objectimous =
Simplify[
Flatten[Union[Map[objvars /. t -> # &, Range[periods]]] /.
x_Symbol[i_Integer /; i < 0] -> 0]];
(*prepare the endogenous constraints*)
endovars = Simplify[endog];
endogenous =
Simplify[
Flatten[Union[Map[endovars /. t -> # &, Range[periods]]] /.
x_Symbol[i_Integer /; i < 0] -> 0]];
(*prepare the independent optimising variables*)
optvars = Simplify[opt];
optimous =
Union[Cases[
Flatten[Map[optvars /. t -> # &, Range[periods]],
1], {x_Symbol[_Integer], _Integer | _Real} | \
{x_Symbol[_Integer], _Integer | _Real, _Integer | _Real} | \
{x_Symbol[_Integer], _Integer | _Real, _Integer | _Real, _Integer | \
_Real}, Infinity]];
optimousvars =
Union[Cases[optimous, x_Symbol[i_Integer] -> x[i], Infinity]];
(*include any additional optimising variables arising from the \
endogenous constraints*)
primafacievars =
Union[Cases[objectimous, x_Symbol[i_Integer] -> x[i], Infinity],
Cases[endogenous, x_Symbol[i_Integer] -> x[i], Infinity]];
extendedvariables =
Union[optimous, Complement[primafacievars, optimousvars]];
(*prepare optimising function*)
Print[Round[(AbsoluteTime[] - starttime), 0.01],
" seconds compilation time, commencing optimisation ..."];
soln = NMinimize[Join[objectimous, endogenous], extendedvariables,
MaxIterations -> maxit];
(*output results*)
Print["Solution: ", soln];
Print[Round[(AbsoluteTime[] - starttime)/60, 0.01],
" minutes from start to completion"];
Print[Round[N[MaxMemoryUsed[]*10^-6], 0.1], " Mb memory used"];
**** Notebook using Hold ****
(*AppendTo[$Echo,"stdout"];*)
(*SetOptions["stdout",PageWidth->110] ; \
*)
(* Nordhaus Brief Climate Change Policy Model May 2008 *)
(* \
Stuart Nettleton July 2008 *)
starttime = AbsoluteTime[];
periods = 11; (* projection periods *)
maxit = 5000; (* maximum \
iterations *)
(* objective function *)
(* program always minimises, so negative for \
maximisation *)
obj = Hold[-cumu[periods]];
objvars = Apply[List, Map[Hold, obj]];
(* preferences for 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 *)
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;
(* exogenous variables *)
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)];
a[t_] := a[t] = a[t - 1]/(1 - ga[t - 1]);
g\[Sigma][t_] :=
g\[Sigma][t] =
g\[Sigma][0]*
Exp[-d\[Sigma]1*10*(t - 1) - d\[Sigma]2*10*(t - 1)^2];
\[Sigma][t_] := \[Sigma][t] = \[Sigma][t - 1]/(1 - g\[Sigma][t]);
\[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);
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)]];
\[CapitalPi][t_] := \[CapitalPi][t] = \[Kappa][t]^(1 - \[Theta]);
(*s[t_]:= s[t] =sr;*)
(* initial values of exogenous variables *)
ga[0] = 0.092;
a[1] = a[0] = 0.02722;
g\[Sigma][0] = -0.0730;
\[Sigma][1] = \[Sigma][0] = 0.13418;
eland[0] = 11;
\[Kappa][1] = \[Kappa][0] = 0.25372;
(* endogenous parameters *)
\[Alpha] = 2.0;
\[Gamma] = .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;
c0 = 30;
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];
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];
matav[t_] := (mat[t] + mat[t + 1])/2;
mlo[t_] := \[Phi]23*mup[t - 1] + \[Phi]33*mlo[t - 1];
mup[t_] := \[Phi]12*mat[t - 1] + \[Phi]22*mup[t - 1] + \[Phi]32*
mlo[t - 1];
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]);
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));
inv[t_] := (y[t] + 0.001)*s[t];
(*k[t_]:=10*inv[t-1]+((1-\[Delta])^10)*k[t-1];*)
ri[t_] := \[Gamma]*y[t]/k[t] - (1 - (1 - \[Delta])^10)/10;
c[t_] := y[t] - inv[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;
cpc[t_] := c[t]*1000/l[t];
pcy[t_] := y[t]*1000/l[t];
(* initial values of endogenous variables *)
cumu[1] = cumu[0] = cumu0;
ceind[1] = ceind[0] = ceind0;
mat[1] = mat[0] = mat0;
mlo[1] = mlo[0] = mlo0;
mup[1] = mup[0] = mup0;
tat[1] = tat[0] = tat0;
tlo[1] = tlo[0] = tlo0;
y[0] = y0;
k[1] = k[0] = k0;
c[0] = c0;
\[Mu][1] = \[Mu][0] = \[Mu]0;
(* endogenous inequality constraints*)
endog = Hold[
k[t] <= 10*inv[t - 1] + ((1 - \[Delta])^10)*k[t - 1],
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,
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 ... *)
(* prepare the objective function *)
objectimous =
Simplify[
Flatten[Apply[List,
Map[ReleaseHold,
Union[Map[objvars /. t -> # &, Range[periods ]]]]] /.
x_Symbol[i_Integer /; i < 0] -> 0]];
(* prepare the endogenous constraints *)
endogenous =
Simplify[
Flatten[Apply[List,
Map[ReleaseHold,
Union[Map[endovars /. t -> # &, Range[periods ]]]]] /.
x_Symbol[i_Integer /; i < 0] -> 0]];
(* prepare the independent optimising variables *)
optimous =
Union[Cases[
Apply[List,
Map[ReleaseHold,
Flatten[Map[optvars /. t -> # &, Range[periods]]]]],
{x_Symbol[_Integer], _Integer | _Real} |
{x_Symbol[_Integer], _Integer | _Real, _Integer | _Real} |
{x_Symbol[_Integer], _Integer | _Real, _Integer | _Real, \
_Integer | _Real},
Infinity]];
optimousvars =
Union[Cases[optimous, x_Symbol[i_Integer] -> x[i], Infinity]];
(* include any additional optimising variables arising from the \
endogenous constraints *)
primafacievars =
Union[Cases[objectimous, x_Symbol[i_Integer] -> x[i], Infinity],
Cases[endogenous, x_Symbol[i_Integer] -> x[i], Infinity]] ;
extendedvariables =
Union[optimous, Complement[primafacievars, optimousvars]];
(* prepare optimising function *)
Print[Round[(AbsoluteTime[] - starttime), 0.01],
" seconds compilation time, commencing optimisation ..."];
soln = NMinimize[Join[objectimous, endogenous], extendedvariables,
MaxIterations -> maxit];
(* output results *)
Print["Solution: ", soln];
Print[Round[(AbsoluteTime[] - starttime)/60, 0.01],
" minutes from start to completion"];
Print[Round[N[MaxMemoryUsed[]*10^-6], 0.1], " Mb memory used"];
--
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.