(101c) Efficient Global Solutions to Optimal Control Problems Via Polynomial Optimization and Sum-of-Squares Relaxations

Efficient solution methods for optimal control problems (OCPs) are useful for optimization-based state and parameter estimation, experimental design, and model-based control, among other tasks. Direct methods are a popular solution approach for OCPs wherein the original infinite-dimensional OCP is approximated as a finite-dimensional one via discretization of the time-varying functions [1]. However, direct methods only seek local optimality, rather than global optimality. The local optima attained by these optimization algorithms may be suboptimal with respect to the global optimum by a significant margin [2]. Alternatively, global optimization algorithms can be used. However, with any approach to global optimization, the worst-case complexity scales exponentially with the number of decision variables, which is high in OCPs even after the use of typical discretization methods.

To reduce the number of decision variables in OCPs, the parsimonious input parameterization approach has been proposed [3]. This approach generates a finite set of plausible arc sequences and describes each sequence by a small number of decision variables. Then, for a given arc sequence, one can compute the optimal values of these decision variables via numerical optimization, which can also be solved to global optimality. This contribution aims to extend the parsimonious input parameterization approach for efficient solution to OCPs to global optimality.

We consider the general class of OCPs formulated as the minimization of a terminal cost subject to terminal constraints, with piecewise-continuous inputs and states described by the differential equations of a continuous-time, nonlinear dynamical system, input bounds, and state constraints. The inputs that represent the solution to such OCPs are composed of several arcs. For each input, each arc can be input constraint-seeking, such that it is determined by the input bounds, state constraint-seeking, such that it is determined by an equality related to the state constraints, or sensitivity-seeking, such that it is determined by an equality that stems from the dynamics [3]. Hence, there is a finite number of arc types from which arc sequences can be formed. It follows that the number of plausible sequences is also finite.

The goal is to show how these OCPs can be solved efficiently to global optimality. The proposed approach for global optimality relies on determining: (i) when and how the optimal switching between arcs takes place for a given plausible arc sequence; and (ii) which sequence provides the optimal solution. Once question (i) is addressed for every plausible sequence, for example via parallel computing, it is trivial to answer question (ii). For a given plausible arc sequence, the input is defined by the following decision variables: the switching times between arcs, the final time, and the initial conditions of the sensitivity-seeking arcs. Then, addressing question (i) above consists in computing the optimal values of the decision variables for the given arc sequence.

For a given arc sequence, we describe the input in each interval by defining states for this interval and combining all the states into a vector of extended states. Then, one can eliminate input dependencies and reformulate the OCP in terms of the extended states and the new decision variables, which is convenient for numerical optimization since there are only a few decision variables. We aim to reformulate the OCP for each arc sequence as a polynomial optimization problem (POP) that is amenable to global optimization. This entails expressing each function in the cost and constraints as a polynomial function [4]. To this end, we compute each function and its first-order partial derivatives with respect to the decision variables.

An efficient approach to approximate each function as a polynomial function consists in using multivariate Hermite interpolation to obtain a polynomial of a certain degree that matches the value and the first-order partial derivatives at some sample points. Note that this requires no more than computing the extended states and the extended adjoint variables for every function that correspond to each point. The number of sample points is polynomial in the number of decision variables, which is typically small owing to the parsimonious input parameterization. This yields the polynomial representation of each function. Hence, the OCP for that arc sequence is reformulated as a POP.

The solution to the POP is not exactly the same as the solution to the OCP for the given arc sequence due to the fact that the polynomial functions are not perfect approximations of the functions in the cost and constraints. One can obtain explicit expressions for the approximations of (i) the difference between the global solution to the POP and the corresponding KKT point of the OCP and (ii) the difference between the globally optimal cost of the POP and the cost of the OCP at the KKT point as a function of the approximation errors. The exact KKT point of the OCP can be obtained via local optimization of the OCP with initial guess given by the global solution to the POP. Moreover, one can show that the difference between the costs of the OCP at this KKT point and at the globally optimal solution is bounded if the approximation errors are bounded.

Then, the concept of sum-of-squares (SOS) polynomials is applied to the OCP reformulated as POPs. In general, if a polynomial is strictly positive on a compact basic semi-algebraic set specified by some polynomials and the set satisfies some technical assumptions, then the former polynomial can be represented as a combination of SOS polynomials for some relaxation order that depends on the degree of the polynomials in the problem [5]. This amounts to satisfying a linear matrix inequality (LMI) for each SOS polynomial, which can be done via a convex semidefinite program (SDP).

This result is very useful for the problem of computing the global minimum of the polynomial that represents the cost of the OCP subject to the constraints that stem from the polynomials that represent the constraints of the OCP. Such a problem can be formulated as an SDP. Hence, if the number of decision variables and the degree of the polynomials in the problem are relatively small, the SDP can be solved efficiently since the relaxation order that provides a representation in terms of SOS polynomials is usually not very large. If this representation exists for some relaxation order, a certificate can be obtained upon convergence of the SDP [6]. Suppose that a global optimum is computed and certified for a small relaxation order. Since the complexity of SDPs is polynomial in their input size, it means that a global solution is computed and certified in polynomial time.

The plausible arc sequences under consideration are the sequences with a number of arcs no larger than some upper bound and without consecutive arcs of the same type. However, in many cases, it is reasonable to assume that one seeks a solution that includes a relatively small number of arcs with an upper bound between 3 and 5 since in practice one is not interested in chattering solutions even if they are globally optimal from a mathematical point of view. Hence, while the approach above can guarantee global optimality only if the globally optimal solution does not include more than the specified number of arcs, global optimality among all plausible arc sequences can be guaranteed in any case.

The proposed approach is applied to a simulation example that corresponds to a problem of production maximization in an acetoacetylation reaction system with the inlet flow rate as decision variable. For all the plausible arc sequences, it is possible to extract the unique solution to the POP from the solution to the SDP for the relaxation order 5 and certify the global optimality of the solution. One can determine the globally optimal solution with no more than 3 arcs, which only requires solving 6 problems in parallel in less than 40 s. If we had only used local optimization to compute a local solution to the OCP, we could have obtained a worse local solution and it would not have been possible to provide any guarantee that the local solution is in any way close to the globally optimal solution.

Future work will focus on extending this method so as to provide efficient global solutions in the presence of multiple inputs, state constraints, disturbances, and mismatch between the real dynamical system and its model.

[1] L. T. Biegler, A. M. Cervantes, and A. Wächter. Advances in simultaneous strategies for dynamic process optimization, Chem. Eng. Sci., 57(4):575-593, 2002.

[2] B. Houska and B. Chachuat. Branch-and-lift algorithm for deterministic global optimization in nonlinear optimal control, J. Optim. Theory Appl., 162(1):208-248, 2014.

[3] D. Rodrigues and D. Bonvin. On reducing the number of decision variables for dynamic optimization, Optim. Control Appl. Methods, 41:292-311, 2019.

[4] D. Henrion and J. B. Lasserre. GloptiPoly: Global optimization over polynomials with Matlab and SeDuMi, ACM Trans. Math. Softw., 29(2):165-194, 2003.

[5] J. B. Lasserre. Global optimization with polynomials and the problem of moments, SIAM J. Optim., 11:796-817, 2001.

[6] J. B. Lasserre. A semidefinite programming approach to the generalized problem of moments, Math. Program., 112:65-92, 2008.