Pendulum swing-up

A pendulum swung up to its inverted position by a bang-bang acceleration of its pivot.

• Name: pendulum_swing_up
• Field: engineering/control
• Formalism: ordinary differential equations
• References: Åström and Furuta 2000; Dilão 2023

Factory

@register_model(
    field='engineering/control',
    tags=('control', 'bang-bang', 'pendulum', 'energy', 'ode'),
    references=('Åström and Furuta 2000', 'Dilão 2023'),
    defaults={'until': 100.0, 'dt': 0.1, 'solver': RK4Solver(dt=0.01)},
)
def pendulum_swing_up(
    theta: float = 0.0,
    omega0: float = 0.5,
    *,
    omega: float = 1.5652475842498528,
    mu: float = 0.1,
) -> System:
    """Return a pendulum swung up to its inverted position by a bang-bang acceleration of its pivot.

    The pivot of an acrobat's pendulum moves horizontally with acceleration `f'' = u`, so the pendulum follows
    `θ'' + ω^2 sin θ + u cos θ = 0`. The control `u = -μ sign(θ' cos θ) sign(E_s - E)` pumps the energy
    `E = θ'^2 / 2 - ω^2 cos θ` towards `E_s = ω^2`, the energy of the inverted position, from below or from above.

    Args:
        theta:
            The initial angle from the hanging position.

        omega0:
            The initial angular velocity.

        omega:
            The natural frequency `sqrt(g / l)`, by default for `g = 9.8` and `l = 4`.

        mu:
            The size of the acceleration of the pivot.

    Returns:
        system:
            The model, with the angle `pendulum.theta`, the position of the pivot `pendulum.f`, and their velocities.

    Examples:
        >>> from skcomplex.models.engineering.control import pendulum_swing_up
        >>> pendulum_swing_up(mu=0.2).elements[0].mu
        0.2
    """
    pendulum = SwingingPendulumElement('pendulum', theta=theta, theta_d1=omega0, f=0.0, f_d1=0.0, omega=omega, mu=mu)
    return System(elements=[pendulum])

Components

class SwingingPendulumElement(Element):
    """Pendulum swung up by a bang-bang control of its energy, with `f` integrating the control."""

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        theta: Annotated[object, StateSpec(order=2)] = 0.0,
        theta_d1: float = 0.0,
        f: Annotated[object, StateSpec(order=2)] = 0.0,
        f_d1: float = 0.0,
        omega: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 1.5652475842498528,
        mu: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    ) -> None:
        super().__init__(label, visual=visual)
        self.theta = theta
        self.theta_d1 = theta_d1
        self.f = f
        self.f_d1 = f_d1
        self.omega = omega
        self.mu = mu

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the controlled rate of the angle and the rate of the control integral."""
        return [
            RateLaw(
                s.theta,
                -(s.omega**2) * sin(s.theta) - SwingingPendulumElement._swing_up_control(s) * cos(s.theta),
            ),
            RateLaw(s.f, SwingingPendulumElement._swing_up_control(s)),
        ]

    @staticmethod
    def _pendulum_energy(s: SimpleNamespace) -> Expression:
        """Return the energy per unit inertia of the pendulum, `θ'^2 / 2 - ω^2 cos θ`."""
        return derivative(s.theta, 1) ** 2 / 2 - s.omega**2 * cos(s.theta)

    @staticmethod
    def _swing_up_control(s: SimpleNamespace) -> Expression:
        """Return the bang-bang control `u = -μ sign(θ' cos θ) sign(E_s - E)` of the energy of the pendulum."""
        gap = s.omega**2 - SwingingPendulumElement._pendulum_energy(s)
        return -s.mu * sign(derivative(s.theta, 1) * cos(s.theta)) * sign(gap)

Equations

pendulum.theta: diff(theta_pendulum(t), t, t) == mu_pendulum*cos(theta_pendulum(t))*where(omega_pendulum^2*cos(theta_pendulum(t)) + omega_pendulum^2 - 1/2*diff(theta_pendulum(t), t)^2 > 0, 1, where(omega_pendulum^2*cos(theta_pendulum(t)) + omega_pendulum^2 - 1/2*diff(theta_pendulum(t), t)^2 < 0, -1, 0))*where(cos(theta_pendulum(t))*diff(theta_pendulum(t), t) > 0, 1, where(cos(theta_pendulum(t))*diff(theta_pendulum(t), t) < 0, -1, 0)) - omega_pendulum^2*sin(theta_pendulum(t))
pendulum.f: diff(f_pendulum(t), t, t) == -mu_pendulum*where(cos(theta_pendulum(t))*diff(theta_pendulum(t), t) > 0, 1, where(cos(theta_pendulum(t))*diff(theta_pendulum(t), t) < 0, -1, 0))*where(omega_pendulum^2*cos(theta_pendulum(t)) + omega_pendulum^2 - 1/2*diff(theta_pendulum(t), t)^2 > 0, 1, where(omega_pendulum^2*cos(theta_pendulum(t)) + omega_pendulum^2 - 1/2*diff(theta_pendulum(t), t)^2 < 0, -1, 0))