Another NMinimize w/ NDSolve question...
- To: mathgroup at smc.vnet.net
- Subject: [mg104828] Another NMinimize w/ NDSolve question...
- From: Joe Hays <hays.joe at gmail.com>
- Date: Wed, 11 Nov 2009 04:26:44 -0500 (EST)
All,
I'm struggling to use NDSolve with NMinimize. To start I'm using the simple
single pendulum equations of motion. I've also selected a simple cost
function to try and minimize, basically,
cost = input.inputs + sum[x(i).x(i)]
where both 'inputs' and all 'x(i)' are vectors over some time horizon. (Some
may recognize this as an optimal control style cost function.)
To do this, I am defining 'input' as a piecewise function where a value is
held constant between certain intervals. If I had 5 input updates every .2
seconds from t(i)=0 to T(f)=1 my piecewise function would look something
like:
{\[Piecewise]
myU[1][0] t<0.2
myU[1][1] t<0.4
myU[1][2] t<0.6
myU[1][3] t<0.8
myU[1][4] t<1.
0 True
}
I've placed my NDSolve and cost function calculations in a custom function,
"myFunc", that only executes if all elements of the vector holding my
optimization variables, 'ControlVars', are numeric. I understand this is
necessary to force NMinimize to use numeric values.
I've implemented all the details I could find about using NDSolve with
NMinimize however, even with this simplified example I'm unable to get
NMinimize to solve the optimization problem. I get the errors pasted in
below my code.
I've included the code for this example. I've tried to document it
reasonably well in hopes that someone is willing to wade through it with
me... (apologies in advance. you'll see from my coding style that I've still
got more to learn about Mathematica to make my code more terse.)
Thanks in advance!
================================================
*(**
*NOTE: this 'StateSpaceEOMs' data structure is a bit legacy from my old
codes. But, for this simple example we'll reuse it...*
**)*
StateSpaceEOMs = Table[0, {i, 1, 3}];
StateSpaceEOMs[[1]]={}; *(* ignore this element *)*
StateSpaceEOMs[[2]]={Y1[t], Y2[t]}; * (* these are my state vectors *)*
StateSpaceEOMs[[3]]={Y2[t], U-Sin[Y1[t]]}; *(* these are my RHS equations
*)*
myICs = {Y1[0]== .1, Y2[0]== .1}; *(* these are some dummy initial
conditions. note: Y1 = Y2 = 0 is a stable fixed point *)*
Tf = 1; *(* final time *)*
dt = .2; *(* time step for piecewise function *)*
numActuators = 1; * (* this example has a single actuator into Y1 *)*
numSimpoints = Tf/dt; *(* how many points in the piecewise input function
*)*
tmp = Table[myU[i][j], {i,1,numActuators}, {j,0,numSimpoints}]; *(* build a
list of variables representing the values of each case in the piecewise
input function *)*
ControlVars = {};
For[i=1,i<=numActuators, i++,
ControlVars = Join[ControlVars, tmp[[i]]]
]; *(* this FOR LOOP support piecewise function for multiple actuators
being concatenated, or joined, together *)*
ControlVars
Out[42]= {myU[1][0],myU[1][1],myU[1][2],myU[1][3],myU[1][4],myU[1][5]}
*(**
*This function is intending to package all the data for NDSolve, perform the
solve, and calculate an objective function using the results of NDSolve. The
objective function is then returned. One important characteristic of this
function is that it will not execute if 'ControlVars' is not a vector of
numeric values. This is supposed to force NMinimize to call the function
with numerical values during the optimization process.*
**)*
myFunc[ControlVars_, myICs_, StateSpaceEOMs_, Tf_, dt_, numActuators_:1,
qPenalty_:1, uPenalty_:1]:=Module[{s, UVec, Obj, myQ, mySolve, tmp, numvars,
Sys, dvars},
*(**
*This section creates a system of equality equations for NDSolve from my
'StateSpaceEOMs' data structure.*
**)*
numvars = Dimensions[StateSpaceEOMs[[2]]][[1]];
dvars = D[StateSpaceEOMs[[2]],t];
Sys=Table[,{i, 1, numvars}];
myICs=Table[,{i, 1, numvars}];
For[k=1, k<=numvars, k++,
Sys[[k]]=dvars[[k]]==StateSpaceEOMs[[3,k]];
];
*(**
*This section builds a piecewise function of time for each actuator using
the current guesses for 'ControlVars' NMinimize sends to 'myFunc'. *
**)*
UVec= {};
For[k=1, k<=numActuators, k++,
UVec = Join[UVec,
{Piecewise[Table[{ControlVars[[((k-1)*numSimpoints)+i]],t<i*dt},
{i,1,numSimpoints}] ]} ];
];
*(**
*Now, solve the system...*
**)*
mySolve=NDSolve[Join[Sys/.U->UVec[[1]],myICs],StateSpaceEOMs[[2]],{t,0,Tf}];
*(**
*NOTE: this is a problem... it assumes a specific var name for the control
inputs... need to generalize!!! *
**)*
*(* *
*build a sum of the inner product of all my states over time.*
**)*
myQ = 0;
For[i=1,i<=Dimensions[StateSpaceEOMs[[2]]][[1]], i++,
tmp = Table[StateSpaceEOMs[[2]][[i]]/.mySolve[[1]]/.t->j, {j,0,Tf,dt}]
myQ=myQ+ tmp.tmp;
];
*(**
*Now, define and return my objective function !!!*
**)*
Obj = qPenalty*myQ + uPenalty*ControlVars.ControlVars
] /; VectorQ[ControlVars,NumericQ] * (* <<<<------- I think this is
correct ! *)*
*(**
*Now, minimize my function given my defined ICs and system of eqs...*
**)*
NMinimize[myFunc[ControlVars, myICs, StateSpaceEOMs, Tf, dt],ControlVars]
================================================
Errors received...
During evaluation of In[45]:= Set::write: Tag Equal in Y1[0]==0.1 is
Protected. >>
During evaluation of In[45]:= Set::write: Tag Equal in Y2[0]==0.1 is
Protected. >>
During evaluation of In[45]:= Set::write: Tag Times in 0
{0.1,0.123956,0.154367,0.191243,0.230475,0.280597} is Protected. >>
During evaluation of In[45]:= General::stop: Further output of Set::write
will be suppressed during this calculation. >>
During evaluation of In[45]:= $RecursionLimit::reclim: Recursion depth of
256 exceeded. >>
During evaluation of In[45]:= $RecursionLimit::reclim: Recursion depth of
256 exceeded. >>
During evaluation of In[45]:= $RecursionLimit::reclim: Recursion depth of
256 exceeded. >>
During evaluation of In[45]:= General::stop: Further output of
$RecursionLimit::reclim will be suppressed during this calculation. >>
================================================
- Follow-Ups:
- Re: Another NMinimize w/ NDSolve question...
- From: Joe Hays <hays.joe@gmail.com>
- Re: Another NMinimize w/ NDSolve question...