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>