{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Human Standing Control Parameter Identification with Direct Collocation\n",
"\n",
"Jason K. Moore \n",
"Ton van den Bogert \n",
"\n",
"\n",
"\n",
"July 10, 2015 \n",
"ISB TGCS 2015, Edinburgh, UK"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Parameter Identification\n",
"\n",
"Find the parameters $\\mathbf{p}$ such that the difference between the model simulation, $\\mathbf{y}$, and measurements, $\\mathbf{y}_m$ is minimized.\n",
"\n",
"### Dynamic system\n",
"\n",
"- Equations of Motion: $\\dot{\\mathbf{x}} = \\mathbf{f}(\\mathbf{x}, \\mathbf{p})$\n",
"- Measurement variables: $\\mathbf{y} = \\mathbf{g}(\\mathbf{x}, \\mathbf{p})$\n",
"\n",
"### Objective\n",
"\n",
"$$\\min\\limits_\\mathbf{p} J(\\mathbf{p})$$\n",
"where\n",
"$$J(\\mathbf{p}) = \\int [\\mathbf{y}_m - \\mathbf{y}(\\mathbf{p})]^2 dt$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Parameter Identification By Shooting\n",
"\n",
"- Repeated simulations are computationally costly\n",
"- Systems may be unstable and thus have an ill-defined objective\n",
"- Local minima are inevitable\n",
" - May requires a superb guess\n",
" - May need time intensive global optimization methods"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Local Minima Example: Simple pendulum\n",
"\n",
"> Vyasarayani, Chandrika P., Thomas Uchida, Ashwin Carvalho, and John McPhee.\n",
"\"Parameter Identification in Dynamic Systems Using the Homotopy Optimization\n",
"Approach\". Multibody System Dynamics 26, no. 4 (2011): 411-24.\n",
"\n",
"\n",
"\n",
"## 1-DoF, 1 parameter pendulum equations of motion\n",
"\n",
"$\\dot{\\mathbf{x}} = \\begin{bmatrix} \\dot{\\theta}(t) \\\\ \\dot{\\omega}(t) \\end{bmatrix} = \\begin{bmatrix} \\omega(t) \\\\ -p \\sin{\\theta}(t) \\end{bmatrix}$\n",
"\n",
"## Objective: Minimize least squares\n",
"\n",
"$J(p) = \\min\\limits_{p} \\int_{t_0}^{t_f} [\\theta_m(t) - \\theta(\\mathbf{x}, p, t)]^2 dt$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Direct Collocation\n",
"\n",
"> Betts, J., and W. Huffman. “Large Scale Parameter Estimation Using Sparse Nonlinear Programming Methods.” SIAM Journal on Optimization 14, no. 1 (January 1, 2003): 223–44. doi:10.1137/S1052623401399216.\n",
"\n",
"## Benefits\n",
"\n",
"- Fast computation times\n",
"- Handles unstable systems with ease\n",
"- Less susceptible to local minima\n",
"\n",
"## Disadvantages\n",
"\n",
"- Accurate solution requires large number of nodes\n",
"- Memory management for large sparse matrices and operations\n",
"- Tedious and error prone to form gradients, Jacobians, and Hessians"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Our Direct Collocation Implementation\n",
"\n",
"### Implicit Continous Equations of Motion\n",
"\n",
"No need to solve for $\\dot{\\mathbf{x}}$.\n",
"\n",
"$$\\mathbf{f}(\\dot{\\mathbf{x}}, \\mathbf{x}, \\mathbf{r}, \\mathbf{p}, t) = 0$$\n",
"\n",
"- $\\mathbf{x}, \\dot{\\mathbf{x}}$: states and their derivatives\n",
"- $\\mathbf{r}$: exogenous inputs\n",
"- $\\mathbf{p}$: constant parameters"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Discretization\n",
"\n",
"### First order discrete integration options:\n",
"\n",
"Backward Euler:\n",
"\n",
"$\\mathbf{f}(\\frac{\\mathbf{x}_i - \\mathbf{x}_{i-1}}{h} , \\mathbf{x}_i, \\mathbf{r}_i, \\mathbf{p}, t_i) = 0$\n",
"\n",
"Midpoint Rule:\n",
"\n",
"$\\mathbf{f}(\\frac{\\mathbf{x}_{i + 1} - \\mathbf{x}_i}{h}, \\frac{\\mathbf{x}_i + \\mathbf{x}_{i + 1}}{2}, \\frac{\\mathbf{r}_i + \\mathbf{r}_{i + 1}}{2}, \\mathbf{p}, t_i) = 0$\n",
"\n",
"### Nonlinear Programming Formulation\n",
" \n",
"$$\\min\\limits_{\\theta} J(\\theta)$$\n",
"\n",
"$$\\textrm{where } \\mathbf{\\theta} = [\\mathbf{x}_i, \\dots, \\mathbf{x}_N, \\mathbf{r}_i, \\ldots, \\mathbf{r}_N, \\mathbf{p}]$$\n",
"\n",
"$$\\textrm{subject to } \\mathbf{f}(\\mathbf{x}_i, \\mathbf{x}_{i+1}, \\mathbf{r}_i, \\mathbf{r}_{i + 1}, \\mathbf{p}, t_i) = 0 \\textrm{ and } \\theta_L \\leq \\theta \\leq \\theta_U$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Software Tool: opty\n",
"\n",
"- User specifies continous symbolic:\n",
" - objective\n",
" - equations of motion (explicit or implicit)\n",
" - additional constraints\n",
" - bounds on free variables: $\\mathbf{x}, \\mathbf{r}, \\mathbf{p}$\n",
"- EoMs can be generated with PyDy (http://pydy.org)\n",
"- Effficient just-in-time compiled C code is generated for functions that evaluate:\n",
" - objective and its gradient\n",
" - constraints and its Jacobian\n",
"- NLP problem automatically formed for IPOPT\n",
"- Open source: BSD license\n",
"- Written in Python\n",
"- http://github.com/csu-hmc/opty\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Example Code: 1 DoF, 1 Parameter Pendulum\n",
"\n",
"
\n",
"\n",
"$$\\dot{\\mathbf{x}} = \\begin{bmatrix} \\dot{\\theta}(t) \\\\ \\dot{\\omega}(t) \\end{bmatrix} = \\begin{bmatrix} \\omega(t) \\\\ -p \\sin{\\theta}(t) \\end{bmatrix}$$\n",
"\n",
"Paraphrased from https://github.com/csu-hmc/opty/blob/master/examples/vyasarayani2011.py\n",
"\n",
"### Symbolics\n",
"\n",
"```python\n",
"# Specify symbols for the parameters\n",
"p, t = symbols('p, t')\n",
"\n",
"# Specify the functions of time\n",
"theta, omega, theta_m = symbols('theta, omega, theta_m', cls=Function)\n",
"\n",
"# Specify the symbolic equations of motion\n",
"eom = (Matrix([theta(t), omega(t)]).diff(t) - \n",
" Matrix([omega(t), -p * sin(theta(t))]))\n",
"\n",
"# Specify the symbolic objective function\n",
"obj = Integral((theta_m(t) - theta(t))**2, t)\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Example Code: 1 DoF, 1 Parameter Pendulum\n",
"\n",
"### Numerics\n",
"\n",
"```python\n",
"# Choose discritzation values\n",
"num_nodes = 1000\n",
"interval = 0.01 # seconds\n",
"\n",
"# Form the problem\n",
"prob = Problem(obj, eom, (theta(t), omega(t)), num_nodes, interval,\n",
" known_trajectory_map={y1_m(t): measured_data},\n",
" integration_method='midpoint')\n",
"\n",
"# Set an initial guess\n",
"initial_guess = random(prob.num_free)\n",
"\n",
"# Solve the system\n",
"solution, info = prob.solve(initial_guess)\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Computational Speed\n",
"\n",
"### Example Larger System\n",
"\n",
"
$$\\dot{\\mathbf{x}} = \\mathbf{f}_o(\\mathbf{x}, \\mathbf{r}_c, \\mathbf{r}_k, \\mathbf{p}_k, t)$$
\n", "\n", " | Known | \n", "Identified | \n", "Error | \n", "
---|---|---|---|
$k_{00}$ | \n", "950 | \n", "946 | \n", "-0.4% | \n", "
$k_{01}$ | \n", "175 | \n", "177 | \n", "1.4% | \n", "
$k_{02}$ | \n", "185 | \n", "185 | \n", "-0.2% | \n", "
$k_{03}$ | \n", "50 | \n", "55 | \n", "9.4% | \n", "
$k_{10}$ | \n", "45 | \n", "45 | \n", "1.1% | \n", "
$k_{11}$ | \n", "290 | \n", "289 | \n", "-0.3% | \n", "
$k_{12}$ | \n", "60 | \n", "59 | \n", "-2.1% | \n", "
$k_{13}$ | \n", "26 | \n", "27 | \n", "4.2% | \n", "