Robotics Application of a Method for Analytically Computing Infinitesimal
111
of a basis vector at a sufficiently small magnitude (1E-2–1E-4 are common values [17])
and re-initialized. The system is then integrated after this perturbation for a sufficiently
long time (typically 10 times the period). The value of the iPRC for the phase at which
the perturbation is applied is calculated as the difference between the duration of the
unperturbed and perturbed limit cycles, divided by the magnitude of the perturbation.
Analytically Finding the Infinitesimal Phase Response Curve. For a general smooth
nonlinear system which produces a limit cycle, the iPRC can be found using an adjoint
equation method [19]. The adjoint equation is
˙
z = −(DF(γ ))
T
· z,
(16)
where z is the iPRC and DF(γ ) is the Jacobian matrix evaluated on the limit cycle γ .
Usually, the adjoint equation cannot be used on nonsmooth systems due to the discontinuities present at the boundaries of piecewise functions or sliding conditions. However,
[16, 17] introduce methods to account for these discontinuities.
To account for the discontinuities at the boundaries, jump matrices are calculated at
each boundary crossing. These jump matrices capture the abrupt changes that happen
to the vector field of the system at the boundaries. To use the jump matrices to properly
calculate the iPRC, one evaluates the adjoint equation backwards in time. The adjoint
equation is calculated backwards in time because the jump matrix at the liftoff point for
sliding conditions would not be well defined otherwise, as detailed in Remark 3.9 and
Theorem 3.7 in [16]. When a boundary is encountered, the iPRC just before (in forwards
time) the boundary is calculated by
z
−
= J z
+
,
(17)
where z − is the iPRC just before the boundary, z + is the iPRC just after the boundary
and J is the jump matrix. How to calculate the jump matrix is detailed in Theorem 2.2
for piecewise smooth limit cycles in [17] and Remark 3.10 for LCSCs in [16]. Once z −
is found, one re-initializes the adjoint equation using z − and continues to evaluate the
adjoint equation backwards in time.
3 Results
3.1 Simulation of the Full System
The state equations shown in Table 1 were simulated in MATLAB (The Mathworks,
Natick, MA) using its ODE15s solver with relative and absolute tolerances set at 1 ×
10 −10 , for 5000 ms to produce Fig. 1. Figure 1 shows that the system is rhythmically
oscillating. The flattening of the peaks in plot (c) in Fig. 1 verifies and displays the
locations of the system’s sliding conditions, where the motor’s current draw is limited.
Because the motor’s current draw is a dynamical state variable, this limitation represents
a sliding condition, not a discontinuity. The period of the system was calculated as
T ≈ 545ms.
111
of a basis vector at a sufficiently small magnitude (1E-2–1E-4 are common values [17])
and re-initialized. The system is then integrated after this perturbation for a sufficiently
long time (typically 10 times the period). The value of the iPRC for the phase at which
the perturbation is applied is calculated as the difference between the duration of the
unperturbed and perturbed limit cycles, divided by the magnitude of the perturbation.
Analytically Finding the Infinitesimal Phase Response Curve. For a general smooth
nonlinear system which produces a limit cycle, the iPRC can be found using an adjoint
equation method [19]. The adjoint equation is
˙
z = −(DF(γ ))
T
· z,
(16)
where z is the iPRC and DF(γ ) is the Jacobian matrix evaluated on the limit cycle γ .
Usually, the adjoint equation cannot be used on nonsmooth systems due to the discontinuities present at the boundaries of piecewise functions or sliding conditions. However,
[16, 17] introduce methods to account for these discontinuities.
To account for the discontinuities at the boundaries, jump matrices are calculated at
each boundary crossing. These jump matrices capture the abrupt changes that happen
to the vector field of the system at the boundaries. To use the jump matrices to properly
calculate the iPRC, one evaluates the adjoint equation backwards in time. The adjoint
equation is calculated backwards in time because the jump matrix at the liftoff point for
sliding conditions would not be well defined otherwise, as detailed in Remark 3.9 and
Theorem 3.7 in [16]. When a boundary is encountered, the iPRC just before (in forwards
time) the boundary is calculated by
z
−
= J z
+
,
(17)
where z − is the iPRC just before the boundary, z + is the iPRC just after the boundary
and J is the jump matrix. How to calculate the jump matrix is detailed in Theorem 2.2
for piecewise smooth limit cycles in [17] and Remark 3.10 for LCSCs in [16]. Once z −
is found, one re-initializes the adjoint equation using z − and continues to evaluate the
adjoint equation backwards in time.
3 Results
3.1 Simulation of the Full System
The state equations shown in Table 1 were simulated in MATLAB (The Mathworks,
Natick, MA) using its ODE15s solver with relative and absolute tolerances set at 1 ×
10 −10 , for 5000 ms to produce Fig. 1. Figure 1 shows that the system is rhythmically
oscillating. The flattening of the peaks in plot (c) in Fig. 1 verifies and displays the
locations of the system’s sliding conditions, where the motor’s current draw is limited.
Because the motor’s current draw is a dynamical state variable, this limitation represents
a sliding condition, not a discontinuity. The period of the system was calculated as
T ≈ 545ms.