Benchmarking solution methods for parameter identification in dynamical systems

Accurate models are crucial for simulating, optimizing, and controlling real‐world processes. Parameter identification—the task of estimating the unknown parameters of a dynamical system based on measurements—is challenging and there exist various methods to approach it. Integration‐based methods, such as shooting methods and full discretization, approximate the model output by numerically solving the dynamical system. Gradient matching methods, on the other hand, avoid solving the dynamical system and focus on minimizing the error between the measurement slope and the state derivatives instead. All approaches have advantages and disadvantages and a recipe for which method is best in a particular situation does not exist. In this paper, we present the results of a benchmark comparing single shooting, multiple shooting, full discretization, and gradient matching using a comprehensive database of test problems. We investigate if there is a notable difference in the performance of the various methods and if one approach is superior to the others. From the benchmark, we conclude that the integration‐based methods outperform gradient matching. While full discretization exhibits robustness, single and multiple shooting provide higher precision. Additionally, we observe that finding an optimal configuration of the methods across a diverse set of parameter identification problems remains challenging and often requires fine‐tuning.

guess for the states and model parameters and suffers from numerical instabilities [2].It also tends to find local minima [2].Therefore, many enhanced shooting methods were developed.The most prominent multiple shooting [3] splits the integration into subintervals, which makes it more robust regarding the numerical integration and often avoids local minima [4].Incremental single shooting [5] applies single shooting on a successively enlarged time horizon and thereby overcomes the deficiencies of the classical approach.Another integration-based parameter identification method is full discretization [6].It discretizes the state space and formulates an integration scheme as a constraint to ensure that the discretized states approximate a solution of the dynamical system.It is considered robust regarding local minima and instabilities in the integration [7].
In contrast to the integration-based parameter identification methods, gradient matching [8] takes an approach that avoids solving the dynamical system.Classically, it consists of two stages.The first stage is to approximate the slope of the measurements.A simple way of doing this is via finite differences or computing a spline approximation of the measurements and taking its derivative as the measurement slope [9].In the second stage, the model parameters are identified by minimizing the distance between the measurement slope and the state derivatives obtained from the dynamical system.The main challenge is to find a suitable approximation of the measurement slope.Many enhancements of the classical gradient matching address this problem.Poyton et al. [10] and Ramsay et al. [11] use the model error as a regularization in the first stage and iteratively improve the measurement slope approximation and the model parameters.Varziri et al. [12] go one step further and perform both stages simultaneously by solving a bi-objective optimization problem.Recent approaches use neural networks that incorporate knowledge about the model to learn the measurement slope [13].
Many articles discuss specific parameter identification methods and outline their advantages, but a comprehensive study is still widely missing.In general, it is unclear which method works well for a particular problem.This obstacle makes it considerably more challenging to solve parameter identification problems in practice.It is therefore desirable to get a better insight into the performance of various methods.To approach this, we benchmark four widely used parameter identification methods: single shooting, multiple shooting, full discretization, and gradient matching.We analyze their performance and systematically compare them using a comprehensive database of test problems from Schittkowski [14].We investigate if there is a notable difference in the performance of the methods and if one is superior to the others.Ultimately, the results can be used in the decision-making of choosing a method for a particular parameter identification problem.
Our objective is to determine the model parameters  and initial value  0 ∈ ℝ   that minimize the error between the model output and collected measurements.These measurements consist of states x ∈ ℝ   and corresponding controls ū ∈ ℝ   , observed at discrete time points   ∈ [ 0 ,   ].The task is formulated as the nonlinear parameter identification problem The expression (  ;  0 , ) indicates that the model solution explicitly depends on   , while its dependency on  0 and  is implied by the dynamical system.We use ũ to denote an approximation of ū, which can be evaluated at any time point within the interval [ 0 ,   ].This could for example be achieved via interpolation.The following parameter identification methods are described by means of Problem (1).For the numerical experiments, an extended problem formulation is used (see Section 4).

PARAMETER IDENTIFICATION METHODS
Solving parameter identification problems is challenging for several reasons.They have to cope with measurements that usually come from real-world observations and show noisy behavior.While these measurements must be fitted, it has to be ensured that the model output follows the governing dynamical system, which can make the problem highly nonlinear.In this paper, we focus on four parameter identification methods, each with advantages and disadvantages: single shooting, multiple shooting, full discretization, and gradient matching.

Single shooting
The most natural method to solve parameter identification problems is single shooting [1].We use an iterative optimization procedure to solve the problem and whenever the objective function needs to be evaluated, the dynamical system is solved numerically to compute the model solution (  ;  0 , ).The solution approach may be intuitive but has several disadvantages.First, the objective function of Problem ( 2) can be strongly nonlinear.As a result, the optimization problem is usually nonconvex and has many local minima [2].Finding a good solution to Problem (2) is, therefore, difficult.Second, the solution of the dynamical system is often sensitive to the initial value and the model parameters [2].Furthermore, numerically solving the dynamical system inevitably introduces errors.This further amplifies the effect of poor estimates for the initial states and model parameters.Due to numerical instabilities, the integration might even fail.This problem is especially noticeable when the dynamical system has to be solved over a large time horizon.The multiple shooting method addresses these challenges.

Multiple shooting
In multiple shooting [3], intermediate time points  0 =  1 < ⋯ <    <    , the shooting nodes, are introduced.They divide the time interval [ 0 ,   ] into disjunct subintervals  1 , … ,    .Instead of integrating over the whole time horizon, the integration is performed over these subintervals.Consequently, each shooting node is associated with an initial value [] 0 ∈ ℝ   ( = 1, … ,   ).A continuity condition is formulated as a constraint to ensure that the subtrajectories form a feasible solution.Summarizing, this yields the optimization problem where  [] ( +1 ; [] 0 , ) denotes the model solution on the -th shooting interval.Since multiple shooting handles smaller time horizons, it is more robust regarding the numerical integration [4].Another advantage is that during the optimization process, infeasible solutions are allowed.Consequently, the trajectory may be discontinuous in intermediate iterations.This freedom allows it to stay closer to the measurements and reduces the chance of terminating in a local minimum [4].

Full discretization
In full discretization [6], Problem (1) is discretized by introducing discretization points  0 =  1 < ⋯ <    =    and associated discrete state variables  [] ∈ ℝ   ( = 1, … ,   ).A discretization scheme is formulated as a constraint to ensure that the discrete state variables approximate a solution of the dynamical system.This may be a one-step method Φ or any other discretization scheme like collocation [7].We obtain the optimization problem subject to  [+1] =  [] + ℎΦ (   ,  +1 ,  [] ,  [+1] ,  ) ,  = 1, ⋯ ,   − 1, where x(  ;  [] ) denotes an interpolation of the discrete states that is evaluable at the measurement time points.Full discretization is robust regarding local minima and instabilities in the integration [7].A further advantage is that implicit integration schemes can be used without extra cost, which makes the approach capable of solving stiff systems.

Gradient matching
The previous parameter identification methods aim to compute a solution of the dynamical system and compare it to the measurements.In contrast to these integration-based methods, gradient matching [8] takes an approach that avoids solving the dynamical system.The slope of the measurements-the "measurement derivative"   x-is computed and compared against the model function , which results in the optimization problem The main challenge is to compute an adequate approximation of the measurement derivative.Various numerical methods exist for this purpose [15].The only variables in Problem (3) are the model parameters, which makes it very efficiently solvable.A disadvantage of gradient matching is that all states need to be measured to evaluate the function .Enhancements of the method solve this issue [10][11][12].

BENCHMARK SETUP
For the benchmark, we extend the formulation of Problem (1).First, we relax the assumption that all states have to be measured.We introduce the set of measured state indices  and take only those into account in the objective.Second, measurements of multiple experiments are incorporated.For each experiment, the dynamical system needs to be solved independently.We assume the measurements split across   ∈ ℕ experiments.A superscript [] indicates the dependence on an experiment.While Problem (1) imposes no restrictions on the model parameters, the extended problem formulation allows constraints of the form   ≤ () ≤   .The function  ∶ ℝ   → ℝ   is assumed to be twice continuously differentiable.This formulation allows general nonlinear constraints on the model parameters and simple box bounds.We allow the residuals to be weighted differently to emphasize particular measurement points.Additionally, we add normalization factors to scale the objective to a numerically better range.All these factors are combined as scaling coefficients with the weights The basis of the benchmark is a selection of 298 parameter identification problems that match Problem (4) from the collection of Schittkowski [14].It contains a variety of examples from both academia and real-world applications.We define two error measures to quantify the performance of a parameter identification method on a particular benchmark problem.If reference values  ref ∈ ℝ   for the model parameters are known, the averaged relative model parameter error can be used to quantify the accuracy of a method.Here,  * ∈ ℝ   denotes the identified model parameters.However, the error ( 5) is unsuitable for real-world problems for which the reference model parameters are naturally unknown.Instead, if available, we use the identified model parameters  * and initial states  * [] 0 ∈ ℝ   to solve the model with an adaptive Dormand-Prince method [16].We evaluate the solution at the measurement time points and compute the model simulation  * [] , ( * [] 0 ,  * ).With this, we can define the state-wise averaged weighted measurement relative error We denote (, ) as the error produced by the parameter identification method  on the benchmark problem , where  is a placeholder for   or  x.To compare two methods  1 and  1 , we use the performance ratio introduced by Kuhlmann [17]: where a method is considered to fail, if the error exceeds a threshold.If both methods fail, the performance ratio ( 7) is undefined.Plotting the performance ratios as a bar plot over the benchmark problems gives a clear visualization of the performance comparison of the two methods (see Figure 1).We also use it to assign a total score to the comparison: with the number of benchmark problems   .For simplicity, we exclude failed instances in Equation (8).A positive total score indicates that  1 performs better on all benchmark problems, and a negative total score implies that  2 is better.

NUMERICAL RESULTS
Gradient matching is only applicable to problems where all states are measured (see Section 3.4).This excludes many problems from the benchmark library.To obtain a better statistical basis, we, therefore, first compare only the three integration-based methods single shooting, multiple shooting, and full discretization.We solve all optimization problems with the nonlinear programming solver WORHP [18].As an initial guess of , we use the values provided by the benchmark library.All variables representing states are initialized by interpolation of the measurements.

Comparison of integration-based methods
For the experiments, we need to choose a configuration of the parameter identification methods.We use a fourth-order Runge-Kutta scheme without error control for the numerical integration in single and multiple shooting.For single Bar plots of the performance ratios for all test problems of single shooting (green), multiple shooting (red), and full discretization (blue).The colors illustrate which method performed better.The sections correspond to the different cases in the definition of the performance ratio (7).The numbers count the instances of problems in each section.A dark color indicates that the other method of the comparison failed.The number in the bottom left corner counts how many problems are solved successfully, and the number in the bottom right corner counts how many problems are not solved successfully.The total score ( 8) is given in the top right corner.
shooting, we choose 1000 integration steps.For multiple shooting, we introduce 20 shooting nodes and 50 integration steps per shooting interval.We apply a trapezoidal scheme as the discretization scheme in full discretization and choose 1000 discretization points.We solve all 298 test problems with these configurations and compute the performance ratios (7) for each combination of parameter identification methods as well as the total score (8).As the error measure, we use the state-wise averaged weighted measurement relative error (6) and a success threshold of 0.5.Figure 1 shows bar plots of the performance ratios for each test problem.It also indicates how many problems are solved successfully and shows the total score of the comparison.The comparison of single shooting and multiple shooting (see Figure 1A) has a total score  = −0.01,which does not indicate a notable difference between the methods.However, we observe that single shooting fails in more cases in which multiple shooting succeeds (33 problems) than vice versa (24 problems), which expresses that multiple shooting performs slightly more robust than single shooting.Nonetheless, the 24 problems for which multiple shooting fails while single shooting succeeds are surprisingly high.If both methods find a solution, single shooting is more precise in the majority of the cases (76 problems).
The comparison of full discretization and single shooting clearly shows an advantage for full discretization (see Figure 1B).This is not only reflected in the positive total score  = 0.17 but also in the fact that full discretization solves 61 problems for which single shooting fails.We conclude that full discretization is much more robust than single shooting.However, the experiments show that single shooting is more precise than full discretization if both methods succeed.In that case, single shooting finds a more precise solution for 87 problems, while full discretization is more precise for only 65 problems.
Figure 1C shows the performance comparison between full discretization and multiple shooting.The observations are similar to the previous comparison.On the one hand, full discretization shows more robust behavior as it solves many problems for which multiple shooting fails (49 vs. 8 problems).On the other hand, multiple shooting computes more precise results if both methods succeed (97 vs. 65 problems).
In addition to the previous observations, we notice that many problems are not solved successfully by any of the three methods.We think that this is because the same method configuration is used to solve all test problems.The results demonstrate that finding a configuration that works well for diverse problems is challenging.F I G U R E 2 Bar plots of the performance comparisons of single shooting (green), multiple shooting (red), full discretization (blue), and gradient matching (cyan) for a subset of the test problems.For a detailed explanation of the visualization, see Figure 1.

Comparison of integration-based methods and gradient matching
In the second comparison, we also include gradient matching.Therefore, we consider only the test problems, which provide measurements for all states.As gradient matching does not identify the initial states, we use the averaged relative model parameter error (5) as the error measure.This also excludes many problems from the benchmark library.Eventually, we consider 151 problems for the comparison.We configure the integration-based methods as in Section 5.1.We use finite differences to compute the measurement derivative in gradient matching.Analogously to Section 5.1, we compute and visualize the performance ratios.Figure 2 shows the results.We observe a clear advantage for the integration-based methods.For all three comparisons, the total score is positive, indicating that the integration-based methods perform better than gradient matching.We see that gradient matching is notably less precise than single shooting, multiple shooting, and full discretization.If both methods succeed, gradient matching finds the better solution for only 6, 7, and 11 problems, respectively.We see one reason for the bad performance of gradient matching in the simple approach of calculating the measurement derivative.Finite differences lead to poor approximations for noisy or low-resolution data.Eventually, this affects the accuracy of the parameter identification.
An interesting observation is that gradient matching relatively often finds a solution if the shooting methods fail, as visualized by the dark cyan bars (see Figure 2A,B).Even if this solution might be imprecise, it can serve as an initial guess for an integration-based method and be further improved.

CONCLUSION
In this paper, we presented a benchmark of the parameter identification methods single shooting, multiple shooting, full discretization, and gradient matching to get a better insight into their performance.Our experiments showed that full discretization is the most robust method.However, the shooting methods computed more precise results than full discretization.This suggests using full discretization to identify model parameters, which can be further improved with shooting methods.We also observed that gradient matching was less precise than the integration-based methods.However, it relatively often found a solution when the other methods failed.It could be used to calculate an initial guess for the integration-based methods.Our results must be viewed carefully, as for all test problems, we used the same method configurations, which may not always have been the optimal one.We experienced that many problems of the benchmark library were still not solved successfully.This fact emphasizes the challenging nature of parameter identification problems.Finding a method configuration that works well for diverse problems is difficult.It would be advantageous to involve more test problems in the benchmark with gradient matching to improve its statistical basis.Therefore, more sophisticated gradient matching approaches that overcome the limitations of the simple algorithm need to be explored.It would be interesting to study if advanced gradient matching methods can compete with integration-based parameter identification methods.Another topic of future research should be to increase the number of successfully solved problems in the benchmark.Experiments with different method configurations would be a starting point.Ultimately, it would be desirable to be able to select and configure the optimal method for a particular parameter identification problem.

A
C K N O W L E D G M E N T S Open access funding enabled and organized by Projekt DEAL.O R C I D Marek Wiesner https://orcid.org/0000-0002-6230-5622R E F E R E N C E S