TACO: a toolkit for AMPL control optimization

Abstract

We describe a set of extensions to the AMPL modeling language to conveniently model mixed-integer optimal control problems for ODE or DAE dynamic processes. These extensions are realized as AMPL user functions and suffixes and do not require intrusive changes to the AMPL language standard or implementation itself. We describe and provide TACO, a Toolkit for AMPL Control Optimization that reads AMPL stub.nl files and detects the structure of the optimal control problem. This toolkit is designed to facilitate the coupling of existing optimal control software packages to AMPL. We discuss requirements, capabilities, and the current implementation. Using the example of the multiple shooting code for optimal control MUSCOD-II, a direct and simultaneous method for DAE-constrained optimal control, we demonstrate how the problem information provided by the TACO toolkit is interfaced to the solver. In addition, we show how the MS-MINTOC algorithm for mixed-integer optimal control can be used to efficiently solve mixed-integer optimal control problems modeled in AMPL. We use the AMPL extensions to model three control problem examples and we discuss how those extensions affect the representation of optimal control problems. Solutions to these problems are obtained by using MUSCOD-II and MS-MINTOC inside the AMPL environment. A collection of further AMPL control models is provided on the web site http://mintoc.de. MUSCOD-II and MS-MINTOC have been made available on the NEOS Server for Optimization, using the TACO toolkit to enable input of AMPL models.

Appendices

Appendices

The following sections contain supplementary material intended to guide software developers interested in using the presented TACO toolkit to interface their optimal control software packages with AMPL. Appendix A explains the most important data structures that hold information about the mapping from the optimal control point of view to the AMPL view. Appendix B lists functions available to optimal control software package developers for evaluating AMPL functions. Appendix C explains error codes emitted by the TACO toolkit, and mentions possible remedies.

Appendix A: TACO data structures exposed to optimal control problem solvers

This sections lists TACO data structures exposed to developers of software packages for solving optimal control problems. We discuss several snippets taken from the header file ocp_frontend.h, which should be consulted for additional details.

1.1 A.1: Data structures mapping from AMPL to optimal control

The optimal control frontend provides a collection of fields that hold the optimal control problem interpretation of every AMPL variable passed to the solver. They are laid out as follows:

The field vartypes gives the OCP variable type of an AMPL variable, that is, whether the variable is $t$, $T$, a component of vector $p$, or a component of one of the vector trajectories $x$, $z$, or $u$ and $w$. For the case of it being a vector component, the field varindex holds the index into the OCP variable or trajectory vector. For AMPL constraints and objectives, no mapping information from AMPL to the OCP perspective is provided.

1.2 A.2: Data structures mapping from optimal control to AMPL

The AMPL optimal control frontend provides a collection of structures holding the AMPL perspective for every component of an optimal control problem according to problem class (MIOCP). Starting with information about problem dimensions, the following variables are provided and should be self-explanatory.

Information about both the independent time variable and the end-time variable is held in a structured variable named endtime of the following layout.

For a fixed-endtime scenario, endtime.index is $-1$, the field endtime.init holds the fixed end time. For a variable end-time scenario, endtime.index is nonnegative and the field endtime.init holds the initial guess for the free end time if available, and is set to $1.0$ otherwise.

The following structured variable xstates holds information about differential state trajectory variables, and associated right-hand side functions.

Here, the field fixed holds the AMPL index of the initial value constraint for this ODE state. It is $-1$ if the ODE state’s initial value is free. The initializer init provides two values for linear interpolation (see suffix .interp_to).

The field ffcn_index holds the AMPL constraint index of the right-hand side function associated with a differential state. Even though we currently support explicit ODEs only, AMPL-internal rearrangement of constraint expressions may cause a negative sign on the diff() call. Hence the field rhs_factor is introduced to compensate for AMPL-internal representations of the form $-\dot{x}(t) = f(t,x(t),\ldots )$.

Similar to differential states, the structured variable zstates holds information about algebraic state trajectories.

It is important to keep in mind that DAE constraints are not associated with algebraic state trajectory variables by a one-to-one mapping. We merely keep both in the same array for simplicity, since their numbers must match.

Information about integer and continuous control trajectories is kept in a structured variable named controls with the following layout.

The field type denotes the (solver-dependent) discretization type to be applied to this control trajectory. The field integer is set to $1$ if the control is a binary or integer control and to $0$ if it is a continuous control. For piecewise linear and piecewise cubic discretization types, additional information about slope limits and initial guesses is provided.

In addition, model parameters information is provided in a structured variable named params, with layout as follows. This concludes AMPL variables information.

Information about Mayer-type, Lagrange-type, and integral least-squares-type objective functions (1a) is found in structured variables named mayer, lagrange, and clsq. Their layout is presented below. In this representation, the optimal control problem objective may consist of multiple parts that are summed up, see Eq. (1a), although only one part of each objective type is allowed. Problem solvers may want to support feasibility problems that do not have an objective function at all.

dpnd is a bit field with bit $k$ ($k\ge 0$) set if and only if the objective function’s AMPL expression depends on an AMPL variable with vartypes entry set to value $k$ (see enum vartype_t). This allows for quick dependency checks that may save runtime, for example, in derivative approximation. The field maximize is set to $1$ if the objective function is to be maximized and to $0$ if it is to be minimized. Note that least-squares functions are recognized only if they are to be minimized. The fields nparts and parts hold information about the AMPL expression DAGs associated with the individual least-squares residual expressions of a least-squares objective.

Path constraints (1d) and coupled constraints (1e) information is held in structured variables named pathcon, cpcon_eq, and cpcon_in with the following layout.

Path constraints (1d) always are inequality constraints. For coupled constraints (1e), equality and inequality constraints are stored in separate arrays. For two-sided inequality constraints $l\le c(x)\le u$, the field side indicates which side of the constraint should be evaluated.

Information about decoupled point constraints (1f) and point least-squares objectives (1a) is stored in a structured variable named grid. For each objective or constraint evaluation time, a grid node is introduced and holds information about the associated objective or constraint. Grid nodes are guaranteed to be unique, in other words, no two nodes share the same time point, and they are sorted in ascending order. The grid is guaranteed to contain at least two nodes: the first grid node will always be at time $0$, and the last grid node will always be at time $t_\text{ f }$.

For equality and inequality point constraints (1f) on a grid node, the fields n_eq and n_in hold the dimensions and the fields con_eq and con_in the constraint information, respectively. For coupled constraints (1e), we do not store pointers to constraint information structures, but rather indices into the global lists of coupled constraints cpcon_eq and cpcon_in. For solvers requiring linear separability of coupled constraints, this layout eases the setup of the coupled constraints’ block structure. The field lsq holds information about the point least-squares objective contribution in a node; again index $-1$ indicates that no point least-squares objective is present.

Appendix B: TACO functions exposed to optimal control problem solvers

This sections lists TACO functions exposed to developers of software packages for solving optimal control problems.

Reading of the AMPL model. For reading and verifying the AMPL mode as well as creation of management of the database, the optimal control frontend provides the following functions:

void allocate_mappings()

• allocates memory for the database to be created, prior to reading the stub.nl file provided by AMPL.

void read_mappings()

• calls the AMPL solver library to read the stub.nl file. Afterwards, the DAGs of all AMPL objectives and constraints are examined for appearance of variables, calls to user functions, and nonsmooth operators. The role of AMPL variables, constraints, and objectives in class (MIOCP) is determined, and the database is filled with appropriate information.

void free_mappings()

• frees memory allocated for the database after the optimal control problem has been solved.

A developer wishing to interface an optimal control problem solver with AMPL is provided with a number of functions that infer the optimal control problem’s structure from a stub.nl file, as shown in Fig. 11.

Evaluating problem functions. Once the database has been successfully filled with the optimal control problem representation of the AMPL model, an appropriate optimal control problem solver can be called. This solver will have to evaluate problem functions. For conveniently doing so, the optimal control frontend provides a set of functions that wrap around the conival and objival functions provided by AMPL’s solver library:

double evaluate_ode_rhs(i,y)

• evaluates the ODE right-hand side $\dot{x}_i(t)=f_i(\cdot )$, see Eq. (1b), for a differential state $x_i(t)$, $0 \le i < n_\text{ x }$. The vector y here is the plain vector of current AMPL variable values.

double evaluate_dae_rhs(i,y)

• evaluates the residual of DAE constraint $g_i(\cdot )$, $0\le i < n_\text{ z }$, see Eq. (1c).

double evaluate_scaled_mayer(y)

• evaluates and scales the value of the Mayer objective term, see Eq. (1a), if a Mayer term exists.

double evaluate_scaled_lagrange(y)

• evaluates and scales the value of the Lagrange objective term’s integrand, see Eq. (1a), if a Lagrange objective exists.

void evaluate_scaled_clsq(y,res)

• evaluates and scales the value of the Lagrange objective term’s integrand, see Eq. (1a), if the integrand has least-squares structure and if a Lagrange objective exists. Upon return the vector res holds the least-squares residuals, that is, the individual values of each squared summand.

void evaluate_scaled_node_lsq(n,y,res)

• evaluates one summand of a point least-squares objective term, see Eq. (1a), if one exists. The integer n specified the number $0<n<N_\text{ lsq }$ of the measurement point to be evaluated.

double evaluate_path_con(i,y)

• evaluates and scales the residual of an inequality path constraint, see Eq. (1d),

double evaluate_coupled_eq_con(i,y)

• evaluates and scales the residual of a coupled equality constraint, see Eq. (1e),

double evaluate_coupled_ineq_con(i,y)

• evaluates and scales the residual of a coupled inequality constraint, see Eq. (1e),

double evaluate_point_eq_con(n,i,y)

• evaluates and scales the residual of an equality point constraint, see Eq. (1f),

double evaluate_point_ineq_con(n,i,y)

• evaluates and scales the residual of an inequality point constraint, see Eq. (1f),

Writing solutions from a solver. To write a solution file of the proposed format, the optimal control frontend provides three function calls to be used by a backend interfacing an optimal control problem solver with AMPL:

void alloc_solution(solution, gridlength)

• allocates a solution structure given a discretization grid length.

void write_solution(solution, filename)

• writes a solution file according to the new format and fills it with the solution data found in the solution structure. The solver backend should have queried the solver itself for this solution and copied it over to this solution structure.

void free_solution(solution)

• free a previously allocated solution structure after it has been written to a solution file.

Appendix C: TACO toolkit error messages

TACO performs a range of consistency tests during the analysis and verification of the optimal control problem structure. If option solvername _auxfiles is set to contain the characters “r” and “c”, error messages will show the name of the offending variable, objective, or constraint. Codes 1 to 5 are reserved for internal errors; the remaining error messages are:

• 6. “AMPL error reading stub.nl file filename” AMPL could not read the stub.nl file.

• 7. “AMPL error getting solver options” AMPL could not read the solver options. Make sure your solver options string contains known keywords and valid values only.

• 8. “I/O error writing solution file filename” AMPL could not write the solution file. Make sure that the current working directory is writeable and that the volume is not full.

• 9. “AMPL network variables are not supported” The AMPL model contains network variables. Network modeling is currently not supported.

• 10. “AMPL piecewise linear terms are not supported” The AMPL model contains piecewise linear terms. Nonsmooth models are currently not supported.

• 11. “diff() must be called exactly once for each differential state variable” You called diff() more than once for the same variable. Every differential state variable is associated with exactly one right hand side function. Most likely, this is a typo in your model.

• 12. “Logical operator operator not permitted” The shown constraint or objective make use of a logical operator. Nonsmooth models are currently not supported.

• 13. “Relational operator operator not permitted” The shown constraint or objective make use of a relational operator. Nonsmooth models are currently not supported.

• 14. “Nonsmooth operator operator not permitted” The shown constraint or objective make use of a nonsmooth operator. Nonsmooth models are currently not supported.

• 15. “ODE constraints must be equality constraints” You called diff() in an inequality constraint. All ODE constraints must be equality constraints, though.

• 16. “Cannot determine type of objective functions” The shown objective cannot be determined to be of Mayer, Lagrange, integral least-squares, or point least-squares type. Please refer to Section 2.3 to learn about the structure of supported objective functions.

• 17. “Mayer objective function must not depend on controls” The indicated objective is of Mayer type, and hence must not depend on control trajectory variables.

• 18. “The number of DAE constraints does not match the apparent number of algebraic variab- les” The number of DAE constraints found in your model does not match the number of algebraic variables you have declared by letting suffix .type to value "dae". Make sure you’ve set the suffix of all algebraic state variables, and make sure all your DAE constraints are equality constraints.

• 19. “Invalid value for suffix .type” A variable, constraint, or objective has suffix .type set to an unrecognized value. This should not happen unless you’re using an incompatible OptimalControl.mod file.

• 20. “Missing independent time variable or end time variable” TACO could not determine which AMPL variables represent independent time or end time. Make sure that the independent time variable appears as second argument to all diff() calls, and that the end-time variable appears as second argument to integral() in Lagrange-type objectives, or to eval() in Mayer-type objectives or end-point constraints. If your AMPL model uses a time horizon of fixed length, make sure the end-time value shows up consistently in these places.

• 21. “Nesting a call to function name inside another is not allowed” You nested two calls to functions diff(), eval(), or integral(). This should never be necessary.

• 22. “As of now, diff() must be explicit. Try to write your ODE as diff(var, t) = ...;” The AMPL optimal control frontend currently supported explicit ODEs only, i.e. the permitted formulation is $\dot{x}(t)=f(t,x(t),\ldots )$ instead of the more general $0=f(t,x(t),\dot{x}(t),\ldots )$. Try to rewrite your model to use explicit ODEs only.

• 23. “First argument to eval() must not be a constant” The shown constraint or objective calls eval() with two constant arguments. This is not sensible. Either remove the calls to eval(), or make sure the expression to be evaluated references independent time or at least one state or control trajectory variable.

• 24. “Second argument to eval() must be >= 0.0, or the final time variable” The shown constraint or objective calls eval() with a time point that a) evaluates to a negative value, or b) is not a simple constant and not the end-time variable.

• 25. “First argument to integral() must not be a constant” An objective calls integral() with two constant arguments. This is not sensible. Either remove the call to integral(), or make sure the expression to be integrated references independent time or at least one state or control trajectory variable.

• 26. “Second argument to integral() must be the final time variable” The second argument of a call to integral() must be the final time variable, or the final time itself.

• 27. “Second argument to diff() must be the independent time variable” The second argument of a call to diff() must be the independent time variable.

• 28. “Function diff() expects differential state variable and independent time variable as arguments” The first argument to diff must be a differential state trajectory variable, and the second argument of a call to diff() must be the independent time variable.

• 29. “Calls to integral() not allowed in constraints” The function integral() was called inside a constraint expression. Calls to integral() are only allowed in a Lagrange-type or integral least-squares-type objective.

• 30. “Function call to name must enclose entire expression” Calls to eval() and integral() must enclose the entire expression to be evaluated or integrated.

• 31. “Only controls or parameters can be binary/integer variables” AMPL’s binary or integer restriction may be applied to control trajectory variables and model parameter variables only. Time as well as differential and algebraic state trajectory variables must be continuous.

• 32. “Invalid grid node position, must be within [0,1] (tf free) or [0,tf] (tf fixed)” A call to eval() used an invalid evaluation time. Allowed times are constants $0$ to $1$, being relative times on the time horizon $[0,t_\text{ f }]$ if $t_\text{ f }$ is free, absolute times $0$ to $t_\text{ f }$ if $t_\text{ f }$ is fixed, or the end-time variable $t_\text{ f }$ itself.

• 33. “Duplicate objective of name -type” At most one objective of Mayer-, Lagrange-, and least-squares type is allowed per problem.

• 34. “Calls to eval() not allowed in ODE constraints” The function eval() was called inside an ODE constraint expression. This is not allowed.

• 35. “Calls to diff() not allowed in objective functions” The function diff() was called inside an objective function expression. This is not allowed.