CN103336427B - A kind of dynamic modeling of three-dimensional output probability density function and controller design method - Google Patents

Abstract

The invention discloses dynamic modeling and the controller design method of a kind of three-dimensional output probability density function in the theoretical field of stochastic distribution control.The method comprises the following steps: step 1: the three-dimensional built based on square root B-spline model exports PDF dynamic model; Step 2: utilize the inputoutput data collected in real system to set up the three-dimensional input/output model exporting PDF by recursive least squares; Step 3: select instantaneous square root performance index CONTROLLER DESIGN, by the controlled quentity controlled variable of optimization instantaneous square root performance index CONTROLLER DESIGN, the system that realizes exports the shape that PDF distribution shape tracing preset exports PDF distribution; The present invention devises optimization routine controller, by optimizing square root performance index, realizes exporting the tracking of PDF distribution to given output PDF.The present invention has enriched three-dimensional output PDF control theory, provides new method for having the three-dimensional industrial process exporting distribution character.

Description

A kind of dynamic modeling of three-dimensional output probability density function and controller design method

Technical field

The invention belongs to the theoretical field of stochastic distribution control, particularly relate to a kind of dynamic modeling and controller design method of three-dimensional output probability density function.

Background technology

Stochastic system control theory is one of important branch of control theory and application, mainly because most industrial process is all subject to the interference of random signal, for this practical problems, has formed the stochastic control theory of system.Its early stage achievement in research concentrates on the statistical property of control system variable self, control objectives is average and the variance of system, in existing method, stochastic variable Gaussian distributed in major part hypothesis stochastic system, but the not realistic application of this hypothesis, the distribution of such as, fiber length distribution in paper-making process, the grain in grain processing and boiler flame temperature distribution etc.General random system distribution output PDF (probability density function) probability density function represents, when stochastic variable is gaussian variable, can realize controlling the output PDF of system by the average of control system and variance, stochastic variable is not met to the system of Gaussian distribution, its average and variance can not comprise the full detail of system, can not realize to the control of system average and variance the control this system being exported to PDF.For this type systematic, Wang Hong teaches and proposed in 1998 the method that direct control system exports PDF shape, namely exports PDF and controls.The direct CONTROLLER DESIGN of these class methods exports PDF distribution shape tracing preset PDF distribution shape to make system.

Export PDF to be approached by B-spline neural network, achieve the decoupling zero of original complicated coupling system so to a certain extent, output PDF for this kind of decoupling zero controls, and is referred to as SDC (stochasticdistribution control) stochastic system distributed controll.These class methods breach the limitation of STOCHASTIC CONTROL research, are converted into traditional Stochastic Control Problem and set up advantages of simple model, design the method for efficient suitable control algolithm.The state-space model that the partial differential equation realizing implying descriptive system dynamic perfromance from complexity is converted into decoupling zero describes, and finally reaches to be similar to that certainty annuity is the same can describe the dynamic behaviour of stochastic distribution system with model more accurately.More tally with the actual situation compared with SDC and stochastic control theory in the past, therefore, this theory is dissolved into new industrial circle, will its powerful vitality be given.

For the stochastic distribution system with two-dimensional characteristics, establish comparatively perfect theoretical system.As in system modelling, establish linear B-spline model, Rational B-splines model, square root B-spline model and reasonable square root B-spline model, input and output ARMAX model, neural network PDF model etc.In Controller gain variations, achieve instantaneous optimal track control algorithm, optimal track control algorithm, model reference self-adapting control algorithm, predictive control algorithm, structuring controller algorithm, Iterative Learning Control Algorithm etc.In recent years, scholars have done a large amount of work in the robust control of stochastic distribution control, minimum entropy control, fault diagnosis and design of filter etc.

In sum, export distributed controll problem for two dimension to make great progress, also there is a class three-dimensional in actual industrial process and export distribution problem, as characterized the distributed in three dimensions etc. of material concentration in the three-dimensional temperature field of boiler flame temperature, the fire coal circulating fluid bed boiler in power station.These distributed in three dimensions and industrial process operation conditions closely bound up, this enhances productivity to whole industrial process, reduce the aspects such as environmental pollution all has important using value.Along with the high speed development of Computer Image Processing and sensor technology, the detection exporting distribution situation to three-dimensional obtains and develops rapidly.But the method being obtained system output distribution by advanced technology on-line measurement implements very complicated, required apparatus expensive.

The three-dimensional of the present invention's research exports the important component part that PDF control problem is SDC theory, but imperfection is gone back to the research of the modeling and control problem of three-dimensional PDF, article " Modeling and control ofthe flame temperature distribution using probability density function shaping " exports PDF to three-dimensional and has carried out static modelling and Controller gain variations, select two-dimentional B-spline basis function and establish the three-dimensional static model exporting PDF by least-squares algorithm, optimize quadratic performance index, system local optimum control inputs is obtained by the method for gradient, computer simulation obtains legitimate result.So far, the research for dynamic process three-dimensional output distributed controll aspect rarely has to be delivered, but not direct using the research report of three-dimensional output distribution as control.

In order to improve three-dimensional output PDF control theory further, three-dimensional random distributed controll problem is made to realize becoming possibility, first the present invention establishes the three-dimensional instantaneous square root B-spline model exporting PDF, the three-dimensional that the dynamic change part that instantaneous square root B-spline model basis adds weights constitutes based on square root B-spline model exports PDF dynamic model, achieve the dynamic decoupling between weights, analyze the condition that three-dimensional output PDF dynamic model meets natural sulfur reservoir; Then three-dimensional output PDF input/output model is established according to system inputoutput data by recursive least squares; Finally select instantaneous square root performance index, devise conventional optimal controller.The present invention exports PDF control theory to three-dimensional and carries out perfect, and the control exporting distribution problem for three-dimensional provides new method and thinking.

Summary of the invention

The present invention is directed to the needs that three-dimensional exports being left to be desired of PDF theory and actual industrial process, propose a kind of dynamic modeling and controller design method of three-dimensional output probability density function.

The dynamic modeling of three-dimensional output probability density function and a controller design method, the method comprises the following steps:

Step 1: the three-dimensional built based on square root B-spline model exports PDF dynamic model;

Described structure comprises the following steps based on the three-dimensional output PDF dynamic model of square root B-spline model:

Step S1: build the three-dimensional instantaneous square root B-spline model exporting PDF according to two-dimentional B-spline function;

As follows with two one dimension B-spline function tensor product representation two dimension B-spline functions:

$B_{j,i}(x,r)=B_{j_x{i}_x}\left(x\right)B_{j_r{i}_r}\left(r\right)$

Wherein, obtained by following recursion formula:

$B_{1,i_x}\left(x\right)=\left\{\begin{array}{cc}1,& x\∈[x_{i_x+1},x_{i_x+1})\\ 0,& x\∉[x_{i_x},x_{i_x+1})\end{array}\right.$

$B_{j_x,i_x}\left(x\right)=\frac{x-x_{i_x}}{x_{i_x+j_x-1}-x_{i_x}}{B}_{j_x-1,i_x}\left(x\right)+\frac{x_{i_x+j_x}-x}{x_{i_x+j_x}-x_{i_x+1}}{B}_{j_x-1,i_x+1}\left(x\right)$

obtained by following recursion formula:

$B_{1,i_r}\left(r\right)=\left\{\begin{array}{cc}1,& r\∈[r_{i_r+1},r_{i_r+1})\\ 0,& r\∉[r_{i_r},r_{i_r+1})\end{array}\right.$

$B_{j_r,i_r}\left(r\right)=\frac{r-r_{i_r}}{r_{i_r+j_r-1}-r_{i_r}}{B}_{j_r-1,i_r}\left(r\right)+\frac{r_{i_r+j_r}-r}{r_{i_r+j_r}-r_{i_r+1}}{B}_{j_r-1,i_r+1}\left(r\right)$

Wherein, B _j,i(x, r) is two-dimentional B-spline basis function; for one dimension B-spline basis function; for one dimension B-spline basis function; X, r are respectively the variable spatially defined, x ∈ [a ₁, b ₁], r ∈ [a ₂, b ₂]; a ₁for the lower limit of setting in X-axis setting district; b ₁for the higher limit of setting in X-axis setting district; [a ₁, b ₁] for comprising a ₁and b ₁an interval; a ₂for the lower limit of setting in R axle setting district; b ₂for the higher limit of setting in R axle setting district; [a ₂, b ₂] for comprising a ₂and b ₂an interval; J represents the order of two-dimentional B-spline; I represents two-dimentional B-spline basis function number; j _xfor the order of basis function that X-axis is chosen; i _xfor the number of basis function that X-axis is chosen; j _rfor the order of basis function that R axle is chosen; i _rfor the number of basis function that R axle is chosen;

be 1 rank i-th _xindividual B-spline function; for j _x-1 rank i-th _xindividual B-spline function; for j _x-1 rank i-th _x+ 1 B-spline basis function; for nodal value, and have m _xfor interval [a ₁, b ₁] in effective nodes, j _x-1 counts for the acromere of the interval left and right sides; for comprising but do not comprise an interval;

be 1 rank i-th _rindividual B-spline function; for j _r-1 rank i-th _rindividual B-spline function; for j _r-1 rank i-th _r+ 1 B-spline basis function; for nodal value, and have m _rfor interval [a ₂, b ₂] in effective nodes, j _r-1 counts for the acromere of the interval left and right sides; for comprising but do not comprise an interval;

By two-dimentional B-spline function B _j,i(x, r), omits the order j of B-spline, i.e. two-dimentional B-spline function B _j,i(x, r) is designated as B _i(x, r);

Obtaining the three-dimensional instantaneous square root B-spline model exporting PDF based on two-dimentional square root B-spline function is:

$\sqrt{\mathrm{\γ}(x,r,u_k)}=C_0(x,r)V_k+B_n(x,r){\mathrm{\ω}}_n\left(V_k\right)$

Wherein,

γ (x, r, u _k) be three-dimensional output probability density function;

N is the two-dimentional B-spline function number selected, and k is sampling instant;

C ₀(x, r)=[B ₁(x, r), B ₂(x, r) ..., B _n-1(x, r)], wherein, C ₀(x, r) is that 1 × (n-1) ties up Basis Function transformation vector;

B _i(x, r) is two-dimentional B-spline function;

V _k=[ω ₁(u _k), ω ₂(u _k) ..., ω _n-1(u _k)] ^t, wherein, V _kweight vector is tieed up in (n-1) × 1 corresponding for the k moment;

ω _i(u _k) for depending on u _kweights, u _kfor the control action that the k moment is corresponding.

B _n(x, r) is two-dimentional B-spline function;

ω _n(V _k) be the weights that the n-th basis function is corresponding.

Step S2: on the basis of step S1, adds the dynamic change part of weights, obtains exporting PDF dynamic model based on square root B-spline model three-dimensional;

Assuming that the weights dynamic part added is:

V _k＝AV _k-1+Bu _k-1

Wherein, A is that (n-1) × (n-1) of expression system dynamic relationship ties up parameter matrix, and B is that parameter matrix is tieed up in (n-1) × 1 of expression system dynamic relationship; V _k-1the n-1 corresponding for the k-1 moment ties up weight vector, u _k-1for the controlled quentity controlled variable that the k-1 moment is corresponding;

Three-dimensional based on square root B-spline model exports PDF dynamic model:

$V_k={\mathrm{AV}}_{k-1}+{\mathrm{Bu}}_{k-1}\sqrt{\mathrm{\γ}(x,r,u_k)}=C_0(x,r)V_k+B_n(x,r){\mathrm{\ω}}_n\left(V_k\right)$

Step S3: the weight vector V of system _kthere is nonlinear relationship in the weights corresponding with the n-th basis function, analyzes and obtain the three-dimensional dynamic decoupling exported between PDF dynamic model weights;

The decoupling zero formula that described three-dimensional exports between PDF dynamic model weights is:

$\left[\begin{array}{c}{V}_k\\ {\mathrm{\ω}}_n\left(V_k\right)\end{array}\right]=Q^{-1}\left[\begin{array}{c}{C}_1\\ C_2\end{array}\right]$

Wherein,

wherein, C ₁for the integration of k moment output probability density function root mean square and Basis Function transformation matrix product;

wherein, C ₂for k moment output probability density function root mean square and the n-th two dimensional basis functions B _nthe integration of (x, r);

$Q=\left[\begin{array}{cc}{\mathrm{\Σ}}_0& {{\mathrm{\Σ}}_1}^T\\ {\mathrm{\Σ}}_1& {\mathrm{\Σ}}_2\end{array}\right],$ Wherein, Q is ∑ ₀, ∑ ₁, ∑ ₂variation, when basis function select after and the inputoutput data of known real system time, Q is known quantity.

wherein, ∑ ₀for Basis Function transformation vector C ₀the integration of square value in its field of definition of (x, r);

wherein, ∑ ₁for Basis Function transformation vector C ₀(x, r) and the integration of the n-th basis function product in its field of definition scope;

wherein, ∑ ₂be the n-th basis function B _n(x, r) square is at the integration of its field of definition scope;

A ₁for the lower limit of setting in X-axis setting district; b ₁for the higher limit of setting in X-axis setting district;

A ₂for the lower limit of setting in R axle setting district; b ₂for the higher limit of setting in R axle setting district;

∑ ₁ ^tfor ∑ ₁transposed matrix.

Step S4: analyzing three-dimensional exports PDF dynamic model and meets the condition that natural sulfur reservoir possesses and be:

||V _k|| _∑≤1

Wherein, || V _k|| _∑=V _k ^t∑ V _k,

In above two formulas, V _kthe n-1 corresponding for the k moment ties up weight vector; V _k ^tfor V _ktransposed matrix; ∑ is ∑ ₀, ∑ ₁and ∑ ₂conversion vector;

Step 2: utilize the inputoutput data collected in real system to set up the three-dimensional input/output model exporting PDF by recursive least squares;

The input/output model that the three-dimensional set up exports PDF is:

$f(x,r,u_k)=\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}{a}_{i}f(x,r,u_{k-i})+\underset{j=0}{\overset{n-2}{\mathrm{\Σ}}}{C}_0(x,r)D_j{u}_{k-j-1}$

Wherein, $f(x,r,u_k)=\sqrt{\mathrm{\γ}(x,r,u_k)}-B_n(x,r){\mathrm{\ω}}_n\left(V_k\right);$

In above two formulas, f (x, r, u _k) be the variation of output probability density function corresponding to k moment; a _ifor f (x, r, u that the k-i moment is corresponding _k-i) coefficient; F (x, r, u _k-i) be the variation of output probability density function corresponding to k-i moment; u _k-ifor the control action that the k-i moment is corresponding; u _k-j-1for the control action that the k-j-1 moment is corresponding; D _j=[d _j1..., d _ji..., d _{j (n-1)}] ^tfor needing the parameter of identification; d _jifor with C ₀the coefficient that item in (x, r) is corresponding.

Step 3: select instantaneous square root performance index CONTROLLER DESIGN, by the controlled quentity controlled variable of optimization instantaneous square root performance index CONTROLLER DESIGN, the system that realizes exports the shape that the distribution of PDF distribution shape tracing preset exports PDF distribution.

The instantaneous square root performance index of described selection are:

$J={\∫}_{a_1}^{b_2}{\∫}_{a_1}^{b_1}{(\sqrt{\mathrm{\γ}(x,r,u_{k+1})}-\sqrt{g(x,r)})}^2\mathrm{dxdr}+{{\mathrm{Ru}}_k}^2$

Wherein, J is instantaneous square root performance index value; γ (x, r, u _k+1) be three-dimensional output probability density function; G (x, r) is given three-dimensional output PDF distribution function; R is the constraint constant of control action; a ₁for the lower limit of setting in X-axis setting district; b ₁for the higher limit of setting in X-axis setting district; a ₂for the lower limit of setting in R axle setting district; b ₂for the higher limit of setting in R axle setting district.

Controlled quentity controlled variable is obtained as follows by optimization instantaneous square root performance index:

$u_k=-\frac{{\∫}_{a_2}^{b_2}{\∫}_{a_1}^{b_1}\left(C_0(x,r)D_0\hat{g}(x,r)\right)\mathrm{dxdr}}{{\∫}_{a_2}^{b_2}{\∫}_{a_1}^{b_1}{\left(C_0(x,r)D_0\right)}^2\mathrm{dxdr}+R}$

Wherein,

$\hat{g}(x,r)=\underset{i=2}{\overset{n-1}{\mathrm{\Σ}}}(a_{i}f(x,r,k-i+1)+C_0(x,r)D_{i-1}{u}_{k-i+1})+a_{1}f(x,r,k)+B_n(x,r){\mathrm{\ω}}_n\left(V_k\right)-\sqrt{g(x,r)};$

Wherein, for the variation of known quantity and parameter; F (x, r, k-i+1) is the variation of output probability density function corresponding to k-i+1 moment; D ₀for the parameter value picked out.

Described controlled quentity controlled variable u _kby right middle f (x, r, k), f (x, r, k-1) ... f (x, r, k-n+2), ω _n(V _k) and u _k-1, u _k-2..., u _k-n+2the adjustment of value, the system that realizes exports the shape that PDF distribution shape tracing preset exports PDF distribution.

Beneficial effect of the present invention: 1, the present invention has ensured that system exports the PDF constraint condition that is greater than 1, the system that analyzes exports PDF when its field of definition inner product is divided into 1, the constraint condition that weights should meet; 2, establish the three-dimensional dynamic model exporting PDF according to square root B-spline model in the present invention, then the dynamic model based on square root B-spline function set up is converted, according to the inputoutput data gathered, establish the input/output model of system; 3, the present invention devises optimization routine controller, by optimizing square root performance index, obtains the control action of system, realizes exporting the tracking of PDF distribution shape to given output PDF distribution shape.The present invention has enriched three-dimensional output PDF control theory, provides new method for having the three-dimensional industrial process exporting distribution character.

Accompanying drawing explanation

Fig. 1 is two-dimentional B-spline function image;

Fig. 2 is the initial p DF distribution of Three-Dimensional Dynamic system;

Fig. 3 is that the given output PDF of Three-Dimensional Dynamic system distributes;

Fig. 4 is the output PDF response surface design of Three-Dimensional Dynamic system;

Fig. 5 controls to export PDF to given PDF tracking error the last moment;

Fig. 6 is the response curve of controlled quentity controlled variable in control procedure;

Fig. 7 is performance index change curve in control procedure;

Fig. 8 is overall flow figure of the present invention.

Embodiment

In order to deepen the understanding of the present invention, below in conjunction with accompanying drawing, specific embodiments of the invention are described in further detail.It is emphasized that following explanation is only exemplary, the scope be not meant to limit the present invention and application thereof.

In order to actual industrial process needs, PDF Control Lyapunov functions will be exported in the system with distributed in three dimensions characteristic, to simplify the complicacy adopting mechanism method to set up system model and Controller gain variations to cause.The present invention proposes a kind of dynamic modeling and controller design method of three-dimensional output probability density function.For realizing the tracking to whole output PDF distribution shape.

The present invention is divided into the following steps:

One, the three-dimensional instantaneous square root B-spline model exporting PDF is built, the three-dimensional that the dynamic change part that instantaneous square root B-spline model basis adds weights constitutes based on square root B-spline model exports PDF dynamic model, realize the three-dimensional dynamic decoupling exported between PDF dynamic model weights, analyzing three-dimensional output PDF dynamic model meets the condition that natural sulfur reservoir possesses;

Two, on the basis of step one for the ease of CONTROLLER DESIGN, above-mentioned dynamic model is converted, utilizes the inputoutput data that collects in real system to set up the three-dimensional input/output model exporting PDF by recursive least squares;

Three, on the basis of step 2, select instantaneous square root performance index CONTROLLER DESIGN, the control action obtained by adjustment controller, the system that reaches exports the shape that the distribution of PDF distribution shape tracing preset exports PDF distribution.

Specifically be divided into:

1, two-dimentional B-spline function method for expressing

The tensor product representation of two one dimension B-spline functions of two dimension B-spline function:

Wherein, computing formula is obtained by following recursion formula:

$B_{1,i_x}\left(x\right)=\left\{\begin{array}{cc}1,& x\∈[x_{i_x+1},x_{i_x+1})\\ 0,& x\∉[x_{i_x},x_{i_x+1})\end{array}\right.---\left(2\right)$

$B_{j_x,i_x}\left(x\right)=\frac{x-x_{i_x}}{x_{i_x+j_x-1}-x_{i_x}}{B}_{j_x-1,i_x}\left(x\right)+\frac{x_{i_x+j_x}-x}{x_{i_x+j_x}-x_{i_x+1}}{B}_{j_x-1,i_x+1}\left(x\right)---\left(3\right)$

obtained by following recursion formula:

$B_{1,i_r}\left(r\right)=\left\{\begin{array}{cc}1,& r\∈[r_{i_r+1},r_{i_r+1})\\ 0,& r\∉[r_{i_r},r_{i_r+1})\end{array}\right.---\left(4\right)$

$B_{j_r,i_r}\left(r\right)=\frac{r-r_{i_r}}{r_{i_r+j_r-1}-r_{i_r}}{B}_{j_r-1,i_r}\left(r\right)+\frac{r_{i_r+j_r}-r}{r_{i_r+j_r}-r_{i_r+1}}{B}_{j_r-1,i_r+1}\left(r\right)---\left(5\right)$

Wherein, B _j,i(x, r) is two-dimentional B-spline basis function; for one dimension B-spline basis function; for one dimension B-spline basis function; X, r are respectively the variable spatially defined, x ∈ [a ₁, b ₁], r ∈ [a ₂, b ₂]; a ₁for the lower limit of setting in setting district; b ₁for the higher limit of setting in setting district; [a ₁, b ₁] for comprising a ₁and b ₁an interval; a ₂for the lower limit of setting in setting district; b ₂for the higher limit of setting in setting district; [a ₂, b ₂] for comprising a ₂and b ₂an interval; J represents the order of two-dimentional B-spline; I represents two-dimentional B-spline basis function number; j _xfor the order of basis function that X-axis is chosen; i _xfor the number of basis function that X-axis is chosen; j _rfor the order of basis function that R axle is chosen; i _rfor the number of basis function that R axle is chosen;

be 1 rank i-th _xindividual B-spline function; for j _x-1 rank i-th _xindividual B-spline function; for j _x-1 rank i-th _x+ 1 B-spline basis function; for nodal value, and have m _xfor interval [a ₁, b ₁] in effective nodes, j _x-1 counts for the acromere of the interval left and right sides; for comprising with an interval;

be 1 rank i-th _rindividual B-spline function; for j _r-1 rank i-th _rindividual B-spline function; for j _r-1 rank i-th _r+ 1 B-spline basis function; for nodal value, and have m _rfor interval [a ₂, b ₂] in effective nodes, j _r-1 counts for the acromere of the interval left and right sides; for comprising with an interval;

2, the three-dimensional instantaneous square root model exporting PDF is built

Namely square root model approaches the square root that system exports PDF, to ensure the nonnegativity of output PDF in control procedure of stochastic system with two-dimentional B-spline function.

The three-dimensional discrete form exporting the instantaneous square root model of PDF is expressed as:

$\sqrt{\mathrm{\γ}(x,r,u_k)}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i\left(u_k\right)B_i(x,r)+e_0---\left(6\right)$

Wherein,

γ (x, r, u _k) be output probability density function;

N is the two-dimentional B-spline function number selected, and k is sampling instant;

B _i(x, r) is two-dimentional B-spline function, wherein, eliminates the order j of B-spline;

ω _i(u _k) for depending on u _kweights; u _kfor the control action that the k moment is corresponding.

E ₀for the approximate error of system;

Ignore e under normal circumstances ₀.Then the three-dimensional instantaneous square root model representation exporting PDF is:

$\sqrt{\mathrm{\γ}(x,r,u_k)}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_i\left(u_k\right)B_i(x,r)---\left(7\right)$

Export PDF function for given three-dimensional, formula (7) is unique, gets x ∈ [a ₁, b ₁] r ∈ [a ₂, b ₂] be stochastic variable span, the instantaneous square root model according to (7) formula three-dimensional being exported PDF is expressed as further:

$\sqrt{\mathrm{\γ}(x,r,u_k)}=C_0(x,r)V_k+B_n(x,r){\mathrm{\ω}}_n\left(V_k\right)---\left(8\right)$

Wherein,

C ₀(x, r)=[B ₁(x, r), B ₂(x, r) ..., B _n-1(x, r)], wherein, C ₀(x, r) is the conversion vector that 1 × (n-1) ties up basis function, B _n-1(x, r) is (n-1)th basis function; V _k=[ω ₁(u _k), ω ₂(u _k) ..., ω _n-1(u _k)] ^t, wherein, V _kfor weight vector is tieed up in (n-1) × 1, ω _n-1(u _k) be the weights that (n-1)th basis function is corresponding; [ω ₁(u _k), ω ₂(u _k) ..., ω _n-1(u _k)] ^tfor [ω ₁(u _k), ω ₂(u _k) ..., ω _n-1(u _k)] transposed matrix; ω _n(V _k) be the weights that the n-th weights are corresponding; B _n(x, r) is the n-th basis function.

3, the three-dimensional that the dynamic change part adding weights is formed based on square root B-spline model exports PDF dynamic model

The three-dimensional of above-mentioned design exports PDF square root model and does not relate to weights change, and in a lot of situation, exporting between PDF and input is dynamic relationship.General hypothesis V _kwith control inputs u _kbetween be linear dynamic be correlated with, the dynamic change part of hypothesis weights is expressed as here:

V _k＝AV _k-1+Bu _k-1(9)

In formula (9), A is that (n-1) × (n-1) of expression system dynamic relationship ties up parameter matrix, and B is that parameter matrix is tieed up in (n-1) × 1 of expression system dynamic relationship; V _k-1the n-1 corresponding for the k-1 moment ties up weight vector, u _k-1for the controlled quentity controlled variable that the k-1 moment is corresponding.

So the three-dimensional based on square root B-spline model exports PDF dynamic model expression and is:

V _k＝AV _k-1+Bu _k-1(10)

$\sqrt{\mathrm{\γ}(x,r,u_k)}=C_0(x,r)V_k+B_n(x,r){\mathrm{\ω}}_n\left(V_k\right)---\left(11\right)$

4, the three-dimensional dynamic decoupling exported between PDF dynamic model weights;

Draw between the weights that basis function is corresponding it is nonlinear relationship by (10), (11) formula, in order to address this problem, need do as down conversion, to realize the dynamic decoupling of weights.

By formula (11) both sides with being multiplied by [C ₀ ^t(x, r) B _n(x, r)] ^tcan obtain:

$\left[\begin{array}{c}{C_0}^T(x,r)\\ B_n(x,r)\end{array}\right]\sqrt{\mathrm{\γ}(x,r,u_k)}=\left[\begin{array}{cc}{C_0}^T(x,r)C_0(x,r)& {C_0}^T(x,r)B_n(x,r)\\ B_n(x,r)C_0(x,r)& {B_n}^2(x,r)\end{array}\right]\left[\begin{array}{c}{V}_k\\ {\mathrm{\ω}}_n\left(V_k\right)\end{array}\right]---\left(12\right)$

Get x ∈ [a ₁, b ₁] r ∈ [a ₂, b ₂] be stochastic variable span, above formula both sides integration can obtain:

$\left[\begin{array}{c}{C}_1\\ C_2\end{array}\right]=Q\left[\begin{array}{c}{V}_k\\ {\mathrm{\ω}}_n\left(V_k\right)\end{array}\right]---\left(13\right)$

Wherein,

wherein, C ₁for the integration of k moment output probability density function root mean square and Basis Function transformation matrix product;

wherein, C ₂for k moment output probability density function root mean square and the n-th two dimensional basis functions B _nthe integration of (x, r);

$Q=\left[\begin{array}{cc}{\mathrm{\Σ}}_0& {{\mathrm{\Σ}}_1}^T\\ {\mathrm{\Σ}}_1& {\mathrm{\Σ}}_2\end{array}\right],$ Wherein, Q is ∑ ₀, ∑ ₁, ∑ ₂variation, when basis function select after and the inputoutput data of known real system time, Q is known quantity.

wherein, ∑ ₀for Basis Function transformation vector C ₀the integration of square value in its field of definition of (x, r);

wherein, ∑ ₁for Basis Function transformation vector C ₀(x, r) and the integration of the n-th basis function product in its field of definition scope;

wherein, ∑ ₂be the n-th basis function B _n(x, r) square is at the integration of its field of definition scope;

A ₁for the lower limit of setting in X-axis setting district; b ₁for the higher limit of setting in X-axis setting district;

A ₂for the lower limit of setting in R axle setting district; b ₂for the higher limit of setting in R axle setting district;

∑ ₁ ^tfor ∑ ₁transposed matrix.

When the B-spline obtained is orthogonal, in formula (13), Q inverse of a matrix always exists, and formula (13) is expressed as:

$\left[\begin{array}{c}{V}_k\\ {\mathrm{\ω}}_n\left(V_k\right)\end{array}\right]=Q^{-1}\left[\begin{array}{c}{C}_1\\ C_2\end{array}\right]---\left(14\right)$

Formula (14) achieves and realizes ω _n(V _k) and V _kdynamic decoupling.

4, the three-dimensional PDF of output dynamic model meets natural sulfur reservoir condition and derives as follows:

Because γ is (x, r, u _k) be output probability density function, then meet its field of definition scope planted agent the constraint that integration is:

${\∫}_{a_1}^{b_2}{\∫}_{a_1}^{b_1}\mathrm{\γ}(x,r,u_k)\mathrm{dxdr}={\∫}_{a_2}^{b_2}{\∫}_{a_1}^{b_1}{\left(\sqrt{\mathrm{\γ}(x,r,u_k)}\right)}^2\mathrm{dxdr}=1---\left(15\right)$

Release according to formula (15):

V _k ^T∑ ₀V _k+2∑ ₁V _kω _n(u _k)+∑ ₂ω _n ²(u _k)＝1 (16)

Solution formula (16) can obtain

${\mathrm{\ω}}_n\left(u_k\right)=\frac{1}{{\mathrm{\Σ}}_2}(-{\mathrm{\Σ}}_1\±\sqrt{{V_k}^T{{\mathrm{\Σ}}_1}^T{\mathrm{\Σ}}_1{V}_k-{V_k}^T{\mathrm{\Σ}}_0{\mathrm{\Σ}}_2{V}_k-{\mathrm{\Σ}}_2})---\left(17\right)$

ω is drawn by above formula _n(u _k) and V _kbetween be nonlinear relationship, and ω _n(u _k) can represent with other n-1 free weights, note ω _n(u _k)=h (V _k).In order to ensure that above formula has solution, should following formula be met:

V _k ^T∑ ₁ ^T∑ ₁ ^TV _k-V _k ^T∑ ₀∑ ₂V _k-∑ ₂≥0 (18)

Abbreviation formula (18) can obtain following non-linear constrain

||V _k|| _∑≤1 (19)

Wherein, V _kfor weight vector is tieed up in (n-1) × 1; || V _k|| _∑=V _k ^t∑ V _k, represent V _k∑ norm; ∑ is ∑ ₀, ∑ ₁and ∑ ₂conversion vector.

As long as the system weights that calculating is tried to achieve meet the condition that formula (19) possesses, just can meet and export PDF is divided into 1 natural sulfur reservoir condition in its field of definition inner product.

Obtain the three-dimensional PDF dynamic model that exports by above-mentioned derivation to meet the condition that natural sulfur reservoir possesses and be:

||V _k|| _∑≤1

5, the three-dimensional input/output model exporting PDF

The weights dynamic relationship that above-mentioned (9) formula represents is not easy to obtain in systems in practice, so, need that PDF dynamic model is exported to the three-dimensional based on square root B-spline model set up and carry out a down conversion, the dynamic model expression that (10) (11) formula represents is become input/output format.Order

$f(x,r,u_k)=\sqrt{\mathrm{\γ}(x,r,u_k)}-B_n(x,r)h\left(V_k\right)=C_0(x,r)V_k---\left(20\right)$

Wherein, f (x, r, u _k) for having the output function of equal value of output characteristics.Displacement operator z is introduced to weights dynamic equation (10) ^-1, formula (20) is rewritten as:

f(x,r,u _k)＝C ₀(x,r)(I-Az ^-1) ^-1Bu _k-1(21)

Then according to (I-Az ^-1) ^-1the expansion of B, by above formula abbreviation, namely obtains the three-dimensional input/output model exporting PDF:

$f(x,r,u_k)=\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}{a}_{i}f(x,r,u_{k-i})+\underset{j=0}{\overset{n-2}{\mathrm{\Σ}}}{C}_0(x,r)D_j{u}_{k-j-1}---\left(22\right)$

Wherein, I is (n-1) × (n-1) is unit matrix; N is the basis function number chosen; F (x, r, u _k-i) be output function of equal value corresponding to k-i moment; u _k-ifor the control action that the k-i moment is corresponding; u _k-j-1for the control action that the k-j-1 moment is corresponding; a _ifor f (x, r, u that the k-i moment is corresponding _k-i) coefficient, D _j=[d _j1..., d _ji..., d _{j (n-1)}] ^tfor needing the parameter of identification, d _jifor with C ₀the coefficient that item in (x, r) is corresponding.

D _jand a _ibe all the unknown quantity in (22) formula, adopt the unknown parameter in the input/output model of recursive least squares identification three-dimensional output PDF according to the form of (22) formula.Identification process is as follows:

Definition

θ＝[a ₁,…,a _n-1,d ₀₁,…,d _0(n-1),d ₁₁,…，d _1(n-1),…,d _(n-2)1,…,d _(n-2)(n-1)] ^T(23)

φ(x,r,k)＝[f(x,r,k-1),…,f(x,r,k-n+1),u _k-1C ₀₁(x,r),…,u _k-1C _0(n-1)(x,r),…,

u _k-n+1C ₀₁(x,r),…,u _0(k-n+1)C _0(n-1)(x,r)] ^T(24)

Wherein, θ is the weights needing identification; φ (x, r, k) is known quantity, when the inputoutput data collected in real system is determined.

Get x ∈ [a ₁, b ₁] r ∈ [a ₂, b ₂] be stochastic variable span, within the scope of its field of definition, select N respectively _xand N _rindividual sampled point forms f (x _i, r _j, k):

f(x _i,r _j,k)＝θ ^Tφ(x _i,r _j,k) (25)

X _i, r _jfor the sampled point of X-axis and R axle, i=1,2 ..., N _x, j=1,2 ..., N _r.

Recursive least squares is defined as follows:

$\mathrm{\θ}(i+1,j+1)=\mathrm{\θ}(i,j)+\frac{P(i,j)\mathrm{\φ}(x_i,r_j,k)\mathrm{\ϵ}(i,j)}{1+{\mathrm{\φ}}^T(x_i,r_j,k)P(i,j)\mathrm{\φ}(x_i,r_j,k)}---\left(26\right)$

ε(i,j)＝f(x _i,r _j,k)-θ ^T(i,j)φ(x _i,r _j,k) (27)

$P(i+1,j+1)=(I-\frac{P(i,j)\mathrm{\φ}(x_i,r_j,k){\mathrm{\φ}}^T(x_i,r_j,k)}{1+{\mathrm{\φ}}^T(x_i,r_j,k)P(i,j)\mathrm{\φ}(x_i,r_j,k)})P(i,j)---\left(28\right)$

Wherein, P (i, j) is transformation matrix required during identified parameters; ε (i, j) is Identification Errors; P (1,1)=10 ^3-6i _{n (n-1)}for initial matrix, wherein, I _{n (n-1)}for n × (n-1) is unit matrix; θ (1,1)=θ ₀for initial weight vector.

The input/output model parameter identification step of three-dimensional output PDF is as follows:

(1) select suitable two-dimentional B-spline basis function, calculate basis function B _i(x, r) (i=1,2 ..., value n).

(2) at sampling instant k (k>=n), collection system control inputs { u _k-1..., u _k-n+1and field of definition within the scope of sample point γ (x _i, r _j, u _k-1) ..., γ (x _i, r _j, u _k-n+1) value;

(3) h (V is calculated by formula (14) _k-1) ..., h (V _k-n-1) value, calculate f (x by definition _i, r _j, u _k-1) ..., f (x _i, r _j, u _k-n+1) and φ (x _i, r _j, value k);

(4) according to formula (26)-(28), estimated parameter θ, note θ (N _x, N _r) be the estimated value of sampling instant k;

(5) if k is less than N, then k increases by 1, turns to second step.

After identification, θ is known quantity, then obtain the three-dimensional input/output model exporting PDF:

$f(x,r,u_k)=\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}{a}_{i}f(x,r,u_{k-i})+\underset{j=0}{\overset{n-2}{\mathrm{\Σ}}}{C}_0(x,r)D_j{u}_{k-j-1}$

6, instantaneous optimization Tracking Control Design

Export the input/output model of PDF above according to the three-dimensional set up, can corresponding Controller gain variations be carried out.The object of CONTROLLER DESIGN selects suitable control inputs that the output PDF distribution shape of system reality is approached as much as possible to expect PDF distribution shape, consider the model set up, the weights of hypothesized model meet constraint condition, so select square root quadratic performance index:

$J={\∫}_{a_1}^{b_2}{\∫}_{a_1}^{b_1}{(\sqrt{\mathrm{\γ}(x,r,u_{k+1})}-\sqrt{g(x,r)})}^2\mathrm{dxdr}+{{\mathrm{Ru}}_k}^2---\left(29\right)$

Be rewritten as according to formula (20) (22) above formula:

$J\left(u_k\right)={\∫}_{a_2}^{b_2}{\∫}_{a_1}^{b_1}{(C_0(x,r)D_0{u}_k+\hat{g}(x,r))}^2\mathrm{dxdr}+{{\mathrm{Ru}}_k}^2---\left(30\right)$

Wherein:

$\hat{g}(x,r)=\underset{i=2}{\overset{n-1}{\mathrm{\Σ}}}(a_{i}f(x,r,k-i+1)+C_0(x,r)D_{i-1}{u}_{k-i+1})+a_{1}f(x,r,k)+B_n(x,r)h\left(V_k\right)-\sqrt{g(x,r)}$

Expansion (30):

$J\left(u_k\right)={u_k}^2{\∫}_{a_2}^{b_2}{\∫}_{a_1}^{b_1}{\left(C_0(x,r)D_0\right)}^2\mathrm{dxdr}+{2u}_k{\∫}_{a_2}^{b_2}{\∫}_{a_1}^{b_1}{C}_0(x,r)D_0\hat{g}(x,r)\mathrm{dxdr}+{\∫}_{a_2}^{b_2}{\∫}_{a_1}^{b_1}\hat{g}{(x,r)}^2\mathrm{dxdr}+{{\mathrm{Ru}}_k}^2---\left(31\right)$

Adopt the algorithm of optimal control, to u _kask local derviation, order then obtain:

$u_k=-\frac{{\∫}_{a_2}^{b_2}{\∫}_{a_1}^{b_1}\left(C_0(x,r)D_0\hat{g}(x,r)\right)\mathrm{dxdr}}{{\∫}_{a_2}^{b_2}{\∫}_{a_1}^{b_1}{\left(C_0(x,r)D_0\right)}^2\mathrm{dxdr}+R}---\left(32\right)$

Since then, in the present invention, a kind of dynamic modeling of three-dimensional output probability density function and controller design method complete.

Embodiment is as follows:

Because experiment condition is limited, be difficult to obtain the input and output data of real system, dynamic vector A, B of supposing the system are known, and lower surface construction based on the three-dimensional output PDF dynamic model of square root B-spline model is:

V _k＝AV _k-1+Bu _k-1

$\sqrt{\mathrm{\γ}(x,r,u_k)}=C_0(x,r)V_k+B_n(x,r)h\left(V_k\right)$

Wherein

$A=\left(\begin{array}{ccccc}0.991& 0& 0& 0& 0\\ 0& 0.994& 0& 0& 0\\ 0.704& 0.056& 0.887& 0& 0\\ 0& 0& 0& 0.947& 0\\ 0& 0& 0& 0.057& 0.976\end{array}\right);$

B＝[0.0209 0.0448 0.0246 0.0292 0.0305] ^T；

According to the definition of two-dimentional B-spline basis function, one dimension B-spline basis function is selected to construct two-dimentional B-spline function here.Suppose x, the span of r is x ∈ [0,1] r ∈ [0,1], and the basis function on x and r axle is defined as follows:

(1) basis function in X-axis:

B ₁(x)＝xI _x1+(2-x)I _x2

B ₂(x)＝(x-1)I _x2+(3-x)I _x3

B ₃(x)＝(x-2)I _x3+(4-x)I _x4

Wherein

(2) basis function on R axle

$B_1\left(r\right)=\frac{1}{2}{r}^2{I}_{r1}+(-r^2+3r-\frac{3}{2})I_{r2}+\frac{1}{2}{(r-3)}^2{I}_{r3}$

$B_2\left(r\right)=\frac{1}{2}{(r-1)}^2{I}_{r2}+(-r^2+5r-\frac{11}{2})I_{r3}+\frac{1}{2}{(r-4)}^2{I}_{r4}$

Wherein

The known quantity that Controller gain variations needs is:

Initial weight V ₀=[0.688 2.129 1.551 0.166 0.792] ^t;

The input constraint factor is R=0.0005;

Control action span is u ∈ [0,1];

Initial control action is u ₀=0.3;

The control inputs that desired output PDF is corresponding is u=0.65;

Here the two-dimentional B-spline basis function selected as shown in Figure 1.Fig. 2, the initial output PDF that Fig. 3 sets forth three-dimensional linear system distributes and desired output PDF distributed image, Fig. 4 is the system output PDF response surface design that tracing preset exports distribution, last moment desired output PDF exports the tracking error of PDF as shown in Figure 5 with control, the response curve of control inputs effect in control procedure that what Fig. 6 provided is, obtains control inputs by figure and can to restrain and close to expectation input.Fig. 7 is the change curve of performance index in control procedure.

Claims (1)

1. the dynamic modeling of three-dimensional output probability density function and a controller design method, it is characterized in that, the method comprises the following steps:

Step 1: the three-dimensional built based on square root B-spline model exports PDF dynamic model;

Described structure comprises the following steps based on the three-dimensional output PDF dynamic model of square root B-spline model:

Step S1: build the three-dimensional instantaneous square root B-spline model exporting PDF according to two-dimentional B-spline function;

As follows with two one dimension B-spline function tensor product representation two dimension B-spline functions:

$B_{j,i}(x,r)=B_{j_x,i_x}\left(x\right)B_{j_r,i_r}\left(r\right)$

Wherein, obtained by following recursion formula:

$B_{1,i_x}\left(x\right)=\left\{\begin{array}{cc}1,& x\∈[x_{i_x},x_{i_x+1})\\ 0,& x\∉[x_{i_x},x_{i_x+1})\end{array}\right.$

$B_{j_x,i_x}\left(x\right)=\frac{x-x_{i_x}}{x_{i_x+j_x-1}-x_{i_x}}{B}_{j_x-1,i_x}\left(x\right)+\frac{x_{i_x+j_x}-x}{x_{i_x+j_x}-x_{i_x+1}}{B}_{j_x-1,i_x+1}\left(x\right)$

obtained by following recursion formula:

$B_{1,i_r}\left(r\right)=\left\{\begin{array}{cc}1,& r\∈[r_{i_r},r_{i_r+1})\\ 0,& r\∉[r_{i_r},r_{i_r+1})\end{array}\right.$

$B_{j_r,i_r}\left(r\right)=\frac{r-r_{i_r}}{r_{i_r+j_r-1}-r_{i_r}}{B}_{j_r-1,i_r}\left(r\right)+\frac{r_{i_r+j_r}-r}{r_{i_r+j_r}-r_{i_r+1}}{B}_{j_r-1,i_r+1}\left(r\right)$

Wherein, B _j,i(x, r) is two-dimentional B-spline basis function; for one dimension B-spline basis function; for one dimension B-spline basis function; X, r are respectively the variable spatially defined, x ∈ [a ₁, b ₁], r ∈ [a ₂, b ₂]; a ₁for the lower limit of setting in X-axis setting district; b ₁for the higher limit of setting in X-axis setting district; [a ₁, b ₁] for comprising a ₁and b ₁an interval; a ₂for the lower limit of setting in R axle setting district; b ₂for the higher limit of setting in R axle setting district; [a ₂, b ₂] for comprising a ₂and b ₂an interval; J represents the order of two-dimentional B-spline; I represents two-dimentional B-spline basis function number; j _xfor the order of basis function that X-axis is chosen; i _xfor the number of basis function that X-axis is chosen; j _rfor the order of basis function that R axle is chosen; i _rfor the number of basis function that R axle is chosen;

be 1 rank i-th _xindividual B-spline function; for j _x-1 rank i-th _xindividual B-spline function; for j _x-1 rank i-th _x+ 1 B-spline basis function; for nodal value, wherein i _x=-j _x+ 1 ..., m _x+ j _x, j _x>1, and have $x_{j_x+1}<...<a_1=x_0<...<x_{m_x+1}=b_1<...<x_{m_x+j_x};$ M _xfor interval [a ₁, b ₁] in effective nodes, j _x-1 counts for the acromere of the interval left and right sides; for comprising but do not comprise an interval;

be 1 rank i-th _rindividual B-spline function; for j _r-1 rank i-th _rindividual B-spline function; for j _r-1 rank i-th _r+ 1 B-spline basis function; for nodal value, wherein i _r=-j _r+ 1 ..., m _r+ j _r, j _r>1, and have $r_{j_r+1}<...<a_2=r_0<...<r_{m_r+1}=b_2<...<r_{m_r+j_r};$ M _rfor interval [a ₂, b ₂] in effective nodes, j _r-1 counts for the acromere of the interval left and right sides; for comprising but do not comprise an interval;

By two-dimentional B-spline function B _j,i(x, r), omits the order j of B-spline, i.e. two-dimentional B-spline function B _j,i(x, r) is designated as B _i(x, r);

Obtaining the three-dimensional instantaneous square root B-spline model exporting PDF based on two-dimentional square root B-spline function is:

$\sqrt{\mathrm{\&gamma;}(x,r,u_k)}=C_0(x,r)V_k+B_n(x,r){\mathrm{\&omega;}}_n\left(V_k\right)$

Wherein,

γ (x, r, u _k) be three-dimensional output probability density function;

N is the two-dimentional B-spline function number selected, and k is sampling instant;

C ₀(x, r)=[B ₁(x, r), B ₂(x, r) ..., B _n-1(x, r)], wherein, C ₀(x, r) is that 1 × (n-1) ties up Basis Function transformation vector;

B _i(x, r) is two-dimentional B-spline function;

V _k=[ω ₁(u _k), ω ₂(u _k) ..., ω _n-1(u _k)] ^t, wherein, V _kweight vector is tieed up in (n-1) × 1 corresponding for the k moment;

ω _i(u _k) for depending on u _kweights, u _kfor the control action that the k moment is corresponding;

B _n(x, r) is two-dimentional B-spline function;

ω _n(V _k) be the weights that the n-th basis function is corresponding;

Step S2: on the basis of step S1, adds the dynamic change part of weights, obtains exporting PDF dynamic model based on square root B-spline model three-dimensional;

Assuming that the weights dynamic part added is:

V _k＝AV _k-1+Bu _k-1

Wherein, A is that (n-1) × (n-1) of expression system dynamic relationship ties up parameter matrix, and B is that parameter matrix is tieed up in (n-1) × 1 of expression system dynamic relationship; V _k-1the n-1 corresponding for the k-1 moment ties up weight vector, u _k-1for the controlled quentity controlled variable that the k-1 moment is corresponding;

Three-dimensional based on square root B-spline model exports PDF dynamic model:

V _k＝AV _k-1+Bu _k-1

$\sqrt{\mathrm{\&gamma;}(x,r,u_k)}=C_0(x,r)V_k+B_n(x,r){\mathrm{\&omega;}}_n\left(V_k\right)$

Step S3: the weight vector V of system _kthere is nonlinear relationship in the weights corresponding with the n-th basis function, analyzes and obtain the three-dimensional dynamic decoupling exported between PDF dynamic model weights;

The decoupling zero formula that described three-dimensional exports between PDF dynamic model weights is:

$\left[\begin{array}{c}{V}_k\\ {\mathrm{\&omega;}}_n\left(V_k\right)\end{array}\right]=Q^{-1}\left[\begin{array}{c}{C}_1\\ C_2\end{array}\right]$

Wherein,

wherein, C ₁for the integration of k moment output probability density function root mean square and Basis Function transformation matrix product;

wherein, C ₂for k moment output probability density function root mean square and the n-th two dimensional basis functions B _nthe integration of (x, r);

$Q=\left[\begin{array}{cc}{\mathrm{\&Sigma;}}_0& {{\mathrm{\&Sigma;}}_1}^T\\ {\mathrm{\&Sigma;}}_1& {\mathrm{\&Sigma;}}_2\end{array}\right],$ Wherein, Q is Σ ₀, Σ ₁, Σ ₂variation, when basis function select after and the inputoutput data of known real system time, Q is known quantity;

wherein, Σ ₀for Basis Function transformation vector C ₀the integration of square value in its field of definition of (x, r);

wherein, Σ ₁for Basis Function transformation vector C ₀(x, r) and the integration of the n-th basis function product in its field of definition scope;

wherein, Σ ₂be the n-th basis function B _n(x, r) square is at the integration of its field of definition scope;

A ₁for the lower limit of setting in X-axis setting district; b ₁for the higher limit of setting in X-axis setting district;

A ₂for the lower limit of setting in R axle setting district; b ₂for the higher limit of setting in R axle setting district;

Σ ₁ ^tfor Σ ₁transposed matrix;

Step S4: analyzing three-dimensional exports PDF dynamic model and meets the condition that natural sulfur reservoir possesses and be:

||V _k|| _Σ≤1

Wherein, || V _k|| _Σ=V _k ^tΣ V _k, $\mathrm{\&Sigma;}={\mathrm{\&Sigma;}}_0-\frac{{{\mathrm{\&Sigma;}}_1}^T{\mathrm{\&Sigma;}}_1}{{\mathrm{\&Sigma;}}_2};$

In above two formulas, V _kthe n-1 corresponding for the k moment ties up weight vector; V _k ^tfor V _ktransposed matrix; Σ is Σ ₀, Σ ₁and Σ ₂conversion vector;

Step 2: utilize the inputoutput data collected in real system to set up the three-dimensional input/output model exporting PDF by recursive least squares;

The input/output model that the three-dimensional set up exports PDF is:

$f(x,r,u_k)=\underset{i=1}{\overset{n-1}{\mathrm{\&Sigma;}}}{a}_{i}f(x,r,u_{k-i})+\underset{j=0}{\overset{n-2}{\mathrm{\&Sigma;}}}{C}_0(x,r)D_j{u}_{k-j-1}$

Wherein, $f(x,r,u_k)=\sqrt{\mathrm{\&gamma;}(x,r,u_k)}-B_n(x,r){\mathrm{\&omega;}}_n\left(V_k\right);$

In above two formulas, f (x, r, u _k) be the variation of output probability density function corresponding to k moment; a _ifor f (x, r, u that the k-i moment is corresponding _k-i) coefficient; F (x, r, u _k-i) be the variation of output probability density function corresponding to k-i moment; u _k-ifor the control action that the k-i moment is corresponding; u _k-j-1for the control action that the k-j-1 moment is corresponding; D _j=[d _j1..., d _ji..., d _{j (n-1)}] ^tfor needing the parameter of identification; d _jifor with C ₀the coefficient that item in (x, r) is corresponding;

Step 3: select instantaneous square root performance index CONTROLLER DESIGN, by the controlled quentity controlled variable of optimization instantaneous square root performance index CONTROLLER DESIGN, the system that realizes exports the shape that the distribution of PDF distribution shape tracing preset exports PDF distribution;

The instantaneous square root performance index of described selection are:

$J={\&Integral;}_{a_2}^{b_2}{\&Integral;}_{a_1}^{b_1}{(\sqrt{\mathrm{\&gamma;}(x,r,u_{k+1})}-\sqrt{g(x,r)})}^2\mathrm{dxdr}+{{\mathrm{Ru}}_k}^2$

Wherein, J is instantaneous square root performance index value; γ (x, r, u _k+1) be three-dimensional output probability density function; G (x, r) is given three-dimensional output PDF distribution function; R is the constraint constant of control action; a ₁for the lower limit of setting in X-axis setting district; b ₁for the higher limit of setting in X-axis setting district; a ₂for the lower limit of setting in R axle setting district; b ₂for the higher limit of setting in R axle setting district;

Controlled quentity controlled variable is obtained as follows by optimization instantaneous square root performance index:

$u_k=-\frac{{\&Integral;}_{a_2}^{b_2}{\&Integral;}_{a_1}^{b_1}\left(C_0(x,r)D_0\hat{g}(x,r)\right)\mathrm{dxdr}}{{\&Integral;}_{a_2}^{b_2}{\&Integral;}_{a_1}^{b_1}{\left(C_0(x,r)D_0\right)}^2\mathrm{dxdr}+R}$

Wherein,

$\hat{g}(x,r)=\underset{i=2}{\overset{n-1}{\mathrm{\&Sigma;}}}(a_{i}f(x,r,k-i+1)+C_0(x,r)D_{i-1}{u}_{k-i+1})+a_{1}f(x,r,k)+B_n(x,r){\mathrm{\&omega;}}_n\left(V_k\right)-\sqrt{g(x,r)};$

Wherein, for the variation of known quantity and parameter; F (x, r, k-i+1) is the variation of output probability density function corresponding to k-i+1 moment; D ₀for the parameter value picked out;

Described controlled quentity controlled variable u _kby right middle f (x, r, k), f (x, r, k-1) ..., f (x, r, k-n+2), ω _n(V _k) and u _k-1, u _k-2..., u _k-n+2the adjustment of value, the system that realizes exports the shape that PDF distribution shape tracing preset exports PDF distribution.