qf = [] #vector containing the state to compute the objective function
tis0 = ... # containing initialization of t_1,t_2, ...

#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:-1] = 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)