F-8 aircraft (AMPL)

This page contains a discretized version of the MIOCP F-8 aircraft in AMPL format. You should be aware of the comments regarding discretization made on the AMPL overview page. Note that you will need to include two generic support AMPL files, ampl_general.mod and ampl_general.dat.

Model

The main difficulty in calculating a time-optimal solution for the F-8 aircraft problem is the determination of the correct switching structure and of the switching points. If we want to formulate a MINLP, we have to slightly modify this problem. Our aim is not a minimization of the overall time, but now we want to get as close as possible to the origin (0,0,0) in a prespecified time $t_f = 3.78086$ on an equidistant time grid. As this time grid is not a superset of the one used for the time-optimal solution, one can not expect to reach the target state exactly.

AMPL

The model in AMPL code for a fixed control discretization grid with a collocation method. We need a model file f8_ampl.mod,

# --------------------------------------------------------------
# F-8 aircraft control problem with collocation (implicit Euler)
# (c) Sebastian Sager
# --------------------------------------------------------------
var x {I, 1..nx};
param xi > 0;
 
minimize Deviation: sum {i in 1..3} x[nt,i]*x[nt,i]; 
 
subj to ODE_DISC_1 {i in I diff {0}}:
x[i,1] = x[i-1,1] + (dt[uidx[i]]/ntperu) * (
 - 0.877*x[i,1] + x[i,3] - 0.088*x[i,1]*x[i,3] + 0.47*x[i,1]*x[i,1] 
 - 0.019*x[i,2]*x[i,2]
 - x[i,1]*x[i,1]*x[i,3] + 3.846*x[i,1]*x[i,1]*x[i,1]
 + 0.215*xi - 0.28*x[i,1]*x[i,1]*xi + 0.47*x[i,1]*xi^2 - 0.63*xi^2
 - 2*w[uidx[i]] * ( 0.215*xi - 0.28*x[i,1]*x[i,1]*xi - 0.63*xi^3 )
);
 
subj to ODE_DISC_2 {i in I diff {0}}:
   x[i,2] = x[i-1,2] + (dt[uidx[i]]/ntperu) * x[i,3];
 
subj to ODE_DISC_3 {i in I diff {0}}:
x[i,3] = x[i-1,3] + (dt[uidx[i]]/ntperu) * (
 - 4.208*x[i,1] - 0.396*x[i,3] - 0.47*x[i,1]*x[i,1] 
 - 3.564*x[i,1]*x[i,1]*x[i,1]
 + 20.967*xi - 6.265*x[i,1]*x[i,1]*xi + 46*x[i,1]*xi^2 - 61.4*xi^3
 - 2*w[uidx[i]] * (20.967*xi - 6.265*x[i,1]*x[i,1]*xi - 61.4*xi^3)
);

a data file f8_ampl.dat,

# Algorithmic parameters
param ntperu := 500;
param nu := 60;
param nt := 30000;
param nx := 3;
param fix_w := 0;
param fix_dt := 1;
 
# Problem parameters
param xi := 0.05236;
param T := 8;
 
# Initial values differential states
let x[0,1] := 0.4655;
let x[0,2] := 0.0;
let x[0,3] := 0.0;
for {i in 1..3} { fix x[0,i]; }
 
# Initial values control
let {i in U} w[i] := 0.0;
for {i in 0..(nu-1) / 2} { let w[i*2] := 1.0; }
let {i in U} dt[i] := 3.78086 / nu;

and a running script f8_ampl.run,

model ampl_general.mod;
model ampl_f8.mod;
data ampl_f8.dat;
data ampl_general.dat;
 
option presolve_eps 1e-10;
option solver ...;
 
solve;
 
param myt;
let myt := 0;
for {i in I} { 
	if ( dt[uidx[i]] > 1e-6 ) then {
		print myt, w[uidx[i]], x[i,1], x[i,2] > resultInt.txt;
	}
	let myt := myt + (dt[uidx[i]]/ntperu);
}
 
display w;

Results

The solution calculated by Bonmin (subversion revision number 1453, default settings, 3 GHz, Linux 2.6.28-13-generic, with ASL(20081205)) has an objective function value of $\Phi= 0.023405$ , while the optimum of the relaxation is $\Phi=0.023079$ . Bonmin needs 85702 iterations and 7031 nodes (64282 seconds). Strong branching is done 45 times (1212 iterations), with 0 fathomed nodes and 4 fixed variables. The intervals on the equidistant grid on which $w(t) = 1$ holds, counting from 0 to 59, are 0, 1, 31, 32, 42, 52, and 54.