Enabeling opty to estimate input time shifts.

[chris-konrad]

This is work in progress addressing #224. The code currently does not run without exceptions and added functions are untested.

Motivation

For a dynamic system $\dot{x}(t) = f(x(t)) + u(t-\tau)$, one may want to estimate the unknown input delay $\tau$. Currently, opty does not support timeshift symbols and hence can't explicitly estimate this delay (See #224).

Method

A potential way to enable opty to find the desired $\tau$ directly is outlined as follows (credits to @moorepants)

• Substitute $u(t-\tau) = u'(t)$ and consider $u'(t)$ an unknown trajectory and $\tau$ an unknown free parameter.
• Consider $u(t)$ a known trajectory.
• Create instance constraints $c_i(u'[i], \tau) = u'[i] - u[i-\frac{\tau}{\Delta t}]$ for all collocation nodes enforcing that.

Calculating the constrain Jacobian

The timeshift constraints have to be added to the constraint jacobian. For the constraints $c_i$ as defined above, the Jacobian is

• for midpoint integration:
  $G_i = [\frac{\partial c_i}{\partial u[i]}, \frac{\partial c_i}{\partial \tau}] = [1, \frac{u'[i-\frac{\tau}{\Delta t}+1] - u'[i-\frac{\tau}{\Delta t}-1]}{2\Delta t} \cdot \frac{1}{\Delta t}]$

• for backwards euler integration:
  $G_i = [\frac{\partial c_i}{\partial u[i]}, \frac{\partial c_i}{\partial \tau}] = [1, \frac{u'[i-\frac{\tau}{\Delta t}] - u'[i-\frac{\tau}{\Delta t}-1]}{\Delta t} \cdot \frac{1}{\Delta t}]$

Status

Development is ongoing in branch input-time-shifts. The current code does not run without exceptions and the added functions are untested.

TODO List

• Add a function to identify and substitute trajectories with a timeshift (16cc611, 89b7db3)
• Verify that timeshift trajectories are given in the correct format (6f79482)
• Require unshifted input to be a known trajectory (11ab07d)
• Generate instance constraints (89b7db3)
• Enable _find_closest_free_index() to create the indices of constraints with timeshift (b52a96f, 1e9a710, 11ab07d)
• Enable _instance_constraints_func() to handle timeshift constraints. (b52a96f, 1e9a710, 11ab07d)
• Enable _instance_constraints_jacobian_indices() to handle timeshift constraints (6f79482).
• Enable jacobian_indices() to find the dynamic indices of timeshift constraints
• Test if constraints_jacobian() needs to be adressed
• Make sure that existing unittests succeed
• Add bounds to timeshift parameters and process known trajectories with different lengths
• Write unittests for new functionality.
• Test new functionalitiy

[moorepants]

See chris-konrad#1 for status and plans.

This would be nice if were an issue here, because forks are often temporary or just may disappear at some point.

[chris-konrad]

This would be nice if were an issue here, because forks are often temporary or just may disappear at some point.

I added it to the first message in this PR

[chris-konrad]

I added the new indices of the new timeshift contraint partials to the jacobian indices. This assumes the following layout of the Jacobian. Can you confirm that this is correct, @moorepants? Green is the existing Jacobian, yellow is the extension with timeshift constraints.

The constraint definition is as above:
$G_i = [\frac{\partial c_i}{\partial u[i]}, \frac{\partial c_i}{\partial \tau}]$

Creating a new problem with Problem() now runs through without exceptions. I did not test it any further, so far.

[Peter230655]

@chris-konrad : Curiosity: Once this works will opty be able to handle time delayed eoms, like $\dfrac{dx(t)}{dt} = f(x(t - \tau), u(x - \tau), parameters)$ ?
Thanks!

[chris-konrad]

@Peter230655: No, unfortunately this will not directly be possible. The implemented approach requires that the timeshifted trajectory is a known trajectory.

[Peter230655]

@Peter230655: No, unfortunately this will not directly be possible. The implemented approach requires that the timeshifted trajectory is a known trajectory.

Thanks!
Still good for me: I am mostly a user of sympy mechanics and of opty, playing around with them to pass time in retirement.
This, when implemented will enlarge the playing field for me. :-)

[chris-konrad]

All previously existing unittests now run fine but it seems to break some examples with variable times. This is little surprising as I did not pay any attention to that during implementation. Before addressing these issues, I will focus on testing if the new functions work at all.

[moorepants]

This assumes the following layout of the Jacobian. Can you confirm that this is correct, @moorepants?

Are you treating $\tau$ the same as any unknown parameter? i.e. is it in the list of unknown_parameters? This would have some bearing on the construction of the jacobian (order of rows).

[moorepants]

Why would u_shift not be added here instead of u(t-t_d)? Shouldn't the eom's not have t_d present?

[chris-konrad]

This is the eom that I pass to the constructor of Problem. It has t_d present in u(t-t_d). The ConstraintCollocator then looks through this eom, finds all expressions of the form U(time_symbol +/- symbol) and substitutes them with U_shift(time_symbol (see: ConstraintCollocator._substitute_timeshift_trajectories(self))

[moorepants]

Why did you choose adding a new attribute here, instead of passing u_shift - u(t - t_d) directly to instance constraints?

[chris-konrad]

After substituting the timeshift trajectories, ConstraintCollocator automatically appends instance constraints representing u_shift(t) - u(t - t_d) to the other instance constraints (see: ConstraintCollocator._generate_timeshift_constraints()) , I just didn't write a test for this yet. The attribute ConstraintCollocator.timeshift_trajectory_substitutes, that I check in this test, is just a map from the substitutes to the original trajectory and the timeshift parameter that I occasionally need throughout the code to look up which functions and parameters relate to the timeshift.

[chris-konrad]

This assumes the following layout of the Jacobian. Can you confirm that this is correct, @moorepants?

Are you treating τ the same as any unknown parameter? i.e. is it in the list of unknown_parameters? This would have some bearing on the construction of the jacobian (order of rows).

Thats what I do. Indeed there seems to be a problem in the Jacobian index function related to this. I will have a look.

[moorepants]

Should be safe to call squeeze() without the if statement, as it will only affect arrays with an empty dimension.

[Peter230655]

I may have found a simple hack to do this:
In the time delayed driving force $u(t - \tau)$ I replaced t with a variable T, and in the eom I set $\dfrac{d}{dt} T(t) = 1.0$, so T has always the same value as t. (with instanc_constraint $T(t_0) = 0$.
In my simple example, a pendulum where u is the torque applied to it, $\tau$ is found easily, if I give reasonable bounds for it.
I set u(t) = 0.0 for t $\leq 0.0.$
@moorepants: If my initial results prove right, would this be an example for examples_gallery/beginner ?

[moorepants]

Yes, this sounds like a useful example to have.

[moorepants]

Are you doing u(T(t) - tau)?

[Peter230655]

Are you doing u(T(t) - tau)?

Yes, this is exactly what I do. So far looks o.k.

[moorepants]

I don't think the internals of opty support anything but t in the function Function('u')(t).

[Peter230655]

Of course only you know this best! I did not use a 'general function'. I did not think that given a bunch of measurements, opty could determine whether the driving for was a sine function or an exponential function, or whatever.
I used $sm.sin(T(t) - \tau)$. My simple thinking was this: opty does accept functions of state variables and treats them correctly.
So if I introduce a state valiable with the same numerical values as t it might work.
So far looks o.k., but I will play around more.
Actually I used torque = $F \cdot \sin(\omega \cdot (T(t) - \tau)$ and let it estimate F, $\omega, \tau$.

[moorepants]

s i n ( T ( t ) − τ ) will work, but not if the function is a state variable, e.g. x(T(t) - t) or input variable r(T(t) - t). I think all of the instance constraint code only works with x(t) or r(t) and should break (or give incorrect results) otherwise.

[moorepants]

Feel free to prove me wrong though!

[Peter230655]

I do not think, you are wrong! :-)
I also do not think, a state variable will work with anything but t, but this is not what I need - if I understood Chris' problem correctly.

Anyway, let me play around a bit more, hopefully I can push it tomorrow, then you will see if useful or not so good.

[moorepants]

I also do not think, a state variable will work with anything but t, but this is not what I need - if I understood Chris' problem correctly.

Christoph precisely needs an unknown input to have a delay, so r(t - tau) instead of r(t).

[Peter230655]

Then I missunderstood, was too easy to be true! 😊😊
So is opty expected to estimate the delay and what function r is?

[Peter230655]

@chris-konrad
I am still confused how this works. If I understood Jason correctly, you want opty to find r and $\tau$ in $r(t - \tau)$
But a function $\hat r(t) := r(( t - \tau)$ would be indistinguishable from r (?). How can opty find r and $\tau$ instead of $\hat r$