Finite dimensional optimization problem over dynamical system

Question

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$.