肌肉骨骼模型
骨骼的剛體動力學模型
首先基于剛體動力學對于人體骨骼系統進行模組化。首先确定系統的廣義坐标 q q q，對于圖中所示的二維模型為例，廣義坐标包含軀幹的水準、數值方向位移和姿态角，以及每條腿三個關節角，一共 3 + 2 ∗ 3 = 9 3+2*3=9 3+2∗3=9個自由度。之後通過這些廣義變量和廣義變量的導數，表示出系統的動能和勢能，并構造拉格朗日函數：
L = T − V L=T-V L=T−V
考慮骨骼系統僅收到地反力和肌肉力的作用，通過拉格朗日方程得到系統的動力學模型：
d d t ( ∂ L ∂ q ˙ j ) − ∂ L ∂ q j = τ c o n t a c t j + τ m u s c l e j \frac{\mathrm{d}}{\mathrm{d} t}\left(\frac{\partial L}{\partial \dot{q}_{j}}\right)-\frac{\partial L}{\partial q_{j}}=\tau^{j}_{contact}+\tau^{j}_{muscle} dtd(∂q˙j∂L)−∂qj∂L=τcontactj+τmusclej
M ( q ) ⋅ q ¨ + B ( q , q ˙ ) = J ( q ) T F c o n t a c t + R ( q ) F m u s c l e \mathbf{M}(\mathbf{q}) \cdot \ddot{\mathbf{q}}+\mathbf{B}(\mathbf{q}, \dot{\mathbf{q}})=\mathbf{J(q)}^T\mathbf{F}_{contact}+\mathbf{R(q)}\mathbf{F}_{muscle} M(q)⋅q¨+B(q,q˙)=J(q)TFcontact+R(q)Fmuscle
其中 M ( q ) \mathbf{M(q)} M(q)為品質矩陣， B ( q , q ˙ ) \mathbf{B}(\mathbf{q}, \dot{\mathbf{q}}) B(q,q˙)包含重力、離心力和科氏力， J ( q ) \mathbf{J(q)} J(q)表示雅克比矩陣，将地反力轉化為對每個廣義坐标的等效力矩， R ( q ) \mathbf{R(q)} R(q)表示肌肉力對廣義坐标的等效力矩。下面我們主要讨論等式右側的肌肉力與地反力。
肌肉動力學
肌肉激活動力學
盡管各個文獻裡的描述不盡一緻，但一般來說excitation用來描述神經信号，是肌肉的控制輸入（整流濾波歸一化的EMG信号就是這個），而activation用來表示肌肉的活躍度，是肌肉的内部狀态。在靜态狀态下，這兩者可能相等，但在動态變化中，兩者一般不同，存在由excitation到activation動态過程，即肌肉激活動力學。各方學者建立了不同的激活動力學模型，一般包含一個微分方程和一些非線性變換，其中Winters開發的一種激活動力學如下：
f = 0.5 tanh ( b ( e − a ) ) f=0.5 \tanh (b(e-a)) f=0.5tanh(b(e−a))
d a d t = [ 1 τ a ( 0.5 + 1.5 a ) ( f + 0.5 ) + 0.5 + 1.5 a τ d ( − f + 0.5 ) ] ( e − a ) \frac{d a}{d t}=\left[\frac{1}{\tau_{a}(0.5+1.5 a)}(f+0.5)+\frac{0.5+1.5 a}{\tau_{d}}(-f+0.5)\right](e-a) dtda=[τa(0.5+1.5a)1(f+0.5)+τd0.5+1.5a(−f+0.5)](e−a)
肌肉收縮動力學
肌肉收縮的特性可以用如上的Hill肌肉模型來描述，它由肌肉收縮單元CE與被動彈性單元PE并聯，并與一個肌腱彈性單元T串聯，肌肉與肌腱之間還存在一個角度成為羽狀角。由于根據上面的模型，肌肉實際産生的力不僅受到收縮單元CE的激活度的影響，還受到肌肉的長度、收縮速度的影響。由肌肉纖維産生的力具有如下關系：
F m t ( t ) = F t = F max [ a ( t ) f ( l ) f ( v ) + f p ( l ) ] cos ( ϕ ( t ) ) \begin{aligned} F^{\mathrm{mt}}(t) &=F^{\mathrm{t}} \\ &=F^{\max }\left[a(t)f(l) f(v) +f_{\mathrm{p}}(l)\right] \cos (\phi(t)) \end{aligned} Fmt(t)=Ft=Fmax[a(t)f(l)f(v)+fp(l)]cos(ϕ(t))
F max F^{\max } Fmax表示最大等長收縮力， f ( l ) , f ( v ) , f ( p ) f(l),f(v),f(p) f(l),f(v),f(p)分别表示力-長度、力-速度、被動力-長度的關系。這些函數由下圖所示的實驗曲線所描述，一般可以通過多項式函數拟合來近似表示。
注意到上式公式中，實際輸入為激活度 a a a、歸一化肌肉長度 l l l、和歸一化收縮速度 v v v，而肌肉長度與收縮速度之間存在關系，收縮動力學的狀态量僅有 a a a和 l l l，且上述方程描述的是一個微分方程。我們将整個肌肉動力學表示為如下形式，但這種表示并不嚴謹，實際上每塊肌肉包含着兩個微分方程：
F m u s c l e = F m u s c l e ( q , a , a ˙ , L C E , L ˙ C E , e ) \mathbf{F}_{muscle} = \mathbf{F}_{muscle}(\mathbf{q},\mathbf{a} , \mathbf{\dot a},\mathbf{L}_{CE},\mathbf{\dot L}_{CE},\mathbf{e}) Fmuscle=Fmuscle(q,a,a˙,LCE,L˙CE,e)
地反力模型
OpenSim中采用Hunt-Crossley模型來表示接觸力，接觸力由剛度項、阻尼項和摩擦項組合而成，其表達式為：
f contact = f stiffness + f dissipation + f friction \mathbf{f}_{\text {contact}}=\mathbf{f}_{\text {stiffness}}+\mathbf{f}_{\text {dissipation}}+\mathbf{f}_{\text {friction}} fcontact=fstiffness+fdissipation+ffriction
接觸力模型中的第一項表示彈性力。如上圖所示，我們可以用兩個彈性小球來描述接觸，當小球接觸時會發生彈性變形，變形的程度與各自的彈性模量 E E E有關，總的變形為 x x x。當變形發生時，每個小球都會受到一個彈性力：
f stiffness = f H z = ( 4 3 σ R 1 / 2 E ∗ ) x 3 / 2 f_{\text {stiffness}}=f_{H z}=\left(\frac{4}{3} \sigma R^{1 / 2} E^{*}\right) x^{3 / 2} fstiffness=fHz=(34σR1/2E∗)x3/2
R R R表示合成相對曲率， E ∗ E^* E∗表示合成彈性模量，它們都可以通過兩個小球的參數計算出來； σ \sigma σ表示偏心系數，當 σ = 1 \sigma=1 σ=1的時候為圓接觸， x x x為總的變形量。
接觸力模型中的第二項表示損耗力，或是阻尼力。碰撞的發生往往伴随着能量的損耗，單純的彈性力無法表示這一點，是以引入Hunt-Crossley耗散力：
f d i s s i p a t i o n = 3 2 f s t i f f n e s s c ∗ x ˙ f_{dissipation}=\frac{3}{2} f_{stiffness} c^{*} \dot{x} fdissipation=23fstiffnessc∗x˙
c ∗ c^{*} c∗表示等效耗散系數，而 x ˙ \dot{x} x˙表現出耗散力的阻尼特性。剛度力和耗散力都垂直接觸面。
f normal = f stiffness + f dissipation f_{\text {normal}}=f_{\text {stiffness}}+f_{\text {dissipation}} fnormal=fstiffness+fdissipation
接觸力模型中的最後一項表示摩擦力，方向平行接觸面：
f friction = μ ( v ) f normal f_{\text {friction}}=\mu(v) f_{\text {normal}} ffriction=μ(v)fnormal
μ ( v ) \mu(v) μ(v)表示等效摩擦因數，是關于相對運動速度的函數。
在使用上述的Hunt-Crossley模型來計算地反力時，地面的剛度認為無窮大，則相對彈性變形 x x x可以通過接觸小球圓心相對地面的距離計算得到，即：
F c o n t a c t = F c o n t a c t ( q , q ˙ ) \mathbf{F}_{contact} = \mathbf{F}_{contact}(\mathbf{q}, \dot{\mathbf{q}}) Fcontact=Fcontact(q,q˙)
綜合模型
綜上所述，肌骨模型的動力學可以統一表示為：
M ( q ) ⋅ q ¨ + B ( q , q ˙ ) = J ( q ) T F c o n t a c t ( q , q ˙ ) + R ( q ) F m u s c l e ( q , a , a ˙ , L C E , L ˙ C E , e ) \mathbf{M}(\mathbf{q}) \cdot \ddot{\mathbf{q}}+\mathbf{B}(\mathbf{q}, \dot{\mathbf{q}})=\mathbf{J(q)}^T \mathbf{F}_{contact}(\mathbf{q}, \dot{\mathbf{q}})+\mathbf{R(q)}\mathbf{F}_{muscle}(\mathbf{q},\mathbf{a} , \mathbf{\dot a},\mathbf{L}_{CE},\mathbf{\dot L}_{CE},\mathbf{e}) M(q)⋅q¨+B(q,q˙)=J(q)TFcontact(q,q˙)+R(q)Fmuscle(q,a,a˙,LCE,L˙CE,e)
定義狀态變量：
x = ( q , q ˙ , L C E , a ) T \mathbf{x}=\left(\mathbf{q}, \dot{\mathbf{q}}, \mathbf{L}_{C E}, \mathbf{a}\right)^{T} x=(q,q˙,LCE,a)T
在定義 u u u為16塊肌肉的控制量（excitations），上面的動力學方程可表示為：
f ( x , x ˙ , u ) = 0 \mathbf{f}(\mathbf{x}, \dot{\mathbf{x}}, \mathbf{u})=0 f(x,x˙,u)=0
肌骨模型的仿真應用
前向仿真
逆向分析
逆向分析是指，在已有運動學、地反力資料後，根據實驗資料推斷出産生相應運動的肌肉活動，也就是OpenSim中的CMC。由于肌肉之間存在和比較嚴重的耦合問題，多塊肌肉連接配接同一個關節，同一塊肌肉連接配接着不同關節，是以逆向分析最重要的是解決肌肉備援 - Solving the Muscle Redundancy Problem。
在這個問題中，已有運動學資料 q k \mathbf{q}_k qk和地反力 F C k \mathbf{F}_{\mathbf{C}k} FCk序列， k = 1 , … , K k=1, \ldots, K k=1,…,K表示離散時間。我們先計算逆動力學得到關節力矩 T I D k T_{\mathbf{ID}k} TIDk，關節力矩應當等于有肌肉收縮的等效力矩：
T I D k = ∑ m = 1 M r m k F T m + e T k T m a x T_{\mathrm{ID} k}=\sum_{m=1}^{M} r_{m k} F_{\mathrm{T} m}+e_{\mathrm{T} k} T_{\mathrm{max}} TIDk=m=1∑MrmkFTm+eTkTmax
m m m表示肌肉序号， r r r表示肌肉對關節等效力矩。有時候關節力矩和肌肉裡橘無法嚴格相等，是以我們添加了公式中的最後一項來表示力矩的殘差，而 e T k e_{\mathrm{T} k} eTk表示殘差系數，在後面我們會設計優化目标是殘差力盡可能的小。
由此建立逆向分析的優化架構：