Semi-probabilistic methods are widely used for structural engineering design and assessment. In such methods, the safety and performance adequacy of structures is typically tied to partial safety factors for load and resistance, and load combination factors. There exist some heuristic strategies for estimating these partial factors, such as the design value method and first order reliability coefficient method. However, when adopted for estimation of load combination factors, these strategies are either inaccurate or have non-unique estimates. Excessive conservatism in factor estimates is not desirable, particularly for the performance assessment of existing structures. In this study we propose a method for estimating load combination factors using a matrix linear algebra approach. Specifically, we develop a closed-form analytical expression to estimate unique load combination factors for typical code calibration problems. A comparative study of the existing heuristic approaches with the new approach is presented and demonstrated on two case studies. It is shown that the proposed method offers a valuable means of deriving load combination factors: one that avoids hueristics, and conservatism. 1. Introduction The basis for modern structural design codes are semi-probabilistic approaches such as load and resistance factor design (LRFD) methods. In these methods, the uncertainties within load and resistance are considered using partial safety factors for load and resistance. The structural loads can be time-variant or time-invariant, and multiple loads can act simultaneously. Structural engineering design involves consideration of these simultaneous loads to establish the most adverse scenario for design. However, the simultaneous occurrence of extremes of these loads is a rare event that must be accounted using a load combination factor in addition to the partial safety factors. Compared to partial safety factor estimation, the estimation of load combination factors has received little attention. The combination of time-variant loads is a complex problem of combining two or more stochastic processes Melchers and Beck (2017). In most code calibration problems, the simplifications introduced by Ferry Borges and Castanheta (1972) and Turkstra and Madsen (1980) are adopted. These simplifications describe strategies for combining load processes. However, they do not specify the estimation strategy for load combination factors for a semi-probabilisticdesignandassessment. Likewise, manystructuralreliabilitytextbooksdescribethedetermination ofpartialfactors,butnotloadcombinationfactors; e.g. NowakandCollins(2013);MelchersandBeck(2017);Haldar and Mahadevan (1999). One of the earliest descriptions of both partial and load combination factor calibration is provided in Thoft-Cristensen and Baker (1982). It adopts a heuristic design value method for a simplified estimation ∗Corresponding author Email addresses: colin.caprani@monash.edu (C.C. Caprani), shihab.khan@monash.edu (Mohammad S. Khan) Preprint submitted to Elsevier May 21, 2023 of partial and combination factors. Notably, Sørensen (2004) contrasts the design value method with a coefficient based estimation of the load combination factor using the design points obtained using first order reliability method (FORM). A founding technical document for structural reliability, ISO2394 (ISO (2015)), takes basis in Sørensen (2004) and Faber and Sørensen (2003) for the estimation of partial and combination factors. However, these strategies for the estimation of load combination factors have some deficiencies. Firstly, and as will be shown, the coefficient approach explained in Sørensen (2004) yields non-unique estimates of load combination factors; one set per load case. The solution typically adopted is then consider the maximum of the estimates for a given load, across all load cases. The second deficiency is then that such consideration adds unnecessary conservatism, which is often undesirable. Finally, these approaches are heuristic in nature and lack an analytical basis for the estimation of unique combination factors. To address these problems in load combination factor calibration, this study proposes a new approach for estimating load combination factors. We use matrix inversion to obtain a unique least-squares estimate of the load combination factors satisfying several loadcase limit states simultaneously. For the standard code calibration problem, an extremely useful closed-form expression for estimating load combination factors is derived. Considering two examples, the results from the existing heuristic approaches and the proposed approach are compared. It is shown that the proposed approach provides an optimum means for estimating load combination factors, particularly for problems involving three or more time-varying loads. The paper is structured as follows. Firstly, Section 2 presents the mathematical background for the code calibration problem and load combination factor estimation using FORM as a basis. Secondly, Section 3 presents the development for estimating unique load combination factors. Then Section 4 demonstrates the new approach on two example problems with varying time-variant loads and compares the results with those obtained using the existing heuristic methods. Finally, Section 5 discusses the results and concludes the study. 2. Background 2.1. FORM Basis This section briefly introduces the aspects of FORM relevant to the calibration of partial and load combination factors. Readers are referred to structural reliability textbooks for more details (e.g. Nowak and Collins (2013) and Melchers and Beck (2017)). At the most probable failure point in physical coordinates, the limit state function (LSF) calibrated to a chosen target reliability index, βT, is g(X∗ 1, X∗ 2, . . .) = 0such that Pr[g≤0] =pf= Φ(−βT) (1) where X∗ iis the design point of the ith random variable, Xi. Here, Equation (1) is generic and does not make any assumption regarding the structure of the LSF, which maybe linear or non-linear with respect to the ran- dom variables. In general the distributions of the random variables are non-normal. Therefore, isoprobabilistic transformations are used (e.g. Nataf or Rosenblatt), such as for uncorrelated variables: X∗ i=F−1 Xi[Φ(u∗ i)] 2 in which u∗ i=αiβis the corresponding design value in standard normal (or u-) space, and αiis its direction cosine. For design purposes, a nominal or characteristic value of the variable Xis used, Xk. Since the design values are intended to satisfy the LSF (Equation (1)), they will yield a design with reliability, β≥βT. Thus, the partial factor to apply to the characteristic value of the variable, to yield the design value corresponding to βT, is γi=X∗ i Xi,k. (2) 2.2. Code Format The prototypical equation to be satisfied in structural design is ϕRk≥X iγSiSki, (3) in which Rkis the characteristic value of the resistance variable, and Skithat of the ith loading variable. Each variable has its respective partial factor, γ, often denoted ϕin the case of the resistance variable. Equation (3) is linear with respect to the load variables. The linearity of the design equation is independent of the LSF equation, but it is not restricted to LSFs that are linear functions of the random variables. Indeed, a wide class of LSFs can be considered for design in the format of Equation (3) and so this study is limited to such design equations. The loading variables, Si, can be classified as permanent or varying. That is, permanent loads, such as dead loads—although uncertain—have just one realization from its distribution for the lifetime of the structure. Con- versely, varying loads, such as traffic or wind loads, have realizations occurring continuously. Of course, this can include periods where the value of loading is zero, such as when there is no traffic on the bridge (i.e. zero-inflated distributions). For a designer to satisfy Equation (3), the combination of permanent and varying loads should allow for each varying load achieving a maximum value while the others remain at their point-in-time values, within the time frame being considered. As a result, Equation (3) should be checked ntimes, as each varying loading variable achieves a maximum. As a special (but important) case, consider there to be a permanent dead load, G, a varying live load Q, and a varying wind load W. Equation (3) then becomes: ϕR≥γGGk+γQQk+γWWk. (4) The difficulty here is that the live and wind loads are varying, and (most likely) uncorrelated. So to satisfy Equation (4) when either varying load reaches a maximum within the period being considered (e.g. 1-year), the designer will then seek to satisfy two load cases: LC1: ϕR≥γGGk+γQmaxQmax k+γWpitWpit k(5) LC2: ϕR≥γGGk+γQpitQpit k+γWmaxWmax k (6) where the live and wind loads achieve their maximum ( max) distribution respectively, while the other remains at its point-in-time (pit) value. In Equation (5), each varying load is now characterized by four values rather than two: its characteristic point- in-time and maximum values, and their corresponding partial factors. This significant increase in the number of 3 parameters (i.e. 4n) required to represent each stochastic variable is clearly undesirable for practical structural design codes. To rationalize the number of parameters required, most code formats have adopted a model which can be expressed as ϕR≥γGGk+γSdSd,k+ψγSaSa,k, (7) in which dis the dominating varying load, ais the accompanying varying load, and ψis a load combination factor. This representation means that each load is parameterized by just a single partial factor γand a characteristic value Sk. It is the load combination factor ψthat then accounts for the reduced probability of both loads achieving their maximum value simultaneously. With this rationalization, Equation (5) can be expressed as: LC1: ϕR≥γGGk+γQQk+ψWγWWk, (8) LC2: ϕR≥γGGk+ψQγQQk+γWWk. (9) Some code formats will merge the load combination and partial factor for the accompanying load and specify a reduced (composite) partial factor (i.e. ψγ) when used in combination. Other code formats will retain the partial factor but provide for a reduced characteristic value (i.e. ψSk) when used in combination. Whatever means is used to represent the problem, the implication for code calibration is that the load combination factor must be estimated. 2.3. Current Methods 2.3.1. Design Value Method Thoft-Cristensen and Baker (1982) gives an empirical estimate of load combination factors. When the variables are ranked in order of their directional cosines αi, the ith variable can be considered to have a relative direction cosine of: ˜α=√ i−√ i−1. (10) This recommendation can be coupled with the αR=−0.8andαS= 0.7estimate (from ISO2394 (ISO (2015)) to yield estimates of partial and combination factors (Equations (2) and (11)). Although approximate, Thoft- Cristensen and Baker (1982) suggest that this approach gives reasonable results. 2.3.2. Coefficient Method Sørensen (2004) explains the coefficient method as the obvious extension of the FORM result for a single limit state function (or loadcase), to that for two loadcases. Consider the two limit state functions of Equation (8). For calibration, the design values for both loadcases are found for a single target reliability index, βT. This yields a set of partial factors (following Equation (2)) for each variable, for each limit state. That is, for example, the wind load will have γW,1andγW,2, corresponding to the calibration for each loadcase. Rationalization of these partial factors is then done such that a single value of the partial factor for each variable results, e.g. γW. With the single set of partial factors and the characteristic value known, the load combination factor for the secondary varying loading at the design point for its point-in-time loadcase becomes: ψS=S∗ pit γSSk. (11) 4 Extension of this approach to problems with n >2time-varying loads is then based on the identity: ψi,jγiSki=S∗ i,j for loadcase iand design variable j. When i=jvariable iis at its maximum, and so we define ψi,i= 1.0from which: γi=S∗ i,i Ski(12) as usual. But when i̸=jwe then have: ψi,j=S∗ i,j γiSki. (13) The difficulty here of course is that for n >2there are multiple values of ψi,j, not a single unique value of ψjas needed. That is, the coefficient method results in unique estimates for partial load factors γbut estimates of ψand ϕthat vary per load case. In such cases, there are too many degrees of freedom in the problem. Thus, a trial and error estimation is required, or more simply and conservatively, the maximum is chosen, ψj= max iψi,j. 3. General Method for Estimating Load Combination Factors 3.1. Matrix Solution Here we propose an revision to the coefficient method in which a unique set of load combination factors is found; a single one for each varying load. To do this, we use an orthogonal projection from the higher dimensionality of the coefficient method to the preferred dimension of the problem; namely, a least-squares projection. Consistentwithpractice, itissupposedthatasingleloadcombinationfactorisdesiredforeachvariable, andthat a set of partial factors has been determined. We consider that each time-varying load must be considered acting at its maximum (i.e. dominating), and accompanying varying loads are acting at their arbitrary point-in-time values. That is, we adopt a form of Turkstra’s rule for load combinations (Turkstra and Madsen, 1980). Thus, for nsuch loads, there will be nlimit state functions to be satisfied. To introduce the method, we consider an example with three time-varying loads, and then show how it can extend to the general case: ϕR≥γGG+γ1S1+ψ2γ2S2+ψ3γ3S3 (14a) ϕR≥γGG+ψ1γ1S1+γ2S2+ψ3γ3S3 (14b) ϕR≥γGG+ψ1γ1S1+ψ2γ2S2+γ3S3 (14c) in which the subscript kis dropped for clarity, it being understood that these variables are at their characteristic values for these design equations. Moving any non-varying load terms to the resistance side of the equation, and swapping, we write: 0ψ1+γ2S2ψ2+γ3S3ψ3≤ϕR−γGG−γ1S1 γ1S1ψ1+ 0ψ2+γ3S3ψ3≤ϕR−γGG−γ2S2 γ1S1ψ1+γ2S2ψ2+ 0ψ3≤ϕR−γGG−γ3S3 5 Now rewrite this in matrix form, separating the unknown ψifrom the remaining known values as follows:  0 γ2S2γ3S3 γ1S1 0 γ3S3 γ1S1γ2S2 0   ψ1 ψ2 ψ3  ≤  ϕR−Q1 ϕR−Q2 ϕR−Q3  (16) in which the total factored permanent, non-varying, and dominating characteristic design loads are denoted: Qi=γGG+γiSi. (17) For the limiting equality case, clearly Equation (16) can then be easily solved for the ψvector by matrix inversion to yield a single set of load combination factors, ψi, that satisfy Equation (14). The Appendix describes the detailed consideration of this method for generic limit-state functional forms, 3.2. Mathematical Basis In this section we develop some general mathematical results that will be used later. Consider the n-dimensional hollow (Gentle, 2017) square matrix Apartly formed as the outer product of a vector of ones 1and a row vector a= [a1, a2, . . . , a n]: A=1aT−diag(a) = 0a2. . . a n a10. . . a n ............ a1a2. . . 0 (18) Where the inverse of Aexists, it will be of interest to solve the following linear equation: Ax=b, (19) such that x=A−1b. Writing ˜a= [1/a1,1/a2, . . . , 1/an]as the vector of component-wise reciprocols of a, we claim for n≥2: A−1= (n−1)−1˜a1T−diag( ˜a) = −(n−2) (n−1)a11 (n−1)a1. . .1 (n−1)a1 1 (n−1)a2−(n−2) (n−1)a2. . .1 (n−1)a2 ............ 1 (n−1)an1 (n−1)an. . .−(n−2) (n−1)an (20) so that the solution of the linear equation becomes: xj=b1 (n−1)aj+b2 (n−1)aj+··· −(n−2)bj (n−1)aj+···+bn (n−1)aj giving: xj=P nbi−(n−1)bj (n−1)aj(21) where 1≤j≤n. This is the main result. 6 To verify Equation (20), we demonstrate that AA−1=I. To show this, we will show that the diagonal entry is1and the off-diagonal entries are 0. Consider then, row iofAmultiplied with column iofA−1: Ai,i= [a1, . . . , a i−1,0, ai+1, . . . , a n]×  1 (n−1)a1 ... 1 (n−1)ai−1 −(n−2) (n−1)ai 1 (n−1)ai+1 ... 1 (n−1)an   =i−1 n−1+ 0 +n−(i+ 1) + 1 n−1 =n−1 n−1 giving Ai,i= 1as required. Next, consider row iofAtimes column j̸=iofA−1for which the result should be 0. We take the case i < jbut the argument is the same considering i > j: Ai,j= [a1, . . . , a i−1,0, ai+1, . . . , a j−1, aj, aj+1, . . . , a n]×  1 (n−1)a1 ... 1 (n−1)ai−1 1 (n−1)ai 1 (n−1)ai+1 ... 1 (n−1)aj−1 −(n−2) (n−1)aj 1 (n−1)aj+1 ... 1 (n−1)an   =i−1 n−1+ 0 +j−1−(i+ 1) + 1 n−1+−(n−2) n−1+n−(j+ 1) + 1 n−1 =i−1 +j−1−i−1 + 1−n+ 2 + n−j−1 + 1 n−1 = 0 Thus Equation (21) is shown to be the solution to the linear equation Equation (19) when the coefficient matrix is Aas defined in Equation (18). 3.3. Closed-Form Solution Referring to the solution Equation (21), for the reliability problem of Equation (16), we see that xj≡ψj, the vector bi≡ϕR−γgG−γSiSi, and the matrix entries are aj≡γSjSj. Hence, when ϕRis the same (a typical 7 structural design) for each row (loadcase), we have X bi=n(ϕR−γgG)−X γSiSi (22) giving: ψj=n(ϕR−γgG)−PγSiSi−(n−1)(ϕR−γgG−γSjSj) (n−1)γSjSj and simplifying this we obtain the main result of this work: ψj= 1−PγSiSi−(ϕR−γgG) (n−1)γSjSj. (23) whichisaclosed-formexpressionforthe nloadcombinationfactorsgivencalibratedpartialfactorsandcharacteristic design values. 3.4. Discussion 3.4.1. Heuristics Equation (23) exhibits some noteworthy aspects. The numerator term in parentheses ( ϕR−γgG) represents the available capacity to carry the varying loads. Moreover, it is clear that the fraction must be ≤1for it to be valid, and so the numerator must be smaller than the denominator, in spite of the (n−1)multiplying term. Indeed, the numerator represents the excess loading demand, should all loads achieve their maximum values simultaneously. Later we denote this excess loading S+. This excess loading is then shared across the (n−1)loads which are not at their maximum, in proportion to their contribution to the excess loading. This fractional component of the excess loading is then subtracted from 1.0 to yield the component’s load combination factor. Together, these combination factors then serve to eliminate the excess loading demand, returning the problem to one where the loads are not simultaneously at their maxima. 3.4.2. Resistance Design Points In a reliability calibration (see later examples), each loadcase in Equation (16) will be separately calibrated. Consequently, there will (in general) be a different design point for the resistance variable, R∗ i=ϕiRin each row. Of course, for structural design purposes, the maximum calibrated resistance must be used in order to satisfy all loadcases. Therefore it is reasonable that in Equation (22) we consider that ϕRis the same for each equation, leading to the closed-form expression of Equation (23). Nevertheless, for some code calibrations it could be of interest to consider the resistance design points ϕRwithin each load case. In such cases, eq. (21) can be referred to directly, or inversion of the matrix in Equation (16) can be done with different design points in each row. 4. Examples Throughout the following examples, results are presented to five significant digits. This is to allow readers to confirm their own calculations with ours. We do not suggest that partial or combination factors should be any more than perhaps three significant digits in a design code intended for practice, nor do we suggest that the results are meaningful to five orders of magnitude. Indeed, in a real code calibration the results in these examples would be rounded off in a safe direction towards more memorable quantities, typically to just two decimal places, such as 8 γ= 1.15. However, we do not do this rationalization here so as to assess the performance of the various approaches, without the influence of this important heuristic stage of code calibration. That is, the following are numerical examples intended to highlight the various numerical approaches. The examples are implemented using The PySTRA Developers (2023), and an open-source jupyter-lab notebook containing the calculations is provided at Caprani and Khan (2023). 4.1. Two Varying Loads 4.1.1. Description Here we consider an example from Sørensen (2004, p. 190) with two time-varying load sources. The limit state function is g=zR−(0.4G+ 0.6Q+ 0.3W), (24) in which zis a scalar design parameter (e.g. a deterministic nominal strength design to achieve the target reliability) and the other variables have their usual interpretations. The basic variables are given in Table 1. The distributions specified for the varying loads are the annual maximum distributions. From the number of occurrences in the 1-year reference period, it can be found that the point-in-time distributions for the imposed and wind loads have expected values of 0.89 and 0.77 respectively (the CoV remains unchanged). Table 1: Random variables of Example 1; τis the occurrence rate of the varying load, and ris the number of such of such occurrences within the reference 1-year period. (*Note: although 360 recurrences per year is reasonable and sated in Sørensen (2004, p. 190), the problem was actually worked out using 2, and this ‘error’ is retained here for comparison with the source.) DistributionAnnual Max. Point-In-Time Characteristic Design Points E[X] CoV τ r Percentile xk x∗ 1 x∗ 2 z- - - - - - - 3.0431 3.0477 RLognormal 1.0 0.15 - - 0.05 0.7738 0.6553 0.6550 GNormal 1.0 0.10 - - 0.50 1.0000 1.0371 1.0371 QGumbel 1.0 0.20 0.5 years 2 0.98 1.5185 1.6235 1.5129 WGumbel 1.0 0.40 1 day 2 (360*) 0.98 2.0369 2.0171 2.2458 Equation(24)isusedtocalibratethedesignparameter zfortwoloadcases: whentheimposedloadhasitsannual maximum distribution, and the wind load is at its point-in-time distribution, and vice-versa. A target reliability index of βT= 4.3is used for the calibration. The resulting design points and design variable values for the two calibration loadcases are shown in Table 1. Further, it is considered that the characteristic values of the imposed and wind loads are given by the 98-percentiles of the annual maximum distributions, and so are Qk= 1.5185and Wk= 2.0369. With these values, the partial factors are: γQ=Q∗ 1 Qk=1.6236 1.5185= 1.0692, γW=W∗ 2 Wk=2.2458 2.0369= 1.1026. 9 4.1.2. Design Value Method We begin by considering the design value method described in Section 2.3.1 which gives (see Equation (10)): α= 0.7×(√ 2−√ 1) = 0 .29 and this compares to the actual FORM estimated values of α= 0.54(which are the same for both the imposed and wind loads). Nevertheless, when dominating, the design points in the physical space corresponding to βT= 4.3 for the imposed and wind loads are 1.945and2.891; and when secondary 1.251and1.502, respectively. From these design point estimates, the partial factors are 1.9454/1.5185 = 1 .2812for imposed load and 2.8908/2.0365 = 1.4192for wind load. Finally, the load combination factors become 1.2509/1.9454 = 0 .6430for imposed and 1.5019/2.8908 = 0 .5195for wind. The efficacy of the estimated factors can be evaluated by initiating a structural design with these estimates. For the example problem in consideration, such a design involves using the factors to estimate the corresponding scalar design parameter, z, and then conducting a reliability analysis with the estimated zand specified random variable statistics. Proceeding to a structural design with these values gives a design parameter of z= 2.6549yields: β= [3.7082,3.7011] which are significantly less than the intended βT= 4.3. Consequently, for the other numerical examples in this study, we will not consider the heuristic design value method further. 4.1.3. Coefficient Method Consider next the common coefficient method for estimating partial factors, described in Section 2.3.2. The load combination factors can be worked out using this method as: ψQ=Q∗ 2 γQQk=1.5185 1.6236= 0.9318, ψW=W∗ 1 γWWk=2.0171 2.2458= 0.8982. Using these values in a structural design yields design parameter values of z= [3.0443,3.0477]. Proceeding with the maximum desgin parameter, i.e. z= max izi= 3.0477, gives design reliabilities of: β= [4.3065,4.3000] which demonstrates a successful calibration to βT= 4.3. 4.1.4. Matrix Inversion Now we utilize the proposed matrix inversion approach as described in Section 3.1. Using the individual cali- brations for each load case (see Table 1), renders Equation (16) in the following form:  0 0 .6737 0.9741 0   ψ1 ψ2  =  0.6051 0.9077   from which the solution for the load combination factors is easily found to be: ψ= [0.9318 0 .8982] 10 which are precisely the same as for the coefficient method. This is expected for n= 2. Of course, the reliabilities achieved in a structural design using these values will also be the same. 4.1.5. Closed-Form Equation To calculate the result for the closed-form equation (Equation (23)) we consider the usual case of a single structural design. Consequently, from the design points for resistance and dead load in each loadcase, we take the critical values of R∗= 0.6553andG∗= 1.0371respectively. These give the partial factors as ϕ= 0.6553/0.7737 = 0.847andγG= 1.0371/1.0000 = 1 .0371. The partial factors for the remaining variables are the same as for the coefficient method. Using Equation (23), and taking the design parameter as z= max( z1, z2) = 3 .0477we establish the excess load term: S+= 0.6×1.62360 .3×2.2458 + 0.4×1.0371−3.0477×0.6553 = 0.06639 from which the load combination factors are: ψQ= 1−S+ 1×0.6×1.6236= 0.9318 ψW= 1−S+ 1×0.3×2.2458= 0.9015. These are practically the same as those from the coefficient and matrix methods, and match the results given in Sørensen (2004, p. 191) to two decimal places. Of course, this is not surprising, as the n= 2case has just one additional degree of freedom. The slight difference is because we have worked with a single resistance design point, as explained in ??. As before, we perform a structural design with these values. Interestingly, in contrast with that for the coefficient method, this yields z= 3.0477for both loadcases, indicating an optimum result. The achieved reliabilities are: β= [4.3065,4.3000] which indicates a very well calibrated design to βT= 4.3. 4.2. Three Varying Loads 4.2.1. Description Here we extend the previous example and demonstrate that the coefficient approach does not yield a single set of load combination factors; a problem resolved using the new method. The limit state function for the new problem is: g=zR−(0.2G+ 0.6Q1+ 0.35Q2+ 0.25Q3) (25) and the variables are described in Table 2. For the calibration of the design parameter z, a target reliability of βT= 4.8is used. The design points for the three load cases in which each time-varying load are considered as a combination in turn are shown in Table 2. 11 Table 2: Random variables of Example 2. DistributionAnnual Max. Point-In-Time Characteristic Design Points E[X] CoV E[ X] CoV Percentile xk x∗ 1 x∗ 2 x∗ 3 z- - - - - - - 3.5045 3.4546 3.3951 RLognormal 1.0 0.15 - - 0.05 0.7738 0.6194 0.6137 0.6124 GNormal 1.0 0.10 - - 0.50 1.0000 1.0194 1.0202 1.0207 Q1Gumbel 1.0 0.20 0.887 0.183 0.95 1.3732 1.8722 1.4497 1.5489 Q2Gumbel 1.0 0.30 0.828 0.278 0.95 1.5597 1.2591 1.1270 1.3686 Q3Gumbel 1.0 0.40 0.802 0.416 0.90 1.5218 1.6108 1.7667 1.8671 4.2.2. Coefficient Method Equations (12) and (13) can be used directly and individually. However, similar to Equation (16), we can write the coefficients of the time-varying loads as a matrix. Dividing the design point values by their characteristic values yields: γ1 ψ2,1γ2ψ3,1γ3 ψ1,2γ1 γ2 ψ3,2γ3 ψ1,3γ1ψ2,3γ2γ3,3  from which the partial factors are found as the diagonal elements to be: γQ= [1.3634 1 .1072 1 .2269] Consequently, dividing each column by its partial factor yields the coefficient method estimate of the load combi- nation factors:  1.0000 0 .7291 0 .8627 0.7743 1 .0000 0 .9463 0.8273 0 .7925 1 .0000  Clearly there is not a single solution for the load combination factor for each time-varying loads. Although not evident from the literature, presumably past practice has then yielded the load combination factors from these results as ψj= max iψi,jgiving: ψ= [0.8273 0 .7925 0 .9463] Proceeding with these values in a structural design, and finding the design parameter z= max izi= 3.6709yields achieved reliabilities of β= [5.0028 5 .0708 5 .1493] which are all quite in excess of βT= 4.8. The root mean squared error (RMSE) between the estimated design and target reliabilities is 0.2807. 12 4.2.3. Matrix Inversion Applying Equation (16) to the problem, and using the individual resistance design points in each row, yields:  0. 0.6044 0 .4668 1.1233 0 . 0.4668 1.1233 0 .6044 0 .   ψ1 ψ2 ψ3  =  0.8434 1.3114 1.4084  , from which the unique solution for the load combination factors is: ψ= [0.8353 0 .7777 0 .7993] . As before we undertake a structural design and find z= 3.5442. Using this design parameter as before, with these combination factors gives achieved reliabilities for each load case as: β= [4.8494 4 .9144 4 .9925] . Although these still exceed the target reliability (as desired), they do so by a lower amount than the (empirical forn > 2) coefficient method. In this case, RMSE for the estimated design and target reliabilities is 0.1324, a substantial reduction. 4.2.4. Closed-Form Equation Finally, consideringthemostcommoncaseforstructuraldesign, weselectthelargestcalibrateddesignparameter value from Table 2, and so z= max izi= 3.5045, which along with its associated resistance variable R= 0.6194, gives the resistance design point as ϕzR = 2.1463. Hence, we get the excess load term as: S+= 0.6×1.8722 + 0 .35×1.7270 + 0 .25×1.8671 −(2.1463−0.2×1.0207) = 0.2524 which gives the individual load combination factors as: ψ1= 1−S+ 2×0.6×1.8722= 0.8876 ψ2= 1−S+ 2×0.35×1.7270= 0.7912 ψ3= 1−S+ 2×0.25×1.8671= 0.7296. As before, we perform a structural design with these values and find the interesting result that zi= 3.5045 for each loadcase. That is, a least-squares optimum load combination factor has resulted (compare with Table 2). Using this, with the above combination factors gives achieved reliabilities for each load case as: β= [4.8000 4 .8641 4 .9421] . It is notable that this design yields precisely the target reliability index for the governing loadcase. Further, the remaining loadcases again exceed the target reliability as desired, but to a smaller degree than the empirical choices in the coefficient method. Indeed, the RMSE here is just 0.0900. 13 5. Discussion & Conclusions 5.1. Discussion The scientific literature on the calibration of partial and load combination factors has focussed on just one or two time-varying variables. This is certainly appropriate for code-calibration for many generic limit state functions, typical of those in building or other simple structures. Bridges on the other hand—especially railway bridges—must cope with a wide variety of time-varying loading variables and loadcases. For example, rail bridges must be designed or assessed for wind loads, train loads, braking forces, and track nosing forces; all of which vary in time. Little attention seems to have been given to the calibration of load combination factors for such cases of many n > 2 time-varying loads. For design codes this has not necessarily caused a problem, since it is relatively inexpensive to build in some conservatism at the construction stage. However, as the population of bridges around the world ages dramatically, probability-based bridge assessment offers a rational and quantitative means of ensuring safety. But in doing so, it is no longer acceptable to rely on ad-hoc methods of load combination factor calibration, such as the traditional coefficient method for n >2. The new method developed here addresses this problem. 5.2. Conclusions The semi-probabilistic methods for structural design and assessment rely on partial safety and load combination factors to capture the underlying uncertainties. The existing heuristics for estimating load combination factors yield estimates that vary per load case, thereby leading to ad-hoc measures that add significant conservatism to the corresponding design. Such unnecessary conservatism is not always desirable and is particularly problematic for the assessment of in-situ structures. These issues are more pronounced in problems with multiple time-variant loads such as design and assessment of bridges. This study proposes a novel matrix-based approach for the simultaneous satisfaction of nloadcases when there arentime-varying load components in a limit state function. We derive a closed-form expression for the calculation of load combination factors for the typical code calibration problem with a single resistance parameter. The proposed approach is generic and valid for typical code calibration problems involving a code-format that is linear with respect to the load variables. The proposed approach works effectively for limit state functions that are both linear and non-linear in their random variables. The study also proposes an alternate formulation in terms of limit state function evaluations, which is advantageous for programmed implementations of the proposed approach. Based on a comparative study of the existing and proposed method for estimating load combination factors, it is shown that the factors estimated with the proposed method yield a design closer to the target reliability index, than the current ad-hoc coefficient approach. For a case study problem with n= 3, it is found that the factors estimated with the proposed method are 52% more efficient than the existing methods. Consequently, the proposed method expressions can be used towards a more rational calibration of load combination factors with multiple time-varying load components. The incorporation of the proposed approach in structural reliability documents like the Joint Committee of Structural Safety’s Probabilistic Model Code and ISO2394 (2015) will make probabilistic methods more accessible for the structural engineering industry. 14 Appendix A. Computer Implementation This appendix describes the implementation of the method for complex limit-state functional forms. For most cases in code-calibration, the direct use of Equation (23) will be adequate. Structural reliability software such as PySTRA ((The PySTRA Developers, 2023)) often allows users to specify their own limit-state function, which must then be considered as a “black box” function. This is very powerful, but consideration must be given to the provision of potentially non-linear function of random variables. For such cases, the direct matrix solution of Equation (16) is more suited, but requires adaptation to facilitate a black-box LSF. For this purpose, we identify the set of LSF inputs, X, with the partitioning: X={R,G,S,Ω} (A.1) where the set Rrepresents resistance variables, Grepresents time-invariant load variables, Srepresents time-variant load variables, and Ωrepresents accompanying random variables (i.e. random variables that accompany the load and resistance variables as a multiplier, such as model errors). Based upon the partitioning of X, we can estimate the partial safety factors as an evaluation of the LSF as γi=−g S∗ i,i,Ω∗,X\Si,Ω= 0 Ski(A.2) where X\Si,Ω= 0indicates that all random variables except SiandΩare set to be zero, Si=S∗ i,i, and Ω = Ω∗. Note the similarity between Equation (A.2) and Equation (12). Substituting, Equation (A.2) in Equation (17), we get, Qi=−g G∗, S∗ i,i,Ω∗,X\G,Si,Ω= 0 (A.3) Using Equations (A.2) and (A.3), the matrix representation of the linear simultaneous equations, such as that in Equation (16), can be constructed and the load combination factors can be evaluated. For example, consider the following non-linear LSF with three time-variant loads g=ωRR−ωS(G+S1+S2+S3) (A.4) Here,R={R},G={G},S={S1, S2, S3}, and Ω ={ωR, ωS}. Using Equations (A.2) and (A.3) and including the ϕRterm in the LSF evaluation, the linear simultaneous design equations can be written as,  0 −g S∗ 2,2,Ω∗,X\S2,Ω= 0 −g S∗ 3,3,Ω∗,X\S3,Ω= 0 −g S∗ 1,1,Ω∗,X\S1,Ω= 0 0 −g S∗ 3,3,Ω∗,X\S3,Ω= 0 −g S∗ 1,1,Ω∗,X\S1,Ω= 0 −g S∗ 2,2,Ω∗,X\S2,Ω= 0 0   ψ1 ψ2 ψ3  ≤   g R∗,G∗, S∗ 1,1,Ω∗,S\S1= 0 g R∗,G∗, S∗ 2,2,Ω∗,S\S2= 0 g R∗,G∗, S∗ 3,3,Ω∗,S\S3= 0  (A.5) Such a representation based on LSF evaluations remains agnostic to the user specification of the LSF. This implementation is included in the open-source Python package PySTRA (The PySTRA Developers (2023)). 15