External Constraint Handling for Solving Optimal Control Problems with Simultaneous Approaches and Interior Point Methods

TL;DR

This paper introduces a novel method for efficiently handling inactive constraints in optimal control problems, significantly reducing computational effort when using interior point methods within a simultaneous approach framework.

Contribution

The paper proposes a systematic constraint removal strategy tailored for interior point solvers in optimal control, improving computational efficiency during mesh refinement.

Findings

Over 40% reduction in computation time in example problems

Effective removal of inactive constraints without compromising solution feasibility

Enhanced efficiency in interior point-based optimal control solvers

Abstract

Inactive constraints do not contribute to the solution of an optimal control problem, but increase the problem size and burden the numerical computations. We present a novel strategy for handling inactive constraints efficiently by systematically removing the inactive and redundant constraints. The method is designed to be used together with simultaneous approaches under a mesh refinement framework, with mild assumptions that the original problem has feasible solutions, and the initial solve of the problem is successful. The method is tailored for interior point-based solvers, which are known to be very sensitive to the choice of initial points in terms of feasibility. In the example problem shown, the proposed scheme achieves more than a 40% reduction in computation time.

Tables

Table 1. TABLE I: External constraint handling history: Aircraft Flight Profile ( t 0 = subscript 𝑡 0 absent t_{0}= 0 s, t f = 7475 subscript 𝑡 𝑓 7475 t_{f}=7475  s, β = 0 𝛽 0 \beta=0  s)

	Constraint Activation Intervals Implemented in the OCP [$s$]
	MR Iteration 1	MR Iteration 2	MR Iteration 3	MR Iteration 4	MR Iteration 5	MR Iteration 6
	(K = 40)	(K = 81)	(K = 136)	(K = 176)	(K = 207)	(K = 256)
NFZ 1	[$t_0$, $t_f$]	[6900 7092]	[6996 7114]	[7056 7162]	[7071 7140]	[7074 7150]
NFZ 2/3/5	[$t_0$, $t_f$]	$\varnothing$	$\varnothing$	$\varnothing$	$\varnothing$	$\varnothing$
NFZ 4	[$t_0$, $t_f$]	[671 863]	[743 942]	[767 875]	[774 880]	[774 878]

Table 2. TABLE II: Computational performance comparison

	Standard	With ECH	With ECH
	Solve	($\beta =0$ s)	($\beta =747.5$ s)
Total Comp.	130.54	92.14	78.27
Time [s]		(29% lower)	(40% lower)
MR Iterations	6	6	6
Re-comp.	21.82	13.50	10.78
Time [s]		(38% lower)	(50% lower)
Fuel Used [kg]	1787.6	1787.6	1787.6

Equations

$$\underset{x,u,p,t_0,t_f}{\min }\Phi (x(t_0),t_0,x(t_f),t_f,p)+{\int }_{t_0}^{t_f}L(x(t),u(t),t,p)dt$$

$$\displaystyle\dot{x}(t)=f(x(t),u(t),t,p),\$$

$$\displaystyle c(x(t),u(t),t,p)\leq 0,\$$

$$\displaystyle\phi(x(t_{0}),t_{0},x(t_{f}),t_{f},p)=0,\$$

$$x^{(k)}(\tau )\approx {\bar{x}}^{(k)}(\tau ):=\sum _{j=1}^{N^{(k)}}{\mathcal{X}}_j^{(k)}{\mathcal{B}}_j^{(k)}(\tau ),$$

$$\begin{gathered} \underset{X,U,p,t_0,t_f}{\min }\Phi (X_1^{(1)},t_0,X_f^{(K)},t_f,p) \\ +\sum _{k=1}^K\sum _{i=1}^{N^{(k)}}{w}_i^{(k)}L(X_i^{(k)},U_i^{(k)},{\tau }_i^{(k)},t_0,t_f,p) \end{gathered}$$

$$\sum _{j=1}^{N^{(k)}}{\mathcal{A}}_{ij}^{(k)}{X}_j^{(k)}+{\mathcal{D}}_i^{(k)}f(X_i^{(k)},U_i^{(k)},{\tau }_i^{(k)},t_0,t_f,p)=$$

$$c(X_i^{(k)},U_i^{(k)},{\tau }_i^{(k)},t_0,t_f,p)\le$$

$$\varphi (X_1^{(1)},t_0,X_f^{(K)},t_f,p)=$$

$$y^{\ast }\in \arg \underset{y}{\min }\Phi (y)\text{subject to}c(y)\le 0,$$

$$\begin{array}{rl}\frac{\partial \Phi (y)}{\partial y}{|}_{y=y^{\ast }}+\lambda \frac{\partial c(y)}{\partial y}{|}_{y=y^{\ast }}& =0,\\ c(y^{\ast })\le 0,\lambda \ge 0,\lambda \circ c& (y^{\ast })=0,\end{array}$$

$$\begin{array}[]{ll}\text{if }\bar{\Lambda}_{\mathcal{T}_{i}}\geq\zeta&\text{constraint potentially active in interval }\mathcal{T}_{i}\\ \text{otherwise}&\text{constraint potentially inactive in interval }\mathcal{T}_{i}\\ \end{array}$$

$$J^{\ast }:=\underset{X,U,p,t_0,t_f}{\min }\sum _{l=1}^{n_g}{s}_l$$

$$\sum _{j=1}^{\tilde{N}}{\mathcal{A}}_{ij}^{(k)}{X}_j^{(k)}+{\mathcal{D}}_i^{(k)}f(X_i^{(k)},U_i^{(k)},{\tau }_i^{(k)},t_0,t_f,p)=$$

$$c(X_i^{(k)},U_i^{(k)},{\tau }_i^{(k)},t_0,t_f,p)\le s,\text{ with }s\ge$$

$$\varphi (X_1^{(1)},t_0,X_K^{(K)},t_f,p)=$$

$$\dot{h}(t)=$$

$$\dot{PO{S}_N}(t)=$$

$$\dot{PO{S}_E}(t)=$$

$$\begin{array}{rl}{\dot{v}}_T(t)=& \frac{1}{m(t)}(T(v_C(t),h(t),\Gamma (t))\\ & -D(v_T(t),h(t),\alpha (t))-m(t)g\sin (\gamma (t)))\end{array}$$

$$\begin{array}{rl}\dot{\gamma }(t)=& \frac{1}{m(t)v_T(t)}(L(v_T(t),h(t),\alpha (t))\cos (\varphi (t))\\ & -m(t)g\cos (\gamma (t)))\end{array}$$

$$\dot{\chi }(t)=$$

$$\dot{m}(t)=$$

$$\begin{array}{r}(PO{S}_N(t)-PO{S}_{N_N}{)}^2+(PO{S}_E(t)-PO{S}_{E_N}{)}^2\ge r_N^2\end{array}$$

Full text

©2019 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.

I Introduction

Optimal control has been very popular for a wide range of applications, thanks to its ability to handle various types of constraints systematically. When formulating the optimal control problem (OCP), it is common practice to impose a large number of constraints to ensure all mission specifications are fulfilled. However, for the solution obtained, it is often the case that only a small subset of the imposed inequality constraints will actually be active. Furthermore, even for the ones in this small subset, the duration for which each constraint is active is generally much shorter than the time dimension of the OCP. One example would be for the design of flight control systems: although all limits of the flight envelope need to be specified in the problem formulation for safety requirements, only in rare (abnormal) situations is it the case that some limits will be reached.

In numerical optimal control, the OCPs are transcribed into sparse nonlinear programming (NLP) problems. A distinction can be made here between simultaneous and sequential approaches, depending on whether all or just the control trajectories are discretized as decision variables [1]. For this work, we focus on the simultaneous approach.

The main computational overheads for solving the NLP problems are directly related to the number of decision variables and constraints. Thus there exist significant computational benefits to exclude inactive constraints in the problem formulation. One possibility is to only include the constraints that are determined to be active. Based on this, an external strategy for the handling of path constraints with active-set based NLP solvers has been proposed in [2]. The idea is to first solve the unconstrained problem and determine which constraints are likely to be active based on constraint violations. These constraints are then added in the OCP and the problem is repetitively solved until all original constraints are satisfied. However, a fundamental problem arises when implementing the same idea on interior point method (IPM) based solvers, since good performance hinges on the initial point to be feasible, or at least close to feasible [3].

Another option is to remove constraints that are inactive. Removal of constraints for model predictive control (MPC) has been studied to accelerate computations for linear MPC [4], tube-based robust linear MPC [5] and recently nonlinear MPC [6], with computational benefits clearly demonstrated. However, all of these are based on a quadratic regulation cost, making their application specifically aimed at receding horizon control of regulation tasks.

In this paper, we introduce an external constraint handling (ECH) strategy that is tailored to IPM-based NLP solvers for solving a variety of OCPs. Implemented together with mesh refinement (MR) schemes, constraints that do not contribute to the solution are systematically removed in the problem formulation. Special attention is paid to ensure feasibility of the initial point. As a result, significant computational savings can be achieved. Section II gives an introduction to numerical optimal control with direct collocation. This is followed by a discussion of the proposed ECH strategy in Section III. A flight control example is presented in Section IV to demonstrate the computational benefits.

II Numerical Optimal Control

Generally speaking, optimization-based control requires the solution of OCPs expressed in the general Bolza form:

[TABLE]

with $x:\mathbb{R}\rightarrow\mathbb{R}^{n}$ is the state trajectory of the system, $u:\mathbb{R}\rightarrow\mathbb{R}^{m}$ is the control input trajectory, $p\in\mathbb{R}^{s}$ are static parameters, $t_{0}\in\mathbb{R}$ and $t_{f}\in\mathbb{R}$ are the initial and terminal time. $\Phi$ is the Mayer cost functional ($\Phi$: $\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}^{s}\to\mathbb{R}$), $L$ is the Lagrange cost functional ($L:\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}\times\mathbb{R}^{s}\to\mathbb{R}$), $f$ is the dynamic constraint ($f:\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}\times\mathbb{R}^{s}\to\mathbb{R}^{n}$), $c$ is the path constraint ($c:\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}\times\mathbb{R}^{s}\to\mathbb{R}^{n_{g}}$) and $\phi$ is the boundary condition ($\phi:\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}^{s}\to\mathbb{R}^{n_{q}}$).

In practice, most optimal control problems formulated as (1) need to be solved with numerical schemes. Compared to sequential methods, simultaneous methods have some advantages with regards to computational efficiency, as well as in the treatment of path constraints and unstable dynamics [1]. In this paper, we will demonstrate our proposed ECH strategy with the simultaneous approach of direct collocation.

II-A Direct collocation methods

Direct collocation methods can be categorized into fixed-order $h$ methods [7], and variable-order $p$/$hp$ methods [8, 9]. Here, we only provide a high level overview. For a mesh of size $N:=\sum_{k=1}^{K}N^{(k)}$, the states can be approximated as

[TABLE]

within mesh interval $k$ $\in$ $\{1,\ldots,K\}$, where $N^{(k)}$ denotes the number of collocation points for interval $k$ and $\mathcal{B}_{j}^{(k)}(\cdot)$ are basis functions. For typical $h$ methods, $\tau\in\mathbb{R}^{N}$ takes values on the interval $[0,1]$ representing $[t_{0},t_{f}]$, and $\mathcal{B}_{j}^{(k)}(\cdot)$ are elementary B-splines of various orders. For $p$/$hp$ methods, $\tau$ $\in$ $[-1,1]$ and $\mathcal{B}_{j}^{(k)}(\cdot)$ are Lagrange interpolating polynomials. We use $X_{j}^{(k)}$ and $U_{j}^{(k)}$ to represent the approximated states and inputs at collocation points, e.g. $X_{j}^{(k)}=\bar{x}^{(k)}(\tau_{j}^{(k)})\in\mathbb{R}^{n}$, where $\tau_{j}^{(k)}$ is the $j^{\text{th}}$ collocation point in mesh interval $k$.

Consequently, the OCP (1) can be approximated by

[TABLE]

where $w_{j}^{(k)}$ are the quadrature weights for the chosen discretization, $\mathcal{A}^{(k)}$ is the numerical differentiation matrix with element $(i,j)$ denoted by $\mathcal{A}_{ij}^{(k)}$ and $\mathcal{D}_{i}^{(k)}$ is a row vector.

The discretized problem can then be solved with off-the-shelf NLP solvers. The NLP solver generates a discretized solution $Z\coloneqq(X,U,p,\tau,t_{0},t_{f})$ as sampled data points. Interpolating splines may be used to construct an approximation of the continuous-time optimal trajectory $t\mapsto\tilde{z}(t)\coloneqq(\tilde{x}(t),\tilde{u}(t),t,p)$. The quality of the interpolated solution needs to be assured through error analysis, assessing the level of accuracy and constraint satisfaction at a much higher resolution than the discretization mesh.

If necessary, appropriate modifications must be made to the discretization mesh, until the solutions obtained with the new mesh fulfills all predefined error tolerance levels (e.g. the absolute local error $\eta_{tol}$ and the absolute local constraint violation $\varepsilon_{c_{tol}}$). This process of MR is crucial in solving large-scale problems efficiently. For instance, it took 6 MR iterations for the example problem to be solved to a specific tolerance level. This level of accuracy was not achievable with any uniform mesh using the same desktop computer.

III External Constraint Handling

III-A Active and inactive constraints

A constraint (1c) is considered active if its presence influences the solution $z^{\ast}(\cdot)\coloneqq(x^{\ast}(\cdot),u^{\ast}(\cdot),p^{\ast},t_{0}^{\ast},t_{f}^{\ast})$. A constraint is inactive if it can be removed without affecting the solution. To clarify, consider a simplified problem:

[TABLE]

where we need to identify conditions such that constraints can be determined to be active. The most obvious criteria is when the solution $y^{\ast}$ is at the boundary of $c_{i}(y)\leq 0$, i.e. $c_{i}(y^{\ast})=0$. Additionally, consider the Lagrangian $\mathcal{L}:=\Phi(y)+\lambda^{T}c(y)$ and the necessary optimality conditions (Karush-Kuhn-Tucker (KKT) conditions):

[TABLE]

with $\circ$ the Hadamard product. From the included complementary slackness condition, we know that for strictly positive Lagrange multipliers ($\lambda_{i}>0$), the corresponding solution will have $c_{i}(y^{\ast})=0$, i.e. the constraint is active.

III-B Identifying active constraints in optimal control

Theoretically, the above-mentioned analysis applies only to the continuous OCP formulation (1). Additional challenges will arise in practice when solving the discretized problem (2) numerically: the NLP solver will only return the values of the discretized state $X$, input $U$, and Lagrange multipliers $\Lambda$ at collocation points.

To estimate the constraint activation status in-between collocation points, a criteria can be introduced based on the interpolated continuous trajectory $\tilde{z}$. By definition, inequality constraints are active if the magnitude of the differences between the actual constraint $c_{l}(\tilde{z}(\cdot))$ and the user-defined constraint bounds are zero. Due to numerical inaccuracies, however, there will always be a remainder. Thus, we consider a constraint to be potentially active if this difference is smaller than the constraint violation tolerance $\epsilon_{c_{tol}}$.

Note that the word potentially is used to emphasise that, for numerical schemes under limited machine precision, no concrete determination of constraint active status can be made. On the other hand, we know that only if the identified inactive constraints are truly inactive, then can they be removed from the OCP without affecting the solutions. Thus, it would be much more preferable to erroneously identify inactive constraints as active, than the opposite situation.

For this reason, we also use the multiplier information to enforce a larger (more conservative) selection of potentially active constraints. Here, a similar numerical challenge arises: with limited machine precision, even when the corresponding constraints are inactive, the multiplier values are rarely truly equal to zero. To identify the regions where the constraints are likely to be active, the numerical multiplier data $\Lambda$ is first normalized between 0 and 1 for each constraint $c_{l}(Z)\leq 0$. Signal processing algorithms can be used to identify different intervals where the behaviour of Lagrange multipliers have significant changes, for example using the Matlab findchangepts function.

For each identified interval $\mathcal{T}_{i}$, the mean value of the normalized multipliers ($\bar{\Lambda}_{\mathcal{T}_{i}}$) is calculated and compared based on the following criteria:

[TABLE]

with $\zeta$ a threshold parameter.

To sum up, the following definition is used to determine whether the constraints are potentially active or potentially inactive at different collocation points.

Definition 1

A constraint $c_{l}(X_{i},U_{i},\tau_{i},t_{0},t_{f},p)\leq 0$ is potentially active at time $t_{i}:=t_{i}(\tau_{i},t_{0},t_{f})$ if one of the following criteria is met:

• •

Between adjacent collocation points ($t\in[t_{i-1},t_{i+1}]$), $c_{l}(\tilde{x}(t),\tilde{u}(t),t,p)\geq-\epsilon_{c_{tol}}$ holds, with $\epsilon_{c_{tol}}>0$.

• •

$\bar{\Lambda}_{\mathcal{T}_{i}}\geq\zeta$, with $t_{i}\in\mathcal{T}_{i}$.

Otherwise, a constraint is potentially inactive at time $t_{i}$.

In addition to identifying the time instances at which certain constraints may be potentially inactive, it is also preferable to determine the sets of constraints that never become active at all times.

Definition 2

A constraint $c_{l}(X_{i},U_{i},\tau_{i},t_{0},t_{f},p)\leq 0$ is potentially redundant if for all $t_{i}:=t_{i}(\tau_{i},t_{0},t_{f})$ with $i=1,\ldots,N$, the constraints $c_{l}(X_{i},U_{i},\tau_{i},t_{0},t_{f},p)\leq 0$ are potentially inactive. Otherwise, this set of constraints is potentially enforced.

III-C Initialization for interior point methods

Interior point methods (IPMs) for solving NLPs were introduced in the early 1960s [10, 11, 12] and have became very popular in numerical optimal control. The idea is to augment the objective function with barrier functions of constraints in order to enforce their satisfaction. Potential solutions will iterate only in the feasible region following the so-called central path, resulting in a very efficient algorithm.

Standard interior point methods are sensitive to the choice of a starting point. To ensure that the initial guess is strictly feasible with respect to constraints, various initialization methods have been developed (e.g. [13] and the collective study in [7]) and implemented in modern solvers.

To ensure reliable and efficient computation of the initialization algorithm, as well as the subsequent NLP iterations, several criteria [3] can be formulated regarding ideal initial points for IPMs. The ideal initial point should:

• •

satisfy or be close to primal and dual feasibility,

• •

be close to the central path,

• •

be as close to optimality as possible.

Because of these characteristics, external constraint handling schemes developed for active-set based SQP solvers, such as [2], are not suitable for IPM-based solvers. By first solving the unconstrained problem and gradually adding constraints based on the constraint violation error, the solution of previous solves will all be infeasible for the new OCP formulation, and the solution may undergo drastic changes as well. For IPM-based NLP solvers, this would lead to a higher computational overhead for initialization, as well as higher chances for the iterations to frequently enter the slow and unreliable feasibility restoration phase.

III-D Proposed Scheme for Constraint Handling

Based on the criteria presented in Section III-B and the characteristics of IPMs as discussed in Section III-C, a strategy for efficiently handling constraints in OCPs solved with IPM-based NLP solvers is proposed, with the work-flow presented in Figure 1. The approach is called external, since the modifications to the OCP are made at the MR iteration level, instead of during the NLP iterations.

The unmodified OCP is first solved on the initial coarse mesh. Even with all constraint equations included, the computation time will still be quite low at this stage due to the small problem size. Once the solution is obtained, potentially inactive constraints and potentially redundant constraint sets can be identified, based on Definitions 1 and 2, with potentially redundant constraints directly excluded from the OCP formulation. Furthermore, if the problem has a fixed terminal time, i.e. the time instance corresponding to a mesh point will not change, then potentially inactive constraints in the potentially enforced constraint sets may also be removed.

Recall that it is preferable to erroneously identify an inactive constraint as potentially active, rather than the opposite. It is therefore often a good idea in practice to enlarge the intervals with potential constraint activation by an interval of length $\beta$ in each direction, with $\beta$ either fixed or adapting during the MR process. This adaptation also guarantees the convergence of the overall scheme, i.e. in the worst case, $\beta$ can be sufficiently large to impose the constraints for the whole trajectory, with the original problem recovered.

In-between MR iterations, special attention must be made to constraints and constraint sets that were determined to be potentially inactive or redundant in the previous solves. If they never become active or enforced again, the constraint removal process may continue until MR is converged. But there will be the chance, after refining the mesh, the constraint violation error analysis dictates that certain constraints and constraint sets that have been removed earlier may become potentially active or enforced again. If this happens, they need to be included again to ensure that the solution of the modified OCP is equivalent to the unmodified problem.

Note that the previous solution will no longer be a feasible initial guess for the new problem formulation. To assist the subsequent solve of NLPs, the following auxiliary feasibility problem (AFP) can be solved before proceeding:

[TABLE]

with $s\in\mathbb{R}^{n_{g}}$ slack variables. The initial guess for the AFP will be $\tilde{Z}\coloneqq(\tilde{X},\tilde{U},p,\tilde{\tau},t_{0},t_{f})$, the values of the interpolated solution $\tilde{z}$ at the collocation points of the refined mesh.

III-E Properties of the external constraint handling strategy

It is possible to derive proofs of feasibility and optimality invariance for the removal of constraints on a given discretization mesh. However, with the size of the mesh changing throughout the refinement process, the analysis of errors, and identification of constraint activation status (Section III-B) are all subject to considerable uncertainties. When a constraint or constraint set must be included again in the problem, it will be challenging to ensure that the subsequent OCP solve can be supplied with a feasible initial guess. The introduction of the AFP is the answer to this challenge, with its solution guaranteed to be a feasible point for the corresponding original OCP.

Proposition 1

If the original OCP (2) has feasible points, then a solution to the auxiliary feasibility problem (3) will be a feasible point of (2) on the same discretization mesh, and (3) will have corresponding objective value $J^{\ast}=0$.

Proof:

If $Z\coloneqq(X,U,\tau,t_{0},t_{f},p)$ is a feasible point of (2), then (2c) must hold. With $s\geq 0$, the solution for the AFP (3) will be the situation where $\sum s_{l}=0$, and (2c) guarantees the existence of such a solution. ∎

Now, for the very same reason, we need to obtain a suitable initial guess for the slack variables $s$ in the AFP. One possible way is by calculating the constraint violation errors of the interpolated solutions on the refined mesh.

Proposition 2

Define $\tilde{s}\in\mathbb{R}^{N\times n_{g}}$ as the absolute local constraint violation error $\epsilon_{c}(t)$ calculated at the collocation points of the refined mesh, with the updated initial guess $\tilde{Z}\coloneqq(\tilde{X},\tilde{U},p,\tilde{\tau},t_{0},t_{f})$. For any set $\{\bar{s}\in\mathbb{R}^{n_{g}}\mid\bar{s_{l}}\geq\max_{i=1,\dots,N}(\tilde{s}_{i,l}),l=1,\ldots,n_{g}\}$ implemented as the initial guess for $s$, the AFP (3) will have a strictly feasible initial point with respect to the constraints (3c).

Proof:

$\tilde{s}_{i,l}:=|\min(-c_{l}(\tilde{X}_{i}$, $\tilde{U}_{i}$, $p$, $\tilde{\tau}_{i}$, $t_{0}$, $t_{f}),0)|$ by definition, for all $i=1,\ldots,N$ and $l=1,\ldots,n_{g}$,

• •

if $c_{l}(\tilde{X}_{i}$, $\tilde{U}_{i}$, $p$, $\tilde{\tau}_{i}$, $t_{0}$, $t_{f})<0$, i.e. the constraint is satisfied and the solution is not on the boundary, then $\tilde{s}_{i,l}=0$, thus $c_{l}(X_{i},U_{i},\tau_{i},t_{0},t_{f},p)<\tilde{s}_{i,l}$ holds.

• •

if $c_{l}(\tilde{X}_{i}$, $\tilde{U}_{i}$, $p$, $\tilde{\tau}_{i}$, $t_{0}$, $t_{f})=0$ (constraint satisfied and solution is on the boundary) or $c_{l}(\tilde{X}_{i}$, $\tilde{U}_{i}$, $p$, $\tilde{\tau}_{i}$, $t_{0}$, $t_{f})>0$ (constraint violation occurs), then $c_{l}(X_{i},U_{i},\tau_{i},t_{0},t_{f},p)=\tilde{s}_{i,l}$ holds.

Therefore, $c_{l}(X_{i},U_{i},\tau_{i},t_{0},t_{f},p)\leq\tilde{s}_{i,l}$ will always be true. From $\bar{s}_{l}\geq\max_{i=1,\dots,N}(\tilde{s}_{i,l})$, it can be concluded that $c_{l}(X_{i},U_{i},\tau_{i},t_{0},t_{f},p)\leq\bar{s}_{l}$ holds. ∎

We can now show that, except for the initial solve, all subsequent solves will have feasible initial guesses.

Proposition 3

If the unmodified OCP has feasible points, and the initial solve of the discretized OCP has been successful, then all subsequent solves of OCPs and AFPs with mesh refinement schemes and the proposed external constraint handling method will have a feasible initial point with respect to the constraints (2c) or (3c).

Proof:

For any interpolated solution $\tilde{Z}\coloneqq(\tilde{X},\tilde{U},p,\tilde{\tau},t_{0},t_{f})$ on the new mesh, if (2c) is not satisfied, then the corresponding AFP will be solved and Proposition 2 ensures that AFP will have a feasible initial guess. From Proposition 1 the solution of the AFP will be a feasible initial guess for the subsequent OCP solve. ∎

III-F A practically more efficient alternative implementation

The proposed external constraint handling scheme can guarantee feasible initial points under conditions stated in Proposition 3. Nevertheless, it is not efficient in practice. The frequent solve of AFPs are not only time consuming, but are often not necessary.

Recall the conditions regarding ideal initial guesses for IPM methods. It is not necessary to satisfy primal and dual feasibility — rather, one only needs to be close to fulfillment. In practice, the computational performance of a modern IPM that uses near-feasible initial guesses is very much comparable to using feasible initial points. In addition, constraint satisfaction for simple bounds can be computationally much easier to achieve by the NLP solver, thus there is no need to enforce those through the solve of an AFP.

Thus, a practically more efficient version of the external handling scheme can be formulated, by restricting the conditions for solving the AFP to the MR iteration when a potentially redundant path constraint set turns into a potentially enforced path constraint set.

IV Example

To demonstrate the computational benefits of the proposed ECH scheme, we show an problem that is relatively large in the horizon length. The task involves finding a fuel-optimal flight path of a commercial aircraft where authorities have identified five non-flight zones (NFZ) for the aircraft to avoid.

From simple flight mechanics with a flat earth assumption, both the longitudinal and lateral motion of the aircraft can be described by the dynamic equations

[TABLE]

with $h$ the altitude [m], $POS_{N}$ and $POS_{E}$ the north and east position [m], $v_{T}$ the true airspeed [m/s], $\gamma$ the flight path angle [rad], $\chi$ the tracking angle [rad], and $m$ the mass [kg]. $T$, $L$, $D$ are the thrust, lift and drag forces. $FF$ is the fuel flow model, requiring an input of calibrated airspeed $v_{C}$, which can be related to $v_{T}$ via a conversion.

Additionally, $g=9.81$ m/s2 is gravitational acceleration. We have three control inputs, the roll angle $\phi$ in [rad], the throttle settings $\Gamma$ normalized between 0 and 1, and the angle of attack $\alpha$ in [rad]. Further details of the modelling of a Fokker 50 aircraft can be obtained from [14].

The avoidance of NFZs can be implemented with the following path constraints

[TABLE]

with $POS_{N_{N}}$ and $POS_{E_{N}}$ the north and east position of the center of the non-flight zones, and $r_{N}$ the radius.

The problem will have the boundary cost $\Phi(x(t_{0}),t_{0},x(t_{f}),t_{f},p)=-m(t_{f})$ (maximize the mass at the end of the flight, with fixed $t_{f}=7475$ s), subject to the dynamics and path constraints. Furthermore, variable simple bounds are imposed together with the boundary conditions.

The OCP is transcribed using the optimal control software ICLOCS2 [15] with Hermite-Simpson discretization, and solved with IPM-based NLP solver IPOPT [16] (version 3.12.4). All computation results shown were obtained on an Intel Core i7-4770 computer with 16G of RAM. Figure 2 illustrates the results solved to a user-defined tolerance.

Using the proposed external constraint handling scheme, we first solved the problem with a worst-case buffer interval setting of $\beta=0$ s. The history for constraint activation intervals implemented in the OCP are demonstrated in Table I. It can be seen that in the initial solve (MR iteration 1), all constraint sets are enforced and all constraints are treated as potentially active. Based only on the solution from this coarse grid, the ECH method correctly identified that the constraint sets related to NFZ 2, 3 and 5 are all potentially redundant. It also determined that constraints related to NFZ 1 are only potentially active near the end of the flight, whereas for NFZ 4 they are at the beginning of the mission.

In later iterations of the MR, these intervals had only some minor adjustments. It can be seen that without implementing any buffer interval, we do see occasions where active constraints got erroneously identified as inactive for finer meshes. However, the constraint violation error analysis in the MR process correctly identified these situations and made corrections accordingly.

Table II compares the computational performance of the standard solve, as well as solves with external constraint handling using the alternative ECH implementation (allowing near-feasible initial guesses) described in Section III-F, with two different buffer interval settings. With the worst-case setting of $\beta=0$, the total computation time saw a 29% reduction, while the number of MR iterations remained the same. Choosing a much more conservative buffer interval setting of $\beta=0.1(t_{f}-t_{0})=747.5$ s further improved this time reduction to 40%, due to the fact that initial guesses were feasible for all later MR iterations.

For real-time applications, it is useful to consider the re-computation time for solving the OCP problem again with the final (refined) discretization mesh, using the obtained solutions as initial guesses. For the ECH method with $\beta=747.5$ s, the time taken was only half compared to the standard solve. Therefore, the benefits of the proposed scheme can be seen for both off-line and online applications.

V Conclusions

A strategy has been developed to systematically identify and handle inactive constraints and redundant constraint sets for numerically solving optimal control problems together with mesh refinement schemes. Unlike previous work that would always result in infeasible initial guesses for intermediate steps, the proposed scheme is capable of providing guarantees on the feasibility of initial points in mesh refinement iterations. The method only requires some mild conditions — the original OCP to have feasible points, and the initial solve of the discretized OCP to be successful, making it particularly suitable for OCP toolboxes that utilize IPM-based NLP solvers in lowering the computational cost.

Due to limitations in time and space, we only illustrated the proposed external constraint handling method with direct collocation. We note that similar benefits should be obtainable, after some adaptations, with other simultaneous methods, such as direct multiple shooting. This might also be possible for sequential methods, such as direct single shooting. Moreover, it would be interesting to test a slightly altered version of our proposed constraint removal method with active-set based NLP solvers — this could be compared under a mesh refinement scheme against earlier work [2], which exploits constraint addition instead.