Puma560 krok-po-kroku: warunek Piepera, środek nadgarstka, q₁…q₆ z geometrii i algebry.
Chcemy znaleźć kąty wszystkich sześciu przegubów Pumy560 tak, aby końcówka narzędzia (TCP) znalazła się w zadanym punkcie z zadaną orientacją. Zrobimy to analitycznie — czyli podamy jawne wzory algebraiczne, bez iteracji i bez szukania numerycznego. Dla Pumy jest to wyjątkowo proste dzięki szczególnej własności geometrycznej: osie trzech ostatnich przegubów (, , ) przecinają się w jednym punkcie — środku nadgarstka. Ta własność jest jedną z form warunku Piepera (Pieper 1968).
Warunek Piepera — wystarczający, nie konieczny
W literaturze (i często w skryptach uczelnianych) pada teza, że „rozwiązanie analityczne istnieje wtedy i tylko wtedy, gdy spełniony jest warunek Piepera". To nieprawda. Pieper pokazał jedynie, że jeśli trzy kolejne osie (a) przecinają się w jednym punkcie albo (b) są wzajemnie równoległe, to istnieje rozwiązanie zamknięte z dekompozycji 3+3.
Kontrprzykłady:
Praktycznie każdy 6-DOF stosowany dziś w przemyśle (Puma, Stäubli, KUKA, ABB, Fanuc, UR) ma analityczne IK — bo producenci celowo projektują geometrię, by upraszczała wyprowadzenie. Dlatego analityczne IK nie jest egzotyką ani ograniczone do archaicznych konstrukcji.
Wracając do Pumy: warunek Piepera (forma A) sprawia, że 6-wymiarowe zadanie rozpada się na dwa łatwiejsze podproblemy 3-wymiarowe — najpierw wyznaczamy , żeby środek nadgarstka trafił w odpowiednie miejsce w przestrzeni, a potem , żeby narzędzie miało żądaną orientację. Cały rachunek sprowadza się do kilku zastosowań twierdzenia cosinusów i funkcji atan2.

Pomocna metafora — Puma jako dwa „podroboty" zmontowane szeregowo
Zamiast od razu liczyć, warto na minutę spojrzeć na manipulator jako na dwa logicznie niezależne urządzenia:
To rozdzielenie nie jest tylko opisem — jest matematycznym faktem: pochodna pozycji wrist centre względem jest zerowa, bo obroty wokół osi przechodzących przez punkt nie przesuwają tego punktu. Dlatego można najpierw rozwiązać podproblem pozycji (3 niewiadome, 3 równania skalarne), a dopiero potem podproblem orientacji (kolejne 3 niewiadome z warunku, że pełna macierz rotacji ma być taka jak zadano). Bez tego rozdzielenia musielibyśmy jednocześnie zmagać się z układem 6 równań nieliniowych, dla którego analityczne rozwiązanie istnieje, ale jest dramatycznie bardziej skomplikowane (równanie 16-tego stopnia w metodzie Raghavana–Rotha dla geometrii nie-pieperowskich).
Do opisu manipulatora używamy zmodyfikowanej konwencji Denavita–Hartenberga (Craig, Introduction to Robotics, wyd. 3, rozdz. 3). Każdemu ogniwu przypisujemy własny układ współrzędnych (w literaturze angielskiej: frame — dalej krócej „układ {i}"). Przypisanie nie jest dowolne: oś musi pokrywać się z osią obrotu przegubu i, a oś leży wzdłuż wspólnej normalnej łączącej osie przegubów i−1 i i.
Cztery parametry na ogniwo — , , , — w pełni wyznaczają transformację z układu {i−1} do {i}. Intuicja geometryczna każdego z nich:
Odległość między osiami przegubów i−1 oraz i, mierzona wzdłuż ich wspólnej normalnej.
Kąt obrotu osi ẑi−1 do ẑi, obracany wokół wspólnej normalnej.
Przesunięcie początku układu {i} wzdłuż osi ẑi od wspólnej normalnej.
Kąt obrotu osi x̂i−1 do x̂i wokół ẑi. Dla przegubu obrotowego jest to zmienna konfiguracji.
Transformację z układu {i−1} do {i} można rozbić na cztery kolejne operacje wzdłuż lub wokół pojedynczej osi:. Każdy krok wprowadza jeden pośredni układ pomocniczy:
Złożenie tych czterech operacji daje pełną transformację 4×4. W praktyce piszemy od razu macierz w postaci zamkniętej (niżej), ale rozbicie na cztery kroki bywa pomocne przy derywacji i debugowaniu.
W literaturze funkcjonują dwie różne wersje konwencji Denavita–Hartenberga. Różnią się miejscem przypinania układu współrzędnych do ogniwa oraz kolejnością elementarnych transformacji. Liczbowo te same parametry, ale inny indeks ogniwa — to najczęstsza pułapka przy porównywaniu wzorów z różnych książek.


| Aspekt | Klasyczna DH (1955) | Zmodyfikowana DH Craiga (1986) |
|---|---|---|
| Pozycja układu {i} | Na końcu ogniwa i (po przegubie i+1) | Na początku ogniwa i (oś = oś przegubu i) |
| Oś | Pokrywa się z osią przegubu i+1 | Pokrywa się z osią przegubu i |
| Transformacja | ||
| Indeksy parametrów | Wszystkie mają indeks i () | Mieszane: i |
| Kolejność operacji | Najpierw rotacja wokół () | Najpierw rotacja wokół () |
| Intuicja indeksu | Parametry „opisują" ogniwo i−1→i jako całość | charakteryzują ogniwo stałe, — przegub zmienny |
| Obsługa manipulatorów rozgałęzionych | Trudna (indeks i zbiera parametry dwóch różnych ogniw) | Naturalna (parametry jednoznacznie przypisane do linku) |
Jak to wpływa na IK? Dla tego samego fizycznego robota liczby są takie same, ale ich przypisanie do indeksów różni się o jeden. Jeśli wyprowadzasz rozwiązanie analityczne, pilnuj której książki używasz — wzory Craiga (które stosujemy w tym module) nie pasują dosłownie do tabeli DH z Siciliano/Khatib bez reindeksacji.
W tym module i w całej aplikacji używamy zmodyfikowanej konwencji Craiga — głównie dlatego, że dla Pumy560 ma ona najkrótsze i najczytelniejsze wzory (tabela α_{i-1} kończy się zerami i , co od razu upraszcza iloczyny macierzy). Wybór ma jednak charakter notacyjny — to nie jest „jedyna poprawna" konwencja.
| 1 | ||||
| 2 | ||||
| 3 | ||||
| 4 | ||||
| 5 | ||||
| 6 |
to długość ramienia (upper arm), to długość przedramienia (forearm). Dodatkowo w Pumie występują dwa niezerowe odsadzenia: (odsadzenie przedramienia od osi ramienia, ok. 12,5 cm) oraz (mały bok łokcia, ok. 2 cm). Te odsadzenia są źródłem większości komplikacji w wyprowadzeniu — bez nich Puma sprowadzałaby się do klasycznego planarnego manipulatora 2R.
Transformacja ogniwa zapisana w konwencji zmodyfikowanej:
gdzie , itd. Mnożąc kolejne takie macierze, dostaniemy — pełne przekształcenie z bazy do układu efektora.
Zanim zanurzymy się w algebrę, ułóżmy w głowie obraz tego, jak ramię „pracuje". Każdy przegub Pumy ma jedną dokładnie wyznaczoną funkcję mechaniczną — i dla każdej istnieje codzienna analogia, która ułatwia zapamiętanie. To krótkie wprowadzenie zwraca się z nawiązką, gdy później pojawi się pytanie „dlaczego kontroluje wysokość, a nie".
Decyduje wyłącznie o azymucie — w którą stronę świata zwrócone jest ramię. Nie zmienia ani wysokości, ani wysunięcia TCP, jedynie obraca całą resztę manipulatora. Z perspektywy IK: ten staw widzimy w „rzucie z góry" — bierzemy docelowej pozycji i pytamy: w którą stronę ramię musi się skręcić, żeby w ogóle „móc patrzeć" na cel.
Wraz z kontroluje punkt 2D, w który sięga końcówka ramienia (wysokość + wysunięcie radialne) — w płaszczyźnie pionowej wyznaczonej już przez . Czysto „pionowy" stopień swobody barku.
Decyduje o „rozłożeniu" całego ramienia: im bardziej wyciągnięte, tym dalej sięga; im bardziej zgięte, tym wyżej/bliżej. Geometrycznie odpowiada za długość przekątnej bark→wrist center (poprzez prawo cosinusów — patrz krok 4). Tutaj pojawia się słynna dwuznaczność „elbow up / elbow down".
Pierwszy z trzech stawów nadgarstka. Nie zmieniapozycji środka nadgarstka, bo oś obrotu przechodzi dokładnie przez ten punkt. Wpływa za to na orientację — to właśnie ten staw „skręca" całą głowę nadgarstka wokół przedramienia.
Kąt ataku narzędzia względem osi przedramienia. Gdy , narzędzie jest „wyciągnięte" w jednej linii z przedramieniem — wówczas występuje gimbal lock: osie i pokrywają się i tracą niezależność.
Końcowy stopień swobody — orientuje narzędzie wokół jego własnej osi po tym, jak ustawiły kierunek, w który narzędzie „patrzy". Bez robot mógłby tylko wskazywać kierunek; z nim — może też ustawić orientację gripper'a, szczęk czy dyszy.
Dlaczego trzy + trzy? Trzy pierwsze stawy () ustalają pozycję 3-wymiarowego punktu (środek nadgarstka) — bo 3 niewiadome wystarczają, żeby opisać dowolny punkt w . Trzy ostatnie () ustalają orientację w przestrzeni — bo grupa obrotów jest również 3-wymiarowa. To nie przypadek — to inżynierski projekt optymalizujący liczbę stopni swobody pod konkretne zadanie: ustawić narzędzie w pełnej dowolności w przy możliwie najprostszej geometrii.
Po lewej: manipulator, którym sterujesz ręcznie sześcioma suwakami — służy do generowania dowolnych poz testowych. Po prawej: pole pozy docelowej z przyciskiem „zrzut z kontrolera", który kopiuje aktualne robota jako cel IK. Wszystkie wartości pośrednie w kolejnych krokach przeliczane są na bieżąco — zmiana pozy celu od razu aktualizuje liczby poniżej.
Zadana poza efektora to macierz transformacji jednorodnej — górne to orientacja (macierz rotacji ), a ostatnia kolumna to pozycja :
Zbiór to zbiór wszystkich takich macierzy — grupa sztywnych przesunięć i obrotów w przestrzeni. Szukamy wektora kątów spełniającego , gdzie to odwzorowanie FK (kinematyka prosta).
Obserwacja kluczowa: obroty w przegubach 4, 5, 6 dzieją się wokół osi przecinających się w jednym punkcie (środku nadgarstka ). Obrót wokół osi przechodzącej przez punkt nie przesuwa tego punktu. Wobec tego położenie nie zależy od — zależy tylko od .
W konwencji Craiga dla Pumy zachodzi , więc środek nadgarstka pokrywa się z początkiem układu {6}. Jeśli dodatkowo nie ma offsetu narzędzia, środek nadgarstka to po prostu pozycja efektora:
Gdy jest narzędzie (pewien stały offset między układem {6} a układem TCP), cofamy się o ten offset:
Reasumując plan: najpierw znajdziemy z pozycji środka nadgarstka (3 równania, 3 niewiadome), potem z warunku, że orientacja pozostałych trzech przegubów ma domknąć zadane .
const T06 = options.toolOffset
? mul4(target, invSE3(options.toolOffset))
: target;
const R = extractRotation(T06); // górne 3×3
const [px, py, pz] = extractPosition(T06); // ostatnia kolumnaPozycja początku układu {4} (czyli środka nadgarstka) w bazie to ostatnia kolumna iloczynu czterech macierzy . Wstawiając parametry z tabeli DH (patrz wyprowadzenie poniżej), otrzymujemy:
gdzie , , , .
Każda w konwencji Craiga ma postać podaną wyżej. Dla Pumy 560 podstawiamy z tabeli DH:
(, więc tylko obrót wokół z).
().
().
Przy szukaniu pozycji środka nadgarstka (pierwsze 3 elementy ostatniej kolumny ) nie potrzebujemy : wystarczy „adres" początku układu {4}, czyli punkt wyrażony w bazie. Liczymy więc krócej:
gdzie ostatnia transformacja sprowadza się do obrotu wokół x i przesunięcia o w x i w z.
Wymnażając kolejno: najpierw daje przesunięcie łokcia o obrócone przez , potem obraca to wokół osi x o i dodaje przegub , finalnie obraca całość wokół pionowej osi bazy. W rezultacie dostajemy podane wyżej — z efektywną długością przedramienia rozłożoną na składowe w „radialnym" kierunku oraz w pionie.
Zauważmy, że w pierwszych dwóch równaniach powtarza się to samo wyrażenie w nawiasie. Wprowadźmy skrót:
Geometrycznie to radialna odległość środka nadgarstka od pionowej osi bazy w obróconym o układzie — innymi słowy, „jak daleko" wrist centre jest od pionowej osi obrotu barku.
Z dwóch pierwszych równań mamy i . Są to składowe obrotu o kąt wektora w płaszczyźnie XY:
Mnożąc obie strony tej równości przez (transpozycja = inwersja, bo macierz obrotu jest ortogonalna) i odczytując:
Podnosząc oba równania do kwadratu i dodając, lewa strona daje , a prawa = po rozwinięciu i uproszczeniu (krzyżowe wyrazy się znoszą, ) . Stąd kluczowa identyczność:
Ta równość jest fundamentem całej dekompozycji. Lewa strona zależy wyłącznie od zadanych (wiemy je), prawa strona zawiera tylko . Dlatego wyliczamy natychmiast — bez znajomości .
# T_target: macierz 4x4 (poza efektora), np. NumPy array.
# R = T_target[:3, :3], p = T_target[:3, 3]
import numpy as np
R = T_target[:3, :3]
p = T_target[:3, 3]
# Środek nadgarstka: cofamy się o d6 wzdłuż osi z układu 6.
# Dla Pumy d6 = 0, więc p_wc == p, ale piszemy ogólnie:
p_wc = p - D6 * R[:, 2]
px_w, py_w, pz_w = p_wcIntuicja: cień TCP na podłodze
Wyobraź sobie, że nad TCP świeci pionowa lampa rzucająca cień na poziomą płaszczyznę podstawy. Ten cień ma w bazie współrzędne . Naturalna pierwsza myśl: obrót podstawy musi po prostu skierować ramię w stronę cienia — czyli .
Dla manipulatora bez odsadzenia tak właśnie by było. Puma ma jednak — boczne odsunięcie osi przedramienia od osi ramienia. Skutkiem tego, gdy bark obraca się o kąt , przedramię (a w konsekwencji środek nadgarstka) nie znajduje się ściśle przed barkiem — leży lekko obok, w bok od osi ramienia o stałą odległość .
Z perspektywy „cienia na podłodze" oznacza to, że cień nie leży dokładnie na linii barku — jest przesunięty w bok. Stąd potrzeba drugiego członu: . Pierwszy atan2 wskazuje na cień, drugi koryguje o ten boczny offset (dla typowych pozycji ~7°). Te 7° to dokładnie błąd, jaki popełniają studenckie kody, w których pominięto drugi człon.
Geometrycznie: wektor () w obróconym układzie ramienia to ten sam wektor co () w bazie — tylko obrócony o . Wektor () ma nachylenie względem osi ramienia, a wektor () ma nachylenie względem osi X bazy. Różnica tych dwóch nachyleń to właśnie kąt obrotu między układami, czyli .
Z identyczności kroku 1 wyznaczamy , ale z dokładnością do znaku (bo pierwiastek kwadratowy daje dwie wartości):
Jeśli , zadanie jest nieosiągalne — cel leży w „zakazanym" cylindrze o promieniu wokół pionowej osi bazy. W pozostałych przypadkach znak wyznacza, czy bark Pumy zwraca się w kierunku celu (shoulder right, ), czy odwrotnie, o obrócony (shoulder left, ). Obie konfiguracje dają identyczną pozycję wrist centre — stąd „gałąź" w algorytmie.
Mając i , rozwiązujemy układ dwóch liniowych równań w z poprzedniej strony:
Formalnie rozwiązujemy go przez wyznacznik, ale ładniej zapisuje się to przez różnicę dwóch funkcji atan2 — jest różnicą kąta do celu () i kąta „w bok" o odsadzeniu :
# ρ — promień rzutu wrist centre na płaszczyznę xy bazy.
# Dwa znaki dają shoulder LEFT i RIGHT.
r_xy_sq = px_w**2 + py_w**2
disc_q1 = r_xy_sq - D3**2
if disc_q1 < 0:
return [] # cel poza "cylindrem zakazanym"
rho_abs = np.sqrt(max(0.0, disc_q1))
phi = np.arctan2(py_w, px_w)
solutions = []
for shoulder_sign in (+1, -1): # +1 = right, -1 = left
rho = shoulder_sign * rho_abs
q1 = phi - np.arctan2(D3, rho)
# ... dalej q3, q2, q4-q6 dla tej gałęziDwuargumentowa funkcja atan2(y, x) — w odróżnieniu od zwykłego arctan(y/x) — zwraca kąt na pełnym przedziale , uwzględniając znaki obu argumentów. To ważne: zwykły arctan traci informację o ćwiartce, przez co „widzi" kąty z góry i z dołu jako identyczne.
Częsty błąd w kodach studenckich: użycie atan2(py, px) bez drugiego członu atan2(d₃, ρ). Dla Pumy daje to stały błąd rzędu 7° — odpowiadający kątowi „w bok" wymuszonemu przez odsadzenie .
const phi = Math.atan2(py, px); // kąt do celu w XY
const rho = rhoSign * Math.sqrt(px*px + py*py - D3*D3); // ±
const q1 = phi - Math.atan2(D3, rho); // q₁ na obu gałęziachPo wyznaczeniu obracamy cały układ o ten kąt wokół pionowej osi bazy. W tak obróconym układzie całe ramię robota leży w jednej pionowej płaszczyźnie — oznaczmy jej współrzędne jako , gdzie to odległość od pionowej osi bazy (dodatnia po „stronie" celu), a — wysokość nad płaszczyzną podstawy.
W tej płaszczyźnie zostało nam zadanie 2-wymiarowe: połączyć bark ze środkiem nadgarstka za pomocą dwóch ogniw — ramienia o długości i przedramienia o pewnej efektywnej długości , z łokciem gdzieś pomiędzy.
Czym jest ? Początek układu {4} (wrist centre) nie leży na prostym przedłużeniu ramienia — jest przesunięty względem łokcia o dwa składowe: wzdłuż osi X układu {3} (mały bok, ok. 2 cm) i wzdłuż osi Z układu {3} (długi, ok. 43 cm). Skąd ta konstrukcja? Fizycznie wynika z tego, że silnik napędzający przedramię musi się gdzieś zmieścić obok łokcia, więc geometrycznie przedramię jest „schowane" za mały boczny segment.
W płaszczyźnie ramienia wektor od łokcia do środka nadgarstka ma długość:
i jest odchylony od osi X układu {3} o stały kąt:
O tym odchyleniu pamiętamy przy obliczaniu w następnym kroku.
Diagram przedstawia planarny podproblem 2R: bark (niebieski) i środek nadgarstka (czerwony) leżą w jednej płaszczyźnie pionowej, obróconej o q₁ wokół osi bazy. Dwie gałęzie — zielona (elbow up) i pomarańczowa (elbow down) — pokazują alternatywne położenia łokcia (pełna linia = ramię a₂, linia przerywana = efektywne przedramię L). Promień D = odległość bark–wrist (fioletowa przerywana) wchodzi do prawa kosinusów w kroku 4.
Dlaczego q₃ daje się wyznaczyć osobno — niezmiennik odległości
Kluczowa obserwacja, która sprawia, że można w ogóle w tym momencie szukać w izolacji od reszty: odległość |bark → wrist center| nie zależy od ani — zależy wyłącznie od .
Intuicyjnie: i obracają cały trójkąt bark-łokieć-wrist jako sztywne ciało — pierwszy wokół pionowej osi bazy, drugi wokół osi barku. Obroty sztywne nie zmieniają długości żadnego boku trójkąta. To, co zmienia wewnętrzny kształt trójkąta — a w szczególności długość przekątnej bark↔wrist — to wyłącznie zgięcie w łokciu, czyli .
Stąd ważny wniosek metodologiczny: mając zadaną pozycję wrist centre (a więc znaną przekątną ), możemy natychmiast wyłuskać z jednego równania skalarnego, nie wiedząc jeszcze ile wynosi . Dopiero potem wrócimy do dwuwymiarowego problemu i wyznaczymy , traktując kształt trójkąta jako już ustalony. To dlatego algorytm liczy stawy w kolejności , a nie w kolejności indeksów.
Spojrzenie geometryczne. W płaszczyźnie ramienia bark , łokieć i środek nadgarstka tworzą trójkąt z trzema znanymi długościami: (ramię), (efektywne przedramię) i (przekątna bark→wrist). Każdy trójkąt o danych trzech bokach ma jednoznacznie wyznaczone wszystkie kąty wewnętrzne — w szczególności kąt przy łokciu:
Trzy znane długości — a₂ (ramię), L (efektywne przedramię, L = √(a₃² + d₄²)), D (odległość bark↔nadgarstek, D = √(ρ² + p_z²)) — domykają trójkąt. Wewnętrzny kąt γ przy łokciu wynika wprost z prawa cosinusów. Drugie rozwiązanie (elbow ↑ vs ↓) odpowiada odbiciu trójkąta względem boku D — znak sin γ przy rekonstrukcji γ = atan2(±sin γ, cos γ).
Obraz pomocniczy — naprężona linka i dwa lustrzane trójkąty
Wyobraź sobie naprężoną linkę o długości , której końce są przymocowane w barku i w środku nadgarstka . Łokieć jest swobodnym przegubem o ustalonej odległości od barku i od wrist — może „opaść" po dowolnej stronie linki. Daje to dwa identycznie wymierzone trójkąty będące względem siebie lustrzanym odbiciem przez bok :
Obie konfiguracje dają dokładnie tę samą pozycję wrist centre, ale wymagają różnych . Stąd binarny wybór znaku w wyrażeniu na — to przełącznik między dwoma lustrzanymi obrazami tego samego trójkąta.
Z prawa cosinusów (zastosowanego do boku przeciwległego kątowi ):
Pełny kąt wewnętrzny trójkąta wyciąga się standardowo przez atan2 (a nie acos), żeby zachować informację o znaku :
Znak przed odpowiada wyborowi elbow up (łokieć powyżej linii bark↔wrist) lub elbow down (poniżej) — dwa odbicia trójkąta względem boku .
Powiązanie z q₃. Kąt przy łokciu nie jest tożsamy z . Gdyby przedramię biegło dokładnie wzdłuż osi (czyli ), wówczas . W Pumie jednak przedramię jest „odchylone" od osi o stały kąt (z kroku 3), więc:
Biorąc dwie zależności z kroku 1, i , podnosimy oba do kwadratu i dodajemy. Wyrazy z i składają się w sumę kątów , podobnie . Po uproszczeniu:
Oznaczając lewą stronę (znaną) przez :
Wyrażenie to liniowa kombinacja sinusa i cosinusa o tym samym argumencie — można ją zapisać jako jeden cosinus z przesuniętym kątem:
(sprawdź: rozwijając prawą stronę przez i podstawiając , — obie strony się zgadzają).
Stąd , więc , i:
Spójność z wyprowadzeniem geometrycznym: (różnica znaku wynika z definicji kąta), więc , czyli — ten sam wynik.
Postać końcowa (wykorzystywana w kodzie — jeden atan2 zamiast łańcucha kątów):
Znak odpowiada dwóm gałęziom geometrycznym — łokieć „w górę" lub „w dół". Gdy , cel leży poza zasięgiem ramienia dla danej gałęzi barku — kombinacja shoulder+elbow nie ma rozwiązania (trójkąt o takich bokach nie istnieje, bo nie spełnia nierówności trójkąta).
# Wewnątrz pętli po shoulder_sign:
K = (rho**2 + pz_w**2 - A2**2 - A3**2 - D4**2) / (2 * A2)
L = np.sqrt(A3**2 + D4**2)
beta = np.arctan2(D4, A3)
disc = L**2 - K**2
if disc < 0:
continue # cel poza zasięgiem
sqrt_d = np.sqrt(disc)
for elbow_sign in (+1, -1): # +1 = up, -1 = down
q3 = np.arctan2(elbow_sign * sqrt_d, K) - beta
# ... dalej q2, q4-q6 dla tej gałęziMając już znane , wracamy do dwóch równań pozycji z kroku 1 — tym razem traktując jako jedyną niewiadomą. Składowe radialna i pionowa to:
Stosujemy wzory na sumę kątów: oraz . Wstawiając do pierwszego równania:
Grupujemy wyrazy przy i osobno:
Tak samo dla : . Wprowadzając skróty na powtarzające się grupy:
(zauważ: zależą tylko od znanego i stałych DH), oba równania zapisują się zwarto jako układ liniowy 2×2:
Wyznacznik macierzy współczynników: . Z reguły Cramera (lub bezpośrednio przez inwersję 2×2):
Wyznacznik znika tylko, gdy jednocześnie — niemożliwe dla fizycznego manipulatora. Mając i jednocześnie, kąt wyznacza atan2:
# Wewnątrz pętli po elbow_sign (q3 znane):
c3, s3 = np.cos(q3), np.sin(q3)
M = A2 + A3 * c3 - D4 * s3
N = A3 * s3 + D4 * c3
denom = M**2 + N**2
c2 = (M * rho - N * pz_w) / denom
s2 = (-M * pz_w - N * rho) / denom
q2 = np.arctan2(s2, c2)Dlaczego nie q₂ = arcsin(s₂) lub arccos(c₂)? Każda z tych funkcji zwraca wartość tylko z połowy okręgu — tracimy informację o znaku drugiego składnika. atan2(s, c) używa znaków obu argumentów, dając jednoznaczny kąt na całym okręgu . Klasyczny błąd implementacyjny: „policz arctan(s/c) i zgadnij ćwiartkę po znakach" — działa tylko po stronie poprawnie pamiętających autorów.
Do tego momentu zajmowaliśmy się pozycją środka nadgarstka. Teraz przechodzimy do orientacji: trzy ostatnie przeguby muszą obrócić nadgarstek tak, aby pełna macierz zgodziła się z zadanym .
Część rotacyjną macierzy — czyli — wyliczamy ze znanych już . Kluczowa obserwacja: w konwencji Craiga i obracają się wokół tej samej (poziomej) osi , bo w tabeli DH (między układami {2} i {3} nie ma skręcenia osi). Obroty wokół wspólnej osi zwyczajnie się dodają: . Razem:
Krok 1 — pomnóżmy :
(mnożenie wiersz·kolumna; ma , stąd druga i trzecia kolumna).
Krok 2 — pomnóżmy z wynikiem powyżej:
Pierwszy wiersz: .
Po uproszczeniu wszystkich trzech wierszy:
Sprawdzenie ortogonalności: każdy wiersz ma normę , a iloczyn skalarny pierwszych dwóch wierszy: . Zgodnie z oczekiwaniem dla macierzy obrotu.
Wynik jawny (przepisany dla wygody — używamy go bezpośrednio w kodzie):
Dalej: jeśli , to nadgarstek musi dostarczyć „brakującą" rotację:
Równość wynika z ortogonalności macierzy rotacji (wiersze są ortonormalne). W implementacji liczymy więc po prostu iloczyn transponowanej macierzy przez — bez wywoływania numerycznej inwersji.
# Wewnątrz pętli po elbow (q2, q3 znane):
c1, s1 = np.cos(q1), np.sin(q1)
c23 = np.cos(q2 + q3)
s23 = np.sin(q2 + q3)
R03 = np.array([
[ c1*c23, -c1*s23, -s1],
[ s1*c23, -s1*s23, c1],
[-s23, -c23, 0],
])
R36 = R03.T @ R # ortogonalność: T zamiast invTrzy ostatnie ogniwa Pumy to klasyczny nadgarstek typu ZYZ: dwa skręcenia osi między kolejnymi przegubami sprawiają, że osie , , są wzajemnie prostopadłe i schodzą się w jednym punkcie (środek nadgarstka). Z tabeli DH (, , ), iloczyn:
Skąd dokładnie te pięć macierzy? — rozbicie iloczynu na cegiełki DH
Część rotacyjna każdej macierzy DH (Craig) ma postać — najpierw skręcamy poprzednią oś o stały kąt (geometria ogniwa), potem obracamy o zmienną wokół nowej osi (kąt przegubu). Sklejając trzy ostatnie ogniwa, dostajemy:
Pierwsza cegiełka to czyste bo skręcenie zostało już „pochłonięte" w macierzy . Pozostałe dwie pary biorą skręcenia z tabeli DH — (oś obraca się o 90° aby pokryć się z osią prostopadłą do przedramienia) oraz (oś wraca o 90° w przeciwną stronę, aby oś wskazywała wzdłuż narzędzia). Stąd „kanapka" z–x–z–x–z: trzy obroty przegubowe przeplatane dwoma stałymi „skręceniami montażowymi" .
Dlaczego to się nazywa „typ ZYZ"? Bo gdy odpowiednio przenieść stałe macierze przez obroty (korzystając z tożsamości ), cała kanapka redukuje się do — klasycznych kątów Eulera ZYZ. Konwencja DH wymusza zapis przez i , ale geometria pod spodem to dokładnie te same trzy obroty z–y–z, których uczy się w kinematyce ciała sztywnego.
ma klasyczną postać macierzy ZYZ: środkowy obraca wokół osi prostopadłej do osi i , więc środkowy element całej macierzy to , a środkowy wiersz/kolumna kodują . Po wymnożeniu:
Element środkowego wiersza daje nam z dokładnością do znaku . Pozostałe elementy średniej kolumny i wiersza rozszyfrowują i :
Ponieważ występuje w każdej równości osobno, a może być dodatni lub ujemny, dostajemy dwie gałęzie nadgarstka (no-flip i flip). Dla :
Skąd ten pierwiastek? — jedynka trygonometryczna w działaniu
Mamy już — to wystarczyłoby do policzenia przez acos, ale dwa istotne powody każą sięgnąć po atan2: (1) tracimy informację o znaku sinusa, więc nie odróżnimy gałęzi no-flip od flip; (2) acos jest źle uwarunkowany w pobliżu (małe błędy zaokrąglenia → duże błędy kąta). Dlatego rekonstruujemy także .
Z dwóch równości z pierwszej kolumny środkowego wiersza:
Podnosimy do kwadratu i dodajemy, korzystając z jedynki trygonometrycznej :
Stąd — znak jest fundamentalnie nierozstrzygalny z tych równań (układ trzech równań na trzy niewiadome ma dwie gałęzie rozwiązań). Mając parę , składamy:
Bonus geometryczny: w nadgarstku typu ZYZ to długość rzutu osi narzędzia na płaszczyznę prostopadłą do . Gdy ta długość maleje do zera, osie i stają się współliniowe — pojawia się osobliwość nadgarstka opisana niżej.
(górne znaki — gałąź no-flip, dolne — flip).
Gdy , dwa końcowe obroty nadgarstka zachodzą wokół tej samej osi. Matematycznie wyzerowuje się cała średnia kolumna i średni wiersz (poza ), a w macierzy pozostaje tylko rotacja sumaryczna . Tracimy stopień swobody: pojedynczy kąt albo z osobna jest nieokreślony, wyznaczona jest tylko ich suma.
W kodzie detekujemy to przez obliczenie . Gdy ta liczba jest mniejsza niż ustalona tolerancja, wybieramy umownie i resztę dopasowujemy tak, by orientacja docelowa była zachowana.
# R36 znane z kroku 6. Indeksy (row, col), 0-based.
sq5_abs = np.hypot(R36[1, 0], R36[1, 1])
if sq5_abs < EPS: # singularność nadgarstka
# tylko (q4 + q6) jednoznaczne — wybieramy q4 = 0
q5 = np.arctan2(0.0, R36[1, 2])
q4 = 0.0
q6 = np.arctan2(-R36[0, 1], R36[0, 0])
solutions.append((q1, q2, q3, q4, q5, q6))
else:
for wrist_sign in (+1, -1): # +1 = noflip, -1 = flip
q5 = np.arctan2( wrist_sign * sq5_abs, R36[1, 2])
q4 = np.arctan2( wrist_sign * R36[2, 2], -wrist_sign * R36[0, 2])
q6 = np.arctan2(-wrist_sign * R36[1, 1], wrist_sign * R36[1, 0])
solutions.append((q1, q2, q3, q4, q5, q6))W trzech miejscach podejmowaliśmy decyzję „znak czy ":
Trzy niezależne wybory binarne dają kombinacji — dlatego typowo dla Pumy mamy do ośmiu różnych konfiguracji osiągających tę samą pozę efektora. W praktyce bywa mniej: jeśli dyskryminant łokcia wyjdzie ujemny, odrzucamy paręshoulder+elbow; jeśli natrafimy na singularność nadgarstka, dwie gałęzie wrist zlewają się w jedną.
Który z tych ośmiu wyników wybrać? Sam algorytm IK tego nie rozstrzyga — selekcja to osobny krok planowania ruchu. Najczęściej stosowane kryteria:
Poniżej miniatury wszystkich rozwiązań dla aktualnej pozy . Kliknięcie ładuje konfigurację do głównego kontrolera powyżej — efektor pozostaje w tym samym miejscu, a ramię robota przybiera drastycznie różne kształty. To daje intuicję, o co walczy wybór gałęzi.
Powyżej każdy krok zawierał już wzór symboliczny. Tutaj wszystkie kroki sklejone w jednym przejściu, dla konkretnej pozy celu. Liczby są stałe — można je przepisać na kalkulator/notebook i sprawdzić własną implementację linijka po linijce, w izolacji od interaktywnego playgroundu wyżej.
[ 0 0 1 ]
[ 0 1 0 ]
[-1 0 0 ]Wszystkie liczby powyżej zostały wyznaczone analitycznie z wzorów wyprowadzonych w krokach 1–7. Można je odtworzyć ręcznie kalkulatorem — ten panel służy jako test kontrolny dla studenta implementującego algorytm we własnym kodzie. Stałe DH: , , , m.
Cała ściąga, którą zobaczysz za chwilę, to litera algorytmu — to, co trzeba wpisać do kodu. Ale za każdym wzorem stoi pomysł algebraiczny: jakaś sprytna kombinacja, która pozwala wyłuskać pojedynczą niewiadomą z układu wielu równań. Jeśli zapamiętasz pomysły, formuły zrekonstruujesz sam; jeśli zapamiętasz tylko formuły, przy pierwszym robocie o innej geometrii będziesz musiał zaczynać od zera.
| Staw | Główny pomysł | Technika algebraiczna |
|---|---|---|
| Rzut na płaszczyznę poziomą eliminuje wszystkie wyższe stawy | Dwa atan2 w różnicy — kąt do celu w XY minus kąt korygujący offset . | |
| Odległość bark↔wrist jest niezmiennikiem obrotów — zależy tylko od . | Prawo cosinusów na trójkącie bark-łokieć-wrist; jedno równanie skalarne . | |
| Po wyznaczeniu układ pozycji staje się liniowy w (). | Pogrupowanie wyrazów pozycji wg i daje układ liniowy 2×2; reguła Cramera zwraca obie wartości i atan2(s_2, c_2) dopina kąt. | |
| Mając , „brakująca" rotacja jest jawnie znana — można z niej odczytać trzy kąty Eulera. | , a następnie porównanie konkretnych elementów macierzy z postacią symboliczną iloczynu trzech rotacji wokół osi nadgarstka. |
Wspólny wzorzec: każdy staw odkrywa swój kąt dzięki jednej szczęśliwej obserwacji, która eliminuje pozostałe niewiadome. Dla tym kluczem jest niezależność rzutu poziomego od ; dla — niezmiennicza długość przekątnej; dla — liniowość po wyłuskaniu ; dla orientacji — rozkład grupowy . Te cztery obserwacje to cała architektura algorytmu; reszta to mechaniczne stosowanie atan2.
Ten sam paradygmat („szukaj kombinacji eliminującej pozostałe niewiadome") sprawdza się też w innych analitycznych IK — dla robotów UR działa równoległość trzech środkowych osi (forma B Piepera), dla robotów Stäubli — inna sekwencja offsetów. Pomysł zostaje, mechanika się zmienia.
Wszystkie 7 snippetów Python z kroków 1–7 powyżej, połączone w jedną funkcję solve_puma560_ik(T_target). Zwraca listę do 8 rozwiązań (każde jako krotka 6 kątów), zgodnie z rozszczepieniem na gałęzie shoulder × elbow × wrist.
import numpy as np
# Parametry DH Puma560 (modified Craig):
A2, A3 = 0.4318, 0.0203
D3, D4 = 0.1500, 0.4331
D6 = 0.0
EPS = 1e-9
def solve_puma560_ik(T_target):
"""Analityczne IK dla Puma560. Zwraca listę krotek (q1..q6)."""
R = T_target[:3, :3]
p = T_target[:3, 3]
# === Krok 1: środek nadgarstka ===
p_wc = p - D6 * R[:, 2]
px_w, py_w, pz_w = p_wc
# === Krok 2: q1 (dwie gałęzie shoulder) ===
r_xy_sq = px_w**2 + py_w**2
disc_q1 = r_xy_sq - D3**2
if disc_q1 < 0:
return []
rho_abs = np.sqrt(max(0.0, disc_q1))
phi = np.arctan2(py_w, px_w)
L = np.sqrt(A3**2 + D4**2)
beta = np.arctan2(D4, A3)
solutions = []
for shoulder_sign in (+1, -1):
rho = shoulder_sign * rho_abs
q1 = phi - np.arctan2(D3, rho)
# === Krok 4: q3 (dwie gałęzie elbow) ===
K = (rho**2 + pz_w**2 - A2**2 - A3**2 - D4**2) / (2 * A2)
disc = L**2 - K**2
if disc < 0:
continue
sqrt_d = np.sqrt(disc)
for elbow_sign in (+1, -1):
q3 = np.arctan2(elbow_sign * sqrt_d, K) - beta
# === Krok 5: q2 (układ 2x2 → atan2) ===
c3, s3 = np.cos(q3), np.sin(q3)
M = A2 + A3 * c3 - D4 * s3
N = A3 * s3 + D4 * c3
denom = M**2 + N**2
c2 = (M * rho - N * pz_w) / denom
s2 = (-M * pz_w - N * rho) / denom
q2 = np.arctan2(s2, c2)
# === Krok 6: R36 = R03.T @ R ===
c1, s1 = np.cos(q1), np.sin(q1)
c23, s23 = np.cos(q2 + q3), np.sin(q2 + q3)
R03 = np.array([
[ c1*c23, -c1*s23, -s1],
[ s1*c23, -s1*s23, c1],
[-s23, -c23, 0],
])
R36 = R03.T @ R
# === Krok 7: q4, q5, q6 (dwie gałęzie wrist + singularność) ===
sq5_abs = np.hypot(R36[1, 0], R36[1, 1])
if sq5_abs < EPS:
q5 = np.arctan2(0.0, R36[1, 2])
q4 = 0.0
q6 = np.arctan2(-R36[0, 1], R36[0, 0])
solutions.append((q1, q2, q3, q4, q5, q6))
else:
for wrist_sign in (+1, -1):
q5 = np.arctan2( wrist_sign * sq5_abs, R36[1, 2])
q4 = np.arctan2( wrist_sign * R36[2, 2], -wrist_sign * R36[0, 2])
q6 = np.arctan2(-wrist_sign * R36[1, 1], wrist_sign * R36[1, 0])
solutions.append((q1, q2, q3, q4, q5, q6))
return solutions
# Weryfikacja: forward kinematics * inverse = identity
if __name__ == "__main__":
from numpy.testing import assert_allclose
# załóż że masz forward_kinematics(q) → T 4×4
q_true = np.array([0.3, -1.2, 1.6, 0.4, 0.5, -0.6])
T = forward_kinematics(q_true) # noqa: F821
sols = solve_puma560_ik(T)
dists = [np.linalg.norm(np.array(s) - q_true) for s in sols]
print(f"#solutions = {len(sols)}, best |q - q*| = {min(dists):.2e}")
# → #solutions = 8, best |q - q*| ≈ 1e-15 (precyzja maszynowa)Ile to linii: ~70 linii czystego kodu, plus ~15 linii parametrów i weryfikacji. Cała aplikacja przemysłowa IK Pumy mieści się w jednym pliku, który czytasz w 5 minut — to jest właśnie wartość zamkniętego rozwiązania względem metod iteracyjnych.
Wszystkie kluczowe wzory algorytmu IK Pumy 560 zebrane w jednym miejscu — przydatne jako kompaktowa referencja przy implementacji albo powtórce przed egzaminem.
Wybór znaków w trzech miejscach (, dyskryminant , ) generuje kombinacji. Część może odpaść jako poza zasięgiem () lub poza limitami przegubów; w singularności nadgarstka () gałęzie wrist zlewają się.
atan2(p_y, p_x) bez − atan2(d₃, ρ) daje błąd ~7° wynikający z odsadzenia .arcsin lub arccos tam gdzie wystarcza atan2 — gubi informację o ćwiartce, błąd zniknie tylko dla niektórych pozycji.W poprzednich krokach kilkukrotnie wspominaliśmy konfiguracje, w których algorytm staje wobec sytuacji nietypowej — osobliwości (w literaturze angielskojęzycznej i wielu polskich tłumaczeniach: singularności). Nie są to błędy algorytmu, tylko miejsca w przestrzeni konfiguracyjnej, w których manipulator fizycznie traci jeden stopień swobody (kierunek, w którym nie może już poruszyć narzędziem, bez zmiany położenia paru przegubów na raz). Algebraicznie objawiają się one tym, że formuła traci jednoznaczność: układ równań ma albo nieskończenie wiele rozwiązań, albo nie ma żadnego.
| Typ osobliwości | Warunek matematyczny | Co się dzieje fizycznie | Konsekwencja dla IK |
|---|---|---|---|
| Osobliwość osi 1 (nad bazą) | (cel na granicy „cylindra zakazanego") | TCP znajduje się dokładnie nad pionową osią bazy — obrót wokół tej osi nie zmienia położenia TCP. | staje się nieokreślony (każdy kąt jest równie dobry). Numerycznie pojawia się atan2(0, 0), które Python zwraca jako 0 — to arbitralny wybór, nie poprawne rozwiązanie. |
| Osobliwość łokcia — wewnętrzna | , czyli — przekątna bark↔wrist osiąga ekstremum. | Ramię jest maksymalnie wyprostowane () lub maksymalnie zgięte (). Trójkąt bark-łokieć-wrist degeneruje się do odcinka. | Dyskryminant — dwie gałęzie elbow up i elbow down zlewają się w jedno rozwiązanie. Algorytm zwróci konkretną liczbę, ale „rezerwa" manewrowa łokcia spadła do zera. |
| Osobliwość łokcia — zewnętrzna (poza zasięgiem) | Cel leży poza fizycznym zasięgiem ramienia — trójkąt o bokach nie spełnia nierówności trójkąta (suma dwóch krótszych boków jest mniejsza od trzeciego). | → wyjątek numeryczny. Trzeba albo odrzucić pozę celu, albo wybrać drugą gałąź barku (czasem cel jest osiągalny tylko z drugiej strony manipulatora). | |
| Osobliwość nadgarstka (gimbal lock) | , czyli . | Osie i stają się współliniowe (oba są tą samą prostą) — obrót wokół jednej jest nieodróżnialny od obrotu wokół drugiej. Klasyczna blokada Eulera/Cardano. | Tylko suma (przy ) lub różnica (przy ) jest wyznaczona; sam i mają nieskończenie wiele par. Trzeba przyjąć umowę (np. ) i dopasować do reszty. |
Wspólny mianownik wszystkich czterech przypadków: dwie kolumny Jakobianu stają się liniowo zależne, więc wyznacznik . Geometrycznie znaczy to, że istnieje kierunek w przestrzeni (kierunek prędkości efektora), w którym żadna kombinacja nie potrafi przesunąć końcówki — manipulator traci jeden stopień swobody w tym kierunku. Dokładniejszą analizą zajmuje się moduł 7; tutaj wystarczy zapamiętać, że nasz solver analityczny nie zawodzi numerycznie w osobliwościach (poza zewnętrzną łokcia, gdzie cel jest fizycznie nieosiągalny) — zwraca konkretne liczby, ale są to liczby arbitralne w obrębie nieskończonej rodziny rozwiązań i nie należy ich traktować jako „tej jedynej poprawnej" konfiguracji.
Cały walkthrough jest realizowany w TypeScript (natywnie w przeglądarce). Poniżej ten sam algorytm, linijka po linijce przeniesiony do Pythona, uruchamiamy przez Pyodide w osobnym web workerze. Oba środowiska startują z identycznej pozy docelowej i powinny zwrócić identyczny zbiór rozwiązań, z różnicą nie większą niż precyzja arytmetyki zmiennoprzecinkowej (~10⁻¹⁵). Pyodide jest ładowany leniwie — pierwsze uruchomienie zajmie kilka sekund (pobranie ~10 MB runtime).
Python runtime (Pyodide) uruchamiany jest w web workerze i wymaga jednorazowego pobrania ~10 MB zasobów. Kliknij, by go zainicjować.
W module 2 bierzemy tę samą matematykę i przenosimy ją w tryb interaktywnego eksperymentu: wszystkie 8 rozwiązań jednocześnie na jednej scenie, sterowanie gałęziami checkboxami, animacja trajektorii pokazująca, jak gałęzie „znikają" przy przekraczaniu granic osiągalności. Moduły 3 i 4 pokazują alternatywne strategie rozwiązania — numeryczne (Jakobian, optymalizacja) — które są wolniejsze i mniej dokładne, ale za to uniwersalne: nie wymagają ręcznego wyprowadzania wzorów dla każdej rodziny robotów i działają dla dowolnej geometrii (także tych, dla których zamknięta forma byłaby trudna do wyprowadzenia).
Mając już z IK, kolejny naturalny krok to: różniczkowanie po czasie i pytanie o moment napędowy w przegubach. Tym zajmuje się moduł 9 (Dynamika odwrotna) — algorytm Newton-Euler na trajektorii (q, q̇, q̈) → siły i momenty napędowe τ. Następnie moduł 10 dodaje model elektromechaniczny silnika DC i przekładni harmonicznej, prowadząc do energii cyklu transportowego.