hp-psuedospectral method¶
To solve the integral constraints within the optimal control problem we employs the hp-pseudospectral method. The hp-psuedospectral method is an form of Gaussian Quadrature, which uses multi-interval collocation points.
Single Phase Optimal Control¶
Find:
- The state: \(\mathbf{x}(t)\)
- The control: \(\mathbf{u}(t)\)
- The integrals: \(\mathbf{q}\)
- The initial time: \(t_0\)
- The final time: \(t_f\)
- To Minimize:
- \[J = \Phi(\mathbf{x}(t_0),\mathbf{x}(t_f),\mathbf{q},t_0,t_f)\]
That Satisfy the Following Constraints:
- Dynamic Constraints:
\[\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t} = \mathbf{\psi}(\mathbf{x}(t),\mathbf{u}(t),t)\]
- Inequality Path Constraints:
\[\mathbf{c}_{min} <= \mathbf{c}(\mathbf{x}(t),\mathbf{u}(t),t) <= \mathbf{c}_{max}\]
- Integral Constraints:
\[q_i = \int_{t_0}^{t_f} \Upsilon_i(\mathbf{x}(t),\mathbf{u}(t),t)\, \mathrm{d}t,\;\;(i=1,....,n_q)\]
- Event Constraints:
\[\mathbf{b}_{min} <= \mathbf{b}(\mathbf{x}(t_0),\mathbf{x}(t_f),t_f,\mathbf{q}) <= \mathbf{b}_{max}\]
Change of Interval¶
To can change the limits of the integration (in order to apply Quadrature), we introduce \(\tau \in [-1,+1]\) as a new independent variable and perform a change of variable for \(t\) in terms of \(\tau\), by defining:
\[t = \frac{t_f - t_0}{2}\tau + \frac{t_f + t_0}{2}\]
The optimal control problem defined above (TODO: figure out equation references), is now redefined in terms of \(\tau\) as:
Find:
- The state: \(\mathbf{x}(\tau)\)
- The control: \(\mathbf{u}(\tau)\)
- The integrals: \(\mathbf{q}\)
- The initial time: \(t_0\)
- The final time: \(t_f\)
- To Minimize:
- \[J = \Phi(\mathbf{x}(-1),\mathbf{x}(+1),\mathbf{q},t_0,t_f)\]
That Satisfy the Following Constraints:
- Dynamic Constraints:
\[\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}\tau} = \frac{t_f-t_0}{2} \mathbf{\psi}(\mathbf{x}(\tau),\mathbf{u}(\tau),\tau,t_0,t_f)\]
- Inequality Path Constraints:
\[\mathbf{c}_{min} <= \mathbf{c}(\mathbf{x}(\tau),\mathbf{u}(\tau),\tau,t_0,t_f) <= \mathbf{c}_{max}\]
- Integral Constraints:
\[q_i = \frac{t_f-t_0}{2} \int_{-1}^{+1} \Upsilon_i(\mathbf{x}(\tau),\mathbf{u}(\tau),\tau,t_0,t_f)\, \mathrm{d}\tau,\;\;(i=1,....,n_q)\]
- Event Constraints:
\[\mathbf{b}_{min} <= \mathbf{b}(\mathbf{x}(-1),\mathbf{x}(+1),t_f,\mathbf{q}) <= \mathbf{b}_{max}\]
Divide The Interval \(\tau \in [-1,+1]\)¶
- The interval \(\tau \in [-1,+1]\) is now divided into a mesh of K mesh intervals as:
- \[[T_{k-1},T_k], k = 1,...,T_K\]
- with \((T_0,...,T_K)\) being the mesh points; which satisfy:
- \[-1 = T_0 < T_1 < T_2 < T_3 < ........... < T_{K-1} < T_K = T_f = +1\]
Rewrite the Optimal Control Problem using the Mesh¶
Find:
- The state : \(\mathbf{x}^{(k)}(\tau)\) in mesh interval k
- The control: \(\mathbf{u}^{(k)}(\tau)\) in mesh interval k
- The integrals: \(\mathbf{q}\)
- The initial time: \(t_0\)
- The final time: \(t_f\)
- To Minimize:
- \[J = \Phi(\mathbf{x}^{(1)}(-1),\mathbf{x}^{(K)}(+1),\mathbf{q},t_0,t_f)\]
That Satisfy the Following Constraints:
- Dynamic Constraints:
\[\frac{\mathrm{d}\mathbf{x}^{(k)}(\tau^{(k)})}{\mathrm{d}\tau^{(k)}} = \frac{t_f-t_0}{2} \mathbf{\psi}(\mathbf{x}^{(k)}(\tau^{(k)}),\mathbf{u}^{(k)}(\tau^{(k)}),\tau^{(k)},t_0,t_f),\;\;(k=1,...,K)\]
- Inequality Path Constraints:
- Integral Constraints:
- Event Constraints:
State Continuity
Also, we must now constrain the state to be continuous at each interior mesh point \((T_1,...T_{k-1})\) by enforcing:
\[\mathbf{y}^{k}(T_k) = \mathbf{y}^{k+1}(T_k)\]
Optimal Control Problem Approximation¶
The optimal control problem will now be approximated using the Radau Collocation Method as which follows the description provided by [BGar11]. In collocation methods, the state and control are discretized at particular points within the selected time interval. Once this is done the problem can be transcribed into a nonlinear programming problem (NLP) and solved using standard solvers for these types of problems, such as IPOPT or KNITRO.
- For each mesh interval \(k\in[1,..,K]\):
- \begin{eqnarray} \mathbf{x}^{(k)}(\tau)&\approx\mathbf{X}^{(k)}(\tau)=\displaystyle\sum_{j=1}^{N_k+1}\mathbf{X}_j^{(k)}\frac{\mathrm{d}\mathcal{L}_j^{k}(\tau)}{\mathrm{d}\tau}\\ where,\;\;&\\ \mathcal{L}_j^{k}(\tau)&=\prod_{\substack{l=1 \\ l\neq j}}^{N_k+1}\frac{\tau-\tau_l^{(k)}}{\tau_j^{(k)}-\tau_l^{(k)}}\\ and,\;\;&\\ &D_{ki}=\dot{\mathcal{L}}_i(\tau_k)=\frac{\mathrm{d}\mathcal{L}_j^{k}(\tau)}{\mathrm{d}\tau} \end{eqnarray}
- also,
\(\mathcal{L}_j^{(k)}(\tau),\;\;(j=1,...,N_k+1)\) is a basis of Lagrange polynomials
\((\tau_1^{k},.....,\tau_{N_k}^{(k)})\) are the Legendre-Gauss-Radau collocation points in mesh interval k
- defined on the subinterval \(\tau^{(k)}\in[T_{k-1},T_k]\)
- \(\tau_{N_k+1}^{(k)}=T_k\) is a noncollocated point
A basic description of Lagrange Polynomials is presented in Lagrange Interpolating Polynomials
- The \(\mathbf{D}\) matrix:
- Has a size \(= [N_c]\times[N_c+1]\)
- with \((1<=k<=N_c), (1<=i<=N_c+1)\)
- this non-square shape because the state approximation uses the \(N_c+1\) points:
- \((\tau_1,...\tau_{N_c+1})\)
- but collocation is only done at the \(N_c\) LGR points:
- \((\tau_1,...\tau_{N_c})\)
If we define the state matrix as:
The dynamics are collocated at the \(N_c\) LGR points using:
\(\mathbf{D}_k\mathbf{X}^{LGR} = \frac{(t_f-t_0)}{2}\mathbf{f}(\mathbf{X}_k,\mathbf{U}_k,\tau,t_0,t_f)\;\;for\;\;k = {1,...,Nc}\)
- with,
- \(\mathbf{D}_k\) being the \(k^{th}\) row of the \(\mathbf{D}\) matrix.
References
| [BGar11] | Divya Garg. Advances in global pseudospectral methods for optimal control. PhD thesis, University of Florida, 2011. |