2/15 MITの安定化制御則(5)


今回も[1]の内容を見ていく. ページで言うと, 18/47~26/47

4 Controller Design
フィードバック制御へのアプローチは次のようにまとめることができる.
動力学の縮小モデル(例えば、重心動力学)についての計画された軌道が与えられると, 軌道の周りの時変線形化を計算し, 時変線形二次レギュレータ(TV− LQR)設計を使用して安定化コントローラを導出する. しかし, LQRのclosed-form optimal controller(閉形式最適コントローラ)を使用するのではなく, フルウォーキングシステムの瞬間的なダイナミクス, 入力, および接触の制約を追加的に取り込む2次計画(QP)を定式化して解く.

closed-form optimal controllerについては[2]を参照

4.1 General Formulation

ダイナミック, \(\dot{\mathbf x}=f(\mathbf{ x, u})\), また, \(t\in [0,t_f]\)に対する目標の軌跡\(\mathbf x_d(t), \mathbf u_d(t)\)を考え, 次式のように軌跡に沿て線形化する.
$$\dot{\bar{\mathbf x}}(t)=\mathbf A(t)\bar{\mathbf x}(t)+\mathbf B(t)\bar{\mathbf u}(t)$$
相対座標を使用して, \(\bar x(t)=x(t)-x_d(t)\)および\(\bar u(t)= u(t) – u_d(t)\)である.
次に, 二次コスト関数を定義する.
$$g(\mathbf x(t), \mathbf u(t), t)=\bar{\mathbf x}^T(t)\mathbf Q\bar{\mathbf x}(t)+\bar{\mathbf u}^T(t)\mathbf R\bar{\mathbf u}(t)$$
そして制約付き最小化問題を解く.

この問題に対する解決策は, 連続Riccati微分方程式を解くことによって与えられ, それはoptimal cost-to-goをもたす[3].([3]はMITで使われているロボット工学の教材です. めちゃくちゃおすすめのサイトです)
$$J(\mathbf x,t)=\bar{\mathbf x}^T(t)\mathbf S(t)\bar{\mathbf x}(t)\mathbf s_1(t)+s_2(t)$$
Hamilton-Jacobi-Bellman (HJB) equation [4]より, the optimal controllerは次式を満足する.

すべての\(\bar{\mathbf x}\)について, 最小化から定数項を削除した. \(\bar{\mathbf u}^*(t)\)は, 単に\(\bar{\mathbf u}\)に関する\(l(\bar{\mathbf x}, \bar{\mathbf u}, t)\)の導関数を取り, それをゼロに等しく設定して解くことによって, 時変線形ポリシーとして書くことができる.
ここで重要な点は, 物理システムでは, トルク制限や地面反力への制約など, 無視できない不等式制約があることである. これらの制約条件が瞬時に線形不等式として表現できる限り, \(\bar{\mathbf x}\)で与えられる各制御ステップで入力を計算するために凸 QPを定式化して解くことができる.
minimize \(_{\bar{\mathbf u}}\:l(\bar{\mathbf x}, \bar{\mathbf u},t)\)
subject to  \(\mathbf M\bar{\mathbf u}\geq\mathbf b\hspace{20mm}(19)\)
言い換えれば, この最適化は, システムの制約に従って, cost-to-go subjectの最急降下方向を見つける.  このQPを解くことによって計算された入力は, 一般に, closed-form LQRポリシーの出力をしきい値処理して計算されたものと等しくないことに注意.  さらに、以下に示すように, このフレームワークを使用することで, 追加の制約を追加の計算コストをほとんどかけずにQPに追加することができる.
※cost-to-goについては[5][6]([5]は動的計画法を参照. wikipedia曰く, cost-to-goの訳は, 残余コスト)

4.2 COM and COP Stabilization
歩行計画が, 所望の足跡計画に基づいて導出された区分的多項式の圧力中心(COP)軌跡, \(\mathbf c_d(t),t\in [0,t_f]\)の形態である場合を考える. 現在の実装では, \(\mathbf c_d(t)\)は, スイング速度やダブルサポートタイムなどの高レベルの歩行パラメータによって管理されるタイミングで計画された足跡の重心間を補間する区分的線形軌道である.
局所的に平坦な地形上の脚付きシステムでは, 重心のダイナミクスはCOP[7],\(\mathbf c\)に関して次式のように書き直すことができる.
$$\dot{\mathbf k}=(\mathbf c-\mathbf r)\times\mathbf \lambda_C+\mathbf{\tau}_N$$
$$\dot{\mathbf l}=m\mathbf g+\mathbf \lambda_c$$
ここで, cはCOP, λcはCOPにおける正味の外力, τnは垂直接触モーメントである. ロボットの重心角運動量\(\dot{\mathbf k}=0, \mathbf k=0\),および法線モーメントτn = 0とすると, 重心ダイナミクスは次のように単純化される.
\(\dot{\mathbf x}_{CM}=\mathbf{Ax}_{CM}+\mathbf{Bu}_{CM}\equiv\begin{bmatrix}0_{2\times 2}&\mathbf I_{2\times 2}\\0_{2\times 2}&0_{2\times 2}\end{bmatrix}\mathbf x_{CM}+\begin{bmatrix}0_{2\times 2}\\ \mathbf I_{2\times 2}\end {bmatrix}\mathbf u_{CM}\)
$$\mathbf c=\begin{bmatrix}\mathbf I_{2\times 2}&0_{2\times 4}\end{bmatrix}\mathbf x_{CM}-\frac{r_Z}{\ddot r_Z+g}\mathbf I_{2\times 2}\mathbf u_{CM}$$
ここで, \(\mathbf x_{CM}=[r_x,r_y,\dot r_x,\dot r_y]^T\)はCOMのe ground projection, また\(\mathbf u_{CM}=[\ddot r_x,\ddot r_y]^T\)である.
出力\(\mathbf c\)は一般に非線形であるが, 2回微分可能なCOMの高さの軌跡を与えられると線形時変, またはCOMの高さが一定のままであると仮定すると線形時不変である(よく研究された線形倒立振子力学をもたらす).
この実装では, COMの高さを固定し, 線形形式の出力を使用する.
LQR問題を次のように定式化する.

ここで, \(\bar{\mathbf x}_{CM}=\mathbf x_{CM}(t)-[\mathbf c_d^T(t_f),0,0]^T\)また\(\bar{\mathbf u}_{CM}(t)=\mathbf u_{CM}(t)\)である. ここでコストは\(\bar{\mathbf c}(t)=\mathbf c(t)-\mathbf c_d(t)\)に関して次のように定義される.
\(g(\mathbf x_{CM}(t),\mathbf u_{CM}(t),t)=\bar{\mathbf c}^T(t)\mathbf Q_c\bar{\mathbf c}(t)+\bar{\mathbf u}_{CM}^T(t)\mathbf R\bar{\mathbf u}_{CM}(t)\)
ここで, 右辺は, 線形出力方程式に代入することによって\(x_{CM}(t),\bar{\mathbf u}_{CM}(t)\)および\(c_d(t)\)で表すことができる.  この問題についてRiccati方程式を導出すると, Sには時間依存項がないことがわかる( continuous algebraic Riccati equation:連続代数Riccati方程式に相当).
次に、線形および affine cost-to-go係数は, 軌跡に沿った単一のbackward passで(すなわち, 数値的に積分する必要なしに)解くことができる線形微分方程式になる. これは実際には\( l(\bar{\mathbf x}_{CM},\bar{\mathbf u}_{CM},t)\)は後退地平線の足跡計画に適した時間スケールで非常に効率的に計算できることを意味する.

4.3 Stabilizing whole-body plans
3.2節で説明した計画アルゴリズムは, 全身運動学, 重心運動量, および外力の軌跡の形で, より豊富な入力をコントローラに提供する. そのため, LQR問題の設計方法に柔軟性がある.  ただし, COP LQRの定式化の単純さに合わせて, 単純な線形COMダイナミクスを使用して3D COM軌道を安定させる. COMの軌跡\(r_{des}(t),\dot r_{des}(t), \ddot r_{des}(t)\)が与えられると, 2次コストは次式のように定義される.
\(g(\mathbf x_{CM}(t),\mathbf u_{CM}(t),t)=\bar{\mathbf x}^T(t)\mathbf Q_c\bar{\mathbf x}(t)+\bar{\mathbf u}_{CM}^T(t)\mathbf R\bar{\mathbf u}_{CM}(t)\)
ここで, \(\bar{\mathbf x}_{CM}(t)=\mathbf x_{CM}(t)-[\mathbf r_{des}^T(t),\dot{\mathbf r}_{des}^T(t)]^T, \bar{\mathbf u}_{CM}(t)=\mathbf u_{CM}(t)-\ddot {\mathbf r}_{des}(t)\)である.
ロボットの運動が空中段階を含む場合, ロボットのCOM運動が完全に重力によって決定されるという事実を説明するためにLQR問題を修正しなければならない. 空中の所望の運動を追跡するために, ロボットは離陸時に所望のCOM状態を達成しなければならない. このように, ロボットCOMのダイナミクスはハイブリッドになる. それは, サポートフェーズではスムーズなダブルインテグレータダイナミクスを持ち, 空中フェーズでは離陸から着陸状態への離散的な遷移を持つ.
ジャンプリカッチ方程式[8]
を使用して, このハイブリッドシステムのcost-to-goを計算する. 直感的には, このハイブリッドLQRアプローチは単純に離陸地点までおよび着陸後にCOM軌跡を追跡するという目標をエンコードする.

4.4 Additional costs and constraints

実際には, QP定式化(19)は, 満足のいく性能を達成するために追加のコスト項と線形の制約で増やされなければならない. まず最初に, 上で導出したLQRダイナミクスは完全なロボットダイナミクスではないので, ロボットへの実行可能な入力を計算するためにQPに追加の制約を追加する必要がある. ロボットの現在の状態\(\mathbf q, \mathbf v\)が与えられると, 運動方程式を計算することができ,
$$\mathbf H(\mathbf q)\dot{\mathbf v}+\mathbf C(\mathbf{q,v})=\mathbf{Bτ}+\mathbf J^T\mathbf λ$$
$$\begin{bmatrix}\mathbf H_f\\ \mathbf H_a\end{bmatrix}\dot{\mathbf v}+\begin{bmatrix}\mathbf C_f\\ \mathbf C_a\end{bmatrix}=\begin{bmatrix}\mathbf 0\\ \mathbf B_a\end{bmatrix}\mathbf τ+\begin{bmatrix}\mathbf J^T_f\\ \mathbf J^T_a\end{bmatrix}\mathbf λ$$
ここで, \(\mathbf H(\mathbf q)\)は質量行列, ベクトル\(\mathbf C(\mathbf q, \mathbf v)\)は重力とコリオリの項, \(\mathbf B\)はcontrol input map, \(\mathbf J^T\)は外力λを一般化力に変換する.  力学を明示的に非作動DOFと作動DOFに区分した.
典型的なヒューマノイドにとって, 唯一の動かされていない自由度は浮動ベースであり, \(\mathbf B_a\)はフルランクである.
この分解を使用して, 入力トルクτを決定変数として除去し, トルクの制約をλと\(\dot{\mathbf v}_υ\)の線形不等式制約として表すことができる. ロボットの歩行の場合, ベクトル\(\mathbf λ= [\mathbf λ^T_1…\mathbf λ^T_{N_c}]^T\)は\(N_c\)個の接触点に作用する接地力を含む.  接触力が摩擦コーンの内側にとどまるようにするために, 3.2で説明したのと同じ多面体近似(\(N_d = 4\))を使用する.
摩擦の制約を尊重することに加えて, コントローラは, 外力変数にゼロ以外の値を代入できる時期を検出する必要がある. 例えば, スイングフェイズ中またはスイングフェイズに入るとき, QPはその足に正の地面反力を割り当ててはいけない. このために, 最適な接触力変数を決定するために単純なロジックを使用した. ロボットが(足跡計画に基づいて)接触していると予想する場合, または足ひずみゲージが足に有意な値を報告している場合, その足のために正の値を割り当てることができる. これに対する例外は, 歩行中に足が接触を破っているときであり, そこでは計画は足が接触しているか否かを決定するために排他的に使用される. 実際には, 計画された接触状態がアクティブで力が測定されていない場合は, 摩擦係数を小さな一定値に設定する(事実上, その物体に通常の接触力のみを許可する), これは, 接触前の足の望ましくない接線方向の動きを減らすのに役立つ.
LQRソリューションはロボットの全身の動きを捉えないため, 足や骨盤などのロボットの身体に取り付けられた一連のフレームの望ましい加速度として追加のモーションゴールを指定する.  時間\(t\)におけるボディフレームの望ましい加速度\(\ddot{\mathbf p}(t)\)は,  PD規則を用いて次式のように計算される.
\(\ddot{\mathbf p}_{ref}(t)=\mathbf K_p(\mathbf p_d(t)-\mathbf p(t))+\mathbf K_d(\dot{\mathbf p}_d(\dot{\mathbf p}_d(t)-\dot{\mathbf p}(t))+\ddot{\mathbf p}_d(t))\)
ここで, \(\mathbf p_d(t)\)はフレームに対する滑らかな2回微分可能な所望の軌跡である. 典型的な歩行では, 滑らかな足の軌跡が足跡計画者によって計算され, 望ましい骨盤の動きが足跡計画から計算される. 所望の骨盤ヨーは平均足ヨーに等しく, 骨盤の高さは足の上の一定の高さに維持される.  全身の動的な動きについては, 計画で指定された滑らかな足, 骨盤, 胴体, および手の軌跡を追跡する.
身体加速度に加えて, またはその代わりに, 所望の一般化された加速度\(\dot{\mathbf v}_d\)を入力し, コストまたは制約としてQPに組み込むことができる. 例えば, 一般化された加速度拘束は, 歩行中および物体を運んでいる間に固定された上半身構成を維持するために有用である. 実際には, 制約のさまざまな組み合わせを試したり, 制約をコストに変更したりすること(およびその逆)をすばやく試すことができると非常に役立つ. [9]はこれらの変更を加えるためのさまざまなオプションをサポートする.
セクション6の実験で使用した完全なQPは, Quadratic Progmam 1.で説明される.
Quadratic Program 1.

ここで, 集合Cは環境と接触している身体の指標を含んでいる.  接触している物体に対しては, 「滑らない」拘束(”no slip” constraint)を適用する. それは, 例えば, 接触点の加速度が0またはそれらの速度と反対の方向にゼロ以外の値であることを要求し得る.  実行不可能性を回避するために, これらの制約の制限された違反を許容するスラック変数εを組み込んでいる(追加costを条件として). 二重サポートにおけるAtlasバランスの場合, QPには90個の決定変数, 30個の等式制約, および112個の不等式制約がある.  \(\bar{\mathbf x}_{CM}\)と\(\bar{\mathbf u}_{CM}\)は, ロボットの状態\((\mathbf q,\mathbf v)\)が与えられたときの決定変数のaffine 関数として次式のように表すことができる.
$$\bar{\mathbf x}_{CM}=\begin{bmatrix}COM(\mathbf q)\\ \mathbf J_{CM}(\mathbf q)\mathbf v\end{bmatrix}-\ddot{\mathbf x}_{CM,des}$$
$$\bar{\mathbf u}_{CM}=\dot{\mathbf J}_{CM}(\mathbf q)\mathbf v+\mathbf J_{CM}(\mathbf q)\dot{\mathbf v}-\bar{\mathbf u}_{CM,des}$$
ここで, \(\mathbf J_{CM}(\mathbf q)\)はCOM jacobianである.
シミュレーション[1011121314]とハードウェア[151617, 18]の両方でQPを使用して2足歩行システムを制御することを最近研究している研究者もいる.  我々の定式化のように, これらの最適化は通常, 重心(COM)加速度[10]、線形および角運動量の変化率[12, 16, 13], ゼロモーメント点(ZMP)[17], またはcanonical 歩行関数[15]などの低次元動的量を使用する.  ここで, QPコストは, 単純なPD制御規則を使用して計算された基準値からの重み付き距離からなる[10, 12, 16].
Johnson et al. [18]は, QPベースの運動量コントローラを実装し, Atlasロボットとのバランスをとりながら歩くことを実現し, 2013年のDRC試験で成功裏に競った.  Herzog et al. [16, 19]は, 油圧式ヒューマノイド下半身のバランス制御にQPコントローラを使用した.  我々のアプローチは, 特に制約の定義において, これらのアプローチといくつかの特徴を共有しているが, 軌跡を安定させQP目的を構築するための時変LQRの使用において異なる.
※[19]に関する動画

※[19]の詳細は[20]. [20]に関する動画

目的として最適なcost-to-goを使用することで, 制御Lyapunov技術との関連性が生まれる.  Ames et al. [15, 21]は, \(l(\mathbf x, \mathbf u,t)\) のnegativityに対する制約を満たしながら, 入力ノルム|| u ||を最小化するQPを解くことによって歩行のために制御Lyapunov関数を使用した.  離散時間設定では, Wang and Boyd [22]は, 状態数が入力数の2乗にほぼ等しい場合に, アクティブ集合の明示的な列挙を使用して制御Lyapunovポリシーを迅速に評価するためのアプローチを説明している.

4.5 Efficient QP Solver

ハードウェア上での実施は, 十分に高い制御速度が達成されることを要求する.  それは短時間でQuadratic Program 1を定式化し解くことができなければならないことを意味する.  重要な観察は, 典型的な動作中に, 不等式制約のアクティブセットが連続する制御ステップ間で非常にまれに変化することである.  効率的なアクティブセットソルバーを使用することでこれを利用できる.  いくつかの汎用QPソルバーに対するタイミング結果を含む詳細については, [23]を参照.
Quadratic Program 1は標準形式で次式のように書ける.
minimize z   \(\frac{1}{2}\mathbf z^T\mathbf{Wz}+\mathbf g^T\mathbf z\)
subject to    \(\mathbf{Az}=\mathbf b\hspace{20mm}(30)\)
\(\hspace{23mm}\mathbf{Mz}\geq\mathbf f\)
ここで, 不等式は\(\mathbf M =(\mathbf m_1, …, \mathbf m_n)^T\), また\(\mathbf f=(f_1, …, f_n)^T\)により定義される. この問題を解決するために, active setと呼ばれる部分集合A⊆{1…n}における各iに対する最適性において, \(\mathbf m^T_i\mathbf z=f_i\)であると仮定される. 時間ステップ\(k> 0\)の場合, このサブセットは時間ステップ\(k-1\)からのactive 不等式のインデックスに等しくなる. この仮定では, QPのKKT条件は\(\mathbf z,\mathbf γ, \mathbf α\)に関して次のように記述できる.
$$-\mathbf g=\mathbf{Wz}+\mathbf A^T\mathbf \alpha+\sum_{i\in A}\gamma_i\mathbf m_i$$

\(\hspace{20mm}\mathbf{Az}=\mathbf b\hspace{30mm}(31)\)
\(\hspace{20mm}\mathbf m^T_i\mathbf z=f_i\;\;\;∀_i\in A\)
\(\hspace{20mm}\gamma_i=0\;\;\;∀_i≠A\)
\(\hspace{20mm}\mathbf{Mz}\leq\mathbf f\)
\(\hspace{20mm}\gamma_i\geq 0\;\;\;∀_i\in A\hspace{30mm}(32)\)
我々の方法は, 線形方程式(31)を解き, そして解\((\mathbf z, \mathbf γ, \mathbf α)\)が不等式(32)を満たすかどうかをチェックする.  不等式が満たされる場合, zはQPを解き, アルゴリズムは終了する. そうでなければ, アルゴリズムは, \(\mathbf m^T_i\mathbf z>f_i\)である場合にインデックスiをAに追加し, \(γ_i<0\)である場合にインデックスiを削除し, (31)を解く.  アルゴリズムは, 不等式(32)が満たされるまで, または指定された最大反復回数に達するまで, このプロセスを繰り返す. この方法は下のAlgorithm 1で概説される.
解が見つからないというまれなケースでは, アルゴリズムはより信頼性の高い(しかし平均的に遅い)内点ソルバーにフェイルオーバーする.  このアルゴリズムでは有限の終了を保証できないため, この対応が必要である. 実際には, QP 1のインスタンスはほとんどの場合1回の反復で解決される. 通常の歩行中にAtlasのQPを解くのに約0.2 msかかります(QPセットアップ時間を含む1 ms)[23]. 足跡の軌跡を評価し, 身体が接触しているかどうかを判断し, ロボットとの間でやり取りするメッセージを処理するコンポーネントなど, すべての追加のコントローラソフトウェアコンポーネントを含め, コントローラは約800 Hzの速度で動作する.


Algorithm 1: 凸QPを解くためのアクティブセット法. 時間ステップ\(k\)でアルゴリズムに渡された集合Aは, 時間\(k – 1\)の最適性でアクティブな制約の集合に等しくなる.

モデル予測制御(MPC)のためのアクティブセット法の他の用途は, MPCにおいて生じるQP間の時間的関係を利用している.  Bartlettらは, MPCのアクティブセット戦略と内点戦略を比較した[24]. 記載されているのは, アクティブセットに変更が行われた後にKKT条件を効率的に解決するためのSchur補完に基づくアクティブセットアプローチである.  Ferreauら[25]は, コスト関数と動的制約が各時間ステップで同じであるMPC問題を考えている.  すなわち, 反復時に解かれるQPは, 初期条件を強制する単一の制約によってのみ異なる.  彼らは初期状態を以前の状態から現在の状態に滑らかに変化させることによって, 最適解によってトレースされた区分的線形パスを追跡することができた. 我々が検討した制御装置は変化するコストと制約構造を持っていたので, この方法は適用するのが困難であっただろう.

4.6 Joint-level control

シミュレーションからハードウェアへの移行には, ロボットの油圧アクチュエータを命令することによってQPからの参照を追跡する, 適切な低レベルの制御が必要です. ロボット上の関節レベルループは1000 Hzで動作し, サーボバルブコマンド, m∈[-10、10] mAを計算する.  下半身制御では, ロボットの動作状況と状態に基づいて選択的に実行する2つのアプローチがある.  1つ目はトルクのみのフィードバックを使用し. 2つ目はトルクと速度のフィードバックを組み合わせたものである.
地形とのコンプライアントな相互作用が望まれる場合は, 2つの足首関節でトルクのみの制御を行う.  たとえば, 地面の傾斜を正確に推定できない場合や, ボードや岩などの小さな障害物がある場合は, 足首の自然なコンプライアンスを利用して足の下の地形に合わせる. トルクのみの制御では, ジョイントで検出されたトルクのフィードバックゲイン(固定の伝達曲線を通してマッピングされたアクチュエータ圧力の測定値を介して計算)とフィードフォワード速度ゲインを組み合わせて, アクチュエータのダイナミクスをキャンセルする.  単純な区分的線形摩擦モデルをデータから当てはめると、適合制御は可能であるが, モデル誤差およびセンサの制限(例えば, 関節トルクセンサの欠如)のために正確に肢を位置決めする能力が制限される.
代わりに, トルクと速度のフィードバックを組み合わせて使用することで, 逆運動学ベースのコントローラーと併用すると非常に優れた位置追跡結果を得ることができるが, ジョイントの自然なコンプライアンスが低下する.  このモードでは, バルブサーボコマンドは次の形式になる.
$$m=K^T_p(\tau_{ref}-\tau)+K^v_p(v_{ref}-v)$$
\(v_{ref}\)を計算するために, QPブロックによって解かれる関節トルクに加えて一般化された加速度を出力し, それらを経時的に積分する. すなわち,
$$v_{ref}(T)=\int_{t_c}^{T}\dot v_{ref}(t)dt$$
ここで, \(t_c\)は, 接触状態が最後に変化した時間である(すなわち, 足が地形と接触して接触しなくなると, 積分速度がリセットされる).

【参考文献】
[1]https://dspace.mit.edu/openaccess-disseminate/1721.1/110533
[2]https://nxr.northwestern.edu/sites/default/files/publications/Sub2015TROAnMu.pdf
[3]http://underactuated.csail.mit.edu/underactuated.html
[4]https://www.researchgate.net/publication/224773123_Dynamic_Programming_and_Optimal_Control
[5]https://ja.wikipedia.org/wiki/%E5%BE%AE%E5%88%86%E5%8B%95%E7%9A%84%E8%A8%88%E7%94%BB%E6%B3%95
[6]http://digital.csic.es/bitstream/10261/133133/1/cost-to-go.pdf
[7]https://pdfs.semanticscholar.org/cf32/008efe00f0748518a2c3321b3b1ff4d2acb8.pdf
[8]https://arxiv.org/pdf/1010.2247.pdf
[9]https://wiki.qut.edu.au/download/attachments/173981546/drake.pdf?version=1&modificationDate=1403592695000&api=v2
[10]http://people.csail.mit.edu/yeuhi/papers/abe-2007-mcf.pdf
[11]https://www.researchgate.net/publication/224401305_Dynamic_balance_control_of_humanoids_for_multiple_grasps_and_non_coplanar_frictional_contacts
[12]http://graphics.cs.ucr.edu/papers/macchietto:2009:mcb.pdf
[13]https://www.researchgate.net/publication/256422371_Summary_of_team_IHMC’s_Virtual_Robotics_Challenge_entry
[14]http://www.cs.cmu.edu/afs/cs/user/sfeng/www/sf_hum13.pdf
[15]http://ames.caltech.edu/ICRA_Ames_2012.pdf
[16]https://www.researchgate.net/publication/236661674_Momentum-based_Balance_Control_for_Torque-controlled_Humanoids
[17]http://projects.laas.fr/gepetto/uploads/Publications/2012-saab-itro.pdf
[18]https://www.researchgate.net/publication/308340634_Team_IHMC’s_Lessons_Learned_from_the_DARPA_Robotics_Challenge_Finding_Data_in_the_Rubble
[19]https://arxiv.org/pdf/1305.2042.pdf
[20]https://arxiv.org/pdf/1410.7284.pdf
[21]http://ames.caltech.edu/hscc42kn.pdf
[22]http://web.stanford.edu/~boyd/papers/pdf/fast_clf.pdf
[23]https://groups.csail.mit.edu/robotics-center/public_papers/Kuindersma13.pdf
[24]http://Active set vs. interior point strategies for model predictive control
[25]http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.491.3821&rep=rep1&type=pdf