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):

#function to compute the objective
def objfun(tis):
    y0[n:n+N] = tis
    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$.

  • 1
    $\begingroup$ Could you post the actual equation for f? $\endgroup$
    – Richard
    May 24, 2019 at 15:22
  • $\begingroup$ I edited my post $\endgroup$
    – Smilia
    May 24, 2019 at 15:35
  • 1
    $\begingroup$ 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. $\endgroup$
    – Richard
    May 24, 2019 at 21:52


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.