I am interested in solving numerically the following mathematical problem Consider an ode of the form $$ \dot q(t) = f(q(t),t_1,\ldots, t_N),\qquad t\in [0,T], $$ where $q\in \mathbb{R}^n$ is the state and $t_i,\, i=1,\ldots N$ are real parameters such that $0\le t_1 \le \ldots \le t_N\le T$.
Moreover $f$ is defined from the heaviside function: $H(t-t_i)$ (so I don't want to use auto differentiation)
I would like to solve: $$ \min_{t_1,\ldots,t_N}\; \lVert q(T)\rVert $$
Here is what I have done so far. I am using the scipy.optimize.minimize() function together with the 'SLSQP' method and I am integrating the ode using the 'dopri5' method from the scipy ode integrator.
Here some "pseudo code":
n = 2 #dimension of q
N = 10 # numbers of ti
qf = [] #vector containing the state to compute the objective function
tis0 = ... # containing initialization of t1,t2, ...
#init integrator
solver = ode(dyn).set_integrator('dopri5')
#function to compute the dynamcis
def dyn(t,z,n,N):
#compute dq/dt
dqdt = [0. for i in range(0,n+N)]
dqdt[...] = ... #heaviside function applied to (t-ti) appearing here
dqdt[n:-2] = 0. #dynamics for the (constants) ti
return dqdt
#compute the state q
def solout(t,q):
qf.append(q)
#function to compute the objective
def objfun(tis):
solver.set_solout(solout)
y0[n:n+N] = tis
solver.set_initial_value(y0,0.).set_f_params(n,N)
solver.integrate(T)
return norm(qf[-1])
#argument for the optimization routines : bounds, constraints...
bnds = tuple([(0.,T) for i in range(0,N)])
cons = ({'type':'ineq', 'fun': lambda x: x[0] < x[1]},...)
#call to the optimization routine
res = minimize(objfun,tis0,method='SLSQP',bounds=bnds,constraints=cons)
I am not convinced that this is a suitable approach for this problem. In fact I am observing that the algorithm is converging but the optimal solution is very close to the initialization vector tis0 ... so tis' variables don't move. Maybe I shouldn't define the tis as constants and don't put them in the dynamics ?
Also I wonder whether SLSQP is the right method to tackle this problem.
Any thoughts even not related to this solution is pretty welcome.
EDIT: the dynamics is given by: \begin{array}{l} &\dot x(t) = y(t)\, \gamma(z(t)) -x(t)\, \beta(z(t)) \\ &\dot y(t) = x(t) - y(t) + 1 \end{array} where $\gamma(z(t))=\frac{z(t)}{z(t)+1}$, $\beta(z(t))=\frac{1}{1+\gamma(z(t))}$ and $$ z(t) = \sum_{i=1}^N R_i(t_i)\,(t-t_i)\, \exp(-(t-t_i))\, H(t-t_i) $$ where $R_i(t_i) = 1+\exp(-(t_i-t_{i-1}))$ and with the conviention $t_0=-\infty$.
f
? $\endgroup$