# Finite dimensional optimization problem over dynamical system

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

• Could you post the actual equation for f? – Richard May 24 '19 at 15:22
• I edited my post – Smilia May 24 '19 at 15:35
• Thanks! Looking at the function, I think it might be useful if you're able to flesh out your pseudocode into actual code. That'll really lower the barrier to entry for folks who might help you. It'll also help identify whether the convergence is a code issue or a maths issue. I'm other words: a minimum, working example demonstrating the issue would be useful. – Richard May 24 '19 at 21:52