Re: Another NMinimize w/ NDSolve question...
- To: mathgroup at smc.vnet.net
- Subject: [mg104865] Re: [mg104828] Another NMinimize w/ NDSolve question...
- From: Joe Hays <hays.joe at gmail.com>
- Date: Thu, 12 Nov 2009 06:02:50 -0500 (EST)
- References: <200911110926.EAA29286@smc.vnet.net>
All,
After a nights rest and revisiting this problem fresh this morning... I
believe I've now successfully fixed my little test example. I'm including it
here for future users of NDSolve and NMinimize hoping it might save them
time...
%%%%%%%%%%%%%%%%%%%
START OF CODE
%%%%%%%%%%%%%%%%%%%
In[1]:= *(**
*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 *)
In[8]:= 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,1,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 supports piecewise functions for multiple actuators
being concatenated, or joined, together *)
In[11]:= *(**
*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[{UVec, Obj, myQ, mySolve, tmp, numvars,
Sys, dvars, k, j, i, numSimpoints},
*(**
*This section creates a system of equality equations for NDSolve from my
'StateSpaceEOMs' data structure.*
**)*
numvars = Dimensions[StateSpaceEOMs[[2]]][[1]];
numSimpoints = Tf/dt;
dvars = D[StateSpaceEOMs[[2]],t];
Sys=Table[0,{i, 1, numvars}];
(*myICs=Table[0,{i, 1, numvars}];*)
For[k=1, k<=numvars, k++,
Sys[[k]]=dvars[[k]]==StateSpaceEOMs[[3,k]];
(*myICs[[k]]=StateSpaceEOMs[[2,k]]==ICs[[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,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 ! *)
In[13]:= *(**
*Now, minimize my function given my defined ICs and system of eqs...*
**)*
NMinimize[myFunc[ControlVars, myICs, StateSpaceEOMs, Tf, dt],ControlVars]
Out[13]=
{0.111733,{myU[1][1]->-0.0693213,myU[1][2]->-0.0354595,myU[1][3]->-0.0112234,myU[1][4]->0.00268734,myU[1][5]->0.0062375}}
%%%%%%%%%%%%%%%%%%%
END OF CODE
%%%%%%%%%%%%%%%%%%%
On Wed, Nov 11, 2009 at 4:26 AM, Joe Hays <hays.joe at gmail.com> wrote:
> 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. >>
>
> ================================================
>
>
>
- References:
- Another NMinimize w/ NDSolve question...
- From: Joe Hays <hays.joe@gmail.com>
- Another NMinimize w/ NDSolve question...