Thread: [python-control-discuss] forced response problem for large dt

Apparently I don't know which email list is subscribed.  Sorry if this
comes through more than once.

I am teaching a classical controls class for mechanical engineering
undergraduates.  Up till now, I have tried to gloss over state-space.  My
students need to do an initial condition simulation for a system that
includes actuator saturation.  I have them do this by integrating for one
time step at a time using control.forced_response inside a for loop.  We
are essentially doing a continuous time approximation of ZOH with the input
held constant for each time step.  As the time step gets too large, sp.
linalg.expm has to use a higher order pade approximation and eventually
throws this error:

ValueError                                Traceback (most recent call last)

/Users/rkrauss/git/python-control/zumo_PID_simulation_modified_ss.py
in <module>()

    103     t0 = dt*(i-1)

    104     t1 = dt*i

--> 105
  to, yo, xo = control.forced_response(G_int, [t0,t1], [v[i],v[i]], X0=X0)

    106     X0 = xo[:,-1]#<-- save for next time through for loop

    107     x[i] = squeeze(X0)


/Users/rkrauss/git/python-control/control/timeresp.pyc in forced_response(sys,
T, U, X0, transpose)

    374                          [np.zeros((n_inputs, n_states + 2 *
n_inputs))]])

    375             print('M=' + str(M))

--> 376             expM = sp.linalg.expm(M)

    377             Ad = expM[:n_states, :n_states]

    378             Bd1 = expM[:n_states, n_states+n_inputs:]


/usr/local/lib/python2.7/site-packages/scipy/linalg/matfuncs.pyc in expm(A,
q)

    258     # Input checking and conversion is provided by
sparse.linalg.expm().

    259     import scipy.sparse.linalg

--> 260     return scipy.sparse.linalg.expm(A)

    261

    262


/usr/local/lib/python2.7/site-packages/scipy/sparse/linalg/matfuncs.pyc
in expm(A)

    580

    581     """

--> 582     return _expm(A, use_exact_onenorm='auto')

    583

    584


/usr/local/lib/python2.7/site-packages/scipy/sparse/linalg/matfuncs.pyc
in _expm(A,
use_exact_onenorm)

    635     if structure == UPPER_TRIANGULAR:

    636         # Invoke Code Fragment 2.1.

--> 637         X = _fragment_2_1(X, h.A, s)

    638     else:

    639         # X = r_13(A)^(2^s) by repeated squaring.


/usr/local/lib/python2.7/site-packages/scipy/sparse/linalg/matfuncs.pyc
in _fragment_2_1(X,
T, s)

    753     exp_diag = np.exp(scale * diag_T)

    754     for k in range(n):

--> 755         X[k, k] = exp_diag[k]

    756

    757     for i in range(s-1, -1, -1):


ValueError: setting an array element with a sequence.


This is probably ultimately a problem for the scipy people, but my
students' project is due in 6 days.  Any suggestions to quickly get the
attached simulation code to work for 60Hz simulation, i.e. dt = 1.0/60?


Thanks,


Ryan

--
Ryan Krauss, Ph.D.
Associate Professor
Mechanical Engineering
Southern Illinois University Edwardsville

My hackish solution to the immediate problem is to integrate 3 times at
dt/3 within the for loop:

Ndt = 3#factor to divide dt by
dt_sim = dt/Ndt

for i in range(1,N):
    e[i] = r[i] - y[i-3]#<-- one step time delay on the measurement
    #v[i] = Gc*e[i]#<-- this is only doing P control
    esum += e[i]
    v[i] = Kp*e[i] + Ki*esum + Kd*(e[i]-e[i-1]) #PID control
    if include_sat:
        v[i] = mysat(v[i])
    t0 = dt*(i-1)

    #integrate Ndt times to get to the next dt
    for q in range(Ndt):
        t1 = t0 + dt_sim
        to, yo, xo = control.forced_response(G_int, [t0,t1], [v[i],v[i]],
X0=X0)
        X0 = xo[:,-1]#<-- save for next time through for loop
        t0 = t1#<-- for next pass

    x[i] = squeeze(X0)
    y[i] = yo[-1]

This keep my students moving forward for now.


--
Ryan Krauss, Ph.D.
Associate Professor
Mechanical Engineering
Southern Illinois University Edwardsville

On Fri, Dec 4, 2015 at 3:21 PM, Ryan Krauss <rk...@si...> wrote:

> Apparently I don't know which email list is subscribed.  Sorry if this
> comes through more than once.
>
> I am teaching a classical controls class for mechanical engineering
> undergraduates.  Up till now, I have tried to gloss over state-space.  My
> students need to do an initial condition simulation for a system that
> includes actuator saturation.  I have them do this by integrating for one
> time step at a time using control.forced_response inside a for loop.  We
> are essentially doing a continuous time approximation of ZOH with the input
> held constant for each time step.  As the time step gets too large, sp.
> linalg.expm has to use a higher order pade approximation and eventually
> throws this error:
>
> ValueError                                Traceback (most recent call last)
>
>
> /Users/rkrauss/git/python-control/zumo_PID_simulation_modified_ss.py in <module>()
>
>     103     t0 = dt*(i-1)
>
>     104     t1 = dt*i
>
> --> 105
>   to, yo, xo = control.forced_response(G_int, [t0,t1], [v[i],v[i]], X0=X0)
>
>     106     X0 = xo[:,-1]#<-- save for next time through for loop
>
>     107     x[i] = squeeze(X0)
>
>
> /Users/rkrauss/git/python-control/control/timeresp.pyc in forced_response(sys,
> T, U, X0, transpose)
>
>     374                          [np.zeros((n_inputs, n_states + 2 *
> n_inputs))]])
>
>     375             print('M=' + str(M))
>
> --> 376             expM = sp.linalg.expm(M)
>
>     377             Ad = expM[:n_states, :n_states]
>
>     378             Bd1 = expM[:n_states, n_states+n_inputs:]
>
>
> /usr/local/lib/python2.7/site-packages/scipy/linalg/matfuncs.pyc in expm(A,
> q)
>
>     258     # Input checking and conversion is provided by
> sparse.linalg.expm().
>
>     259     import scipy.sparse.linalg
>
> --> 260     return scipy.sparse.linalg.expm(A)
>
>     261
>
>     262
>
>
>
> /usr/local/lib/python2.7/site-packages/scipy/sparse/linalg/matfuncs.pyc in expm(A)
>
>     580
>
>     581     """
>
> --> 582     return _expm(A, use_exact_onenorm='auto')
>
>     583
>
>     584
>
>
> /usr/local/lib/python2.7/site-packages/scipy/sparse/linalg/matfuncs.pyc in _expm(A,
> use_exact_onenorm)
>
>     635     if structure == UPPER_TRIANGULAR:
>
>     636         # Invoke Code Fragment 2.1.
>
> --> 637         X = _fragment_2_1(X, h.A, s)
>
>     638     else:
>
>     639         # X = r_13(A)^(2^s) by repeated squaring.
>
>
> /usr/local/lib/python2.7/site-packages/scipy/sparse/linalg/matfuncs.pyc in _fragment_2_1(X,
> T, s)
>
>     753     exp_diag = np.exp(scale * diag_T)
>
>     754     for k in range(n):
>
> --> 755         X[k, k] = exp_diag[k]
>
>     756
>
>     757     for i in range(s-1, -1, -1):
>
>
> ValueError: setting an array element with a sequence.
>
>
> This is probably ultimately a problem for the scipy people, but my
> students' project is due in 6 days.  Any suggestions to quickly get the
> attached simulation code to work for 60Hz simulation, i.e. dt = 1.0/60?
>
>
> Thanks,
>
>
> Ryan
>
> --
> Ryan Krauss, Ph.D.
> Associate Professor
> Mechanical Engineering
> Southern Illinois University Edwardsville
>

At one point in time we were doing the simulations in control.forced_response differently, using the scipy.integrate.odeint() function.  A brief discussion is here:

  https://github.com/python-control/python-control/issues/48

plus look at the changes here

  https://github.com/python-control/python-control/commit/d7d278ba6072fce1ef28402b7580ffa698424f76#diff-e216ba2d66c950242b1475a045e3cd33

So one other possibility, separate from you quick fix, would be to call scipy.integrate.odeint().


-richard

> On 4 Dec 15, at 13:21, Ryan Krauss <rk...@si...> wrote:
> 
> Apparently I don't know which email list is subscribed.  Sorry if this comes through more than once.
> 
> I am teaching a classical controls class for mechanical engineering undergraduates.  Up till now, I have tried to gloss over state-space.  My students need to do an initial condition simulation for a system that includes actuator saturation.  I have them do this by integrating for one time step at a time using control.forced_response inside a for loop.  We are essentially doing a continuous time approximation of ZOH with the input held constant for each time step.  As the time step gets too large, sp.linalg.expm has to use a higher order pade approximation and eventually throws this error:
> 
> ValueError                                Traceback (most recent call last)
> 
> /Users/rkrauss/git/python-control/zumo_PID_simulation_modified_ss.py in <module>()
> 
>     103     t0 = dt*(i-1)
> 
>     104     t1 = dt*i
> 
> --> 105     to, yo, xo = control.forced_response(G_int, [t0,t1], [v[i],v[i]], X0=X0)
> 
>     106     X0 = xo[:,-1]#<-- save for next time through for loop
> 
>     107     x[i] = squeeze(X0)
> 
> 
> 
> /Users/rkrauss/git/python-control/control/timeresp.pyc in forced_response(sys, T, U, X0, transpose)
> 
>     374                          [np.zeros((n_inputs, n_states + 2 * n_inputs))]])
> 
>     375             print('M=' + str(M))
> 
> --> 376             expM = sp.linalg.expm(M)
> 
>     377             Ad = expM[:n_states, :n_states]
> 
>     378             Bd1 = expM[:n_states, n_states+n_inputs:]
> 
> 
> 
> /usr/local/lib/python2.7/site-packages/scipy/linalg/matfuncs.pyc in expm(A, q)
> 
>     258     # Input checking and conversion is provided by sparse.linalg.expm().
> 
>     259     import scipy.sparse.linalg
> 
> --> 260     return scipy.sparse.linalg.expm(A)
> 
>     261 
> 
>     262 
> 
> 
> 
> /usr/local/lib/python2.7/site-packages/scipy/sparse/linalg/matfuncs.pyc in expm(A)
> 
>     580 
> 
>     581     """
> 
> --> 582     return _expm(A, use_exact_onenorm='auto')
> 
>     583 
> 
>     584 
> 
> 
> 
> /usr/local/lib/python2.7/site-packages/scipy/sparse/linalg/matfuncs.pyc in _expm(A, use_exact_onenorm)
> 
>     635     if structure == UPPER_TRIANGULAR:
> 
>     636         # Invoke Code Fragment 2.1.
> 
> --> 637         X = _fragment_2_1(X, h.A, s)
> 
>     638     else:
> 
>     639         # X = r_13(A)^(2^s) by repeated squaring.
> 
> 
> 
> /usr/local/lib/python2.7/site-packages/scipy/sparse/linalg/matfuncs.pyc in _fragment_2_1(X, T, s)
> 
>     753     exp_diag = np.exp(scale * diag_T)
> 
>     754     for k in range(n):
> 
> --> 755         X[k, k] = exp_diag[k]
> 
>     756 
> 
>     757     for i in range(s-1, -1, -1):
> 
> 
> 
> ValueError: setting an array element with a sequence.
> 
> 
> 
> This is probably ultimately a problem for the scipy people, but my students' project is due in 6 days.  Any suggestions to quickly get the attached simulation code to work for 60Hz simulation, i.e. dt = 1.0/60?
> 
> 
> 
> Thanks,
> 
> 
> 
> Ryan
> 
> 
> --
> Ryan Krauss, Ph.D.
> Associate Professor
> Mechanical Engineering
> Southern Illinois University Edwardsville
> <kp_0_3_rotate_only_less_delay.csv><zumo_PID_simulation_modified_ss.py>------------------------------------------------------------------------------
> Go from Idea to Many App Stores Faster with Intel(R) XDK
> Give your users amazing mobile app experiences with Intel(R) XDK.
> Use one codebase in this all-in-one HTML5 development environment.
> Design, debug & build mobile apps & 2D/3D high-impact games for multiple OSs.
> http://pubads.g.doubleclick.net/gampad/clk?id=254741911&iu=/4140_______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss

Thanks for the suggestion.


--
Ryan Krauss, Ph.D.
Associate Professor
Mechanical Engineering
Southern Illinois University Edwardsville

On Sat, Dec 5, 2015 at 12:47 PM, Richard Murray <mu...@cd...>
wrote:

> At one point in time we were doing the simulations in
> control.forced_response differently, using the scipy.integrate.odeint()
> function.  A brief discussion is here:
>
>   https://github.com/python-control/python-control/issues/48
>
> plus look at the changes here
>
>
> https://github.com/python-control/python-control/commit/d7d278ba6072fce1ef28402b7580ffa698424f76#diff-e216ba2d66c950242b1475a045e3cd33
>
> So one other possibility, separate from you quick fix, would be to call
> scipy.integrate.odeint().
>
>
> -richard
>
> On 4 Dec 15, at 13:21, Ryan Krauss <rk...@si...> wrote:
>
> Apparently I don't know which email list is subscribed.  Sorry if this
> comes through more than once.
>
> I am teaching a classical controls class for mechanical engineering
> undergraduates.  Up till now, I have tried to gloss over state-space.  My
> students need to do an initial condition simulation for a system that
> includes actuator saturation.  I have them do this by integrating for one
> time step at a time using control.forced_response inside a for loop.  We
> are essentially doing a continuous time approximation of ZOH with the input
> held constant for each time step.  As the time step gets too large, sp.
> linalg.expm has to use a higher order pade approximation and eventually
> throws this error:
>
> ValueError                                Traceback (most recent call last)
>
>
> /Users/rkrauss/git/python-control/zumo_PID_simulation_modified_ss.py in <module>()
>
>     103     t0 = dt*(i-1)
>
>     104     t1 = dt*i
>
> --> 105
>   to, yo, xo = control.forced_response(G_int, [t0,t1], [v[i],v[i]], X0=X0)
>
>     106     X0 = xo[:,-1]#<-- save for next time through for loop
>
>     107     x[i] = squeeze(X0)
>
>
> /Users/rkrauss/git/python-control/control/timeresp.pyc in forced_response(sys,
> T, U, X0, transpose)
>
>     374                          [np.zeros((n_inputs, n_states + 2 *
> n_inputs))]])
>
>     375             print('M=' + str(M))
>
> --> 376             expM = sp.linalg.expm(M)
>
>     377             Ad = expM[:n_states, :n_states]
>
>     378             Bd1 = expM[:n_states, n_states+n_inputs:]
>
>
> /usr/local/lib/python2.7/site-packages/scipy/linalg/matfuncs.pyc in expm(A,
> q)
>
>     258     # Input checking and conversion is provided by
> sparse.linalg.expm().
>
>     259     import scipy.sparse.linalg
>
> --> 260     return scipy.sparse.linalg.expm(A)
>
>     261
>
>     262
>
>
>
> /usr/local/lib/python2.7/site-packages/scipy/sparse/linalg/matfuncs.pyc in expm(A)
>
>     580
>
>     581     """
>
> --> 582     return _expm(A, use_exact_onenorm='auto')
>
>     583
>
>     584
>
>
> /usr/local/lib/python2.7/site-packages/scipy/sparse/linalg/matfuncs.pyc in _expm(A,
> use_exact_onenorm)
>
>     635     if structure == UPPER_TRIANGULAR:
>
>     636         # Invoke Code Fragment 2.1.
>
> --> 637         X = _fragment_2_1(X, h.A, s)
>
>     638     else:
>
>     639         # X = r_13(A)^(2^s) by repeated squaring.
>
>
> /usr/local/lib/python2.7/site-packages/scipy/sparse/linalg/matfuncs.pyc in _fragment_2_1(X,
> T, s)
>
>     753     exp_diag = np.exp(scale * diag_T)
>
>     754     for k in range(n):
>
> --> 755         X[k, k] = exp_diag[k]
>
>     756
>
>     757     for i in range(s-1, -1, -1):
>
>
> ValueError: setting an array element with a sequence.
>
>
> This is probably ultimately a problem for the scipy people, but my
> students' project is due in 6 days.  Any suggestions to quickly get the
> attached simulation code to work for 60Hz simulation, i.e. dt = 1.0/60?
>
>
> Thanks,
>
>
> Ryan
>
> --
> Ryan Krauss, Ph.D.
> Associate Professor
> Mechanical Engineering
> Southern Illinois University Edwardsville
> <kp_0_3_rotate_only_less_delay.csv><zumo_PID_simulation_modified_ss.py>
> ------------------------------------------------------------------------------
> Go from Idea to Many App Stores Faster with Intel(R) XDK
> Give your users amazing mobile app experiences with Intel(R) XDK.
> Use one codebase in this all-in-one HTML5 development environment.
> Design, debug & build mobile apps & 2D/3D high-impact games for multiple
> OSs.
>
> http://pubads.g.doubleclick.net/gampad/clk?id=254741911&iu=/4140_______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>
>
>
>
> ------------------------------------------------------------------------------
> Go from Idea to Many App Stores Faster with Intel(R) XDK
> Give your users amazing mobile app experiences with Intel(R) XDK.
> Use one codebase in this all-in-one HTML5 development environment.
> Design, debug & build mobile apps & 2D/3D high-impact games for multiple
> OSs.
> http://pubads.g.doubleclick.net/gampad/clk?id=254741911&iu=/4140
> _______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>
>