A Comparative Study of Sensitivity Computations in ESDIRK-Based Optimal Control Problems

In this paper, we compare the impact of iterated and direct approaches to sensitivity computation in fixed-step explicit singly diagonally-implicit Runge-Kutta (ESDIRK) methods when applied to optimal control problems (OCPs). We use the principle of internal numerical differentiation (IND) strictly for the iterated approach, i.e., reusing the iteration matrix factorizations, the number of Newton-type iterations, and Newton iterates, to compute the sensitivities. The direct method computes the sensitivities without using the Newton schemes. We compare the impact of the iterated and direct sensitivity computations in OCPs for the quadruple tank system. We benchmark the iterated and direct approaches with a base case. This base case is an OCP that applies an ESDIRK method that refactorizes the iteration matrix in every Newton iteration and uses a direct approach for sensitivity computations. In these OCPs, we vary the number of integration steps between control intervals and we evaluate the performance based on the number of SQP and QPs iterations, KKT violations, and the total number of function evaluations, Jacobian updates, and iteration matrix factorizations. The results indicate that the iterated approach outperforms the direct approach but yields similar performance to the base case.


I. INTRODUCTION
Efficient computation of integrator sensitivities is essential for gradient-based optimization algorithms applied to systems described by ordinary differential equations (ODEs).Optimal control problems (OCPs) that are discretized using a multiple shooting approach belong to this category.They require the numerical solution of the system dynamics and the corresponding integrator sensitivities.One-step methods such as Runge-Kutta (RK) integration schemes have proven to be particularly effective when dealing with frequent discontinuities inherent in OCPs that implement zero-order hold parameterizations for the inputs [1].
For systems described by stiff ODEs or differentialalgebraic equations, a multiple shooting-based OCP algorithm may require implicit Runge-Kutta (IRK) methods.Explicit singly diagonally-implicit Runge-Kutta (ESDIRK) integrators are such types of methods and we characterize them by a diagonal structure with identical coefficients on the diagonal and with the first stage being explicit.We can implement these methods such that they reuse the matrix factorizations for all Newton-type iterations within a single integration step.We refer to this matrix as the iteration matrix, and its reuse leads to a computationally efficient method.ESDIRK-based multiple shooting OCPs have been Jørgensen (E-mail: jbjo@dtu.dk).jbjo@dtu.dkapplied to the continuous stirred tank reactor (CSTR) and the quadruple tank system (QTS) in [2].
There exist numerous approaches for sensitivity computations such as the Staggered direct method, the Simultaneous corrector method, and the Staggered corrector method [1].The Staggered direct method computes the sensitivity solutions directly from the iteration matrix after the system of nonlinear equations has converged.However, it requires frequent iteration matrix factorization to produce accurate integrator sensitivity information.The Simultaneous corrector method solves a larger combined state and sensitivity equations system.This approach utilizes a special Jacobian structure and allows the factored iteration matrix to be reused for multiple steps while keeping a high accuracy of the sensitivities.Finally, the Staggered corrector method applies a separate Newton-type scheme to obtain sensitivities, after converging the system of nonlinear equations [3].
The principle of internal numerical differentiation (IND) provides a methodology for developing sensitivities that are closely related to the integrator code [4], [5].The principle of IND involves computing the sensitivities by directly differentiating the discretization scheme generated adaptively from the integrator.A strict implementation of IND for IRK methods results in reusing not only the step-sizes of the integrator but also the iteration matrix factorizations, number of Newton-type iterations, and the sequence of Newton iterates to compute sensitivities.Such sensitivity computations are applied to backward differentiation formula (BDF) methods and [5] refers to it as the iterated IND.This iterated IND method is similar to the Staggered corrector method.
In an alternative IND approach to sensitivity computations, we may assume that we solve the system of nonlinear equations exactly at each integration step.This allows for direct computation of the sensitivities based on the iteration matrix factorizations without using the Newton-type scheme.This is referred to as the direct IND approach in [5] and is identical to the Staggered direct method.
Different sensitivity computations for ESDIRK methods have already been shown, e.g., [1] applies the direct approach for sensitivity computation for ESDIRK34 and suggests reusing the Jacobians for all the internal stages for one integrator iteration.The Staggered corrector approach is applied to ESDIRK12, ESDIRK23, and ESDIRK34 methods in [2].The direct method, unlike the iterated method, avoids using the Newton-type scheme for sensitivity computation, making it a more computationally efficient option.However, this efficiency comes at the cost of obtaining only approximate sensitivities.A computational overview of the iterated and direct approaches is shown in [5] and [6] compares various computational aspects such as LU factorizations and Jacobian updates for iterated and direct IND methods.However, there has been no clear demonstration of the impact of these sensitivities on gradient-based optimization.
In this paper, we compare the computational performance of using either an iterated or a direct approach for sensitivity computations in ESDIRK-based multiple shooting OCPs.We do this using a basic implementation of a sequential quadratic program (SQP) solver and compare the performances in terms of the number of SQP and QP iterations, function evaluations, Jacobian updates, and the number of iteration matrix factorizations.Our comparison is based on repeatedly solving OCPs for a model of the QTS with varying numbers of integration steps between control intervals.Additionally, we evaluate both approaches against a base case, where a direct approach to sensitivity computation is employed for an ESDIRK method using an exact Newton scheme, i.e., the iteration matrix being refactorized in all Newton-type iterations.The study shows that the iterated method is similar to the base case in terms of computational performance while avoiding some Jacobian updates and refactorizations in the integrator.We also show that the direct approach may only achieve convergence in the optimizer if we perform many integration steps for a small sampling time of the controller.
The rest of the paper is organized as follows.Section II introduces the ESDIRK integration method for ODEs.Section III shows how we compute sensitivities for the ES-DIRK methods using the iterated and direct IND approaches.In Section IV we present the optimal control problem.Section V shows numerical examples using the ESDIRKbased multiple shooting OCP with either a direct or iterated IND approach and the base case applied to the QTS.Finally, Section VI presents conclusions.

II. THE ESDIRK METHODS FOR ODES
We consider the initial value problem (IVP) where x(t) ∈ R nx is the state vector, u(t) ∈ R nu is the vector of inputs, and d(t) ∈ R n d is the vector of disturbances.
The ESDIRK integration scheme for solving (1) is where T i and X i are the internal nodes and stages at iteration k for i = 1, . . ., s, with s being the number of stages, and ESDIRK12 x k and x k+1 are the steps computed at t k and t k+1 = t k + h, respectively.We may compute the integration step-size h based on the local error estimate e k+1 = x k+1 − xk+1 , using the embedded method in (2d).However, for simplicity, we only consider a fixed integration step-size selection in this paper.Here, h is defined as h = N , where N is the number of integration steps.
Table I presents the Butcher tableaus for ESDIRK12, ESDIRK23, and ESDIRK34.We apply the numerical values for these tableaus from [7].For the remaining part of the paper, we apply the simplified notation f For each integration step, we solve the systems of nonlinear equations for i = 2, . . ., s, with We do this by applying the inexact Newton method for l = 0, . . ., w i,k − 1, with w i,k being the number of iterations required for stage i at iteration k to satisfy a chosen convergence criteria.M k is the iteration matrix at iteration k defined as where I is the identity matrix.We solve (5) using the LU factorizations of M k and we update these in every integration step.We choose the convergence criteria for the inexact Newton method as < τ, (7) with abs and rel being absolute and relative tolerances, respectively, and τ = 0.1, for i = 2, . . ., s [2].As the ESDIRK methods are stiffly accurate, i.e., a sj = b j , we avoid the computations in (2c) and obtain the next step directly as A. Computation of initial guesses for Newton iterations We apply stage value predictors (SVPs) to generate initial guesses of the Newton iteration schemes.The SVPs use information from the previously converged step, x k−1 , and previous converged stages, Xj = Xj [w j,k−1 −1] for j = 2, . . .s, to construct guesses on the form For fixed integration step-size r = 1.We construct these predictors using the order conditions in [8].The SVPs for ESDIRK12, (α 12 , β 12 ), and for ESDIRK23, (α 23 , β 23 ) are .
We don't show the SVPs for ESDIRK34 due to the lack of space.For the first integration step (k = 0), we use the trivial predictor, i.e., X [0] i = x k .We note, that due to (8), (9) applies information of

A. Iterated IND
In the iterated IND, we differentiate the adaptively generated discretization scheme, including the operations in the Newton-type scheme.The adaptive components for the fixed step-size ESDIRK methods are the LU factorizations of the iteration matrices, M k , the number of Newton-type iterations for all stages, w i,k , and the sequence of Newton iterates, X [l] i for l = 0, . . ., w i,k − 1 and for i = 2, . . ., s.We express the implementation of the ESDIRK methods as a combination of the elementary operations where (11a) represent the SVPs for stage i in ( 9) with X = X[w 2,k−1 −1] 2 ; . . .; x k , (11b) represent the Newton-type iterations in (5), and (11c) represent the mapping of the final stage to the next step in (8).
1) State sensitivity: We construct the state sensitivities by computing the partial derivatives of (11) with respect to the initial state, x 0 , i.e., for i = 2, . . ., s, with (12b) being We compute the partial derivative of the residual function as with where for j = 1, . . ., i − 1 in (15b) are the stages converged after w j,k -number of Newton-type iterations at integrator iteration k.For (12a), we compute the state sensitivities of the SVPs as for i = 2, . . ., s and for (12c) we write 2) Input sensitivity: We construct sensitivities with respect to the input as for i = 2, . . ., s.We compute (18b) as where the partial derivative of the residual function is ∂Ri ∂X [l] i and ∂Ri ∂ψi are the same as in (15a) and we calculate ∂ψi ∂u and ∂Ri ∂u as The contribution to the sensitivities for the SVPs in (18a) is for i = 2, . . ., s.For (18c) the input sensitivity at the next iteration is directly obtained as We initialize the state and input sensitivities at integration iteration k = 0 as 3) Implementations details: The ESDIRK integration schemes are implemented such that at least one Newton-type iteration is performed for each stage.This ensures that all elementary operations in (11) are included when computing the sensitivities.

B. Direct IND
In the direct IND approach, we assume we solve (3) for i = 2, . . ., s, exactly.Using the implicit function theorem, this enables the direct computation of state and input sensitivities for (2) as where J i = ∂Ri(Xi) ∂Xi , and ∂ψi ∂x0 and ∂ψi ∂u are described in (15b) and (21a), respectively.We apply the approximation in (6) to reuse the factorizations of ∂Ri(Xi) ∂Xi in (25) Identically to the iterated approach, we obtained ∂x k+1 ∂x0 and ∂x k+1 ∂u as ( 17) and ( 23), respectively, and we apply the same initial values as in (24).We refer to (26) as the direct IND approach for the ESDIRK methods (2).

C. Comments on iteration matrix refactorization strategies
Some ODE/DAE solvers improve computational efficiency by reusing iteration matrix factorizations through multiple consecutive integration steps when solving (1).Refactorization strategies similar to the strategies in [9] have successfully been implemented for ESDIRK methods, as demonstrated in [10].However, integrators that apply such strategies in gradient-based optimization may experience convergence issues, especially those using the direct IND method.This is because factorizations are not only used for solving nonlinear equations but also for directly obtaining the sensitivities.In contrast, as suggested in [5], we can equip integrators using the iterated approach with such refactorization strategies.This makes the iterated IND approach an attractive choice for IRK methods, where these strategies significantly enhance efficiency.However, it is important to note that the iterated method still requires frequent Jacobian updates, which can potentially eliminate advantages gained from such a strategy.

IV. OPTIMAL CONTROL PROBLEM
We formulate an OCP as a weighted least-squares problem that minimizes the deviation between outputs and setpoints and penalizes the input rate of movement.We formulate this OCP from t 0 to t f = t 0 + T , where T is the prediction horizon.We apply a zero-order hold parameterization of the inputs, i.e., u(t) = u n for t n ≤ t ≤ t n+1 = t n + T s with n = 0, . . ., N c − 1, and T s = T /N c being the interval between inputs (sampling time).We define the input rate of movement as ∆u n = u n − u n−1 and we also impose input bound constraints.This constitutes the OCP with N = 0, . . ., N c − 1, and with the objectives ϕ z penalizes the difference between the output z(t) and the setpoint z(t) and ϕ ∆u penalizes the input rate of movement.Q z and Q∆u = Q ∆u /T s represent the weight matrices for ϕ z and ϕ ∆u , respectively.We transcribe the OCP ( 27)-(28) into a nonlinear programming (NLP) problem using a multiple shooting approach.We solve an IVP between each input interval using the ESDIRK methods with N steps.To solve the NLP, we apply an SQP solver with a backtracking line-search algorithm and BFGS to update the Lagrangian Hessian, developed in [11].

V. NUMERICAL EXPERIMENTS
We compare the computational performance of the iterated and direct approaches in ESDIRK-based OCP.We do this by repeatedly solving OCPs with a varying number of integration steps between control intervals for the model of the QTS presented in [12].We compare these two approaches to a base case.This base case applies an ESDIRK method with an exact Newton method, i.e., refactorizing the iteration matrix at every Newton-type iteration.This base case applies a direct approach which due to the frequent refactorizations, achieves accurate sensitivity information.Table II presents the model parameters for the QTS.We use the initial conditions x(t 0 ) = [7602.7,11404.0,1000.0,1000.0]′ g and apply the constant vector of disturbances d(t) = [0, 0, 100, 100] ′ cm 3 /s for all experiments.The outputs of the QTS are the water levels in the two bottom tanks, i.e., z(t) = 1 ρA1 x 1 (t), For the OCP, we choose the sampling time T s = 10 s.We use the weight matrices Q z = diag( [10,10]) and Q ∆u = diag([0.1,0.1]), and the number of steps in the horizon N c = 40.The bounds are between 0 cm 3 /s and 500 cm 3 /s for both inputs.For the experiments, we construct the time-varying setpoints z(t) = 20, 30 ′ cm for 0 ≤ t < T /2 and z(t) = 30, 20 ′ cm for T /2 ≤ t ≤ T .We initialize all states and inputs for the SQP as a vector with all entries set to 300.For the SQP algorithm, we choose the SQP tolerances as 10 −3 , the QP tolerance 10 −8 , the SQP step length tolerance 10 −8 , and abs and rel in (7) both as 10 −8 .Fig. 1 presents one solution of the OCP using the ESDIRK12 with 10 integrator steps between control intervals and iterated IND sensitivity computations.

A. Varying integration step-size between control intervals
We solve the OCPs repeatedly while increasing the number of integration steps between control intervals.We apply the sequence N = [5,10,15,20,25,30,35,40,45,50]    integration steps between controls.We also do this with the base case.To measure the performance of the three cases, we store the number of SQP iterations, the KKT violations, and the number of QP iterations in the SQP solver.We also store the total number of function and Jacobian evaluations in the ESDIRK integration scheme and LU factorizations of the iteration matrix.Fig. 2, 3, and 4 present the results using ESDIRK12, ESDIRK23, and ESDIRK34, respectively.The dotted black lines are the SQP tolerance representing the allowed KKT violations.The crosses in these figures represent experiments that do not converge.The non-convergence is due to the step length, computed using the backtracking line-search algorithm, exceeding the step tolerance.

B. Low SQP tolerance and short control intervals
We also demonstrate the performance when we lower the tolerances for the SQP to 10 −6 , the QP to 10 −10 , and abs and rel to 10 −10 .We also choose a smaller control interval at T s = 2.0 s.We keep the same step tolerance.Fig. 5 shows the results.Crosses represent experiments that do not converge to the SQP tolerance of 10 −6 .tolerance of 10 −3 in all experiments.The direct approach only achieves convergence for ESDIRK23 when using 10 integration steps between control intervals.Additionally, the iterated approach performs similarly the base case while requiring fewer total Jacobian updates and LU factorizations.The base case requires more LU factorizations compared to the iterated approach as it refactorizes the iteration matrix at every Newton iteration.Both the base case and the iterated approach update the Jacobians once per Newton iteration, but the base case requires one additional Jacobian update before every Newton scheme.It requires this, as it computes the Jacobian for the explicit state and also for constructing the iteration matrix based on the SVPs before every Newton scheme.In contrast, the iterated approach uses the same Jacobians for both the explicit step and the iteration matrix.

C. Comments on the numerical experiments
Fig. 5 demonstrates that the direct method can achieve convergence for ESDIRK23 and ESDIRK34 as long as we maintain an SQP tolerance of 10 −3 and a short control interval.However, it still cannot reach the low tolerance requirement at 10 −6 .We assume that the direct method fails to converge due to insufficient gradient information at SQP iterations close to convergence.We assume this because all sensitivity methods in SQP result in full step lengths in the initial phases, while the line-search algorithm drastically reduces this length near the solution when using the direct method.This reduction of the step length suggests a lack of accurate gradient information, which hinders the SQP method's ability to navigate to the optimal solution.VI.CONCLUSIONS In this paper, we compare the iterated and direct approaches to sensitivity computation for ESDIRK-based OCPs using multiple shooting discretization.The iterated approach strictly applies the principle of internal numerical differentiation, and the direct method approximates sensitivities without employing the Newton-type schemes from the ES-DIRK methods.We evaluate both methods in terms of some computational performance metrics by repeatedly solving ESDIRK-based OCPs for the QTS with a varying number of integration steps between control intervals.The computational performance metrics are the number of SQP and QP iterations, KKT violations, the total number of function evaluations, Jacobian updates, and iteration matrix factorizations.We compare these approaches to an ESDIRK method that applies an exact Newton scheme and employs a direct sensitivity computation approach.We refer to this as the base case and the study shows that the iterated method performs similarly to the base case but requires fewer Jacobian updates and matrix factorizations.Furthermore, the iterated approach outperforms the direct method for longer sampling times.
A. H. D. Christensen and J. B. Jørgensen are with the Department of Applied Mathematics and Computer Science, Technical University of Denmark, DK-2800 Kgs.Lyngby, Denmark.Corresponding author: J. B. of

Fig. 2 -
Fig.2-4demonstrate that ESDIRK methods using the iterated sensitivity approach converge to the desired SQP

Fig. 5 .
Fig. 5. Comparison of ESDIRK-based OCP algorithms with direct and iterated sensitivity computations.N = 10 for all methods.

TABLE I BUTCHER
TABLEAU'S FOR SOME ESDIRK METHODS.