Numerical Simulations of a Rolling Ball Robot Actuated by Internal Point Masses

TL;DR

This paper presents a numerical approach for controlling a rolling ball robot actuated by internal point masses, enabling trajectory tracking and obstacle avoidance through continuation methods.

Contribution

It introduces a numerical framework for controlling a ball robot with internal masses along arbitrary rails, advancing motion planning techniques.

Findings

Successful trajectory tracking demonstrated

Effective obstacle avoidance achieved

Numerical solutions validated for complex internal mass configurations

Abstract

The controlled motion of a rolling ball actuated by internal point masses that move along arbitrarily-shaped rails fixed within the ball is considered. The controlled equations of motion are solved numerically using a predictor-corrector continuation method, starting from an initial solution obtained via a direct method, to realize trajectory tracking and obstacle avoidance maneuvers.

Figures

Tables

Table 1. Table 2.1: Initial condition parameter values for the rolling disk. Refer to ( 2.10 ) and ( 2.11 ).

Parameter	Value
${𝜽}_a$
${\dot{𝜽}}_a$
${\varphi }_a$	$0$
$z_a$	$0$
${\dot{z}}_a$	$0$

Table 2. Table 2.2: Final condition parameter values for the rolling disk. Refer to ( 2.11 ).

Parameter	Value
${\dot{𝜽}}_b$
$z_b$	$1$
${\dot{z}}_b$	$0$

Table 3. Table 2.3: Integrand cost function coefficient values for the rolling disk when predictor-corrector continuation is performed in α 𝛼 \alpha . Refer to ( 2.15 ).

Parameter	Value
$\alpha \left(\mu \right)$	$.1+\frac{.95-\mu }{.95-.00001}\left(5000-.1\right)$
${\gamma }_1={\gamma }_2={\gamma }_3={\gamma }_4$	$.1$

Table 4. Table 3.4: Initial condition parameter values for the rolling ball. Refer to ( 3.10 ).

Parameter	Value
${𝜽}_a$
${\dot{𝜽}}_a$
${𝔮}_a$
${𝛀}_a$
${𝒛}_a$

Table 5. Table 3.5: Final condition parameter values for the rolling ball. Refer to ( 3.11 ).

Parameter	Value
${\dot{𝜽}}_b$
${𝛀}_b$
${𝒛}_b$

Table 6. Table 3.6: Integrand cost function coefficient values for the rolling ball when predictor-corrector continuation is performed in the obstacle heights. Refer to ( 3.15 ).

Parameter	Value
${\gamma }_1={\gamma }_2={\gamma }_3$	$10$
$h_1\left(\mu \right)=h_2\left(\mu \right)$	$\frac{.95-\mu }{.95-.00001}\left(1000\right)$
${𝒗}_1$
${𝒗}_2$
${\rho }_1={\rho }_2$	$.282$

Table 7. Table 3.7: Integrand cost function coefficient values for the rolling ball when a second round of predictor-corrector continuation is performed in the control coefficients. Refer to ( 3.15 ).

Parameter	Value
${\gamma }_1\left(\mu \right)={\gamma }_2\left(\mu \right)={\gamma }_3\left(\mu \right)$	$10+\frac{.95-\mu }{.95-.00001}\left(-1000-10\right)$
$h_1=h_2$	$7.846\mathrm{e8}$
${𝒗}_1$
${𝒗}_2$
${\rho }_1={\rho }_2$	$.282$

Table 8. Table B.8: Explanation of shorthand notation for zeroth and first derivatives of 𝐟 ^ ^ 𝐟 \hat{\mathbf{f}} and first and second derivatives of H ^ ^ 𝐻 \hat{H} used in ( B.5 ) and ( B.6 ).

Shorthand	$|$	Extended Shorthand	$|$	Normalized	$|$	Un-Normalized
$\widehat{𝐟}$	=	${\widehat{𝐟}|}_{(s,\widetilde{𝒛}\left(s\right),\mu )}$	=	$\widehat{𝐟}(t\left(s\right),\widetilde{𝒙}\left(s\right),\widetilde{𝝀}\left(s\right),\mu )$	=	$\widehat{𝐟}(t\left(s\right),𝒙\left(t\left(s\right)\right),𝝀\left(t\left(s\right)\right),\mu )$
${\widehat{𝐟}}_{𝝀}$	=	${{\widehat{𝐟}}_{𝝀}|}_{(s,\widetilde{𝒛}\left(s\right),\mu )}$	=	${\widehat{𝐟}}_{𝝀}(t\left(s\right),\widetilde{𝒙}\left(s\right),\widetilde{𝝀}\left(s\right),\mu )$	=	${\widehat{𝐟}}_{𝝀}(t\left(s\right),𝒙\left(t\left(s\right)\right),𝝀\left(t\left(s\right)\right),\mu )$
${\widehat{𝐟}}_{𝒙}$	=	${{\widehat{𝐟}}_{𝒙}|}_{(s,\widetilde{𝒛}\left(s\right),\mu )}$	=	${\widehat{𝐟}}_{𝒙}(t\left(s\right),\widetilde{𝒙}\left(s\right),\widetilde{𝝀}\left(s\right),\mu )$	=	${\widehat{𝐟}}_{𝒙}(t\left(s\right),𝒙\left(t\left(s\right)\right),𝝀\left(t\left(s\right)\right),\mu )$
${\widehat{𝐟}}_t$	=	${{\widehat{𝐟}}_t|}_{(s,\widetilde{𝒛}\left(s\right),\mu )}$	=	${\widehat{𝐟}}_t(t\left(s\right),\widetilde{𝒙}\left(s\right),\widetilde{𝝀}\left(s\right),\mu )$	=	${\widehat{𝐟}}_t(t\left(s\right),𝒙\left(t\left(s\right)\right),𝝀\left(t\left(s\right)\right),\mu )$
${\widehat{𝐟}}_{\mu }$	=	${{\widehat{𝐟}}_{\mu }|}_{(s,\widetilde{𝒛}\left(s\right),\mu )}$	=	${\widehat{𝐟}}_{\mu }(t\left(s\right),\widetilde{𝒙}\left(s\right),\widetilde{𝝀}\left(s\right),\mu )$	=	${\widehat{𝐟}}_{\mu }(t\left(s\right),𝒙\left(t\left(s\right)\right),𝝀\left(t\left(s\right)\right),\mu )$
${\hat{H}}_{𝒙}^{𝖳}$	=	${{\hat{H}}_{𝒙}^{𝖳}|}_{(s,\widetilde{𝒛}\left(s\right),\mu )}$	=	${\hat{H}}_{𝒙}^{𝖳}(t\left(s\right),\widetilde{𝒙}\left(s\right),\widetilde{𝝀}\left(s\right),\mu )$	=	${\hat{H}}_{𝒙}^{𝖳}(t\left(s\right),𝒙\left(t\left(s\right)\right),𝝀\left(t\left(s\right)\right),\mu )$
${\hat{H}}_{𝒙𝒙}$	=	${{\hat{H}}_{𝒙𝒙}|}_{(s,\widetilde{𝒛}\left(s\right),\mu )}$	=	${\hat{H}}_{𝒙𝒙}(t\left(s\right),\widetilde{𝒙}\left(s\right),\widetilde{𝝀}\left(s\right),\mu )$	=	${\hat{H}}_{𝒙𝒙}(t\left(s\right),𝒙\left(t\left(s\right)\right),𝝀\left(t\left(s\right)\right),\mu )$
${\hat{H}}_{𝒙t}$	=	${{\hat{H}}_{𝒙t}|}_{(s,\widetilde{𝒛}\left(s\right),\mu )}$	=	${\hat{H}}_{𝒙t}(t\left(s\right),\widetilde{𝒙}\left(s\right),\widetilde{𝝀}\left(s\right),\mu )$	=	${\hat{H}}_{𝒙t}(t\left(s\right),𝒙\left(t\left(s\right)\right),𝝀\left(t\left(s\right)\right),\mu )$
${\hat{H}}_{𝒙\mu }$	=	${{\hat{H}}_{𝒙\mu }|}_{(s,\widetilde{𝒛}\left(s\right),\mu )}$	=	${\hat{H}}_{𝒙\mu }(t\left(s\right),\widetilde{𝒙}\left(s\right),\widetilde{𝝀}\left(s\right),\mu )$	=	${\hat{H}}_{𝒙\mu }(t\left(s\right),𝒙\left(t\left(s\right)\right),𝝀\left(t\left(s\right)\right),\mu )$

Table 9. Table B.9: Explanation of shorthand notation for 𝐟 ^ ^ 𝐟 \hat{\mathbf{f}} and first derivatives of H ^ ^ 𝐻 \hat{H} evaluated at a 𝑎 a used in ( B.32 ), ( B.33 ), and ( B.34 ). Note that H ^ 𝝀 | a = H ^ 𝝀 ​ ( a , 𝒙 ​ ( a ) , 𝝀 ​ ( a ) , μ ) = 𝐟 ^ 𝖳 | a evaluated-at subscript ^ 𝐻 𝝀 𝑎 subscript ^ 𝐻 𝝀 𝑎 𝒙 𝑎 𝝀 𝑎 𝜇 evaluated-at superscript ^ 𝐟 𝖳 𝑎 \left.\hat{H}_{\boldsymbol{\lambda}}\right|_{a}=\hat{H}_{\boldsymbol{\lambda}}\left(a,{\boldsymbol{x}}(a),\boldsymbol{\lambda}(a),\mu\right)=\left.\hat{\mathbf{f}}^{\mathsf{T}}\right|_{a} .

Shorthand	$|$	Meaning	$|$	Simplification
${{\hat{H}}_{𝒙}|}_a$	=	${\hat{H}}_{𝒙}(a,𝒙\left(a\right),𝝀\left(a\right),\mu )$	=	$H_{𝒙}(a,𝒙\left(a\right),𝝀\left(a\right),𝝅(a,𝒙\left(a\right),𝝀\left(a\right),\mu ),\mu )$
${{\widehat{𝐟}}^{𝖳}|}_a$	=	${\widehat{𝐟}}^{𝖳}(a,𝒙\left(a\right),𝝀\left(a\right),\mu )$	=	${𝐟}^{𝖳}(a,𝒙\left(a\right),𝝀\left(a\right),𝝅(a,𝒙\left(a\right),𝝀\left(a\right),\mu ),\mu )$
${{\hat{H}}_t|}_a$	=	${\hat{H}}_t(a,𝒙\left(a\right),𝝀\left(a\right),\mu )$	=	$H_t(a,𝒙\left(a\right),𝝀\left(a\right),𝝅(a,𝒙\left(a\right),𝝀\left(a\right),\mu ),\mu )$
${{\hat{H}}_{\mu }|}_a$	=	${\hat{H}}_{\mu }(a,𝒙\left(a\right),𝝀\left(a\right),\mu )$	=	$H_{\mu }(a,𝒙\left(a\right),𝝀\left(a\right),𝝅(a,𝒙\left(a\right),𝝀\left(a\right),\mu ),\mu )$

Table 10. Table B.10: Explanation of shorthand notation for first derivatives of 𝝈 𝝈 \boldsymbol{\sigma} used in ( B.32 ), ( B.33 ), and ( B.34 ).

Shorthand	$|$	Meaning
${𝝈}_{𝒙\left(a\right)}$	=	${𝝈}_{𝒙\left(a\right)}(a,𝒙\left(a\right),\mu )$
${𝝈}_a$	=	${𝝈}_a(a,𝒙\left(a\right),\mu )$
${𝝈}_{\mu }$	=	${𝝈}_{\mu }(a,𝒙\left(a\right),\mu )$

Table 11. Table B.11: Explanation of shorthand notation for 𝐟 ^ ^ 𝐟 \hat{\mathbf{f}} and first derivatives of H ^ ^ 𝐻 \hat{H} evaluated at b 𝑏 b used in ( B.32 ), ( B.33 ), and ( B.34 ). Note that H ^ 𝝀 | b = H ^ 𝝀 ​ ( b , 𝒙 ​ ( b ) , 𝝀 ​ ( b ) , μ ) = 𝐟 ^ 𝖳 | b evaluated-at subscript ^ 𝐻 𝝀 𝑏 subscript ^ 𝐻 𝝀 𝑏 𝒙 𝑏 𝝀 𝑏 𝜇 evaluated-at superscript ^ 𝐟 𝖳 𝑏 \left.\hat{H}_{\boldsymbol{\lambda}}\right|_{b}=\hat{H}_{\boldsymbol{\lambda}}\left(b,{\boldsymbol{x}}(b),\boldsymbol{\lambda}(b),\mu\right)=\left.\hat{\mathbf{f}}^{\mathsf{T}}\right|_{b} .

Shorthand	$|$	Meaning	$|$	Simplification
${{\hat{H}}_{𝒙}|}_b$	=	${\hat{H}}_{𝒙}(b,𝒙\left(b\right),𝝀\left(b\right),\mu )$	=	$H_{𝒙}(b,𝒙\left(b\right),𝝀\left(b\right),𝝅(b,𝒙\left(b\right),𝝀\left(b\right),\mu ),\mu )$
${{\widehat{𝐟}}^{𝖳}|}_b$	=	${\widehat{𝐟}}^{𝖳}(b,𝒙\left(b\right),𝝀\left(b\right),\mu )$	=	${𝐟}^{𝖳}(b,𝒙\left(b\right),𝝀\left(b\right),𝝅(b,𝒙\left(b\right),𝝀\left(b\right),\mu ),\mu )$
${{\hat{H}}_t|}_b$	=	${\hat{H}}_t(b,𝒙\left(b\right),𝝀\left(b\right),\mu )$	=	$H_t(b,𝒙\left(b\right),𝝀\left(b\right),𝝅(b,𝒙\left(b\right),𝝀\left(b\right),\mu ),\mu )$
${{\hat{H}}_{\mu }|}_b$	=	${\hat{H}}_{\mu }(b,𝒙\left(b\right),𝝀\left(b\right),\mu )$	=	$H_{\mu }(b,𝒙\left(b\right),𝝀\left(b\right),𝝅(b,𝒙\left(b\right),𝝀\left(b\right),\mu ),\mu )$

Table 12. Table B.12: Explanation of shorthand notation for first derivatives of 𝝍 𝝍 \boldsymbol{\psi} used in ( B.32 ), ( B.33 ), and ( B.34 ).

Shorthand	$|$	Meaning
${𝝍}_{𝒙\left(b\right)}$	=	${𝝍}_{𝒙\left(b\right)}(b,𝒙\left(b\right),\mu )$
${𝝍}_b$	=	${𝝍}_b(b,𝒙\left(b\right),\mu )$
${𝝍}_{\mu }$	=	${𝝍}_{\mu }(b,𝒙\left(b\right),\mu )$

Table 13. Table B.13: Equality between Jacobians of two-point boundary condition functions in normalized and un-normalized coordinates.

Normalized	$|$	Un-Normalized
${\widetilde{𝚼}}_{\widetilde{𝒛}\left(0\right)}(\widetilde{𝒛}\left(0\right),\widetilde{𝒛}\left(1\right),\mu )$	=	${𝚼}_{𝒛\left(a\right)}(𝒛\left(a\right),𝒛\left(b\right),\mu )$
${\widetilde{𝚼}}_{1,\widetilde{𝒛}\left(0\right)}(\widetilde{𝒛}\left(0\right),\widetilde{𝒛}\left(1\right),\mu )$	=	${𝚼}_{1,𝒛\left(a\right)}(𝒛\left(a\right),𝒛\left(b\right),\mu )$
${\widetilde{𝚼}}_{2,\widetilde{𝒛}\left(0\right)}(\widetilde{𝒛}\left(0\right),\widetilde{𝒛}\left(1\right),\mu )$	=	${𝚼}_{2,𝒛\left(a\right)}(𝒛\left(a\right),𝒛\left(b\right),\mu )$
${\widetilde{𝚼}}_{\widetilde{𝒛}\left(1\right)}(\widetilde{𝒛}\left(0\right),\widetilde{𝒛}\left(1\right),\mu )$	=	${𝚼}_{𝒛\left(b\right)}(𝒛\left(a\right),𝒛\left(b\right),\mu )$
${\widetilde{𝚼}}_{1,\widetilde{𝒛}\left(1\right)}(\widetilde{𝒛}\left(0\right),\widetilde{𝒛}\left(1\right),\mu )$	=	${𝚼}_{1,𝒛\left(b\right)}(𝒛\left(a\right),𝒛\left(b\right),\mu )$
${\widetilde{𝚼}}_{2,\widetilde{𝒛}\left(1\right)}(\widetilde{𝒛}\left(0\right),\widetilde{𝒛}\left(1\right),\mu )$	=	${𝚼}_{2,𝒛\left(b\right)}(𝒛\left(a\right),𝒛\left(b\right),\mu )$
${\widetilde{𝚼}}_{\mu }(\widetilde{𝒛}\left(0\right),\widetilde{𝒛}\left(1\right),\mu )$	=	${𝚼}_{\mu }(𝒛\left(a\right),𝒛\left(b\right),\mu )$
${\widetilde{𝚼}}_{1,\mu }(\widetilde{𝒛}\left(0\right),\widetilde{𝒛}\left(1\right),\mu )$	=	${𝚼}_{1,\mu }(𝒛\left(a\right),𝒛\left(b\right),\mu )$
${\widetilde{𝚼}}_{2,\mu }(\widetilde{𝒛}\left(0\right),\widetilde{𝒛}\left(1\right),\mu )$	=	${𝚼}_{2,\mu }(𝒛\left(a\right),𝒛\left(b\right),\mu )$

Equations

$${\mathbf{e}}_1={\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^{\mathsf{T}},{\mathbf{e}}_2={\left[\begin{array}{ccc}0& 1& 0\end{array}\right]}^{\mathsf{T}},\mathrm{and}{\mathbf{e}}_3={\left[\begin{array}{ccc}0& 0& 1\end{array}\right]}^{\mathsf{T}}.$$

$${\mathbf{E}}_1={\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^{\mathsf{T}},{\mathbf{E}}_2={\left[\begin{array}{ccc}0& 1& 0\end{array}\right]}^{\mathsf{T}},\mathrm{and}{\mathbf{E}}_3={\left[\begin{array}{ccc}0& 0& 1\end{array}\right]}^{\mathsf{T}}.$$

$$\begin{split}\dot{\boldsymbol{\Omega}}&=\left[\sum_{i=0}^{n}m_{i}\widehat{\mathbf{s}_{i}}^{2}-\mathbb{I}\right]^{-1}\Bigg{[}\boldsymbol{\Omega}\times\mathbb{I}\boldsymbol{\Omega}+r\tilde{\boldsymbol{\Gamma}}\times\boldsymbol{\Gamma}\\ &\hphantom{=}+\sum_{i=0}^{n}m_{i}\mathbf{s}_{i}\times\left\{g\boldsymbol{\Gamma}+\boldsymbol{\Omega}\times\left(\boldsymbol{\Omega}\times\boldsymbol{\zeta}_{i}+2\dot{\theta}_{i}\boldsymbol{\zeta}_{i}^{\prime}\right)+\dot{\theta}_{i}^{2}\boldsymbol{\zeta}_{i}^{\prime\prime}+\ddot{\theta}_{i}\boldsymbol{\zeta}_{i}^{\prime}\right\}\Bigg{]},\\ \dot{\Lambda}&=\Lambda\widehat{\boldsymbol{\Omega}},\\ \dot{\boldsymbol{z}}&=\left(\Lambda\boldsymbol{\Omega}\times r\mathbf{e}_{3}\right)_{12},\end{split}$$

$${\hat{\mathbf{v}}}^2=\left[\begin{array}{ccc}-(v_2^2+v_3^2)& v_1{v}_2& v_1{v}_3\\ v_1{v}_2& -(v_1^2+v_3^2)& v_2{v}_3\\ v_1{v}_3& v_2{v}_3& -(v_1^2+v_2^2)\end{array}\right]$$

$${\mathbf{v}}_{12}={\left[\begin{array}{cc}{v}_1& v_2\end{array}\right]}^{\mathsf{T}}\in {\mathbb{R}}^2.$$

$$N=Mg+⟨\sum _{i=0}^n{m}_i\left[\dot{\mathbf{\Omega }}\times {\mathbf{s}}_i+\mathbf{\Omega }\times \left(\mathbf{\Omega }\times {\mathbf{\zeta }}_i+2{\dot{\theta }}_i{\mathbf{\zeta }}_i^{\prime }\right)+{\dot{\theta }}_i^2{\mathbf{\zeta }}_i^{\prime \prime }+{\ddot{\theta }}_i{\mathbf{\zeta }}_i^{\prime }\right],\mathbf{\Gamma }⟩-F_{\mathrm{e},3}$$

$$-f_{\mathrm{s}}\mathbf{\sigma }=\left[\begin{array}{c}{\left(\Lambda {\sum }_{i=0}^n{m}_i\left[\dot{\mathbf{\Omega }}\times {\mathbf{s}}_i+\mathbf{\Omega }\times \left(\mathbf{\Omega }\times {\mathbf{\zeta }}_i+2{\dot{\theta }}_i{\mathbf{\zeta }}_i^{\prime }\right)+{\dot{\theta }}_i^2{\mathbf{\zeta }}_i^{\prime \prime }+{\ddot{\theta }}_i{\mathbf{\zeta }}_i^{\prime }\right]-{\mathbf{F}}_{\mathrm{e}}\right)}_{12}\\ 0\end{array}\right].$$

$$\ddot{\varphi }=\frac{-r{F}_{\mathrm{e},1}+{\sum }_{i=0}^n{m}_i{K}_i}{d_2+{\sum }_{i=0}^n{m}_i\left[{\left(r\sin \varphi +{\zeta }_{i,1}\right)}^2+{\left(r\cos \varphi +{\zeta }_{i,3}\right)}^2\right]}\equiv \kappa \left(t,\mathbf{\theta },\dot{\mathbf{\theta }},\varphi ,\dot{\varphi },\ddot{\mathbf{\theta }}\right),$$

$$\begin{array}{rl}{K}_i& \equiv \left(g+r{\dot{\varphi }}^2\right)\left({\zeta }_{i,3}\sin \varphi -{\zeta }_{i,1}\cos \varphi \right)+\left(r\cos \varphi +{\zeta }_{i,3}\right)\left(-2\dot{\varphi }{\dot{\theta }}_i{\zeta }_{i,3}^{\prime }+{\dot{\theta }}_i^2{\zeta }_{i,1}^{\prime \prime }+{\ddot{\theta }}_i{\zeta }_{i,1}^{\prime }\right)\\ & -\left(r\sin \varphi +{\zeta }_{i,1}\right)\left(2\dot{\varphi }{\dot{\theta }}_i{\zeta }_{i,1}^{\prime }+{\dot{\theta }}_i^2{\zeta }_{i,3}^{\prime \prime }+{\ddot{\theta }}_i{\zeta }_{i,3}^{\prime }\right).\end{array}$$

$$z=z_a-r\left(\varphi -{\varphi }_a\right),$$

$$\begin{split}N&=Mg+\sum_{i=0}^{n}m_{i}\Big{[}\left(-\ddot{\phi}\zeta_{i,3}-{\dot{\phi}}^{2}\zeta_{i,1}-2\dot{\phi}{\dot{\theta}}_{i}\zeta_{i,3}^{\prime}+{\dot{\theta}}_{i}^{2}\zeta_{i,1}^{\prime\prime}+{\ddot{\theta}}_{i}\zeta_{i,1}^{\prime}\right)\sin\phi\\ &\hphantom{=Mg+\sum_{i=0}^{n}m_{i}\Big{[}}+\left(\ddot{\phi}\zeta_{i,1}-{\dot{\phi}}^{2}\zeta_{i,3}+2\dot{\phi}{\dot{\theta}}_{i}\zeta_{i,1}^{\prime}+{\dot{\theta}}_{i}^{2}\zeta_{i,3}^{\prime\prime}+{\ddot{\theta}}_{i}\zeta_{i,3}^{\prime}\right)\cos\phi\Big{]}-F_{\mathrm{e},3}\end{split}$$

$$\begin{split}-f_{\mathrm{s}}\boldsymbol{\sigma}&=-\Big{\{}Mr\ddot{\phi}+\sum_{i=0}^{n}m_{i}\Big{[}\left(\ddot{\phi}\zeta_{i,3}+{\dot{\phi}}^{2}\zeta_{i,1}+2\dot{\phi}{\dot{\theta}}_{i}\zeta_{i,3}^{\prime}-{\dot{\theta}}_{i}^{2}\zeta_{i,1}^{\prime\prime}-{\ddot{\theta}}_{i}\zeta_{i,1}^{\prime}\right)\cos\phi\\ &\hphantom{=-\Big{\{}Mr\ddot{\phi}+\sum_{i=0}^{n}m_{i}\Big{[}}+\left(\ddot{\phi}\zeta_{i,1}-{\dot{\phi}}^{2}\zeta_{i,3}+2\dot{\phi}{\dot{\theta}}_{i}\zeta_{i,1}^{\prime}+{\dot{\theta}}_{i}^{2}\zeta_{i,3}^{\prime\prime}+{\ddot{\theta}}_{i}\zeta_{i,3}^{\prime}\right)\sin\phi\Big{]}+F_{\mathrm{e},1}\Big{\}}\mathbf{e}_{1}.\end{split}$$

$${\mathbf{\zeta }}_i\left({\theta }_i\right)=r_i\left[\begin{array}{c}\cos {\theta }_i\\ 0\\ \sin {\theta }_i\end{array}\right].$$

$$z_{\mathrm{d}}(t)\equiv \left[z_{a}w(t)+{\tilde{z}}_{\mathrm{d}}(t)\left(1-w(t)\right)\right]\left(1-y(t)\right)+z_{b}y(t),$$

$$S(t)\equiv \frac{1}{2}\left[1+\tanh \left(\frac{-t}{ϵ}\right)\right],$$

$$w(t)\equiv S\left(t-a\right),$$

$$y(t)\equiv S\left(-t+b\right),$$

$${\tilde{z}}_{\mathrm{d}}(t)\equiv \left[z_a+\left(z_b-z_a\right)\frac{t-a}{b-a}\right]\sin \left(\frac{9\pi }{2}\frac{t-a}{b-a}\right),$$

$$\min_{\boldsymbol{u}}J\mbox{\, s.t. \,}\left\{\begin{array}[]{ll}\dot{{\boldsymbol{x}}}=\mathbf{f}\left({\boldsymbol{x}},\boldsymbol{u}\right),\\ \boldsymbol{\sigma}\left({\boldsymbol{x}}(a)\right)=\mathbf{0},\\ \boldsymbol{\psi}\left({\boldsymbol{x}}(b)\right)=\mathbf{0}.\end{array}\right.$$

$$\mathbf{x}\equiv \left[\begin{array}{c}\mathbf{\theta }\\ \dot{\mathbf{\theta }}\\ \varphi \\ \dot{\varphi }\end{array}\right]\mathrm{and}\mathbf{u}\equiv \ddot{\mathbf{\theta }},$$

$$\dot{\mathbf{x}}=\left[\begin{array}{c}\dot{\mathbf{\theta }}\\ \ddot{\mathbf{\theta }}\\ \dot{\varphi }\\ \ddot{\varphi }\end{array}\right]=\mathbf{f}\left(\mathbf{x},\mathbf{u}\right)\equiv \left[\begin{array}{c}\dot{\mathbf{\theta }}\\ \mathbf{u}\\ \dot{\varphi }\\ \kappa \left(\mathbf{x},\mathbf{u}\right)\end{array}\right],$$

$$\mathbf{\sigma }\left(\mathbf{x}(a)\right)\equiv \left[\begin{array}{c}\mathbf{\theta }(a)-{\mathbf{\theta }}_a\\ \dot{\mathbf{\theta }}(a)-{\dot{\mathbf{\theta }}}_a\\ \varphi (a)-{\varphi }_a\\ -r\dot{\varphi }(a)-{\dot{z}}_a\end{array}\right]=0,$$

$$\mathbf{\psi }\left(\mathbf{x}(b)\right)\equiv \left[\begin{array}{c}\Pi \left(\tilde{\Lambda }\left(\varphi (b)\right)\left[\frac{1}{M}{\sum }_{i=0}^4{m}_i{\mathbf{\zeta }}_i\left({\theta }_i(b)\right)\right]\right)\\ \dot{\mathbf{\theta }}(b)-{\dot{\mathbf{\theta }}}_b\\ z_a-r\left(\varphi (b)-{\varphi }_a\right)-z_b\\ -r\dot{\varphi }(b)-{\dot{z}}_b\end{array}\right]=0.$$

$$\frac{1}{M}\sum _{i=0}^4{m}_i{\mathbf{\zeta }}_i\left({\theta }_i(t)\right)$$

$$\tilde{\Lambda }(\varphi (t))=\Lambda (t)=\left[\begin{array}{ccc}\cos \varphi (t)& 0& -\sin \varphi (t)\\ 0& 1& 0\\ \sin \varphi (t)& 0& \cos \varphi (t)\end{array}\right]$$

$$J\equiv {\int }_a^{b}L\left(t,\mathbf{x},\mathbf{u},\mu \right)\mathrm{d}t={\int }_a^b\left[\frac{\alpha (\mu )}{2}{\left(z_a-r\left(\varphi -{\varphi }_a\right)-z_{\mathrm{d}}\right)}^2+\sum _{i=1}^4\frac{{\gamma }_i}{2}{\ddot{\theta }}_i^2\right]\mathrm{d}t,$$

$$L\left(t,\mathbf{x},\mathbf{u},\mu \right)\equiv \frac{\alpha (\mu )}{2}{\left(z_a-r\left(\varphi -{\varphi }_a\right)-z_{\mathrm{d}}\right)}^2+\sum _{i=1}^4\frac{{\gamma }_i}{2}{\ddot{\theta }}_i^2,$$

$$\begin{array}{rl}\dot{\mathbf{x}}& ={\hat{H}}_{\mathbf{\lambda }}^{\mathsf{T}}\left(t,\mathbf{x},\mathbf{\lambda },\mu \right)=\hat{\mathbf{f}}\left(\mathbf{x},\mathbf{\lambda }\right)\equiv \mathbf{f}\left(\mathbf{x},\mathbf{\pi }\left(\mathbf{x},\mathbf{\lambda }\right)\right),\\ \dot{\mathbf{\lambda }}& =-{\hat{H}}_{\mathbf{x}}^{\mathsf{T}}\left(t,\mathbf{x},\mathbf{\lambda },\mu \right)=-H_{\mathbf{x}}^{\mathsf{T}}\left(t,\mathbf{x},\mathbf{\lambda },\mathbf{\pi }\left(\mathbf{x},\mathbf{\lambda }\right),\mu \right),\\ {\mathbf{\lambda }|}_{t=a}& =-G_{\mathbf{x}(a)}^{\mathsf{T}},G_{\mathbf{\xi }}^{\mathsf{T}}=\mathbf{\sigma }\left(\mathbf{x}(a)\right)=0,\\ {\mathbf{\lambda }|}_{t=b}& =G_{\mathbf{x}(b)}^{\mathsf{T}},G_{\mathbf{\nu }}^{\mathsf{T}}=\mathbf{\psi }\left(\mathbf{x}(b)\right)=0.\end{array}$$

$$\begin{array}{rl}G\left(\mathbf{x}(a),\mathbf{\xi },\mathbf{x}(b),\mathbf{\nu }\right)& \equiv {\mathbf{\xi }}^{\mathsf{T}}\mathbf{\sigma }\left(\mathbf{x}(a)\right)+{\mathbf{\nu }}^{\mathsf{T}}\mathbf{\psi }\left(\mathbf{x}(b)\right)\\ & ={\mathbf{\xi }}^{\mathsf{T}}\left[\begin{array}{c}\mathbf{\theta }(a)-{\mathbf{\theta }}_a\\ \dot{\mathbf{\theta }}(a)-{\dot{\mathbf{\theta }}}_a\\ \varphi (a)-{\varphi }_a\\ -r\dot{\varphi }(a)-{\dot{z}}_a\end{array}\right]+{\mathbf{\nu }}^{\mathsf{T}}\left[\begin{array}{c}\Pi \left(\tilde{\Lambda }\left(\varphi (b)\right)\left[\frac{1}{M}{\sum }_{i=0}^4{m}_i{\mathbf{\zeta }}_i\left({\theta }_i(b)\right)\right]\right)\\ \dot{\mathbf{\theta }}(b)-{\dot{\mathbf{\theta }}}_b\\ z_a-r\left(\varphi (b)-{\varphi }_a\right)-z_b\\ -r\dot{\varphi }(b)-{\dot{z}}_b\end{array}\right],\end{array}$$

$$\begin{array}{rl}H\left(t,\mathbf{x},\mathbf{\lambda },\mathbf{u},\mu \right)& \equiv L\left(t,\mathbf{x},\mathbf{u},\mu \right)+{\mathbf{\lambda }}^{\mathsf{T}}\mathbf{f}\left(\mathbf{x},\mathbf{u}\right)\\ & =\frac{\alpha (\mu )}{2}{\left(z_a-r\left(\varphi -{\varphi }_a\right)-z_{\mathrm{d}}\right)}^2+\sum _{i=1}^4\frac{{\gamma }_i}{2}{\ddot{\theta }}_i^2+{\mathbf{\lambda }}^{\mathsf{T}}\left[\begin{array}{c}\dot{\mathbf{\theta }}\\ \mathbf{u}\\ \dot{\varphi }\\ \kappa \left(\mathbf{x},\mathbf{u}\right)\end{array}\right],\end{array}$$

Full text

Keywords: optimal control, rolling ball robots, trajectory tracking, obstacle avoidance, predictor-corrector continuation

1. 1 Introduction

2. 1.1 Overview

3. 1.2 Rolling Ball

4. 1.3 Rolling Disk

5. 1.4 Numerical Methods

6. 2 Trajectory Tracking for the Rolling Disk

7. 2.1 Optimal Control Problem and Controlled Equations of Motion

8. 2.2 Numerical Solutions: Trajectory Tracking

9. 2.2.1 Turning Points

10. 3 Obstacle Avoidance for the Rolling Ball

11. 3.1 Optimal Control Problem and Controlled Equations of Motion

12. 3.2 Numerical Solutions: Sigmoid Obstacle Avoidance

13. 3.3 Numerical Solutions: $\operatorname{ReLU}$ Obstacle Avoidance

14. 4 Summary, Discussion, and Future Work

15. Acknowledgements

16. A Optimal Control: Variational Pontryagin’s Minimum Principle

17. B Implementation Details for Solving the ODE TPBVP for a Regular Optimal Control Problem

18. B.1 Normalization and ODE Velocity Function

19. B.2 Two-Point Boundary Condition Function

20. B.3 Final Details

21. C Predictor-Corrector Continuation Method for Solving an ODE TPBVP

22. C.1 Introduction

23. C.2 A Hilbert Space

24. C.3 The Fréchet Derivative and Newton’s Method

25. C.4 The Davidenko ODE IVP

26. C.5 Construct the Tangent

27. C.6 Normalize the Tangent

28. C.7 Construct the Tangent Predictor

29. C.8 Construct the Corrector

30. C.9 Polish the Corrector

31. C.10 Pseudocode for Predictor-Corrector Continuation

32. D Sweep Predictor-Corrector Continuation Method for Solving an ODE TPBVP

33. D.1 Introduction

34. D.2 Construct the Tangent

35. D.3 Determine the Tangent Direction

36. D.4 Sweep along the Tangent

37. D.5 Pseudocode for Sweep Predictor-Corrector Continuation

1 Introduction

@fb@secFB

1.1 Overview

This paper is a continuation of [1], providing numerical solutions of the controlled equations of motion for several special cases of the rolling disk and ball actuated by moving internal point masses. The paper [1] invokes Pontryagin’s minimum principle to derive the theoretical background for the optimal control of the rolling disk and ball having general performance indexes. This paper implements the theory derived in [1] to solve several practical examples, such as trajectory tracking for the rolling disk and obstacle avoidance for the rolling ball. The key contributions of this paper are listed below.

• •

The controlled equations of motion, for a rolling ball actuated by internal point masses that move along arbitrarily-shaped rails fixed within the ball, are solved numerically by a predictor-corrector continuation method, starting from an initial solution provided by a direct method.

• •

Jacobians of the ordinary differential equations (ODEs) and boundary conditions (BCs) which constitute the controlled equations of motion (i.e. an ordinary differential equation two-point boundary value problem (ODE TPBVP)) corresponding to the optimal control of a dynamical system are derived. These Jacobians are useful for numerically solving the controlled equations of motion.

• •

Algorithms for solving an ODE TPBVP by predictor-corrector continuation are developed and were implemented in LAB to numerically solve the controlled equations of motion for the rolling ball. There are not very many predictor-corrector continuation methods publicly available for solving dynamical systems. The idea of using a monotonic continuation ODE TPBVP solver in conjunction with a predictor-corrector continuation method to advance (or “sweep”) as far along the tangent as possible is new, and this novel technique was used to obtain all the numerical results in this paper.

The paper is organized as follows. Subsections 1.2 and 1.3 review the specific types of rolling disk and ball considered, define coordinate systems and notation used to describe this rolling disk and ball, and present the uncontrolled equations of motion for this rolling disk and ball derived earlier in [2, 3]. In Sections 2 and 3, the controlled equations of motion for the rolling disk and ball are formulated and solved numerically via a predictor-corrector continuation method, starting from an initial solution provided by a direct method. Subsection 1.4 provides details of the numerical methods used to solve the controlled equations of motion. Section 4 summarizes the results of the paper and discusses topics for future work. The background material for this paper is contained in several appendices. In particular, Appendix A reviews the theory of optimal control needed to derive the controlled equations of motion for a generic dynamical system given initial and final conditions, given a performance index to be minimized, and in the absence of path inequality constraints. Appendix B provides details for numerically solving the controlled equations of motion. Appendices C and D develop two predictor-corrector continuation algorithms which numerically solve an ODE TPBVP, the latter of which is utilized to numerically solve the controlled equations of motion for the rolling disk and ball.

@fb@secFB

1.2 Rolling Ball

Consider a rigid ball of radius $r$ containing some static internal structure as well as $n$ point masses, where either $n$ is a positive integer denoting the number of moving masses or $n=0$ if no moving masses are used and the structure of the ball is static. This ball rolls without slipping on a horizontal surface in the presence of a uniform gravitational field of magnitude $g$, as illustrated in Figure 1.1. The ball with its static internal structure has mass $m_{0}$ and the $i^{\mathrm{th}}$ point mass has mass $m_{i}$ for $1\leq i\leq n$. Let $M=\sum_{i=0}^{n}m_{i}$ denote the mass of the total system. The total mechanical system consisting of the ball with its static internal structure and the $n$ point masses is referred to as the ball or the rolling ball, the ball with its static internal structure but without the $n$ point masses may also be referred to as $m_{0}$, and the $i^{\mathrm{th}}$ point mass may also be referred to as $m_{i}$ for $1\leq i\leq n$. Note that this system is the Chaplygin ball [4] equipped with point masses.

Two coordinate systems, or frames of reference, will be used to describe the motion of the rolling ball, an inertial spatial coordinate system and a body coordinate system in which each particle within the ball is always fixed. For brevity, the spatial coordinate system will be referred to as the spatial frame and the body coordinate system will be referred to as the body frame. These two frames are depicted in Figure 1.1. The spatial frame has orthonormal axes $\mathbf{e}_{1}$, $\mathbf{e}_{2}$, $\mathbf{e}_{3}$, such that the $\mathbf{e}_{1}$-$\mathbf{e}_{2}$ plane is parallel to the horizontal surface and passes through the ball’s geometric center (i.e. the $\mathbf{e}_{1}$-$\mathbf{e}_{2}$ plane is a height $r$ above the horizontal surface), such that $\mathbf{e}_{3}$ is vertical (i.e. $\mathbf{e}_{3}$ is perpendicular to the horizontal surface) and points “upward” and away from the horizontal surface, and such that $\left(\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{e}_{3}\right)$ forms a right-handed coordinate system. For simplicity, the spatial frame axes are chosen to be

[TABLE]

The acceleration due to gravity in the uniform gravitational field is $\mathfrak{g}=-g\mathbf{e}_{3}=\begin{bmatrix}0&0&-g\end{bmatrix}^{\mathsf{T}}$ in the spatial frame.

The body frame’s origin is chosen to coincide with the position of $m_{0}$’s center of mass. The body frame has orthonormal axes $\mathbf{E}_{1}$, $\mathbf{E}_{2}$, and $\mathbf{E}_{3}$, chosen to coincide with $m_{0}$’s principal axes, in which $m_{0}$’s inertia tensor $\mathbb{I}$ is diagonal, with corresponding principal moments of inertia $d_{1}$, $d_{2}$, and $d_{3}$. That is, in this body frame the inertia tensor is the diagonal matrix $\mathbb{I}=\mathrm{\textbf{diag}}\left(\begin{bmatrix}d_{1}&d_{2}&d_{3}\end{bmatrix}\right)$. Moreover, $\mathbf{E}_{1}$, $\mathbf{E}_{2}$, and $\mathbf{E}_{3}$ are chosen so that $\left(\mathbf{E}_{1},\mathbf{E}_{2},\mathbf{E}_{3}\right)$ forms a right-handed coordinate system. For simplicity, the body frame axes are chosen to be

[TABLE]

In the spatial frame, the body frame is the moving frame $\left(\Lambda\left(t\right)\mathbf{E}_{1},\Lambda\left(t\right)\mathbf{E}_{2},\Lambda\left(t\right)\mathbf{E}_{3}\right)$, where $\Lambda\left(t\right)\in SO(3)$ defines the orientation (or attitude) of the ball at time $t$ relative to its reference configuration, for example at some initial time.

For $1\leq i\leq n$, it is assumed that $m_{i}$ moves along its own 1-d rail. It is further assumed that the $i^{\mathrm{th}}$ rail is parameterized by a 1-d parameter $\theta_{i}$, so that the trajectory $\boldsymbol{\zeta}_{i}$ of the $i^{\mathrm{th}}$ rail, in the body frame translated to the ball’s geometric center, as a function of $\theta_{i}$ is $\boldsymbol{\zeta}_{i}(\theta_{i})$. Refer to Figure 1.1 for an illustration. Therefore, the body frame vector from the ball’s geometric center to $m_{i}$’s center of mass is denoted by $\boldsymbol{\zeta}_{i}(\theta_{i}(t))$. Since $m_{0}$ is stationary in the body frame and to be consistent with the positional notation for $m_{i}$ for $1\leq i\leq n$, $\boldsymbol{\zeta}_{0}\equiv\boldsymbol{\zeta}_{0}(\theta_{0})\equiv\boldsymbol{\zeta}_{0}(\theta_{0}(t))$ is the constant (time-independent) vector from the ball’s geometric center to $m_{0}$’s center of mass for any scalar-valued, time-varying function $\theta_{0}(t)$. In addition, suppose a time-varying external force $\mathbf{F}_{\mathrm{e}}(t)$ acts at the ball’s geometric center.

Let $\mathbf{z}_{i}(t)$ denote the position of $m_{i}$’s center of mass in the spatial frame so that the position of $m_{i}$’s center of mass in the spatial frame is $\mathbf{z}_{i}(t)=\mathbf{z}_{0}(t)+\Lambda(t)\left[\boldsymbol{\zeta}_{i}(\theta(t))-\boldsymbol{\zeta}_{0}\right]$. In general, a particle with position $\mathbf{w}(t)$ in the body frame has position $\mathbf{z}(t)=\mathbf{z}_{0}(t)+\Lambda(t)\mathbf{w}(t)$ in the spatial frame and has position $\mathbf{w}(t)+\boldsymbol{\zeta}_{0}$ in the body frame translated to the ball’s geometric center.

For conciseness, the ball’s geometric center is often denoted GC, $m_{0}$’s center of mass is often denoted CM, and the ball’s contact point with the surface is often denoted CP. The GC is located at $\mathbf{z}_{\mathrm{GC}}(t)=\mathbf{z}_{0}(t)-\Lambda(t)\boldsymbol{\zeta}_{0}$ in the spatial frame, at $-\boldsymbol{\zeta}_{0}$ in the body frame, and at $\mathbf{0}=\begin{bmatrix}0&0&0\end{bmatrix}^{\mathsf{T}}$ in the body frame translated to the GC. The CM is located at $\mathbf{z}_{0}(t)$ in the spatial frame, at $\mathbf{0}$ in the body frame, and at $\boldsymbol{\zeta}_{0}$ in the body frame translated to the GC. The CP is located at $\mathbf{z}_{\mathrm{CP}}(t)=\mathbf{z}_{0}(t)-\Lambda(t)\left[r\boldsymbol{\Gamma}(t)+\boldsymbol{\zeta}_{0}\right]$ in the spatial frame, at $-\left[r\boldsymbol{\Gamma}(t)+\boldsymbol{\zeta}_{0}\right]$ in the body frame, and at $-r\boldsymbol{\Gamma}(t)$ in the body frame translated to the GC, where $\boldsymbol{\Gamma}(t)\equiv\Lambda^{-1}(t)\mathbf{e}_{3}$. Since the third spatial coordinate of the ball’s GC is always [math] and of the ball’s CP is always $-r$, only the first two spatial coordinates of the ball’s GC and CP, denoted by $\boldsymbol{z}(t)$, are needed to determine the spatial location of the ball’s GC and CP.

For succintness, the explicit time dependence of variables is often dropped. That is, the orientation of the ball at time $t$ is denoted simply $\Lambda$ rather than $\Lambda(t)$, the position of $m_{i}$’s center of mass in the spatial frame at time $t$ is denoted $\mathbf{z}_{i}$ rather than $\mathbf{z}_{i}(t)$, the position of $m_{i}$’s center of mass in the body frame translated to the GC at time $t$ is denoted $\boldsymbol{\zeta}_{i}$ or $\boldsymbol{\zeta}_{i}(\theta_{i})$ rather than $\boldsymbol{\zeta}_{i}(\theta_{i}(t))$, the spatial $\mathbf{e}_{1}$-$\mathbf{e}_{2}$ position of the ball’s GC at time $t$ is denoted $\boldsymbol{z}$ rather than $\boldsymbol{z}(t)$, and the external force is denoted $\mathbf{F}_{\mathrm{e}}$ rather than $\mathbf{F}_{\mathrm{e}}(t)$.

As shown in [2, 3], the uncontrolled equations of motion for this rolling ball are

[TABLE]

where $\mathbf{s}_{i}\equiv r\boldsymbol{\Gamma}+\boldsymbol{\zeta}_{i}$ is the body frame vector from the CP to $m_{i}$ for $0\leq i\leq n$, $\boldsymbol{\Omega}\equiv\left(\Lambda^{-1}\dot{\Lambda}\right)^{\vee}$ is the ball’s body angular velocity, $\boldsymbol{\Gamma}\equiv\Lambda^{-1}\mathbf{e}_{3}$ is the spatial unit normal expressed in the body frame, and $\tilde{\boldsymbol{\Gamma}}\equiv\Lambda^{-1}\mathbf{F}_{\mathrm{e}}$ is the external force expressed in the body frame. For $\mathbf{v}=\begin{bmatrix}v_{1}&v_{2}&v_{3}\end{bmatrix}^{\mathsf{T}}\in\mathbb{R}^{3}$, $\widehat{\mathbf{v}}^{2}=\widehat{\mathbf{v}}\widehat{\mathbf{v}}$ is the symmetric matrix given by

[TABLE]

and $\mathbf{v}_{12}$ is the projected vector consisting of the first two components of $\mathbf{v}$ so that

[TABLE]

Let $N$ denote the magnitude of the normal force acting at the ball’s CP. Let $f_{\mathrm{s}}$ and $\boldsymbol{\sigma}$ denote the magnitude of and unit-length direction antiparallel to the static friction acting at the ball’s CP, respectively. As shown in [3], the magnitude of the normal force is

[TABLE]

and the static friction is

[TABLE]

The dynamics encapsulated by (1.3) are valid only if the ball does not detach from the surface ($N>0$) and rolls without slipping ($\mu_{\mathrm{s}}N\geq f_{\mathrm{s}}$), where $\mu_{\mathrm{s}}$ denotes the coefficient of static friction between the ball and the surface.

@fb@secFB

1.3 Rolling Disk

Now suppose that the ball’s inertia is such that one of the ball’s principal axes, say the one labeled $\mathbf{E}_{2}$, is orthogonal to the plane containing the GC and CM. Also assume that all the point masses move along 1-d rails which lie in the plane containing the GC and CM. Moreover, suppose that the ball is oriented initially so that the plane containing the GC and CM coincides with the $\mathbf{e}_{1}$-$\mathbf{e}_{3}$ plane and that the external force $\mathbf{F}_{\mathrm{e}}$ acts in the $\mathbf{e}_{1}$-$\mathbf{e}_{3}$ plane. Then for all time, the ball will remain oriented so that the plane containing the GC and CM coincides with the $\mathbf{e}_{1}$-$\mathbf{e}_{3}$ plane and the ball will only move in the $\mathbf{e}_{1}$-$\mathbf{e}_{3}$ plane, with the ball’s rotation axis always parallel to $\mathbf{e}_{2}$. Note that the dynamics of this system are equivalent to that of the Chaplygin disk [4], equipped with point masses, rolling in the $\mathbf{e}_{1}$-$\mathbf{e}_{3}$ plane, and where the Chaplygin disk (minus the point masses) has polar moment of inertia $d_{2}$. Therefore, henceforth, this particular ball with this special inertia, orientation, and placement of the rails and point masses, may be referred to as the disk or the rolling disk. Figure 1.2 depicts the rolling disk. Let $\phi$ denote the angle between $\mathbf{e}_{1}$ and $\mathbf{E}_{1}$, measured counterclockwise from $\mathbf{e}_{1}$ to $\mathbf{E}_{1}$. Thus, if $\dot{\phi}>0$, the disk rolls in the $-\mathbf{e}_{1}$ direction and $\boldsymbol{\Omega}$ has the same direction as $-\mathbf{e}_{2}$, and if $\dot{\phi}<0$, the disk rolls in the $\mathbf{e}_{1}$ direction and $\boldsymbol{\Omega}$ has the same direction as $\mathbf{e}_{2}$.

As shown in [2, 3], the uncontrolled equation of motion for this rolling disk is

[TABLE]

where

[TABLE]

In (1.8), $\kappa$ is a function that depends on time ($t$) through the possibly time-varying external force $F_{\mathrm{e},1}(t)$, on the point mass parameterized positions ($\boldsymbol{\theta}$), velocities ($\dot{\boldsymbol{\theta}}$), and accelerations ($\ddot{\boldsymbol{\theta}}$), and on the disk’s orientation angle ($\phi$) and its time-derivative ($\dot{\phi}$). The spatial $\mathbf{e}_{1}$ position $z$ of the disk’s GC is given by

[TABLE]

where $z_{a}$ is the spatial $\mathbf{e}_{1}$ position of the disk’s GC at time $t=a$ and $\phi_{a}$ is the disk’s angle at time $t=a$. Let $N$ denote the magnitude of the normal force acting at the disk’s CP. Let $f_{\mathrm{s}}$ and $\boldsymbol{\sigma}$ denote the magnitude of and unit-length direction antiparallel to the static friction acting at the disk’s CP, respectively. As shown in [3], the magnitude of the normal force is

[TABLE]

and the static friction is

[TABLE]

The dynamics encapsulated by (1.8) are valid only if the disk does not detach from the surface ($N>0$) and rolls without slipping ($\mu_{\mathrm{s}}N\geq f_{\mathrm{s}}$), where $\mu_{\mathrm{s}}$ denotes the coefficient of static friction between the disk and the surface.

@fb@secFB

1.4 Numerical Methods

In Sections 2 and 3, the motions of the rolling disk and ball are simulated in LAB R2019b by numerically solving the controlled equations of motion (2.16) and (3.23) corresponding to the optimal control problems (2.7) and (3.4) for the rolling disk and ball, respectively. Subsection 2.2 simulates the rolling disk, while Subsections 3.2 and 3.3 simulate the rolling ball. Because the controlled equations of motion have a very small radius of convergence [5, 6, 7], a direct method, namely the LAB toolbox PS-II [8] version 2.5, is first used to construct a good initial guess. In these simulations, PS-II is configured to use the NLP solver SNOPT [9, 10] version 7.6.0, though PS-II can also be configured to use the NLP solver IPOPT [11]. For the rolling disk, the direct method is used to solve the rolling disk optimal control problem (2.7). When using the direct method to solve the rolling ball optimal control problem, the differential-algebraic equation (DAE) formulation (3.19) is solved first. The direct method solution to the DAE formulation is then used as an initial guess to solve the ODE formulation (3.4), which is consistent with the controlled equations of motion (3.23) for the rolling ball, by the direct method. The LAB automatic differentiation toolbox Gator [12, 13] version 1.5 is used to supply vectorized first derivatives (i.e. Jacobians) to the direct method solver PS-II, since SNOPT accepts first, but not second, derivatives.

Starting from the initial guess provided by the direct method, the controlled equations of motion (2.16) and (3.23) are solved by predictor-corrector continuation in the parameter $\mu$, utilizing the algorithm described in Appendix D. The predictor-corrector continuation method uses the LAB global method ODE TPBVP solvers p [14] version 1.0 or twp [15] version 1.0. By vectorized automatic differentiation of $H_{{\boldsymbol{x}}}$, $\boldsymbol{\pi}$, and $\hat{\mathbf{f}}$, Gator is used to numerically construct the Jacobians of the normalized ODE velocity function (B.5) and (B.6). By non-vectorized automatic differentiation of the Hamiltonian $H$, the initial condition function $\boldsymbol{\sigma}$, the final condition function $\boldsymbol{\psi}$, and the endpoint function $G$, Gator is used to numerically construct the normalized BC function (B.40) and the Jacobians of the normalized BC function (B.42), (B.43), and (B.44). These functions are needed by the ODE TPBVP solvers p and twp to solve the controlled equations of motion (2.16) and (3.23) by predictor-corrector continuation in the parameter $\mu$.

In contrast to the direct method, the controlled equations of motion obtained via the indirect method have a very small radius of convergence [5, 6, 7]. Therefore, the direct method is needed to initialize the predictor-corrector continuation of the controlled equations of motion. Predictor-corrector continuation is used in conjunction with the indirect, rather than direct, method, because a predictor-corrector continuation direct method requires a predictor-corrector continuation NLP solver. Even though predictor-corrector continuation NLP solver algorithms are provided in [16, 17], there do not seem to be any publicly available predictor-corrector continuation NLP solvers.

2 Trajectory Tracking for the Rolling Disk

@fb@secFB

2.1 Optimal Control Problem and Controlled Equations of Motion

In the next subsection, numerical solutions of the controlled equations of motion for the rolling disk are presented, where the goal is to move the disk between a pair of points while the disk’s GC tracks a prescribed trajectory. A rolling disk of mass $m_{0}=4$, radius $r=1$, polar moment of inertia $d_{2}=1$, and with the CM coinciding with the GC (i.e. $\boldsymbol{\zeta}_{0}=\mathbf{0}$) is simulated. These physical parameters for the disk are consistent with the necessary and sufficient conditions stipulated by Inequality 3.2 in [18]; in that inequality, note that $d_{1}+d_{3}=d_{2}=1$ since the disk is a planar distribution of mass. There are $n=4$ control masses, each of mass $1$ so that $m_{1}=m_{2}=m_{3}=m_{4}=1$, located on concentric circular control rails centered on the GC of radii $r_{1}=.9$, $r_{2}=.6\overline{3}$, $r_{3}=.3\overline{6}$, and $r_{4}=.1$, as shown in Figure 2.1. For $1\leq i\leq 4$, the position of $m_{i}$ in the body frame centered on the GC is

[TABLE]

The total mass of the system is $M=8$, and gravity is $g=9.81$. There is no external force acting on the disk’s GC, so that $F_{\mathrm{e},1}=0$ in (1.8). The initial time is fixed to $a=0$ and the final time is fixed to $b=12$. The disk’s GC starts at rest at $z_{a}=0$ at time $a=0$ and stops at rest at $z_{b}=1$ at time $b=12$. Table 2.1 shows parameter values used in the rolling disk’s initial conditions (2.10) and final conditions (2.11). Since the initial orientation of the disk is $\phi_{a}=0$ and since the initial configurations of the control masses are given by \boldsymbol{\theta}_{a}=\begin{bmatrix}\scalebox{0.75}[1.0]{-}\frac{\pi}{2}&\scalebox{0.75}[1.0]{-}\frac{\pi}{2}&\scalebox{0.75}[1.0]{-}\frac{\pi}{2}&\scalebox{0.75}[1.0]{-}\frac{\pi}{2}\end{bmatrix}^{\mathsf{T}}, all the control masses are initially located directly below the GC. In order for the disk to start and stop at rest, $\dot{\boldsymbol{\theta}}_{a}=\dot{\boldsymbol{\theta}}_{b}=\begin{bmatrix}0&0&0&0\end{bmatrix}^{\mathsf{T}}$ and $\dot{z}_{a}=\dot{z}_{b}=0$. Table 2.2 shows parameter values used in the rolling disk’s final conditions (2.11).

The desired GC path $z_{\mathrm{d}}$ is depicted by the red curve in Figures 2.2a and 2.2b. $z_{\mathrm{d}}$ encourages the disk’s GC to track a sinusoidally-modulated linear trajectory connecting $z(0)=0$ with $z(12)=1$. That is, the disk is encouraged to roll right, then left, then right, then left, and finally to the right, with the amplitude of each successive roll increasing from the previous one. Specifically, $z_{\mathrm{d}}$ is given by

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

with $\epsilon=.01$ in (2.3). The reader is referred to [1] for further details about the properties and construction of (2.2).

The optimal control problem for the rolling disk is

[TABLE]

In (2.7), the system state ${\boldsymbol{x}}$ and control $\boldsymbol{u}$ are

[TABLE]

where $\boldsymbol{\theta},\dot{\boldsymbol{\theta}},\ddot{\boldsymbol{\theta}}\in\mathbb{R}^{4}$ and $\phi,\dot{\phi}\in\mathbb{R}$. In (2.7), the system dynamics defined for $a\leq t\leq b$ are

[TABLE]

where $\kappa\left({\boldsymbol{x}},\boldsymbol{u}\right)$ is given by the right-hand side of (1.8). In (2.9), the time-dependence of $\kappa$ is dropped since $F_{\mathrm{e},1}=0$ in (1.8) for these simulations. In (2.7), the prescribed initial conditions at time $t=a$ are

[TABLE]

and the prescribed final conditions at time $t=b$ are

[TABLE]

In (2.11),

[TABLE]

is the total system CM expressed in the body frame translated to the disk’s GC at time $t$,

[TABLE]

is the rotation matrix that maps the body to spatial frame at time $t$, and $\Pi$ is the projection onto the first component. Therefore, the first constraint in (2.11) ensures that the total system CM is above or below the disk’s GC in the spatial frame at the final time $t=b$, so that, in conjunction with the final condition parameter values given in Table 2.2, the disk stops at rest. In (2.7), the performance index is

[TABLE]

where the integrand cost function is

[TABLE]

for positive coefficients $\alpha(\mu)$ and $\gamma_{i}$, $1\leq i\leq 4$. The first summand $\frac{\alpha(\mu)}{2}\left(z_{a}-r(\phi-\phi_{a})-z_{\mathrm{d}}\right)^{2}$ in $L$ encourages the disk’s GC to track the desired spatial $\mathbf{e}_{1}$ path $z_{\mathrm{d}}$, and $\mu$ is a scalar continuation parameter used to construct a sequence of optimal control problems. The next $4$ summands $\frac{\gamma_{i}}{2}{\ddot{\theta}}_{i}^{2}$, $1\leq i\leq 4$, in $L$ limit the magnitude of the acceleration of the $i^{\mathrm{th}}$ control mass parameterization. Table 2.3 shows the values set for the integrand cost function coefficients in (2.15).

As explained in Appendix A, the controlled equations of motion for the rolling disk’s optimal control problem (2.7) are encapsulated by the ODE TPBVP:

[TABLE]

Subappendix A.1 of [1] derives the formulas for constructing $H_{{\boldsymbol{x}}}^{\mathsf{T}}$. In (2.16), $G$ is the endpoint function

[TABLE]

where $\boldsymbol{\xi}\in\mathbb{R}^{10}$ and $\boldsymbol{\nu}\in\mathbb{R}^{7}$ are constant Lagrange multiplier vectors enforcing the initial and final conditions, (2.10) and (2.11), respectively. In (2.16), $H$ is the Hamiltonian

[TABLE]

where $\boldsymbol{\lambda}\in\mathbb{R}^{10}$ is a time-varying Lagrange multiplier vector enforcing the dynamics (2.9). In (2.16), $\boldsymbol{\pi}$ is an analytical formula expressing the control $\boldsymbol{u}$ as a function of the state ${\boldsymbol{x}}$ and the costate $\boldsymbol{\lambda}$. The components of $\boldsymbol{\pi}$ are given by

[TABLE]

for $1\leq i\leq 4$. In (2.16), $\hat{H}$ is the regular Hamiltonian

[TABLE]

The reader is referred to [1] for a more general description of the rolling disk’s optimal control problem (2.7) and the associated controlled equations of motion (2.16).

@fb@secFB

2.2 Numerical Solutions: Trajectory Tracking

The direct method solver PS-II is used to solve the optimal control problem (2.7) when the integrand cost function coefficient is $\alpha=.1$. Predictor-corrector continuation is then used to solve the controlled equations of motion (2.16), starting from the direct method solution. The continuation parameter is $\mu$, which is used to adjust $\alpha$ according to the linear homotopy given in Table 2.3, so that $\alpha=.1$ when $\mu=.95$ and $\alpha\approx 272$ when $\mu\approx.8983$. The predictor-corrector continuation begins at $\mu=.95$, which is consistent with the direct method solution obtained at $\alpha=.1$.

For the direct method, PS-II is run using the NLP solver SNOPT. The PS-II mesh error tolerance is 1\mathrm{e}\scalebox{0.75}[1.0]{-}6 and the SNOPT error tolerance is 1\mathrm{e}\scalebox{0.75}[1.0]{-}7. In order to encourage convergence of SNOPT, a constant $C=50$ is added to the integrand cost function $L$ in (2.15). The sweep predictor-corrector continuation method discussed in Appendix D is used by the indirect method. For the sweep predictor-corrector continuation method, the maximum tangent steplength $\sigma_{\mathrm{max}}$ is adjusted according to Figure 2.4d over the course of $6$ predictor-corrector steps, the maximum tangent steplength in each step is $\sigma_{\mathrm{max}}=\begin{bmatrix}40&40&40&5&1&1\end{bmatrix}$, the direction of the initial unit tangent is determined by setting d=\scalebox{0.75}[1.0]{-}2 to force the continuation parameter $\mu$ to initially decrease, the relative error tolerance is 1\mathrm{e}\scalebox{0.75}[1.0]{-}8, the unit tangent solver is bvpc_m, and the monotonic “sweep” continuation solver is cc. The numerical results are shown in Figures 2.2, 2.3, and 2.4. As $\mu$ decreases from $.95$ down to $.8983$ during continuation (see Figure 2.4a), $\alpha$ increases from $.1$ up to $272$ (see Figure 2.4b). Since $\alpha$ is ratcheted up during continuation, thereby increasing the penalty in the integrand cost function (2.15) for deviation between the disk’s GC and $z_{\mathrm{d}}$, by the end of continuation, the disk’s GC tracks $z_{\mathrm{d}}$ very accurately (compare Figures 2.2a vs 2.2b), at the expense of more serpentine control mass trajectories (compare Figures 2.2c vs 2.2d) and larger magnitude controls (compare Figures 2.2e vs 2.2f). The disk does not detach from the surface since the magnitude of the normal force is always positive (see Figures 2.3a and 2.3b). The disk rolls without slipping if the coefficient of static friction $\mu_{\mathrm{s}}$ is at least $\hat{\mu}_{\mathrm{s}}\approx.07799$ for the direct method solution (see Figure 2.3c) and if $\mu_{\mathrm{s}}$ is at least $\hat{\mu}_{\mathrm{s}}\approx.3502$ for the indirect method solution (see Figure 2.3d). That is, the indirect method solution requires a much larger coefficient of static friction.

2.2.1 Turning Points

The previous sweep predictor-corrector continuation indirect method is repeated, but this time with 22 steps, where the maximum tangent steplength in each step is

[TABLE]

Note that the first 6 maximum tangent steplengths in (2.21) agree with those used in the previous simulation, so that the two simulations agree for the first 7 solutions. Figure 2.5 shows the evolution of the continuation parameter $\mu$, GC path weighting factor $\alpha$, performance index $J$, GC tracking error $\left\|z-z_{\mathrm{d}}\right\|^{2}$, and tangent steplength $\sigma$ over the course of the 23 solutions (the first solution is initialized by the direct method) constructed by the sweep predictor-corrector continuation indirect method. Note the turning points (local maxima or minima) at solutions 7, 10, and 18 in Figures 2.5a-2.5d. Figure 2.5d shows that the GC tracking error realizes a minimum at solution 7, which explains how the stopping point for the previous simulation was selected.

3 Obstacle Avoidance for the Rolling Ball

@fb@secFB

3.1 Optimal Control Problem and Controlled Equations of Motion

In the next two subsections, numerical solutions of the controlled equations of motion for the rolling ball are presented, where the goal is to move the ball between a pair of points while the ball’s GC avoids a pair of obstacles. A rolling ball of mass $m_{0}=4$, radius $r=1$, principal moments of inertia $d_{1}=d_{2}=d_{3}=1$, and with the CM coinciding with the GC (i.e. $\boldsymbol{\zeta}_{0}=\mathbf{0}$) is simulated. These physical parameters for the ball are consistent with the necessary and sufficient conditions stipulated by Inequalities 3.1 and 3.2 in [18]. There are $n=3$ control masses, each of mass $1$ so that $m_{1}=m_{2}=m_{3}=1$, located on circular control rails centered on the GC of radii $r_{1}=.95$, $r_{2}=.9$, and $r_{3}=.85$, oriented as shown in Figure 3.1. For $1\leq i\leq 3$, the position of $m_{i}$ in the body frame centered on the GC is

[TABLE]

where $\mathcal{B}_{i}\left(\mathbf{n}\right)\in SO(3)$ is a rotation matrix whose columns are the right-handed orthonormal basis constructed from the unit vector $\mathbf{n}\in\mathbb{R}^{3}$ based on the algorithm given in Section 4 and Listing 2 of [19], $\boldsymbol{\varsigma}\colon\mathbb{R}^{3}\to\mathbb{R}^{3}$ maps spherical coordinates to Cartesian coordinates:

[TABLE]

and

[TABLE]

are spherical coordinates of unit vectors in $\mathbb{R}^{3}$.

The total mass of the system is $M=7$, and gravity is $g=9.81$. There is no external force acting on the ball’s GC, so that $\tilde{\boldsymbol{\Gamma}}\equiv\Lambda^{-1}\mathbf{F}_{\mathrm{e}}=\mathbf{0}$ in (1.3). The initial time is fixed to $a=0$ and the final time is fixed to $b=5$. The ball’s GC starts at rest at $\boldsymbol{z}_{a}=\begin{bmatrix}0&0\end{bmatrix}^{\mathsf{T}}$ at time $a=0$ and stops at rest at $\boldsymbol{z}_{b}=\begin{bmatrix}1&1\end{bmatrix}^{\mathsf{T}}$ at time $b=5$. The ball’s GC should avoid $2$ circular obstacles, depicted in Figures 3.2a and 3.2b. The obstacles each have radius $\rho_{1}=\rho_{2}=.282$ and are centered at $\boldsymbol{v}_{1}=\begin{bmatrix}.2&.2\end{bmatrix}^{\mathsf{T}}$ and $\boldsymbol{v}_{2}=\begin{bmatrix}.8&.8\end{bmatrix}^{\mathsf{T}}$.

The ODE formulation of the optimal control problem for the rolling ball is

[TABLE]

In (3.4), the system state ${\boldsymbol{x}}$ and control $\boldsymbol{u}$ are

[TABLE]

where $\boldsymbol{\theta},\,\dot{\boldsymbol{\theta}},\,\ddot{\boldsymbol{\theta}}\in\mathbb{R}^{3}$, $\mathfrak{q}\in\mathscr{S}\cong\mathbb{S}^{3}\subset\mathbb{R}^{4}$, $\boldsymbol{\Omega}\in\mathbb{R}^{3}$, and $\boldsymbol{z}\in\mathbb{R}^{2}$. In the system state defined in (3.5), the rolling ball’s orientation matrix $\Lambda\in SO(3)$ is parameterized by $\mathfrak{q}\in\mathscr{S}$, where $\mathscr{S}$ denotes the set of versors (i.e. unit quaternions) [4, 20, 21, 22]. The properties of versors and the notation used to manipulate versors are explained in Appendix D of [2]. Recall from [2] that given a column vector $\boldsymbol{v}\in\mathbb{R}^{3}$, $\boldsymbol{v}^{\sharp}$ is the quaternion

[TABLE]

and given a quaternion $\mathfrak{p}\in\mathbb{H}$, $\mathfrak{p}^{\flat}\in\mathbb{R}^{3}$ is the column vector such that

[TABLE]

As explained in [2], the transformation of a body frame vector $\mathbf{Y}\in\mathbb{R}^{3}$ into the spatial frame by the ball’s orientation matrix $\Lambda\in SO(3)$ can be realized using the versor $\mathfrak{q}\in\mathscr{S}$ via the Euler-Rodrigues formula

[TABLE]

In (3.4), the system dynamics defined for $a\leq t\leq b$ are

[TABLE]

where $\boldsymbol{\kappa}\left({\boldsymbol{x}},\boldsymbol{u}\right)$ is given by the right-hand side of the formula for $\dot{\boldsymbol{\Omega}}$ in (1.3). In (3.9), the time-dependence of $\boldsymbol{\kappa}$ is dropped since $\tilde{\boldsymbol{\Gamma}}\equiv\Lambda^{-1}\mathbf{F}_{\mathrm{e}}=\mathbf{0}$ in (1.3) for these simulations. In (3.4), the prescribed initial conditions at time $t=a$ are

[TABLE]

and the prescribed final conditions at time $t=b$ are

[TABLE]

Table 3.4 shows the parameter values used in the rolling ball’s initial conditions (3.10). The initial configurations of the control masses are selected so that the total system CM in the spatial frame is initially located above or below the ball’s GC. Hence, in conjunction with the other initial condition parameter values given in Table 3.4, the ball starts at rest. In (3.11),

[TABLE]

is the total system CM in the body frame translated to the ball’s GC at time $t$,

[TABLE]

is the total system CM in the spatial frame translated to the ball’s GC at time $t$, and $\boldsymbol{\Pi}$ is the projection onto the first two components. Therefore, the first two constraints in (3.11) ensure that the total system CM in the spatial frame is above or below the ball’s GC at the final time $t=b$. Hence, in conjunction with the final condition parameter values given in Table 3.5, the ball stops at rest. In (3.4), the performance index is

[TABLE]

where the integrand cost function is

[TABLE]

for positive coefficients $\gamma_{i}(\mu)$, $1\leq i\leq 3$, and $h_{j}(\mu)$, $1\leq j\leq 2$. The first $3$ summands $\frac{\gamma_{i}(\mu)}{2}{\ddot{\theta}}_{i}^{2}$, $1\leq i\leq 3$, in $L$ limit the magnitude of the acceleration of the $i^{\mathrm{th}}$ control mass parameterization and the final $2$ summands $h_{j}(\mu)\,S\left(\left|\boldsymbol{z}-\boldsymbol{v}_{j}\right|-\rho_{j}\right)$, $1\leq j\leq 2$, in $L$ encourage the ball’s GC to avoid the pair of obstacles. For the obstacle avoidance function in $L$, $S$ is either the time-reversed sigmoid function (2.3) or the $C^{2}$ cutoff function

[TABLE]

In (3.16), $\operatorname{ReLU}$ is the rectified linear unit function frequently used in the machine learning literature [23]. In (3.15), the coefficients $\gamma_{i}(\mu)$, $1\leq i\leq 3$, and $h_{j}(\mu)$, $1\leq j\leq 2$, depend on the scalar continuation parameter $\mu$ so that a sequence of optimal control problems may be constructed. Note that a solution obtained by the optimal control procedure, that minimizes (2.14) for the rolling disk or (3.14) for the rolling ball, is a “compromise” between several, often conflicting, components, where some components of the performance index can be made more prominent by making their coefficients appropriately larger. The minimization of the performance index does not guarantee the minimization of each component individually.

There is also a DAE formulation of the optimal control problem for the rolling ball which explicitly enforces the algebraic versor constraint on $\mathfrak{q}$ and which is mathematically equivalent to (3.4). In the DAE formulation, the first component, $q_{0}$, of the versor $\mathfrak{q}$ is moved from the state ${\boldsymbol{x}}$ to the control $\boldsymbol{u}$ and an imitator state, $\tilde{q}_{0}$, is used to replace $q_{0}$ in ${\boldsymbol{x}}$. $\tilde{q}_{a,0}=q_{a,0}$, so that with perfect integration (i.e. no numerical integration errors), $\tilde{q}_{0}(t)=q_{0}(t)$ for $a\leq t\leq b$. Defining

[TABLE]

then with perfect integration,

[TABLE]

for $a\leq t\leq b$. $\tilde{q}_{0}$ is added to the state since the final conditions require knowledge of $q_{0}$, which is unavailable if it has been moved to the control since the final conditions are not a function of the control. The DAE formulation of the rolling ball’s optimal control problem is

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Even though the DAE formulation (3.19) is mathematically equivalent to the ODE formulation (3.4), the DAE formulation (3.19) tends to be numerically more stable to solve than the ODE formulation (3.4), as explained in Example 6.12 “Reorientation of an Asymmetric Rigid Body” of [7].

As explained in Appendix A, the controlled equations of motion for the ODE formulaton (3.4) of the rolling ball’s optimal control problem are encapsulated by the ODE TPBVP:

[TABLE]

Subappendix A.2 of [1] derives the formulas for constructing $H_{{\boldsymbol{x}}}^{\mathsf{T}}$. In (3.23), $G$ is the endpoint function

[TABLE]

where $\boldsymbol{\xi}\in\mathbb{R}^{15}$ and $\boldsymbol{\nu}\in\mathbb{R}^{10}$ are constant Lagrange multiplier vectors enforcing the initial and final conditions, (3.10) and (3.11), respectively. In (3.23), $H$ is the Hamiltonian

[TABLE]

where $\boldsymbol{\lambda}\in\mathbb{R}^{15}$ is a time-varying Lagrange multiplier vector enforcing the dynamics (3.9). In (3.23), $\boldsymbol{\pi}$ is an analytical formula expressing the control $\boldsymbol{u}$ as a function of the state ${\boldsymbol{x}}$ and the costate $\boldsymbol{\lambda}$. The components of $\boldsymbol{\pi}$ are given by

[TABLE]

for $1\leq i\leq 3$ and where $\boldsymbol{\lambda}_{\boldsymbol{\Omega}}\equiv\begin{bmatrix}\lambda_{11}&\lambda_{12}&\lambda_{13}\end{bmatrix}^{\mathsf{T}}$. In (3.23), $\hat{H}$ is the regular Hamiltonian

[TABLE]

The reader is referred to [1] for a more general description of the ODE and DAE formulations, (3.4) and (3.19), of the rolling ball’s optimal control problem and the controlled equations of motion (3.23) correpsonding to (3.4). The DAE TPBVP encapsulating the controlled equations of motion corresponding to the DAE formulation (3.19) of the rolling ball’s optimal control problem were not investigated since a robust DAE TPBVP solver is not readily available in LAB. COLDAE is a robust DAE TPBVP solver that uses collocation [24]; however, COLDAE is only available in Fortran and thus was not used in our calculations.

@fb@secFB

3.2 Numerical Solutions: Sigmoid Obstacle Avoidance

The controlled equations of motion (3.23) for the rolling ball are solved numerically to move the ball between the pair of points while avoiding the pair of obstacles, where the obstacle avoidance function $S$ in (3.15) is realized via the time-reversed sigmoid function (2.3) with $\epsilon=.01$. The direct method solver PS-II is used to solve the DAE formulation (3.19) of the optimal control problem, where the obstacle heights appearing in the integrand cost function (3.15) are $h_{1}=h_{2}=0$ and where the values for the other integrand cost function coefficients are given in Table 3.6. Using this direct method solution as an initial guess, PS-II is used again to solve the ODE formulation (3.4) of the same optimal control problem. Predictor-corrector continuation is then used to solve the controlled equations of motion (3.23), starting from the second direct method solution. The continuation parameter is $\mu$, which is used to adjust $h_{1}=h_{2}$ according to the linear homotopy shown in Table 3.6, so that $h_{1}=h_{2}=0$ when $\mu=.95$ and $h_{1}=h_{2}=1{,}000$ when $\mu=.00001$. The predictor-corrector continuation begins at $\mu=.95$, which is consistent with the direct method solution obtained at $h_{1}=h_{2}=0$.

For the direct method, PS-II is run using the NLP solver SNOPT. The PS-II mesh error tolerance is 1\mathrm{e}\scalebox{0.75}[1.0]{-}6 and the SNOPT error tolerance is 1\mathrm{e}\scalebox{0.75}[1.0]{-}7. In order to encourage convergence of SNOPT, a constant $C=50$ is added to the integrand cost function $L$ in (3.15). The sweep predictor-corrector continuation method discussed in Appendix D is used by the indirect method. For the sweep predictor-corrector continuation method, there are $4$ predictor-corrector steps, the maximum tangent steplength in each step is $\sigma_{\mathrm{max}}=500$, the direction of the initial unit tangent is determined by setting d=\scalebox{0.75}[1.0]{-}2 to force the continuation parameter $\mu$ to initially decrease, the relative error tolerance is 1\mathrm{e}\scalebox{0.75}[1.0]{-}6, the unit tangent solver is bvpc_m, and the monotonic “sweep” continuation solver is cc. The numerical results are shown in Figures 3.2, 3.3, and 3.4. As $\mu$ decreases from $.95$ down to $.9406$ during continuation (see Figure 3.4a), $h_{1}=h_{2}$ increases from [math] up to $9.93$ (see Figure 3.4b). Since $h_{1}=h_{2}$ is ratcheted up during continuation, thereby increasing the penalty in the integrand cost function (3.15) when the GC intrudes into the obstacles, by the end of continuation, the ball’s GC avoids both obstacles while veering smartly around the first obstacle (compare Figures 3.2a vs 3.2b), at the expense of slightly larger magnitude controls (compare Figures 3.2e vs 3.2f). The ball does not detach from the surface since the magnitude of the normal force is always positive (see Figures 3.3a and 3.3b). The ball rolls without slipping if the coefficient of static friction $\mu_{\mathrm{s}}$ is at least $\hat{\mu}_{\mathrm{s}}\approx.1055$ for the direct method solution (see Figure 3.3c) and if $\mu_{\mathrm{s}}$ is at least $\hat{\mu}_{\mathrm{s}}\approx.0988$ for the indirect method solution (see Figure 3.3d). As shown in Figures 3.4a-3.4c, the sweep predictor-corrector continuation indirect method encounters turning points at solutions 3 and 4.

@fb@secFB

3.3 Numerical Solutions: $\operatorname{ReLU}$ Obstacle Avoidance

The controlled equations of motion (3.23) for the rolling ball are solved numerically again to move the ball between the pair of points while avoiding the pair of obstacles, but this time the obstacle avoidance function $S$ in (3.15) is realized via the $C^{2}$ cutoff function (3.16). With the obstacle heights appearing in the integrand cost function (3.15) set to $h_{1}=h_{2}=0$ and with the values for the other integrand cost function coefficients set according to Table 3.6, the same double pass (DAE formulation followed by ODE formulation) direct method, with all the same settings, parameters, and initial conditions, is used as in the previous subsection to generate the same initial solution. Starting from the direct method solution, two rounds of the same sweep predictor-corrector continuation method that was used in the previous subsection are again used here to solve the controlled equations of motion (3.23), with all the same settings and parameters, except for the number of predictor-corrector steps and maximum tangent steplengths used. In the first round, the continuation parameter $\mu$ is used to adjust $h_{1}=h_{2}$ according to the linear homotopy shown in Table 3.6, so that $h_{1}=h_{2}=0$ when $\mu=.95$ and $h_{1}=h_{2}=1{,}000$ when $\mu=.00001$. This predictor-corrector continuation begins at $\mu=.95$, which is consistent with the direct method solution obtained at $h_{1}=h_{2}=0$. In the first round, two predictor-corrector steps are made with the maximum tangent steplengths $\sigma_{\mathrm{max}}=\begin{bmatrix}20000&100000\end{bmatrix}$. Starting from the predictor-corrector continuation solution obtained in the first round, a second predictor-corrector continuation is used to solve the controlled equations of motion (3.23), where the continuation parameter is $\mu$, which is now used to adjust $\gamma_{1}=\gamma_{2}=\gamma_{3}$ according to the linear homotopy shown in Table 3.7, so that $\gamma_{1}=\gamma_{2}=\gamma_{3}=10$ when $\mu=.95$ and \gamma_{1}=\gamma_{2}=\gamma_{3}=\scalebox{0.75}[1.0]{-}1{,}000 when $\mu=.00001$. The second predictor-corrector continuation begins at $\mu=.95$, which is consistent with the first predictor-corrector continuation solution obtained at $\gamma_{1}=\gamma_{2}=\gamma_{3}=10$. Moreover, during the second predictor-corrector continuation, the obstacle heights $h_{1}=h_{2}$ are fixed at $7.846\mathrm{e}{8}$, which is consistent with the final obstacle heights obtained by the first predictor-corrector continuation solution. In the second continuation, four predictor-corrector steps are made with the maximum tangent steplengths $\sigma_{\mathrm{max}}=\begin{bmatrix}5000&200&1&1\end{bmatrix}$.

The numerical results are shown in Figures 3.5, 3.6, and 3.7. As $\mu$ decreases from $.95$ down to \scalebox{0.75}[1.0]{-}7.453\mathrm{e}{5} during the first round of predictor-corrector continuation, $h_{1}=h_{2}$ increases from [math] up to $7.846\mathrm{e}{8}$ (see Figure 3.7a). Since $h_{1}=h_{2}$ is ratcheted up during continuation, thereby increasing the penalty in the integrand cost function (3.15) when the GC intrudes into the obstacles, by the end of continuation, the ball’s GC completely exits the second obstacle and approaches the boundary of the first obstacle (compare Figures 3.5a vs 3.5b). As $\mu$ decreases from $.95$ down to $.9406$ during the second round of predictor-corrector continuation, $\gamma_{1}=\gamma_{2}=\gamma_{3}$ decreases from $10$ down to $3.602\mathrm{e}{-5}$ (see Figure 3.7b). Since $\gamma_{1}=\gamma_{2}=\gamma_{3}$ is ratcheted down during continuation, thereby decreasing the penalty in the integrand cost function (3.15) for large magnitude accelerations of the control mass parameterizations, and since the obstacle heights are held fixed at $7.846\mathrm{e}{8}$, by the end of continuation, the ball’s GC avoids both obstacles while veering smartly around the first obstacle (compare Figures 3.5b vs 3.5c). Figure 3.7c shows that the performance index $J$ increases from $284.2$ up to $289.7$ as the obstacle heights are ramped up in the first round of predictor-corrector continuation; Figure 3.7d shows that the performance index $J$ then decreases down to $250$ as the control coefficients are ramped down and the ball’s GC fully departs the tall obstacles in the second round of predictor-corrector continuation.

4 Summary, Discussion, and Future Work

The controlled equations of motion for the rolling disk and ball were solved numerically using predictor-corrector continuation, starting from an initial solution obtained via a direct method, to solve trajectory tracking problems for the rolling disk and obstacle avoidance problems for the rolling ball. These optimal control maneuvers were achieved by performing predictor-corrector continuation in weighting factors that scale penalty functions in the integrand cost function of the performance index.

This paper focused on the indirect, rather than direct, method to numerically solve the optimal control problems. Because the indirect and direct methods only converge to a local minimum solution near the initial guess, a robust continuation algorithm capable of handling turning points is needed to obtain indirect and direct method solutions of complicated, nonconvex optimal control problems. A continuation indirect method requires a continuation ODE or DAE TPBVP solver, while a continuation direct method requires a continuation NLP solver. Predictor-corrector continuation ODE TPBVP algorithms were presented in Appendices C and D and implemented in LAB to realize the continuation indirect method used to solve the rolling disk and ball optimal control problems. Even though predictor-corrector continuation NLP solver algorithms are provided in the literature (e.g. see [16, 17]), there do not seem to be any publicly available predictor-corrector continuation NLP solvers, which inhibited the use of a continuation direct method in this paper. When compared against the direct method, the indirect method suffers from two major deficiencies:

Unlike the direct method, the indirect method has a very small radius of convergence and therefore requires a very accurate initial solution guess [7, 6, 5]. Moreover, unlike the direct method, the indirect method requires a guess of the costates, which are unphysical. 2. 2.

Unlike the direct method, the indirect method is unable to construct the switching structure (i.e. the times when the states and/or controls enter and exit the boundary) of an optimal control problem having path inequality constraints.

Since predictor-corrector continuation was used in this paper, the first deficiency in the indirect method only applied when constructing the solution of the initial ODE TPBVP, and this deficiency was circumvented by using a direct method to solve the optimal control problem corresponding to that initial ODE TPBVP. To circumvent the second deficiency in the indirect method, path inequality constraints were incorporated into the optimal control problems as soft constraints through penalty functions in the integrand cost functions.

The predictor-corrector continuation methods presented in Appendices C and D work if the control can be expressed analytically as a function of the state and costate, e.g. if the Hamiltonian is quadratic in the control. That is, those methods perform continuation only in the state and costate after the original optimal control DAE TPBVP is transformed into an ODE TPBVP through elimination of the control. For more complicated Hamiltonians (such as when penalty functions are added to the integrand cost function to softly enforce bounded body frame accelerations of the control masses, the no-detachment constraint, and the no-slip constraint), numerical methods (such as Newton’s method) must be used to construct the control numerically from the state, costate, and a good initial guess of the control. In these cases, the predictor-corrector continuation method of Appendix C must be extended to perform continuation in the state, costate, and control by solving the optimal control DAE TPBVP. This will be investigated in subsequent work.

In future work, instead of using LAB, the simulation code could be reimplemented in the higher performance programming languages Julia or C++, while relying on Fortran routines like COLNEW [25], COLMOD [26], TWPBVP(C) [27, 28], TWPBVPL(C) [29, 30], and ACDC(C) [26, 15] to solve the underlying ODE TPBVPs, to obtain faster numerical results. Julia and C++ feature several mature and efficient automatic differentiation libraries [31] capable of constructing the Jacobians and Hessians needed by the ODE TPBVP solvers. In addition, a more efficient and robust predictor-corrector adaptive tangent steplength algorithm, such as described in [32, 33], could be implemented.

Another avenue for future investigation is to use a neighboring extremal optimal control (NEOC) method [34], which constructs a homotopy between the controlled equations of motion and their linearization about a nominal solution; however, the NEOC method in [34] could be made more robust by using predictor-corrector, rather than monotonic, continuation in the homotopy parameter. Yet another avenue for future investigation is to perform predictor-corrector continuation in a weighting factor that scales a term in the endpoint cost function measuring the deviation between the actual and prescribed final conditions.

Additionally, throughout the paper we have kept the initial and final times fixed. It would be interesting to perform additional studies when the time duration is free, as outlined in Appendix A, especially regarding problems of navigation over complex terrains and slippery surfaces. Complex terrains will affect both the uncontrolled equations of motion, due to gravity, and the performance index, for example, through potential energy penalty functions that discourage the ball from ascending steep slopes. If the substrate is slippery, for example, due to the presence of moisture, one can imagine situations where the valleys are wet and the slopes are dry and thus less slippery. Then, one can introduce an additional term in the performance index penalizing motion through the valleys where the possibility of a slip is high. The interplay between the terms in the performance index discouraging and encouraging motion up the slopes will give a very interesting control system. These and other interesting questions will be considered in future work.

Acknowledgements

We are indebted to our colleagues A.M. Bloch, D.M. de Diego, F. Gay-Balmaz, D.D. Holm, M. Leok, A. Lewis, T. Ohsawa, and D.V. Zenkov for useful and fruitful discussions. M.J. Weinsten provided copious advice on using the LAB automatic differentiation toolbox Gator and fixed numerous bugs in Gator that were revealed in the course of this research. A.V. Rao provided a free license to use the LAB direct method optimal control solver PS-II for some of this research. P. Tallapragada observed that the reaction forces exerted on the ball by the accelerating internal point masses may cause the ball to detach from the surface, which prompted the inclusion of plots depicting the normal force and minimum coefficient of static friction. G.M. Rozenblat pointed out the necessary and sufficient conditions that must be satisfied by a ball’s physical parameters in [18].

This research was partially supported by the NSERC Discovery Grant, the University of Alberta Centennial Fund, and the Alberta Innovates Technology Funding (AITF) which came through the Alberta Centre for Earth Observation Sciences (CEOS). S.M. Rogers also received support from the University of Alberta Doctoral Recruitment Scholarship, the FGSR Graduate Travel Award, the IGR Travel Award, the GSA Academic Travel Award, the AMS Fall Sectional Graduate Student Travel Grant, Target Corporation, and the Institute for Mathematics and its Applications at the University of Minnesota, Twin Cities.

Appendix A Optimal Control: Variational Pontryagin’s Minimum Principle

This appendix presents necessary conditions, called the variational Pontryagin’s minimum principle, which a solution to an optimal control problem lacking path inequality constraints must satisfy; there is a more general version of Pontryagin’s minimum principle that applies to optimal control problems possessing path inequality constraints. In this paper, these necessary conditions, in the context of describing the optimal control of the rolling ball, are referred to as the controlled equations of motion. In the literature, application of Pontryagin’s minimum principle to solve an optimal control problem is called the indirect method. Let $n,m\in\mathbb{N}$. Let $a$ be a prescribed or free initial time and let $k_{1}\in\mathbb{N}^{0}$ be such that $0\leq k_{1}\leq n$ if $a$ is prescribed and $1\leq k_{1}\leq n+1$ if $a$ is free. Let $b$ be a prescribed or free final time and let $k_{2}\in\mathbb{N}^{0}$ be such that $0\leq k_{2}\leq n$ if $b$ is prescribed and $1\leq k_{2}\leq n+1$ if $b$ is free. Suppose a dynamical system has state ${\boldsymbol{x}}\in\mathbb{R}^{n}$ and control $\boldsymbol{u}\in\mathbb{R}^{m}$ and the control $\boldsymbol{u}$ is sought that minimizes the performance index

[TABLE]

subject to the system dynamics defined for $a\leq t\leq b$

[TABLE]

the prescribed initial conditions at time $t=a$

[TABLE]

and the prescribed final conditions at time $t=b$

[TABLE]

$p$ is a scalar-valued function called the endpoint cost function, $L$ is a scalar-valued function called the integrand cost function, ${\boldsymbol{x}}$ and $\mathbf{f}$ are $n\times 1$ vector-valued functions, $\boldsymbol{u}$ is an $m\times 1$ vector-valued function, $\boldsymbol{\sigma}$ is a $k_{1}\times 1$ vector-valued function, and $\boldsymbol{\psi}$ is a $k_{2}\times 1$ vector-valued function. $\mu$ is a prescribed scalar parameter which may be exploited to numerically solve this problem via continuation. More concisely, this optimal control problem may be stated as

[TABLE]

Observe that the optimal control problem encapsulated by (A.5) ignores path inequality constraints such as $\mathbf{D}\left(t,{\boldsymbol{x}},\boldsymbol{u},\mu\right)\leq\mathbf{0}$, where $\mathbf{D}$ is an $r\times 1$ vector-valued function for $r\in\mathbb{N}^{0}$. Path inequality constraints can be incorporated into (A.5) as soft constraints through penalty functions in the integrand cost function $L$ or the endpoint cost function $p$. By omitting hard path inequality constraints from (A.5), a solution of (A.5) does not lie on the boundary of a compact set and the calculus of variations may be applied to derive necessary conditions, called the variational Pontryagin’s minimum principle, which a solution of (A.5) must satisfy.

Define the endpoint function $G$ and the Hamiltonian $H$ by

[TABLE]

where $\boldsymbol{\xi}$ is a $k_{1}\times 1$ constant Lagrange multiplier vector, $\boldsymbol{\nu}$ is a $k_{2}\times 1$ constant Lagrange multiplier vector, and $\boldsymbol{\lambda}$ is an $n\times 1$ time-varying Lagrange multiplier vector. In the literature, the time-varying Lagrange multiplier vector used to adjoin the system dynamics to the integrand cost function is often called the adjoint variable or the costate. Henceforth, the time-varying Lagrange multiplier vector is referred to as the costate and the elements in this vector are referred to as the costates. The necessary conditions [35] on ${\boldsymbol{x}}$, $\boldsymbol{\lambda}$, and $\boldsymbol{u}$ which a solution of (A.5) must satisfy are the DAEs defined for $a\leq t\leq b$

[TABLE]

the left boundary conditions defined at time $t=a$

[TABLE]

and the right boundary conditions defined at time $t=b$

[TABLE]

If the initial time $a$ is prescribed, then the left boundary condition $\left.H\right|_{t=a}=G_{a}$ is dropped. If the final time $b$ is prescribed, then the right boundary condition $\left.H\right|_{t=b}=-G_{b}$ is dropped. The necessary conditions (A.7), (A.8), and (A.9) constitute a differential-algebraic equation two-point boundary value problem (DAE TPBVP).

If $H_{\boldsymbol{u}\boldsymbol{u}}$ is nonsingular, then the optimal control problem is said to be regular or nonsingular; otherwise if $H_{\boldsymbol{u}\boldsymbol{u}}$ is singular, then the optimal control problem is said to be singular. If $H_{\boldsymbol{u}\boldsymbol{u}}$ is nonsingular, then by the implicit function theorem, the algebraic equation $H_{\boldsymbol{u}}=\mathbf{0}_{1\times m}$ in (A.7) guarantees the existence of a unique function, say $\boldsymbol{\pi}$, for which

[TABLE]

If $H_{\boldsymbol{u}\boldsymbol{u}}$ is nonsingular, it may be possible to solve the algebraic equation $H_{\boldsymbol{u}}=\mathbf{0}_{1\times m}$ in (A.7) analytically for $\boldsymbol{u}$ in terms of $t$, ${\boldsymbol{x}}$, $\boldsymbol{\lambda}$, and $\mu$ to construct $\boldsymbol{\pi}$ explicitly in (A.10); otherwise, the value $\boldsymbol{u}$ of $\boldsymbol{\pi}$ in (A.10) may be constructed numerically in an efficient and accurate manner (with quadratic convergence) via a few iterations of Newton’s method applied to $H_{\boldsymbol{u}}=\mathbf{0}_{1\times m}$ starting from an initial guess $\boldsymbol{u}_{0}$ of $\boldsymbol{u}$:

[TABLE]

Using (A.10), the Hamiltonian may be re-expressed as a function of $t$, ${\boldsymbol{x}}$, $\boldsymbol{\lambda}$, and $\mu$ via the regular or reduced Hamiltonian

[TABLE]

Note that by construction of $\boldsymbol{\pi}$,

[TABLE]

By using the definition (A.12) of the regular Hamiltonian $\hat{H}$, invoking (A.13), and defining

[TABLE]

it follows from the chain rule that

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

By using (A.10) to eliminate the algebraic equation $H_{\boldsymbol{u}}=\mathbf{0}_{1\times m}$ from (A.7), by plugging (A.17) and (A.16) into the right hand sides of the ODEs in (A.7), and by plugging the definition (A.12) into the left and right boundary conditions (A.8) and (A.9), the necessary conditions on ${\boldsymbol{x}}$ and $\boldsymbol{\lambda}$ which a solution of (A.5) must satisfy are the ODEs defined for $a\leq t\leq b$

[TABLE]

the left boundary conditions defined at time $t=a$

[TABLE]

and the right boundary conditions defined at time $t=b$

[TABLE]

If the initial time $a$ is prescribed, then the left boundary condition $\left.\hat{H}\right|_{t=a}=G_{a}$ is dropped. If the final time $b$ is prescribed, then the right boundary condition $\left.\hat{H}\right|_{t=b}=-G_{b}$ is dropped. The necessary conditions (A.19), (A.20), and (A.21) constitute an ODE TPBVP. Appendix B provides implementation details for numerically solving the ODE TPBVP (A.19), (A.20), and (A.21).

A solution of the DAE TPBVP (A.7), (A.8), and (A.9) or of the ODE TPBVP (A.19), (A.20), and (A.21) is said to be an extremal solution of the optimal control problem (A.5). Note that an extremal solution only satisfies necessary conditions for a minimum of the optimal control problem (A.5), so that an extremal solution is not guaranteed to be a local minimum of (A.5).

Since the DAE TPBVP (A.7), (A.8), and (A.9) and the ODE TPBVP (A.19), (A.20), and (A.21) have small convergence radii [5, 6, 7], a continuation method (performing continuation in the parameter $\mu$) is often required to numerically solve them starting from a solution to a simpler optimal control problem [36]. The solution to the simpler optimal control problem might be obtained via analytics, the gradient method [5, 6], the method of successive approximations [37, 38, 39, 40], or the direct method [7, 8]. For example in [41], the continuation parameter $\mu$ is used to vary integrand cost function coefficients in $L$ in order to numerically solve the optimal control ODE TPBVPs for Suslov’s problem via monotonic continuation, starting from an analytical solution to a singular optimal control problem. In Sections 2 and 3, the continuation parameter $\mu$ is used to vary integrand cost function coefficients in $L$ in order to numerically solve the optimal control ODE TPBVPs for the rolling disk and ball via predictor-corrector continuation, starting from a direct method solution to a simpler optimal control problem. Appendices C and D describe predictor-corrector continuation methods for solving ODE TPBVPs and which are used to solve the optimal control ODE TPBVPs for the rolling disk and ball in Sections 2 and 3.

Appendix B Implementation Details for Solving the ODE TPBVP for a Regular Optimal Control Problem

Details for numerically solving the ODE TPBVP (A.19), (A.20), and (A.21) associated with the indirect method solution of a regular optimal control problem are presented here. There are two general methods, initial value and global, for numerically solving an ODE TPBVP. An initial value method, such as single or multiple shooting, subdivides the integration interval $[a,b]$ into a fixed, finite mesh and integrates the ODE on each mesh subinterval using a guess of the unknown initial conditions at one endpoint in each mesh subinterval. A root-finder is used to iteratively adjust the guesses of the unknown initial conditions until the solution segments are continuous at the internal mesh points and until the boundary conditions at the endpoints $a$ and $b$ are satisfied. A global method, such as a Runge-Kutta, collocation, or finite-difference scheme, subdivides the integration interval $[a,b]$ into a finite, adaptive mesh and solves a large nonlinear system of algebraic equations obtained by imposing the ODE constraints at a finite set of points in each mesh subinterval, by imposing continuity of the solution at internal mesh points, and by imposing the boundary conditions at the endpoints $a$ and $b$. By estimating the error in each mesh subinterval, a global method iteratively refines or adapts the mesh until a prescribed error tolerance is satisfied. Because initial value methods cannot integrate unstable ODEs, global methods are preferred [36, 42, 43].

@fb@secFB

B.1 Normalization and ODE Velocity Function

There are many solvers available to numerically solve the ODE TPBVP (A.19), (A.20), and (A.21). For example, 4c [44], 5c [45], p [14], and twp [15] (which encapsulates bvp_m, bvpc_m, bvp_l, bvpc_l, c, and cc) are LAB Runge-Kutta or collocation ODE TPBVP solvers, while COLSYS [46], COLNEW [25], COLMOD [26], COLCON [32], BVP_M-2 [47], TWPBVP [27], TWPBVPC [28], TWPBVPL [29], TWPBVPLC [30], ACDC [26], and ACDCC [15] are Fortran Runge-Kutta or collocation ODE TPBVP solvers. The reader is referred to the Appendix in [41] for a comprehensive list of ODE TPBVP solvers. In order to numerically solve the ODE TPBVP (A.19), (A.20), and (A.21), many solvers (such as the global method LAB and Fortran solvers just listed) require that the ODE TPBVP be defined on a fixed time interval and any unknown parameters, such as $\boldsymbol{\xi}$, $\boldsymbol{\nu}$, $a$, and $b$, must often be modeled as dummy constant dependent variables with zero derivatives. In addition, to aid convergence, many solvers can exploit Jacobians of the ODE velocity function and of the two-point boundary condition function. Thus, (A.19) is redefined on the normalized time interval $[0,1]$ through the change of independent variable $s\equiv\frac{t-a}{T}$, where $T\equiv b-a$. Note that $t(s)=Ts+a$. Define the normalized state $\tilde{\boldsymbol{x}}(s)\equiv{\boldsymbol{x}}(t(s))$ and normalized costate $\tilde{\boldsymbol{\lambda}}(s)\equiv\boldsymbol{\lambda}(t(s))$. Define the expanded un-normalized ODE TPBVP dependent variable vector

[TABLE]

Defining $\tilde{\boldsymbol{z}}(s)\equiv\boldsymbol{z}(t(s))$, the expanded normalized ODE TPBVP dependent variable vector is

[TABLE]

By the chain rule, (A.19), and since $\frac{\mathrm{d}t(s)}{\mathrm{d}s}=T$,

[TABLE]

Define $\tilde{\boldsymbol{\Phi}}\left(s,\tilde{\boldsymbol{z}}(s),\mu\right)$ to be the right-hand side of (B.3), i.e. the normalized ODE velocity function, so that

[TABLE]

The Jacobian of $\tilde{\boldsymbol{\Phi}}$ with respect to $\tilde{\boldsymbol{z}}(s)$ is

[TABLE]

and the Jacobian of $\tilde{\boldsymbol{\Phi}}$ with respect to $\mu$ is

[TABLE]

In (B.5) and (B.6), shorthand notation is used for conciseness and all zeroth and first derivatives of $\hat{\mathbf{f}}$ and all first and second derivatives of $\hat{H}$ are evaluated at $\left(s,\tilde{\boldsymbol{z}}(s),\mu\right)$. An explanation of the meaning of the shorthand notation used to express all zeroth and first derivatives of $\hat{\mathbf{f}}$ and all first and second derivatives of $\hat{H}$ is given in Table B.8. In rows $n+1$ through $2n$ and columns $n+1$ through $2n$ of (B.5), Clairaut’s Theorem was used to obtain $\hat{H}_{{\boldsymbol{x}}\boldsymbol{\lambda}}=\hat{H}_{\boldsymbol{\lambda}{\boldsymbol{x}}}^{\mathsf{T}}=\hat{\mathbf{f}}_{{\boldsymbol{x}}}^{\mathsf{T}}$, recalling from (A.17) that $\hat{H}_{\boldsymbol{\lambda}}=\hat{\mathbf{f}}^{\mathsf{T}}$.

Recall that $\hat{\mathbf{f}}$ is defined in (A.14) in terms of $\mathbf{f}$ and $\boldsymbol{\pi}$. By using the chain rule, the first derivatives of $\hat{\mathbf{f}}$ that appear in (B.5), (B.6), and Table B.8 may be computed from first derivatives of $\mathbf{f}$ and $\boldsymbol{\pi}$ as follows:

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Recall that by construction of $\boldsymbol{\pi}$, $H_{\boldsymbol{u}}\left(t,{\boldsymbol{x}},\boldsymbol{\lambda},\boldsymbol{\pi}\left(t,{\boldsymbol{x}},\boldsymbol{\lambda},\mu\right),\mu\right)=\mathbf{0}_{1\times m}$, as stated previously in (A.13). Differentiating $H_{\boldsymbol{u}}\left(t,{\boldsymbol{x}},\boldsymbol{\lambda},\boldsymbol{\pi}\left(t,{\boldsymbol{x}},\boldsymbol{\lambda},\mu\right),\mu\right)=\mathbf{0}_{1\times m}$ with respect to $\boldsymbol{\lambda}$, ${\boldsymbol{x}}$, $t$, and $\mu$, in turn, and using the chain rule gives

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

(B.11), (B.12), (B.13), and (B.14) may be solved for $\boldsymbol{\pi}_{\boldsymbol{\lambda}}$, $\boldsymbol{\pi}_{{\boldsymbol{x}}}$, $\boldsymbol{\pi}_{t}$, and $\boldsymbol{\pi}_{\mu}$, respectively:

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

In (B.15), Clairaut’s Theorem was used to obtain $H_{\boldsymbol{u}\boldsymbol{\lambda}}=H_{\boldsymbol{\lambda}\boldsymbol{u}}^{\mathsf{T}}={\mathbf{f}}_{\boldsymbol{u}}^{\mathsf{T}}$, since $H_{\boldsymbol{\lambda}}=\mathbf{f}^{\mathsf{T}}$. As will be stated again later, (B.15), (B.16), (B.17), and (B.18) are especially useful if the value of $\boldsymbol{\pi}$ is constructed numerically via Newton’s method as in (A.11); if $\boldsymbol{\pi}$ is given analytically, then it should be possible to construct $\boldsymbol{\pi}_{\boldsymbol{\lambda}}$, $\boldsymbol{\pi}_{{\boldsymbol{x}}}$, $\boldsymbol{\pi}_{t}$, and $\boldsymbol{\pi}_{\mu}$ via manual, symbolic, or automatic differentiation of the analytical formula for $\boldsymbol{\pi}$.

Since Clairaut’s Theorem guarantees that

[TABLE]

and

[TABLE]

(B.12) may be solved for $H_{{\boldsymbol{x}}\boldsymbol{u}}\left(t,{\boldsymbol{x}},\boldsymbol{\lambda},\boldsymbol{\pi}\left(t,{\boldsymbol{x}},\boldsymbol{\lambda},\mu\right),\mu\right)$:

[TABLE]

By differentiating (A.16) with respect to ${\boldsymbol{x}}$, $t$, and $\mu$, using the chain rule, and exploiting (B.21), the second derivatives of $\hat{H}$ that appear in (B.5), (B.6), and Table B.8 may be computed from first derivatives of $\mathbf{\boldsymbol{\pi}}$ and second derivatives of $H$ as follows:

[TABLE]

[TABLE]

and

[TABLE]

If the value of $\boldsymbol{\pi}$ is constructed numerically via Newton’s method as in (A.11) rather than analytically, then (B.15), (B.16), (B.17), and (B.18) should be used to evaluate $\boldsymbol{\pi}_{\boldsymbol{\lambda}}$, $\boldsymbol{\pi}_{{\boldsymbol{x}}}$, $\boldsymbol{\pi}_{t}$, and $\boldsymbol{\pi}_{\mu}$, which appear in the formulas for $\hat{\mathbf{f}}_{\boldsymbol{\lambda}}$ given in (B.7), $\hat{\mathbf{f}}_{{\boldsymbol{x}}}$ given in (B.8), $\hat{\mathbf{f}}_{t}$ given in (B.9), $\hat{\mathbf{f}}_{\mu}$ given in (B.10), $\hat{H}_{{\boldsymbol{x}}{\boldsymbol{x}}}$ given in (B.22), $\hat{H}_{{\boldsymbol{x}}t}$ given in (B.23), and $\hat{H}_{{\boldsymbol{x}}\mu}$ given in (B.24). The second equation in (B.22), (B.23), and (B.24) is given because it may be more computationally efficient than the first equation if $\boldsymbol{\pi}$ is given analytically, so that (B.15), (B.16), (B.17), and (B.18) need not be used to evaluate $\boldsymbol{\pi}_{\boldsymbol{\lambda}}$, $\boldsymbol{\pi}_{{\boldsymbol{x}}}$, $\boldsymbol{\pi}_{t}$, and $\boldsymbol{\pi}_{\mu}$.

@fb@secFB

B.2 Two-Point Boundary Condition Function

Now the boundary conditions (A.20)-(A.21) are considered. Letting

[TABLE]

[TABLE]

and

[TABLE]

the boundary conditions (A.20)-(A.21) in un-normalized dependent variables are given by the two-point boundary condition function

[TABLE]

The Jacobians of $\boldsymbol{\Upsilon}$ with respect to $\boldsymbol{z}(a)$, $\boldsymbol{z}(b)$, and $\mu$ are

[TABLE]

[TABLE]

and

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

In equations (B.32), (B.33), and (B.34), $\hat{\mathbf{f}}$ and all first derivatives of $\hat{H}$ in row $1$ are evaluated at $\left(a,{\boldsymbol{x}}(a),\boldsymbol{\lambda}(a),\mu\right)$ as shown in Table B.9, all first derivatives of $\boldsymbol{\sigma}$ in rows $n+2$ through $n+1+k_{1}$ are evaluated at $\left(a,{\boldsymbol{x}}(a),\mu\right)$ as shown in Table B.10, $\hat{\mathbf{f}}$ and all first derivatives of $\hat{H}$ in row $n+2+k_{1}$ are evaluated at $\left(b,{\boldsymbol{x}}(b),\boldsymbol{\lambda}(b),\mu\right)$ as shown in Table B.11, and all first derivatives of $\boldsymbol{\psi}$ in rows $2n+3+k_{1}$ through $2n+2+k_{1}+k_{2}$ are evaluated at $\left(b,{\boldsymbol{x}}(b),\mu\right)$ as shown in Table B.12. Since $\hat{H}_{\boldsymbol{\lambda}}=\hat{\mathbf{f}}^{\mathsf{T}}$, $\left.\hat{H}_{\boldsymbol{\lambda}}\right|_{a}=\left.\hat{\mathbf{f}}^{\mathsf{T}}\right|_{a}$ in row $1$ and columns $n+1$ through $2n$ of (B.32) and $\left.\hat{H}_{\boldsymbol{\lambda}}\right|_{b}=\left.\hat{\mathbf{f}}^{\mathsf{T}}\right|_{b}$ in row $n+2+k_{1}$ and columns $n+1$ through $2n$ of (B.33). In equations (B.35), (B.36), and (B.37), all second derivatives of $G$ are evaluated at $\left(a,{\boldsymbol{x}}(a),\boldsymbol{\xi},b,{\boldsymbol{x}}(b),\boldsymbol{\nu},\mu\right)$, while all first derivatives of $\boldsymbol{\sigma}$ and $\boldsymbol{\psi}$ are evaluated at $\left(a,{\boldsymbol{x}}(a),\mu\right)$ and $\left(b,{\boldsymbol{x}}(b),\mu\right)$, respectively, as shown in Tables B.10 and B.12. To simplify (B.35) and (B.36), Clairaut’s Theorem, $G_{\boldsymbol{\xi}}^{\mathsf{T}}=\boldsymbol{\sigma}$, and $G_{\boldsymbol{\nu}}^{\mathsf{T}}=\boldsymbol{\psi}$ are used to get $G_{a\boldsymbol{\xi}}=G_{\boldsymbol{\xi}a}^{\mathsf{T}}=\boldsymbol{\sigma}_{a}^{\mathsf{T}}$, $G_{a\boldsymbol{\nu}}=G_{\boldsymbol{\nu}a}^{\mathsf{T}}=\boldsymbol{\psi}_{a}^{\mathsf{T}}=\mathbf{0}_{1\times k_{2}}$, $G_{{\boldsymbol{x}}(a)\boldsymbol{\xi}}=G_{\boldsymbol{\xi}{\boldsymbol{x}}(a)}^{\mathsf{T}}=\boldsymbol{\sigma}_{{\boldsymbol{x}}(a)}^{\mathsf{T}}$, $G_{{\boldsymbol{x}}(a)\boldsymbol{\nu}}=G_{\boldsymbol{\nu}{\boldsymbol{x}}(a)}^{\mathsf{T}}=\boldsymbol{\psi}_{{\boldsymbol{x}}(a)}^{\mathsf{T}}=\mathbf{0}_{n\times k_{2}}$, $G_{b\boldsymbol{\xi}}=G_{\boldsymbol{\xi}b}^{\mathsf{T}}=\boldsymbol{\sigma}_{b}^{\mathsf{T}}=\mathbf{0}_{1\times k_{1}}$, $G_{b\boldsymbol{\nu}}=G_{\boldsymbol{\nu}b}^{\mathsf{T}}=\boldsymbol{\psi}_{b}^{\mathsf{T}}$, $G_{{\boldsymbol{x}}(b)\boldsymbol{\xi}}=G_{\boldsymbol{\xi}{\boldsymbol{x}}(b)}^{\mathsf{T}}=\boldsymbol{\sigma}_{{\boldsymbol{x}}(b)}^{\mathsf{T}}=\mathbf{0}_{n\times k_{1}}$, and $G_{{\boldsymbol{x}}(b)\boldsymbol{\nu}}=G_{\boldsymbol{\nu}{\boldsymbol{x}}(b)}^{\mathsf{T}}=\boldsymbol{\psi}_{{\boldsymbol{x}}(b)}^{\mathsf{T}}$.

To express the boundary conditions (A.20)-(A.21) in terms of normalized dependent variables, let $\tilde{\boldsymbol{\Upsilon}}_{1}\left(\tilde{\boldsymbol{z}}(0),\tilde{\boldsymbol{z}}(1),\mu\right)\equiv\boldsymbol{\Upsilon}_{1}\left(\boldsymbol{z}(a),\boldsymbol{z}(b),\mu\right)$, $\tilde{\boldsymbol{\Upsilon}}_{2}\left(\tilde{\boldsymbol{z}}(0),\tilde{\boldsymbol{z}}(1),\mu\right)\equiv\boldsymbol{\Upsilon}_{2}\left(\boldsymbol{z}(a),\boldsymbol{z}(b),\mu\right)$, and $\tilde{\boldsymbol{\Upsilon}}\left(\tilde{\boldsymbol{z}}(0),\tilde{\boldsymbol{z}}(1),\mu\right)\equiv\boldsymbol{\Upsilon}\left(\boldsymbol{z}(a),\boldsymbol{z}(b),\mu\right)$. Thus

[TABLE]

[TABLE]

and

[TABLE]

and the boundary conditions (A.20)-(A.21) in normalized dependent variables are given by the normalized two-point boundary condition function

[TABLE]

The Jacobians of $\tilde{\boldsymbol{\Upsilon}}$ with respect to $\tilde{\boldsymbol{z}}(0)$, $\tilde{\boldsymbol{z}}(1)$, and $\mu$ are

[TABLE]

[TABLE]

and

[TABLE]

where the equality between the Jacobians of $\tilde{\boldsymbol{\Upsilon}}$, $\tilde{\boldsymbol{\Upsilon}}_{1}$, and $\tilde{\boldsymbol{\Upsilon}}_{2}$ with respect to $\tilde{\boldsymbol{z}}(0)$, $\tilde{\boldsymbol{z}}(1)$, and $\mu$ and the Jacobians of $\boldsymbol{\Upsilon}$, $\boldsymbol{\Upsilon}_{1}$, and $\boldsymbol{\Upsilon}_{2}$ with respect to $\boldsymbol{z}(0)$, $\boldsymbol{z}(1)$, and $\mu$ is given in Table B.13.

Special care must be taken when implementing the Jacobians (B.42) and (B.43). Since the unknown constants $\boldsymbol{\xi}$, $\boldsymbol{\nu}$, $a$, and $b$ appear at the end of both $\tilde{\boldsymbol{z}}(0)$ and $\tilde{\boldsymbol{z}}(1)$, the unknown constants from only one of $\tilde{\boldsymbol{z}}(0)$ and $\tilde{\boldsymbol{z}}(1)$ are actually used to construct each term in $\tilde{\boldsymbol{\Upsilon}}$ involving $\boldsymbol{\xi}$, $\boldsymbol{\nu}$, $a$, and $b$. The trailing columns in (B.42) are actually the Jacobian of $\tilde{\boldsymbol{\Upsilon}}$ with respect to $\boldsymbol{\xi}$, $\boldsymbol{\nu}$, $a$, and $b$ in $\tilde{\boldsymbol{z}}(0)$, while the trailing columns in (B.43) are actually the Jacobian of $\tilde{\boldsymbol{\Upsilon}}$ with respect to $\boldsymbol{\xi}$, $\boldsymbol{\nu}$, $a$, and $b$ in $\tilde{\boldsymbol{z}}(1)$. Thus, the trailing columns in (B.42) and (B.43) corresponding to the Jacobian of $\tilde{\boldsymbol{\Upsilon}}$ with respect to $\boldsymbol{\xi}$, $\boldsymbol{\nu}$, $a$, and $b$ should not coincide in a software implementation. For example, if the unknown constants are extracted from $\tilde{\boldsymbol{z}}(0)$ to construct $\tilde{\boldsymbol{\Upsilon}}$, $\tilde{\boldsymbol{\Upsilon}}_{\tilde{\boldsymbol{z}}(0)}$ is as shown in (B.42) while the trailing columns in (B.43) corresponding to the Jacobian of $\tilde{\boldsymbol{\Upsilon}}$ with respect to the unknown constants in $\tilde{\boldsymbol{z}}(1)$ should be all zeros. Alternatively, if the unknown constants are extracted from $\tilde{\boldsymbol{z}}(1)$ to construct $\tilde{\boldsymbol{\Upsilon}}$, $\tilde{\boldsymbol{\Upsilon}}_{\tilde{\boldsymbol{z}}(1)}$ is as shown in (B.43) while the trailing columns in (B.42) corresponding to the Jacobian of $\tilde{\boldsymbol{\Upsilon}}$ with respect to the unknown constants in $\tilde{\boldsymbol{z}}(0)$ should be all zeros.

@fb@secFB

B.3 Final Details

In equations (B.1), (B.2), (B.3), (B.4), and (B.6), the second to last row is needed only if the initial time $a$ is free and the last row is needed only if the final time $b$ is free. In equation (B.5), the second to last row and column are needed only if the initial time $a$ is free and the last row and column are needed only if the final time $b$ is free.

In equations (B.25), (B.26), (B.27), (B.28), (B.31), (B.34), (B.37), (B.38), (B.39), (B.40), (B.41), and (B.44) the first row is needed only if the initial time $a$ is free and row $n+k_{1}+2$ is needed only if the final time $b$ is free. In equations (B.29), (B.30), (B.32), (B.33), (B.35), (B.36), (B.42), and (B.43) the first row and second to last column are needed only if the initial time $a$ is free and row $n+k_{1}+2$ and the last column are needed only if the final time $b$ is free.

In order to numerically solve the ODE TPBVP (A.19), (A.20), and (A.21) without continuation or with a monotonic continuation solver (such as c or cc), the solver should be provided (B.4), (B.5), (B.40), (B.42), and (B.43). In order to numerically solve the ODE TPBVP (A.19), (A.20), and (A.21) with a non-monotonic continuation solver (such as the predictor-corrector methods discussed in Appendices C and D), the solver should be provided (B.4), (B.5), (B.6), (B.40), (B.42), (B.43), and (B.44).

The first and second derivatives required to construct (B.4), (B.5), (B.6), (B.40), (B.42), (B.43), and (B.44) are generally quite tedious to derive manually. Instead, symbolic differentiation [48], complex/bicomplex step differentiation [49, 50, 51, 52], dual/hyper-dual numbers [53, 54, 55, 56], and automatic differentiation [57, 58] are computational alternatives. If fact, it may be shown that the use of dual/hyper-dual numbers to compute first and second derivatives is equivalent to automatic differentiation [56]. While symbolic differentiation suffers from expression explosion and complex/bicomplex step differentiation only applies to real analytic functions, dual/hyper-dual numbers and automatic differentiation are more robust and broadly-applicable.

Therefore, while (B.4), (B.5), (B.6), (B.40), (B.42), (B.43), and (B.44) are complicated, they may be readily constructed numerically through automatic differentiation of $H$, $\boldsymbol{\pi}$, $\hat{\mathbf{f}}$, $G$, $\boldsymbol{\sigma}$, and $\boldsymbol{\psi}$ if $\boldsymbol{\pi}$ is given analytically and of $H$, $\mathbf{f}$, $G$, $\boldsymbol{\sigma}$, and $\boldsymbol{\psi}$ if the value of $\boldsymbol{\pi}$ is constructed numerically via Newton’s method as in (A.11). There are many free automatic differentiation toolboxes available [31], such as the LAB automatic differentiation toolbox Gator [12, 13]. Moreover, Gator is able to construct vectorized automatic derivatives, which is extremely useful for realizing the vectorized version of (B.4), (B.5), and (B.6), as the non-vectorized version of these equations execute too slowly in LAB to solve the ODE TPBVP (A.19), (A.20), and (A.21) in a timely manner.

Appendix C Predictor-Corrector Continuation Method for Solving an ODE TPBVP

@fb@secFB

C.1 Introduction

Suppose it is desired to solve the ODE TPBVP:

[TABLE]

where $a,b\in\mathbb{R}$ are prescribed with $a<b$, $s\in\left[a,b\right]\subset\mathbb{R}$ is the independent variable, $n\in\mathbb{N}$ is the prescribed number of dependent variables in $\mathbf{y}$, $\mathbf{y}\colon\left[a,b\right]\to\mathbb{R}^{n}$ is an unknown function which must be solved for, $\lambda\in\mathbb{R}$ is a prescribed scalar parameter, $\mathbf{F}\colon\left[a,b\right]\times\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n}$ is a prescribed ODE velocity function defining the velocity of $\mathbf{y}$, and $\mathbf{G}\colon\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n}$ is a prescribed two-point boundary condition function. Observe that if $n=1$, $\mathbf{y}$, $\mathbf{F}$, and $\mathbf{G}$ are scalar-valued functions, while if $n>1$, $\mathbf{y}$, $\mathbf{F}$, and $\mathbf{G}$ are vector-valued functions. The Jacobian of $\mathbf{F}$ with respect to $\mathbf{y}$ is $\mathbf{F}_{\mathbf{y}}\colon\left[a,b\right]\times\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n\times n}$ and the Jacobian of $\mathbf{F}$ with respect to $\lambda$ is $\mathbf{F}_{\lambda}\colon\left[a,b\right]\times\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n\times 1}$. The Jacobian of $\mathbf{G}$ with respect to $\mathbf{y}(a)$ is $\mathbf{G}_{\mathbf{y}(a)}\colon\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n\times n}$, the Jacobian of $\mathbf{G}$ with respect to $\mathbf{y}(b)$ is $\mathbf{G}_{\mathbf{y}(b)}\colon\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n\times n}$, and the Jacobian of $\mathbf{G}$ with respect to $\lambda$ is $\mathbf{G}_{\lambda}\colon\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n\times 1}$. If $\mathbf{F}$ is linear in $\mathbf{y}$ and $\mathbf{G}$ is linear in $\mathbf{y}(a)$ and $\mathbf{y}(b)$, then (C.1) is said to be a linear ODE TPBVP; otherwise, (C.1) is said to be a nonlinear ODE TPBVP.

Note that a solution $\mathbf{y}$ to (C.1) depends on the given value of the scalar parameter $\lambda$, so a solution to (C.1) will be denoted by the pair $\left(\mathbf{y},\lambda\right)$. Usually it is not possible to solve (C.1) analytically. Instead, a numerical method such as a shooting, finite-difference, or Runge-Kutta method (collocation is a special kind of Runge-Kutta method) must be utilized to construct an approximate solution to (C.1). All such numerical methods require an initial solution guess and convergence to a solution is guaranteed only if the initial solution guess is sufficiently near the solution. Thus, solving (C.1) numerically requires construction of a good initial solution guess.

One way to construct a good initial solution guess for (C.1) is through continuation in the scalar parameter $\lambda$. If $\left(\mathbf{y}_{\mathrm{I}},\lambda_{\mathrm{I}}\right)$ solves (C.1) and it is desired to solve (C.1) for $\lambda=\lambda_{\mathrm{F}}$, it may be possible to construct a finite sequence of solutions $\left\{\left(\mathbf{y}_{j},\lambda_{j}\right)\right\}_{j=1}^{J+1}$ starting at the known solution $\left(\mathbf{y}_{1},\lambda_{1}\right)=\left(\mathbf{y}_{\mathrm{I}},\lambda_{\mathrm{I}}\right)$ and ending at the desired solution $\left(\mathbf{y}_{J+1},\lambda_{J+1}\right)=\left(\mathbf{y}_{\mathrm{F}},\lambda_{\mathrm{F}}\right)$, using the previous solution $\left(\mathbf{y}_{j},\lambda_{j}\right)$ as an initial solution guess for the numerical solver to obtain the next solution $\left(\mathbf{y}_{j+1},\lambda_{j+1}\right)$, $1\leq j\leq J$, in the sequence. $J\in\mathbb{N}$ denotes the number of solutions in the sequence after the known solution.

This appendix describes a particular such continuation method, called predictor-corrector continuation, for solving (C.1). The treatment given here follows [59]. In the literature, predictor-corrector continuation is also called path-following [60], predictor-corrector path-following [61], and differential path-following [62]. AUTO [63], COLCON [32], and the algorithm presented in [64] are Fortran predictor-corrrector continuation codes, while suite1.1 [65], bfun’s lowpath [59], and O [66] are LAB predictor-corrrector continuation codes. All these codes rely on global methods for solving ODE BVPs (e.g. Runge-Kutta, collocation, and finite-difference schemes), which are more robust than initial value methods for solving ODE BVPs (i.e. single and multiple shooting) because initial value methods cannot integrate unstable ODEs [36, 42, 43].

Before delving into the details, some functional analysis is reviewed which is necessary to understand how the predictor-corrector continuation method is applied to solve (C.1).

@fb@secFB

C.2 A Hilbert Space

Let $\mathcal{H}=\left\{\left(\mathbf{y},\lambda\right):\mathbf{y}\in L^{2}\left(\left[a,b\right],\mathbb{R}^{n}\right),\lambda\in\mathbb{R}\right\}$. $\mathcal{H}$ is a Hilbert space over $\mathbb{R}$. If $\alpha,\beta\in\mathbb{R}$ and $\left(\mathbf{y},\lambda\right),\left(\tilde{\mathbf{y}},\tilde{\lambda}\right)\in\mathcal{H}$, then

[TABLE]

the inner product on $\mathcal{H}$ is

[TABLE]

and the norm on $\mathcal{H}$, induced by the inner product, is

[TABLE]

$\left(\mathbf{y},\lambda\right)\in\mathcal{H}$ and $\left(\tilde{\mathbf{y}},\tilde{\lambda}\right)\in\mathcal{H}$ are said to be orthogonal iff $\left<\left(\mathbf{y},\lambda\right),\left(\tilde{\mathbf{y}},\tilde{\lambda}\right)\right>=0$, and $\left(\mathbf{y},\lambda\right)\in\mathcal{H}$ is said to be of unit length iff $\left\|\left(\mathbf{y},\lambda\right)\right\|=1$.

@fb@secFB

C.3 The Fréchet Derivative and Newton’s Method

Given a function $\mathbf{F}:\mathbb{R}^{n}\to\mathbb{R}^{m}$, recall that ordinary vector calculus defines the Jacobian of $\mathbf{F}$ as the function $\mathbf{F}^{\prime}:\mathbb{R}^{n}\to\mathbb{R}^{m\times n}$ such that $\mathbf{F}^{\prime}\left({\boldsymbol{x}}\right)$ is the linearization of $\mathbf{F}$ at ${\boldsymbol{x}}\in\mathbb{R}^{n}$. Given normed vector spaces $V$ and $W$ and an open subset $U$ of $V$, the Fréchet derivative is an extension of the Jacobian to an operator $\mathcal{F}:U\to W$. Before giving the definition of the Fréchet derivative, recall that $L\left(V,W\right)$ denotes the space of continuous linear operators from $V$ to $W$. Now for the definition of the Fréchet derivative, which comes from Definition 2.2.4 of [59].

Definition C.3.1

Suppose that $V$ and $W$ are normed vector spaces, and let $U$ be an open subset of $V$. Then the operator $\mathcal{F}:U\to W$ is said to be Fréchet differentiable at $u\in U$ if and only if there exists an operator $\mathcal{L}\in L\left(V,W\right)$ such that

[TABLE]

The operator $\mathcal{L}$ is then called the Fréchet derivative of $\mathcal{F}$ at $u$, often denoted by $\mathcal{F}^{\prime}(u)$. If $\mathcal{F}$ is Fréchet differentiable at all points in $U$, $\mathcal{F}$ is said to be Fréchet differentiable in $U$.

Given a function $\mathbf{H}:\mathbb{R}^{m}\to\mathbb{R}^{m}$, Newton’s method is an algorithm to solve $\mathbf{H}\left({\boldsymbol{x}}\right)=\mathbf{0}$ for ${\boldsymbol{x}}\in\mathbb{R}^{m}$ and $\mathbf{0}\in\mathbb{R}^{m}$ when $\mathbf{H}$ satisfies certain mild conditions. Starting from an initial solution guess ${\boldsymbol{x}}_{0}\in\mathbb{R}^{m}$ sufficiently close to a solution, Newton’s method converges to a solution of $\mathbf{H}\left({\boldsymbol{x}}\right)=\mathbf{0}$ by iteratively solving the equations

[TABLE]

starting at $k=0$, where $\mathbf{H}^{\prime}$ denotes the Jacobian of $\mathbf{H}$ and ${\boldsymbol{x}}_{k},\delta{\boldsymbol{x}}_{k}\in\mathbb{R}^{m}$ for $k\geq 0$. The iteration in (C.6) continues until $\mathbf{H}\left({\boldsymbol{x}}_{k}\right)\approx\mathbf{0}$ (or $\delta{\boldsymbol{x}}_{k}\approx\mathbf{0}$) or until $k$ exceeds a maximum iteration threshold.

Now consider an operator $\mathcal{H}:U\subset V\to W$, where $V$ and $W$ are Banach spaces and $U$ is an open subset of $V$. Kantorovich [67] provided an extension of Newton’s method to solve $\mathcal{H}(u)=0$ for $u\in U$ and $0\in W$ when $\mathcal{H}$ satisfies certain mild conditions. Starting from an initial solution guess $u_{0}\in U$ sufficiently close to a solution, Kantorovich’s extension of Newton’s method converges to a solution of $\mathcal{H}(u)=0$ by iteratively solving the equations

[TABLE]

starting at $k=0$, where $\mathcal{H}^{\prime}$ denotes the Fréchet derivative of $\mathcal{H}$ and $u_{k},\delta u_{k}\in U$ for $k\geq 0$. The iteration in (C.7) continues until $\mathcal{H}\left(u_{k}\right)\approx 0$ (or $\delta u_{k}\approx 0$) or until $k$ exceeds a maximum iteration threshold.

@fb@secFB

C.4 The Davidenko ODE IVP

To motivate the predictor-corrector continuation method, the Davidenko ODE IVP is first presented. Let $\mathcal{C}=\left\{\left(\mathbf{y},\lambda\right):\left(\mathbf{y},\lambda\right)\mathrm{\,solves\,}\eqref{eqn_ode_tpbvp}\right\}$ denote the solution manifold of (C.1). Suppose the solution manifold $\mathcal{C}$ is parameterized by arclength $\nu$, so that an element of $\mathcal{C}$ is $\left(\mathbf{y}(\nu),\lambda(\nu)\right)$, the tangent $\left(\mathbf{v}(\nu),\tau(\nu)\right)$ to $\mathcal{C}$ at $\left(\mathbf{y}(\nu),\lambda(\nu)\right)$ satisfies $\left\|\left(\mathbf{v}(\nu),\tau(\nu)\right)\right\|^{2}=\int_{a}^{b}\mathbf{v}^{\mathsf{T}}(s,\nu)\mathbf{v}(s,\nu)\mathrm{d}s+\left[\tau(\nu)\right]^{2}=1$ (i.e. $\left(\mathbf{v}(\nu),\tau(\nu)\right)$ is a unit tangent), and the solution manifold $\mathcal{C}$ can be described as a solution curve. With this arclength parameterization, $\mathbf{y}\colon\left[a,b\right]\times\mathbb{R}\to\mathbb{R}^{n}$, $\lambda\colon\mathbb{R}\to\mathbb{R}$, $\mathbf{v}\colon\left[a,b\right]\times\mathbb{R}\to\mathbb{R}^{n}$, $\tau\colon\mathbb{R}\to\mathbb{R}$, $\mathbf{y}(\nu)$ is shorthand for $\mathbf{y}(\cdot,\nu)\colon\left[a,b\right]\to\mathbb{R}^{n}$, and $\mathbf{v}(\nu)$ is shorthand for $\mathbf{v}(\cdot,\nu)\colon\left[a,b\right]\to\mathbb{R}^{n}$. Note that the components of the unit tangent $\left(\mathbf{v}(\nu),\tau(\nu)\right)$ to $\mathcal{C}$ at $\left(\mathbf{y}(\nu),\lambda(\nu)\right)$ are given explicitly by $\mathbf{v}(s,\nu)=\frac{\partial\mathbf{y}(s,\nu)}{\partial\nu}$ and $\tau(\nu)=\frac{\mathrm{d}\lambda(\nu)}{\mathrm{d}\nu}$.

The Fréchet derivative of the ODE TPBVP (C.1) with respect to $\nu$ about the solution $\left(\mathbf{y}(\nu),\lambda(\nu)\right)$, in conjunction with the arclength constraint and the initial condition $\left(\mathbf{y}_{\mathrm{I}},\lambda_{\mathrm{I}}\right)$, gives the nonlinear ODE IVP in the independent arclength variable $\nu$:

[TABLE]

which must be solved for $\left(\mathbf{y}(\nu),\lambda(\nu)\right)$ starting at $\nu_{0}$ from an initial solution $\left(\mathbf{y}_{\mathrm{I}},\lambda_{\mathrm{I}}\right)$ of (C.1). (C.8) is called the Davidenko ODE IVP and its solution is called the Davidenko flow [68, 69, 70, 33, 59, 71]. The first two equations in (C.8) constitute the Fréchet derivative of the ODE TPBVP (C.1), the third equation is the arclength constraint, and the final equation is the initial condition. By introducing a dummy scalar-valued function $w$ to represent the integrand of the arclength constraint, (C.8) can be re-written:

[TABLE]

Again, letting $\nu$ vary, (C.9) is a nonlinear ODE IVP which must be solved for $\left(\mathbf{y}(\nu),\lambda(\nu)\right)$ (i.e. $\mathbf{y}\colon\left[a,b\right]\times\mathbb{R}\to\mathbb{R}^{n}$ and $\lambda\colon\mathbb{R}\to\mathbb{R}$) starting at $\nu_{0}$ from an initial solution $\left(\mathbf{y}_{\mathrm{I}},\lambda_{\mathrm{I}}\right)$ of (C.1). However, for a fixed $\nu$, (C.9) is a nonlinear ODE TPBVP which must be solved for $\mathbf{v}(\cdot,\nu)\colon\left[a,b\right]\to\mathbb{R}^{n}$, $\tau(\nu)\in\mathbb{R}$, and $w(\cdot,\nu)\colon\left[a,b\right]\to\mathbb{R}$ and where the independent variable is $s\in\left[a,b\right]$.

As explained in Chapter 5 of [33], it is inadvisable to integrate the Davidenko ODE IVP (C.8), or equivalently (C.9). Instead, a predictor-corrector continuation method, depicted in Figure C.1 and explained in detail in the following subappendices, is used to generate a solution sequence $\left\{\left(\mathbf{y}_{j},\lambda_{j}\right)\right\}_{j=1}^{J}$ which is a discrete subset of the Davidenko flow such that $\left(\mathbf{y}_{1},\lambda_{1}\right)=\left(\mathbf{y}_{\mathrm{I}},\lambda_{\mathrm{I}}\right)$.

@fb@secFB

C.5 Construct the Tangent

Given a solution $\left(\mathbf{y}_{j},\lambda_{j}\right)$ to (C.1) and a unit tangent $\left(\mathbf{v}_{j-1},\tau_{j-1}\right)$ to the previous solution $\left(\mathbf{y}_{j-1},\lambda_{j-1}\right)$ to (C.1), we seek to construct a tangent $\left(\mathbf{v}_{j},\tau_{j}\right)$ to the solution curve $\mathcal{C}$ at $\left(\mathbf{y}_{j},\lambda_{j}\right)$ which is roughly of unit length. The arclength constraint is

[TABLE]

which is nonlinear in the tangent $\left(\mathbf{v}_{j},\tau_{j}\right)$. An alternative constraint, the pseudo-arclength constraint, is

[TABLE]

which, in contrast to the arclength constraint (C.10), is linear in the tangent $\left(\mathbf{v}_{j},\tau_{j}\right)$. The linearization (i.e. Fréchet derivative) of the ODE TPBVP (C.1) about the solution $\left(\mathbf{y}_{j},\lambda_{j}\right)$, in conjunction with the pseudo-arclength condition (C.11), gives the linear ODE TPBVP:

[TABLE]

which must be solved for $\mathbf{v}_{j}\colon\left[a,b\right]\to\mathbb{R}^{n}$, $\tau_{j}\in\mathbb{R}$, and $w\colon\left[a,b\right]\to\mathbb{R}$ and where $\left(\mathbf{v}_{j},\tau_{j}\right)$ is a tangent to $\mathcal{C}$ at $\left(\mathbf{y}_{j},\lambda_{j}\right)$. Note that the first, second, and third equations in (C.12) are the ODEs, while the fourth, fifth, and sixth equations constitute the boundary conditions. The first, second, and fourth equations in (C.12) are the linearization (i.e. Fréchet derivative) of (C.1) about the solution $\left(\mathbf{y}_{j},\lambda_{j}\right)$ and ensure that a tangent is produced, while the third, fifth, and sixth equations in (C.12) enforce the pseudo-arclength condition (C.11). The initial solution guess to solve (C.12) is $\left(\mathbf{v}_{j},\tau_{j}\right)=\left(\mathbf{v}_{j-1},\tau_{j-1}\right)$ and $w(s)=\int_{a}^{s}\mathbf{v}_{j-1}^{\mathsf{T}}(\tilde{s})\mathbf{v}_{j-1}(\tilde{s})\mathrm{d}\tilde{s}$, $s\in\left[a,b\right]$, for $j\geq 1$. For $j=1$, define $\left(\mathbf{v}_{0},\tau_{0}\right)=\left(\mathbf{0}_{n\times 1},1\right)$. Note that the construction of the initial guess for $w$ can be realized efficiently via the LAB routine trapz.

Note that the linear ODE TPBVP (C.12) can be solved numerically via the LAB routines p or twp, which offers 4 algorithms: bvp_m, bvpc_m, bvp_l, and bvpc_l; moreover, p and twp have special algorithms to solve linear ODE TPBVPs. Since $\mathbf{y}_{j}$ and $\mathbf{v}_{j-1}$ are usually only known at a discrete set of points in $\left[a,b\right]$, the values of these functions at the other points in $\left[a,b\right]$ must be obtained through interpolation in order to numerically solve (C.12). The LAB routine erp1 performs linear, cubic, pchip, makima, and spline interpolation and may be utilized to interpolate $\mathbf{y}_{j}$ and $\mathbf{v}_{j-1}$ while solving (C.12).

Because the numerical solvers usually converge faster when provided Jacobians of the ODE velocity function and of the two-point boundary condition function, these are computed below. Let

[TABLE]

The ODE velocity function in (C.12) is

[TABLE]

The Jacobian of the ODE velocity function $\mathbf{H}^{\mathrm{t}}$ with respect to $\mathbf{x}$ is

[TABLE]

The two-point boundary condition in (C.12) is

[TABLE]

where $\mathbf{K}^{\mathrm{t}}$ is the two-point boundary condition function

[TABLE]

The Jacobians of the two-point boundary condition function $\mathbf{K}^{\mathrm{t}}$ with respect to $\mathbf{x}(a)$ and $\mathbf{x}(b)$ are

[TABLE]

and

[TABLE]

Special care must be taken when implementing the Jacobians (C.18) and (C.19). Since the unknown constant $\tau_{j}$ appears as the second to last element of both $\mathbf{x}(a)$ and $\mathbf{x}(b)$, $\tau_{j}$ from only one of $\mathbf{x}(a)$ and $\mathbf{x}(b)$ is actually used to construct each term in $\mathbf{K}^{\mathrm{t}}$ involving $\tau_{j}$. The middle column of (C.18) is actually the derivative of $\mathbf{K}^{\mathrm{t}}$ with respect to the $\tau_{j}$ in $\mathbf{x}(a)$, while the middle column of (C.19) is actually the derivative of $\mathbf{K}^{\mathrm{t}}$ with respect to the $\tau_{j}$ in $\mathbf{x}(b)$. Thus, the middle columns in (C.18) and (C.19) corresponding to the derivative of $\mathbf{K}^{\mathrm{t}}$ with respect to $\tau_{j}$ should not coincide in a software implementation. For example, if $\mathbf{K}^{\mathrm{t}}$ is constructed from the $\tau_{j}$ in $\mathbf{x}(a)$, $\mathbf{K}_{\mathbf{x}(a)}^{\mathrm{t}}$ is as shown in (C.18) while the middle column of (C.19) corresponding to the derivative of $\mathbf{K}^{\mathrm{t}}$ with respect to the $\tau_{j}$ in $\mathbf{x}(b)$ is all zeros. Alternatively, if $\mathbf{K}^{\mathrm{t}}$ is constructed from the $\tau_{j}$ in $\mathbf{x}(b)$, $\mathbf{K}_{\mathbf{x}(b)}^{\mathrm{t}}$ is as shown in (C.19) while the middle column of (C.18) corresponding to the derivative of $\mathbf{K}^{\mathrm{t}}$ with respect to the $\tau_{j}$ appearing in $\mathbf{x}(a)$ is all zeros.

@fb@secFB

C.6 Normalize the Tangent

The tangent $\left(\mathbf{v}_{j},\tau_{j}\right)$ at $\left(\mathbf{y}_{j},\lambda_{j}\right)$ obtained by solving (C.12) in the previous step is only roughly of unit length. A unit tangent at $\left(\mathbf{y}_{j},\lambda_{j}\right)$ is obtained from $\left(\mathbf{v}_{j},\tau_{j}\right)$ through normalization:

[TABLE]

where

[TABLE]

The integration operator to construct the normalization scalar $\kappa$ in (C.21) can be realized via the LAB routine pz.

@fb@secFB

C.7 Construct the Tangent Predictor

The unit tangent $\left(\mathbf{v}_{j},\tau_{j}\right)$ constructed in (C.20) is used to obtain a guess (the so-called “tangent predictor”) $\left(\mathbf{y}_{1}^{\mathrm{c}},\lambda_{1}^{\mathrm{c}}\right)$ for the next solution $\left(\mathbf{y}_{j+1},\lambda_{j+1}\right)$ as follows:

[TABLE]

where $\sigma\in\left[\sigma_{\mathrm{min}},\sigma_{\mathrm{max}}\right]$ is a steplength and where $0<\sigma_{\mathrm{min}}\leq\sigma_{\mathrm{max}}$. Concretely, $\sigma_{\mathrm{min}}$ might be $.0001$ and $\sigma_{\mathrm{max}}$ might be $\frac{1}{2}$. $\sigma$ is adapted during the predictor-corrector continuation method based on the corrector step, discussed in the next subappendix. Initially, the value of $\sigma$ is set to $\sigma_{\mathrm{init}}\in\left[\sigma_{\mathrm{min}},\sigma_{\mathrm{max}}\right]$. The notation $\left(\mathbf{y}_{1}^{\mathrm{c}},\lambda_{1}^{\mathrm{c}}\right)$ is used to denote the tangent predictor in (C.22) because, as discussed in the next subappendix, the tangent predictor is used as the initial corrector in an iterative Newton’s method that projects the tangent predictor onto $\mathcal{C}$.

@fb@secFB

C.8 Construct the Corrector

Since the tangent predictor $\left(\mathbf{y}_{1}^{\mathrm{c}},\lambda_{1}^{\mathrm{c}}\right)$ constructed in (C.22) does not necessarily lie on $\mathcal{C}$, $\left(\mathbf{y}_{1}^{\mathrm{c}},\lambda_{1}^{\mathrm{c}}\right)$ must be projected onto $\mathcal{C}$ to obtain the next solution (the so-called “corrector”) $\left(\mathbf{y}_{j+1},\lambda_{j+1}\right)$. This projection process is the corrector step. In order to perform the projection efficiently, the difference between the next solution and the tangent predictor, $\left(\mathbf{y}_{j+1},\lambda_{j+1}\right)-\left(\mathbf{y}_{1}^{\mathrm{c}},\lambda_{1}^{\mathrm{c}}\right)$, should be orthogonal to the unit tangent $\left(\mathbf{v}_{j},\tau_{j}\right)$. That is, the orthogonality constraint is

[TABLE]

The tangent predictor $\left(\mathbf{y}_{1}^{\mathrm{c}},\lambda_{1}^{\mathrm{c}}\right)$ can be iteratively corrected by applying Newton’s method to (C.1), while enforcing the orthogonality constraint (C.23), to generate a sequence of correctors $\left\{\left(\mathbf{y}_{k}^{\mathrm{c}},\lambda_{k}^{\mathrm{c}}\right)\right\}_{k=1}^{K+1}$. Applying Newton’s method to the ODE TPBVP (C.1) about the current corrector $\left(\mathbf{y}_{k}^{\mathrm{c}},\lambda_{k}^{\mathrm{c}}\right)$, in conjunction with the orthogonality constraint (C.23), gives the linear ODE TPBVP:

[TABLE]

which must be solved for $\delta\mathbf{y}_{k}^{\mathrm{c}}\colon\left[a,b\right]\to\mathbb{R}^{n}$, $\delta\lambda_{k}^{\mathrm{c}}\in\mathbb{R}$, and $w\colon\left[a,b\right]\to\mathbb{R}$ and where $\left(\delta\mathbf{y}_{k}^{\mathrm{c}},\delta\lambda_{k}^{\mathrm{c}}\right)$ represents a correction to the current corrector $\left(\mathbf{y}_{k}^{\mathrm{c}},\lambda_{k}^{\mathrm{c}}\right)$. Note that the first, second, and third equations in (C.24) are the ODEs, while the fourth, fifth, and sixth equations constitute the boundary conditions. The first, second, and fourth equations in (C.24) are the result of applying Newton’s method to (C.1) about the current corrector $\left(\mathbf{y}_{k}^{\mathrm{c}},\lambda_{k}^{\mathrm{c}}\right)$, while the third, fifth, and sixth equations in (C.24) enforce the orthogonality constraint (C.23). (C.24) must be solved iteratively for at most $K$ iterations, so that $1\leq k\leq K$. The initial solution guess to solve (C.24) at the beginning of each iteration is $\left(\delta\mathbf{y}_{k}^{\mathrm{c}},\delta\lambda_{k}^{\mathrm{c}}\right)=\left(\mathbf{0}_{n\times 1},0\right)$ and $w(s)=0$, $s\in\left[a,b\right]$. The initial corrector about which Newton’s method is applied in the first iteration is the tangent predictor $\left(\mathbf{y}_{1}^{\mathrm{c}},\lambda_{1}^{\mathrm{c}}\right)$. At the end of each iteration, the corrector about which Newton’s method is applied for the next iteration is updated via $\left(\mathbf{y}_{k+1}^{\mathrm{c}},\lambda_{k+1}^{\mathrm{c}}\right)=\left(\mathbf{y}_{k}^{\mathrm{c}},\lambda_{k}^{\mathrm{c}}\right)+\left(\delta\mathbf{y}_{k}^{\mathrm{c}},\delta\lambda_{k}^{\mathrm{c}}\right)$. At the end of each iteration, convergence to $\mathcal{C}$ should be tested via:

[TABLE]

where $\gamma$ is a small threshold such as $.001$. Since Newton’s method enjoys quadratic convergence near a solution, only a few (say $K=5$) iterative solves of (C.24) should be attempted. If convergence has not been attained in $K$ iterations, the steplength $\sigma$ should be reduced:

[TABLE]

where $\sigma_{\mathrm{r}}$ is a reduction scale factor such as $\frac{1}{4}$ and the corrector step should be restarted at the new tangent predictor $\left(\mathbf{y}_{1}^{\mathrm{c}},\lambda_{1}^{\mathrm{c}}\right)=\left(\mathbf{y}_{j},\lambda_{j}\right)+\sigma\left(\mathbf{v}_{j},\tau_{j}\right)$, based on the updated value of $\sigma$ realized in (C.26). If, as a result of the reduction realized in (C.26), $\sigma<\sigma_{\mathrm{min}}$, the algorithm should halt and predictor-corrector continuation failed. However, if convergence has been achieved in $k\leq K$ iterations, the next solution can be taken to be $\left(\mathbf{y}_{j+1},\lambda_{j+1}\right)=\left(\mathbf{y}_{k+1}^{\mathrm{c}},\lambda_{k+1}^{\mathrm{c}}\right)$ or the corrector can be further polished as explained in the next subappendix. Moreover, if convergence has been achieved rapidly in no more than $k_{\mathrm{fast}}$ iterations, where $1\leq k_{\mathrm{fast}}\leq K$ and, concretely, $k_{\mathrm{fast}}$ might be 3, then the steplength $\sigma$ may be increased:

[TABLE]

where $\sigma_{\mathrm{i}}$ is an increase scale factor such as $2$.

Note that the linear ODE TPBVP (C.24) can be solved numerically via the LAB routines p or twp, which offers 4 algorithms: bvp_m, bvpc_m, bvp_l, and bvpc_l; moreover, p and twp have special algorithms to solve linear ODE TPBVPs. Since $\mathbf{y}_{k}^{\mathrm{c}}$, $\frac{\mathrm{d}}{\mathrm{d}s}\mathbf{y}_{k}^{\mathrm{c}}$, and ${\mathbf{v}}_{j}$ are usually only known at a discrete set of points in $\left[a,b\right]$, the values of these functions at the other points in $\left[a,b\right]$ must be obtained through interpolation in order to numerically solve (C.24). The LAB routine erp1 performs linear, cubic, pchip, makima, and spline interpolation and may be utilized to interpolate $\mathbf{y}_{k}^{\mathrm{c}}$, $\frac{\mathrm{d}}{\mathrm{d}s}\mathbf{y}_{k}^{\mathrm{c}}$, and ${\mathbf{v}}_{j}$ while solving (C.24).

Because the numerical solvers usually converge faster when provided Jacobians of the ODE velocity function and of the two-point boundary condition function, these are computed below. Let

[TABLE]

The ODE velocity function in (C.24) is

[TABLE]

The Jacobian of the ODE velocity function $\mathbf{H}^{\mathrm{c}}$ with respect to $\mathbf{x}$ is

[TABLE]

The two-point boundary condition in (C.24) is

[TABLE]

where $\mathbf{K}^{\mathrm{c}}$ is the two-point boundary condition function

[TABLE]

The Jacobians of the two-point boundary condition function $\mathbf{K}^{\mathrm{c}}$ with respect to $\mathbf{x}(a)$ and $\mathbf{x}(b)$ are

[TABLE]

and

[TABLE]

Special care must be taken when implementing the Jacobians (C.33) and (C.34). Since the unknown constant $\delta\lambda_{k}^{\mathrm{c}}$ appears as the second to last element of both $\mathbf{x}(a)$ and $\mathbf{x}(b)$, $\delta\lambda_{k}^{\mathrm{c}}$ from only one of $\mathbf{x}(a)$ and $\mathbf{x}(b)$ is actually used to construct each term in $\mathbf{K}^{\mathrm{c}}$ involving $\delta\lambda_{k}^{\mathrm{c}}$. The middle column of (C.33) is actually the derivative of $\mathbf{K}^{\mathrm{c}}$ with respect to the $\delta\lambda_{k}^{\mathrm{c}}$ in $\mathbf{x}(a)$, while the middle column of (C.34) is actually the derivative of $\mathbf{K}^{\mathrm{c}}$ with respect to the $\delta\lambda_{k}^{\mathrm{c}}$ in $\mathbf{x}(b)$. Thus, the middle columns in (C.33) and (C.34) corresponding to the derivative of $\mathbf{K}^{\mathrm{c}}$ with respect to $\delta\lambda_{k}^{\mathrm{c}}$ should not coincide in a software implementation. For example, if $\mathbf{K}^{\mathrm{c}}$ is constructed from the $\delta\lambda_{k}^{\mathrm{c}}$ in $\mathbf{x}(a)$, $\mathbf{K}_{\mathbf{x}(a)}^{\mathrm{c}}$ is as shown in (C.33) while the middle column of (C.34) corresponding to the derivative of $\mathbf{K}^{\mathrm{c}}$ with respect to the $\delta\lambda_{k}^{\mathrm{c}}$ in $\mathbf{x}(b)$ is all zeros. Alternatively, if $\mathbf{K}^{\mathrm{c}}$ is constructed from the $\delta\lambda_{k}^{\mathrm{c}}$ in $\mathbf{x}(b)$, $\mathbf{K}_{\mathbf{x}(b)}^{\mathrm{c}}$ is as shown in (C.34) while the middle column of (C.33) corresponding to the derivative of $\mathbf{K}^{\mathrm{c}}$ with respect to the $\delta\lambda_{k}^{\mathrm{c}}$ appearing in $\mathbf{x}(a)$ is all zeros.

@fb@secFB

C.9 Polish the Corrector

The final corrector $\left(\mathbf{y}_{k+1}^{\mathrm{c}},\lambda_{k+1}^{\mathrm{c}}\right)$ from the previous step can be further polished by finding $\left(\mathbf{y}_{j+1},\lambda_{j+1}\right)$ that solves (C.1) while satisfying the orthogonality constraint (C.23). This yields the ODE TPBVP:

[TABLE]

which must be solved for $\mathbf{y}_{j+1}\colon\left[a,b\right]\to\mathbb{R}^{n}$, $\lambda_{j+1}\in\mathbb{R}$, and $w\colon\left[a,b\right]\to\mathbb{R}$. Note that the first, second, and third equations in (C.35) are the ODEs, while the fourth, fifth, and sixth equations constitute the boundary conditions. The first, second, and fourth equations in (C.35) ensure that the solution lies on $\mathcal{C}$ (i.e. satisfies (C.1)), while the third, fifth, and sixth equations in (C.35) enforce the orthogonality constraint (C.23). The initial solution guess to solve (C.35) is the final corrector $\left(\mathbf{y}_{k+1}^{\mathrm{c}},\lambda_{k+1}^{\mathrm{c}}\right)$ from the previous step and $w(s)=0$, $s\in\left[a,b\right]$.

Note that the ODE TPBVP (C.35) can be solved numerically via the LAB routines p or twp, which offers 4 algorithms: bvp_m, bvpc_m, bvp_l, and bvpc_l. Since $\mathbf{y}_{1}^{\mathrm{c}}$ and ${\mathbf{v}}_{j}$ are usually only known at a discrete set of points in $\left[a,b\right]$, the values of these functions at the other points in $\left[a,b\right]$ must be obtained through interpolation in order to numerically solve (C.35). The LAB routine erp1 performs linear, cubic, pchip, makima, and spline interpolation and may be utilized to interpolate $\mathbf{y}_{1}^{\mathrm{c}}$ and ${\mathbf{v}}_{j}$ while solving (C.35).

Because the numerical solvers usually converge faster when provided Jacobians of the ODE velocity function and of the two-point boundary condition function, these are computed below. Let

[TABLE]

The ODE velocity function in (C.35) is

[TABLE]

The Jacobian of the ODE velocity function $\mathbf{H}^{\mathrm{p}}$ with respect to $\mathbf{x}$ is

[TABLE]

The two-point boundary condition in (C.35) is

[TABLE]

where $\mathbf{K}^{\mathrm{p}}$ is the two-point boundary condition function

[TABLE]

The Jacobians of the two-point boundary condition function $\mathbf{K}^{\mathrm{p}}$ with respect to $\mathbf{x}(a)$ and $\mathbf{x}(b)$ are

[TABLE]

and

[TABLE]

Special care must be taken when implementing the Jacobians (C.41) and (C.42). Since the unknown constant $\lambda_{j+1}$ appears as the second to last element of both $\mathbf{x}(a)$ and $\mathbf{x}(b)$, $\lambda_{j+1}$ from only one of $\mathbf{x}(a)$ and $\mathbf{x}(b)$ is actually used to construct each term in $\mathbf{K}^{\mathrm{p}}$ involving $\lambda_{j+1}$. The middle column of (C.41) is actually the derivative of $\mathbf{K}^{\mathrm{p}}$ with respect to the $\lambda_{j+1}$ in $\mathbf{x}(a)$, while the middle column of (C.42) is actually the derivative of $\mathbf{K}^{\mathrm{p}}$ with respect to the $\lambda_{j+1}$ in $\mathbf{x}(b)$. Thus, the middle columns in (C.41) and (C.42) corresponding to the derivative of $\mathbf{K}^{\mathrm{p}}$ with respect to $\lambda_{j+1}$ should not coincide in a software implementation. For example, if $\mathbf{K}^{\mathrm{p}}$ is constructed from the $\lambda_{j+1}$ in $\mathbf{x}(a)$, $\mathbf{K}_{\mathbf{x}(a)}^{\mathrm{p}}$ is as shown in (C.41) while the middle column of (C.42) corresponding to the derivative of $\mathbf{K}^{\mathrm{p}}$ with respect to the $\lambda_{j+1}$ in $\mathbf{x}(b)$ is all zeros. Alternatively, if $\mathbf{K}^{\mathrm{p}}$ is constructed from the $\lambda_{j+1}$ in $\mathbf{x}(b)$, $\mathbf{K}_{\mathbf{x}(b)}^{\mathrm{p}}$ is as shown in (C.42) while the middle column of (C.41) corresponding to the derivative of $\mathbf{K}^{\mathrm{p}}$ with respect to the $\lambda_{j+1}$ appearing in $\mathbf{x}(a)$ is all zeros.

@fb@secFB

C.10 Pseudocode for Predictor-Corrector Continuation

Below is pseudocode that realizes the predictor-corrector continuation method.

Appendix D Sweep Predictor-Corrector Continuation Method for Solving an ODE TPBVP

@fb@secFB

D.1 Introduction

In this appendix, an alternative predictor-corrector continuation method is presented that exploits a monotonic continuation ODE TPBVP solver, such as twp’s c or cc, to monotonically increase (i.e. sweep) the tangent steplength $\sigma$ from [math] up until a maximum threshold $\sigma_{\mathrm{max}}$ is reached or until the next turning point is reached.

@fb@secFB

D.2 Construct the Tangent

Given a solution $\left(\mathbf{y}_{j},\lambda_{j}\right)$ to (C.1), we seek to construct a unit tangent $\left(\mathbf{v}_{j},\tau_{j}\right)$ to the solution curve $\mathcal{C}$ at $\left(\mathbf{y}_{j},\lambda_{j}\right)$. Recall the arclength constraint

[TABLE]

The linearization (i.e. Fréchet derivative) of the ODE TPBVP (C.1) about the solution $\left(\mathbf{y}_{j},\lambda_{j}\right)$, in conjunction with the arclength constraint (D.1), gives the nonlinear ODE TPBVP:

[TABLE]

which must be solved for $\mathbf{v}_{j}\colon\left[a,b\right]\to\mathbb{R}^{n}$, $\tau_{j}\in\mathbb{R}$, and $w\colon\left[a,b\right]\to\mathbb{R}$ and where $\left(\mathbf{v}_{j},\tau_{j}\right)$ is a unit tangent to $\mathcal{C}$ at $\left(\mathbf{y}_{j},\lambda_{j}\right)$. Note that the first, second, and third equations in (D.2) are the ODEs, while the fourth, fifth, and sixth equations constitute the boundary conditions. The first, second, and fourth equations in (D.2) are the linearization (i.e. Fréchet derivative) of (C.1) about the solution $\left(\mathbf{y}_{j},\lambda_{j}\right)$ and ensure that a tangent is produced, while the third, fifth, and sixth equations in (D.2) enforce the arclength constraint (D.1) ensuring that the tangent is of unit length. The initial solution guess to solve (D.2) is $\left(\mathbf{v}_{j},\tau_{j}\right)=\left(\mathbf{0}_{n\times 1},1\right)$ and $w(s)=0$, $s\in\left[a,b\right]$.

Note that the ODE TPBVP (D.2) can be solved numerically via the LAB routines p or twp, which offers 4 algorithms: bvp_m, bvpc_m, bvp_l, and bvpc_l. Since $\mathbf{y}_{j}$ is usually only known at a discrete set of points in $\left[a,b\right]$, the values of this function at the other points in $\left[a,b\right]$ must be obtained through interpolation in order to numerically solve (D.2). The LAB routine erp1 performs linear, cubic, pchip, makima, and spline interpolation and may be utilized to interpolate $\mathbf{y}_{j}$ while solving (D.2).

Because the numerical solvers usually converge faster when provided Jacobians of the ODE velocity function and of the two-point boundary condition function, these are computed below. Let

[TABLE]

The ODE velocity function in (D.2) is

[TABLE]

The Jacobian of the ODE velocity function $\mathbf{H}^{\mathrm{t}}$ with respect to $\mathbf{x}$ is

[TABLE]

The two-point boundary condition in (D.2) is

[TABLE]

where $\mathbf{K}^{\mathrm{t}}$ is the two-point boundary condition function

[TABLE]

The Jacobians of the two-point boundary condition function $\mathbf{K}^{\mathrm{t}}$ with respect to $\mathbf{x}(a)$ and $\mathbf{x}(b)$ are

[TABLE]

and

[TABLE]

Special care must be taken when implementing the Jacobians (D.8) and (D.9). Since the unknown constant $\tau_{j}$ appears as the second to last element of both $\mathbf{x}(a)$ and $\mathbf{x}(b)$, $\tau_{j}$ from only one of $\mathbf{x}(a)$ and $\mathbf{x}(b)$ is actually used to construct each term in $\mathbf{K}^{\mathrm{t}}$ involving $\tau_{j}$. The middle column of (D.8) is actually the derivative of $\mathbf{K}^{\mathrm{t}}$ with respect to the $\tau_{j}$ in $\mathbf{x}(a)$, while the middle column of (D.9) is actually the derivative of $\mathbf{K}^{\mathrm{t}}$ with respect to the $\tau_{j}$ in $\mathbf{x}(b)$. Thus, the middle columns in (D.8) and (D.9) corresponding to the derivative of $\mathbf{K}^{\mathrm{t}}$ with respect to $\tau_{j}$ should not coincide in a software implementation. For example, if $\mathbf{K}^{\mathrm{t}}$ is constructed from the $\tau_{j}$ in $\mathbf{x}(a)$, $\mathbf{K}_{\mathbf{x}(a)}^{\mathrm{t}}$ is as shown in (D.8) while the middle column of (D.9) corresponding to the derivative of $\mathbf{K}^{\mathrm{t}}$ with respect to the $\tau_{j}$ in $\mathbf{x}(b)$ is all zeros. Alternatively, if $\mathbf{K}^{\mathrm{t}}$ is constructed from the $\tau_{j}$ in $\mathbf{x}(b)$, $\mathbf{K}_{\mathbf{x}(b)}^{\mathrm{t}}$ is as shown in (D.9) while the middle column of (D.8) corresponding to the derivative of $\mathbf{K}^{\mathrm{t}}$ with respect to the $\tau_{j}$ appearing in $\mathbf{x}(a)$ is all zeros.

@fb@secFB

D.3 Determine the Tangent Direction

The unit tangent $\left(\mathbf{v}_{j},\tau_{j}\right)$ at $\left(\mathbf{y}_{j},\lambda_{j}\right)$ obtained by solving (D.2) must be scaled so that the sweep predictor-corrector continuation method does not reverse direction. As shown in [60], the correct direction for the unit tangent is obtained via:

[TABLE]

where $\kappa$ is the inner product of the previous and current unit tangents:

[TABLE]

The integration operator to construct the inner product $\kappa$ in (D.11) can be realized via the LAB routine pz. With the sign direction selected by (D.10), the inner product of the previous and current unit tangents is positive:

[TABLE]

@fb@secFB

D.4 Sweep along the Tangent

By monotonically increasing (or sweeping) the tangent steplength $\sigma$ from [math], the current solution $\left(\mathbf{y}_{j},\lambda_{j}\right)$ and its unit tangent $\left(\mathbf{v}_{j},\tau_{j}\right)$ can be used to find the next solution $\left(\mathbf{y}_{j+1},\lambda_{j+1}\right)$ that solves (C.1) while satisfying the orthogonality constraint:

[TABLE]

This yields the ODE TPBVP:

[TABLE]

which must be solved for $\mathbf{y}_{j+1}\colon\left[a,b\right]\to\mathbb{R}^{n}$, $\lambda_{j+1}\in\mathbb{R}$, and $w\colon\left[a,b\right]\to\mathbb{R}$ by monotonically increasing (or sweeping) $\sigma$. Note that the first, second, and third equations in (D.14) are the ODEs, while the fourth, fifth, and sixth equations constitute the boundary conditions. The first, second, and fourth equations in (D.14) ensure that the solution lies on $\mathcal{C}$ (i.e. satisfies (C.1)), while the third, fifth, and sixth equations in (D.14) enforce the orthogonality constraint (D.13). The initial solution guess to solve (D.14) is the current solution $\left(\mathbf{y}_{j},\lambda_{j}\right)$ and $w(s)=0$, $s\in\left[a,b\right]$. $\sigma$ starts at [math], since the initial solution guess for $\left(\mathbf{y}_{j+1},\lambda_{j+1}\right)$ is $\left(\mathbf{y}_{j},\lambda_{j}\right)$, and increases monotonically until the maximum threshold $\sigma_{\mathrm{max}}$ is reached or until the ODE TPBVP solver halts (due to reaching a turning point).

Note that the ODE TPBVP (D.14) can be solved numerically via the LAB routine twp, which offers 2 continuation algorithms: c and cc. The continuation algorithms c and cc assume that the continuation parameter (in this case $\sigma$) is monotonically increasing or decreasing, so that they will halt at a turning point in the continuation parameter. Since $\mathbf{y}_{j}$ and ${\mathbf{v}}_{j}$ are usually only known at a discrete set of points in $\left[a,b\right]$, the values of these functions at the other points in $\left[a,b\right]$ must be obtained through interpolation in order to numerically solve (D.14). The LAB routine erp1 performs linear, cubic, pchip, makima, and spline interpolation and may be utilized to interpolate $\mathbf{y}_{j}$ and ${\mathbf{v}}_{j}$ while solving (D.14).

Because the numerical solvers usually converge faster when provided Jacobians of the ODE velocity function and of the two-point boundary condition function, these are computed below. Let

[TABLE]

The ODE velocity function in (D.14) is

[TABLE]

The Jacobian of the ODE velocity function $\mathbf{H}^{\mathrm{q}}$ with respect to $\mathbf{x}$ is

[TABLE]

The two-point boundary condition in (D.14) is

[TABLE]

where $\mathbf{K}^{\mathrm{q}}$ is the two-point boundary condition function

[TABLE]

The Jacobians of the two-point boundary condition function $\mathbf{K}^{\mathrm{q}}$ with respect to $\mathbf{x}(a)$ and $\mathbf{x}(b)$ are

[TABLE]

and

[TABLE]

Special care must be taken when implementing the Jacobians (D.20) and (D.21). Since the unknown constant $\lambda_{j+1}$ appears as the second to last element of both $\mathbf{x}(a)$ and $\mathbf{x}(b)$, $\lambda_{j+1}$ from only one of $\mathbf{x}(a)$ and $\mathbf{x}(b)$ is actually used to construct each term in $\mathbf{K}^{\mathrm{q}}$ involving $\lambda_{j+1}$. The middle column of (D.20) is actually the derivative of $\mathbf{K}^{\mathrm{q}}$ with respect to the $\lambda_{j+1}$ in $\mathbf{x}(a)$, while the middle column of (D.21) is actually the derivative of $\mathbf{K}^{\mathrm{q}}$ with respect to the $\lambda_{j+1}$ in $\mathbf{x}(b)$. Thus, the middle columns in (D.20) and (D.21) corresponding to the derivative of $\mathbf{K}^{\mathrm{q}}$ with respect to $\lambda_{j+1}$ should not coincide in a software implementation. For example, if $\mathbf{K}^{\mathrm{q}}$ is constructed from the $\lambda_{j+1}$ in $\mathbf{x}(a)$, $\mathbf{K}_{\mathbf{x}(a)}^{\mathrm{q}}$ is as shown in (D.20) while the middle column of (D.21) corresponding to the derivative of $\mathbf{K}^{\mathrm{q}}$ with respect to the $\lambda_{j+1}$ in $\mathbf{x}(b)$ is all zeros. Alternatively, if $\mathbf{K}^{\mathrm{q}}$ is constructed from the $\lambda_{j+1}$ in $\mathbf{x}(b)$, $\mathbf{K}_{\mathbf{x}(b)}^{\mathrm{q}}$ is as shown in (D.21) while the middle column of (D.20) corresponding to the derivative of $\mathbf{K}^{\mathrm{q}}$ with respect to the $\lambda_{j+1}$ appearing in $\mathbf{x}(a)$ is all zeros.

@fb@secFB

D.5 Pseudocode for Sweep Predictor-Corrector Continuation

Below is pseudocode that realizes the sweep predictor-corrector continuation method.