Non-Holonomic Constraints

jupyter nbconvert --to html notebook.ipynb

Vincoli in meccanica analitica

Un $vincolo$ limita le configurazioni o i moti possibili di un sistema meccanico. Siano $ q_1,\dots,q_n $ le coordinate generalizzate, $ \dot q_i $ le velocità e $ t $ il tempo.

Vincoli olonomi

Definizione

Un vincolo è detto $olonomo$ se può essere espresso come $$ f(q_1,\dots,q_n,t) = 0, $$ cioè dipende esclusivamente dalle coordinate generalizzate (e, al più, dal tempo).

Esempi

• Punto vincolato su una circonferenza:

$$\large x^2 + y^2 - R^2 = 0 $$

• Pendolo rigido:

$$\large x^2 + y^2 - l^2 = 0 $$

Trattamento lagrangiano

Nel caso di vincoli olonomi esistono due metodi equivalenti.

Riduzione delle coordinate Si scelgono coordinate generalizzate che soddisfano automaticamente il vincolo. In questo modo si riduce il numero di gradi di libertà e il vincolo non compare esplicitamente nelle equazioni del moto.

Moltiplicatori di Lagrange Si introduce un moltiplicatore $ \lambda $ e si definisce il Lagrangiano esteso $$\large \mathcal{L}' = L(q,\dot q,t) + \lambda f(q,t). $$

Le equazioni di Lagrange assumono la forma $$\large \frac{d}{dt}\frac{\partial L}{\partial \dot q_i}\frac{\partial L}{\partial q_i}= \lambda \frac{\partial f}{\partial q_i}. $$

Il termine proporzionale a $ \lambda $ rappresenta la forza vincolare.

Vincoli non olonomi

Definizione

Un vincolo è detto $non olonomo$ se non è integrabile in una relazione che coinvolga soltanto le coordinate. Tipicamente assume la forma $$\large \sum_i a_i(q,t)\,\dot q_i + a_0(q,t) = 0. $$

Tali vincoli dipendono dalle velocità e limitano lo spazio delle velocità ammissibili, ma non quello delle configurazioni.

Esempi

• Rotolamento senza strisciamento
• Disco che rotola su un piano
• Veicolo che può muoversi solo nella direzione delle ruote

Difficoltà concettuali

I vincoli non olonomi non permettono, in generale, l'eliminazione di coordinate. Inoltre, non possono essere trattati semplicemente aggiungendo il vincolo al Lagrangiano, poiché il principio variazionale standard presuppone variazioni indipendenti $ \delta q_i $.

Equazioni di Lagrange con vincoli non olonomi

Il trattamento corretto si basa sul principio di d'Alembert--Lagrange. Le equazioni del moto si scrivono $$\large \frac{d}{dt}\frac{\partial L}{\partial \dot q_i}\frac{\partial L}{\partial q_i}= \sum_\alpha \lambda_\alpha a_{\alpha i}(q,t), $$ soggette ai vincoli $$\large \sum_i a_{\alpha i}(q,t),\dot q_i + a_{\alpha 0}(q,t) = 0. $$

I moltiplicatori $ \lambda_\alpha $ non derivano da un principio variazionale puro, ma rappresentano le forze vincolari necessarie a garantire il rispetto dei vincoli.

Principio di d'Alembert--Lagrange

Il principio fondamentale è $$\large \sum_i (Q_i - m_i \ddot q_i)\,\delta q_i = 0, $$ dove le variazioni $ \delta q_i $ sono compatibili con i vincoli imposti.

Riepilogo

from tabella import ddf; df = ddf(); df

	Tipo di vincolo	Dipendenza	Riduzione dei g.d.l.	Trattamento
0	Olonomo	q, t	Sì	Coordinate o moltiplicatori
1	Non olonomo	q̇	No (in generale)	d'Alembert–Lagrange
2	Integrabile	—	Sì	Diventa olonomo

Ruota volvente senza attrito in un piano - imposto moto quasi circolare

Dinamica di $Cetaev$

Contesto

Nel caso di vincoli olonomi, il principio variazionale standard di Hamilton è applicabile poiché le variazioni virtuali $ \delta q_i $ sono indipendenti (a meno dei vincoli sulle coordinate).

Nel caso di vincoli non olonomi, che agiscono sulle velocità, non tutte le variazioni virtuali sono ammissibili e il principio di Hamilton non è direttamente applicabile. La dinamica di Cetaev fornisce una regola precisa per selezionare tali variazioni.

Vincoli non olonomi lineari nelle velocità

Si considerino vincoli della forma \begin{equation} \sum_{i=1}^n a_{\alpha i}(q,t)\,\dot q_i + a_{\alpha 0}(q,t) = 0, \qquad \alpha = 1,\dots,m . \end{equation}

Questi vincoli sono tipici di problemi di rotolamento senza strisciamento e di contatto meccanico.

Condizione di Cetaev

La regola di Cetaev impone che le variazioni virtuali ammissibili soddisfino \begin{equation} \sum_{i=1}^n a_{\alpha i}(q,t)\,\delta q_i = 0, \qquad \alpha = 1,\dots,m . \end{equation}

Tale condizione si ottiene sostituendo formalmente $$ \dot q_i \longrightarrow \delta q_i $$ nel vincolo, eliminando i termini indipendenti dalle velocità.

Principio variazionale

Applicando il principio di d'Alembert--Lagrange con variazioni compatibili con la condizione di Cetaev, si richiede \begin{equation} \delta \int_{t_1}^{t_2} L(q,\dot q,t)\,dt = 0, \end{equation} soggetto ai vincoli sulle variazioni \begin{equation} \sum_i a_{\alpha i}(q,t)\,\delta q_i = 0 . \end{equation}

Introducendo moltiplicatori di Lagrange $ \lambda_\alpha $, si ottengono le equazioni del moto \begin{equation} \frac{d}{dt}\frac{\partial L}{\partial \dot q_i}\frac{\partial L}{\partial q_i}= \sum_{\alpha=1}^m \lambda_\alpha a_{\alpha i}(q,t) . \end{equation}

Significato fisico

Le forze vincolari associate ai moltiplicatori $ \lambda_\alpha $ non compiono lavoro virtuale sulle variazioni ammissibili. Il moto risultante rispetta i vincoli cinematici ed è coerente con l'evidenza fisica dei sistemi reali.

Osservazioni teoriche

La sostituzione formale $$\large \dot q_i \rightarrow \delta q_i $$ non deriva direttamente dal calcolo delle variazioni, ma costituisce un postulato fisico. Per questo motivo esistono formulazioni alternative, come la dinamica vakonomica o il formalismo di Gibbs--Appell, che conducono a equazioni del moto differenti.

Sintesi

La dinamica di Cetaev costituisce una formulazione della meccanica lagrangiana specifica per vincoli non olonomi lineari nelle velocità, basata su una scelta mirata delle variazioni virtuali e fondamentale per la corretta descrizione di sistemi con rotolamento o contatto.

Confronto tra approcci

from tabella import *; dg = ddg(); dg

	Approccio	Principio	Risultato
0	Cetaev	d'Alembert--Lagrange	Dinamica fisicamente corretta
1	Vakonomico	Hamilton con vincoli	Moto differente
2	Olonomo	Hamilton standard	Risultato unico

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Parameters
R = 1.0    # wheel radius
m = 1.0    # mass (translational)
I = 0.5    # moment of inertia (rotational)

def A_matrix(theta):
    """
    Constraint matrix A(q) (2x3) for constraints A dq = 0:
    [ 1, 0, -R*cos(theta) ]
    [ 0, 1, -R*sin(theta) ]
    """
    return np.array([
        [1.0, 0.0, -R * np.cos(theta)],
        [0.0, 1.0, -R * np.sin(theta)]
    ])

def dA_dt(theta, omega):
    """
    Time derivative of A: dA/dt = (dA/dtheta) * theta_dot
    dA/dtheta =
      [ 0, 0,  R*sin(theta) ]
      [ 0, 0, -R*cos(theta) ]
    So dA/dt = dA/dtheta * omega
    """
    return np.array([
        [0.0, 0.0,  R * np.sin(theta) * omega],
        [0.0, 0.0, -R * np.cos(theta) * omega]
    ])

def dynamics_chetaev(t, state):
    """
    State: [x, y, theta, vx, vy, omega]
    Returns derivative of state (first order):
      dx/dt = vx
      dy/dt = vy
      dtheta/dt = omega
      dv/dt = ddq (computed solving Chetaev linear system)
    """
    x, y, theta, vx, vy, omega = state
    qdot = np.array([vx, vy, omega])  # dq

    # Mass matrix M (3x3)
    M = np.diag([m, m, I])

    # Constraint matrix and its time derivative
    A = A_matrix(theta)           # 2x3
    dA = dA_dt(theta, omega)      # 2x3

    # Build the KKT block matrix:
    # [ M  -A^T ]
    # [ A   0  ]
    KKT_top = np.hstack((M, -A.T))   # 3 x (3+2)
    KKT_bot = np.hstack((A, np.zeros((2,2))))  # 2 x (3+2)
    KKT = np.vstack((KKT_top, KKT_bot))        # 5x5

    # RHS:
    # top RHS is zero (no external generalized forces)
    rhs_top = np.zeros(3)
    # bottom RHS is - dA/dt * qdot
    rhs_bot = - dA.dot(qdot)

    rhs = np.concatenate((rhs_top, rhs_bot))

    # Solve for [ddq; lambda]
    sol = np.linalg.solve(KKT, rhs)
    ddq = sol[:3]
    lambdas = sol[3:]

    # pack derivatives
    dxdt = vx
    dydt = vy
    dthetadt = omega
    dvxdt, dvydt, domegadt = ddq

    return np.array([dxdt, dydt, dthetadt, dvxdt, dvydt, domegadt])

# Initial conditions: choose omega nonzero and set vx,vy consistent with constraints
theta0 = 0.3  # radians initial orientation
omega0 = 1.0  # initial angular speed
vx0 = R * omega0 * np.cos(theta0)
vy0 = R * omega0 * np.sin(theta0)
state0 = np.array([0.0, 0.0, theta0, vx0, vy0, omega0])

# Integrate
t_span = (0.0, 10.0)
t_eval = np.linspace(t_span[0], t_span[1], 400)
sol = solve_ivp(dynamics_chetaev, t_span, state0, t_eval=t_eval, rtol=1e-9, atol=1e-9)

# Extract
t = sol.t
x = sol.y[0]
y = sol.y[1]
theta = sol.y[2]
vx = sol.y[3]
vy = sol.y[4]
omega = sol.y[5]

# Compute constraint residuals over time: should be (near) zero
c1 = vx - R * omega * np.cos(theta)
c2 = vy - R * omega * np.sin(theta)

# Plot 1: Trajectory x-y
plt.figure(figsize=(6,5))
plt.plot(x, y)
plt.xlabel("x")
plt.ylabel("y")
plt.title("Traiettoria della ruota (metodo di Chetaev)")
plt.grid(True)
plt.show()

# Plot 2: Orientation theta vs time
plt.figure(figsize=(6,3))
plt.plot(t, theta)
plt.xlabel("t")
plt.ylabel("theta (rad)")
plt.title("Orientazione theta(t)")
plt.grid(True)
plt.show()

# Plot 3: Velocities vx, vy, omega vs time (each in its own figure per rules)
plt.figure(figsize=(6,3))
plt.plot(t, vx)
plt.xlabel("t")
plt.ylabel("vx")
plt.title("vx(t)")
plt.grid(True)
plt.show()

plt.figure(figsize=(6,3))
plt.plot(t, vy)
plt.xlabel("t")
plt.ylabel("vy")
plt.title("vy(t)")
plt.grid(True)
plt.show()

plt.figure(figsize=(6,3))
plt.plot(t, omega)
plt.xlabel("t")
plt.ylabel("omega")
plt.title("omega(t)")
plt.grid(True)
plt.show()

# Plot 4: Constraint residuals
plt.figure(figsize=(6,3))
plt.plot(t, c1)
plt.xlabel("t")
plt.ylabel("c1 (vx - R*omega*cosθ)")
plt.title("Residuo vincolo c1 nel tempo")
plt.grid(True)
plt.show()

plt.figure(figsize=(6,3))
plt.plot(t, c2)
plt.xlabel("t")
plt.ylabel("c2 (vy - R*omega*sinθ)")
plt.title("Residuo vincolo c2 nel tempo")
plt.grid(True)
plt.show()

# Print quick diagnostics (first/last values)
print("Initial state:", state0)
print("Final state (t=%.2f):" % t[-1], sol.y[:, -1])
print("Max absolute constraint residuals:", np.max(np.abs(c1)), np.max(np.abs(c2)))

Initial state: [0.         0.         0.3        0.95533649 0.29552021 1.        ]
Final state (t=10.00): [-1.06320602  1.59616291 10.3        -0.64082642 -0.76768581  1.        ]
Max absolute constraint residuals: 2.630888451538027e-09 1.8532565482232144e-09

Ruota volvente senza attrito in un piano - imposto moto quasi circolare, ma con frizione nella rotazione

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Parameters
R = 1.0     # wheel radius
m = 1.0     # mass
I = 0.5     # inertia
k_fric = 2.0   # rotational friction coefficient (tunable)

def A_matrix(theta):
    return np.array([
        [1.0, 0.0, -R * np.cos(theta)],
        [0.0, 1.0, -R * np.sin(theta)]
    ])

def dA_dt(theta, omega):
    return np.array([
        [0.0, 0.0,  R*np.sin(theta)*omega],
        [0.0, 0.0, -R*np.cos(theta)*omega]
    ])

def dynamics_chetaev_friction(t, state):
    x, y, theta, vx, vy, omega = state
    qdot = np.array([vx, vy, omega])

    M = np.diag([m, m, I])
    A = A_matrix(theta)
    dA = dA_dt(theta, omega)

    # friction torque τ = -k_fric * ω
    Q = np.array([0.0, 0.0, -k_fric * omega])

    KKT_top = np.hstack((M, -A.T))
    KKT_bot = np.hstack((A, np.zeros((2,2))))
    KKT = np.vstack((KKT_top, KKT_bot))

    rhs_top = Q
    rhs_bot = - dA.dot(qdot)
    rhs = np.concatenate((rhs_top, rhs_bot))

    sol = np.linalg.solve(KKT, rhs)
    ddq = sol[:3]

    dxdt = vx
    dydt = vy
    dthetadt = omega
    dvxdt, dvydt, domegadt = ddq

    return np.array([dxdt, dydt, dthetadt, dvxdt, dvydt, domegadt])


theta0 = 0.3
omega0 = 2.0   # faster to see friction effect
vx0 = R * omega0 * np.cos(theta0)
vy0 = R * omega0 * np.sin(theta0)
state0 = np.array([0.0, 0.0, theta0, vx0, vy0, omega0])

t_span = (0, 10)
t_eval = np.linspace(0, 10, 400)

sol = solve_ivp(dynamics_chetaev_friction, t_span, state0, t_eval=t_eval, rtol=1e-9, atol=1e-9)

t = sol.t
x, y, theta, vx, vy, omega = sol.y

c1 = vx - R * omega * np.cos(theta)
c2 = vy - R * omega * np.sin(theta)

plt.figure(figsize=(6,5))
plt.plot(x, y)
plt.xlabel("x")
plt.ylabel("y")
plt.title("Traiettoria con attrito rotazionale (Chetaev)")
plt.grid(True)
plt.show()

plt.figure(figsize=(6,3))
plt.plot(t, omega)
plt.xlabel("t")
plt.ylabel("omega")
plt.title("Velocità angolare ω(t) con attrito")
plt.grid(True)
plt.show()

plt.figure(figsize=(6,3))
plt.plot(t, theta)
plt.xlabel("t")
plt.ylabel("theta")
plt.title("Orientazione θ(t)")
plt.grid(True)
plt.show()

plt.figure(figsize=(6,3))
plt.plot(t, c1)
plt.xlabel("t")
plt.ylabel("c1")
plt.title("Residuo c1")
plt.grid(True)
plt.show()

plt.figure(figsize=(6,3))
plt.plot(t, c2)
plt.xlabel("t")
plt.ylabel("c2")
plt.title("Residuo c2")
plt.grid(True)
plt.show()