Question: Option remember

Hi,

I am trying to implement some numerical algorithm for an optimal control problem (still).

It involves double cycle. In the inner one two DE systems are solved sequentially and the result is used to update values of some objective function.

This part uses remeber option thanks to pagan  and is ok.

However, in the outer cycle I use the values of this objective function to obtain solution for yet another DE and update ICs for the inner cycle.

In this part option remeber seems of no help and I have no idea as to how deal with it.

Here is the code:

restart:
> read "BaU":
> sys1:=STATC:
> sys2:=COSTATC:
>
> systot:=STATG;
> GRAD:=subs(u=FF(t),GRADNUM):
> RECC:=FF(t)+alpha*GRAD:
> REC:=piecewise(RECC>1,1,RECC<1,RECC,RECC=1,1):
> F[0]:=t->0.1;
> UU[0]:=t->0.1:
> G[0]:=t->0.1:
> GR[0]:=0:
> N:=5:
> alpha:=0.005:
> TAUS[0]:=0.7307:
> KS[0]:=733.2:
> M[0]:=808.9:x:=10:
> KS[0];
> kk[0]:=t->733.2:

#the outer cycle starts here

for j from x to 220 by x do
> sys1:=subs(K[0]=KS[j-x],tau[s]=TAUS[j-x],STATC):
> sys2:=subs(K[0]=KS[j-x],tau[s]=TAUS[j-x],COSTATC):
> GRAD:=subs(K[0]=KS[j-x],tau[s]=TAUS[j-x],subs(u=FF(t),GRADNUM)):
> RECC:=FF(t)+alpha*GRAD:
> REC:=piecewise(RECC>1,1,RECC<1,RECC,RECC=1,1):
> print(j);
> F[0][j]:=UU[j-x]:
> #print(F[0][j]);

#the inner cycle starts here

> for i from 1 to N do
>
> print(i);
> newsys1:=subs(u(t)=F[i-1][j](t),sys1):
> #print(newsys1);
>
> #print(F[i-1]);

#first DE system computation

>  sol1[i][j]:=dsolve(newsys1, numeric,range=j-x+0.1..j,output=listprocedure, known=[F[i-1][j]]):
> #print(i);
> KK[i][j]:=subsop(3=remember,eval(K(t),sol1[i][j])):
>
> #KKS[i]:=piecewise(t<j-x,kk[j-x](t),t>j-x,KK[i](t),t=j-x,kk[j-x](t)):
> #print(KKS[i]);

#second DE system computation

> newsys2:=subs(u=F[i-1][j],K(t)=KK[i][j](t),sys2):
>
> #print(newsys2);
> sol2[i][j]:=dsolve(newsys2, numeric, output=listprocedure, range=j-x+0.1..j,known=[F[i-1][j],KK[i][j]]):
>
> PSI1[i][j]:=subsop(3=remember,eval(psi[1](t),sol2[i][j])):

#update of an objective function

> F[i][j]:=unapply('evalf'(subs(psi[1](t)=subsop(3=remember,eval(psi[1](t),sol2[i][j]))(t),K(t)=subsop(3=remember,eval(K(t),sol1[i][j]))(t),
> FF=F[i-1][j],REC)),t,numeric,proc_options=[remember]):
> #print(F[i]);
> GR[i][j]:=unapply('evalf'(subs(psi[1](t)=subsop(3=remember,eval(psi[1](t),sol2[i][j]))(t),K(t)=subsop(3=remember,eval(K(t),sol1[i][j]))(t),FF=F[i-1][j],                              GRAD)),t,numeric,proc_options=[remember]):
> end do:

#end of inner cycle

#outer DE system

aftersys:=subs(k[0]=KS[j-x],tau[0]=TAUS[j-x],m[0]=M[j-x],u(t)=F[N][j](t),systot);
>
> #print(aftersys);
> aftersol[j]:=dsolve(aftersys, numeric,range=j-x+0.1..j,output=listprocedure, known=[F[N][j]]):
> kk[j]:=subsop(3=remember,eval(k(t),aftersol[j])):
> mm[j]:=subsop(3=remember,eval(m(t),aftersol[j])):
> tt[j]:=subsop(3=remember,eval(tau(t),aftersol[j])):

#updates of ICs for the inner DE systems newsys1 and newsys2

> KS[j]:=evalf(kk[j](j)):
> TAUS[j]:=evalf(tt[j](j)):
> M[j]:=evalf(mm[j](j)):
> UU[j]:=unapply('evalf'(subs(psi[1](t)=subsop(3=remember,eval(psi[1](t),sol2[N][j]))(t),K(t)=subsop(3=remember,eval(K(t),sol1[N][j]))(t),
> FF=F[N-1][j],REC)),t,numeric,proc_options=[remember]):
>
> print(UU[j](1));
> end do:

BaU is just the file with exact formulae for DEs and functions I use.

I have two problems.

First, the algorithm takes around 2 hours to compute until j=120 with step x=10, which is long. I think there might be some shortcut with remeber trick, as I use in the inner cycle (with i steps), but I do not know how to make it.

Second, when I evaluate variables from the output of sol1[i][j], it evaluates variable not in the range j-x..j as desired, but  on the whole available range, given by boundary conditions. Can I push Maple to evaluate variables only in some restricted (smaller) range. This may boost performance, I think...

Please, any ideas might be useful

AB