Input Shaping

Rather the thumbnail to be the sensitivity plot

In the 1950s, Dr. Otto J. M. Smith conceived his Posicast¹ control method. He describes the method, “This is what happens when a fisherman drops his fly in the water at the maximum-position and zero velocity instant. Hence the descriptive name positive-cast or Posicast for this kind of system.”. The basic idea is outlined below in the figures I redrew from his publication.

Let’s say you have a simple gantry crane with the circle indicating the controllable hoist, Figure (a). You want move the payload, W, from the starting position to the final position f. You could push the forward button for the hoist until the payload arrives at f, but when you get there the payload would be swinging. This is highly undesirable. An uncontrolled swinging payload is the last thing a crane operator wants. Dr. Smith’s idea is to break up the command into two commands. The first command would move the hoist to position m, Figure (b). As the payload swings due to the motion, the hoist pauses until⊕“…one-half cycle of the natural transient period of the pendulum.” —Dr. Smith the payload reaches the final position. The second command will move the hoist quickly to the final position and the payload will experience little to zero sway. This concept is very simple in theory, but very powerful in controlling unwanted motion. Lets dig into the concept a little more and see how this idea became Input Shaping. Let’s begin by replicating the gantry crane model in Sympy:

from sympy import symbols, init_printing
import sympy
import sympy.physics.mechanics as me
from pydy.system import System
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
init_printing()

Let’s assign some arbitrary variables.

Length = 0.5 #meters
payload = 0.01 #kg
hoist = 1.0 #kg

# Create the Main Frame
A = me.ReferenceFrame('A')

# Create the symbols
x, y, theta, f = me.dynamicsymbols('x y theta f')
x_dot, y_dot, theta_dot = me.dynamicsymbols('x_dot y_dot theta_dot')

M, m, g, t = sympy.symbols('M m g t')
L1 = sympy.symbols('L1')

# Orient the other frames
B = A.orientnew('B', 'Axis', [theta, A.z])

# Create the Origin
Origin = me.Point('Origin')
Origin.set_pos(Origin, 0)
Origin.set_vel(A, 0)

# Create the hoist and payload points
P0 = me.Point('P0')
P1 = me.Point('P1')

# Set the hoist and payload positions
P0.set_pos(Origin, x * A.x)
P1.set_pos(P0, -L1 * B.y)

P0.set_vel(A, x_dot * A.x)
P1.v2pt_theory(P0, A, B)

I am not describing the hoist with a vertical dimension. It will only move in the X directions. The speed of the payload in Frame A: \[ x_{dot}\mathbf{\hat{a}_x} + L_{1} \dot{\theta}\mathbf{\hat{b}_x} \]

# Set up the kinematic differential equations
kde = [x_dot - x.diff(t),
       theta_dot - theta.diff(t)]

# Create the Particles
Pa0 = me.Particle('Pa0', P0, M)
Pa1 = me.Particle('Pa1', P1, m)

payload_gravity = (P1, -m * g * A.y)
payload_normal = (P1, m * g * B.y)
payload_on_hoist = (P0, -m * g * B.y)

forcing = (P0, f * A.x)

loads = [payload_gravity,
         payload_normal,
         payload_on_hoist,
         forcing]

This is how we will control the crane. This is a force (impulse) in Newtons that is applied to the hoist. The hoist has no friction so we need to apply a force in the opposite direction to stop it. The function allows us to specify the impulse amplitudes, the times applied, and the duration for each (although they need to be the same duration). Just to demonstrate what the force looks like, I’ll plot bang([5,-5],[1,3],1, sim_time).

def bang(Amps, Times, Interval, t):
    '''
    This function will issue impulses at the times given
    '''
    n = np.size(Amps)
    F = []
    for i in range(n):
        Fi = Amps[i]
        Ti = Times[i]
        fi = Fi * (t >= Ti) * (t <= Ti + Interval)
        F.append(fi)
    f = sum(F)
    return f

plt.plot(sim_time,bang([5,-5],[1,3],1, sim_time))

# Setting up the coordinates, speeds, and creating KanesMethod
coordinates = [x, theta]
speeds = [x_dot, theta_dot]
kane = me.KanesMethod(A, coordinates, speeds, kde)

# Creating Fr and Fr_star
fr, frstar = kane.kanes_equations([Pa0, Pa1], loads)

# Creating the PyDy System
sys = System(kane)

At first, we can tilt the payload and let it freely swing. Because there is no rotational damping crane’s oscillations will continue forever. Let’s run the simulation and look at the crane response.

# Assigning all the constants
runtime = 15

sys.constants = {m: payload,
                 M: hoist,
                 g: 9.81,
                 L1: Length}
sys.initial_conditions = {x:0, theta:0.78}
sys.specifieds={f:lambda y, t:bang([5,-5],[1,3],1, t)}

sys.times = np.linspace(0.0, runtime, runtime*30)
sys.generate_ode_function(generator='cython')
sim_time = np.linspace(0.0, runtime, runtime*30)
resp = sys.integrate()

The pendulum response will have a period very close to a simple pendulum. However, because the hoist can move—the pendulum’s motion will be influenced by the hoist. This is called a Multiple Degree of Freedom System. For a simple pendulum with a length of 0.5m, the natural frequency is \(\omega = \sqrt{\frac{gravity}{length}} \)Rad/s or 0.705Hz.

zeros = sim_time[np.argwhere(np.isclose(resp[:,1], 0, atol=0.1))]
x_resp = resp[:,0]
theta_resp = resp[:,1]
fig, ax = plt.subplots()

ax.scatter(zeros, np.zeros_like(zeros))
plt.plot(sim_time, np.rad2deg(resp[:,1]), label='Unshaped')
plt.xlabel(r'Time (s)')
plt.ylabel(r'Angle (deg)')

On the plot, the dots are plotted when the angle is zero. With these zeros we can calculate the period. The first dot occurs at 0.367s. A full period later, on the third dot, the time is 1.837s. This makes the period 1.47s or 0.68Hz. This is very close to our expected frequency of 0.7Hz. However, because of the dynamic coupling the payload will never exactly match a simple pendulum. Now, let’s rerun the simulation with the payload at rest, but let’s hit the hoist with our previous bang function. Below is the response of the hoist to the input we gave.

It is interesting to see the displacement curve form through the mathematics we outlined when creating the model. There are essentially 5 sections to this curve. Initially, the hoist is at rest at 0 distance. At 1s the force is applied and the distance increases with a slight curve—this indicates a constant acceleration, a linear increase in velocity, and a quadratic increase in displacement. At 2s, the force is stopped and the hoist continues at a constant velocity. At 3s, a force in the opposite direction is given and the hoist decelerates before finally resting at 10m.

The payload response is shown with the same colored sections. The payload is clearly affected by the gantry’s motion. After the gantry has stopped moving at 4s, the payload swings freely at almost 60deg off center. This is obviously unacceptable. A payload swinging freely like this would be an extreme hazard to anyone near the crane not to mention damaging to the crane itself. We will not be attempting to control the motion of the payload during the commands. The command shaping technique we are using is only concerned with the motion at the end of the commands. For this example, we want to stop the oscillation after t = 5s. What happens between t = 1 and t = 5 does not concern this technique.

After t =5s, the payload response is that of an undamped second-order system. Second-order systems are very well known and their responses have similar characteristics. In fact, we can generalize the payload response into a broad “second-order system” response. We can then use this information to make instructions on how to alter our bang command to stop the payload oscillations.

Command Shaping

The command given to the gantry can be viewed as two impulses. One to start the movement and the second to stop the movement. The inputs given have no understanding of the system dynamics. From our example, after our command is finished the payload continues to swing. Inputs only provide information to get from one state to another. Input shaping is a command shaping technique that will “encode” system dynamics within our inputs and stop the swinging with its own motion. Basically, you will give a command to the input shaper and it will spit out an altered or “shaped” command that cancels out it’s own oscillations. So let’s create the input shaper.

If an impulse is applied to a second-order system, \(A_0\), at time \(t_0\) the response will be of the form²:

\[ y(t)=\left[\frac{A_0 \omega}{\sqrt{1-\zeta^2}}e^{-\zeta\omega(t-t_0)}\right]\sin\left(\omega\sqrt{1-\zeta^2}(t-t_0)\right) \]

Plotted below is the general response of a second-order system where \(A_0=1\) \(\omega=0.68\)Hz \(\zeta=0\) and \(t_0=1\)s. The response has the same natural frequency as our system, but it does not exactly replicate our payload response. Explain that we have not included damping and how the impulse needs to be decreased by a certain amount for the next impulse due to the damping

In order to control our response we require an additional impulse. This is one of the drawbacks of the input shaping technique. Because we require an additional impulse, our command time will be increased. From our gantry example, we gave an impulse of 1s. Once applied, our shaper will increase this time by a small amount. For a second-order system, the response to a series of impulses is:

\[ y_\Sigma(t)=\sum_{i=1}^{n}\left[\frac{A_i \omega}{\sqrt{1-\zeta^2}}e^{-\zeta\omega(t-t_i)}\right]\sin\left(\omega\sqrt{1-\zeta^2}(t-t_i)\right) \]

This is the superposition principle. Given two responses, the result is a summation of the two responses. This leads to the “big idea” resulting from wave interference. So let create another plot to understand where the input shaping technique is headed. Let’s plot two impulses \(A_i=[1,-1]\) where the first applied at \(t_1=1\) and the second, an impulse in the opposite direction, a full period later at \(t_2=0.68*2*\pi\). Additionally, the resultant from the summation of the two responses is plotted.

The two responses are directly out of phase and when summed produce zero. The response to a series of impulses gives us the oscillatory nature of second-order systems. The equation also provides \(A_\Sigma\), the amplitude of residual vibration. Residual vibration refers to the motion of the system after an input is given. This is exactly what we are trying to stop, vibration after our command. What follows is a rearrangement of the terms to solve for \(A_\Sigma\).

Using the following trigonometric identity:

\[ \sum_{i=1}^{n} B_i \sin(\omega t + \phi_i) = A_\Sigma \sin (\omega t + \psi) \]

Where the amplitude of residual vibration is

\[ A_\Sigma = \sqrt{\left( \sum_{i=1}^{n} B_i \cos(\phi_i)\right)^2 + \left( \sum_{i=1}^{n} B_i \sin(\phi_i)\right)^2} \]

The coefficients from the response to a series of impulses using the trig identity are:

\[ B_i = \frac{A_i \omega}{\sqrt{1-\zeta^2}}e^{-\zeta\omega(t-t_i)} \]

Using \(B_i\) and moving the constants out of the square root, \(A_\Sigma\) evaluated at the time of the last impulse, \(t=t_n\), is:

\[ A_\Sigma = \frac{\omega}{\sqrt{1-\zeta^2}}e^{-\zeta\omega t_n} \sqrt{\left[ C(\omega,\zeta)\right]^2 + \left[ S(\omega,\zeta)\right]^2} \]

where,

\[ C(\omega, \zeta) = \sum_{i=1}^{n} A_i e^{\zeta \omega t_i}\cos(\omega_d t_i) \]

\[ S(\omega, \zeta) = \sum_{i=1}^{n} A_i e^{\zeta \omega t_i}\sin(\omega_d t_i) \]

and

\[ \omega_d = \omega \sqrt{1 - \zeta^2} \]

The residual vibration amplitude is given with units corresponding to the system. However, if we wanted a dimensionless value which relays the amount of vibration, \(A_\Sigma\) needs to be normalized. The resulting value, \(V(\omega, \zeta)\), the percentage residual vibration “gives us the ratio of vibration with input shaping to that without input shaping”². Applying a single impulse of 1 units (unity magnitude) at time 0, \(A_\Sigma\) will reduce to:

\[ A_\uparrow = \frac{\omega}{\sqrt{1-\zeta^2}} \]

Which allows \(A_\Sigma\) to be normalized and creates the percent residual vibration:

\[ V(\omega, \zeta) = \frac{A_\Sigma}{A_\uparrow} = e^{-\zeta\omega t_n} \sqrt{\left[ C(\omega,\zeta)\right]^2 + \left[ S(\omega,\zeta)\right]^2} \]

If \(V(\omega, \zeta)\) is set to zero, which means we want zero residual vibration, the amplitudes and times can be solved for. However, this is totally open ended so we could end up with infinite amplitudes and shapers which continue forever. Applying constraints stops this behavior. The first limit is to make sure our shaped input and unshaped input reach the same set-point. That is, we do not want the shaped input to end up in a different position from the unshaped input. This requires the total amplitudes created by the shaper to add to one³:

\[ \sum_{i=1}^{n} A_i = 1 \]

and we require that the shaper is positive:

\[ A_i > 0, \quad i = 1,2,…,n \]

For example, to solve for a zero vibration, two-impulse sequence set \(V(\omega, \zeta)=0\). This requires \(C(\omega,\zeta)\) and \(S(\omega,\zeta)\) to equal 0. Additionally, the first impulse will be applied at the same time the input begins at \(t_1=0\). This will simplify to:

\[ C(\omega, \zeta) = 0 = A_1 + A_2 e^{\zeta \omega t_2}\cos(\omega_d t_2) \]

\[ S(\omega, \zeta) = 0 = A_2 e^{\zeta \omega t_2}\sin(\omega_d t_2) \]

and three unknowns \(A_1 \), \(A_2\), and \(t_2\)

Because the system is oscillatory, \(S(\omega,\zeta)\) is zero for every multiple of the damped half period of oscillation³, \(\tau_d\).

\[ \omega_d t_2 = n\pi,\ \Rightarrow t_2 = \frac{n\pi}{\omega_d} = \frac{n\tau_d}{2}, \quad n=1,2,… \]

There are an infinite number of solutions. In order to select the shortest shaper duration, the last impulse time is minimized, \(min(t_n)\). For this example, the smallest value of \(t_2\) occurs at \(\frac{\tau_d}{2}\). Constraining the amplitudes to sum to one and the minimized value of \(t_2\), \(C(\omega,\zeta)\) can now be solved for \(A_1\) and \(A_2\) to find:

\[ A_1 = \frac{\exp{\left(\frac{\zeta \pi}{\sqrt{1-\zeta^2}}\right)}}{1 + \exp{\left(\frac{\zeta \pi}{\sqrt{1-\zeta^2}}\right)}} \]

and

\[ A_2 = \frac{1}{1 + \exp{\left(\frac{\zeta \pi}{\sqrt{1-\zeta^2}}\right)}} \]

Defining \(K=\exp{\Big(\frac{-\zeta \pi}{\sqrt{1-\zeta^2}}}\Big)\), the closed-form solution of the Zero-Vibration ZV shaper is then⁴:

\[ \begin{bmatrix} A_i\\ t_i \end{bmatrix}= \begin{bmatrix} \frac{1}{1+K} & \frac{K}{1+K} \\ 0 & \frac{\tau_d}{2} \end{bmatrix} \]

Easy enough, let’s create the shaper based on our gantry model: \(\omega=4.396\), \(\zeta=0.0456\), and \(\tau_d=1.43\)

\[ \begin{bmatrix} A_i\\ t_i \end{bmatrix}= \begin{bmatrix} 0.536 & 0.464 \\ 0 & 1.43 \end{bmatrix} \]

Solve for the response Show the percent error plot rerun the values and show the response (consider adding damping) Hint at other shapers Hint at multi mode shapers Conclude

… to be continued