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:
\[\begin{split}\mathbf{c}_{min} <= \mathbf{c}(\mathbf{x}(t),\mathbf{u}(t),t) <= \mathbf{c}_{max}\end{split}\]
  • 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:
\[\begin{split}\mathbf{b}_{min} <= \mathbf{b}(\mathbf{x}(t_0),\mathbf{x}(t_f),t_f,\mathbf{q}) <= \mathbf{b}_{max}\end{split}\]

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:
\[\begin{split}\mathbf{c}_{min} <= \mathbf{c}(\mathbf{x}(\tau),\mathbf{u}(\tau),\tau,t_0,t_f) <= \mathbf{c}_{max}\end{split}\]
  • 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:
\[\begin{split}\mathbf{b}_{min} <= \mathbf{b}(\mathbf{x}(-1),\mathbf{x}(+1),t_f,\mathbf{q}) <= \mathbf{b}_{max}\end{split}\]

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:
\[\begin{split}-1 = T_0 < T_1 < T_2 < T_3 < ........... < T_{K-1} < T_K = T_f = +1\end{split}\]

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:
\[\begin{split}\mathbf{c}_{min} <= \mathbf{c}(\mathbf{x}^{(k)}(\tau^{(k)}),\mathbf{u}^{(k)}(\tau^{(k)}),\tau^{(k)},t_0,t_f) <= \mathbf{c}_{max},\;\;(k=1,...,K)\end{split}\]
  • Integral Constraints:
\[q_i = \frac{t_f-t_0}{2} \displaystyle\sum_{k=1}^{K} \int_{T_{k-1}}^{T_k} \Upsilon_i(\mathbf{x}^{(k)}(\tau^{(k)}),\mathbf{u}^{(k)}(\tau^{(k)}),\tau,t_0,t_f)\, \mathrm{d}\tau,\;\;(i=1,....,n_q, k=1,...,K)\]
  • Event Constraints:
\[\begin{split}\mathbf{b}_{min} <= \mathbf{b}(\mathbf{x}^{(1)}(-1),\mathbf{x}^{(K)}(+1),t_f,\mathbf{q}) <= \mathbf{b}_{max}\end{split}\]

  • 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. 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}\ell_j^{k}(\tau)}{\mathrm{d}\tau}\\ where,\;\;&\\ \ell_j^{k}(\tau)&=\prod_{\substack{l=1 \\ l\neq j}}^{N_k+1}\frac{\tau-\tau_l^{(k)}}{\tau_j^{(k)}-\tau_l^{(k)}} \end{eqnarray}\]
also,

  • \(\ell_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