CN104882884A - System harmonic probability evaluating method based on Markov chain Monte Carlo method - Google Patents

Abstract

The invention discloses a system harmonic probability evaluating method based on a Markov chain Monte Carlo method. The system harmonic probability evaluating method comprises the steps of: s1, establishing an h-time harmonic current model generated by individual harmonic with total number of m in a certain class composite load and a general current model of the h-time harmonic of the certain class composite load injected into a PCC point; s2, analyzing uncertain parameters of the current model according to a Bayes formula in a classical theory of statistics, acquiring more objective and accurate unknown quantity posteriori distribution after adjustment through derivation and combining with rough prior distribution of the uncertain parameters; s3, acquiring the Markov chain Monte Carlo method; and s4, arriving at a probability statistics characteristic value and a probability density curve of percentage of each harmonic current according to the Markov chain Monte Carlo method. The results obtained by adopting the system harmonic probability evaluating method based on the Markov chain Monte Carlo method are more comprehensive, more objective and more similar to the actual condition of a power distribution network.

Description

Based on the system harmonics probability evaluation method of failure of Markov chain-Monte Carlo Method

Technical field

The present invention relates to system for distribution network of power distributed power generation quality of power supply field, particularly relate to a kind of system harmonics probability evaluation method of failure based on Markov chain-Monte Carlo Method.

Background technology

At present, the energy is the important substance basis of human survival and development, and the consumption of mankind's nearly decades to the energy particularly electric energy reaches unprecedented level, and traditional extensive centrally connected power supply pattern can not meet the demand of people gradually.Distributed power generation is as a kind of generation technology based on regenerative resource, and cleanliness without any pollution, and can utilize decentralized resource fully efficiently, has broad application prospects.2011, China's grid-connected generation of electricity by new energy amount 933.55 hundred million kilowatt hour, accounted for 2% of gross generation, and by the end of the year 2012, the grid-connected generation of electricity by new energy installed capacity of China will reach 5,159 ten thousand kilowatts, accounts for 4.89% of total installation of generating capacity.

Along with regenerative resource distributed power source fast development and more and more access electric power system, the power quality problem of distribution system is just more outstanding, containing in the distribution system of distributed power source, distributed power source is all connected with electrical network by electronic power convertor equipment mostly, these equipment all belong to non-linear element, therefore become another group of source point of system harmonics pulse current injectingt.Although the harmonic characterisitic of each independently harmonic source is fixed value in distribution system, but the feature such as stochastic volatility exported due to access load and distributed power source, this just makes to add a lot of non-linear element in power distribution network, increases the power quality problem of harmonic pollution.In order to harmonic levels in accurate and comprehensive assessment power distribution network, how to adopt more science and this power quality index of harmonic injection electric current at effective method accurate evaluation points of common connection (PCC) place, be planning and design of power system and run the major issue controlled.

For stochastic volatility and the time variation of harmonic voltage, electric current in power distribution network, need to adopt the method for probabilistic statistical characteristics to assess percent harmonic distortion.The suggestion of IEEE 519 standard should adopt the gap providing probability curve (comprising block diagram and probability distribution curve) and come between comparative assessment percent harmonic distortion grade and defined limit value.And IEC 61000-3-6 provides the compatibility of the harmonic voltage of intermediate and low electric network and plan degree, and the compatibility rule of thumb compatible grade of demand fulfillment 95%, that is usual satisfied compatible in the time situation of 95%.And in the random process appraisal procedure of harmonic problem, Monte-Carlo simulation is widely used in a large amount of existing documents, tradition Monte Carlo Method also exists the problem being difficult to realize sampling from higher-dimension probability distribution, therefore, existing harmonic current analyzing evaluation method, many employings DATA REASONING or experience are estimated, the mathematical tool lacking more science utilizes, and the result thus produced seems not enough in accuracy, promptness.Especially for a large amount of distribution system using the distributed power source access of new forms of energy, because these power supplys export the random fluctuation characteristic with environmental change had, more result in the seriousness of harmonic pollution, thus the research of problem is become more complicated, required analyzing evaluation method also should science more.Also have some application such as: a kind of method that have employed certainty and randomness mixed problem for nonlinear-load is analyzed, and carries out the method such as predicting according to 95% compatibility of a kind of semiempirical formula to harmonic current.Method involved but, too much relies on historical measurement data and micro-judgment result, and complete objective reaction does not contain the harmonic characterisitic in the power distribution network of distributed power source and multiple nonlinear-load now.

Summary of the invention

The object of this invention is to provide a kind of system harmonics probability evaluation method of failure based on Markov chain-Monte Carlo Method, result can be made more comprehensively objective, closer to actual power distribution network situation.

The present invention adopts following technical proposals: a kind of system harmonics probability evaluation method of failure based on Markov chain-Monte Carlo Method, comprises the following steps:

Step 1: the h subharmonic current model that the individual harmonic being m by total quantity in certain class synthetic load produces is:

$I_h^{\left(AGG\right)}={\mathrm{df}}_h\·\underset{i=1}{\overset{m}{\Σ}}(A_i\·I_{h\left(i\right)})---\left(1\right)$

I _{h (i)}for certain h subharmonic current size that independent harmonic source produces, A _ifor the power number in the harmonic source of single nonlinear-load place, df _hfor kind factor;

The current model that the h subharmonic of injection PCC point class synthetic load is total is:

$I_h^{\left(PCC\right)}=K_E\·I_h^{\left(AGG\right)}=K_E\·{\mathrm{df}}_h\·\underset{i=1}{\overset{m}{\Σ}}(A_i\·I_{h\left(i\right)})---\left(2\right)$

K _efor the participant factor of the power demand of nonlinear-load in synthetic load, usually all below 30%;

Step 2: according to the uncertain parameter A of the Bayes formula in classical theory of statistics to the current model of formula (1) _ianalyze, by conjunction with its comparatively rough prior distribution, by deriving the more objective and accurate unknown quantity Posterior distrbutionp after being adjusted;

Step 3:, ergodic homogeneity according to the time of Markov chain, in conjunction with Monte Carlo Calculation definite integral, builds the random function under two kinds of methods combining, and obtains Markov chain-Monte Carlo Method according to the Posterior distrbutionp reaching Stationary Distribution;

Step 4: according to Markov chain-Monte Carlo Method, adopts Gibbs sampling plan to carry out sampled analog to the described model of the formula set up (2), thus draws probability statistics characteristic value and the probability density curve of individual harmonic current containing ratio.

Described step (2) specifically comprises the following steps:

(201), prior distribution: in the classical theory of statistics, if the information containing unknown parameter θ in sample, then unknown parameter θ is thought a stochastic variable, go to describe by probability or probability distribution, and be just called prior distribution by the probability distribution that the prior information of unknown parameter θ is determined; For prior information extraction and determine, determine the sample space of parameter mainly through subjective probability judgment, expertise and historical data, for discrete variable, ask for its probability according to this sample space; To continuous variable, then construct the priori probability density function π (θ) of unknown parameter θ according to sample space;

(202), Posterior distrbutionp:

Under unknown parameter θ condition, the probability of stochastic variable X is designated as: p (x| θ), and the sample { x|x of X ₁, x ₂..., x _ngeneration need first to extract an observed value θ ' according to the prior distribution of unknown parameter θ, and the conditional probability of stochastic variable X is under this observed value θ ':

$p\left(x\right|{\mathrm{\θ}}^{\′})=\underset{i=1}{\overset{n}{\Π}}p\left(x_i\right|{\mathrm{\θ}}^{\′})---\left(3\right)$

The prior distribution of the sample information that formula (3) produces and unknown parameter θ carries out comprehensively, thus obtains the joint probability density function of x and θ:

f(x,θ)＝p(x|θ)π(θ) (4)

After obtaining the sample of x, need to do θ according to f (x, θ) to infer, therefore, f (x, θ) needs to decompose as follows:

f(x,θ)＝π(θ|x)m(x) (5)

Wherein, π (θ | x) be the conditional probability density function of θ under sample x condition; The marginal probability density function that m (x) is x;

By the joint probability density function of x and θ, marginal probability density function m (x) that can obtain x is as follows:

m(x)＝∫ _Θf(x,θ)dθ＝∫ _Θp(x|θ)π(θ)dθ (6)

Wherein, Θ is the sample space of θ;

Comprehensive above each formula, the probability density function expression formula π of Bayes formula (θ | x) be expressed as:

$\mathrm{\π}\left(\mathrm{\θ}\right|x)=\frac{f(x,\mathrm{\θ})}{m\left(x\right)}=\frac{p\left(x\right|\mathrm{\θ})\mathrm{\π}\left(\mathrm{\θ}\right)}{{\∫}_{\mathrm{\Θ}}p\left(x\right|\mathrm{\θ})\mathrm{\π}\left(\mathrm{\θ}\right)d\mathrm{\θ}}---\left(7\right)$

The distribution of this θ under sample conditions is just called Posterior distrbutionp, also claims Bayes distribution;

(203), the calculating of Posterior distrbutionp:

From formula (7): for including the information of θ in denominator m (x), saved in Practical Calculation, the θ Posterior distrbutionp probability density after being simplified is calculated as:

π(θ|x)∝p(x|θ)π(θ) (8)

Wherein, " ∝ " expression " is proportional to ", in Practical Calculation, with on the right side of above formula form replace Posterior distrbutionp, thus simplify Posterior distrbutionp π (θ | calculating x).

Described step (3) specifically comprises the following steps:

(301), Markov chain:

If one of stochastic variable X group of time dependent sequence x _n, its state space is S, works as t=t _n+1the quantity of state in moment is only by t _nthe quantity of state in moment determines, and and t _nstate before moment has nothing to do, then claim this process to be Markov process, and all discrete Markov process of time and state is called Markov chain; And for a _i, a _j∈ S, n>0, if x is at t=t _nthe quantity of state in moment is a _i, then x is at t=t _n+1moment quantity of state is a _jconditional probability can be expressed as:

p _i,j＝P _ij(1)＝P{x _n+1＝a _j|x _n＝a _i} (9)

Wherein, p _i,jfor a step transition probability of Markov chain, referred to as transition kernel;

When a step transition probability is only relevant with i, j, and when haveing nothing to do with initial time, then claim this transition probability to have stationarity, claim this Markov chain to be temporal homogeneous simultaneously;

For temporal homogeneous Markov chain, if no matter met from any state a _iset out (namely i is arbitrary value), arrive state a _jprobability all level off to a certain value τ _j, then claim this chain to have ergodic, claim τ simultaneously _jthe Limit Distribution of chain for this reason, for Limit Distribution, has following relational expression:

p _j(n)＝P{X _n＝a _j}＝τ _j＝limP _ij(n)j＝1,2,... (10)

Wherein, p _jn () is arbitrary n moment X _n=a _jprobability, formula (10) represents that the distribution p (n) in arbitrary n moment is all consistent with Limit Distribution t, therefore also known as making Stationary Distribution;

(302), Monte Carlo Method

When needs calculation of complex integration time, if g (x) is the probability density function on interval (a, b), then integration type is rewritten as:

${\∫}_a^{b}h\left(x\right)dx={\∫}_a^b\frac{h\left(x\right)}{g\left(x\right)}g\left(x\right)dx=E\[\frac{h\left(x\right)}{g\left(x\right)}\]---\left(11\right)$

Wherein, the E [] desired value that is function;

If parameter a, b are finite value, get being uniformly distributed on g (x)=1/ (b-a), then the unbiased esti-mator of above-mentioned integration type is expressed as:

$\hat{Y}=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}\frac{h\left(x\right)}{g\left(x\right)}=\frac{b-a}{n}\underset{i=1}{\overset{n}{\Σ}}h\left(x_i\right)---\left(12\right)$

Integrating step (203), above-mentioned integration may be used for calculating the Posterior distrbutionp function in Bayes statistics, for calculating integration type Y=∫ p (x| θ) π (θ) the d θ of Posterior distrbutionp function, is approximated by:

$\hat{Y}=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}p\left(x\right|{\mathrm{\θ}}_i)---\left(13\right)$

Wherein, θ _ifor pressing the probability sampling of density function π (θ);

(303), Markov chain-Monte Carlo Method is derived:

According to the step of Bayes statistic law, obtain Posterior distrbutionp probability density function ζ (θ)=π (θ | x) after ∝ p (x| θ) π (θ), need to utilize Posterior distrbutionp to calculate some statistics, certain function h (θ) can be expressed as about the desired value of θ Posterior distrbutionp ζ (θ):

$E\[h\left(\mathrm{\θ}\right)\]={\∫}_{-\mathrm{\∞}}^{x}h\left(\mathrm{\θ}\right)\mathrm{\ζ}\left(\mathrm{\θ}\right)d\mathrm{\θ}---\left(14\right)$

About quantity of state original samples θ ₁if it meets distribution ζ (θ), then by θ arbitrary after the feature of Stationary Distribution _iedge distribution be also ζ (θ), but in reality be often difficult to ensure initial state θ ₁distribution be exactly ζ (θ), need, after n ' step transfer above, just can reach Stationary Distribution ζ (θ), therefore, the individual quantity of state of n-n ' is below used for estimation function desired value, and formula (14) is rewritten into:

$\hat{E}\[h\left(\mathrm{\θ}\right)\]=\frac{1}{n-n^{\′}}\underset{i=n^{\′}+1}{\overset{n}{\Σ}}h\left({\mathrm{\θ}}_i\right)---\left(15\right).$

Described step (4) specifically comprises the following steps:

(401), Gibbs sampling plan:

First parameter θ is divided into B block, piecemeal condition demand fulfillment is from each conditional probability density f (θ _b| θ ₁, θ ₂..., θ _b-1, θ _b+1..., θ _b) in can sample, and then extract required sample from each distribution, detailed process is as follows:

(1) before iteration, setting initial point θ ⁽⁰⁾=(θ ₁ ⁽⁰⁾, θ ₂ ⁽⁰⁾..., θ _b ⁽⁰⁾..., θ _b ⁽⁰⁾);

(2) the 1st iteration, extract as follows from conditional probability density:

θ ₁ ⁽¹⁾～f(θ ₁|θ ₂ ⁽⁰⁾,…,θ _B ⁽⁰⁾)

θ ₂ ⁽¹⁾～f(θ ₂|θ ₁ ⁽¹⁾,θ ₃ ⁽⁰⁾,…,θ _B ⁽⁰⁾)

……

θ _b ⁽¹⁾～f(θ _b|θ ₁ ⁽¹⁾,θ ₂ ⁽¹⁾,…,θ _b-1 ⁽¹⁾,θ _b+1 ⁽⁰⁾,…,θ _B ⁽⁰⁾)

……

θ _B ⁽¹⁾～f(θ _B|θ ₁ ⁽¹⁾,θ ₂ ⁽¹⁾,…,θ _b ⁽¹⁾,…,θ _B-1 ⁽⁰⁾)

(3) i-th iteration (i >=2):

θ ₁ ⁽ⁱ⁾～f(θ ₁|θ ₂ ^(i-1),…,θ _B ^(i-1))

θ ₂ ⁽ⁱ⁾～f(θ ₂|θ ₁ ⁽ⁱ⁾,θ ₃ ^(i-1),…,θ _B ^(i-1))

……

θ _b ⁽ⁱ⁾～f(θ _b|θ ₁ ⁽ⁱ⁾,θ ₂ ⁽ⁱ⁾,…,θ _b-1 ⁽ⁱ⁾,θ _b+1 ^(i-1),…,θ _B ^(i-1))

……

θ _B ⁽ⁱ⁾～f(θ _B|θ ₁ ⁽ⁱ⁾,θ ₂ ⁽ⁱ⁾,…,θ _b ⁽ⁱ⁾,…,θ _B-1 ⁽ⁱ⁾)

Through the iteration of above process, just create sequence θ ⁽⁰⁾, θ ⁽¹⁾..., θ ⁽ⁱ⁾... construct a Markov chain, parameter θ is I _h ^(AGG);

(402) history of harmonic current containing ratio or statistics, is utilized to sample as prior distribution;

(403), I is being obtained _h ^(AGG)prior distribution after, then formula (2) the described model in integrating step 1, samples according to the Gibbs sampling plan described in step (401);

(404), after maximum times j sampling iteration, obtained conditional probability density stably and, as Posterior distrbutionp, secondary for j ' the above quantity of state being considered to not reach Stationary Distribution has been removed, just obtains with I _h ^(AGG)the Stationary Distribution of Markov chain that becomes of the posterior information sample architecture of taking out, and according to probability statistics principle, utilize desired value finally to try to achieve the probability statistics characteristic value of PCC point harmonic current containing ratio, and obtain probability density curve.

The present invention is a kind of method for parameter estimation combined with Monte-Carlo step based on the statistical method of Bayes formulae discovery, analyzing fully on the basis injecting PCC point harmonic current feature containing distributed power source distribution system, proposing a kind of based on Markov chain-Monte Carlo Method first.According to Bayes is theoretical, the prior distribution that historical statistical data is formed is processed, then homogeneity according to Markov chain, finally determine after revising and reach stable Posterior distrbutionp.Finally, carried out case verification at 15 node medium voltage distribution network test macros, result is more comprehensively objective, closer to actual power distribution network situation.

Accompanying drawing explanation

Fig. 1 is the schematic diagram of flow process of the present invention;

Fig. 2 presses radial pattern distribution system line chart in 15 nodes;

Fig. 3 is PCC point harmonic current containing ratio probability density in the 1st period;

Fig. 4 is PCC point harmonic current containing ratio probability density in the 2nd period;

Fig. 5 is PCC point harmonic current containing ratio probability density in the 3rd period.

Embodiment

Below in conjunction with accompanying drawing, preferred embodiment is elaborated.It should be emphasized that following explanation is only exemplary, instead of in order to limit the scope of the invention and apply.

As shown in Figure 1, the invention discloses a kind of system harmonics probability evaluation method of failure based on Markov chain-Monte Carlo Method, specifically comprise the following steps:

Step 1: the h subharmonic current model that the individual harmonic being m by total quantity in certain class synthetic load produces is:

$I_h^{\left(AGG\right)}={\mathrm{df}}_h\·\underset{i=1}{\overset{m}{\Σ}}(A_i\·I_{h\left(i\right)})---\left(1\right)$

I _{h (i)}for certain h subharmonic current size that independent harmonic source produces, A _ifor the power number in the harmonic source of single nonlinear-load place, df _hfor kind factor, it is determined value;

The current model that the h subharmonic of injection PCC point class synthetic load is total is:

$I_h^{\left(PCC\right)}=K_E\·I_h^{\left(AGG\right)}=K_E\·{\mathrm{df}}_h\·\underset{i=1}{\overset{m}{\Σ}}(A_i\·I_{h\left(i\right)})---\left(2\right)$

K _efor the participant factor of the power demand of nonlinear-load in synthetic load, usually all below 30%, can value as required;

Step 2: according to the uncertain parameter A of the Bayes formula in classical theory of statistics to the current model of formula (1) _ianalyze, by conjunction with its comparatively rough prior distribution, by deriving the more objective and accurate unknown quantity Posterior distrbutionp after being adjusted;

Step 3:, ergodic homogeneity according to the time of Markov chain, in conjunction with Monte Carlo Calculation definite integral, builds the random function under two kinds of methods combining, and obtains Markov chain-Monte Carlo Method according to the Posterior distrbutionp reaching Stationary Distribution;

Step 4: according to Markov chain-Monte Carlo Method, adopts Gibbs sampling plan to carry out sampled analog to the described model of the formula set up (2), thus draws probability statistics characteristic value and the probability density curve of individual harmonic current containing ratio.

Described step (2) specifically comprises the following steps:

(201), prior distribution: in the classical theory of statistics, if the information containing unknown parameter θ in sample, then unknown parameter θ is thought a stochastic variable, go to describe by probability or probability distribution, and be just called prior distribution by the probability distribution that the prior information of unknown parameter θ is determined; For prior information extraction and determine, determine the sample space of parameter mainly through subjective probability judgment, expertise and historical data, for discrete variable, ask for its probability according to this sample space; To continuous variable, then construct the priori probability density function π (θ) of unknown parameter θ according to sample space;

(202), Posterior distrbutionp:

Under unknown parameter θ condition, the probability of stochastic variable X is designated as: p (x| θ), and the sample { x|x of X ₁, x ₂..., x _ngeneration need first to extract an observed value θ ' according to the prior distribution of unknown parameter θ, and the conditional probability of stochastic variable X is under this observed value θ ':

$p\left(x\right|{\mathrm{\θ}}^{\′})=\underset{i=1}{\overset{n}{\Π}}p\left(x_i\right|{\mathrm{\θ}}^{\′})---\left(3\right)$

The prior distribution of the sample information that formula (3) produces and unknown parameter θ carries out comprehensively, thus obtains the joint probability density function of x and θ:

f(x,θ)＝p(x|θ)π(θ) (4)

After obtaining the sample of x, need to do θ according to f (x, θ) to infer, therefore, f (x, θ) needs to decompose as follows:

f(x,θ)＝π(θ|x)m(x) (5)

Wherein, π (θ | x) be the conditional probability density function of θ under sample x condition; The marginal probability density function that m (x) is x;

By the joint probability density function of x and θ, marginal probability density function m (x) that can obtain x is as follows:

m(x)＝∫ _Θf(x,θ)dθ＝∫ _Θp(x|θ)π(θ)dθ (6)

Wherein, Θ is the sample space of θ;

Comprehensive above each formula, the probability density function expression formula π of Bayes formula (θ | x) be expressed as:

$\mathrm{\π}\left(\mathrm{\θ}\right|x)=\frac{f(x,\mathrm{\θ})}{m\left(x\right)}=\frac{p\left(x\right|\mathrm{\θ})\mathrm{\π}\left(\mathrm{\θ}\right)}{{\∫}_{\mathrm{\Θ}}p\left(x\right|\mathrm{\θ})\mathrm{\π}\left(\mathrm{\θ}\right)d\mathrm{\θ}}---\left(7\right)$

The distribution of this θ under sample conditions is just called Posterior distrbutionp, also claims Bayes distribution;

(203), the calculating of Posterior distrbutionp:

From formula (7): for including the information of θ in denominator m (x), saved in Practical Calculation, the θ Posterior distrbutionp probability density after being simplified is calculated as:

π(θ|x)∝p(x|θ)π(θ) (8)

Wherein, " ∝ " expression " is proportional to ", in Practical Calculation, with on the right side of above formula form replace Posterior distrbutionp, thus simplify Posterior distrbutionp π (θ | calculating x).

Described step (3) specifically comprises the following steps:

(301), Markov chain:

If one of stochastic variable X group of time dependent sequence x _n, its state space is S, works as t=t _n+1the quantity of state in moment is only by t _nthe quantity of state in moment determines, and and t _nstate before moment has nothing to do, then claim this process to be Markov process, and all discrete Markov process of time and state is called Markov chain; And for a _i, a _j∈ S, n>0, if x is at t=t _nthe quantity of state in moment is a _i, then x is at t=t _n+1moment quantity of state is a _jconditional probability can be expressed as:

p _i,j＝P _ij(1)＝P{x _n+1＝a _j|x _n＝a _i} (9)

Wherein, p _i,jfor a step transition probability of Markov chain, referred to as transition kernel;

When a step transition probability is only relevant with i, j, and when haveing nothing to do with initial time, then claim this transition probability to have stationarity, claim this Markov chain to be temporal homogeneous simultaneously;

For temporal homogeneous Markov chain, if no matter met from any state a _iset out (namely i is arbitrary value), arrive state a _jprobability all level off to a certain value τ _j, then claim this chain to have ergodic, claim τ simultaneously _jthe Limit Distribution of chain for this reason, for Limit Distribution, has following relational expression:

p _j(n)＝P{X _n＝a _j}＝τ _j＝limP _ij(n)j＝1,2,... (10)

Wherein, p _jn () is arbitrary n moment X _n=a _jprobability, formula (10) represents that the distribution p (n) in arbitrary n moment is all consistent with Limit Distribution t, therefore also known as making Stationary Distribution;

(302), Monte Carlo Method

When needs calculation of complex integration time, if g (x) is the probability density function on interval (a, b), then integration type is rewritten as:

${\∫}_a^{b}h\left(x\right)dx={\∫}_a^b\frac{h\left(x\right)}{g\left(x\right)}g\left(x\right)dx=E\[\frac{h\left(x\right)}{g\left(x\right)}\]---\left(11\right)$

Wherein, the E [] desired value that is function;

If parameter a, b are finite value, get being uniformly distributed on g (x)=1/ (b-a), then the unbiased esti-mator of above-mentioned integration type is expressed as:

$\hat{Y}=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}\frac{h\left(x\right)}{g\left(x\right)}=\frac{b-a}{n}\underset{i=1}{\overset{n}{\Σ}}h\left(x_i\right)---\left(12\right)$

Integrating step (203), above-mentioned integration may be used for calculating the Posterior distrbutionp function in Bayes statistics, for calculating integration type Y=∫ p (x| θ) π (θ) the d θ of Posterior distrbutionp function, is approximated by:

$\hat{Y}=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}p\left(x\right|{\mathrm{\θ}}_i)---\left(13\right)$

Wherein, θ _ifor pressing the probability sampling of density function π (θ);

(303), Markov chain-Monte Carlo Method is derived:

According to the step of Bayes statistic law, obtain Posterior distrbutionp probability density function ζ (θ)=π (θ | x) after ∝ p (x| θ) π (θ), need to utilize Posterior distrbutionp to calculate some statistics, certain function h (θ) can be expressed as about the desired value of θ Posterior distrbutionp ζ (θ):

$E\[h\left(\mathrm{\θ}\right)\]={\∫}_{-\mathrm{\∞}}^{x}h\left(\mathrm{\θ}\right)\mathrm{\ζ}\left(\mathrm{\θ}\right)d\mathrm{\θ}---\left(14\right)$

About quantity of state original samples θ ₁if it meets distribution ζ (θ), then by θ arbitrary after the feature of Stationary Distribution _iedge distribution be also ζ (θ), but in reality be often difficult to ensure initial state θ ₁distribution be exactly ζ (θ), need, after n ' step transfer above, just can reach Stationary Distribution ζ (θ), therefore, the individual quantity of state of n-n ' is below used for estimation function desired value, and formula (14) is rewritten into:

$\hat{E}\[h\left(\mathrm{\θ}\right)\]=\frac{1}{n-n^{\′}}\underset{i=n^{\′}+1}{\overset{n}{\Σ}}h\left({\mathrm{\θ}}_i\right)---\left(15\right).$

Described step (4) specifically comprises the following steps:

(401), Gibbs sampling plan:

First parameter θ is divided into B block, piecemeal condition demand fulfillment is from each conditional probability density f (θ _b| θ ₁, θ ₂..., θ _b-1, θ _b+1..., θ _b) in can sample, and then extract required sample from each distribution, detailed process is as follows:

(1) before iteration, setting initial point θ ⁽⁰⁾=(θ ₁ ⁽⁰⁾, θ ₂ ⁽⁰⁾..., θ _b ⁽⁰⁾..., θ _b ⁽⁰⁾);

(2) the 1st iteration, extract as follows from conditional probability density:

θ ₁ ⁽¹⁾～f(θ ₁|θ ₂ ⁽⁰⁾,…,θ _B ⁽⁰⁾)

θ ₂ ⁽¹⁾～f(θ ₂|θ ₁ ⁽¹⁾,θ ₃ ⁽⁰⁾,…,θ _B ⁽⁰⁾)

……

θ _b ⁽¹⁾～f(θ _b|θ ₁ ⁽¹⁾,θ ₂ ⁽¹⁾,…,θ _b-1 ⁽¹⁾,θ _b+1 ⁽⁰⁾,…,θ _B ⁽⁰⁾)

……

θ _B ⁽¹⁾～f(θ _B|θ ₁ ⁽¹⁾,θ ₂ ⁽¹⁾,…,θ _b ⁽¹⁾,…,θ _B-1 ⁽⁰⁾)

(3) i-th iteration (i >=2):

θ ₁ ⁽ⁱ⁾～f(θ ₁|θ ₂ ^(i-1),…,θ _B ^(i-1))

θ ₂ ⁽ⁱ⁾～f(θ ₂|θ ₁ ⁽ⁱ⁾,θ ₃ ^(i-1),…,θ _B ^(i-1))

……

θ _b ⁽ⁱ⁾～f(θ _b|θ ₁ ⁽ⁱ⁾,θ ₂ ⁽ⁱ⁾,…,θ _b-1 ⁽ⁱ⁾,θ _b+1 ^(i-1),…,θ _B ^(i-1))

……

θ _B ⁽ⁱ⁾～f(θ _B|θ ₁ ⁽ⁱ⁾,θ ₂ ⁽ⁱ⁾,…,θ _b ⁽ⁱ⁾,…,θ _B-1 ⁽ⁱ⁾)

Through the iteration of above process, just create sequence θ ⁽⁰⁾, θ ⁽¹⁾..., θ ⁽ⁱ⁾... construct a Markov chain, parameter θ is I _h ^(AGG);

(402) history of harmonic current containing ratio or statistics, is utilized to sample as prior distribution;

(403), I is being obtained _h ^(AGG)prior distribution after, then formula (2) the described model in integrating step 1, samples according to the Gibbs sampling plan described in step (401);

(404), after maximum times j sampling iteration, obtained conditional probability density stably and, as Posterior distrbutionp, secondary for j ' the above quantity of state being considered to not reach Stationary Distribution has been removed, just obtains with I _h ^(AGG)the Stationary Distribution of Markov chain that becomes of the posterior information sample architecture of taking out, and according to probability statistics principle, utilize desired value finally to try to achieve the probability statistics characteristic value (average, standard deviation) of PCC point harmonic current containing ratio, and obtain probability density curve.

The present invention adopts in certain 15 node and presses radial pattern distribution system data, and its topological structure as shown in Figure 2.

This distribution system PCC point voltage is 11kV, and three-phase shortcircuit capacity is 300MVA.4 loads in Fig. 2 represent the comprehensive harmonic-producing load of 4 class, and use A respectively, B, C, D represent, all contain independent nonlinear-load and the distributed power source of some quantity in every class synthetic load.I _h ^(AGG)prior distribution be considered to normal distribution, and the different harmonic current I that every class synthetic load to produce respectively from underload in high load capacity period period _h ^(AGG)with fundamental current value I ₁ratio and harmonic current containing ratio, its dependent probability statistical nature data (average and standard deviation) are in table 1.First according to I _h ^(AGG)prior distribution sample, about I _h ^(AGG)the acquisition of prior distribution, adopt the data in table 1, this packet contains the probability statistics feature of 5 times and 7 times two kinds of comparatively typical harmonic currents in power distribution network, is by obtaining the measurement in harmonic current magnetic field.

5 times and 7 subharmonic current average and the variances of table 1 high load capacity phase and low-load period

According to actual conditions needs, the change of load in a day is divided into 3 periods, and in each period, 4 class synthetic loads are in different high load capacities and low load condition.Concrete condition is shown in Table 2.

Every class synthetic load state in 3 periods of table 2

Carry out emulation sampling assessment by MATLAB2008a coding, total frequency in sampling is 30000 times, and give up 4000 times above, sample reaches Stationary Distribution.Finally can obtain the probability level of certain percent harmonic distortion of PCC point in each period, as shown in table 3; Dependent probability distribution is shown in shown in Fig. 3-Fig. 5, wherein Fig. 3 is PCC point harmonic current containing ratio probability density 5 subharmonic current containing ratio and 7 subharmonic current containing ratios in the 1st period, wherein Fig. 4 is PCC point harmonic current containing ratio probability density 5 subharmonic current containing ratio and 7 subharmonic current containing ratios in the 2nd period, and wherein Fig. 5 is PCC point harmonic current containing ratio probability density 5 subharmonic current containing ratio and 7 subharmonic current containing ratios in the 3rd period.

Table 4-3 PCC point harmonic current containing ratio probability level

Invention introduces Markov Chain Monte Carlo method (MCMC) method, the flow process solving the Mathematical Modeling containing the Distribution Network Harmonics electric current evaluation problem using the access of the distributed power source of new forms of energy and deal with problems.In conjunction with concrete example, adopt in actual power distribution network and measure the harmonic current statistics of the dissimilar synthetic load obtained as prior distribution by Harmonic Field, utilize in Bayes theory and there is with markovian homogeneity structure one the harmonic current Posterior distrbutionp after the adjustment of Stationary Distribution carry out sampling and emulate, finally obtain statistical property data and the probability distribution of each harmonic Injection Current of power distribution network PCC point, this result is more clear compared with the normal distribution curve formed by prior distribution, also demonstrate and utilize this algorithm to carry out the validity revised.Have employed comparatively typical 5 times and 7 subharmonic data, and the impact that the ratio considering linear load and nonlinear-load in Different periods is brought, make result more comprehensively objective, closer to actual power distribution network situation.

Claims (4)

1., based on a system harmonics probability evaluation method of failure for Markov chain-Monte Carlo Method, it is characterized in that: comprise the following steps:

Step 1: the h subharmonic current model that the individual harmonic being m by total quantity in certain class synthetic load produces is:

$I_h^{\left(\mathrm{AGG}\right)}={\mathrm{df}}_h\·\underset{i=1}{\overset{m}{\mathrm{\Σ}}}(A_i\·I_{h\left(i\right)})---\left(1\right)$

I _{h (i)}for certain h subharmonic current size that independent harmonic source produces, A _ifor the power number in the harmonic source of single nonlinear-load place, df _hfor kind factor;

The current model that the h subharmonic of injection PCC point class synthetic load is total is:

$I_h^{\left(\mathrm{PCC}\right)}=K_E\·I_h^{\left(\mathrm{AGG}\right)}=K_E\·{\mathrm{df}}_h\·\underset{i=1}{\overset{m}{\mathrm{\Σ}}}(A_i\·I_{h\left(i\right)})---\left(2\right)$

K _efor the participant factor of the power demand of nonlinear-load in synthetic load, usually all below 30%;

Step 2: according to the uncertain parameter A of the Bayes formula in classical theory of statistics to the current model of formula (1) _ianalyze, by conjunction with its comparatively rough prior distribution, by deriving the more objective and accurate unknown quantity Posterior distrbutionp after being adjusted;

Step 3:, ergodic homogeneity according to the time of Markov chain, in conjunction with Monte Carlo Calculation definite integral, builds the random function under two kinds of methods combining, and obtains Markov chain-Monte Carlo Method according to the Posterior distrbutionp reaching Stationary Distribution;

Step 4: according to Markov chain-Monte Carlo Method, adopts Gibbs sampling plan to carry out sampled analog to the described model of the formula set up (2), thus draws probability statistics characteristic value and the probability density curve of individual harmonic current containing ratio.

2. a kind of system harmonics probability evaluation method of failure based on Markov chain-Monte Carlo Method according to claim 1, is characterized in that: described step (2) specifically comprises the following steps:

(201), prior distribution: in the classical theory of statistics, if the information containing unknown parameter θ in sample, then unknown parameter θ is thought a stochastic variable, go to describe by probability or probability distribution, and be just called prior distribution by the probability distribution that the prior information of unknown parameter θ is determined; For prior information extraction and determine, determine the sample space of parameter mainly through subjective probability judgment, expertise and historical data, for discrete variable, ask for its probability according to this sample space; To continuous variable, then construct the priori probability density function π (θ) of unknown parameter θ according to sample space;

(202), Posterior distrbutionp:

Under unknown parameter θ condition, the probability of stochastic variable X is designated as: p (x| θ), and the sample { x|x of X ₁, x ₂..., x _ngeneration need first to extract an observed value θ ' according to the prior distribution of unknown parameter θ, and the conditional probability of stochastic variable X is under this observed value θ ':

$p\left(x\right|{\mathrm{\θ}}^{\′})=\underset{i=1}{\overset{n}{\mathrm{\Π}}}p\left(x_i\right|{\mathrm{\θ}}^{\′})---\left(3\right)$

The prior distribution of the sample information that formula (3) produces and unknown parameter θ carries out comprehensively, thus obtains the joint probability density function of x and θ:

f(x,θ)＝p(x|θ)π(θ) (4)

After obtaining the sample of x, need to do θ according to f (x, θ) to infer, therefore, f (x, θ) needs to decompose as follows:

f(x,θ)＝π(θ|x)m(x) (5)

Wherein, π (θ | x) be the conditional probability density function of θ under sample x condition; The marginal probability density function that m (x) is x;

By the joint probability density function of x and θ, marginal probability density function m (x) that can obtain x is as follows:

m(x)＝∫ _Θf(x,θ)dθ＝∫ _Θp(x|θ)π(θ)dθ (6)

Wherein, Θ is the sample space of θ;

Comprehensive above each formula, the probability density function expression formula π of Bayes formula (θ | x) be expressed as:

$\mathrm{\π}\left(\mathrm{\θ}\right|x)=\frac{f(x,\mathrm{\θ})}{m\left(x\right)}=\frac{p\left(x\right|\mathrm{\θ})\mathrm{\π}\left(\mathrm{\θ}\right)}{{\∫}_{\mathrm{\Θ}}\left(x\right|\mathrm{\θ})\mathrm{\π}\left(\mathrm{\θ}\right)\mathrm{d\θ}}---\left(7\right)$

The distribution of this θ under sample conditions is just called Posterior distrbutionp, also claims Bayes distribution;

(203), the calculating of Posterior distrbutionp:

From formula (7): for including the information of θ in denominator m (x), saved in Practical Calculation, the θ Posterior distrbutionp probability density after being simplified is calculated as:

π(θ|x)∝p(x|θ)π(θ) (8)

Wherein, " ∝ " expression " is proportional to ", in Practical Calculation, with on the right side of above formula form replace Posterior distrbutionp, thus simplify Posterior distrbutionp π (θ | calculating x).

3. a kind of system harmonics probability evaluation method of failure based on Markov chain-Monte Carlo Method according to claim 2, is characterized in that: described step (3) specifically comprises the following steps:

(301), Markov chain:

If one of stochastic variable X group of time dependent sequence x _n, its state space is S, works as t=t _n+1the quantity of state in moment is only by t _nthe quantity of state in moment determines, and and t _nstate before moment has nothing to do, then claim this process to be Markov process, and all discrete Markov process of time and state is called Markov chain; And for a _i, a _j∈ S, n>0, if x is at t=t _nthe quantity of state in moment is a _i, then x is at t=t _n+1moment quantity of state is a _jconditional probability can be expressed as:

p _i,j＝P _ij(1)＝P{x _n+1＝a _j|x _n＝a _i} (9)

Wherein, p _i,jfor a step transition probability of Markov chain, referred to as transition kernel;

When a step transition probability is only relevant with i, j, and when haveing nothing to do with initial time, then claim this transition probability to have stationarity, claim this Markov chain to be temporal homogeneous simultaneously;

For temporal homogeneous Markov chain, if no matter met from any state a _iset out (namely i is arbitrary value), arrive state a _jprobability all level off to a certain value τ _j, then claim this chain to have ergodic, claim τ simultaneously _jthe Limit Distribution of chain for this reason, for Limit Distribution, has following relational expression:

p _j(n)＝P{X _n＝a _j}＝τ _j＝limP _ij(n)j＝1,2,... (10)

Wherein, p _jn () is arbitrary n moment X _n=a _jprobability, formula (10) represents that the distribution p (n) in arbitrary n moment is all consistent with Limit Distribution t, therefore also known as making Stationary Distribution;

(302), Monte Carlo Method

When needs calculation of complex integration time, if g (x) is the probability density function on interval (a, b), then integration type is rewritten as:

${\∫}_a^{b}h\left(x\right)\mathrm{dx}={\∫}_a^b\frac{h\left(x\right)}{g\left(x\right)}g\left(x\right)\mathrm{dx}=E\left[\frac{h\left(x\right)}{g\left(x\right)}\right]---\left(11\right)$

Wherein, the E [] desired value that is function;

If parameter a, b are finite value, get being uniformly distributed on g (x)=1/ (b-a), then the unbiased esti-mator of above-mentioned integration type is expressed as:

$\hat{Y}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{h\left(x\right)}{g\left(x\right)}=\frac{b-a}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}h\left(x_i\right)---\left(12\right)$

Integrating step (203), above-mentioned integration may be used for calculating the Posterior distrbutionp function in Bayes statistics, for calculating integration type Y=∫ p (x| θ) π (θ) the d θ of Posterior distrbutionp function, is approximated by:

$\hat{Y}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}p\left(x\right|{\mathrm{\θ}}_i)---\left(13\right)$

Wherein, θ _ifor pressing the probability sampling of density function π (θ);

(303), Markov chain-Monte Carlo Method is derived:

According to the step of Bayes statistic law, obtain Posterior distrbutionp probability density function ζ (θ)=π (θ | x) after ∝ p (x| θ) π (θ), need to utilize Posterior distrbutionp to calculate some statistics, certain function h (θ) can be expressed as about the desired value of θ Posterior distrbutionp ζ (θ):

$E\left[h\left(\mathrm{\θ}\right)\right]={\∫}_{-\∞}^{x}h\left(\mathrm{\θ}\right)\mathrm{\ζ}\left(\mathrm{\θ}\right)\mathrm{d\θ}---\left(14\right)$

About quantity of state original samples θ ₁if it meets distribution ζ (θ), then by θ arbitrary after the feature of Stationary Distribution _iedge distribution be also ζ (θ), but in reality be often difficult to ensure initial state θ ₁distribution be exactly ζ (θ), need, after n ' step transfer above, just can reach Stationary Distribution ζ (θ), therefore, the individual quantity of state of n-n ' is below used for estimation function desired value, and formula (14) is rewritten into:

$\hat{E}\left[h\left(\mathrm{\θ}\right)\right]=\frac{1}{n-n^{\′}}\underset{i=n^{\′}+1}{\overset{n}{\mathrm{\Σ}}}h\left({\mathrm{\θ}}_i\right)---\left(15\right).$

4. a kind of system harmonics probability evaluation method of failure based on Markov chain-Monte Carlo Method according to claim 3, is characterized in that: described step (4) specifically comprises the following steps:

(401), Gibbs sampling plan:

First parameter θ is divided into B block, piecemeal condition demand fulfillment is from each conditional probability density f (θ _b| θ ₁, θ ₂..., θ _b-1, θ _b+1..., θ _b) in can sample, and then extract required sample from each distribution, detailed process is as follows:

(1) before iteration, setting initial point θ ⁽⁰⁾=(θ ₁ ⁽⁰⁾, θ ₂ ⁽⁰⁾..., θ _b ⁽⁰⁾..., θ _b ⁽⁰⁾);

(2) the 1st iteration, extract as follows from conditional probability density:

θ ₁ ⁽¹⁾～f(θ ₁|θ ₂ ⁽⁰⁾,…,θ _B ⁽⁰⁾)

θ ₂ ⁽¹⁾～f(θ ₂|θ ₁ ⁽¹⁾,θ ₃ ⁽⁰⁾,…,θ _B ⁽⁰⁾)

……

θ _b ⁽¹⁾～f(θ _b|θ ₁ ⁽¹⁾,θ ₂ ⁽¹⁾,…,θ _b-1 ⁽¹⁾,θ _b+1 ⁽⁰⁾,…,θ _B ⁽⁰⁾)

……

θ _B ⁽¹⁾～f(θ _B|θ ₁ ⁽¹⁾,θ ₂ ⁽¹⁾,…,θ _b ⁽¹⁾,…,θ _B-1 ⁽⁰⁾)

(3) i-th iteration (i >=2):

θ ₁ ⁽ⁱ⁾～f(θ ₁|θ ₂ ^(i-1),…,θ _B ^(i-1))

θ ₂ ⁽ⁱ⁾～f(θ ₂|θ ₁ ⁽ⁱ⁾,θ ₃ ^(i-1),…,θ _B ^(i-1))

……

θ _b ⁽ⁱ⁾～f(θ _b|θ ₁ ⁽ⁱ⁾,θ ₂ ⁽ⁱ⁾,…,θ _b-1 ⁽ⁱ⁾,θ _b+1 ^(i-1),…,θ _B ^(i-1))

……

θ _B ⁽ⁱ⁾～f(θ _B|θ ₁ ⁽ⁱ⁾,θ ₂ ⁽ⁱ⁾,…,θ _b ⁽ⁱ⁾,…,θ _B-1 ⁽ⁱ⁾)

Through the iteration of above process, just create sequence θ ⁽⁰⁾, θ ⁽¹⁾..., θ ⁽ⁱ⁾... construct a Markov chain, parameter θ is I _h ^(AGG);

(402) history of harmonic current containing ratio or statistics, is utilized to sample as prior distribution;

(403), I is being obtained _h ^(AGG)prior distribution after, then formula (2) the described model in integrating step 1, samples according to the Gibbs sampling plan described in step (401);

(404), after maximum times j sampling iteration, obtained conditional probability density stably and, as Posterior distrbutionp, secondary for j ' the above quantity of state being considered to not reach Stationary Distribution has been removed, just obtains with I _h ^(AGG)the Stationary Distribution of Markov chain that becomes of the posterior information sample architecture of taking out, and according to probability statistics principle, utilize desired value finally to try to achieve the probability statistics characteristic value of PCC point harmonic current containing ratio, and obtain probability density curve.