Feature request: Time delay parameters

[chrismo-schmidt]

For our work on human control, we would be interested in introducing EOMs with an unknown time shift parameter $t_d$ that should be a free parameter for the optimization.

For example an unknown input delay, aka "reaction time":
$\dot x(t) = f(x(t)) + g(u(t-t_d))$

If one introduces this in the symbolic expressions of the EOM, the problem does not compile. It seems like this feature is currently not supported:
eom = x(t).diff - (f(t) + g(t-t_d))

I am trying workarounds with explicit time variables or an allpass. But generally it would be great if opty could support time delay out of the box.

Any other ideas for workarounds are of course also greatly appreciated.

Regards,
Christoph

[tjstienstra]

Could you perhaps write up code of a simple problem you would like to see working? That would make it easier to provide suggestions for workarounds and for implementations.

I expect that it would be something like this:

import sympy as sm
import sympy.physics.mechanics as me
from opty import Problem, create_objective_function

N = 101
dt= 0.1

t = me.dynamicsymbols._t
td = sm.Symbol("td")
x, y, u = me.dynamicsymbols("x y u")
u_delayed = u.subs({t: t - td})
eoms = sm.Matrix([
  x.diff() - u * y + u_delayed,
  y - sm.sin(x),
])

obj, obj_grad = create_objective_function(u**2, (x, y), (u,), (), N, dt)
prob = Problem(obj, obj_grad, eoms, (x, y), N, dt, instance_constraints=[x.subs({t: 0.0}) - 3])
prob.solve()

Note, it is important to choose a good representative equation of motion, e.g. x.diff() - u_delayed could just be solved by offsetting the solution of x.diff() - u.

[moorepants]

opty assumes all time varying symbols are given as x(t) and then it discretizes that to x(t_i). If you provide x(t - 5.0) as a symbol in an expression, that continuous expression needs to be replace internally with something like x(t_i - 5.0//h) where h is the constant time step and that would then map to x(t_{i-j} where j=5.0/h. Once it made it to that form, the derivatives wrt these shifted variables could then be taken as it is already done internally and things would generally work. We'd have to deal with not having values in the past before the initial time in some way.

[moorepants]

@chrismo-schmidt is your need only to shift a control (specified) trajectory relative to the output trajectories? Timo is right that we can shift things after the fact if we have the simpliest case. In a mimo system where different controls have different delays, that wouldn't be the case.

Also, you are doing parameter identification against data so, you can apply the delay shifts in your objective function or are you trying to let the optimizer also discover the delay as an unknown parameter?

[chrismo-schmidt]

Thanks for the quick reply!

Below is a toy example of what I want to achieve. Indeed my application is a system id against data and I want the optimizer to estimate an unknown input time delay as well as other parameters within the A matrix. For simplicity, I did not include the identification of any parameters in A in the example below. The only free "parameter" is the unknown delay $t_d$ from $u(t-t_d)$. The system dynamics themselves are choses arbitrarily and do not really matter to show what I want to do.

I believe just offsetting the solution as @tjstienstra suggested is not an option because this disregards that any estimated parameters of A would change with the delay (if they were included in the optimisation). Applying the delay in the objective is also not an option because it is unknown.

# -*- coding: utf-8 -*-
"""
Created on Mon Sep  9 15:52:02 2024

Toy example for input time-delay identification with sympy

@author: Christoph M. Schmidt
"""

import sympy as sm
import numpy as np

from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

import opty


#%% Simulate some data

T = 0.1 

# state space system
A = np.array([[0, 1], [-1, -1]]).T  
B = np.array([[0], [1]])  

def state_space_system(t, x):      
    
    # time-dealyed input
    u = np.sin(t-T)
    
    dxdt = A @ x[:,np.newaxis] + B * u
    return dxdt.flatten()

# initial conditions
x0 = [0, 0]  

# evaluation time
t_end = 10
steps = 1000
t_span = (0, t_end) 
t_eval = np.linspace(0, t_end, steps)

# Solve the system numerically
sol = solve_ivp(state_space_system, t_span, x0, t_eval=t_eval)

# output 
C = np.array([[0,1]])
y = C @ sol.y

# plot
fig, ax = plt.subplots(2,1)
ax[0].plot(sol.t, np.sin(sol.t), label='u(t)')
ax[0].plot(sol.t, np.sin(sol.t-T), label='u(t-T)')
ax[0].set_title("input")
ax[0].legend()
ax[1].plot(sol.t, sol.y[1], label='y(t)')
ax[1].set_title("output")
ax[1].legend()
fig.suptitle('Some time-delayed input data to a toy system')

#%% Attempt to identify time delay using sympy 

#time variables
t, td = sm.symbols("t t_d")

#input
u = sm.symbols("u", cls=sm.Function)

#states
x1, x2 = sm.symbols("x_1 x_2", cls=sm.Function)
x = sm.Matrix([[x1(t)],[x2(t)]])

#equations of motion
eoms = x.diff() - sm.Matrix(A) * x - sm.Matrix(B) * u(t-td)


#objective function and gradient
N = len(t_eval)
n = 2
q = 0

def obj(free):

    states, specified_values, constant_values = opty.parse_free(free, n, q, N) 
    
    return (states[1,:] - y)**2


def obj_grad(free):

    states, specified_values, constant_values = opty.parse_free(free, n, q, N)

    dsse = np.r_[
        np.zeros(N),
        2  * (states[1, :] - y)
    ]

    return dsse

# problem definition
prob = opty.Problem(obj, obj_grad, eoms, (x1(t), x2(t)), len(t_eval), 
                    t_end/steps, known_trajectory_map = {u(t): np.sin(t_eval)})

# solution
prob.solve()

This fails with
ValueError: dict_keys([u(t)]) are not in the provided equations of motion.

[Peter230655]

From what little I know I think that numerically integrating time delayed differential equations is highly non-trivial. I believe scipy's solve_ivp cannot do it. Would the same and more not apply to opty - or is my thinking way off?
Thanks!
The reason for my question is this:
If opty would have the same very non trivial issues like regular numerical integration, there is no point in my playing around with this. If opty, for whatever reason does not suffer from this, I'd love to play around with this issue.

[tjstienstra]

Below is a toy example of what I want to achieve.

Thank you

I believe just offsetting the solution as @tjstienstra suggested is not an option

That known_parameter map is indeed making things more complex. A problem that will remain even if timeshifts are supported is that opty would be required to extrapolate u(t) because u(0-td) with td>0 is not defined.

Another option would be to make u(t) a free variable and force it to a shape in the objective function. The disadvantage is that this makes it a, possibly stiff, multi-objective problem.
Another approach would be to define "an additional EoM" that constrains u(t) in a certain way.

[moorepants]

Maybe the instance constraints should support doing something like u(t) = u_shifted(t - 5.0). That might be the easiest to implement. Or we just have a "shift_constraint" kwarg to pass in such things.

[chrismo-schmidt]

A problem that will remain even if timeshifts are supported is that opty would be required to extrapolate u(t) because u(0-td) with td>0 is not defined.

If we require td to be bounded, we could crop all other trajectories according to the given bounds.

Another approach would be to define "an additional EoM" that constrains u(t) in a certain way.

Yesterday, I tried to develop a time-delaying filter (e.g. Pade approximation allpass, Bessel filter). That may be a feasible workaround for small delays (in the order of 10-100 ms), but delays in the orders of seconds require very high filter orders to maintain constant group delay over a certain bandwidth. So in most cases, this is probably not a suitable approach.

[Peter230655]

There is a repository on Github called ddeint. Apparently this can integrate time delay differential equations based on odeint.
I have never tried it. Maybe this can give ideas?

[moorepants]

I think we know how to solve this conceptually, but it is more about how we solve it given the design of opty. An ODE that looks like: $\dot{x} = x(t - t_d)$ maps to (with backward Euler) $\frac{x_i - x_{i-1}}{h} = x_{i - \textrm{round}(t_d/h)}$. We can create a constraint like that internally but then we need to handle the x_i where i < t_d/h period and we need to be able to compute the Jacobian with respect to these i < t_d/h variables. All that is possible, but would need some major edits in the internals to support simply adding $x(t - t_d)$ being placed directly into the equations of motion. Opty treats instance constraints separately from the EoM constraints and manually constructs the Jacobian based on knowing $i$ of each $x_i$, so adopting it to handle $x_{i - td/h}$ might be the easiest path to get opty to handle such a thing. But instance constraints are designed to be at an instant, and this is valid at all time instances, so it is a bit clunky to use instance constraints for something that should be a "path constraint" (valid at all time).