Differential dynamic programming

Algorithm for trajectory optimization

Differential dynamic programming (DDP) is an optimal control algorithm of the trajectory optimization class. The algorithm was introduced in 1966 by Mayne[1] and subsequently analysed in Jacobson and Mayne's eponymous book.[2] The algorithm uses locally-quadratic models of the dynamics and cost functions, and displays quadratic convergence. It is closely related to Pantoja's step-wise Newton's method.[3][4]

## Finite-horizon discrete-time problems

[edit]

The dynamics

x i + 1 = f ( x i , u i ) {\displaystyle \mathbf {x} _{i+1}=\mathbf {f} (\mathbf {x} _{i},\mathbf {u} _{i})} |  | 1

|
|

describe the evolution of the state  x {\displaystyle \textstyle \mathbf {x} } given the control  u {\displaystyle \mathbf {u} } from time  i {\displaystyle i} to time  i + 1 {\displaystyle i+1} . The total cost J 0 {\displaystyle J_{0}} is the sum of running costs  ℓ {\displaystyle \textstyle \ell } and final cost  ℓ f {\displaystyle \ell _{f}} , incurred when starting from state  x {\displaystyle \mathbf {x} } and applying the control sequence  U ≡ { u 0 , u 1 … , u N − 1 } {\displaystyle \mathbf {U} \equiv \\{\mathbf {u} _{0},\mathbf {u} _{1}\dots ,\mathbf {u} _{N-1}\\}} until the horizon is reached:

J 0 ( x , U ) = ∑ i = 0 N − 1 ℓ ( x i , u i ) + ℓ f ( x N ) , {\displaystyle J_{0}(\mathbf {x} ,\mathbf {U} )=\sum _{i=0}^{N-1}\ell (\mathbf {x} _{i},\mathbf {u} _{i})+\ell _{f}(\mathbf {x} _{N}),}

where  x 0 ≡ x {\displaystyle \mathbf {x} _{0}\equiv \mathbf {x} } , and the  x i {\displaystyle \mathbf {x} _{i}} for  i > 0 {\displaystyle i>0} are given by Eq. 1. The solution of the optimal control problem is the minimizing control sequence  U ∗ ( x ) ≡ argmin U ⁡ J 0 ( x , U ) . {\displaystyle \mathbf {U} ^{*}(\mathbf {x} )\equiv \operatorname {argmin} _{\mathbf {U} }J_{0}(\mathbf {x} ,\mathbf {U} ).} Trajectory optimization means finding  U ∗ ( x ) {\displaystyle \mathbf {U} ^{*}(\mathbf {x} )} for a particular  x 0 {\displaystyle \mathbf {x} _{0}} , rather than for all possible initial states.

## Dynamic programming

[edit]

Let  U i {\displaystyle \mathbf {U} _{i}} be the partial control sequence  U i ≡ { u i , u i + 1 … , u N − 1 } {\displaystyle \mathbf {U} _{i}\equiv \\{\mathbf {u} _{i},\mathbf {u} _{i+1}\dots ,\mathbf {u} _{N-1}\\}} and define the cost-to-go J i {\displaystyle J_{i}} as the partial sum of costs from  i {\displaystyle i} to  N {\displaystyle N} :

J i ( x , U i ) = ∑ j = i N − 1 ℓ ( x j , u j ) + ℓ f ( x N ) . {\displaystyle J_{i}(\mathbf {x} ,\mathbf {U} _{i})=\sum _{j=i}^{N-1}\ell (\mathbf {x} _{j},\mathbf {u} _{j})+\ell _{f}(\mathbf {x} _{N}).}

The optimal cost-to-go or value function at time  i {\displaystyle i} is the cost-to-go given the minimizing control sequence:

V ( x , i ) ≡ min U i J i ( x , U i ) . {\displaystyle V(\mathbf {x} ,i)\equiv \min _{\mathbf {U} _{i}}J_{i}(\mathbf {x} ,\mathbf {U} _{i}).}

Setting  V ( x , N ) ≡ ℓ f ( x N ) {\displaystyle V(\mathbf {x} ,N)\equiv \ell _{f}(\mathbf {x} _{N})} , the dynamic programming principle reduces the minimization over an entire sequence of controls to a sequence of minimizations over a single control, proceeding backwards in time:

V ( x , i ) = min u [ ℓ ( x , u ) + V ( f ( x , u ) , i + 1 ) ] . {\displaystyle V(\mathbf {x} ,i)=\min _{\mathbf {u} }[\ell (\mathbf {x} ,\mathbf {u} )+V(\mathbf {f} (\mathbf {x} ,\mathbf {u} ),i+1)].} |  | 2

|
|

This is the Bellman equation.

## Differential dynamic programming

[edit]

DDP proceeds by iteratively performing a backward pass on the nominal trajectory to generate a new control sequence, and then a forward-pass to compute and evaluate a new nominal trajectory. We begin with the backward pass. If

ℓ ( x , u ) + V ( f ( x , u ) , i + 1 ) {\displaystyle \ell (\mathbf {x} ,\mathbf {u} )+V(\mathbf {f} (\mathbf {x} ,\mathbf {u} ),i+1)}

is the argument of the  min [ ⋅ ] {\displaystyle \min[\cdot ]} operator in Eq. 2, let  Q {\displaystyle Q} be the variation of this quantity around the  i {\displaystyle i} -th  ( x , u ) {\displaystyle (\mathbf {x} ,\mathbf {u} )} pair:

Q ( δ x , δ u ) ≡ ℓ ( x + δ x , u + δ u ) + V ( f ( x + δ x , u + δ u ) , i + 1 ) − ℓ ( x , u ) − V ( f ( x , u ) , i + 1 ) {\displaystyle {\begin{aligned}Q(\delta \mathbf {x} ,\delta \mathbf {u} )\equiv &\ell (\mathbf {x} +\delta \mathbf {x} ,\mathbf {u} +\delta \mathbf {u} )&&{}+V(\mathbf {f} (\mathbf {x} +\delta \mathbf {x} ,\mathbf {u} +\delta \mathbf {u} ),i+1)\\\\-&\ell (\mathbf {x} ,\mathbf {u} )&&{}-V(\mathbf {f} (\mathbf {x} ,\mathbf {u} ),i+1)\end{aligned}}}

and expand to second order

≈ 1 2 [ 1 δ x δ u ] T [ 0 Q x T Q u T Q x Q x x Q x u Q u Q u x Q u u ] [ 1 δ x δ u ] {\displaystyle \approx {\frac {1}{2}}{\begin{bmatrix}1\\\\\delta \mathbf {x} \\\\\delta \mathbf {u} \end{bmatrix}}^{\mathsf {T}}{\begin{bmatrix}0&Q_{\mathbf {x} }^{\mathsf {T}}&Q_{\mathbf {u} }^{\mathsf {T}}\\\Q_{\mathbf {x} }&Q_{\mathbf {x} \mathbf {x} }&Q_{\mathbf {x} \mathbf {u} }\\\Q_{\mathbf {u} }&Q_{\mathbf {u} \mathbf {x} }&Q_{\mathbf {u} \mathbf {u} }\end{bmatrix}}{\begin{bmatrix}1\\\\\delta \mathbf {x} \\\\\delta \mathbf {u} \end{bmatrix}}} |  | 3

|
|

The  Q {\displaystyle Q} notation used here is a variant of the notation of Morimoto where subscripts denote differentiation in denominator layout.[5] Dropping the index  i {\displaystyle i} for readability, primes denoting the next time-step  V ′ ≡ V ( i + 1 ) {\displaystyle V'\equiv V(i+1)} , the expansion coefficients are

Q x = ℓ x + f x T V x ′ Q u = ℓ u + f u T V x ′ Q x x = ℓ x x + f x T V x x ′ f x + V x ′ ⋅ f x x Q u u = ℓ u u + f u T V x x ′ f u + V x ′ ⋅ f u u Q u x = ℓ u x + f u T V x x ′ f x + V x ′ ⋅ f u x . {\displaystyle {\begin{alignedat}{2}Q_{\mathbf {x} }&=\ell _{\mathbf {x} }+\mathbf {f} _{\mathbf {x} }^{\mathsf {T}}V'_{\mathbf {x} }\\\Q_{\mathbf {u} }&=\ell _{\mathbf {u} }+\mathbf {f} _{\mathbf {u} }^{\mathsf {T}}V'_{\mathbf {x} }\\\Q_{\mathbf {x} \mathbf {x} }&=\ell _{\mathbf {x} \mathbf {x} }+\mathbf {f} _{\mathbf {x} }^{\mathsf {T}}V'_{\mathbf {x} \mathbf {x} }\mathbf {f} _{\mathbf {x} }+V_{\mathbf {x} }'\cdot \mathbf {f} _{\mathbf {x} \mathbf {x} }\\\Q_{\mathbf {u} \mathbf {u} }&=\ell _{\mathbf {u} \mathbf {u} }+\mathbf {f} _{\mathbf {u} }^{\mathsf {T}}V'_{\mathbf {x} \mathbf {x} }\mathbf {f} _{\mathbf {u} }+{V'_{\mathbf {x} }}\cdot \mathbf {f} _{\mathbf {u} \mathbf {u} }\\\Q_{\mathbf {u} \mathbf {x} }&=\ell _{\mathbf {u} \mathbf {x} }+\mathbf {f} _{\mathbf {u} }^{\mathsf {T}}V'_{\mathbf {x} \mathbf {x} }\mathbf {f} _{\mathbf {x} }+{V'_{\mathbf {x} }}\cdot \mathbf {f} _{\mathbf {u} \mathbf {x} }.\end{alignedat}}}

The last terms in the last three equations denote contraction of a vector with a tensor. Minimizing the quadratic approximation (3) with respect to  δ u {\displaystyle \delta \mathbf {u} } we have

δ u ∗ = argmin δ u ⁡ Q ( δ x , δ u ) = − Q u u − 1 ( Q u + Q u x δ x ) , {\displaystyle {\delta \mathbf {u} }^{*}=\operatorname {argmin} \limits _{\delta \mathbf {u} }Q(\delta \mathbf {x} ,\delta \mathbf {u} )=-Q_{\mathbf {u} \mathbf {u} }^{-1}(Q_{\mathbf {u} }+Q_{\mathbf {u} \mathbf {x} }\delta \mathbf {x} ),} |  | 4

|
|

giving an open-loop term  k = − Q u u − 1 Q u {\displaystyle \mathbf {k} =-Q_{\mathbf {u} \mathbf {u} }^{-1}Q_{\mathbf {u} }} and a feedback gain term  K = − Q u u − 1 Q u x {\displaystyle \mathbf {K} =-Q_{\mathbf {u} \mathbf {u} }^{-1}Q_{\mathbf {u} \mathbf {x} }} . Plugging the result back into (3), we now have a quadratic model of the value at time  i {\displaystyle i} :

Δ V ( i ) = − 1 2 Q u T Q u u − 1 Q u V x ( i ) = Q x − Q x u Q u u − 1 Q u V x x ( i ) = Q x x − Q x u Q u u − 1 Q u x . {\displaystyle {\begin{alignedat}{2}\Delta V(i)&=&{}-{\tfrac {1}{2}}Q_{\mathbf {u} }^{T}Q_{\mathbf {u} \mathbf {u} }^{-1}Q_{\mathbf {u} }\\\V_{\mathbf {x} }(i)&=Q_{\mathbf {x} }&{}-Q_{\mathbf {xu} }Q_{\mathbf {u} \mathbf {u} }^{-1}Q_{\mathbf {u} }\\\V_{\mathbf {x} \mathbf {x} }(i)&=Q_{\mathbf {x} \mathbf {x} }&{}-Q_{\mathbf {x} \mathbf {u} }Q_{\mathbf {u} \mathbf {u} }^{-1}Q_{\mathbf {u} \mathbf {x} }.\end{alignedat}}}

Recursively computing the local quadratic models of  V ( i ) {\displaystyle V(i)} and the control modifications  { k ( i ) , K ( i ) } {\displaystyle \\{\mathbf {k} (i),\mathbf {K} (i)\\}} , from  i = N − 1 {\displaystyle i=N-1} down to  i = 1 {\displaystyle i=1} , constitutes the backward pass. As above, the Value is initialized with  V ( x , N ) ≡ ℓ f ( x N ) {\displaystyle V(\mathbf {x} ,N)\equiv \ell _{f}(\mathbf {x} _{N})} . Once the backward pass is completed, a forward pass computes a new trajectory:

x ^ ( 1 ) = x ( 1 ) u ^ ( i ) = u ( i ) + k ( i ) + K ( i ) ( x ^ ( i ) − x ( i ) ) x ^ ( i + 1 ) = f ( x ^ ( i ) , u ^ ( i ) ) {\displaystyle {\begin{aligned}{\hat {\mathbf {x} }}(1)&=\mathbf {x} (1)\\\\{\hat {\mathbf {u} }}(i)&=\mathbf {u} (i)+\mathbf {k} (i)+\mathbf {K} (i)({\hat {\mathbf {x} }}(i)-\mathbf {x} (i))\\\\{\hat {\mathbf {x} }}(i+1)&=\mathbf {f} ({\hat {\mathbf {x} }}(i),{\hat {\mathbf {u} }}(i))\end{aligned}}}

The backward passes and forward passes are iterated until convergence. If the Hessians  Q x x , Q u u , Q u x , Q x u {\displaystyle Q_{\mathbf {x} \mathbf {x} },Q_{\mathbf {u} \mathbf {u} },Q_{\mathbf {u} \mathbf {x} },Q_{\mathbf {x} \mathbf {u} }} are replaced by their Gauss-Newton approximation, the method reduces to the iterative Linear Quadratic Regulator (iLQR).[6]

## Regularization and line-search

[edit]

Differential dynamic programming is a second-order algorithm like Newton's method. It therefore takes large steps toward the minimum and often requires regularization and/or line-search to achieve convergence.[7][8] Regularization in the DDP context means ensuring that the  Q u u {\displaystyle Q_{\mathbf {u} \mathbf {u} }} matrix in Eq. 4 is positive definite. Line-search in DDP amounts to scaling the open-loop control modification  k {\displaystyle \mathbf {k} } by some  0 < α < 1 {\displaystyle 0<\alpha <1} .

## Monte Carlo version

[edit]

Sampled differential dynamic programming (SaDDP) is a Monte Carlo variant of differential dynamic programming.[9][10][11] It is based on treating the quadratic cost of differential dynamic programming as the energy of a Boltzmann distribution. This way the quantities of DDP can be matched to the statistics of a multidimensional normal distribution. The statistics can be recomputed from sampled trajectories without differentiation.

Sampled differential dynamic programming has been extended to Path Integral Policy Improvement with Differential Dynamic Programming.[12] This creates a link between differential dynamic programming and path integral control,[13] which is a framework of stochastic optimal control.

## Constrained problems

[edit]

Interior Point Differential dynamic programming (IPDDP) is an interior-point method generalization of DDP that can address the optimal control problem with nonlinear state and input constraints.[14]

## See also

[edit]

* Optimal control

## References

[edit]

1. ^ Mayne, D. Q. (1966). "A second-order gradient method of optimizing non-linear discrete time systems". Int J Control. 3: 85–95\. doi:10.1080/00207176608921369.
2. ^ Mayne, David Q.; Jacobson, David H. (1970). Differential dynamic programming. New York: American Elsevier Pub. Co. ISBN 978-0-444-00070-5.
3. ^ de O. Pantoja, J. F. A. (1988). "Differential dynamic programming and Newton's method". International Journal of Control. 47 (5): 1539–1553\. doi:10.1080/00207178808906114. ISSN 0020-7179.
4. ^ Liao, L. Z.; C. A Shoemaker (1992). "Advantages of differential dynamic programming over Newton's method for discrete-time optimal control problems". Cornell University. hdl:1813/5474.
5. ^ Morimoto, J.; G. Zeglin; C.G. Atkeson (2003). "Minimax differential dynamic programming: Application to a biped walking robot". Intelligent Robots and Systems, 2003.(IROS 2003). Proceedings. 2003 IEEE/RSJ International Conference on. Vol. 2\. pp. 1927–1932.
6. ^ Baumgärtner, K. (2023). A Unified Local Convergence Analysis of Differential Dynamic Programming, Direct Single Shooting, and Direct Multiple Shooting. 2023 European Control Conference (ECC). pp. 1–7\. doi:10.23919/ECC57647.2023.10178367.
7. ^ Liao, L. Z; C. A Shoemaker (1991). "Convergence in unconstrained discrete-time differential dynamic programming". IEEE Transactions on Automatic Control. 36 (6): 692. doi:10.1109/9.86943.
8. ^ Tassa, Y. (2011). Theory and implementation of bio-mimetic motor controllers (PDF) (Thesis). Hebrew University. Archived from the original (PDF) on 2016-03-04. Retrieved 2012-02-27.
9. ^ "Sampled differential dynamic programming". 2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). doi:10.1109/IROS.2016.7759229. S2CID 1338737.
10. ^ Rajamäki, Joose; Hämäläinen, Perttu (June 2018). Regularizing Sampled Differential Dynamic Programming - IEEE Conference Publication. 2018 Annual American Control Conference (ACC). pp. 2182–2189\. doi:10.23919/ACC.2018.8430799. S2CID 243932441. Retrieved 2018-10-19.
11. ^ Rajamäki, Joose (2018). Random Search Algorithms for Optimal Control. Aalto University. ISBN 978-952-60-8156-4. ISSN 1799-4942.
12. ^ Lefebvre, Tom; Crevecoeur, Guillaume (July 2019). "Path Integral Policy Improvement with Differential Dynamic Programming". 2019 IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM). pp. 739–745\. doi:10.1109/AIM.2019.8868359. hdl:1854/LU-8623968. ISBN 978-1-7281-2493-3. S2CID 204816072.
13. ^ Theodorou, Evangelos; Buchli, Jonas; Schaal, Stefan (May 2010). "Reinforcement learning of motor skills in high dimensions: A path integral approach". 2010 IEEE International Conference on Robotics and Automation. pp. 2397–2403\. doi:10.1109/ROBOT.2010.5509336. ISBN 978-1-4244-5038-1. S2CID 15116370.
14. ^ Pavlov, Andrei; Shames, Iman; Manzie, Chris (2020). "Interior Point Differential Dynamic Programming". arXiv:2004.12710 [math.OC].

## External links

[edit]

* A Python implementation of DDP
* A MATLAB implementation of DDP

* The open-source software framework acados provides an efficient and embeddable implementation of DDP.