Identification of mechanical properties of arteries with certification of global optimality

In this study, we consider identification of parameters in a non-linear continuum-mechanical model of arteries by fitting the models response to clinical data. The fitting of the model is formulated as a constrained non-linear, non-convex least-squares minimization problem. The model parameters are directly related to the underlying physiology of arteries, and correctly identified they can be of great clinical value. The non-convexity of the minimization problem implies that incorrect parameter values, corresponding to local minima or stationary points may be found, however. Therefore, we investigate the feasibility of using a branch-and-bound algorithm to identify the parameters to global optimality. The algorithm is tested on three clinical data sets, in each case using four increasingly larger regions around a candidate global solution in the parameter space. In all cases, the candidate global solution is found already in the initialization phase when solving the original non-convex minimization problem from multiple starting points, and the remaining time is spent on increasing the lower bound on the optimal value. Although the branch-and-bound algorithm is parallelized, the overall procedure is in general very time-consuming.


Introduction
The initiation and development of cardiovascular diseases has been associated with changes in the mechanical properties of the underlying arterial tissue [1,2]. Measures reflecting arterial stiffness, such as the pulse wave velocity (PWV) [3], the stiffness index (β) [4], and the pressure-strain elastic modulus (Ep) [5], are often used to assess patients in the clinic [2,6,7]. The amount of information obtainable from these popular clinical measures is limited, however, since they typically average stiffness over the artery or assume a constant stiffness despite the distinctive non-linear stiffening behavior of the arterial wall [8]. Several research groups have tried to address these shortcomings and proposed methods to identify more realistic patient specific mechanical properties of arteries in vivo (in the living human) [9-B Jan-Lucas Gade jan-lucas.gade@liu.se 1 Linköping University, Linköping, Sweden 13]. These studies use non-linear continuum-mechanical models and identify the material parameters by fitting the models to clinical data. The fitting is performed by solving an optimization problem; typically minimization of the non-linear least-squares differences of measured and predicted pressure-radius response.
Although these new methods offer more insight into the mechanical properties of arteries in vivo, they suffer a fundamental problem: the non-linear optimization problem to be solved is non-convex and generally possesses local solutions which are not the global solution. This becomes a problem when the new information is used to develop diagnostic tools and intervention criteria. Intrasubject comparison requires a set of parameters for each patient which is uniquely relatable to a 'normal' set for healthy individuals, and the global solution is a natural candidate.
The issue of non-convexity is commonly addressed by using a gradient-based optimization algorithm which finds a local solution in combination with multiple starting points [11][12][13]. The best local solution is then taken to be the global solution. Unfortunately, there is no guarantee that this heuristic method identifies the global solution. Furthermore, no estimate can be provided for the difference between the best local and the true global solution.
In biomedicine, simulated annealing, particle swarm and genetic algorithms are frequently used to solve optimization problems [14][15][16][17], but these algorithms are also of heuristic nature and do not necessarily identify the global solution. Deterministic global optimization methods such as branch-and-bound (B&B) on the other hand, guarantee that in finite time, a solution is found whose optimal value is at worst a user-specified value higher than that of the global solution [18].
B&B algorithms have been used previously to solve least-squares minimization problems to global optimality [19,20]. These studies show that the algorithm's efficiency strongly depends on the reformulation of the original optimization problem and benefits from a low number of unknown parameters and simple functional expressions. The parameter identification method in Gade et al. [13] exhibits these characteristics, thus making B&B a suitable candidate when searching for the global solution. General purpose global optimization solvers exist that can be used as black-box solvers [21][22][23][24], but efficient global optimization requires exploitation of problem-specific characteristics, something which is difficult without access to source code.
In this paper, we propose a B&B-type algorithm for the parameter identification method presented in Gade et al. [13] and explore the feasibility of computing global solutions for the mechanical properties of arteries in vivo. In particular, we compute global solutions for the clinical data of three subjects. For each data set we analyze four successively larger parameter regions around a candidate global solution to study the algorithms efficiency.

Parameter identification method
The parameter identification method is taken from Gade et al. [13]. For a thorough description and discussion, the reader is referred to the original paper.
In the parameter identification method, an artery is treated as a homogeneous, incompressible, thin-walled cylinder which consists of an isotropic matrix with embedded collagen fibers. In the unloaded state, i.e. outside the human body, the artery has an inner radius R i , a wall thickness H , and a length L, see Fig. 1a. In vivo, the artery is subjected to the blood pressure P from the inside and stretched to a length l which is invariant with respect to the changing blood pressure [25], see Fig. 1b. The stretch in the axial direction, i.e. λ z = l/L,

Unloaded state
Loaded state (a) (b) Fig. 1 Unloaded (a) and loaded (b) state of an artery. The coordinates r , θ , and z of the cylindrical coordinate system are associated with the radial, circumferential, and axial direction, respectively. The thin dark-grey lines represent the collagen fibers embedded in the isotropic matrix which is colored in light grey results in an axial reaction force F. Furthermore, the inner radius and wall thickness in the loaded state are denoted with lower case letters, i.e. r i and h, and the cross-sectional area of the arterial wall A = 2πr i h + π h 2 is constant because the arterial tissue is assumed to be incompressible and the axial stretch is constant. In this loaded state, two sets of stresses are calculated for an artery: equilibrium stresses and constitutively determined stresses. By stating equilibrium in the circumferential and axial direction, the corresponding equilibrium stresses are calculated as, i.e. Laplace laws, respectively. The axial reaction force cannot be measured in vivo and is estimated by assuming that the ratio between the axial and circumferential stresses is γ = 0.59 at the mean arterial blood pressureP = 13.3 kPa [9]. Using Eqs. (1) and (2), the axial force is estimated as where the inner radiusr i and the wall thicknessh are associated withP. The second set of stresses is determined using the Holzapfel-Gasser-Ogden (HGO) constitutive model [26]. Following Gade et al. [13], the constitutively determined stresses are: in the circumferential direction and in the axial direction where c > 0 is a material constant describing the isotropic matrix, k 1 , k 2 > 0 are material constants associated with the embedded collagen fibers, and β ∈ [0, π/2] is the pitch angle of the symmetrically arranged collagen fibers relative to the circumferential direction, see Fig. 1a. Furthermore, I 4 is the squared stretch along the collagen fibers calculated as and λ θ is the stretch in the circumferential direction defined as the ratio of the loaded to the unloaded mid-wall circumference, i.e.
where H has been replaced by R i , λ z , r i , and h using that the wall volume is constant due to incompressibility. The equilibrium stresses in Eqs. (1) and (2) are fully determined given a data set comprising time-resolved blood pressure (P) and inner radius (r i ) measurements sampled n times and information about the cross-sectional area (A) of the arterial wall. The constitutively determined stresses in Eqs. (4) and (5), however, contain six unknown model parameters: the unloaded inner radius R i , the axial stretch λ z , the material constants c, k 1 , k 2 , and the angle β. These parameters are determined by solving the following non-convex weighted least-squares minimization problem: where UP denotes upper problem in anticipation of the B&B algorithm, κ = R i , λ z , c, k 1 , k 2 , β is the parameter vector, j = 1, . . . , n indicates a data point sample, and the superscripts L and U denote lower and upper bound, respectively. The weighting factor ψ = 0.99 is used to let the objective function be dominated by the circumferential stresses [13]. The non-linear inequality constraints in (UP) are not present in the original optimization problem and are introduced to reduce the parameter space to obtain physiologically reasonable values, see the Discussion. The fitting ranges and the limits on the non-linear constraints are based on experimentally observed values [27][28][29][30] and summarized in Table 1.
The parameter identification method in Gade et al. [13] has been numerically validated and the 95% limits of agreement have been determined for each of the six parameters, see Table  2. These limits represent the interval around the identified value in which the true parameter is lying and are used to create the vicinity regions around a candidate global solution later on in the Sect. 5. For parameter k 2 the difference of the identified and the correct value increases as the correct value increases. To compensate for this systematic error, the 95% limits of agreement are based on the ratio instead [13].

Global optimization approach
A B&B algorithm is used to solve the problem (UP) to global optimality, meaning that the relative difference between the optimal values of the estimated and the true global solution is less than the tolerance ε. The basic idea is to generate a sequence of non-increasing upper bounds UB and a sequence of non-decreasing lower bounds L B on the global solution. By successive subdivision of the parameter space along the B&B tree, ε-convergence is Lower limit (L) 1 1 1 · 10 −3 1 · 10 −3 1 · 10 −3 0 0.5 1 0 a Upper limit (U) 20 1.5 1 · 10 3 1 · 10 3 1 · 10 3 90 2 2 20 a This constraint is implicitly enforced by the constraint on I 4, j  [18]. In each region of the parameter space, the original non-convex problem, i.e. the upper problem (UP), is solved to local optimality and the upper bound is updated in case this local solution provides a lower value on the upper bound. Then the convex relaxation of the original problem, which will be referred to as the lower problem (LBP), is solved to local, hence global optimality to establish a new lower bound on the global solution in the current region. If the lower bound is above the upper bound, the current region cannot contain the global solution and can therefore be excluded. Otherwise the current region is added to a list of active problems A which possibly contain the global optimum and the process is repeated until the relative difference between the upper and lower bound is less than the specified tolerance ε.

Construction of convex relaxation
To facilitate the construction of the convex relaxation, the model equations (4) and (5) are not introduced in the objective function but enforced as constraints instead [19]. This makes the objective function convex in the auxiliary variables σ mod θθ, j , σ mod zz, j which represent the model stresses and are introduced as additional unknowns to the optimization problem. The added constraints enforcing the model equations are non-linear due to the presence of bilinear, fractional, and componentwise convex terms. To facilitate the creation of convex relaxations of the model equations, the terms r i, j , R i , λ z , β are replaced by scaled counterparts according to where s = 1000 is a scaling factor chosen such that the magnitude ofR is similar to the other unknown parameters. After substitution of the scaled counterparts in the model equations, the following auxiliary variables are introduced: where → means an auxiliary variable replacing the term following the arrow. The auxiliary variables w 1 , w 2 , w 3 , w 5 , w 6 , w 7 , w 8 , w 9 , w 12, j , w 13, j replace bilinear terms, w 4 a fractional term, and w 14, j componentwise convex terms. For each of the bilinear, fractional, and componentwise convex auxiliary variables in Eq. (9), four inequality constraints are added to the optimization problem, see "Appendix A". The auxiliary variables w 10 and w 11, j are enforced as linear equality constraints and are introduced to decrease the total number of auxiliary variables. The resulting intermediate optimization problem is where u = κ, σ mod θθ , σ mod zz , w ,κ = R ,λ z , c, k 1 , k 2 ,β , and w contains the auxiliary variables w a . Due to the introduction of auxiliary variables, the degree of non-linearity in the model equations has been greatly reduced. However, the model equations in (INP) still have componentwise convex terms of the form x exp(y) which are addressed according to "Appendix A.4". For these terms, instead of introducing new auxiliary variables, the relaxation is incorporated directly into the model equations. The complete convex relaxation of (UP) now reads constraints for auxiliary variables in Eq. (9), The non-linear inequality constraints present in (UP) are enforced as simple box constraints on the corresponding (auxiliary) variablesR, w 11, j , and w 14, j . The bounds on auxiliary variables w 11, j , w 12, j , and w 13, j are non-trivial due to the non-linear dependence onR,λ z , andβ, see "Appendix B". Furthermore, tight bounds on the auxiliary variables representing the model stresses are added, see Sect. 4.1.
The lower bounding problem (LBP) has in total 16+6n variables and 37+21n constraints. A comparison of the size of (UP) and (LBP) is provided in Table 3.

Branching method
A crucial choice for the efficiency of a B&B algorithm is the branching method, i.e. the selection of the parameter to branch on. The proposed strategy follows the ideas in Esposito and Floudas [19] and considers the deviations between the values of the auxiliary variables w a ∈ w and the non-convex terms they replace expressed in terms of the six model parameters p ∈κ. The deviations are thus calculated as where w * a denotes the value of the ath auxiliary variable at the solution of (LBP) and w a κ * is the ath auxiliary variable expressed as an explicit function, see Eq. (9), of the model parameters and evaluated at the solution of (LBP). The deviation of auxiliary variable w 1 is for example defined as In order to account for w 14, j being an exponent, the associated deviation is calculated using instead. Given the deviations in Eqs. (11) and (13), a region-reduction measure is introduced for every model parameter according to where p U , p L are the upper and lower bounds of parameter p in the current region and p U,orig , p L,orig are the original upper and lower bounds. Finally, a total deviation δ p is calculated for every model parameter p by summing up all δ a 's in which this parameter is present and multiplying by the corresponding region reduction measure r p . The parameter which has the highest δ p is chosen to branch on and the corresponding region is bisected.

Proposed branch-and-bound algorithm
The rudimentary B&B algorithm described in the beginning of Sect. 3 is adapted to the properties of the problem at hand. In Gade et al. [13], the non-convex problem is solved from several starting points and the solution with the lowest error function value is taken to be the global solution. The same methodology is applied here when solving the upper problem for the first time. The supposedly good estimate of the global solution allows for not solving the upper problem until the lower bound reaches a user-specified threshold t. The solution of both (UP) and (LBP) can be accelerated by using a good starting point. Since the upper problem is solved prior to the lower problem, little information about the current region is available. Therefore, when solving (UP) the only criteria for choosing a starting point is whether it is feasible or not. For that purpose ten feasible or nearly feasible starting points are generated using Latin Hypercube Sampling [31] in the current parameter region. The starting point of the lower problem is determined by using the solution of the upper problem if available, otherwise it is generated using the same method as for the upper problem.
To speed up the B&B algorithm even further it is parallelized. Since it is not essential for the B&B algorithm that the sequence of lower bounds is always non-decreasing, it is possible to split the list of active problems A into A 1 , . . . , A v , where v is the number of available CPUs, and assess them individually. 1 The individual assessment is done for a certain time T and then A 1 , . . . , A v are merged again into one list of active problems and the process is repeated until in less than 1000 active regions the lower bound is below a user-specified threshold t. From then on, the B&B algorithm is executed serially to make sure that the sequence of lower bounds is non-decreasing considering all active regions. The whole B&B algorithm is summarized in Algorithm 1.

Implementation
The proposed B&B algorithm is implemented in Matlab R2019b [32]. The upper problem (UP) is solved using the function fmincon with the interior-point optimization algorithm. For numerical efficiency the model parameters c, k 1 , k 2 , and β are replaced in (UP) by scaled counterparts [13], c = ec, k 1 = ek 1 , k 2 = ek 2 , and β = arcsin β , where the superscribed tilde indicates a scaled parameter. The analytical gradient of the objective function, the analytical gradient of the non-linear inequality constraints, and the analytical Hessian of the Lagrangian are determined with Maple 2019.1 [33] and supplied to fmincon. To solve the upper problem from 10 starting points, the MultiStart-class feature is used. The lower problem (LBP) is solved using the Interior Point Optimizer (IPOPT ) software package [34] which is interfaced to Matlab as mex code. Again, the analytical gradient of the objective function, the analytical gradient of the constraints, and the analytical Hessian of the Lagrangian are provided.
Although both optimization algorithms are based on the interior-point method, fmincon solves (UP) and IPOPT (LBP) more efficiently. The same stopping criteria are used for determine u L , u U for region κ L , κ U and u SP based on κ * solve (LBP) in region u L , u U using u SP to get L B iter BB i else execute Algorithm 2 and determine lowest lower bound L B of A iter BB ← iter BB + 1 end end both algorithms: absolute change of objective function value must be less than 10 −10 , and constraint violation must be less than 10 −8 .

Material
The material for this study comes from the abdominal aorta of three healthy, non-smoking Caucasian males, representative for each age group. Subjects I, II, and III are 25, 41, and 69 years of age, respectively. The pressure-radius data is taken from Sonesson et al. [35] and for details about the data acquisition the reader is referred to that paper. From the raw data we create a pressure-radius curve consisting of n = 18 equidistant data points for each subject according to Stålhand [36], see Fig. 2. The total number of data points is chosen as small as possible, see the Discussion. Neither the deformed wall thickness nor the deformed wall cross-sectional area were recorded in Sonesson et al. [35]. The deformed wall cross-sectional area A is therefore estimated according to A = 19.6+0.8 · age, where A is in mm 2 and age is in years [36].

Bounds on the auxiliary variables representing model stresses
The pressure-radius measurements show distinctive hysteresis, see Fig. 2. When the input data is converted into the Laplace stresses according to Eqs. (1) and (2), the hysteresis is still present, see Fig. 3. The constitutive model in Sect. 2 assumes a purely elastic material, however, and neglects hysteresis accordingly. The curves describing the model stresses in the  circumferential and axial direction must, therefore, be (mostly) contained within the corresponding hysteresis loops if they are to be the best solution to the least-squares optimization problem. Individual lower and upper bounds on each auxiliary variable representing the model stresses are therefore introduced based on the hysteresis loops of the Laplace stresses. Afterwards the average difference between the upper and lower bounds, excluding the endpoints, is calculated and the stress ranges smaller than this average difference are symmetrically enlarged. The lower and upper limits for the endpoints are created by using twice the average stress difference. The resulting bounds on the auxiliary variables representing the circumferential and axial model stresses are shown in Fig. 3.

Results
The proposed B&B algorithm is applied to the clinical data and four successively increasing regions in the parameter space around a candidate global solution for each of the three subjects is investigated. The candidate global solutions are determined by solving (UP) with 100 starting points generated with Latin Hypercube Sampling over the entire parameter space defined by the fitting ranges in Table 1. For subjects I, II, and III, respectively, 81%, 79%, and 64% of all starting points end up at the same solution with the lowest objective function value. The fact that the same solution is not always found demonstrates the non-convexity of the optimization problem and the associated existence of local solutions which are not the global solution. The candidate global solutions are summarized in Table 4. Around these candidate solutions four vicinity regions are established: 1%, 10%, 100%, and 500% of the limits of agreement for each parameter, see Table 2. The vicinity regions respect the parameter boundaries in Table 1 and the explicit parameter spaces are summarized in Table 5. For each subject and each vicinity region, a stopping-tolerance in the range ε = 0.01−0.095 is required for the B&B algorithm, with larger tolerances for larger vicinity regions, see Table 6. The threshold t until the B&B algorithm is executed in parallel is put to 90% of the objective function value of the respective candidate global solution.
During the initialization step of the B&B algorithm, the respective candidate global solution is identified in every case and a better solution for (UP) is never found. The qualitative progress of the difference between the upper and lower bound is similar for each subject in the same vicinity region, see Fig. 4. For many iterations, the relative difference is unity, i.e. L B = 0, followed by a comparably rapid drop until the lower bound has reached approximately 90% of the upper bound and the progress is substantially slowed down. The total amount of B&B iterations to achieve the required tolerance is large, especially for the two largest vicinity regions. Although each B&B iteration requires only a fraction of a second, the total CPU solution times are huge and are only manageable through parallelization, see Table 6.

Discussion
In all cases tested, the candidate global solution is identified during the initialization step of the B&B algorithm and no better solution is found, even after solving (UP) several million times in the larger parameter regions. This is not an unusual behavior of deterministic global optimization methods which often require many iterations to verify that a promising solution found early in the process is indeed the global solution [18]. That no better solution of (UP) is found highlights the successful combination of using multiple starting points, selection of starting points and the rescaling when solving (UP) in the initialization step.
Initially a tolerance of ε = 0.01 was required for each subject and vicinity region. It was, however, soon realized that ε = 0.01 is not feasible for every optimization run, since the required computational power exceeded our available resources. An in-house cluster with 108 2.6 GHz cores and 576 GB RAM, as well as two personal computers, were running for approximately four months to solve all optimization runs with a combined CPU time of 5.4 years. 2 Especially the optimization runs for subject III lasted for a long time because the parameter space within each vicinity region is larger compared to the other two subjects. The increased parameter spaces for subject III are due to the higher k 2 -value of the candidate global solution and the relative confidence interval for that parameter, see Sect. 2. The required tolerance is therefore lowered on an individual basis such that (UP) is at least solved a few million times during which no better solution must be found.
In its current form, the proposed B&B algorithm does not utilize parallelization once the lower bound has been lifted to 90% of the upper bound. This threshold is chosen for two reasons. First, the upper problem is solved again when the threshold has been passed. It is assumed that (UP) becomes practically convex in the small parameter space of the   The reported times correspond to CPU time, i.e. the time a single core would have required, and size(A) denotes the number of active regions at the end of the respective optimization run respective active region and that a better solution is found in this state. The best solution of (UP), i.e. the lowest upper bound, should be available throughout the B&B algorithm to exclude regions as early as possible and communication between the active regions split into parts is not provided. Second, when the list of active regions is split into several parts, the sequence of lower bounds will be different in each individually processed part. The more parts the list of active regions is split into and the longer these parts are processed individually, the more the lower bounds will diverge and the iteration number to achieve a certain tolerance will be overestimated. This overestimation has a step-like appearance where the tolerance progress is shifted towards higher B&B iterations. This is especially pronounced for subject III whose optimization runs have been parallelized the most, see Fig.  4. In order to get the correct tolerance progress towards the end of each optimization run, the B&B iterations are performed serially once the threshold is passed. Both aspects are of minor concern and the algorithm should be parallelized throughout if the goal is solely to identify mechanical properties of arteries with certification of global optimality. Non-linear inequality constraints are added into the original optimization problem (UP). Their main purpose is to reduce the parameter space to the physiological range. Circumferential stretches of 0.8 < λ θ < 1.6 are reported for the abdominal aorta [27,37] and we take λ L θ = 0.5 and λ U θ = 2. The invariant I 4 is a macroscopic quantity representing the squared stretch of collagen fibers. Hence, the upper bound on I 4 can be related to the squared failurestretch of collagen fibers and is estimated to be I U 4 = 2 [38,39]. The lower limit on I 4 is introduced to facilitate the construction of the lower bounding problem. Collagen fibers are, generally, considered to support tensile loads but buckle in compression [26]. Hence, the contribution of collagen in the constitutive model, cf. Sect. 2, should be omitted if I 4 < 1. Due the difficulty in constructing a valid underestimator which is smooth at the transition from non-recruited (I 4 < 1) to recruited collagen fibers (I 4 ≥ 1), collagen fibers are assumed to be always in extension. In order to test the validity of this assumption, the model parameters of the three subjects are identified with the original parameter identification method in Gade et al. [13] which does not restrict I 4 . Collagen fibers are recruited throughout the cardiac cycle in all cases, thus I 4 ≥ 1 appears to be a valid assumption. Lastly, the upper bound of the exponent in the exponential term appears explicitly in some constraints of (LBP), i.e. exp w U 14 = exp k U 2 I U 4 − 1 2 . Due to the high upper bound on k 2 in Table 1, the exponential expression might become very large. To prevent numerical problems, the exponent is limited to 20. The major determinant for the efficiency of a B&B-type algorithm is the creation of the convex underestimator. The convex underestimator must be as close as possible to the original function while the number of additional auxiliary variables and constraints should be kept to a minimum. Several alternative formulations of (LBP) were examined, but the presented one turned out to be most efficient. In particular, the newly proposed underestimator for componentwise convex terms [40] gave a substantial speed-up compared to the more general α-estimator [41]. Furthermore, the introduction of physiologically inspired auxiliary variables and the tight bounds on the auxiliary variables representing the model stresses based on the hysteresis loop, cf. Sect. 4.1, decrease the required computational time. With respect to the branching method, other variants were tested, in particular strategies 1-4 proposed in Adjiman et al. [42]: (S1) bisect the parameter which has the largest region reduction measure defined in Eq. (14); (S2) determine the auxiliary variable whose underestimator possesses the highest separation distance to the actual term it replaces. Bisect the involved parameter according to S1; (S3) Determine the auxiliary variable which differs the most from the actual term it replaces at the solution of (LBP), cf. Eq. (11). Bisect the involved parameter according to S1; and (S4) calculate the deviation of each auxiliary variable according to S3. For every parameter, sum up all the deviations that the parameter is involved in and bisect the parameter with the highest sum. The branching method described in Sect. 3.2 was found to be superior, however. Branching on auxiliary variables is not used since updating all lower and upper bounds is non-trivial for branching on a member of u compared to a member ofκ.
From an implementation point of view, one point is worth discussing. By looking at the total amount of variables and constraints of (LBP), see Table 3, it is obvious that the amount of samples, i.e. the number of pressure-radius pairs, must be as low as possible. For the three investigated subjects in this study, n = 18 is found to be the lower limit to get a solution which is independent from the number of samples. For other subjects this number might be different.
Considering the clinical application, the proposed B&B algorithm appears to be impractical. The solution time is orders of magnitude too large, even if only the vicinity around a candidate global solution is explored and the process is massively parallelized. Furthermore, the smaller the relative tolerance is required to be, the longer the identification process will take since the B&B algorithm slows down substantially towards smaller relative differences between lower and upper bound. However, no better solution than the one identified during the initialization is found, even after solving (UP) several million times. This indicates that the heuristic approach of using multiple starting points works well for this kind of optimization problem. The original parameter identification method in Gade et al. [13] with its comparably small computational requirements, therefore, provides a good estimate of the global solution if it is not even obtained.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Code availability All necessary information for code development is included in this published article.
Availability of data and material All data analyzed during this study are included in this published article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A.1 Bilinear term
A bilinear term of the form x y, (x, y) ∈ x L , x U × y L , y U ⊂ R 2 , is replaced by the auxiliary variable w xy and the following four linear inequality constraints are added to replace the non-linear equality constraint w xy = x y: The first two constraints represent the convex hull [43] and the latter two the concave hull [44].

A.2 Fractional term
A fractional term of the form x y −1 , (x, y) ∈ x L , x U × y L , y U ⊂ R 2 >0 , is replaced by the auxiliary variable w xy and the following four inequality constraints are added to replace the non-linear equality constraint w xy = x y −1 : The first two convex constraints represent the convex hull [45] and the latter two linear constraints the concave hull [46].

A.3 Componentwise convex term xy 2
A componentwise convex term of the form x y 2 , (x, y) ∈ x L , x U × y L , y U ⊂ R 2 ≥0 , is replaced by the auxiliary variable w xy and the following four inequality constraints are added to replace the non-linear equality constraint w xy = x y 2 : x − x L y L 2 + y 2 x L − w xy ≤ 0, x − x U y U 2 + y 2 x U − w xy ≤ 0, −x y L 2 − y − y L x U y U + y L + w xy ≤ 0, −x y U 2 − y − y U x L y U + y L + w xy ≤ 0. (18) The first two convex constraints represent the convex hull [40] and the latter two linear constraints the concave hull, whose creation is inspired by the ideas for a bilinear term [44].

A.4 Componentwise convex term x exp y
A componentwise convex term of the form x exp y, (x, y) ∈ x L , x U × y L , y U ⊂ R 2 ≥0 , is replaced by the auxiliary variable w xy and the following four inequality constraints are added to replace the non-linear equality constraint w xy = x exp y: x − x L exp y L + x L exp y − w xy ≤ 0, The first two convex constraints represent the convex hull [40] and the latter two linear constraints the concave hull, whose creation is inspired by the ideas for a bilinear term [44].

B Bounds on auxiliary variables
In order to calculate the tightest possible bounds on the auxiliary variables, they are determined by writing them as explicit functions of the original model parameters κ. For most auxiliary variables the bounds are trivial and in the following only the non-trivial bounds are shown. The bounds on auxiliary variable w 11, j =r jR 1 −β +λ zβ − 1 are: w U 11, j = ⎧ ⎨ ⎩ w 11, j R U ,λ U z ,β U , ifλ U z ≥r jR U , In addition, I L 4, j − 1 ≤ w L 11, j , w U 11, j ≤ I U 4, j − 1 must hold.