Stationary non-gaussian wind pressure simulation method and system based on maximum entropy method and moment conversion function

Abstract

This invention discloses a method and system for simulating stationary non-Gaussian wind pressure based on the maximum entropy method and moment transformation function. The method calculates the marginal probability density function (PDF) based on the k-th primitive moment of the target standard non-Gaussian wind pressure using the maximum entropy method; and uses the marginal probability density function (PDF) f Xj (x) Determine the transfer function; calculate the cross-correlation matrix based on the target cross-power spectral density matrix; divide the transformation relationship into region I and region II based on the error boundary line in the non-Gaussian correlation coefficient estimation; calculate the cross-correlation matrix ρ of the latent standard non-Gaussian vector for region I and region II respectively. G (τ) and the latent standard non-Gaussian vector cross power spectral density (PSD) matrix S G (ω); based on the latent Gaussian cross-power spectral density matrix S G (ω) Obtain Gaussian samples Z(t), and then obtain samples of non-Gaussian wind pressure vector X(t) based on the Gaussian samples Z(t). This method can provide better simulation accuracy, especially for wind pressure with strong non-Gaussian characteristics; since this method derives the Gaussian correlation function as a function of the non-Gaussian correlation function, it has higher simulation efficiency and better simulation accuracy.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of non-Gaussian wind pressure process simulation, and particularly to a stationary non-Gaussian wind pressure simulation method and system based on maximum entropy method and moment conversion function. BACKGROUND

[0002] It is well known that wind pressure simulation is helpful for time domain analysis of wind-induced structural response, and thus is of great significance to structural wind resistance design, especially to the design of curtain walls and glass curtain walls of buildings. Many field and wind tunnel measurements of wind pressure show that Gaussian distribution can be used to describe wind pressure in the windward region, while it cannot describe wind pressure in the flow separation region on the surface of a building, which usually exhibits mild to strong non-Gaussian characteristics. This non-Gaussian characteristic makes wind pressure simulation a challenging and necessary task.

[0003] At present, many scholars have solved the problem of non-Gaussian wind pressure or process simulation and proposed some methods, such as time series method, Karhunen-Loeve expansion method, correlation distortion method, and high-order spectral representation method. The correlation distortion method proposed by Grigoriu (1998) attracts more attention in these methods due to its conceptual directness and is widely used in non-Gaussian process simulation. The simulation process of this method is summarized as follows.

[0004] Step 1: Use the marginal probability distribution function (PDF) or cumulative distribution function (CDF) of each component of the target non-Gaussian process to determine the translation function, which relates the non-Gaussian variable to the latent Gaussian variable.

[0005] Step 2: In addition to the determined translation function, the relationship between the cross-correlation function (CCF) matrix of the target non-Gaussian process and the latent Gaussian process vector, referred to as the correlation distortion relationship, also needs to be determined.

[0006] Step 3: Obtain the cross-power spectral density (PSD) matrix of the latent Gaussian process using the correlation distortion relationship, and then simulate the Gaussian samples through spectral representation method (SRM) (Deodatis 1996). Non-Gaussian samples can be obtained by inputting these generated samples into the translation function.

[0007] In Step 1, the demand is a monotonic translation function usually determined according to the mapping of the CDF of non-Gaussian processes and the basic Gaussian process (e.g., Deodatis and Micaletti 2001; Shields and Deodatis. 2013). Such CDF-based translation functions are usually implicit and thus can be inconvenient in non-Gaussian simulation and extreme value estimation. In addition, it can not be available when only statistical moments are known. For this purpose, moment-based models, such as Hermite Polynomial Model (HPM) (Winterstein 1988), Johnson Transform Model (JTM) (Ma et al. 2016), and Unified HPM (UHPM) (Zhang et al. 2019; Lu et al. 2020), are widely used as translation functions in non-Gaussian process simulation.

[0008] The analytical translation functions of HPM, JTM, and UHPM can only be determined using the first four statistical moments of the target non-Gaussian process. Compared with HPM and UHPM, JTM has a larger range of application. Wu et al. (2020) compared the performance of HPM and JTM in non-Gaussian wind pressure simulation. The comparison shows that the simulation accuracy of HPM and JTM is not ideal for strong non-Gaussian wind pressure. This is mainly because the first four moments are not enough to capture the non-Gaussian characteristics of these wind pressures, and better estimation usually requires more high-order moments. Some researchers have proposed a simulation method based on piecewise HPM (PHPM) (Liu et al. 2020; Peng et al. 2020) to improve the accuracy and efficiency of HPM. The PHPM for standard Gaussian values less than 0 and greater than 0 is defined as HPM determined by the new first four moments related to the negative side and the positive side, respectively. The two sets of moments are calculated from the data below and above the median of the non-Gaussian data, respectively (Liu et al. 2017). However, the closed-form expression describing the correlation distortion relationship in this method is only limited to non-Gaussian processes with two new kurtosis greater than 3. To broaden the range of application, Wu et al. (2022) developed a closed-form expression describing the correlation distortion relationship of all non-Gaussian processes using piecewise JTM (PJTM). However, PHPM and PJTM are only applicable to cases where the data or PDF of the target non-Gaussian process is available. For the case where only statistical moments are given, a simulation method that is feasible in terms of accuracy and efficiency should be considered.

[0009] Simulation of non-Gaussian stationary multivariate wind pressure is of great significance to structural wind-resistant design. Moment-based translation processing methods, such as Johnson Transform Model (JTM), have been widely used in non-Gaussian simulation due to their conceptual simplicity. However, they perform poorly under wind pressure with strong non-Gaussian characteristics. In addition, the simulation efficiency is still limited by the iterative calculation of the cross-correlation function (CCF) of the potential Gaussian wind pressure. SUMMARY

[0010] Therefore, the present application aims to provide a stationary non-Gaussian wind pressure simulation method and system based on maximum entropy method and moment conversion function.

[0011] To achieve the above-mentioned purpose, the present application provides the following technical solutions:

[0012] The stationary non-Gaussian wind pressure simulation method based on maximum entropy method and moment conversion function provided by the present application comprises the following steps:

[0013] Step 1: based on the target standard non-Gaussian wind pressure X j (t) of the kth (k = 1, 2, …, M) origin moment μ jk , the marginal probability distribution function PDF is calculated by the maximum entropy method MEM;

[0014] Step 2: using the marginal probability distribution function to determine the translation function;

[0015] Step 3: according to the target cross power spectral density matrix S NG (ω), the cross-correlation coefficient matrix ρ NG (τ) is calculated;

[0016] Wherein, ρ NG (τ) represents the cross-correlation coefficient matrix; or represents the cross-cumulative distribution function CCF matrix of the target standard non-Gaussian vector;

[0017] Step 4: according to the error boundary line in the non-Gaussian correlation coefficient estimation, the transformation relationship is divided into I region and II region;

[0018] The cross-correlation coefficient (CCF) matrix ρ G (τ) of the potential standard non-Gaussian vector and the cross-PSD matrix S G (ω) of the potential standard non-Gaussian vector are obtained by calculating the I region and the II region respectively;

[0019] Wherein, ρ G (τ) represents the potential Gaussian cross-correlation coefficient matrix; S G (ω) represents the potential Gaussian cross power spectral density matrix;

[0020] Step 5: based on the potential Gaussian cross power spectral density matrix S G (ω), the Gaussian sample Z(t) is obtained, and then the sample of the non-Gaussian wind pressure vector X(t) is obtained according to the Gaussian sample Z(t).

[0021] Further, the marginal probability distribution function PDF is calculated by the maximum entropy method MEM, and the specific calculation process is as follows:

[0022]

[0023]

[0024] λ0and λ k (k = 1, 2,..., M) are Lagrange multipliers.

[0025] Further, the translation function is calculated by the following formula:

[0026] determine the translation function based on moments related to x j (t) as follows:

[0027]

[0028] When , and are the mean, standard deviation and kurtosis of a new process , respectively, whose probability density function is symmetric around the median x j (t) of x jm and

[0029] When x ≥ x jm , and are the mean, standard deviation and kurtosis of another new process , respectively, whose probability density function is centered around x jm and and are the moment-based models determined by and , respectively;

[0030] JTM is used as the moment-based model. To obtain these newly defined moments, the median of x j (t) should be estimated first as follows:

[0031]

[0032] Then and can be calculated as:

[0033]

[0034]

[0035]

[0036] Similarly, and can be calculated as

[0037]

[0038]

[0039]

[0040] Further, the cross CCF matrix ρ NG (ω) of the target standard non-Gaussian vectors is obtained by Wiener-Khinchine transformation from the corresponding PSD matrix S

[0041] The cross CCF matrix ρ NG (τ) of the target standard non-Gaussian vectors is obtained by Wiener-Khinchine transformation from the corresponding PSD matrix S NG (ω);

[0042] The cross PSD matrix S G (ω) of the potential standard non-Gaussian vectors is obtained by inverse Wiener-Khinchine transformation from ρ G (τ).

[0043] Further, the cross CCF matrix ρ G (τ) of the potential standard non-Gaussian vectors can be obtained by the closed-form formula proposed for the I region and the interpolation-based scheme proposed for the II region.

[0044] Further, the sample of the non-Gaussian wind pressure vector X(t) is obtained by substituting the Gaussian sample Z(t) into the following formula:

[0045]

[0046] The application provides a stationary non-Gaussian wind pressure simulation system based on the maximum entropy method and the moment conversion function, which comprises a memory, a processor, and a computer program stored in the memory and capable of running on the processor, and the processor implements the above method when executing the program.

[0047] The application has the following advantages:

[0048] The application provides a stationary non-Gaussian wind pressure simulation method based on the maximum entropy method and the moment conversion function, which is suitable for stationary non-Gaussian wind pressure simulation problems in the case of only target wind pressure statistical moments, and can give better simulation accuracy compared with traditional simulation methods based on moments (such as JTM), especially for wind pressure with strong non-Gaussianity; since the method derives the Gaussian correlation function as a function of the non-Gaussian correlation function, it has higher efficiency.

[0049] The method develops a series of closed-form expressions for determining the potential Gaussian wind pressure CCFs for the super large regions in the Pearson diagram, and estimates the Gaussian CCFs for the remaining regions in the Pearson diagram through an interpolation-based scheme. Compared with the traditional moment-based simulation method (such as JTM), the proposed hybrid method has higher simulation efficiency and better simulation accuracy.

[0050] Additional advantages, objects, and features of the application will be set forth in part in the description which follows, and in part will become apparent to those skilled in the art upon examination of the following or can be learned by practice of the application. The objects and other advantages of the application can be realized and attained by the structure particularly pointed out in the description. BRIEF DESCRIPTION OF DRAWINGS

[0051] In order to make the purposes, technical solutions and beneficial effects of the present application clearer, the present application provides the following drawings for illustration:

[0052] Figure 1 are numerical examples of Gaussian-non-Gaussian cross correlation functions A1, A2 and A3.

[0053] Figure 2 is the error of formula (13) when estimating the non-Gaussian correlation coefficient.

[0054] Figure 3 is the function graph of h[ρ Nju (τ)] under case 1.1

[0055] Figure 4 is the function graph of h[ρ under case 1.2

[0056] Figure 5 is the function graph of h[ρ under case 1.2 Nju has three roots.

[0057] Figure 6 is the function graph of h[ρ under case 1.3

[0058] Figure 7 is the function graph of h[ρ under case 1.3 Nju has three roots.

[0059] Figure 8 is the Gaussian correlation coefficient obtained by the interpolation method.

[0060] Figure 9To estimate different Δρ G When the non-Gaussian cross-correlation function is used, the error of formula (13) is calculated for N=10.

[0061] Figure 10 A simulation program for the proposed method.

[0062] Figure 11 This refers to the wind pressure information used in this study.

[0063] Figure 12 Select a translation function for the wind pressure coefficient for multiple variables.

[0064] Figure 13 This represents the power spectral density of the corresponding multivariable Gaussian wind pressure coefficient.

[0065] Figure 14 To simulate the power spectral density of multivariable non-Gaussian wind pressure coefficient.

[0066] Figure 15 To simulate the probability density function of multivariable non-Gaussian wind pressure coefficient. Detailed Implementation

[0067] The present invention will be further described below with reference to the accompanying drawings and specific embodiments, so that those skilled in the art can better understand and implement the present invention. However, the embodiments described are not intended to limit the present invention.

[0068] Example 1

[0069] like Figure 1 As shown in the figure, the multivariate non-Gaussian wind pressure simulation method provided in this embodiment is as follows:

[0070] Y(t)=[Y1(t), Y2(t),...,Y n (t)] T As an n-variable stationary non-Gaussian process vector, its j-th (j = 1, 2, ..., n) component process Y j The mean and standard deviation (STD) of (t) are expressed as follows: and A standardized stationary non-Gaussian process vector X(t) = [X1(t), X2(t), ..., X...] with zero mean and unit standard deviation. n (t)] T It can be easily obtained using the following formula:

[0071] X j (t)=[Y j (t)-μ Yj ] / σ Yj (1)

[0072] Among them, X j(t) is the jth component of X(t); t denotes time;

[0073] Obviously, Y(t) can be easily obtained by equation (1) as long as X(t) and its first two moment vectors (i.e., mean and standard deviation) are known. Without loss of generality, the simulation target in this embodiment is X(t), whose only available information is the kth (k = 1, 2,..., M) raw moment vector μ k = [μ 1k , μ 2k ,..., μ nk ] and the cross-PSD or CF matrix, where, is the kth raw moment of X j (t).

[0074] Reconstruction of non-Gaussian probability density function PDF by maximum entropy method

[0075] Given the kth (k = 1, 2,..., M) raw moment μ j of the component process X jk (t), its corresponding can be developed by maximum entropy method (MEM) through the following optimization problem (Jaynes 1957a, 1957b):

[0076]

[0077] where H is the information entropy of X j (t); denotes the probability density function of X j (t); x denotes the random variable; k denotes the order of the raw moment; and M denotes the Mth order.

[0078] By introducing the Euler-Lagrange equation into equation (2), the maximum entropy PDF can be approximated as:

[0079]

[0080] where λ0and λ k (k = 1, 2,..., M) are the Lagrange multipliers, which can be determined by the following convex unconstrained optimization problem:

[0081]

[0082] In order to more accurately and effectively estimate the marginal probability distribution function PDF, the following convergence mechanism is used to reasonably determine the appropriate moment order compared with the traditional maximum entropy method MEM which usually uses constant low-order moments (Wang et al. 2022)

[0083] ||H m -H m-2|| / ||H m ||≤ε H (5)

[0084] where H m represents the information entropy, which can be calculated by the estimated probability density function in equation (3) ;

[0085] where M = m; ε H is the precision value of the relative error, which is recommended to be 10 -5 ;

[0086] Generally, the integral range of the maximum entropy distribution needs to be known to ensure the automatic convergence of equation (5). In this embodiment, the lower and upper limits of the integral are the real solutions of the following equations with respect to the point of interest C:

[0087]

[0088] m i is the time of linear offset; The mathematical theory behind the process of determining the integral boundary is relatively complex, so its details are omitted here for the sake of brevity and can be found in Racz et al. (2006). Once the integral limits are determined, the Lagrange multipliers λ0and λ k (k = 1, 2,..., M) in equation (3) can be obtained by the Newton method and the modified Gram-Schmidt algorithm;

[0089] Piecewise translation function based on moments

[0090] Based on the theory of translation process, X(t) can be obtained by a nonlinear memoryless transformation of n variables Z(t) under a standard Gaussian process vector, where Z(t) = [Z1(t), Z2(t),..., Z n (t)] T , as follows (Grigoriu 1998)

[0091]

[0092] where x j (t) represents the target standard non-Gaussian wind pressure; z j (t) represents the standard Gaussian wind pressure; g j represents the transfer function;

[0093] is the translation function;

[0094] and Φ are the cumulative distribution functions (CDFs) of the standard non-Gaussian and Gaussian processes, respectively;

[0095] is the inverse function of can be obtained from the translation function g(·) in equation (7) must be monotonic.

[0096] In general, cannot give an explicit expression of x j (t) as a function of z j (t), which makes the simulation of non-Gaussian processes relatively inconvenient. To simulate non-Gaussian processes more directly and conveniently, several moment-based models such as HPM and JTM can be used as the translation function. By using these models, whose parameters can only be determined by the first four moments (mean, standard deviation, skewness, kurtosis) of the target non-Gaussian process, the non-Gaussian process can be explicitly expressed as a function of a latent standard Gaussian process. However, these moment-based models perform poorly in simulating processes with strong non-Gaussian characteristics, which can be attributed to the fact that the first four moments do not well reflect the probabilistic characteristics of the non-Gaussian process. Inspired by the piecewise HPM developed in Liu et al. (2017) mainly for non-Gaussian extreme value estimation, the present embodiment uses a moment-based piecewise transfer function related to x j (t) and gives:

[0097]

[0098] where x j (t) denotes the target standard non-Gaussian wind pressure; denotes the newly defined negative-side mean; denotes the newly defined positive-side mean; P j denotes the moment-based transfer function model; denotes the newly defined negative-side standard deviation; denotes the newly defined positive-side standard deviation; denotes the newly defined negative-side kurtosis; denotes the newly defined positive-side kurtosis;

[0099] When x≤x jm , and are the mean, standard deviation, and kurtosis of a new process , respectively, whose is symmetric about the median x j of x jm (t) and

[0100] When x≥x jm , and are the mean, standard deviation, and kurtosis of another new process , respectively, whose Also around x jm and

[0101] where, and are the moment-based models determined by and respectively;

[0102] x jm denotes the median of x j (t); denotes the probability density function of x≥x jm ; denotes the probability density function of x≤x jm ; denotes the probability density function of x j (t);

[0103] In this embodiment, JTM is used as the moment-based model. To obtain these newly defined moments, the median of x j (t) should be estimated first as follows:

[0104]

[0105] Then and can be calculated as:

[0106]

[0107] Similarly, and can be calculated as

[0108]

[0109] where x jm denotes the median of x j (t); denotes the probability density function of x≥x jm ; denotes the probability density function of x≤x jm ; denotes the probability density function of x j (t).

[0110] It should be noted that the first three moments (i.e. mean, standard deviation and skewness) used to determine or are equal to 0, 1 and 0 respectively.

[0111] Correlation function transformation relationship

[0112] To simulate a stationary non-Gaussian process vector using the moment-based piecewise translation function described above, the CCF matrix or PSD matrix in the potential stationary Gaussian process vector should first be obtained based on the cross-correlation function (CCF) matrix or cross-power spectral density (PSD) matrix of the target non-Gaussian process vector. Therefore, establishing the transformation relationship between non-Gaussian and Gaussian CCFs is extremely important, including the transformation of CCF from a Gaussian process to a non-Gaussian process and vice versa. Simply put, the former and the latter are referred to as Gaussian to non-Gaussian and non-Gaussian to Gaussian CCF transformation relationships, respectively, as follows:

[0113] Gaussian to non-Gaussian CCF transformation relationship

[0114] According to equation (2), x j (t) and x u The cross-correlation function R between (t)(u=1,2,...,n) Nju (τ) (τ is time lag) and z j (t) and z u The cross-correlation function ρ between (t) Gjk (τ) can be expressed by the following expression (e.g., Grigoriu 1998).

[0115]

[0116] Among them, z j (t) represents the j-th variable of standard Gaussian wind pressure; z u (t) represents the u-th component of the standard Gaussian wind pressure; g j Represents the transfer function of the j-th component; g u The transfer function representing the u-th component;

[0117] This represents the joint PDF of the bivariate standard Gaussian variables.

[0118] In this embodiment, R Nju (τ) is equal to the corresponding cross-correlation coefficient ρ. Nju (τ), because x j (t) and x u (t) are all standardized. Therefore, for the sake of brevity, ρ will be used in the following analysis. Nju (τ) instead of R Nju (τ).

[0119] By using Mehler's formula, equation (12) can be approximated as (Fan et al. 2023).

[0120]

[0121] where I l，s (l = j, u) is given by

[0122]

[0123] where N denotes the order; g l denotes the transfer function; H s denotes the s-th order Hermite polynomial; denotes the probability density function of a standard Gaussian variable; μ l denotes the mean value of x l (t);

[0124] μ l = μ l1 and are the mean value and the standard deviation of x l (t), respectively; is the PDF of a standard Gaussian variable; H s (·) is the s-th (s = 0, 1, 2,...) order Hermite polynomial, which can be calculated as

[0125]

[0126] Using the Gauss-Hermite quadrature formula, I l，s in equation (14) can be further approximated as

[0127]

[0128] d is the number of abscissas; z gl，a (a = 1, 2,..., d) and w gl，a are the abscissas and weights of the Gauss-Hermite quadrature formula, respectively.

[0129] Equation (13) makes the calculation of equation (12) in the forward direction effective. The accuracy of equation (13) at N = 3 is investigated here through numerical examples. First, three numerical examples denoted by A1, A2 and A3 are used. In each numerical example, the PDF of the non-Gaussian process is determined by JTM using the first four moments. Table 1 summarizes the first four moments of x j (t) and x u (t) for the three numerical examples. For each numerical example, the non-Gaussian correlation coefficients are compared as a function of the Gaussian correlation coefficients of equation (13) with the corresponding coefficients of the numerical integration in equation (12). Figure 1 The results show that the accuracy of equation (13) at N = 3 is satisfactory in the three cases. To verify its accuracy more comprehensively,

[0130] Figure 2 ​In equation (a), the error of equation (13) is given, where N = 3 is used as a function of skewness and kurtosis to estimate the non-Gaussian correlation coefficient. The error is defined as follows:

[0131]

[0132] in, and It is derived from equations (12) and (13) using ρ respectively. Gju，i =0.01i⁻¹ is the estimated non-Gaussian correlation coefficient; L = 199 is an integer;

[0133] It was observed that equation (13) can provide a sufficiently accurate estimate, where the I region is ε. ρ <2%. For those located at α 4j =(1.205α) 3j ) 1.892 +1.376 (represented by the error boundary line in this embodiment) and A very small area between (limit boundary lines) Figure 2 (a) Region II), the error is greater than 2%. When using N=4, such as Figure 2 As shown in (b) and (c), the corresponding region II has shrunk, and the error boundary line is α. 4j =(1.013α) 3j ) 2.146 +1.376. Note that the error boundary lines used to divide regions I and II depend on an acceptable error of 2% in this embodiment. When using other acceptable errors, the error boundary lines will be changed and need to be redefined.

[0134] Table 1. Numerical Examples A1 to A3

[0135]

[0136] Non-Gaussian to Gaussian CCF Transform Relationship

[0137] Equation (13) makes the calculation of equation (12) valid in the positive direction. However, the inverse calculation of equation (13), i.e., the estimation of the Gaussian CCF based on the given non-Gaussian CCF and translation function, is still not available. In order to simulate the non-Gaussian process more effectively, the inverse function of equation (13) will be given here for region I, and an interpolation scheme will be proposed for region II.

[0138] Area I

[0139] Since the formula (13) with N=4 has a larger I region than the formula (13) with N=3, this embodiment uses the formula (13) with N=4 and rewrites it as follows:

[0140]

[0141] The translation function g in equation (14) l (z l ) is monotonic, when z l tends to negative and positive infinity, g l (z l ) tends to negative and positive infinity, respectively, which indicates that in most sub-intervals of the interval (-∞, ∞),

[0142] Therefore, I j，1 > 0 and I u，1 > 0 and equation (18) can be given by:

[0143]

[0144] where a ju = I j，1 I u，1 > 0; b ju = I j，2 I u，2 / (2I j，1 I u，1 );

[0145] c ju = I j，3 I u，3 / (6I j，1 I u，1 ); d ju = I j，4 I u，4 / (6I j，1 I u，1 ).

[0146] Obviously, p Nju (t) can be a linear, quadratic, cubic or quartic polynomial function of p Gju (t) with different b ju , c ju and d ju combinations. Therefore, it is necessary to develop a monotonic inverse function for h[p Gju (t)] as follows with different b ju , c ju and d ju combinations.

[0147] Case 1: d ju ≠ 0

[0148] In this case, p Nju (t) is p GjuThe first and second derivatives of the quartic polynomial function (τ) can be derived.

[0149]

[0150]

[0151] Similar to h[ρ in case 3] Gju (τ)] and and They are ρ Gju The third-order and second-order functions. Clearly, The symbol can be derived from d ju and discriminant It is determined that the discriminant is given by the following formula.

[0152]

[0153] According to d ju and Specific combinations can be determined within certain specific ranges. and The sign of h[ρ] helps to determine the symbol of h[ρ]. Gju The shape of (τ)].

[0154] Case 1.1: ie,

[0155] In this case, it satisfies the following condition throughout the entire range (-∞, +∞). or therefore It is monotonic. When ρ Gju It tends towards negative infinity and positive infinity (denoted as ρ). Gju →-∞ and ρ Gju When →+∞), for d ju >0, respectively and For d ju <0, respectively and Therefore, equation There is a root, as shown below:

[0156]

[0157] and It can be calculated using the following formula:

[0158]

[0159] Because of d ju >0, in and in the range, and Thus h[ρ Gju (τ)] is monotonically decreasing and increasing, respectively. Similarly, for d ju < 0, h[ρ Gju (τ)] is monotonically increasing and decreasing, respectively, in the two ranges. The shape of h[ρ Gju (τ)] can be explained by Figure 9 It is seen that for d ju > 0 and d ju < 0, h -1 [ρ Nju (τ)] requires h[ρ Gju (τ)] in the range and and the solid line portion of the curve Figure 9 , respectively, and can be obtained using Ferrari's formula (see Appendix C). For d ju > 0 and d ju < 0, h -1 [ρ Nju (τ)] can be expressed as:

[0160]

[0161]

[0162] Δ f1 , Δ f2 and R can be expressed as

[0163]

[0164]

[0165]

[0166] Note that y * is an arbitrary root of the equation

[0167]

[0168] Clearly, y * can be solved using Cardano's formula, whose expression is omitted here. The applicable ranges of (25) and (26) are and two ranges, respectively.

[0169] Case 1.2:

[0170] Since equation It has the following two different solutions

[0171]

[0172] Obviously, Because of d ju <0, in and Within the range, therefore Monotonically decreasing. Within the range, therefore Monotonically increasing. When ρ Gju When →-∞ and +∞, And -∞ are its asymptotes. Additionally, Through the point (0, a) ju ).according to and The combination of these factors leads to the following conclusions: Graphics, such as Figure 4 As shown. It can be seen that, for Figure 4 The cases in (b), (d), (f), and (g) are as follows: There is a root; for Figure 4 The cases in (a), (c) and (e) are as follows: There are three roots.

[0173] for Figure 4 The cases in (b), (d), (f), and (g) are as follows: The root can be calculated using equation (23). and Within the range, The values ​​are >0 and <0 respectively, therefore h[ρ Gju [τ] are monotonically increasing and monotonically decreasing, respectively. Gju The shape of (τ)] and Figure 6 The similarity to (b) in [the text]. It can be seen that h[ρ] Gju [τ] has two roots, which need to be determined. Within the range, h is obtained by equation (26). -1 [ρ Nju (τ)。

[0174] for Figure 4 For cases (a), (c), and (e), Cardano's formula can be used to calculate... The three roots are shown below:

[0175]

[0176] in, Please note,

[0177] In and range, and h[ρ Nju (τ)] can be illustrated by Figure 5 For the cases of (a) and (d) in Figure 5 h[ρ Gju (τ)] has four roots, and h Gju [ρ -1 (τ)] is obtained from h[ρ Nju (τ)] in the range of and This can be obtained using Ferrari's formula, respectively

[0178]

[0179]

[0180] The applicable range of formula (33) and (34) is, respectively:

[0181] and

[0182] For the cases of (b), (c) and (e) in Figure 5 h[ρ Gju (τ)] has two roots, and h Gju [ρ -1 (τ)] is obtained from h[ρ Nju (τ)] in the solid part, which can be expressed as

[0183]

[0184] For the cases given in (b) and (c) in Figure 5 the applicable range of formula (35) is:

[0185] For the case given in (e) in Figure 5 the applicable range of formula (48) is

[0186] Case 1.3:

[0187] Similar to Case 4.2, the two roots of h can be calculated by formula (31). Since d ju > 0, in and Within the range, therefore Increase. In Within the range, therefore Decrease. When ρ Gju When →-∞ and +∞, And +∞. Similarly, according to and The combination can be made by Figure 6 illustrate The shape. It can be seen that, for Figure 6 The cases in (b), (c), (d), and (f) are as follows: There is a root; for Figure 6 The cases in (a), (e), and (g) are as follows: There are three roots.

[0188] for Figure 6 The cases in (b), (c), (d), and (f) can be calculated using equation (23). The root. In and Within the range, and Therefore h[ρ Gju (τ)] Decreases and increases respectively. h[ρ] Gju The shape of (τ)] is similar to Figure 3 (a) shows that h[ρ] Gju [τ] has two roots, and it is necessary to know that in h[ρ] within the range Gju Only then can we obtain h, which can be expressed by equation (25). -1 [ρ Nju (τ)。

[0189] Similarly, for Figure 6 The cases in (a), (e), and (g) are as follows: The three roots and It can be calculated using equation (32). It should be noted that... exist and Within the range, respectively and h[ρ Nju The shape of (τ)] can be determined by Figure 7 illustrate.

[0190] for Figure 7 In cases (a) and (c), h[ρ Gju(τ) has four roots and h and h Gju (τ) in the range -1 [ρ Nju (τ)] can be obtained as follows:

[0191]

[0192]

[0193] The applicable ranges of equation (36) and equation (37) are: and

[0194] For the three cases of (b), (d) and (e) in Figure 7 h Gju (τ) has two roots, and h -1 [ρ Nju (τ)] is obtained using h Gju (τ) represented by the solid line, and h -1 [ρ Nju (τ)] is as follows:

[0195]

[0196] For the two cases of (d) and (e) in Figure 7 the applicable range of equation (38) is:

[0197] and for the case of (b) in Figure 7 the applicable range of equation (38) is

[0198] For ease of application, the equation h -1 [ρ Njk (τ)] for describing the non-Gaussian to Gaussian CCF transformation relationship of the I region in all cases is summarized in Table 2.

[0199] Table 2 inverse function h -1 [ρ Njk (τ)] of equation (19) for I region

[0200]

[0201]

[0202] Note: equation (15) is used to calculate h and Equation (31) computes the Case 4 in Equation (30) and In Case 4, when respectively, use Equation (23) and Equation (32) to compute Equation (32) computes the Case 4 in Equation (30) and

[0203] where, denote the solutions of Equation (20) for different cases, respectively; denote the solutions of Equation (21) for different cases, respectively; denote the discriminant of solving Equation (21); j，1 denote the integrals of the transfer function for j components; u，1 denote the integral of the transfer function for the u-th component; f1 , Δ f2 denote the intermediate parameters when solving Ferrari's equation, respectively; denote the intermediate parameter of Equation (23);

[0204] Region II

[0205] As mentioned above, for the CCF estimation of a non-Gaussian process located in Region II, the estimation performance of Equation (13) with N = 4 is poor. In order to obtain more accurate estimation, N needs to be greater than 4. However, when N > 4, it becomes extremely difficult to obtain the inverse function of Equation (13). The present embodiment proposes an interpolation-based scheme for obtaining the Gaussian CCF from a given non-Gaussian CCF with N > 4. The interpolation-based scheme for estimating the Gaussian CCF is described as follows:

[0206] Step 1: Determine the appropriate Gaussian correlation coefficient increment Δρ G and ρ Gju = -1 + iΔρ G (i = 0, 1, 2,..., L), where L is a positive integer. It is noted that the maximum value of ρ Gju needs to be taken as 1, thus, when 2 / Δρ G are integer and decimal, respectively, L = [2 / Δρ G ] and L = [2 / Δρ G ]+1, where [.] denotes the upward rounding symbol.

[0207] Step 2: Substitute ρ Gju into Equation (13) with a larger N value to obtain the corresponding non-Gaussian correlation coefficient ρ Nju = h(ρ Gju ). In particular, let and For most Region II, equation (13) with N=10 performs well, so it is recommended to use N=10.

[0208] Step 3: Based on the given non-Gaussian correlation coefficient Data [ρ] can be used Nju ,h(ρ Gju And advanced interpolation methods (such as linear interpolation or spline interpolation) are used to estimate the corresponding Gaussian correlation coefficient ρ. Gju，intp .

[0209] In interpolation-based schemes, it is necessary to determine Δρ. G Interpolation methods. Linear interpolation and spline interpolation are widely used for interpolation due to their simplicity and accuracy; these two methods can be considered. The error caused by estimating the Gaussian correlation coefficient using an interpolation-based scheme is defined as:

[0210]

[0211] in, and The i-th Gaussian correlation coefficient is estimated by using the iterative equation (13) and the interpolation scheme proposed in this paper.

[0212] Here, we study the performance of these two interpolation methods through two examples, in which two standard non-Gaussian processes are used, with skewness of 1.92 and 1.3 and kurtosis of 5.2 and 3.0; the probability density function of the non-Gaussian process is determined by JTM using the first four moments.

[0213] Figure 8 Comparison of Δρ G =0.1, The Gaussian correlation coefficients obtained by the two interpolation methods are used as Gaussian correlation coefficient functions. For Figure 8 (a) Compared with linear interpolation, spline interpolation can provide more accurate estimation results, among which... However, in Figure 8 In (b), due to overfitting at the negative tail, the estimation results given by the spline interpolation method are poor. Therefore, this embodiment recommends using linear interpolation.

[0214] To determine Δρ G , Figure 9 The text demonstrates the use of linear interpolation for Δρ. G When the values ​​are 0.01, 0.1, and 0.2 respectively. The changes with skewness and kurtosis. Observation shows that when Δρ G When the value is 0.01, the interpolation-based scheme can provide sufficiently accurate estimation results for almost all regions II, where Therefore, it is suggested to use Δρ G = 0.01.

[0215] As Figure 10 shown, the overall simulation process, the detailed steps of the process of the proposed method are as follows:

[0216] Step 1: based on the target standard non-Gaussian wind pressure X j (t) of the k (k = 1, 2,..., M) original matrix μ jk , that is, the j (j = 1, 2,..., n) component of the vector X(t), the edge probability distribution function PDF is calculated by the maximum entropy method MEM; The specific calculation process is estimated using equations (3) and (4);

[0217] Step 2: use the edge probability density function to determine the transfer function, and the specific process can be used to determine the transfer function by equations (8)-(11);

[0218] Step 3: according to the target cross power spectral density matrix S NG (ω) to calculate the cross correlation coefficient matrix ρ NG (τ);

[0219] Where, ρ NG (τ) represents the cross correlation coefficient matrix;

[0220] In this embodiment, the cross CCF matrix of the target standard non-Gaussian vector ρ NG (τ) is obtained by Wiener-Khinchine conversion from the corresponding PSD matrix S

[0221] ρ G (τ) represents the cross CCF matrix of the potential standard non-Gaussian vector;

[0222] ρ Nju (τ) represents the cross correlation coefficient between z j (t) and z u (t);

[0223] R Nju (τ) represents the cross correlation function between x j (t) and x u (t);

[0224] S NG (ω) represents the target cross power spectral density matrix;

[0225] The cross CCF matrix of the target standard non-Gaussian vector ρ NG (τ) in this embodiment can be obtained by Wiener-Khinchine conversion from the corresponding PSD matrix S NG (ω);

[0226] Step 4: Divide the transform relationship into I region and II region according to the error boundary line in the non-Gaussian correlation coefficient estimation;

[0227] Calculate I region and II region respectively to obtain the cross CCF matrix p G (τ) and the cross PSD matrix S G (ω) of potential standard non-Gaussian vectors;

[0228] The cross CCF matrix p G (τ) of potential standard non-Gaussian vectors in this embodiment can be obtained by the closed-form formula proposed in Table 2 for I region and the interpolation-based scheme for II region.

[0229] The cross PSD matrix S G (ω) of potential standard non-Gaussian vectors in this embodiment can be obtained from p G (τ) by inverse Wiener-Khinchine transform.

[0230] where p G (τ) represents the potential Gaussian cross-correlation coefficient matrix; S G (ω) represents the potential Gaussian cross-power spectral density matrix.

[0231] Step 5: Obtain Gaussian samples Z(t) based on the potential Gaussian cross-power spectral density matrix S G (ω), and then obtain samples of non-Gaussian wind pressure vector X(t) according to Gaussian samples Z(t).

[0232] The details of the spectral representation method SRM in this embodiment can be found in Deodatis (1996);

[0233] In this embodiment, the Gaussian samples Z(t) are substituted into equation (8) to obtain samples of non-Gaussian wind pressure vector x(t).

[0234] Numerical example

[0235] As shown in Figure 11 , wind pressure based on roof wind tunnel test is used to evaluate the method proposed in this embodiment. Figure 11 As shown in (a) of Figure 11 , in the wind tunnel test, the model scale is 1:100, and the sampling frequency is 312.5 Hz. The total sampling time is 55 minutes. There are a total of 265 measurement points on the roof (marked with "+"), as shown in (b). The wind pressure at each measurement point is indicated at 0°, 45°, and 90°. More detailed information about this test can be found in Liu et al. (2017).

[0236] For simplicity, wind pressure with wind direction of 90° is used in the following analysis. To verify the effectiveness of the proposed method, two numerical examples are considered, including simulation of univariate and multivariate non-Gaussian wind pressure. The measurement points used are marked with red circles on the roof as Figure 11 (b). For comparison, the simulation results of JTM are also given. Detailed information about JTM can be found in Wu et al. (2020), which is omitted here for brevity.

[0237] Multivariate wind pressure simulation

[0238] In this section, the wind pressures at Taps 41, 48 and 49 are used as an example of multivariate non-Gaussian simulation. Again, the first ten origin moments of these normalized wind pressure coefficients are shown in Table 3. The results show that the wind pressure coefficients at Taps 41 and 49 are located in the SU system, exhibiting mild and strong non-Gaussianity, respectively, while the wind pressure coefficient at Tap 48 is located in the S B system, exhibiting strong non-Gaussianity.

[0239] Table 3. Origin moments of non-Gaussian wind pressure coefficients

[0240]

[0241] From the data given in Table 3, the translation functions of JTM and the proposed method can be obtained. Figure 12 The translation functions of JTM and the proposed method as well as the translation functions from the data are shown. Similar to the observation in the univariate simulation case, the translation functions of the proposed method are very close to the target values for all wind pressure coefficients. However, for the wind pressure with strong non-Gaussianity, JTM significantly reduces the target values of the positive tail. Using the translation functions, the corresponding Gaussian correlation coefficients can be estimated. Similar to the wind pressures in the univariate simulation case, the wind pressures used in this section are located in the I region, so the Gaussian correlation coefficients can be estimated using the closed-form formula proposed. The corresponding auto-power spectral densities and cross-power spectral densities can be easily obtained. For brevity, only the auto-power spectral density of the corresponding Gaussian wind pressure at Tap 41 and the cross-power spectral density of the corresponding Gaussian processes at Taps 48 and 49 are provided, as Figure 13 shown. This shows that JTM and the proposed method can give accurate estimates. In this multivariate simulation case, a total of 9 Gaussian correlation coefficients at each time lag need to be calculated. To obtain these coefficients, the CPU time is 68 seconds and 1.4 seconds, respectively, when calculated using the iterative and the proposed analytical formula.

[0242] Using the JTM and the proposed method, 100 samples of normalized wind pressure coefficients were simulated. The auto and cross power spectral densities of these samples can be easily estimated. For brevity, only the ensemble averages of the auto power spectral densities of the simulated wind pressure at Tap 41 and the cross power spectral densities of the simulated wind pressure at Taps 48 and 49 are provided, as shown in Figs. 6 and 7, respectively. This shows that both methods yield auto and cross power spectral densities very close to the target values. Figure 14 The probability density functions (PDFs) of the simulated samples by the JTM and the proposed method were compared. Similar to the observations in the univariate simulation case, the proposed method can provide satisfactory estimates of the probability density functions for non-Gaussian wind pressures. However, for wind pressures with strong non-Gaussianity, the JTM generally underestimates the values in the positive tail of the probability density function. Figure 15 The probability density functions (PDFs) of the simulated samples by the JTM and the proposed method were compared. Similar to the observations in the univariate simulation case, the proposed method can provide satisfactory estimates of the probability density functions for non-Gaussian wind pressures. However, for wind pressures with strong non-Gaussianity, the JTM generally underestimates the values in the positive tail of the probability density function.

[0243] The above-described embodiments are merely preferred embodiments of the present application, and the scope of the present application is not limited thereto. Any equivalent alternatives or modifications of the present application made by those skilled in the art based on the technical concept of the present application belong to the scope of the present application. The scope of the present application is defined by the appended claims.

Claims

1. A stationary non-Gaussian wind pressure simulation method based on the maximum entropy method and the matrix transfer function, for simulating the structural wind resistance design of a building, characterized in that: The method comprises the following steps: Step 1: Target criterion based non-Gaussian wind pressure of the first original matrix , The marginal probability distribution function PDF is calculated by the maximum entropy method MEM. Step 2: Using the edge probability density function, PDF determining a transfer function; Step 3: Compute the target cross-power spectral density matrix Compute the cross-correlation coefficient matrix ; wherein denotes the cross-correlation matrix; Step 4: dividing the transform relationship into I region and II region according to the error boundary line in non-Gaussian correlation coefficient estimation; The intercorrelation coefficient CCF matrix of the potential standard non-Gaussian vectors is calculated for region I and region II respectively and the interpower spectral density PSD matrix of the potential standard non-Gaussian vectors ; wherein, represents a potential Gaussian cross-correlation coefficient matrix; represents a potential Gaussian cross-power spectral density matrix; Step 5: Obtain the potential Gaussian cross-power spectral density matrix Obtain Gaussian samples Then obtain non-Gaussian wind pressure vector samples from the Gaussian samples ​ The edge probability density function is calculated by a maximum entropy method (MEM), and the specific calculation process is as follows: wherein, denotes an edge probability distribution function; denotes a random variable; denotes the order of the origin moment; denotes the order; and is a Lagrange multiplier, ; The transfer function is calculated by the following formula: determining a segment translation function based on a matrix related to a matrix-based segment translation function: wherein, denotes the target standard non-Gaussian wind pressure; denotes the newly defined negative side mean; denotes the newly defined positive side mean; denotes the moment-based transfer function model; denotes the newly defined negative side standard deviation; denotes the newly defined positive side standard deviation; denotes the newly defined negative side kurtosis; denotes the newly defined positive side kurtosis; the newly defined negative side mean , negative side standard deviation , negative side kurtosis is calculated according to the following formula: ； The , positive side standard deviation , positive side kurtosis is calculated according to the following formula: ； wherein denotes the median of denotes the probability density function of 2. The stationary non-Gaussian wind pressure simulation method based on the maximum entropy method and the moment conversion function according to claim 1, characterized in that: The translation function is obtained by Wiener-Khinchine transformation from the corresponding PSD matrix obtained; Cross CCF matrix of target standard non-Gaussian vectors from the corresponding PSD matrix by Wiener-Khinchine transform obtained; Cross-PSD matrix of potential standard non-Gaussian vectors obtained by inverse Wiener-Khinchine transform from is obtained.

3. The stationary non-Gaussian wind pressure simulation method based on the maximum entropy method and the moment conversion function according to claim 1, characterized in that: the cross-CCF matrix of the potential standard non-Gaussian vectors It can be obtained by the closed-form formula proposed for region I and the interpolation-based scheme proposed for region II.

4. The method for simulating stationary non-Gaussian wind pressure based on maximum entropy method and moment conversion function according to claim 1, characterized in that: The non-Gaussian wind pressure vector The obtaining of the sample of the Gaussian sample is calculated by substituting the following formula: wherein, denotes the target standard non-Gaussian wind pressure; denotes the newly defined negative side mean; denotes the newly defined positive side mean; denotes the newly defined negative side skewness; and the newly defined positive side skewness; denotes the newly defined negative side kurtosis; denotes the newly defined positive side mean; denotes the newly defined negative side kurtosis; denotes the newly defined positive side kurtosis.

5. A stationary non-Gaussian wind pressure simulation system based on the maximum entropy method and the matrix transfer function, comprising a memory, a processor and a computer program stored in the memory and executable on the processor, characterized in that, The processor implements the method in any one of claims 1 to 4 when executing the program.

Images