METHOD FOR CONSTRUCTING GENERAL PROBABILITY MODEL OF HARMONIC EMISSION LEVEL FOR INDUSTRIAL LOAD

Abstract

A method for constructing a general probability model of a harmonic emission level for an industrial load is provided. The method establishes, based on harmonic data monitored by a power quality monitoring system, a general probability model by combining a parametric estimation method based on a normal distribution function and a lognormal distribution function with a nonparametric estimation method represented by a kernel density estimation method, taking a degree of approximation between the general probability model and an actual probability distribution of each harmonic current as an objective function based on parameters required by the general probability model, and optimizing and solving the parameters of the proposed general probability model by using a multiplier method to determine parameters of the general probability model to finally obtain a general probability model applicable to different industrial loads.

Claims

1. A method for constructing a general probability model of a harmonic emission level for an industrial load, comprising the following steps: step 1: extracting harmonic monitoring data of the industrial load to obtain a harmonic characteristic dataset X of a user: $X=\left[\begin{array}{cccc}{I}_2^1& I_3^1& L& I_{25}^1\\ I_2^2& I_3^2& L& I_{25}^2\\ M& M& M& M\\ I_2^N& I_3^N& L& I_{25}^N\end{array}\right],$ wherein N represents a total quantity of sampling points, each column vector in the X represents a harmonic current monitoring sequence, a subscript of the I represents a harmonic order, and a superscript of the I represents a quantity of sampling sequences; step 2: constructing the general probability model for an h-th harmonic I.sub.h in the harmonic characteristic dataset: $f\left(I_h\right)=\underset{i=1}{\overset{3}{.Math.}}{\lambda }_i{f}_i\left(I_h\right),$ wherein f.sub.i(.) represents a probability density subfunction and λ.sub.i represents a weight coefficient of the probability density subfunction; f.sub.1(.) represents a part that is of the I.sub.h and obeys a normal distribution, f.sub.2(.) represents a part that is of the I.sub.h and obeys a lognormal distribution, and f.sub.3(.) represents a part that is of the I.sub.h and obeys another distribution; and the f.sub.i(.) is expressed as follows: $f_1\left(I_h\right)=\frac{1}{\sqrt{2\pi }{\sigma }_1}{e}^{-\frac{\left(I_h-{\mu }_1\right)}{2{\sigma }_1^2}}{f}_2\left(I_h\right)=\frac{1}{\sqrt{2\pi }{I}_h{\sigma }_2}{e}^{-\frac{{\left({\mathrm{lnI}}_h-{\mu }_2\right)}^2}{2{\sigma }_2^2}}{f}_3\left(I_h\right)=\frac{1}{\mathrm{nb}}\underset{j=1}{\overset{n}{.Math.}}K\left(\frac{I_h-I_h^j}{b}\right)$ wherein μ.sub.1 and μ.sub.2 represent mathematical expectations of the probability density subfunction; σ.sub.1 and σ.sub.2 represent standard deviations of the probability density subfunction; K(.) represents a kernel function, b>0, wherein b represents a smoothing parameter, which is referred to as a window; I.sub.h.sup.j represents a j.sup.th sample of the I.sub.h in each window; and n represents a total quantity of samples in each window; the weight coefficient of the probability density subfunction meets the following formula: $\underset{i=1}{\overset{3}{.Math.}}{\lambda }_i=1\left(0\le {\lambda }_i\le 1\right),$ wherein λ.sub.1=1 represents that the I.sub.h obeys a single normal distribution; and A2=1 represents that the I.sub.h obeys the lognormal distribution; step 3: discretizing the general probability model of the I.sub.h, wherein the I.sub.h is discretized to discretize the f.sub.1(.) and the f.sub.2(.) to obtain the following discretized general probability model: $f\left({\tilde{I}}_h\right)=\underset{i=1}{\overset{3}{.Math.}}{\lambda }_i{f}_i\left({\tilde{I}}_h\right),{\tilde{I}}_h=0:0.01:\max \left(I_h\right),$ wherein max (I.sub.h) represents a maximum current value of the h-th harmonic; step 4: constructing a parameter optimization model of the general probability model, which specifically comprises the following substeps: step 4.1: constructing an objective function, wherein a degree of approximation between the general probability model and an actual probability distribution of the I.sub.h is reflected by a difference between a mathematical expectation calculated by a general optimization model and an actual value, as well as a difference between a standard deviation calculated by the general optimization model and an actual value, wherein the objective function is as follows: $\min y_1={\left(\underset{i=1}{\overset{3}{.Math.}}{\lambda }_i{E}_i\left({\tilde{I}}_h\right)-E_0\left({\tilde{I}}_h\right)\right)}^2\min y_2={\left(\sqrt{\underset{i=1}{\overset{3}{.Math.}}{\lambda }_i^2{\mathrm{Var}}_i\left({\tilde{I}}_h\right)}-\sqrt{{\mathrm{Var}}_0\left({\tilde{I}}_h\right)}\right)}^2$ wherein y.sub.1 and y.sub.2 respectively represent mean square errors of a mathematical expectation and a standard deviation of the general probability model; E.sub.i(Ī.sub.h) and E.sub.0(Ī.sub.h) respectively represent a mathematical expectation that is of the I.sub.h and calculated by the probability density subfunction, and an actual mathematical expectation of the I.sub.h; and Var.sub.i(Ī.sub.h) and Var.sub.0(Ī.sub.h) respectively represent a standard deviation that is of the I.sub.h and calculated by the probability density subfunction, and an actual standard deviation of the I.sub.h; and min y.sub.1 and min y.sub.2 are combined into a minimum objective function, and the combined minimum objective function is as follows: $\min y=\frac{y_1+y_2}{2}$ step 4.2: determining constraints, wherein the constraints comprise an equality constraint and an inequality constraint, wherein 1) an equality constraint for optimizing the weight coefficient λ.sub.i of the probability density subfunction is determined according to the following formula, and is expressed by l: $l=\underset{i=1}{\overset{3}{.Math.}}{\lambda }_i-1=0$ 2) the inequality constraint comprises a value range of the weight coefficient λ.sub.i and a value range that is of the random variable I.sub.h and determined by numerical characteristics (μ.sub.1, σ.sub.1) and (μ.sub.2, σ.sub.2) when the single probability density subfunction takes effect; an inequality constraint of the λ.sub.i is as follows: $\{\begin{array}{c}{g}_1={\lambda }_1\ge 0\\ g_2={\lambda }_2\ge 0\\ g_3={\lambda }_3\ge 0\end{array},$ wherein assuming that 95% confidence intervals of {μ.sub.1, μ.sub.2, σ.sub.1, σ.sub.2} are [{circumflex over (θ)}.sub.1, {circumflex over (θ)}.sub.2], [{circumflex over (θ)}.sub.3, {circumflex over (θ)}.sub.4], [{circumflex over (θ)}.sub.5, {circumflex over (θ)}.sub.6], and [{circumflex over (θ)}.sub.7, {circumflex over (θ)}.sub.8] respectively, inequality constraints of the optimization variables {μ.sub.1, μ.sub.2, σ.sub.1, σ.sub.2} are as follows: $\{\begin{array}{c}{g}_4={\mu }_1-{\hat{\theta }}_1\ge 0\\ g_5={\hat{\theta }}_2-{\mu }_1\ge 0\end{array},\{\begin{array}{c}{g}_6={\mu }_2-{\hat{\theta }}_3\ge 0\\ g_7={\hat{\theta }}_4-{\mu }_2\ge 0\end{array},\{\begin{array}{c}{g}_8={\sigma }_1-{\hat{\theta }}_5\ge 0\\ g_9={\hat{\theta }}_6-{\sigma }_1\ge 0\end{array},\mathrm{and}\{\begin{array}{c}{g}_{10}={\sigma }_2-{\hat{\theta }}_7\ge 0\\ g_{11}={\hat{\theta }}_8-{\sigma }_2\ge 0\end{array},$ wherein g.sub.q represents the inequality constraint and q=1, 2, . . . , 11; step 5: solving parameters {λ.sub.1, λ.sub.2, λ.sub.3, μ.sub.1, μ.sub.2, σ.sub.1, σ.sub.2} of the general probability model, wherein a constrained problem is converted into an unconstrained problem, and a multiplier method is used for solving, that is, an optimization variable set is defined as γ={λ.sub.1, λ.sub.2, λ.sub.3, μ.sub.1, μ.sub.2, σ.sub.1, σ.sub.2}, and an augmented Lagrange function is defined as J and is expressed is as follows: $J\left(\gamma ,\omega ,\nu ,\rho \right)=y\left(\gamma \right)-\nu l\left(\gamma \right)+\frac{\rho }{2}{l}^2\left(\gamma \right)+\frac{1}{2\rho }\underset{q=1}{\overset{11}{.Math.}}\left\{{\left[\max \left(0,{\omega }_q-\rho g_q\left(\gamma \right)\right)\right]}^2-{\omega }_q^2\right\},$ wherein y(γ) represents the objective function, l(γ) represents the equality constraint, g.sub.q(γ) represents the inequality constraint, ω.sub.q represents a Lagrange multiplier of the inequality constraint, and ν represents a Lagrange multiplier of the equality constraint; and for the J(γ, ω, ν, φ, the sufficiently large parameter ρ is taken, and the multipliers ω and ν are continuously corrected to obtain a local optimal solution by minimizing the J(γ, ω, ν, φ, wherein correction formulas of the multipliers ω and ν are as follows: $\{\begin{array}{c}{\omega }_q^{\left(k+1\right)}=\max \left(0,{\omega }_q^{\left(k\right)}-\rho g_q\left({\gamma }^{\left(k\right)}\right)\right),q=1,2,.Math.,11\\ v^{\left(k+1\right)}=v^{\left(k\right)}-\rho l\left({\gamma }^{\left(k\right)}\right)\end{array},$ wherein k in a superscript represents a quantity of corrections; and step 6: obtaining the general probability model of the I.sub.h.

2. The method for constructing the general probability model of the harmonic emission level for the industrial load according to claim 1, wherein in the objective function in step 4.1, $E_1\left({\tilde{I}}_h\right)={\mu }_1,{\mathrm{Var}}_1\left({\tilde{I}}_h\right)={\sigma }_1^2,E_2\left({\tilde{I}}_h\right)=e^{{\mu }_2+{\sigma }_2^2/2},{\mathrm{Var}}_2\left({\tilde{I}}_h\right)=\left(e^{{\sigma }_2^2}-1\right)e^{2{\mu }_2+{\sigma }_2^2},E_3\left({\tilde{I}}_h\right)=\underset{j=1}{\overset{N_1}{.Math.}}{\tilde{I}}_h^{j}p\left({\tilde{I}}_h\right),{\mathrm{Var}}_3\left({\tilde{I}}_h\right)=\underset{j=1}{\overset{N_1}{.Math.}}{\left({\tilde{I}}_h^j-E_3\left({\tilde{I}}_h\right)\right)}^{2}p\left({\tilde{I}}_h^j\right),N_1=\mathrm{length}\left({\tilde{I}}_h\right),\mathrm{and}p\left({\tilde{I}}_h^j\right)=\frac{f_3\left({\tilde{I}}_h^j\right)}{\mathrm{sum}\left(f_3\left({\tilde{I}}_h\right)\right)}.$

3. The method for constructing the general probability model of the harmonic emission level for the industrial load according to claim 1, wherein the multiplier method in step 5 specifically comprises the following steps: step a: setting an initial point γ.sup.(0), initial multiplier vector estimates ω.sup.(1) and ν.sup.(1), a parameter ρ, an allowable error ε (>0), a constant c (>1), β (∈(0,1)), and k (=1); step b: taking γ.sup.(k−1) as an initial point and solving an unconstrained problem shown in the following formula to obtain a solution γ.sup.(k):
min J(γ,ω.sup.(k),ν.sup.(k),ρ); step c: if ∥l(γ.sup.(k))∥<ε, stopping the calculation and obtaining a point γ.sup.(k)); otherwise, performing step d; step d: if ∥l(γ.sup.(k))∥/∥l(γ.sup.(k−1))∥≥β, setting ρ=cρ and performing step e; otherwise, performing step e directly; and step e: using a formula $\{\begin{array}{c}{\omega }_q^{\left(k+1\right)}=\max \left(0,{\omega }_q^{\left(k\right)}-\rho g_q\left({\gamma }^{\left(k\right)}\right)\right),q=1,2,.Math.,11\\ v^{\left(k+1\right)}=v^{\left(k\right)}-\rho l\left({\gamma }^{\left(k\right)}\right)\end{array}$ to correct multipliers ω.sub.q.sup.(k+1) and ν.sup.(k+1), setting k=k+1 and performing step b.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0045] FIG. 1 is a basic flowchart according to the present disclosure;

[0046] FIG. 2 shows measured data of a harmonic current; and

[0047] FIG. 3 is a basic flowchart of a multiplier method.

DETAILED DESCRIPTION OF THE EMBODIMENTS

[0048] The present disclosure will be further described below in conjunction with the accompanying drawings and specific embodiments.

[0049] An industrial load has a large capacity and accounts for a large proportion, which imposes a great impact on the power quality of a power system. To accurately depict a harmonic problem caused by the industrial load to a grid, the present disclosure provides a general probability model of a harmonic emission level for an industrial load. A basic flowchart is shown in FIG. 1, which is divided into six steps, namely, S1 to S6.

[0050] S1: Harmonic monitoring data of an industrial load is extracted to obtain harmonic characteristic dataset X of a user.

[0051] A bus inlet of the industrial load bus is usually equipped with a power quality monitoring device whose sampling interval is generally 3 to 15 min. Monitoring data mainly includes maximum values, minimum values, average values, and the 95% of the maximum value of a fundamental voltage, a fundamental current, active power, reactive power, apparent power, a total harmonic voltage/current distortion rate, a 2 to 25-th harmonic voltage ratio/effective current values, and the like. The present disclosure is intended to use an average value of a 2 to 25-th harmonic current measured at a 3-min sampling interval on Dec. 20, 2021, in a 110 kV steelmaking plant in Taiyuan City, Shanxi Province to perform the analysis. FIG. 2 shows monitoring data of several harmonic currents. Finally, the harmonic characteristic dataset X of the user can be obtained.

[00014]$\begin{array}{cc}X=\left[\begin{array}{cccc}{I}_2^1& I_3^1& L& I_{25}^1\\ I_2^2& I_3^2& L& I_{25}^2\\ M& M& M& M\\ I_2^N& I_3^N& L& I_{25}^N\end{array}\right].& \left(1\right)\end{array}$

[0052] In the above formula, N represents the total quantity of sampling points and is equal to 480 herein. Each column vector in the X represents a harmonic current monitoring sequence, a subscript of I represents a harmonic order, and a superscript of the I represents the quantity of sampling sequences. For example, I.sub.h.sup.m represents an m.sup.th sampling point of an h-th harmonic; j=1, 2, . . . , 480; h represents a harmonic order; and h=2, 3, . . . , 25.

[0053] Step 2: A general probability model is constructed for the h-th harmonic (I.sub.h) in the harmonic characteristic dataset, as shown in the following formula (2):

[00015]$\begin{array}{cc}f\left(I_h\right)=\underset{i=1}{\overset{3}{.Math.}}{\lambda }_i{f}_i\left(I_h\right).& \left(2\right)\end{array}$

[0054] In the above formula, f.sub.i(.) represents a probability density subfunction, and λ.sub.i represents the weight coefficient of the probability density subfunction. The general harmonic probability model is a linear combination of three probability density subfunctions. f.sub.1(.) represents a part that is of the I.sub.h and obeys a normal distribution, f.sub.2(.) represents a part that is of the I.sub.h and obeys a lognormal distribution, and f.sub.3(.) represents a part that is of the I.sub.h and obeys another distribution. Formulas (3) to (5) are mathematical expressions of the f.sub.i(.):

[00016]$\begin{array}{cc}{f}_1\left(I_h\right)=\frac{1}{\sqrt{2\pi }{\sigma }_1}{e}^{-\frac{\left(I_h-{\mu }_1\right)}{2{\sigma }_1^2}},& \left(3\right)\end{array}$$\begin{array}{cc}{f}_2\left(I_h\right)=\frac{1}{\sqrt{2\pi }{I}_h{\sigma }_2}{e}^{-\frac{{\left({\ln I}_h-{\mu }_2\right)}^2}{2{\sigma }_2^2}},& \left(4\right)\end{array}$$\begin{array}{cc}{f}_3\left(I_h\right)=\frac{1}{nb}\underset{{j=1}^{}}{\overset{n}{.Math.}}K\left(\frac{I_h-I_h^j}{b}\right).& \left(5\right)\end{array}$

[0055] In the above formulas, μ.sub.1 and μ.sub.2 represent mathematical expectations of the probability density subfunction, and σ.sub.1 and σ.sub.2 represent standard deviations of the probability density subfunction. K(.) represents a kernel function, and b>0, where b represents a smoothing parameter, which is referred to as a bandwidth or a window. I.sub.h.sup.j represents a j.sup.th sample of the I.sub.h in each window, and n represents the total quantity of samples in each window.

[0056] A probability density function is non-negative and normative, so the weight coefficient of the probability density subfunction should meet the following formula (6).

[0057] In the formula, λ.sub.1=1 or λ.sub.2=1, indicating that the I.sub.h obeys a single normal distribution or the lognormal distribution:

[00017]$\begin{array}{cc}\underset{i=1}{\overset{3}{.Math.}}{\lambda }_i=1\left(0\le {\lambda }_i\le 1\right).& \left(6\right)\end{array}$

[0058] Step 3: The general probability model of the I.sub.h is discretized.

[0059] Since f.sub.1(.) and f.sub.2(.) are continuous functions and f.sub.3(.) is a discrete function, these two types of functions cannot be added directly by using the following formula (7). Instead, it is necessary to discretize f.sub.1(.) and f.sub.2(.). f.sub.1(.) and f.sub.2(.) can be discretized by discretizing the I.sub.h. Therefore, the following discretized general harmonic probability model is obtained:

[00018]$\begin{array}{cc}f\left({\tilde{I}}_h\right)=\underset{i=1}{\overset{3}{.Math.}}{\lambda }_i{f}_i\left({\tilde{I}}_h\right),{\tilde{I}}_h=0:0.01:\max \left({\tilde{I}}_h\right).& \left(7\right)\end{array}$

[0060] In the above formula, max(I.sub.h) represents the maximum current value of the h-th harmonic.

[0061] Step 4: A parameter optimization model of the general probability model is constructed.

[0062] It can be seen from the formulas (3) to (7) that by adjusting a value of a parameter set {λ.sub.1, λ.sub.2, λ.sub.3, μ.sub.1, μ.sub.2, σ.sub.1, σ.sub.2, b}, the above discretized general harmonic probability model can fit and approximate a probability distribution function of any random variable. The smoothing parameter b can be set based on experience, and other parameters need to be solved by constructing the parameter optimization model. The parameter optimization model is mainly divided into the following three parts.

[0063] 1) Determining an Objective Function

[0064] The degree of approximation between the general probability model and the actual probability distribution of the I.sub.h may be visually reflected by the difference between the mathematical expectation calculated by a general optimization model and the actual value, as well as the difference between the standard deviation calculated by the general optimization model and the actual value. Higher accuracy of the parameter optimization model leads to a smaller difference between the mathematical expectations and a smaller difference between the standard deviations. Therefore, the following two objective functions are constructed in this specification:

[00019]$\begin{array}{cc}\min y_1={\left(\underset{i=1}{\overset{3}{.Math.}}{\lambda }_i{E}_i\left({\tilde{I}}_h\right)—E_0\left({\tilde{I}}_h\right)\right)}^2\mathrm{and}& \left(8\right)\end{array}$$\begin{array}{cc}\min y_2=(\sqrt{\underset{i=1}{\overset{3}{.Math.}}{\lambda }_i^2{\mathrm{Var}}_i\left({\tilde{I}}_h\right)}-\sqrt{{{\mathrm{Var}}_0\left({\tilde{I}}_h\right))}^2}.& \left(9\right)\end{array}$

[0065] In the above formulas,

[00020]$\begin{array}{cc}{E}_1\left({\tilde{I}}_h\right)={\mu }_1,{\mathrm{Var}}_1\left({\tilde{I}}_h\right)={\sigma }_1^2,& \left(10\right)\end{array}$$\begin{array}{cc}{E}_2\left({\tilde{I}}_h\right)=e^{{\mu }_2+{\sigma }_2^2/2},{\mathrm{Var}}_2\left({\tilde{I}}_h\right)=\left(e^{{\sigma }_2^2}-1\right)e^{2{\mu }_2+{\sigma }_2^2},& \left(11\right)\end{array}$$\begin{array}{cc}{E}_3\left({\tilde{I}}_h\right)=\underset{j=1}{\overset{N_1}{.Math.}}{\tilde{I}}_h^{j}p\left({\tilde{I}}_h\right),& \left(12\right)\end{array}$$\begin{array}{cc}{\mathrm{Var}}_3\left({\tilde{I}}_h\right)=\underset{j=1}{\overset{N_1}{.Math.}}{\left({\tilde{I}}_h^j-E_3\left({\tilde{I}}_h\right)\right)}^{2}p\left({\tilde{I}}_h^j\right),& \left(13\right)\end{array}$$\begin{array}{cc}{N}_1=\mathrm{length}\left({\tilde{I}}_h\right),\mathrm{and}& \left(14\right)\end{array}$$\begin{array}{cc}p\left({\tilde{I}}_h^j\right)=\frac{f_3\left({\tilde{I}}_h^j\right)}{\mathrm{sum}\left(f_3\left({\tilde{I}}_h\right)\right)}.& \left(15\right)\end{array}$

[0066] A solution to this optimization problem is to combine two minimum objective subfunctions into a minimum objective function, and then an optimization method of a single objective optimization problem is used for solving. The following formula (16) is the combined objective function:

[00021]$\begin{array}{cc}\min y=\frac{y_1+y_2}{2}.& \left(16\right)\end{array}$

[0067] 2) Determining Constraints

[0068] The constraints include an equality constraint and an inequality constraint based on the forms of the constraints.

[0069] Equality constraint: An equality constraint for optimizing the variable λ.sub.i may be determined according to a formula (6) and is expressed by l:

[00022]$\begin{array}{cc}l=\underset{i=1}{\overset{3}{.Math.}}{\lambda }_i-1=0.& \left(17\right)\end{array}$

[0070] Inequality constraint: Based on the optimization parameter set {λ.sub.1, λ.sub.2, λ.sub.3, μ.sub.1, μ.sub.2, σ.sub.1, σ.sub.2, b}, it can be seen that there are mainly two types of inequality constraints. One type of inequality constraint is a value range of the weight coefficient λ.sub.i. The other type of inequality constraint is a value range that is of the random variable I.sub.h and determined by numerical characteristics (μ.sub.1, σ.sub.1) and (μ.sub.2, σ.sub.2) when a single probability density subfunction takes effect.

[0071] An inequality constraint of the λ.sub.i is as follows:

[00023]$\begin{array}{cc}\{\begin{array}{c}{g}_1={\lambda }_1\ge 0\\ g_2={\lambda }_2\ge 0\\ g_3={\lambda }_3\ge 0\end{array}.& \left(18\right)\end{array}$

[0072] Inequality constraints of {μ.sub.1, μ.sub.2, σ.sub.1, σ.sub.2}: This specification uses a maximum likelihood estimation method to evaluate the numerical characteristics of the I.sub.h obeying the single normal distribution or the lognormal distribution and takes the upper and lower confidence limits with a confidence level of 0.95 as value ranges of [μ.sub.1, μ.sub.2, σ.sub.1, σ.sub.2]. Assuming that 95% confidence intervals of {μ.sub.1, —2, σ.sub.1, σ.sub.2} are [θ1, θ2], [θ3, θ4], [θ5, θ6], and [θ7, θ8] respectively, the inequality constraints of the optimization variables {μ.sub.1, μ.sub.2, σ.sub.1, σ.sub.2} can be obtained, namely:

[00024]$\begin{array}{cc}\{\begin{array}{c}{g}_4={\mu }_1-{\hat{\theta }}_1\ge 0\\ g_5={\hat{\theta }}_2-{\mu }_1\ge 0\end{array},& \left(19\right)\end{array}$$\begin{array}{cc}\{\begin{array}{c}{g}_6={\mu }_2-{\hat{\theta }}_3\ge 0\\ g_7={\hat{\theta }}_4-{\mu }_2\ge 0\end{array},& \left(20\right)\end{array}$$\begin{array}{cc}\{\begin{array}{c}{g}_8={\sigma }_1-{\hat{\theta }}_5\ge 0\\ g_9={\hat{\theta }}_6-{\sigma }_1\ge 0\end{array},\mathrm{and}& \left(21\right)\end{array}$$\begin{array}{cc}\{\begin{array}{c}{g}_{10}={\sigma }_2-{\hat{\theta }}_7\ge 0\\ g_{11}={\hat{\theta }}_8-{\sigma }_2\ge 0\end{array}.& \left(22\right)\end{array}$

[0073] The parameter optimization model of the general probability model in the present disclosure is composed of the objective function (16) and the constraint conditions (17) to (22).

[0074] Step 5: Parameters {λ.sub.1, λ.sub.2, λ.sub.3, μ.sub.1, μ.sub.2, σ.sub.1, σ.sub.2} of the general probability model are solved.

[0075] The parameter optimization model in the present disclosure is a constrained optimization problem in an optimization theory. A solution of the parameter optimization model is to convert a constrained problem into an unconstrained problem. Common solving methods include a Lagrange multiplier method, Karush-Kuhn-Tucker (KKT) conditions, a penalty function method, and the like. The parameter optimization model in the present disclosure contains both the equality constraint and the inequality constraint and can be solved directly by using a multiplier method.

[0076] It is assumed that the optimization variable set is expressed by γ, namely, γ={λ.sub.1, λ.sub.2, λ.sub.3, μ.sub.1, μ.sub.2, σ.sub.1, σ.sub.2}, and an augmented Lagrange function is defined as J and is expressed is as follows:

[00025]$\begin{array}{cc}J\left(\gamma ,\omega ,v,\rho \right)=y\left(\gamma \right)-\mathrm{vl}\left(\gamma \right)+\frac{\rho }{2}{l}^2\left(\gamma \right)+\frac{1}{2\rho }\underset{q=1}{\overset{11}{.Math.}}\left\{{\left[\max \left(0,{\omega }_q-\rho g_q\left(\gamma \right)\right)\right]}^2-{\omega }_q^2\right\}.& \left(23\right)\end{array}$

[0077] For the J(γ, ω, ν, ρ), only the sufficiently large parameter ρ needs to be taken, and the multipliers ω and ν are continuously corrected to obtain a local optimal solution in the formula (23) by minimizing the J(γ, ω, ν, ρ), where correction formulas of the multipliers ω and ν are as follows:

[00026]$\begin{array}{cc}\{\begin{array}{c}{\omega }_q^{\left(k+1\right)}=\max \left(0,{\omega }_q^{\left(k\right)}-\rho g_q\left({\gamma }^{\left(k\right)}\right)\right),q=1,2,L,11\\ v^{\left(k+1\right)}=v^{\left(k\right)}-\rho l\left({\gamma }^{\left(k\right)}\right)\end{array}.& \left(24\right)\end{array}$

[0078] FIG. 3 is a basic flowchart of the multiplier method. The basic procedure of the multiplier method is as follows:

[0079] Step 1: Initial point γ.sup.(0), initial multiplier vector estimates ω.sup.(1) and ν.sup.(1), parameter ρ, allowable error ε (>0), constant c (>1), β (ε(0,1)), and k (=1) are set.

[0080] Step 2: γ.sup.(k−1) is taken as an initial point, and an unconstrained problem shown in the following formula (25) is solved to obtain solution γ.sup.(k):

min J(γ,ω.sup.(k),ν.sup.(k),ρ)  (25).

[0081] Step 3: If ∥l(γ.sup.(k))∥<ε, the calculation is stopped, and point γ.sup.(k)) is obtained. Otherwise, step 4 is performed.

[0082] Step 4: If ∥l(γ.sup.(k))∥/∥l(γ.sup.(k−1))∥≥β, ρ=cρ is set, and step 5 is performed. Otherwise, step 5 is performed directly.

[0083] Step 5: Multipliers ω.sub.q.sup.(k+1) and ν.sup.(k+1) are corrected by using the formula (24); k=k+1 is set, and step 2 is performed.

[0084] Step 6: The general probability model of the I.sub.h is obtained.