How to use scipy.integrate.ode.set_f_params() to make time dependant parameter changes?

Question

The following script illustrates some behaviour of scipy.integrate.ode.set_f_params() that confuses me.

from scipy.integrate import ode

def f(t,ys,a):
    return a

p = [1]
r = ode(f).set_initial_value(0,0).set_f_params(*p)

s = 0
while r.successful() and s<5:
    r.integrate(r.t+1)
    print r.t, r.y, p

    p[0] += 1
    r = ode(f).set_initial_value(r.y,r.t).set_f_params(*p) ### Option 1
    # r = r.set_f_params(*p) ### Option 2

    s += 1

The output using Option 1 is:

1.0 [ 1.] [1]
2.0 [ 3.] [2]
3.0 [ 6.] [3]
4.0 [10.] [4]
5.0 [ 15.] [5]

which is what I expect. I thought Option 2 should give the same output, yet it gives the following:

1.0 [ 1.] [1]
2.0 [ 2.] [2]
3.0 [ 3.] [3]
4.0 [ 5.64491239] [4]
5.0 [ 9.64491239] [5]

If someone could shed some light on this I would greatly appreciate it.

Answer 1

The reason you get different output is because in Option 2, the scipy integrator decides that it does not need to call your function 'f' very often. (In fact, in Option 2, between t = 1 and t = 3, it doesn't bother to call your 'f' even once.)

To see this, you can modify your function f like so:

def f(t,ys,a):
    print "f(%.17f) is returning %f" % (t, a)
    return a

Then Option 2 produces the output here: (Notice how the integrator, being smart, is jumping the value it probes for 't' by a factor of 10. This causes it to skip .345 all the way to 3.45. So, the fact that f is going to return different isn't noticed by the integrator until it gets up to t 4.0)

f(0.00000000000000000) is returning 1.000000
f(0.00000000000014901) is returning 1.000000
f(0.00000000000038602) is returning 1.000000
f(0.00000000000031065) is returning 1.000000
f(0.00000000310683694) is returning 1.000000
f(0.00000003417209978) is returning 1.000000
f(0.00000034482472820) is returning 1.000000
f(0.00000345135101245) is returning 1.000000
f(0.00003451661385491) is returning 1.000000
f(0.00034516924227954) is returning 1.000000
f(0.00345169552652583) is returning 1.000000
f(0.03451695836898876) is returning 1.000000
f(0.34516958679361798) is returning 1.000000
f(3.45169587103990994) is returning 1.000000
1.0 [ 1.] [1]
2.0 [ 2.] [2]
3.0 [ 3.] [3]
f(34.51695871350283085) is returning 4.000000
f(34.51695871350283085) is returning 4.000000
f(6.55822215528620234) is returning 4.000000
f(6.55822215528620234) is returning 4.000000
f(3.76234849946453931) is returning 4.000000
f(3.76234849946453931) is returning 4.000000
f(3.45169587103990994) is returning 4.000000
f(3.48276113388237274) is returning 4.000000
f(3.51382639672483554) is returning 4.000000
f(3.82447902514946492) is returning 4.000000
f(6.93100530939575776) is returning 4.000000
4.0 [ 5.64491239] [4]
5.0 [ 9.64491239] [5]

In contrast, Option 1 produces this:

f(0.00000000000000000) is returning 1.000000
f(0.00000000000014901) is returning 1.000000
f(0.00000000000038602) is returning 1.000000
f(0.00000000000031065) is returning 1.000000
f(0.00000000310683694) is returning 1.000000
f(0.00000003417209978) is returning 1.000000
f(0.00000034482472820) is returning 1.000000
f(0.00000345135101245) is returning 1.000000
f(0.00003451661385491) is returning 1.000000
f(0.00034516924227954) is returning 1.000000
f(0.00345169552652583) is returning 1.000000
f(0.03451695836898876) is returning 1.000000
f(0.34516958679361798) is returning 1.000000
f(3.45169587103990994) is returning 1.000000
1.0 [ 1.] [1]
f(1.00000000000000000) is returning 2.000000
f(1.00000004712160906) is returning 2.000000
f(1.00004853947319172) is returning 2.000000
f(1.00002426973659575) is returning 2.000000
f(1.24272163569515759) is returning 2.000000
f(3.66969529528077576) is returning 2.000000
2.0 [ 3.] [2]
f(2.00000000000000000) is returning 3.000000
f(2.00000008161702114) is returning 3.000000
f(2.00009034213922021) is returning 3.000000
f(2.00004517106961011) is returning 3.000000
f(2.45175586717085858) is returning 3.000000
f(6.96886282818334202) is returning 3.000000
3.0 [ 6.] [3]
f(3.00000000000000000) is returning 4.000000
f(3.00000009424321812) is returning 4.000000
f(3.00009707894638256) is returning 4.000000
f(3.00004853947319150) is returning 4.000000
f(3.48544327138667454) is returning 4.000000
f(8.33939059052150533) is returning 4.000000
4.0 [10.] [4]
f(4.00000000000000000) is returning 5.000000
f(4.00000010536712125) is returning 5.000000
f(4.00010264848819030) is returning 5.000000
f(4.00005132424409471) is returning 5.000000
f(4.51329376519484793) is returning 5.000000
f(9.64571817470238457) is returning 5.000000
5.0 [ 15.] [5]

Answer 2

There's usually no need to use set_f_params. In Python you can use variables from the outer scope:

def f(t, ys):
    return a

r = ode(f).set_initial_value(0,0)

a = 1
s = 0

while r.successful() and s < 5:
    r.integrate(r.t+1)
    print r.t, r.y, a
    a += 1
    s += 1

Answer 3

I have been having the same problem trying to simulate control problems using sensors that take measurements and actuators that apply forces at regular intervals. I have found a solution that works for both of our applications. The trick is to set the maximum step size to be the same as your timestep. This forces the solver to make a calculation at every timestep, meaning you can set_f_params() in the while loop and know that your changes will have an immediate effect. Heres a little example of a damped pendulum system with a controller that shows it working.

This is the best solution I have come up with, comments welcome.

from numpy import *
import matplotlib.pyplot as plt
from scipy.integrate import ode

def propagate(t, state, k, T, torque_list, integrator_t):
    integrator_t.append(t)
    torque_list.append(T)
    return array([state[1], -k*state[1] + T])


k = .5
#state is angle and angular rate
istate = array([0, 2*pi])
dt = .1

torque_list = []
integrator_t = []
solver = ode(propagate)
solver.set_integrator('dopri5',max_step = dt)
solver.set_initial_value(istate, 0)
solver.set_f_params(k, 0, torque_list, integrator_t )


newstate = []
t = []
theta_target = pi/4
while solver.successful() and solver.t < 14:
    newstate.append(solver.y)
    t.append(solver.t)

    solver.integrate(solver.t + dt)

    T = -2*(solver.y[0]-theta_target) - solver.y[1]
    solver.set_f_params(k, T, torque_list, integrator_t)

torque_list = vstack(torque_list)
integrator_t = vstack(integrator_t)

plt.figure()
plt.plot(t, newstate)
plt.title('States vs Time')
plt.xlabel('Time [s]')
plt.ylabel('Angle, Angle Rate [rad], [rad/s]')
plt.legend(['Angle','Angular Rate'])
#plt.savefig('states.png')

plt.figure()
plt.plot(integrator_t, torque_list)
plt.title('Command Torques vs Time')
plt.xlabel('Time [s]')
plt.ylabel('Command Torque [Nm]')
#plt.savefig('torques.png')
plt.show()

A plot of the systems states of a controlled pendulum. The angle is controlled to pi/4 and the angular rate is controlled to zero.

A plot of the command torques over time. Notice that they regularly update every 0.1 seconds, same as my timestep.