GeezerVison

Coverage decreased (-0.008%) to 88.791% when pulling dd681e2 on sawyerbfuller:dynamics into 9c6c619 on python-control:master.

Following #83 I don't see why anyone won't write

def axplusbu(x, u): return A @ x + B @ u

Instead what you might want to do is to provide a callable so that people can plug it in wherever they want with all bells and whistles included. In english, you return a function which is built out of the parts of the IOsystem and that is independent from the model itself. Otherwise this doesn't seem like doing anything too useful.

What happens if you have a LinearIOSystem, which is both an InputOutputSystem and a StateSpace system? Should the call signature be sys.dynamics(t, x, u) or sys.dynamics(x, u)?

Also, what happens if a system has parameters (allowed for I/O systems)? Should I be able to set those parameter values when I call the (public) dynamics function? If not, how do I change parameter values?

@ilayn I agree, it's not too onerous to write that out. I added these because I have been playing around with such an integrator and it does seem to clean up/improve code appearance to be able to write x[k+1] = P.dynamics(x[k], u[k]) than to write out all of the matrices (in my code, I'm working with a number of systems at once and spelling out all of their respective A, B, C, D matrices becomes unwieldy). You can also use it in scipy.solve_ivp with something like scipy.solve_ivp(lambda t, x: P.dynamics(x, -K@x), Tfinal) (Note: I might have the syntax on that slightly wrong).

can you elaborate more on what you had in mind?

@murrayrm

What happens if you have a LinearIOSystem, which is both an InputOutputSystem and a StateSpace system? Should the call signature be sys.dynamics(t, x, u) or sys.dynamics(x, u)? How would what you're proposing work?

I see 2 options:

1. Make both compatible consistent with scipy's solve_ivp and put the t argument first.
2. put the t argument at the end as an optional keyword argument and possibly raise a warning if it's not None for LTI systems (can LinearIOsystems be time-varying?).

In iosys, for numerical integration using solve_ivp, it creates a separate function ivp_dynamics, partly because solve_ivp doesn't have accomodation for a u. So there's no particular need to have dynamics keep the t argument first. So I think I favor option 2 given the relative importance of t vs. x and u arguments.

Also, what happens if a system has parameters (allowed for I/O systems)? Should I be able to set those parameter values when I call the (public) dynamics function? If not, how do I change parameter values?

NOt sure how params are used in IOSystems, but it looks like the old _rhs method doesn't accept them as currently written. Seems like it should be possible to revise it in the future to accomodate those parameters.
For LTI systems, we could also allow a params dict as a keyword argument for consistency's sake, but raise a warning if it's not empty for LTI systems? (which are immutable).

No matter what, I think we want the InputOutputSystem and StateSpace signatures for dynamics and output to be the same. Otherwise it will be very confusing.

If we leave in t, I think it should be listed first:

• Consistent with SciPy but also with MATLAB (eg, ode45).
• Consistent with output arguments (t, y = step_response(sys)) and input arguments (t, y = initial_response(sys, t, x0), t, y = forced_response(sys, t, u, x0).

Another possibility would be to define dynamics and output for LTI systems and then define these for InputOutputSystems with a signature of the form sys.dynamics(x, u[, t]) where t = 0 by default. This would then give the same signature for both and would make t optional (which is certainly going to be the most common case).

However, I would leave the signatures for _rhs and _out unchanged, so that when when you define a nonlinear system that the argument order would be t, x, u (consistent with MATLAB and SciPy [minus the u argument]). So for

sys = NonlinearIOSystem(updfcn, outfcn)

the signatures of updfcn and outfcn are still f(t, x, u). If we change this order, then we break existing user code and are incompatible with other packages. Of course, we could just redefine things internally, but from a developer point of view it is might be tedious to try to remember to flip things around all the time.

I'll need to look through the params processing. It should be the case that is optional but if specified then it overrides the local parameter dictionary stored in the instance. See examples/steering.ipynb to see examples of how it is used.

FYI: for the iosys module, LinearIOSystem already issues a warning if you pass parameters. We could do the same for state space.

OK good point, it looks like there's a pretty strong convention of leaving t first in the library and elsewhere so I'm in favor of sticking with that, I think it's not too bad.

The latest commit changes Statesp.dynamics and output to have t first, and to ignore it.

With this, the convention now is that the call signature is dynamics(t, x, u) and output(t, x, u) everywhere, including what used to be _rhs and _out. Perhaps remaining to be done, in a future commit, is to allow for a params argument.

You can also use it in scipy.solve_ivp with something like scipy.solve_ivp(lambda t, x: P.dynamics(x, -K@x), Tfinal)

Yes instead of doing that you can actually call the dynamics or the output methods and obtain a function so it becomes

.....
    def dynamics(self):
        def f(x, u):
           return self.a.copy() @ x + self.b.copy() @ u
        return f
....

And then you can ask for that function and use it elsewhere

f = P.dynamics()
solve_ivp(f, t, x0, ....)

So even P is gone f will survive and become independent from the model itself hence the method becomes a constructor. But maybe this is not what you want.

It seems like in the current version, the 't' argument isn't really "optional". In particular, if you don't want to specify 't' I think you would have to call the function as sys.dynamics(x=x, u=u), which looks funny to me.

An alternative would be to define dynamics and output so that they take a variable number of arguments, with the idea that dynamics(x, u) would set t=0 and dynamics(x) might set t=0 and u=0. You could also add a param argument.

If you do that change, I would leave _rhs and _out in place for iosys since you don't want the argument processing to have to be done each time you call the function (which happens a lot in a simulation).

Have decided I should start over and create a new PR.

In the meantime I thought I would do a quick benchmark to find out whether there are notable speed increases using a callable rather than the dynamics method, which has a few more attribute lookups. The results suggest that there's not much of a difference, so I think I'll just stick with having the method:

import control as ct
import time
sys = ct.ss([[-1.37394092,  0.08654488],
        [-0.47735689, -0.43755742]], [[-0.27547913],
        [ 0.        ]], [[-0.03284931, -0.        ]],  [[0.]] )

def dynamics(self, *args):
    if len(args) not in (2, 3): 
        raise ValueError("received"+len(args)+"args, expected 2 or 3")

    if np.size(args[1]) != self.nstates:
        raise ValueError("len(x) must be equal to number of states")
    if len(args) == 2: # received t and x, ignore t
        return self.A.dot(args[1]).reshape((-1,)) # return as row vector
    else: # received t, x, and u, ignore t
        #x = np.reshape(args[-2], (-1, 1)) # x must be column vector
        #u = np.reshape(args[-1], (-1, 1)) # u must be column vector
        if np.size(args[2]) != self.ninputs:
            raise ValueError("len(u) must be equal to number of inputs")
        return (self.A.dot(args[1]) + self.B.dot(args[2])).reshape((-1,)) # row
setattr(ct.StateSpace, 'dynamics', dynamics)

def dynamics_callable(sys): 
    return lambda t, x, u: sys.A.copy()@x + sys.B.copy()@u
dynamics = dynamics_callable(sys)

time1 = time.time()
x = sp.integrate.solve_ivp(sys.dynamics, (0,1000), (0,0),  args=((1,),))
print(time.time() - time1)
time2 = time.time()
x = sp.integrate.solve_ivp(dynamics, (0,1000), (0,0), args=((1,),))
print(time.time() - time2)

prints out:

0.1018521785736084
0.08774900436401367

@murrayrm

It seems like in the current version, the 't' argument isn't really "optional". In particular, if you don't want to specify 't' I think you would have to call the function as sys.dynamics(x=x, u=u), which looks funny to me.

makes sense.

An alternative would be to define dynamics and output so that they take a variable number of arguments, with the idea that dynamics(x, u) would set t=0 and dynamics(x) might set t=0 and u=0. You could also add a param argument.

One possible problem arises is if t is optional and two arguments are received: are they t and x or x and u?

To solve that ambiguity without keyword arguments, it seems like t has to be a required argument. And we just ignore it for LTI systems. Seem reasonable? not sure what else we could do there.

If you do that change, I would leave _rhs and _out in place for iosys since you don't want the argument processing to have to be done each time you call the function (which happens a lot in a simulation).

I can do that. Though it appears that arg processing doesn't actually seem to impact speed much, see above.

For _rhs versus and dynamics and arg processing: the optimal package calls _rhs thousands of times, so it might be worth checking the benchmarks in PR #549 and making sure they aren't affected.

I’ll add a note in the doc strings about the different anticipated use cases for _rhs vs dynamics (Fast vs user friendly) In the new version of this pr, #558.