Moduł 01 · analityczne

Wyprowadzenie analityczne

Puma560 krok-po-kroku: warunek Piepera, środek nadgarstka, q₁…q₆ z geometrii i algebry.

O czym jest ten moduł

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 (q4q_4, q5q_5, q6q_6) 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:

  • UR5/UR10 (Universal Robots) — wrist nie zbiega w punkcie (przesunięcie d5d_5 niezerowe), więc klasyczna „forma A" Piepera nie jest spełniona. Ale q₂, q₃, q₄ są wzajemnie równoległe → spełnia formę B i ma rozwiązanie zamknięte (Hawkins 2013, Kufieta 2014).
  • Manipulatory nie spełniające żadnej z form czasem także mają zamkniętą formę — przez ogólniejsze metody (Raghavan–Roth, redukcja do równania 16. stopnia). Trudniejszą geometrycznie, ale wciąż nie iteracyjną.

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 q1,q2,q3q_1, q_2, q_3, żeby środek nadgarstka trafił w odpowiednie miejsce w przestrzeni, a potem q4,q5,q6q_4, q_5, q_6, żeby narzędzie miało żądaną orientację. Cały rachunek sprowadza się do kilku zastosowań twierdzenia cosinusów i funkcji atan2.

bark (q₂)a₂łokieć (q₃)L = √(a₃²+d₄²)środek nadgarstkaz₄z₅z₆TCPd₆Forma A warunku Piepera (spełniona dla Puma560):Trzy kolejne osie obrotu q₄, q₅, q₆ przecinają się w jednym punkcie — środku nadgarstka.W efekcie: pozycja środka zależy tylko od q₁, q₂, q₃, a orientacja efektora — od q₄, q₅, q₆.
PUMA 560 — rzeczywisty manipulator z nadgarstkiem sferycznym
Forma A warunku Piepera w praktyce: trzy ostatnie osie obrotu nadgarstka manipulatora PUMA 560 przecinają się fizycznie w jednym punkcie wewnątrz obudowy nadgarstka. Dzięki temu obroty q₄, q₅, q₆ zmieniają tylko orientację narzędzia, nie przesuwając środka nadgarstka.
Źródło: Wikimedia Commons · autor: NASA / Dominic Hart · licencja: Public Domain

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:

  • Podrobot pozycjonujący — 3-osiowe ramię (podstawa + bark + łokieć, czyli q1,q2,q3q_1, q_2, q_3), kończące się „kulką" w punkcie przecięcia trzech ostatnich osi. Jego jedynym zadaniem jest umieścić tę kulkę w żądanym punkcie pwc\mathbf{p}_\mathrm{wc} przestrzeni.
  • Podrobot orientujący — 3-osiowy nadgarstek sferyczny (q4,q5,q6q_4, q_5, q_6) zamontowany w tej samej kulce. Jego zadaniem jest obrócić narzędzie do żądanej orientacji nie ruszając kulką, bo wszystkie jego trzy osie obrotu przechodzą dokładnie przez nią.

To rozdzielenie nie jest tylko opisem — jest matematycznym faktem: pochodna pozycji wrist centre względem q4,q5,q6q_4, q_5, q_6 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).

Geometria Pumy560 w konwencji DH (Craig)

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ś z^i\hat{\mathbf{z}}_i musi pokrywać się z osią obrotu przegubu i, a oś x^i1\hat{\mathbf{x}}_{i-1} leży wzdłuż wspólnej normalnej łączącej osie przegubów i−1 i i.

Cztery parametry Denavita–Hartenberga
Cztery parametry DH: skręcenie ogniwa α (zielony), długość ogniwa a (fioletowy), odsadzenie d (czerwony), kąt przegubu θ (niebieski). Osie pokazane na dwóch sąsiednich ogniwach.
Źródło: Wikimedia Commons · autor: Jahobr · licencja: Public Domain

Cztery parametry na ogniwo — αi1\alpha_{i-1}, ai1a_{i-1}, did_i, θi\theta_i — w pełni wyznaczają transformację z układu {i−1} do {i}. Intuicja geometryczna każdego z nich:

aᵢ₋₁długość ogniwa i−1

Odległość między osiami przegubów i−1 oraz i, mierzona wzdłuż ich wspólnej normalnej.

αᵢ₋₁skręcenie ogniwa i−1

Kąt obrotu osi i−1 do i, obracany wokół wspólnej normalnej.

dᵢodsadzenie przegubu i

Przesunięcie początku układu {i} wzdłuż osi i od wspólnej normalnej.

θᵢkąt przegubu i

Kąt obrotu osi i−1 do i wokół i. Dla przegubu obrotowego jest to zmienna konfiguracji.

Cztery elementarne transformacje składające się na T_{i−1}^{i}

Transformację z układu {i−1} do {i} można rozbić na cztery kolejne operacje wzdłuż lub wokół pojedynczej osi:  Ti1i=Rotx(αi1)Transx(ai1)Rotz(θi)Transz(di)\;T_{i-1}^{i} = \mathrm{Rot}_x(\alpha_{i-1})\,\mathrm{Trans}_x(a_{i-1})\,\mathrm{Rot}_z(\theta_i)\,\mathrm{Trans}_z(d_i). Każdy krok wprowadza jeden pośredni układ pomocniczy:

Krok 1 — rotacja wokół osi x o kąt α_{i-1}
Krok 1 — Rot_x(α): obrót wokół osi x̂_{i−1} o kąt skręcenia ogniwa.
Źródło: Wikimedia Commons · autor: Jahobr · licencja: Public Domain
Krok 2 — translacja wzdłuż osi x o długość a_{i-1}
Krok 2 — Trans_x(a): translacja wzdłuż osi x̂ o długość ogniwa.
Źródło: Wikimedia Commons · autor: Jahobr · licencja: Public Domain
Krok 3 — rotacja wokół osi z o kąt θ_i
Krok 3 — Rot_z(θ): obrót wokół nowej osi ẑ o kąt przegubu (zmienna konfiguracji).
Źródło: Wikimedia Commons · autor: Jahobr · licencja: Public Domain
Krok 4 — translacja wzdłuż osi z o odsadzenie d_i
Krok 4 — Trans_z(d): translacja wzdłuż osi ẑ o odsadzenie przegubu.
Źródło: Wikimedia Commons · autor: Jahobr · licencja: Public Domain

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.

Klasyczna vs zmodyfikowana konwencja DH

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.

Klasyczna konwencja DH (Denavit-Hartenberg 1955)
Klasyczna DH: układ {i} przypięty do końca ogniwa i. Transformacja T_{i-1}^i = Rot_z(θᵢ)·Trans_z(dᵢ)·Trans_x(aᵢ)·Rot_x(αᵢ). Wszystkie cztery parametry mają ten sam indeks i.
Źródło: Wikimedia Commons · autor: Pushpendra050 · licencja: CC BY-SA 4.0
Zmodyfikowana konwencja DH (Craig 1986)
Zmodyfikowana DH (Craig): układ {i} przypięty do początku ogniwa i, na osi przegubu i. Transformacja T_{i-1}^i = Rot_x(αᵢ₋₁)·Trans_x(aᵢ₋₁)·Rot_z(θᵢ)·Trans_z(dᵢ). Indeksy α i a to i−1, nie i.
Źródło: Wikimedia Commons · autor: Ollydbg · licencja: GFDL
AspektKlasyczna 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ś z^i\hat{\mathbf{z}}_i = oś przegubu i)
z^i\hat{\mathbf{z}}_iPokrywa się z osią przegubu i+1Pokrywa się z osią przegubu i
Transformacja Ti1iT_{i-1}^iRz(θi)Tz(di)Tx(ai)Rx(αi)R_z(\theta_i)\,T_z(d_i)\,T_x(a_i)\,R_x(\alpha_i)Rx(αi1)Tx(ai1)Rz(θi)Tz(di)R_x(\alpha_{i-1})\,T_x(a_{i-1})\,R_z(\theta_i)\,T_z(d_i)
Indeksy parametrówWszystkie mają indeks i (αi,ai,di,θi\alpha_i, a_i, d_i, \theta_i)Mieszane: αi1,ai1\alpha_{i-1}, a_{i-1} i di,θid_i, \theta_i
Kolejność operacjiNajpierw rotacja wokół zz (θ\theta)Najpierw rotacja wokół xx (α\alpha)
Intuicja indeksuParametry „opisują" ogniwo i−1→i jako całośćα,a\alpha, a charakteryzują ogniwo stałe, θ,d\theta, d — przegub zmienny
Obsługa manipulatorów rozgałęzionychTrudna (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 α,a,d,θ\alpha, a, d, \theta 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 ±π/2\pm\pi/2, co od razu upraszcza iloczyny macierzy). Wybór ma jednak charakter notacyjny — to nie jest „jedyna poprawna" konwencja.

z₁z₂ ⊙q₂a₂q₃d₃(w głąb)a₃d₄wrist centerTCPPuma 560 (Craig, wym. [m])a₂ = 0.4318ramięa₃ = 0.0203offset łokciad₃ = 0.1254odsunięcied₄ = 0.4318przedramięd₆ = 0szare odcinki = długości aᵢbrązowe odcinki = odsunięcia dᵢ⊙ = oś prostopadła do ekranua₃ na rysunku celowo powiększone (~×2.5)
iiαi1\alpha_{i-1}ai1  [m]a_{i-1}\;[\text{m}]di  [m]d_i\;[\text{m}]θi\theta_i
1000000q1q_1
2π/2-\pi/20000q2q_2
300a2=0.4318a_2 = 0.4318d3=0.1254d_3 = 0.1254q3q_3
4π/2-\pi/2a3=0.0203a_3 = 0.0203d4=0.4318d_4 = 0.4318q4q_4
5π/2\pi/20000q5q_5
6π/2-\pi/20000q6q_6

a2a_2 to długość ramienia (upper arm), d4d_4 to długość przedramienia (forearm). Dodatkowo w Pumie występują dwa niezerowe odsadzenia:d3d_3 (odsadzenie przedramienia od osi ramienia, ok. 12,5 cm) oraz a3a_3 (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:

i1Ti=Rotx(αi1)Transx(ai1)Rotz(θi)Transz(di){}^{i-1}T_{i} = \mathrm{Rot}_x(\alpha_{i-1})\,\mathrm{Trans}_x(a_{i-1})\,\mathrm{Rot}_z(\theta_i)\,\mathrm{Trans}_z(d_i)
=[cθsθ0ai1sθcαcθcαsαsαdisθsαcθsαcαcαdi0001]= \begin{bmatrix} c\theta & -s\theta & 0 & a_{i-1} \\ s\theta\,c\alpha & c\theta\,c\alpha & -s\alpha & -s\alpha\,d_i \\ s\theta\,s\alpha & c\theta\,s\alpha & c\alpha & c\alpha\,d_i \\ 0 & 0 & 0 & 1 \end{bmatrix}

gdzie cθ=cosθic\theta = \cos\theta_i, sθ=sinθis\theta = \sin\theta_i itd. Mnożąc kolejne takie macierze, dostaniemy T06T_0^{6} — pełne przekształcenie z bazy do układu efektora.

Co fizycznie robi każdy z sześciu przegubów

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 q3q_3 kontroluje wysokość, a q5q_5 nie".

q₁Swivel podstawy
Oś obrotu: pionowa oś bazy (ẑ₀)
Dolny obrót żurawia portowego — cała wieża obraca się wokół własnej pionowej osi.

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 (px,py)(p_x, p_y) docelowej pozycji i pytamy: w którą stronę ramię musi się skręcić, żeby w ogóle „móc patrzeć" na cel.

q₂Bark — pochylenie ramienia
Oś obrotu: pozioma, prostopadła do ramienia
Podnosisz wyprostowaną rękę w górę lub w dół — bez zginania łokcia. To dokładnie q₂.

Wraz z q3q_3 kontroluje punkt 2D, w który sięga końcówka ramienia (wysokość + wysunięcie radialne) — w płaszczyźnie pionowej wyznaczonej już przez q1q_1. Czysto „pionowy" stopień swobody barku.

q₃Łokieć — zgięcie przedramienia
Oś obrotu: równoległa do osi barku
Zginasz własną rękę w łokciu — dłoń przybliża się do barku, choć ramię się nie obraca.

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".

q₄Roll przedramienia
Oś obrotu: wzdłuż osi przedramienia
Obracasz nadgarstek trzymając śrubokręt — sama dłoń kręci się wokół osi narzędzia.

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.

q₅Pitch nadgarstka
Oś obrotu: prostopadła do z₄
Zginasz dłoń w nadgarstku — od pozycji wyprostowanej (ramię w jednej linii z dłonią) do pełnego zgięcia.

Kąt ataku narzędzia względem osi przedramienia. Gdy q5=0q_5 = 0, narzędzie jest „wyciągnięte" w jednej linii z przedramieniem — wówczas występuje gimbal lock: osie q4q_4 i q6q_6 pokrywają się i tracą niezależność.

q₆Roll TCP
Oś obrotu: wzdłuż osi narzędzia (ẑ₆)
Klucz nasadowy włożony w nakrętkę — kierunek wskazany przez klucz ustalają q₄ i q₅, ale obrót klucza wokół własnej osi to dokładnie q₆.

Końcowy stopień swobody — orientuje narzędzie wokół jego własnej osi po tym, jak q4,q5q_4, q_5 ustawiły kierunek, w który narzędzie „patrzy". Bez q6q_6 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 (q1,q2,q3q_1, q_2, q_3) ustalają pozycję 3-wymiarowego punktu (środek nadgarstka) — bo 3 niewiadome wystarczają, żeby opisać dowolny punkt w R3\mathbb{R}^3. Trzy ostatnie (q4,q5,q6q_4, q_5, q_6) ustalają orientację w przestrzeni — bo grupa obrotów SO(3)SO(3) 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 SE(3)SE(3) przy możliwie najprostszej geometrii.

Laboratorium

Po lewej: manipulator, którym sterujesz ręcznie sześcioma suwakami — służy do generowania dowolnych poz testowych. Po prawej: pole pozy docelowej TT^* z przyciskiem „zrzut z kontrolera", który kopiuje aktualne T06T_0^{6} 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.

Konfiguracja przegubów

θ₁2.8°
θ₂-77.7°
θ₃6.8°
θ₄8.5°
θ₅-19.3°
θ₆172.0°

Poza efektora (T₀⁶)

Pozycja [m]
x = 0.5000
y = 0.1500
z = 0.3000
Orientacja RPY [°]
R = 0.00
P = 90.00
Y = 0.00
φ = atan2(p_y, p_x)= 16.699°ρ = √(p_x² + p_y² − d₃²)= 0.5067 mq₁ = φ − atan2(d₃, ρ)= 2.800°L = √(a₃² + d₄²)= 0.4323 mβ = atan2(d₄, a₃)= 87.308°K= -0.0307L² − K² (elbow discriminant)= 0.1859 q₃ (elbow up)= 6.768°M = a₂ + a₃c₃ − d₄s₃= 0.4011N = a₃s₃ + d₄c₃= 0.4312q₂= -77.699°q₄, q₅, q₆ (z R₃⁶)= -171.49°, 19.27°, -8.04°

Poza docelowa T*

Pozycja [m]
Orientacja RPY [°]
Krok 0

Odseparowanie pozycji od orientacji (dekompozycja wristowa)

Zadana poza efektora to macierz transformacji jednorodnej 4×44\times 4 — górne 3×33\times 3 to orientacja (macierz rotacji RR), a ostatnia kolumna to pozycja p\mathbf{p}:

T=[Rp01]SE(3)T^* = \begin{bmatrix} R & \mathbf{p} \\ \mathbf{0}^\top & 1 \end{bmatrix} \in SE(3)

Zbiór SE(3)SE(3) to zbiór wszystkich takich macierzy — grupa sztywnych przesunięć i obrotów w przestrzeni. Szukamy wektora kątów q=(q1,,q6)q = (q_1, \dots, q_6) spełniającego f(q)=Tf(q) = T^*, gdzie ff 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 pwc\mathbf{p}_\mathrm{wc}). Obrót wokół osi przechodzącej przez punkt nie przesuwa tego punktu. Wobec tego położenie pwc\mathbf{p}_\mathrm{wc} nie zależy odq4,q5,q6q_4, q_5, q_6 — zależy tylko od q1,q2,q3q_1, q_2, q_3.

W konwencji Craiga dla Pumy zachodzi d6=0d_6 = 0, 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:

pwc=pd6Rz^  =d6=0  p\mathbf{p}_\mathrm{wc} = \mathbf{p} - d_6 \cdot R\,\hat{\mathbf{z}} \;\overset{d_6 = 0}{=}\; \mathbf{p}

Gdy jest narzędzie (pewien stały offset TtoolT_\mathrm{tool} między układem {6} a układem TCP), cofamy się o ten offset:

T06  =  TTtool1T_0^{6} \;=\; T^* \cdot T_\mathrm{tool}^{-1}

Reasumując plan: najpierw znajdziemy q1,q2,q3q_1, q_2, q_3 z pozycji środka nadgarstka (3 równania, 3 niewiadome), potem q4,q5,q6q_4, q_5, q_6 z warunku, że orientacja pozostałych trzech przegubów ma domknąć zadane RR.

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 kolumna
Krok 1

Wzór na pozycję środka nadgarstka

Pozycja początku układu {4} (czyli środka nadgarstka) w bazie to ostatnia kolumna iloczynu czterech macierzy T04=T01T12T23T34T_0^4 = T_0^1\,T_1^2\,T_2^3\,T_3^4. Wstawiając parametry z tabeli DH (patrz wyprowadzenie poniżej), otrzymujemy:

px=c1(a2c2+a3c23d4s23)d3s1py=s1(a2c2+a3c23d4s23)+d3c1pz=a2s2a3s23d4c23\begin{aligned} p_x &= c_1\,\bigl(a_2\,c_2 + a_3\,c_{23} - d_4\,s_{23}\bigr) - d_3\,s_1 \\ p_y &= s_1\,\bigl(a_2\,c_2 + a_3\,c_{23} - d_4\,s_{23}\bigr) + d_3\,c_1 \\ p_z &= -a_2\,s_2 - a_3\,s_{23} - d_4\,c_{23} \end{aligned}

gdzie ci=cosqic_i = \cos q_i, si=sinqis_i = \sin q_i, c23=cos(q2+q3)c_{23} = \cos(q_2+q_3), s23=sin(q2+q3)s_{23} = \sin(q_2+q_3).

▸ Pokaż pełne wyprowadzenie iloczynu T04

Każda Ti1iT_{i-1}^i w konwencji Craiga ma postać podaną wyżej. Dla Pumy 560 podstawiamy z tabeli DH:

T01=Rz(q1)=[c1s100s1c10000100001]T_0^1 = R_z(q_1) = \begin{bmatrix} c_1 & -s_1 & 0 & 0 \\ s_1 & c_1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}

(α0=0,a0=0,d1=0\alpha_0 = 0,\,a_0 = 0,\,d_1 = 0, więc tylko obrót wokół z).

T12=Rx(π/2)Rz(q2)=[c2s2000010s2c2000001]T_1^2 = R_x(-\pi/2)\,R_z(q_2) = \begin{bmatrix} c_2 & -s_2 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ -s_2 & -c_2 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}

(α1=π/2,a1=0,d2=0\alpha_1 = -\pi/2,\,a_1 = 0,\,d_2 = 0).

T23=Transx(a2)Rz(q3)Transz(d3)=[c3s30a2s3c300001d30001]T_2^3 = \mathrm{Trans}_x(a_2)\,R_z(q_3)\,\mathrm{Trans}_z(d_3) = \begin{bmatrix} c_3 & -s_3 & 0 & a_2 \\ s_3 & c_3 & 0 & 0 \\ 0 & 0 & 1 & d_3 \\ 0 & 0 & 0 & 1 \end{bmatrix}

(α2=0,a2=a2,d3=d3\alpha_2 = 0,\,a_2 = a_2,\,d_3 = d_3).

T34=Rx(π/2)Transx(a3)Rz(q4)Transz(d4)T_3^4 = R_x(-\pi/2)\,\mathrm{Trans}_x(a_3)\,R_z(q_4)\,\mathrm{Trans}_z(d_4)

Przy szukaniu pozycji środka nadgarstka (pierwsze 3 elementy ostatniej kolumny T04T_0^4) nie potrzebujemy q4q_4: wystarczy „adres" początku układu {4}, czyli punkt (a3,0,d4)(a_3, 0, d_4)^\top wyrażony w bazie. Liczymy więc krócej:

pwc=T01T12T23T34(0,0,0,1)=T03(a3,0,d4,1)\mathbf{p}_\mathrm{wc} = T_0^1\,T_1^2\,T_2^3\,T_3^4\,(0,0,0,1)^\top = T_0^3\,(a_3,\,0,\,d_4,\,1)^\top

gdzie ostatnia transformacja T34T_3^4 sprowadza się do obrotu wokół x i przesunięcia o a3a_3 w x i d4d_4 w z.

Wymnażając kolejno: najpierw T23T34T_2^3\,T_3^4 daje przesunięcie łokcia o (a3,0,d4)(a_3, 0, d_4) obrócone przez Rz(q3)R_z(q_3), potem T12T_1^2 obraca to wokół osi x o π/2-\pi/2 i dodaje przegub q2q_2, finalnie T01T_0^1 obraca całość wokół pionowej osi bazy. W rezultacie dostajemy podane wyżej px,py,pzp_x, p_y, p_z — z efektywną długością przedramienia rozłożoną na składowe a3c23d4s23a_3 c_{23} - d_4 s_{23} w „radialnym" kierunku oraz a3s23d4c23-a_3 s_{23} - d_4 c_{23} w pionie.

Zauważmy, że w pierwszych dwóch równaniach powtarza się to samo wyrażenie w nawiasie. Wprowadźmy skrót:

ρ    a2c2+a3c23d4s23\rho \;\equiv\; a_2\,c_2 + a_3\,c_{23} - d_4\,s_{23}

Geometrycznie ρ\rho to radialna odległość środka nadgarstka od pionowej osi bazy w obróconym o q1q_1 układzie — innymi słowy, „jak daleko" wrist centre jest od pionowej osi obrotu barku.

Z dwóch pierwszych równań mamy px=c1ρd3s1p_x = c_1\,\rho - d_3\,s_1 i py=s1ρ+d3c1p_y = s_1\,\rho + d_3\,c_1. Są to składowe obrotu o kąt q1q_1 wektora (ρ,d3)(\rho, d_3) w płaszczyźnie XY:

[pxpy]=[c1s1s1c1]Rz(q1)[ρd3]\begin{bmatrix} p_x \\ p_y \end{bmatrix} = \underbrace{\begin{bmatrix} c_1 & -s_1 \\ s_1 & c_1 \end{bmatrix}}_{R_z(q_1)} \begin{bmatrix} \rho \\ d_3 \end{bmatrix}

Mnożąc obie strony tej równości przez Rz(q1)R_z(q_1)^{\top} (transpozycja = inwersja, bo macierz obrotu jest ortogonalna) i odczytując:

[ρd3]=[c1s1s1c1][pxpy]    {ρ=pxcosq1+pysinq1d3=pxsinq1+pycosq1\begin{bmatrix} \rho \\ d_3 \end{bmatrix} = \begin{bmatrix} c_1 & s_1 \\ -s_1 & c_1 \end{bmatrix} \begin{bmatrix} p_x \\ p_y \end{bmatrix} \;\Longleftrightarrow\; \begin{cases} \rho = p_x\cos q_1 + p_y\sin q_1 \\ d_3 = -p_x\sin q_1 + p_y\cos q_1 \end{cases}

Podnosząc oba równania do kwadratu i dodając, lewa strona daje ρ2+d32\rho^2 + d_3^2, a prawa (pxc1+pys1)2+(pxs1+pyc1)2(p_x c_1 + p_y s_1)^2 + (-p_x s_1 + p_y c_1)^2 = po rozwinięciu i uproszczeniu (krzyżowe wyrazy się znoszą, c12+s12=1c_1^2 + s_1^2 = 1) px2+py2p_x^2 + p_y^2. Stąd kluczowa identyczność:

  px2+py2  =  ρ2+d32  \boxed{\;p_x^2 + p_y^2 \;=\; \rho^2 + d_3^2\;}

Ta równość jest fundamentem całej dekompozycji. Lewa strona zależy wyłącznie od zadanych px,pyp_x, p_y (wiemy je), prawa strona zawiera tylko ρ\rho. Dlatego ρ\rho wyliczamy natychmiast — bez znajomości q2,q3q_2, q_3.

Python · krok 1implementacja: pozycja wrist centre w bazie
# 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_wc
Krok 2

Wyznaczenie q₁ — dwie gałęzie barku

Intuicja: 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 (px,py,0)(p_x, p_y, 0). Naturalna pierwsza myśl: obrót podstawy musi po prostu skierować ramię w stronę cienia — czyli q1=atan2(py,px)q_1 = \operatorname{atan2}(p_y, p_x).

Dla manipulatora bez odsadzenia tak właśnie by było. Puma ma jednak d312,5cmd_3 \approx 12{,}5\,\mathrm{cm} boczne odsunięcie osi przedramienia od osi ramienia. Skutkiem tego, gdy bark obraca się o kąt q1q_1, 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ść d3d_3.

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: q1=atan2(py,px)atan2(d3,ρ)q_1 = \operatorname{atan2}(p_y, p_x) - \operatorname{atan2}(d_3, \rho). 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 (ρ,d3\rho, d_3) w obróconym układzie ramienia to ten sam wektor co (px,pyp_x, p_y) w bazie — tylko obrócony o q1q_1. Wektor (ρ,d3\rho, d_3) ma nachylenie atan2(d3,ρ)\operatorname{atan2}(d_3, \rho) względem osi ramienia, a wektor (px,pyp_x, p_y) ma nachylenie atan2(py,px)\operatorname{atan2}(p_y, p_x) względem osi X bazy. Różnica tych dwóch nachyleń to właśnie kąt obrotu między układami, czyli q1q_1.

Z identyczności kroku 1 wyznaczamy ρ\rho, ale z dokładnością do znaku (bo pierwiastek kwadratowy daje dwie wartości):

ρ=±px2+py2d32\rho = \pm\sqrt{p_x^2 + p_y^2 - d_3^2}

Jeśli px2+py2<d32p_x^2 + p_y^2 < d_3^2, zadanie jest nieosiągalne — cel leży w „zakazanym" cylindrze o promieniu d3d_3 wokół pionowej osi bazy. W pozostałych przypadkach znak ρ\rho wyznacza, czy bark Pumy zwraca się w kierunku celu (shoulder right, ρ>0\rho > 0), czy odwrotnie, o π\pi obrócony (shoulder left, ρ<0\rho < 0). Obie konfiguracje dają identyczną pozycję wrist centre — stąd „gałąź" w algorytmie.

Mając ρ\rho i d3d_3, rozwiązujemy układ dwóch liniowych równań w (cosq1,sinq1)(\cos q_1, \sin q_1) z poprzedniej strony:

{pxcosq1+pysinq1=ρpxsinq1+pycosq1=d3\begin{cases} p_x\cos q_1 + p_y\sin q_1 = \rho \\ -p_x\sin q_1 + p_y\cos q_1 = d_3 \end{cases}

Formalnie rozwiązujemy go przez wyznacznik, ale ładniej zapisuje się to przez różnicę dwóch funkcji atan2q1q_1 jest różnicą kąta do celu (px,pyp_x, p_y) i kąta „w bok" o odsadzeniu d3d_3:

  q1  =  atan2(py,px)    atan2(d3,ρ)  \boxed{\;q_1 \;=\; \operatorname{atan2}(p_y,\,p_x) \;-\; \operatorname{atan2}(d_3,\,\rho)\;}
Python · krok 2dwie gałęzie shoulder dla q₁
# ρ — 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łęzi

Dwuargumentowa funkcja atan2(y, x) — w odróżnieniu od zwykłego arctan(y/x) — zwraca kąt na pełnym przedziale (π,π](-\pi, \pi], 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 d3d_3.

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łęziach
Krok 3

Przejście do 2-wymiarowej płaszczyzny ramienia

Po wyznaczeniu q1q_1 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 (ρ,z)(\rho, z), gdzie ρ\rho to odległość od pionowej osi bazy (dodatnia po „stronie" celu), a zz — wysokość nad płaszczyzną podstawy.

W tej płaszczyźnie zostało nam zadanie 2-wymiarowe: połączyć bark (0,0)(0, 0) ze środkiem nadgarstka (ρ,z)(\rho, z) za pomocą dwóch ogniw — ramienia o długości a2a_2 i przedramienia o pewnej efektywnej długości LL, z łokciem gdzieś pomiędzy.

Czym jest LL? 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: a3a_3 wzdłuż osi X układu {3} (mały bok, ok. 2 cm) i d4d_4 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ść:

L=a32+d420,4323mL = \sqrt{a_3^2 + d_4^2} \approx 0{,}4323\,\mathrm{m}

i jest odchylony od osi X układu {3} o stały kąt:

β=atan2(d4,a3)87,3\beta = \operatorname{atan2}(d_4, a_3) \approx 87{,}3^\circ

O tym odchyleniu pamiętamy przy obliczaniu q3q_3 w następnym kroku.

ρ [m]z [m]-0.3-0.10.10.30.50.7-0.10.10.3bark (pocz. układu {2})środek nadgarstkaρ = 0.507 m, z = 0.300 mD = 0.589 melbow ↑ (q₂ = -77.7°, q₃ = 6.8°)a₂Lelbow ↓ (q₂ = 16.4°, q₃ = -181.4°)Stałe DH (Craig, Puma560):a₂ = 0.4318 m (ramię)L = √(a₃²+d₄²) = 0.4323 mβ = atan2(d₄, a₃) = 87.31°

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.

Krok 4

Wyznaczenie q₃ — twierdzenie cosinusów, dwie gałęzie łokcia

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ć q3q_3 w izolacji od reszty: odległość |bark → wrist center| nie zależy od q1q_1 ani q2q_2 — zależy wyłącznie od q3q_3.

Intuicyjnie: q1q_1 i q2q_2 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 q3q_3.

Stąd ważny wniosek metodologiczny: mając zadaną pozycję wrist centre (a więc znaną przekątną D=ρ2+pz2D = \sqrt{\rho^2 + p_z^2}), możemy natychmiast wyłuskać q3q_3 z jednego równania skalarnego, nie wiedząc jeszcze ile wynosi q2q_2. Dopiero potem wrócimy do dwuwymiarowego problemu i wyznaczymy q2q_2, traktując kształt trójkąta jako już ustalony. To dlatego algorytm liczy stawy w kolejności q1q3q2q_1 \to q_3 \to q_2, a nie w kolejności indeksów.

Spojrzenie geometryczne. W płaszczyźnie ramienia bark AA, łokieć EE i środek nadgarstka WW tworzą trójkąt z trzema znanymi długościami: a2a_2 (ramię), LL (efektywne przedramię) i D=ρ2+pz2D = \sqrt{\rho^2 + p_z^2} (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 γ\gamma przy łokciu:

A — barkE — łokiećW — środek nadgarstkaa₂LDγPrawo cosinusów dla △ AEW:D² = a₂² + L² − 2·a₂·L·cos γStąd:cos γ = (a₂² + L² − D²) / (2·a₂·L)Kąt q₃ względem γ:q₃ = π − γ − β

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 DD, której końce są przymocowane w barku AA i w środku nadgarstka WW. Łokieć EE jest swobodnym przegubem o ustalonej odległości a2a_2 od barku i LL 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 DD:

  • elbow up — łokieć powyżej linii bark↔wrist, ramię „zgięte w górę" jak u człowieka podnoszącego rękę.
  • elbow down — łokieć poniżej, ramię „złożone w dół" jak u robota chwytającego coś z poziomu stołu.

Obie konfiguracje dają dokładnie tę samą pozycję wrist centre, ale wymagają różnych q2,q3q_2, q_3. Stąd binarny wybór znaku ±\pm\sqrt{\,\cdot\,} w wyrażeniu na sinγ\sin\gamma — to przełącznik między dwoma lustrzanymi obrazami tego samego trójkąta.

Z prawa cosinusów (zastosowanego do boku DD przeciwległego kątowi γ\gamma):

D2=a22+L22a2Lcosγ    cosγ=a22+L2D22a2LD^2 = a_2^2 + L^2 - 2\,a_2\,L\,\cos\gamma \;\Rightarrow\; \cos\gamma = \frac{a_2^2 + L^2 - D^2}{2\,a_2\,L}

Pełny kąt wewnętrzny trójkąta wyciąga się standardowo przez atan2 (a nie acos), żeby zachować informację o znaku sinγ\sin\gamma:

γ=atan2(±1cos2γ,  cosγ)\gamma = \operatorname{atan2}\bigl(\pm\sqrt{1 - \cos^2\gamma},\;\cos\gamma\bigr)

Znak ±\pm przed sinγ\sin\gamma odpowiada wyborowi elbow up (łokieć powyżej linii bark↔wrist) lub elbow down (poniżej) — dwa odbicia trójkąta względem boku DD.

Powiązanie z q₃. Kąt γ\gamma przy łokciu nie jest tożsamy z q3q_3. Gdyby przedramię biegło dokładnie wzdłuż osi x^3\hat{x}_3 (czyli d4=0d_4 = 0), wówczas q3=πγq_3 = \pi - \gamma. W Pumie jednak przedramię jest „odchylone" od osi x^3\hat{x}_3 o stały kąt β=atan2(d4,a3)87,3\beta = \operatorname{atan2}(d_4, a_3) \approx 87{,}3^\circ (z kroku 3), więc:

q3  =  πγβq_3 \;=\; \pi - \gamma - \beta
▸ Pokaż równoważne wyprowadzenie algebraiczne (z dwóch równań pozycji)

Biorąc dwie zależności z kroku 1, ρ=a2c2+a3c23d4s23\rho = a_2 c_2 + a_3 c_{23} - d_4 s_{23} i pz=a2s2a3s23d4c23p_z = -a_2 s_2 - a_3 s_{23} - d_4 c_{23}, podnosimy oba do kwadratu i dodajemy. Wyrazy z c2s23c_2 s_{23} i s2c23s_2 c_{23} składają się w sumę kątów sin(q2(q2+q3))=sinq3\sin(q_2 - (q_2 + q_3)) = -\sin q_3, podobnie c2c23+s2s23=cosq3c_2 c_{23} + s_2 s_{23} = \cos q_3. Po uproszczeniu:

ρ2+pz2  =  a22+a32+d42+2a2(a3c3d4s3)\rho^2 + p_z^2 \;=\; a_2^2 + a_3^2 + d_4^2 + 2a_2\bigl(a_3\,c_3 - d_4\,s_3\bigr)

Oznaczając lewą stronę (znaną) przez KK:

K  =  ρ2+pz2a22a32d422a2  =  a3c3d4s3K \;=\; \frac{\rho^2 + p_z^2 - a_2^2 - a_3^2 - d_4^2}{2 a_2} \;=\; a_3\,c_3 - d_4\,s_3

Wyrażenie a3c3d4s3a_3\,c_3 - d_4\,s_3 to liniowa kombinacja sinusa i cosinusa o tym samym argumencie — można ją zapisać jako jeden cosinus z przesuniętym kątem:

a3cosq3d4sinq3  =  Lcos(q3+β)a_3\cos q_3 - d_4\sin q_3 \;=\; L\,\cos(q_3 + \beta)

(sprawdź: rozwijając prawą stronę przez cos(q3+β)=c3cosβs3sinβ\cos(q_3+\beta) = c_3\cos\beta - s_3\sin\beta i podstawiając Lcosβ=a3L\cos\beta = a_3, Lsinβ=d4L\sin\beta = d_4 — obie strony się zgadzają).

Stąd cos(q3+β)=K/L\cos(q_3 + \beta) = K/L, więc sin(q3+β)=±1K2/L2\sin(q_3 + \beta) = \pm\sqrt{1 - K^2/L^2}, i:

q3  =  atan2 ⁣(±L2K2,  K)    βq_3 \;=\; \operatorname{atan2}\!\bigl(\pm\sqrt{L^2 - K^2},\;K\bigr) \;-\; \beta

Spójność z wyprowadzeniem geometrycznym: K=LcosγK = -L\cos\gamma (różnica znaku wynika z definicji kąta), więc cos(q3+β)=cosγ=cos(πγ)\cos(q_3+\beta) = -\cos\gamma = \cos(\pi - \gamma), czyli q3+β=πγq_3 + \beta = \pi - \gamma — ten sam wynik.

Postać końcowa (wykorzystywana w kodzie — jeden atan2 zamiast łańcucha kątów):

  q3  =  atan2 ⁣(±L2K2,  K)    β  \boxed{\;q_3 \;=\; \operatorname{atan2}\!\bigl(\pm\sqrt{L^2 - K^2},\;K\bigr) \;-\; \beta\;}

Znak ±\pm odpowiada dwóm gałęziom geometrycznym — łokieć „w górę" lub „w dół". Gdy L2<K2L^2 < K^2, 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).

Python · krok 4dwie gałęzie elbow dla q₃
# 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łęzi
Krok 5

Wyznaczenie q₂ — 2×2 układ liniowy w cosinusach

Mając już znane q3q_3, wracamy do dwóch równań pozycji z kroku 1 — tym razem traktując q2q_2 jako jedyną niewiadomą. Składowe radialna i pionowa to:

ρ  =  a2c2+a3c23d4s23,pz  =  a2s2a3s23d4c23\rho \;=\; a_2\,c_2 + a_3\,c_{23} - d_4\,s_{23}, \qquad p_z \;=\; -a_2\,s_2 - a_3\,s_{23} - d_4\,c_{23}

Stosujemy wzory na sumę kątów: c23=c2c3s2s3c_{23} = c_2 c_3 - s_2 s_3 oraz s23=s2c3+c2s3s_{23} = s_2 c_3 + c_2 s_3. Wstawiając do pierwszego równania:

ρ  =  a2c2+a3(c2c3s2s3)d4(s2c3+c2s3)\rho \;=\; a_2 c_2 + a_3(c_2 c_3 - s_2 s_3) - d_4(s_2 c_3 + c_2 s_3)

Grupujemy wyrazy przy c2c_2 i s2s_2 osobno:

ρ  =  (a2+a3c3d4s3)c2    (a3s3+d4c3)s2\rho \;=\; (a_2 + a_3 c_3 - d_4 s_3)\,c_2 \;-\; (a_3 s_3 + d_4 c_3)\,s_2

Tak samo dla pzp_z: (a3s3+d4c3)c2(a2+a3c3d4s3)s2-(a_3 s_3 + d_4 c_3)\,c_2 - (a_2 + a_3 c_3 - d_4 s_3)\,s_2. Wprowadzając skróty na powtarzające się grupy:

M    a2+a3c3d4s3,N    a3s3+d4c3M \;\equiv\; a_2 + a_3\,c_3 - d_4\,s_3, \qquad N \;\equiv\; a_3\,s_3 + d_4\,c_3

(zauważ: M,NM, N zależą tylko od znanego q3q_3 i stałych DH), oba równania zapisują się zwarto jako układ liniowy 2×2:

[MNNM][c2s2]  =  [ρpz]\begin{bmatrix} M & -N \\ -N & -M \end{bmatrix} \begin{bmatrix} c_2 \\ s_2 \end{bmatrix} \;=\; \begin{bmatrix} \rho \\ p_z \end{bmatrix}

Wyznacznik macierzy współczynników: det=M(M)(N)(N)=(M2+N2)\det = M\cdot(-M) - (-N)\cdot(-N) = -(M^2 + N^2). Z reguły Cramera (lub bezpośrednio przez inwersję 2×2):

c2  =  MρNpzM2+N2,s2  =  MpzNρM2+N2c_2 \;=\; \frac{M\,\rho - N\,p_z}{M^2 + N^2}, \qquad s_2 \;=\; \frac{-M\,p_z - N\,\rho}{M^2 + N^2}

Wyznacznik M2+N2M^2 + N^2 znika tylko, gdy jednocześnie a2=a3=d4=0a_2 = a_3 = d_4 = 0 — niemożliwe dla fizycznego manipulatora. Mając c2c_2 i s2s_2 jednocześnie, kąt wyznacza atan2:

  q2  =  atan2(s2,c2)  \boxed{\;q_2 \;=\; \operatorname{atan2}(s_2,\,c_2)\;}
Python · krok 5układ 2×2 dla q₂ + atan2 zachowuje ćwiartkę
# 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 (π,π](-\pi, \pi]. Klasyczny błąd implementacyjny: „policz arctan(s/c) i zgadnij ćwiartkę po znakach" — działa tylko po stronie poprawnie pamiętających autorów.

Krok 6

Przejście od pozycji do orientacji — macierz R₃⁶

Do tego momentu zajmowaliśmy się pozycją środka nadgarstka. Teraz przechodzimy do orientacji: trzy ostatnie przeguby q4,q5,q6q_4, q_5, q_6 muszą obrócić nadgarstek tak, aby pełna macierz R06R_0^6 zgodziła się z zadanym RR.

Część rotacyjną macierzy T03T_0^3 — czyli R03R_0^3 — wyliczamy ze znanych już q1,q2,q3q_1, q_2, q_3. Kluczowa obserwacja: w konwencji Craiga q2q_2 i q3q_3 obracają się wokół tej samej (poziomej) osi z^2=z^3\hat{z}_2 = \hat{z}_3, bo α2=0\alpha_2 = 0 w tabeli DH (między układami {2} i {3} nie ma skręcenia osi). Obroty wokół wspólnej osi zwyczajnie się dodają: Rz(q2)Rz(q3)=Rz(q2+q3)R_z(q_2)\,R_z(q_3) = R_z(q_2 + q_3). Razem:

R03  =  Rz(q1)T01Rx(π/2)T12,  częsˊcˊ rotacyjnaRz(q2+q3)T23,  częsˊcˊ rotacyjnaR_0^3 \;=\; \underbrace{R_z(q_1)}_{T_0^1}\,\underbrace{R_x(-\pi/2)}_{T_1^2,\;\text{część rotacyjna}}\,\underbrace{R_z(q_2 + q_3)}_{T_2^3,\;\text{część rotacyjna}}
▸ Pokaż jawne mnożenie trzech macierzy 3×3

Krok 1 — pomnóżmy Rx(π/2)Rz(q2+q3)R_x(-\pi/2)\,R_z(q_2 + q_3):

[100001010][c23s230s23c230001]=[c23s230001s23c230]\begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & -1 & 0 \end{bmatrix} \begin{bmatrix} c_{23} & -s_{23} & 0 \\ s_{23} & c_{23} & 0 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} c_{23} & -s_{23} & 0 \\ 0 & 0 & 1 \\ -s_{23} & -c_{23} & 0 \end{bmatrix}

(mnożenie wiersz·kolumna; Rx(π/2)R_x(-\pi/2) ma cos(π/2)=0,sin(π/2)=1\cos(-\pi/2) = 0,\,\sin(-\pi/2) = -1, stąd druga i trzecia kolumna).

Krok 2 — pomnóżmy Rz(q1)R_z(q_1) z wynikiem powyżej:

[c1s10s1c10001][c23s230001s23c230]\begin{bmatrix} c_1 & -s_1 & 0 \\ s_1 & c_1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} c_{23} & -s_{23} & 0 \\ 0 & 0 & 1 \\ -s_{23} & -c_{23} & 0 \end{bmatrix}

Pierwszy wiersz: (c1c23+(s1)0+0(s23),  c1(s23)+(s1)0+0(c23),  c10+(s1)1+00)(c_1\cdot c_{23} + (-s_1)\cdot 0 + 0\cdot(-s_{23}),\;c_1\cdot(-s_{23}) + (-s_1)\cdot 0 + 0\cdot(-c_{23}),\;c_1\cdot 0 + (-s_1)\cdot 1 + 0\cdot 0).

Po uproszczeniu wszystkich trzech wierszy:

R03  =  [c1c23c1s23s1s1c23s1s23c1s23c230]R_0^3 \;=\; \begin{bmatrix} c_1\,c_{23} & -c_1\,s_{23} & -s_1 \\ s_1\,c_{23} & -s_1\,s_{23} & c_1 \\ -s_{23} & -c_{23} & 0 \end{bmatrix}

Sprawdzenie ortogonalności: każdy wiersz ma normę c12c232+c12s232+s12=c12+s12=1c_1^2 c_{23}^2 + c_1^2 s_{23}^2 + s_1^2 = c_1^2 + s_1^2 = 1, a iloczyn skalarny pierwszych dwóch wierszy: c1s1c232+c1s1s232s1c1=0c_1 s_1 c_{23}^2 + c_1 s_1 s_{23}^2 - s_1 c_1 = 0. Zgodnie z oczekiwaniem dla macierzy obrotu.

Wynik jawny (przepisany dla wygody — używamy go bezpośrednio w kodzie):

R03  =  [c1c23c1s23s1s1c23s1s23c1s23c230]R_0^3 \;=\; \begin{bmatrix} c_1\,c_{23} & -c_1\,s_{23} & -s_1 \\ s_1\,c_{23} & -s_1\,s_{23} & c_1 \\ -s_{23} & -c_{23} & 0 \end{bmatrix}

Dalej: jeśli R=R06=R03R36R = R_0^6 = R_0^3 \cdot R_3^6, to nadgarstek musi dostarczyć „brakującą" rotację:

R36  =  (R03)1R  =  (R03) ⁣RR_3^6 \;=\; (R_0^3)^{-1}\,R \;=\; (R_0^3)^{\!\top}\,R

Równość (R03)1=(R03)(R_0^3)^{-1} = (R_0^3)^{\top} wynika z ortogonalności macierzy rotacji (wiersze są ortonormalne). W implementacji liczymy więc po prostu iloczyn transponowanej macierzy przez RR — bez wywoływania numerycznej inwersji.

Python · krok 6R₀³ wyliczone z (q₁, q₂, q₃), R₃⁶ = R₀³ᵀ · R
# 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 inv
Krok 7

Ekstrakcja q₄, q₅, q₆ z macierzy R₃⁶

Trzy ostatnie ogniwa Pumy to klasyczny nadgarstek typu ZYZ: dwa skręcenia osi ±π/2\pm\pi/2 między kolejnymi przegubami sprawiają, że osie z^4\hat z_4, z^5\hat z_5, z^6\hat z_6 są wzajemnie prostopadłe i schodzą się w jednym punkcie (środek nadgarstka). Z tabeli DH (α3=π/2\alpha_3 = -\pi/2, α4=π/2\alpha_4 = \pi/2, α5=π/2\alpha_5 = -\pi/2), iloczyn:

R36=Rz(q4)Rx(π/2)Rz(q5)Rx(π/2)Rz(q6)R_3^6 = R_z(q_4)\,R_x(-\pi/2)\,R_z(q_5)\,R_x(\pi/2)\,R_z(q_6)

Skąd dokładnie te pięć macierzy? — rozbicie iloczynu na cegiełki DH

Część rotacyjna każdej macierzy DH (Craig) ma postać Rx(αi1)Rz(θi)R_x(\alpha_{i-1})\,R_z(\theta_i) — najpierw skręcamy poprzednią oś z^\hat z o stały kąt αi1\alpha_{i-1} (geometria ogniwa), potem obracamy o zmienną θi=qi\theta_i = q_i wokół nowej osi z^i\hat z_i (kąt przegubu). Sklejając trzy ostatnie ogniwa, dostajemy:

R36  =  Rz(q4)przegub 4Rx(α4)Rz(q5)ogniwo 4\to5Rx(α5)Rz(q6)ogniwo 5\to6R_3^6 \;=\; \underbrace{R_z(q_4)}_{\text{przegub 4}}\,\underbrace{R_x(\alpha_4)\,R_z(q_5)}_{\text{ogniwo 4{\to}5}}\,\underbrace{R_x(\alpha_5)\,R_z(q_6)}_{\text{ogniwo 5{\to}6}}

Pierwsza cegiełka to czyste Rz(q4)R_z(q_4) bo skręcenie α3=π/2\alpha_3 = -\pi/2 zostało już „pochłonięte" w macierzy R03R_0^3. Pozostałe dwie pary biorą skręcenia z tabeli DH — α4=+π/2\alpha_4 = +\pi/2 (oś z^4\hat z_4 obraca się o 90° aby pokryć się z osią z^5\hat z_5 prostopadłą do przedramienia) oraz α5=π/2\alpha_5 = -\pi/2 (oś z^5\hat z_5 wraca o 90° w przeciwną stronę, aby oś z^6\hat z_6 wskazywała wzdłuż narzędzia). Stąd „kanapka" z–x–z–x–z: trzy obroty przegubowe Rz(q4),Rz(q5),Rz(q6)R_z(q_4), R_z(q_5), R_z(q_6) przeplatane dwoma stałymi „skręceniami montażowymi" Rx(±π/2)R_x(\pm\pi/2).

Dlaczego to się nazywa „typ ZYZ"? Bo gdy odpowiednio przenieść stałe macierze Rx(±π/2)R_x(\pm\pi/2) przez obroty RzR_z (korzystając z tożsamości Rx(π/2)Rz(θ)Rx(π/2)=Ry(θ)R_x(\pi/2)\,R_z(\theta)\,R_x(-\pi/2) = R_y(\theta)), cała kanapka redukuje się do Rz(q4)Ry(q5)Rz(q6+const)R_z(q_4)\,R_y(q_5)\,R_z(q_6 + \mathrm{const}) — klasycznych kątów Eulera ZYZ. Konwencja DH wymusza zapis przez RzR_z i RxR_x, 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 q5q_5 obraca wokół osi prostopadłej do osi q4q_4 i q6q_6, więc środkowy element całej macierzy to cosq5\cos q_5, a środkowy wiersz/kolumna kodują q4±q6q_4 \pm q_6. Po wymnożeniu:

R36=[c4c5c6s4s6c4c5s6s4c6c4s5s5c6s5s6c5s4c5c6c4s6s4c5s6c4c6s4s5]R_3^6 = \begin{bmatrix} c_4 c_5 c_6 - s_4 s_6 & -c_4 c_5 s_6 - s_4 c_6 & -c_4 s_5 \\ s_5 c_6 & -s_5 s_6 & c_5 \\ -s_4 c_5 c_6 - c_4 s_6 & s_4 c_5 s_6 - c_4 c_6 & s_4 s_5 \end{bmatrix}

Element środkowego wiersza R36[1][2]=c5R_3^6[1][2] = c_5 daje nam q5q_5 z dokładnością do znaku sinq5\sin q_5. Pozostałe elementy średniej kolumny i wiersza rozszyfrowują q4q_4 i q6q_6:

c5=R36[1][2],s5cosq6=R36[1][0],s5sinq6=R36[1][1]c_5 = R_3^6[1][2], \qquad s_5\cos q_6 = R_3^6[1][0], \qquad -\,s_5\sin q_6 = R_3^6[1][1]
c4s5=R36[0][2],s4s5=R36[2][2]-\,c_4\,s_5 = R_3^6[0][2], \qquad s_4\,s_5 = R_3^6[2][2]

Ponieważ sinq5\sin q_5 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 sinq50\sin q_5 \neq 0:

q5  =  atan2 ⁣(±R36[1][0]2+R36[1][1]2,  R36[1][2])q_5 \;=\; \operatorname{atan2}\!\left(\pm\sqrt{R_3^6[1][0]^2 + R_3^6[1][1]^2},\;R_3^6[1][2]\right)

Skąd ten pierwiastek? — jedynka trygonometryczna w działaniu

Mamy już cosq5=R36[1][2]\cos q_5 = R_3^6[1][2] — to wystarczyłoby do policzenia q5q_5 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±1\pm 1 (małe błędy zaokrąglenia → duże błędy kąta). Dlatego rekonstruujemy także sinq5\sin q_5.

Z dwóch równości z pierwszej kolumny środkowego wiersza:

R36[1][0]=s5c6,R36[1][1]=s5s6R_3^6[1][0] = s_5\,c_6, \qquad R_3^6[1][1] = -s_5\,s_6

Podnosimy do kwadratu i dodajemy, korzystając z jedynki trygonometrycznej c62+s62=1c_6^2 + s_6^2 = 1:

R36[1][0]2+R36[1][1]2  =  s52c62+s52s62  =  s52(c62+s62)  =  s52R_3^6[1][0]^2 + R_3^6[1][1]^2 \;=\; s_5^2\,c_6^2 + s_5^2\,s_6^2 \;=\; s_5^2\,(c_6^2 + s_6^2) \;=\; s_5^2

Stąd sinq5=±R36[1][0]2+R36[1][1]2\sin q_5 = \pm\sqrt{R_3^6[1][0]^2 + R_3^6[1][1]^2} znak ±\pm jest fundamentalnie nierozstrzygalny z tych równań (układ trzech równań na trzy niewiadome ma dwie gałęzie rozwiązań). Mając parę (sinq5,cosq5)(\sin q_5, \cos q_5), składamy:

q5  =  atan2(sinq5,  cosq5)  =  atan2 ⁣(±R36[1][0]2+R36[1][1]2,  R36[1][2])q_5 \;=\; \operatorname{atan2}(\sin q_5,\;\cos q_5) \;=\; \operatorname{atan2}\!\left(\pm\sqrt{R_3^6[1][0]^2 + R_3^6[1][1]^2},\;R_3^6[1][2]\right)

Bonus geometryczny: w nadgarstku typu ZYZ  sinq5\;|\sin q_5| to długość rzutu osi narzędziaz^6\hat z_6 na płaszczyznę prostopadłą do z^4\hat z_4. Gdy ta długość maleje do zera, osie z^4\hat z_4 i z^6\hat z_6 stają się współliniowe — pojawia się osobliwość nadgarstka opisana niżej.

q4  =  atan2(±R36[2][2],  R36[0][2])q_4 \;=\; \operatorname{atan2}(\pm R_3^6[2][2],\;\mp R_3^6[0][2])
q6  =  atan2(R36[1][1],  ±R36[1][0])q_6 \;=\; \operatorname{atan2}(\mp R_3^6[1][1],\;\pm R_3^6[1][0])

(górne znaki — gałąź no-flip, dolne — flip).

Osobliwość (singularność) nadgarstka (q₅ ≈ 0)

Gdy sinq50\sin q_5 \to 0, dwa końcowe obroty nadgarstka zachodzą wokół tej samej osi. Matematycznie wyzerowuje się cała średnia kolumna i średni wiersz (poza c5=±1c_5 = \pm 1), a w macierzy R36R_3^6 pozostaje tylko rotacja sumaryczna q4+q6q_4 + q_6. Tracimy stopień swobody: pojedynczy kąt q4q_4 albo q6q_6 z osobna jest nieokreślony, wyznaczona jest tylko ich suma.

W kodzie detekujemy to przez obliczenie sinq5=R36[1][0]2+R36[1][1]2|\sin q_5| = \sqrt{R_3^6[1][0]^2 + R_3^6[1][1]^2}. Gdy ta liczba jest mniejsza niż ustalona tolerancja, wybieramy umownie q4=0q_4 = 0 i resztę dopasowujemy tak, by orientacja docelowa była zachowana.

Python · krok 7dwie gałęzie wrist + obsługa singularności q₅ ≈ 0
# 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))
Krok 8

Osiem rozwiązań — enumeracja i selekcja praktyczna

W trzech miejscach podejmowaliśmy decyzję „znak ++ czy -":

  • shoulder — znak ρ\rho (krok 2),
  • elbow — znak L2K2\sqrt{L^2 - K^2} (krok 4),
  • wrist — znak sinq5\sin q_5 (krok 7).

Trzy niezależne wybory binarne dają 23=82^3 = 8 kombinacji — dlatego typowo dla Pumy mamy do ośmiu różnych konfiguracji qq 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:

  • Najbliższy aktualnej konfiguracji — minimalizacja ruchu przegubów, z uwzględnieniem zawinięcia kątów (bo qiq_i i qi+2πq_i + 2\pi to mechanicznie różne położenia, a matematycznie ten sam kąt).
  • W granicach przegubowych — niektóre teoretyczne rozwiązania mogą być fizycznie nieosiągalne.
  • Z dala od singularności — wartość manipulacyjności (patrz moduł 7) jako dodatkowe kryterium.
  • Bez kolizji — sprawdzenie zewnętrznym kolizją checker.

Poniżej miniatury wszystkich rozwiązań dla aktualnej pozy TT^*. 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.

Wszystkie rozwiązania: 8

kliknij, aby wczytać do głównego kontrolera

Wzorzec liczbowy — pełny rachunek dla jednej pozy

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.

Poza testowa

pozycja TCP:(px, py, pz) = (0.5, 0.15, 0.3) morientacja (rpy):(r, p, y) = (0°, 90°, 0°)macierz R = Rz·Ry·Rx:[ 0 0 1 ] [ 0 1 0 ] [-1 0 0 ]

Krok 1–2 · q₁ z dwóch atan2

φ= atan2(p_y, p_x)16.6992°ρ= √(p_x² + p_y² − d₃²)0.5067 mq₁= φ − atan2(d₃, ρ)2.7996°

Krok 3 · płaszczyzna ramienia

L= √(a₃² + d₄²)0.4323 mβ= atan2(d₄, a₃)87.3084°D= √(ρ² + p_z²)0.5889 m

Krok 4 · q₃ z prawa cosinusów

cos γ= (a₂² + L² − D²)/(2 a₂ L)0.0711γ= acos(cos γ)85.9233°K= (D² − a₂² − L²)/(2 a₂)-0.0307L² − K²0.1859q₃ (elbow ↑)= atan2(+√…, K) − β6.7684°

Krok 5 · q₂ z układu 2×2

M= a₂ + a₃c₃ − d₄s₃0.4011N= a₃s₃ + d₄c₃0.4312M² + N²0.3468c₂= (Mρ − Np_z)/(M²+N²)0.2130s₂= (−Mp_z − Nρ)/(M²+N²)-0.9770q₂= atan2(s₂, c₂)-77.6992°

Krok 7 · q₄, q₅, q₆ z R₃⁶ (gałąź wrist=flip)

q₄8.5126°q₅-19.2660°q₆171.9577°

Suma · pełna konfiguracja q

q = (2.80°, -77.70°, 6.77°, 8.51°, -19.27°, 171.96°)Sprawdzenie: FK(q) odtwarza pozycję z dokładnością ≈ 6·10⁻¹⁷ m. Wszystkie qi w granicach przegubów. Pozostałe siedem gałęzi (przy istnieniu) odtwarza tę samą pozę efektora innymi kątami przegubów.

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: a2=0,4318a_2 = 0{,}4318, a3=0,0203a_3 = 0{,}0203, d3=0,1254d_3 = 0{,}1254, d4=0,4318d_4 = 0{,}4318 m.

Klucz przewodni — pomysły algebraiczne, nie wzory

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.

StawGłówny pomysłTechnika algebraiczna
q1q_1Rzut na płaszczyznę poziomą eliminuje wszystkie wyższe stawyDwa atan2 w różnicy — kąt do celu w XY minus kąt korygujący offset d3d_3.
q3q_3Odległość bark↔wrist jest niezmiennikiem obrotów q1,q2q_1, q_2 — zależy tylko od q3q_3.Prawo cosinusów na trójkącie bark-łokieć-wrist; jedno równanie skalarne K=f(q3)K = f(q_3).
q2q_2Po wyznaczeniu q3q_3 układ pozycji staje się liniowy w (c2,s2c_2, s_2).Pogrupowanie wyrazów pozycji wg c2c_2 i s2s_2 daje układ liniowy 2×2; reguła Cramera zwraca obie wartości i atan2(s_2, c_2) dopina kąt.
q4,q5,q6q_4, q_5, q_6Mając q1,q2,q3q_1, q_2, q_3, „brakująca" rotacja R36R_3^6 jest jawnie znana — można z niej odczytać trzy kąty Eulera.R36=(R03) ⁣RR_3^6 = (R_0^3)^{\!\top}\,R, 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 q1q_1 tym kluczem jest niezależność rzutu poziomego od q2,q3q_2, q_3; dla q3q_3 — niezmiennicza długość przekątnej; dla q2q_2 — liniowość po wyłuskaniu q3q_3; dla orientacji — rozkład grupowy R06=R03R36R_0^6 = R_0^3 \cdot R_3^6. 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.

Kompletna funkcja Python — wszystkie kroki sklejone

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.

Python · pełna implementacjakopiuj-wklej do notatnika i uruchom — działa na czystym NumPy
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.

Ściąga formuł

Wszystkie kluczowe wzory algorytmu IK Pumy 560 zebrane w jednym miejscu — przydatne jako kompaktowa referencja przy implementacji albo powtórce przed egzaminem.

Stałe DH (Craig) Pumy 560
a2=0,4318a_2 = 0{,}4318 ma3=0,0203a_3 = 0{,}0203 md3=0,1254d_3 = 0{,}1254 md4=0,4318d_4 = 0{,}4318 m
Krok 0 · Odsprzężenie pozycji od orientacji
T06=TTtool1,pwc=pd6Rz^=d6=0pT_0^{6} = T^* \cdot T_\mathrm{tool}^{-1}, \qquad \mathbf{p}_\mathrm{wc} = \mathbf{p} - d_6\,R\,\hat{\mathbf{z}} \overset{d_6=0}{=} \mathbf{p}
Krok 1 · Wzory na pozycję środka nadgarstka
ρ=a2c2+a3c23d4s23,px2+py2=ρ2+d32\rho = a_2\,c_2 + a_3\,c_{23} - d_4\,s_{23}, \qquad p_x^2 + p_y^2 = \rho^2 + d_3^2
Krok 2 · q₁ — gałąź barku (shoulder)
ρ=±px2+py2d32,q1=atan2(py,px)atan2(d3,ρ)\rho = \pm\sqrt{p_x^2 + p_y^2 - d_3^2}, \qquad q_1 = \operatorname{atan2}(p_y, p_x) - \operatorname{atan2}(d_3, \rho)
Krok 3 · stałe płaszczyzny ramienia
L=a32+d42,β=atan2(d4,a3),D2=ρ2+pz2L = \sqrt{a_3^2 + d_4^2}, \qquad \beta = \operatorname{atan2}(d_4, a_3), \qquad D^2 = \rho^2 + p_z^2
Krok 4 · q₃ — gałąź łokcia (elbow), prawo cosinusów
K=D2a22L22a2=a3c3d4s3    q3=atan2 ⁣(±L2K2,  K)βK = \frac{D^2 - a_2^2 - L^2}{2\,a_2} = a_3\,c_3 - d_4\,s_3 \;\Rightarrow\; q_3 = \operatorname{atan2}\!\bigl(\pm\sqrt{L^2 - K^2},\;K\bigr) - \beta
Krok 5 · q₂ z układu 2×2 w (c₂, s₂)
M=a2+a3c3d4s3,N=a3s3+d4c3M = a_2 + a_3\,c_3 - d_4\,s_3, \qquad N = a_3\,s_3 + d_4\,c_3
c2=MρNpzM2+N2,s2=MpzNρM2+N2,q2=atan2(s2,c2)c_2 = \frac{M\rho - N p_z}{M^2 + N^2}, \quad s_2 = \frac{-M p_z - N\rho}{M^2 + N^2}, \quad q_2 = \operatorname{atan2}(s_2, c_2)
Krok 6 · macierz R₀³ z gotowych kątów q₁, q₂, q₃
R03=[c1c23c1s23s1s1c23s1s23c1s23c230],R36=(R03) ⁣RR_0^3 = \begin{bmatrix} c_1 c_{23} & -c_1 s_{23} & -s_1 \\ s_1 c_{23} & -s_1 s_{23} & c_1 \\ -s_{23} & -c_{23} & 0 \end{bmatrix}, \qquad R_3^6 = (R_0^3)^{\!\top}\,R
Krok 7 · q₄, q₅, q₆ z elementów R₃⁶ (gałąź wrist)
q5=atan2 ⁣(±R10362+R11362,  R1236)q_5 = \operatorname{atan2}\!\bigl(\pm\sqrt{R^{36}_{10}{}^2 + R^{36}_{11}{}^2},\;R^{36}_{12}\bigr)
q4=atan2(±R2236,  R0236),q6=atan2(R1136,  ±R1036)q_4 = \operatorname{atan2}(\pm R^{36}_{22},\;\mp R^{36}_{02}), \qquad q_6 = \operatorname{atan2}(\mp R^{36}_{11},\;\pm R^{36}_{10})
Krok 8 · 8 gałęzi = 2 (shoulder) × 2 (elbow) × 2 (wrist)

Wybór znaków w trzech miejscach (ρ\rho, dyskryminant L2K2\sqrt{L^2 - K^2}, sinq5\sin q_5) generuje 23=82^3 = 8 kombinacji. Część może odpaść jako poza zasięgiem (L2<K2L^2 < K^2) lub poza limitami przegubów; w singularności nadgarstka (sinq50\sin q_5 \approx 0) gałęzie wrist zlewają się.

Pułapki — częste źródła błędów
  • Zapomnienie drugiego członu w q₁: atan2(p_y, p_x) bez − atan2(d₃, ρ) daje błąd ~7° wynikający z odsadzenia d3d_3.
  • Użycie arcsin lub arccos tam gdzie wystarcza atan2 — gubi informację o ćwiartce, błąd zniknie tylko dla niektórych pozycji.
  • Mieszanie konwencji DH klasycznej i Craiga — liczby α,a,d\alpha, a, d są te same, ale przesunięte o jeden indeks. Wzory w tym walkthrough zakładają Craiga.
  • Numeryczna inwersja (R03)1(R_0^3)^{-1} zamiast transpozycji — kosztownie i z gorszą precyzją. Macierze rotacji są ortogonalne: R1=RR^{-1} = R^{\top}.
  • Brak detekcji singularności sinq5<ε|\sin q_5| < \varepsilon przed dzieleniem — kończy się wartościami NaN dla q₄ i q₆.

Osobliwości (singularności) w jednym miejscu

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ściWarunek matematycznyCo się dzieje fizycznieKonsekwencja dla IK
Osobliwość osi 1 (nad bazą)px2+py2=d32p_x^2 + p_y^2 = d_3^2 (cel na granicy „cylindra zakazanego")TCP znajduje się dokładnie nad pionową osią bazy — obrót q1q_1 wokół tej osi nie zmienia położenia TCP.q1q_1 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ętrznaL2=K2L^2 = K^2, czyli D=a2±LD = |a_2 \pm L| — przekątna bark↔wrist osiąga ekstremum.Ramię jest maksymalnie wyprostowane (D=a2+LD = a_2 + L) lub maksymalnie zgięte (D=a2LD = |a_2 - L|). Trójkąt bark-łokieć-wrist degeneruje się do odcinka.Dyskryminant L2K2=0L^2 - K^2 = 0 — 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)L2K2<0L^2 - K^2 < 0Cel leży poza fizycznym zasięgiem ramienia — trójkąt o bokach a2,L,Da_2, L, D nie spełnia nierówności trójkąta (suma dwóch krótszych boków jest mniejsza od trzeciego).ujemna liczba\sqrt{\text{ujemna liczba}} → 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)sinq5=0\sin q_5 = 0, czyli q5{0,π}q_5 \in \{0, \pi\}.Osie z^4\hat z_4 i z^6\hat z_6 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 q4+q6q_4 + q_6 (przy q5=0q_5 = 0) lub różnica q4q6q_4 - q_6 (przy q5=πq_5 = \pi) jest wyznaczona; sam q4q_4 i q6q_6 mają nieskończenie wiele par. Trzeba przyjąć umowę (np. q4=0q_4 = 0) i dopasować q6q_6 do reszty.

Wspólny mianownik wszystkich czterech przypadków: dwie kolumny Jakobianu J(q)J(q) stają się liniowo zależne, więc wyznacznik detJ=0\det J = 0. Geometrycznie znaczy to, że istnieje kierunek w przestrzeni SE(3)SE(3) (kierunek prędkości efektora), w którym żadna kombinacja q˙1,,q˙6\dot q_1, \dots, \dot q_6 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.

Ten sam solver, dwa języki

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).

Porównanie uruchomień: TypeScript vs Python

Python: nie załadowany

Python runtime (Pyodide) uruchamiany jest w web workerze i wymaga jednorazowego pobrania ~10 MB zasobów. Kliknij, by go zainicjować.

TypeScript (natywnie)

środowisko jeszcze się ładuje…

Python 3.12 (Pyodide + NumPy)

środowisko jeszcze się ładuje…

Co dalej

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ż qq 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.