Tutorial for GPOPS-II

The Very Beginning

GPOPS-II is a MATLAB software based on Gauss pseudo-spectral method, mainly used on solving optimal control problems (OCP). This blog is a summary of the documentation of GPOPS-II, including the structure of the program and frequently-used parameters, mainly based on my own projects completed with GPOPS-II. Thanks for kunpeng’s blog and the user’s guide from GPOPS-II official website !

There’s something you should know in advance:

1. GPOPS-II divide the whole time interval into multiple time points. The state and control corresponding to each time point are a row vector, forming the matrix for calculation arranging in rows.
2. Here’s the symbol list:

Symbol	Meaning
$p$	Order of phase
$N^{(p)}$	Length of time sequence for each phase
$M$	Number of guessing points
$g$	Number of eventgroups
$n_b$	Dimension of each eventgroup
$n_y$	Dimension of state vector for each phase
$n_u$	Dimension of control vector for each phase
$n_s$	Number of static parameters for each phase
$n_c$	Number of path constraints for each phase
$n_d$	Dimension of integrands for each phase

1. For numerical examples, please go to user’s guide for help.

Program Structure of GPOPS-II

A complete GPOPS-II solving program is shown below:

%-------------------------------------------------------------%
% ------------------- Variables defination -------------------%
%-------------------------------------------------------------%
% Bounds strucure assignment
% Insert your code ...

% Guess structure assignment
% Insert your code ...

% Auxdata structure assignment
% Insert your code ...


%-------------------------------------------------------------%
% -------------------------- Setup ---------------------------%
%-------------------------------------------------------------%
% Name
setup.name = "OCP-Name";

% Functions
setup.functions.continuous = @OCPContinuous;
setup.functions.endpoint = @OCPEndpoint;

% Bounds, Guess, auxdata
setup.bounds = bounds;
setup.guess = guess;
setup.auxdata = auxdata;

% Other setup parameters
% Insert your code ...


%-------------------------------------------------------------%
% -------------------------- Solve ---------------------------%
%-------------------------------------------------------------%
output = gpops2(setup);


%-------------------------------------------------------------%
% -------------------- Continuous Function -------------------%
%-------------------------------------------------------------%
function output = OCPContinuous(input)
% Defination
% ...
end


%-------------------------------------------------------------%
% --------------------- Endpoint Function --------------------%
%-------------------------------------------------------------%
function output = OCPEndpoint(input)
% Defination
% ...
end

Continuous function: Dynamics equations

function output = OCPContinuous(input)

Components of input:

• input.phase(p).time: time vector for $p^{th}$ phase, $size=N^{(p)}\times 1$
• input.phase(p).state: state matrix for $p^{th}$ phase, $size=N^{(p)}\times n_y^{(p)}$
• input.phase(p).control: control matrix for $p^{th}$ phase, $size=N^{(p)}\times n_u^{(p)}$
• input.parameter: static parameter for $p^{th}$ phase, $size=N^{(p)}\times n_s$
• input.auxdata

Components of output:

• output(p).dynamics: dynamic equations for $p^{th}$ phase, $size=N^{(p)}\times n_y^{(p)}$
• output(p).path: path constraint for $p^{th}$ phase, $size=N^{(p)}\times n_c^{(p)}$
• output(p).integrand: integrand for $p^{th}$ phase, $size=N^{(p)}\times n_d^{(p)}$

Endpoint function:

function output = OCPEndpoint(input)

Components of input:

• input.phase(p).initialtime: initial time for $p^{th}$ phase
• input.phase(p).finaltime: final time for $p^{th}$ phase
• input.phase(p).initialstate: initial state for $p^{th}$ phase, $size=1\times n_y^{(p)}$
• input.phase(p).finalstate: final state for $p^{th}$ phase, $size=1\times n_y^{(p)}$
• input.phase(p).integral: integral for $p^{th}$ phase, $size=1\times n_d^{(p)}$
• input.parameter: static parameter for $p^{th}$ phase, $size=1\times n_s$
• input.auxdata

Components of output:

• output.objective: Objective function

• output.eventgroup(g): event constraints, $size=1\times n_b^{(g)}$

Setup

setup is a structure composed of multiple structures. Therefore, for necessary components, we strongly recommend that you define the substructures in advance and use them to assign values to the components of setup.

Necessary:

% Name
setup.name = "OCP-Name";

% Functions
setup.functions.continuous = @OCPContinuous;
setup.functions.endpoint = @OCPEndpoint;

% Bounds, Guess, auxdata
setup.bounds = bounds;
setup.guess = guess;
setup.auxdata = auxdata;

Bounds:

Each variable declared should have its lower bound and upper bound defined, or GPOPS-II could not recognize the boundary for sprinkling.

Be sure that the boundary not be too strict, or GPOPS-II cannot give a correct solution (SNOPT INFO -31) or you may encounter numerical difficulties (SNOPT INFO -41).

• Time Bounds:

bounds.phase(p).initialtime.lower = t0_lower;
bounds.phase(p).initialtime.upper = t0_upper;
bounds.phase(p).finaltime.lower = tf_lower;
bounds.phase(p).finaltime.upper = tf_upper;
bounds.phase(p).durations.lower = dt_lower;
bounds.phase(p).durations.upper = dt_upper;

• State Bounds: $size=1\times n_y^{(p)}$

bounds.phase(p).initialstate.lower = state0_lower;
bounds.phase(p).initialstate.upper = state0_upper;
bounds.phase(p).finalstate.lower = statef_lower;
bounds.phase(p).finalstate.upper = statef_upper;
bounds.phase(p).state.lower = state_lower;
bounds.phase(p).state.upper = state_upper;

• Control Bounds: $size=1\times n_u^{(p)}$

bounds.phase(p).control.lower = u_lower;
bounds.phase(p).control.upper = u_upper;

• Path Constraint Bounds: $size=1\times n_c^{(p)}$

bounds.phase(p).path.lower = path_lower;
bounds.phase(p).path.upper = path_upper;

• Integral Bounds: $size=1\times n_d^{(p)}$

bounds.phase(p).integral.lower = int_lower;
bounds.phase(p).integral.upper = int_upper;

• Parameter Bounds (if defined): $size=1\times n_s $

bounds.parameter.lower = para_lower;
bounds.parameter.upper = para_upper;

• Eventgroup Bounds (if defined): $size=1\times n_b$

bounds.eventgroup(g).lower = event_lower;
bounds.eventgroup(g).upper = event_upper;

Guess:

Each variable declared should have its guessing value, or GPOPS-II could not iterate based on the guess (except for eventgroup and path).

Attention: For state and control, each row of guess should correspond to the time points you provide: \(X_{guess-i}=X_{guess}(t_i)\)

• Time Guess: $size=M^{(p)}\times1$

guess.phase(p).time = t_guess;

• State Guess: $size=M^{(p)}\times n_y^{(p)}$

guess.phase(p).state = state_guess;

• Control Guess: $size=M^{(p)}\times n_u^{(p)}$

guess.phase(p).control = control_guess;

• Integral Guess: $size=1\times n_d$

guess.integral = integral_guess;

• Parameter Guess: $size=1\times n_s$

guess.parameter = parameter_guess;

Auxdata

auxdata includes all the constants you are going to use during the solution, making it easier for users to import the constants in continuous function and endpoint function without repeatedly defining the same constant.

Optional:

Here we just list some frequently-used parameters.

Derivatives and Cofig:

Field	Possible Values	Default
setup.derivatives.supplier	‘sparseBD’, ‘sparseCD’, ‘sparseFD’	‘sparseFD’
setup.derivatives.derivativelevel	‘first’, ‘second’	‘first’
setup.derivatives.dependencies	‘full’, ‘sparse’, ‘sparseNaN’	‘sparseNaN’
setup.method	‘RPM-Differentiation’, ‘RPM-Integration’	’RPM-Differentiation’

Scales:

Field	Possible Values	Default
setup.scales.method	‘none’, ‘auto-bounds’, ‘auto-guess’	‘none’

Mesh:

Field	Possible Values	Default
setup.mesh.method	‘hp-PattersonRao’, ‘hp-DarbyRao’, ’hp-LiuRao’, ’hp-LiuRao-Legendre’	‘hp-PattersonRao’
setup.mesh.tolerance	$[0, 1]$	$10^{-3}$
setup.mesh.maxiterations	Integer $\geq 0$	$10$
setup.mesh.colpointsmin	Interger $\geq 2$	3
setup.mesh.colpointsmax	Integer $\geq$ colpoints	10
setup.mesh.phase(p).fraction	Length $M\geq 1$ of Positive Numbers $>0$ and $<1$ that Sum to Unity	0.1*ones(1,10)
setup.mesh.phase(p).colpoints	Length $M \geq 1$ of Positive Integers $> 1$ and $< 10$	4*ones(1,10)

NLP:

Field	Possible Values	Default
setup.nlp.solver	‘snopt’, ‘ipopt’	‘ipopt’
setup.nlp.ipoptoptions.linear solver	‘mumps’ or ‘ma57’	‘mumps’
setup.nlp.ipoptoptions.tolerance	$>0$	$10^{-7}$
setup.nlp.ipoptoptions.maxiterations	Integer $>0$	2000
setup.nlp.snoptoptions.tolerance	$>0$	$10^{-6}$
setup.nlp.snoptoptions.maxiterations	Integer $>0$	2000

Output

The most important component of structure output is output.result, which includes:

• output.result.objective: value of your objective

• output.result.solution: optimal solution, includes:

  • output.result.solution.phase(p).time: time sequence

  • output.result.solution.phase(p).state: state matrix, consists of row vectors

  • output.result.solution.phase(p).costate: approximated costate values for checking your results

  • output.result.solution.phase(p).control: optimal control, cosists of row vectors

Possible Errors

Usually, we use SNOPT instead of IPOPT if we can get one, because the performance of SNOPT is better (but IPOPT is free). Normally, you should get the message below if GPOPS-II solves your problem correctly:

SNOPTA EXIT   0 -- finished successfully
SNOPTA INFO   1 -- optimality conditions satisfied

But you may receive different messages if the problem were not solved successfully. For detailed information, you may go to the userguide of SNOPT. Here we just list error messages easiest to meet during your programming.

Infeasible Solution

If you get this message, it means that the solution inexists because the conditions you provided for problem solving are impossible to satisfy:

SNOPTA EXIT -- 10 The problem appears to be infeasible
SNOPTA INFO -- 11 infeasible linear constraints
SNOPTA INFO -- 12 infeasible linear equalities

You should check the conditions you provided and fix them.

Resource Limit Error

If you get this message, it means that the solution exists but it cannot converge to the correct result within iterations:

SNOPTA EXIT -- 30 Resource limit error
SNOPTA INFO -- 31 iteration limit reached
SNOPTA INFO -- 32 major iteration limit reached

You should try to increase the time of iterations.

Numeraical Difficulties

If you get this message, it means that SNOPT cannot improve the result during iterations:

SNOPTA EXIT -- 40 Terminated after numerical difficulties
SNOPTA INFO -- 41 current point cannot be improved
SNOPTA INFO -- 42 singular basis
SNOPTA INFO -- 43 cannot satisfy the general constraints
SNOPTA INFO -- 44 ill-conditioned null-space basis

Firstly, you should check the dynamic equations you defined, because mostly the error is caused by the wrong dynamic equations. Then you should try to reduce the numerical illness of the problem.