6.5. Mixed-integer nonlinear solver: F8 Crusader aircraft

In this example we illustrate the simplicity of the high-level user interface on a mixed-integer nonlinear program. In particular, we use an F8 Crusader aircraft model described by a set of ordinary differential equations (ODEs):

\[\begin{split}\dot{x}_0 = & -0.877 x_0 + x_2 - 0.088 x_0 x_2 + 0.47 x_0^{2} - 0.019 x_1^{2} - x_0^{2} x_2 + 3.846 x_0^{3} \\ & -0.215 w + 0.28 x_0^2w + 0.47 x_0w^{2} + 0.63 w^{3} \\ \dot{x}_1 = & x_2 \\ \dot{x}_2 = & -0.4208 x_0 - 0.396 x_2 - 0.47 x_0^{2} - 3.564 x_0^{3} - 20.967 w \\ & + 6.265 x_0^{2}w + 46 x_0w^{2} + 61.4 w^{3}\end{split}\]

The model is taken from [GarJor77] and consists of three differential states: \(x_0\) the angle of attack in radians, \(x_1\) the pitch angle in radians and \(x_2\) the pitch angle rate in radians per second. There is one control input \(w\), the tail deflection angle in radians. The input is the discrete component of the model, since it can take values within the discrete set \(\{-0.05236, 0.05236\}\). This makes the solution process more complicated in comparison to a nonlinear program, as the different combinations of inputs have to be checked over the control horizon.

The trajectory of the aircraft is to be computed by solving a mixed-integer nonlinear program (MINLP). First, we define the stage variable \(z\) by stacking the input and differential state variables:

\[z = [w, x_0, x_1, x_2]^\top\]

You can download the Matlab code of this example to try it out for yourself by clicking here.

6.5.1. Defining the problem data

6.5.1.1. Objective

Our goal is to minimize the distance of the final state to the origin, which can be translated in the following cost function on the final stage variable:

\[f(z) = 150 x_0^{2} + 5 x_1^{2} + 5 x_2^{2}\]

The terminal cost function is coded in MATLAB as the following function:

model.objectiveN = @(z) 150 * z(2)^2 + 5 * z(3)^2 + 5 * z(4)^2;

Moreover, control inputs are penalized at every stage via the following stage cost function:

model.objective = @(z) 0.1 * z(1)^2;

6.5.1.2. Equality constraints

In this example, the only equality constraints are related to the dynamics. They are provided to FORCES Pro in continuous form. The discretization is then computed internally by the FORCES Pro integrators.

In the code snippet below, it is important to notice that the control input \(w\) is replaced with \(u\) such that

\[w~~\hat{=}~~0.05236~\cdot~(2 u - 1)\]

If \(w\) has values within \(\{-0.05236, 0.05236\}\), then \(u\) lies within the binary set \(\{0, 1\}\).

wa = 0.05236;
wa2 = wa^2;
wa3 = wa^3;
continuous_dynamics = @(x, u) [-0.877 * x(1) + x(3) - 0.088 * x(1) * x(3)...
                               + 0.47 * x(1) * x(1) - 0.019 * x(2) * x(2)...
                               - x(1) * x(1) * x(3)...
                               + 3.846 * x(1) * x(1) * x(1)...
                               - 0.215 * wa * (2 * u(1) - 1)...
                               + 0.28 * x(1) * x(1) * wa * (2 * u(1) - 1)...
                               + 0.47*x(1)*wa2*(2*u(1)-1)*(2*u(1)-1)...
                               + 0.63*wa3*(2*u(1)-1)*(2*u(1)-1)*(2*u(1)-1);
                               x(3);
                               -4.208*x(1) - 0.396 * x(3) - 0.47 * x(1)*x(1)...
                               - 3.564 * x(1) * x(1) * x(1)...
                               - 20.967 * wa * (2 * u(1) - 1)...
                               + 6.265 * x(1) * x(1) * wa * (2 * u(1) -1 )...
                               + 46.0 * x(1)*wa2*(2*u(1)-1)*(2*u(1)-1)...
                               + 61.4*wa3*(2*u(1)-1)*(2*u(1)-1)*(2*u(1)-1)];
model.continuous_dynamics = continuous_dynamics;
model.E = [zeros(3, 1), eye(3)];

6.5.1.3. Inequality constraints

The maneuver is subjected to a set of constraints, involving only the simple bounds:

\begin{align*} 0\,\mathrm{rad}\leq & u \leq 1\,\mathrm{rad} \\ -10\,\mathrm{rad}\leq & x_0 \leq 10\,\mathrm{rad} \\ -10\,\mathrm{rad}\leq & x_1 \leq 10\,\mathrm{rad} \\ -10\,\mathrm{rad/sec}\leq & x_2 \leq 10\,\mathrm{rad/sec} \\ \end{align*}

6.5.1.4. Initial and final conditions

The goal of the maneuver is to steer the aircraft from an initial condition with nose pointing upwards

\[(0.4655, 0, 0)^{T}\]

to the origin.

6.5.2. Defining the MPC problem

With the above de fined MALTAB functions for objective and equality constraints, we can completely define the MINLP formulation in the next code snippet. For this example, the number of stages has been set to \(N = 100\).

In the code snippet below, it is important to notice that the lower and upper bounds are declared as parametric before generating the solver. This needs to be done for generating mixed-integer NLP solvers. Lower and upper bounds are meant to be provided at run-time.

%% Problem dimension
nx = 3;
nu = 1;
nz = nx + nu;
model.N = 100;
model.nvar = nz;
model.neq = nx;

%% Indices of initial state in stage variable
model.xinitidx = nu+1:model.nvar;

%% Lower and upper bound need to be set as parametric for generating an MINLP solver
model.lb = [];
model.ub = [];
model.lbidx{1} = 1 : nu;
model.ubidx{1} = 1 : nu;
for i = 2 : model.N
  model.lbidx{i} = 1 : model.nvar;
  model.ubidx{i} = 1 : model.nvar;
end

%% Dynamics
wa = 0.05236;
wa2 = wa^2;
wa3 = wa^3;
continuous_dynamics = @(x, u) [-0.877 * x(1) + x(3) - 0.088 * x(1) * x(3)...
                               + 0.47 * x(1) * x(1) - 0.019 * x(2) * x(2)...
                               - x(1) * x(1) * x(3)...
                               + 3.846 * x(1) * x(1) * x(1)...
                               - 0.215 * wa * (2 * u(1) - 1)...
                               + 0.28 * x(1) * x(1) * wa * (2 * u(1) - 1)...
                               + 0.47 *x(1)*wa2*(2*u(1)-1)*(2*u(1)-1)...
                               + 0.63*wa3*(2*u(1)-1)*(2*u(1)-1)*(2*u(1)-1);
                                x(3);
                                -4.208 * x(1) - 0.396 * x(3)...
                                - 0.47 * x(1) * x(1)...
                                - 3.564 * x(1) * x(1) * x(1)...
                                - 20.967 * wa * (2 * u(1) - 1)...
                                + 6.265*x(1)*x(1)*wa*(2*u(1)-1)...
                                + 46.0*x(1)*wa2*(2*u(1)-1)*(2*u(1)-1)...
                                + 61.4*wa3*(2*u(1)-1)*(2*u(1)-1)*(2*u(1)-1)];
model.continuous_dynamics = continuous_dynamics;
model.E = [zeros(nx, nu), eye(nx)];

%% Objective
mode.objective =  @(z) 0.1 * z(nu)^2;
model.objectiveN = @(z) 150 * z(nu+1)^2...
                      + 5 * z(nu+2)^2...
                      + 5 * z(nu+3)^2;

%% Indices of integer variables within every stage
for s = 1:model.N
  model.intidx{s} = [1];
end

6.5.3. Generating an MINLP solver

We have now populated model with the necessary fields to generate a mixed-integer solver for our problem. Now we set some options for our solver and then use the function FORCES_NLP to generate a solver for the problem defined by model with the initial state and the lower and upper bounds as a parameters:

%% Set code-generation options
codeoptions = getOptions('F8aircraft');
codeoptions.printlevel = 1;
codeoptions.misra2012_check = 1;
codeoptions.maxit = 2000;
codeoptions.timing = 0;
codeoptions.nlp.integrator.type = 'IRK2';
codeoptions.nlp.integrator.Ts = 0.05;
codeoptions.nlp.integrator.nodes = 20;

%% Generate MINLP solver
FORCES_NLP(model, codeoptions);

In the code snippet above, we have set some integrator options, since the continuous-time dynamics has been provided in the model. The branch-and-bound search can be run on several threads in parallel by setting the run-time parameter numThreadsBnB equal to the number of threads to be used. The default value is 1. Moreover, the maximum number of threads for the branch-and-bound search can be set via the option max_num_threads. By default, max_num_threads = 4.

6.5.4. Calling the generated MINLP solver

Once all parameters have been populated, the MEX interface of the solver can be used to invoke it:

%% Set run-time parameters
problem.(sprintf('lb%02d', 1)) = 0;
problem.(sprintf('ub%02d', 1)) = 1;
for s = 2:99
  problem.(sprintf('lb%02d', s)) = [0, -1e1 * ones(1, 3)]';
  problem.(sprintf('ub%02d', s)) = [1, 1e1 * ones(1, 3)]';
end
problem.(sprintf('lb%02d', 100)) = [0, -1e1 * ones(1, 3)]';
problem.(sprintf('ub%02d', Nstages)) = [1, 1e1 * ones(1, 3)]';

problem.x0 = repmat([0; zeros(3, 1)], 100, 1);
problem.xinit = zeros(3, 1);
problem.xinit(1) = 0.4655;

%% Call MINLP solver
[sol, exitflag, info] = F8aircraft(problem);

6.5.5. Results

The control objective is to drive the angle of attack as close as possible to zero within a five seconds time frame. The control input is the tail deflection angle, which can take values with the set \(\{-0.05236, 0.05236\}\) and the initial state is \((0.4655, 0, 0)^{T}\), where the first component is the angle of attack, the second component is the pitch angle and the third component is the pitch angle rate.

The angle of attack computed by FORCES Pro MINLP solver running on one thread is shown in Figure Figure 6.7 and the input sequence is in Figure Figure 6.8. One can notice the bang-bang behaviour of the solution. When running on three threads the FORCES Pro MINLP solver provides a solution with lower final primal objective. Results are shown on Figures Figure 6.9 and Figure 6.10.

../../_images/f8States_1thread.png

Figure 6.7 Aircraft’s angle of attack over time computed with one thread.

../../_images/f8Input_1thread.png

Figure 6.8 Aircraft’s tail deflection angle over time with one thread.

../../_images/f8States_3threads.png

Figure 6.9 Aircraft’s angle of attack over time computed with three threads.

../../_images/f8Input_3threads.png

Figure 6.10 Aircraft’s tail deflection angle over time with three threads.

[GarJor77]Garrard, W.L.; Jordan, J.M..: Design of Nonlinear Automatic Control Systems. In: Automatica 1977, vol. 13, 497-505.