CN103760931A - Oil-gas-water horizontal type three-phase separator pressure control method optimized through dynamic matrix control - Google Patents

Abstract

The invention discloses an oil-gas-water horizontal type three-phase separator pressure control method optimized through dynamic matrix control. The method comprises the steps that firstly, the model of a pressure object in an oil-gas-water horizontal type three-phase separator is built based on the step response data of the pressure object in the oil-gas-water horizontal type three-phase separator, and basic object characters are excavated out; then, the parameters of a corresponding PI-PD controller are set according to the character of the dynamic matrix control; finally, PI-PD control is carried out on the pressure object in the oil-gas-water horizontal type three-phase separator. The good control performance of the PI-PD control and the good control performance of the dynamic matrix control are combined, and the defects of a traditional control method are effectively overcome.

Description

The oil gas water horizontal three-phase separator compress control method that dynamic matrix control is optimized

Technical field

The invention belongs to technical field of automation, relate to a kind of oil gas water horizontal three-phase separator internal pressure proportional integral-proportion differential (PI-PD) control method of optimizing based on dynamic matrix control (DMC).

Background technology

PID controller architecture is simple, and it is convenient to control, and is widely used in various industrial control systems.But for integration, vibration or unsettled control object, PID is sometimes difficult to meet higher control requirement.For example, when step is inputted, often produce larger hyperharmonic vibration, this may bring potential safety hazard to production.At present, the control of oil gas water horizontal three-phase separator pressure is to adopt PID to control mostly, if can add that PD controls at interior ring, first suppresses its overshoot, and outer shroud adopts PI to control, and will obtain better production performance.Dynamic array control algorithm is as the one of advanced control algorithm, and very low to model requirement, control performance is good simultaneously, if can be by dynamic matrix control and the combination of PI-PD technology, and the efficiency that can further improve oil refining and collect rock gas.

Summary of the invention

The object of the invention is the weak point for existing PID controller, a kind of PI-PD control method of oil gas water horizontal three-phase separator internal pressure of optimizing based on dynamic matrix control is provided, be used for Reducing overshoot, to obtain better working control performance.The method, by controlling in conjunction with dynamic matrix control and PI-PD, has obtained a kind of PI-PD control method with dynamic matrix control performance.The method has not only been inherited the premium properties of dynamic matrix control, the simple needs that also can meet actual industrial process of Simultaneous Forms.

The inventive method first step response data of the pressure object based in oil gas water horizontal three-phase separator is set up the model of oil gas water horizontal three-phase separator internal pressure object, excavates basic plant characteristic; Then according to the characteristic of dynamic matrix control, go the parameter of adjusting corresponding PI-PD controller; Finally the pressure object in oil gas water horizontal three-phase separator being implemented to PI-PD controls.

Technical scheme of the present invention is by data acquisition, set up dynamic matrix, set up forecast model, predict the means such as mechanism, optimization, establish a kind of PI-PD control method of optimizing based on dynamic matrix control, utilized effectively Reducing overshoot improve the stability of system of the method.

The step of the inventive method comprises:

Step (1). by the real-time step response data of process object, set up the model of controlled device, concrete grammar is:

A. give step input signal of controlled device, record the step response curve of controlled device.

B. step response curve a step being obtained carries out filtering processing, then fits to a smooth curve, records step response data corresponding to each sampling instant on smooth curve, and first sampling instant is T _s, adjacent two sampling instant interludes are T _s, sampling instant order is T _s, 2T _s, 3T _sthe step response of controlled device will be at some moment t _nafter=NT, tend to be steady, work as a _i(i > N) and a _nerror and measuring error while having the identical order of magnitude, can think a _nbe approximately equal to the steady-state value of step response.Set up the model vector a of object:

a＝[a ₁,a ₂,…a _N] ^Τ

Wherein Τ is transpose of a matrix symbol, a _ibe the data of process object step response, N is modeling time domain.

Step (2). the PIPD controller of design controlled device, concrete grammar is:

A. set up the dynamic matrix of controlled device

Utilize the model vector a of step (1) b acquisition, set up the dynamic gating matrix of controlled device, its form is as follows:

$A=\left(\begin{array}{cccc}{a}_1& 0& ...& 0\\ a_2& a_1& ...& 0\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ a_P& a_{P-1}& ...& a_{P-M+1}\end{array}\right)$

Wherein, A is P × M rank dynamic matrix of controlled device, the optimization time domain that P is Dynamic array control algorithm, the control time domain that M is Dynamic array control algorithm, M < P < N.

B. calculate the model prediction initial response value y in current k moment of controlled device _m(k)

1.. calculate the model predication value y after k-1 moment access control increment Delta u (k-1) _p(k-1):

y _P(k-1)＝y _M(k-1)+A ₀Δu(k-1)

Wherein,

$y_P(k-1)=\left[\begin{array}{c}{y}_1\left(k\right|k-1)\\ y_1(k+1|k-1)\\ .\\ .\\ .\\ y_1(k+N-1|k-1)\end{array}\right],y_M\left(k\right)=\left[\begin{array}{c}{y}_0\left(k\right|k-1)\\ y_0(k+1|k-1)\\ .\\ .\\ .\\ y_0(k+N-1|k-1)\end{array}\right],A_0=\left[\begin{array}{c}{a}_1\\ a_2\\ .\\ .\\ .\\ a_N\end{array}\right]$

Y ₁(k|k-1), y ₁(k+1|k-1) ..., y ₁(k+N-1|k-1) represent respectively controlled device in the k-1 moment to k, k+1 ..., the model predication value after k+N-1 moment access control increment Delta u (k-1), y ₀(k|k-1), y ₀(k|k-1) ... y ₀(k+N-1|k-1) represent that the k-1 moment is to k, k+1 ..., the initial predicted value in k+N-1 moment, A ₀for the matrix that step response data is set up, Δ u (k-1) is the input control increment in k-1 moment.

2.. calculate the model predictive error value e (k) of k moment controlled device:

ess(k)＝y(k)-y ₁(k|k-1)

Wherein, the real output value of the controlled device that y (k) the expression k moment records, y ₁(k|k-1) represent to have added after controlling increment Δ u (k-1), controlled device is the model predication value to the k moment in the k-1 moment.

3.. calculate the modified value y of k moment model output _cor(k):

y _cor(k)＝y _M(k-1)+h*ess(k)

Wherein,

$y_{\mathrm{cor}}\left(k\right)=\left[\begin{array}{c}{y}_{\mathrm{cor}}\left(k\right|k)\\ y_{\mathrm{cor}}(k+1|k)\\ .\\ .\\ .\\ y_{\mathrm{cor}}(k+N-1|k)\end{array}\right],h=\left[\begin{array}{c}1\\ \mathrm{\&alpha;}\\ .\\ .\\ .\\ \mathrm{\&alpha;}\end{array}\right]$

Y _cor(k|k), y _cor(k+1|k) ... y _cor(k+N-1|k) represent respectively the modified value of controlled device at k moment forecast model, the weight matrix that h is error compensation, α is error correction coefficient.

4.. calculate the initial response value y of the model prediction in k moment _m(k):

y _M(k)＝Sy _cor(k)

Wherein, S is the state-transition matrix on N × N rank,

C. calculate controlled device at M continuous controlling increment Δ u (k) ..., the prediction output valve y under Δ u (k+M-1) _pM, concrete grammar is:

y _PM(k)＝y _p0(k)+AΔu _M(k)

Wherein,

$y_{\mathrm{PM}}\left(k\right)=\left[\begin{array}{c}{y}_M(k+1|k)\\ y_M(k+2|k)\\ .\\ .\\ .\\ y_M(k+P|k)\end{array}\right],y_{P0}\left(k\right)=\left[\begin{array}{c}{y}_0(k+1|k)\\ y_0(k+2|k)\\ .\\ .\\ .\\ y_0(k+P|k)\end{array}\right],{\mathrm{\&Delta;u}}_M\left(k\right)=\left[\begin{array}{c}\mathrm{\&Delta;u}\left(k\right)\\ \mathrm{\&Delta;u}(k+1)\\ .\\ .\\ .\\ \mathrm{\&Delta;u}(k+M-1)\end{array}\right]$

Y _m(k+1|k), y _m(k+2|k) ..., y _m(k+P|k) be the k moment to k+1, k+2 ..., the model prediction output valve in k+P moment, y ₀(k+1|k), y ₀(k+2|k-1) ... y ₀(k+N|k) represent that the k moment is to k+1, k+2 ... the initial predicted value in k+P moment.

D. make the control time domain M=1 of controlled device, choose the objective function J (k) of controlled device, J (k) form is as follows:

$\min J\left(k\right)={{\left|\right|(\mathrm{ref}\left(k\right)-y_{\mathrm{PM}}\left(k\right))\left|\right|}_Q}^2+{\left|\right|\mathrm{\&Delta;u}\left(k\right)\left|\right|}_r^2=Q{(\mathrm{ref}\left(k\right)-y_{P0}\left(k\right)-\mathrm{A\&Delta;u}\left(k\right))}^2+{\mathrm{r\&Delta;u}}^2\left(k\right)$

ref(k)＝[ref ₁(k),ref ₂(k),…,ref _P(k)] ^Τ

Q＝diag(q ₁,q ₂,…q _P)

r＝diag(r ₁,r ₂,…r _M)

ref _i(k)＝β ⁱy(k)+(1-β ⁱ)c(k)

Wherein, Q is error weighting matrix, q ₁, q ₂..., q _pfor the weighting coefficient of weighting matrix; β is softening coefficient, the setting value that c (k) is process object; R is for controlling weighting matrix, r ₁, r ₂... r _mfor controlling the weighting coefficient of weighting matrix, the reference locus that ref (k) is system, ref _i(k) be the value of i reference point in reference locus.

E. controlled quentity controlled variable u (k) is converted:

e(k)＝c(k)-y(k)

u(k)＝u(k-1)+K _p(k)(e(k)-e(k-1))+K _i(k)e(k)-K _f(k)(y(k)-y(k-1)-Kd(y(k)-2y(k-1)+y(k-2))＝u(k-1)+K _p(k)(e(k)-e(k-1))+K _i(k)e(k)-K _f(k)(y(k)-y(k-1)-Kd(y(k)-y(k-1))+Kd(y(k-1)-y(k-2))

U (k) is further processed, can obtain

u(k)＝u(k-1)+w(k) ^ΤE(k)

Wherein,

w(:,k)＝[K _p(k)+K _i(k),-K _p(k),-K _f(k)-K _d(k),K _d(k)] ^Τ

E(k)＝(e(k),e(k-1),y(k)-y(k-1),y(k-1)-y(k-2)) ^Τ

Kp (k), K _i(k), K _f(k), K _d(k) be respectively the ratio of k moment PI-PD controller outer shroud, the integration of outer shroud, the ratio of interior ring, the differential parameter of interior ring, e (k) is the error between k moment reference locus value and real output value, Τ is transpose of a matrix symbol, w (:, k) be four lines k column matrix.

F. objective function u (k) being updated in steps d solves the parameter in PI-PD controller, can obtain:

$w(:,k)=\frac{(\mathrm{ref}\left(k\right)-y_{p0}\left(k\right))\mathrm{QAE}}{(A^T\mathrm{QA}+r)E^{T}E}$

Further can obtain:

K _p(k)＝w(1,k)+w(2,k)

K _i(k)＝-w(2,k)

K _f(k)＝-w(3,k)-w(4,k)

K _d(k)＝w(4,k)

G. obtain the parameter K of PI-PD controller _p(k), K _i(k), K _f(k), K _d(k) after, form controlled quentity controlled variable u (k) and act on controlled device

u(k)＝u(k-1)+K _p(k)(e(k)-e(k-1))+K _i(k)e(k)-K _f(k)(y(k)-y(k-1)-K _d(y(k)-2y(k-1)+y(k-2))。

H. in next moment, according to b, to the step in g, continue to solve the parameter k that PI-PD controller is new _p(k+1), k _i(k+1), k _f(k+1), k _d(k+1) value, successively circulation.

The present invention proposes a kind of PI-PD control method of the oil gas water horizontal three-phase separator internal pressure of optimizing based on dynamic matrix control, the method combines the good control performance of PI-PD control and dynamic matrix control, effectively improve the deficiency of traditional control method, also promoted development and the application of advanced control algorithm simultaneously.

Embodiment

Take the process control of oil gas water horizontal three-phase separator internal pressure as example:

The control of oil gas water horizontal three-phase separator internal pressure is a Delay Process, and regulating measure adopts the aperture of controlling vent valve valve in setting chamber.

Step (1). by the real-time step response data of oil gas water horizontal three-phase separator internal pressure object, set up the model of controlled device, concrete grammar is:

A. step input signal of oil supply air water horizontal three-phase separator, records its step response curve.

B. corresponding step response curve is carried out to filtering processing, then fit to a smooth curve, record step response data corresponding to each sampling instant on smooth curve, first sampling instant is T _s, adjacent two sampling instant interludes are T _s, sampling instant order is T _s, 2T _s, 3T _sresponse will be at some moment t _nafter=NT, tend to be steady, work as a _i(i > N) and a _nerror and measuring error while having the identical order of magnitude, can think a _nbe approximately equal to the steady-state value of step response.Set up the model vector a of oil gas water horizontal three-phase separator internal pressure object:

a＝[a ₁,a ₂,…a _N] ^Τ

Wherein Τ is transpose of a matrix symbol, a _ibe the data of the step response of oil gas water horizontal three-phase separator setting chamber internal pressure, N is modeling time domain.

Step (2). the PI-PD controller of design oil gas water horizontal three-phase separator internal pressure, concrete grammar is:

A. the model vector a that utilizes step (1) b to obtain sets up the dynamic matrix of oil gas water horizontal three-phase separator internal pressure, and its form is as follows:

$A=\left(\begin{array}{cccc}{a}_1& 0& ...& 0\\ a_2& a_1& ...& 0\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ a_P& a_{P-1}& ...& a_{P-M+1}\end{array}\right)$

Wherein, A is P × M rank dynamic matrix of oil gas water horizontal three-phase separator internal pressure, the optimization time domain that P is Dynamic array control algorithm, the control time domain that M is Dynamic array control algorithm, M < P < N.

B. set up the initial model predicted value y of oil gas water horizontal three-phase separator internal pressure in the current k moment _m(k)

1.. calculate the model predication value y of the rear oil gas water horizontal three-phase separator internal pressure of k-1 moment access control increment Delta u (k-1) _p(k-1):

y _P(k-1)＝y _M(k-1)+A ₀Δu(k-1)

Wherein,

$y_P(k-1)=\left[\begin{array}{c}{y}_1\left(k\right|k-1)\\ y_1(k+1|k-1)\\ .\\ .\\ .\\ y_1(k+N-1|k-1)\end{array}\right],A_0=\left[\begin{array}{c}{a}_1\\ a_2\\ .\\ .\\ .\\ a_N\end{array}\right],y_M\left(k\right)=\left[\begin{array}{c}{y}_0\left(k\right|k-1)\\ y_0\left(k\right|k-1)\\ .\\ .\\ .\\ y_0(k+N-1|k-1)\end{array}\right]$

Y ₁(k|k-1), y ₁(k+1|k-1) ..., y ₁(k+N-1|k-1) represent respectively oil gas water horizontal three-phase separator internal pressure in the k-1 moment to k, k+1 ..., the k+N-1 moment adds the model predication value after Δ u (k-1), y ₀(k|k-1), y ₀(k|k-1) ... y ₀(k+N-1|k-1) represent oil gas water horizontal three-phase separator internal pressure in the k-1 moment to k, k+1 ..., the initial predicted value in k+N-1 moment, A ₀for the matrix of being set up by oil gas water horizontal three-phase separator setting chamber internal pressure step response data, Δ u (k-1) is the controlling increment of vent valve valve opening in k-1 moment oil gas water horizontal three-phase separator.

2.. calculate the model predictive error value ess (k) of k moment oil gas water horizontal three-phase separator internal pressure:

ess(k)＝y(k)-y ₁(k|k-1)

Wherein, the real output value of the oil gas water horizontal three-phase separator internal pressure that y (k) the expression k moment records, y ₁(k|k-1) represent to have added after controlling increment Δ u (k-1), oil gas water horizontal three-phase separator internal pressure is the model predication value to the k moment in the k-1 moment.

3.. calculate the modified value y of the pressure model output in k moment oil gas water horizontal three-phase separator _cor(k):

y _cor(k)＝y _M(k-1)+h*ess(k)

Wherein,

$y_{\mathrm{cor}}\left(k\right)=\left[\begin{array}{c}{y}_{\mathrm{cor}}\left(k\right|k)\\ y_{\mathrm{cor}}(k+1|k)\\ .\\ .\\ .\\ y_{\mathrm{cor}}(k+N-1|k)\end{array}\right],h=\left[\begin{array}{c}1\\ \mathrm{\&alpha;}\\ .\\ .\\ .\\ \mathrm{\&alpha;}\end{array}\right]$

Y _cor(k|k), y _cor(k+1|k) ... y _cor(k+N-1|k) represent respectively pressure in the oil gas water horizontal three-phase separator modified value at k moment model, the weight matrix that h is error compensation, α is error correction coefficient.

4.. the pressure in calculating oil gas water horizontal three-phase separator is at the model prediction initial response value y in k moment _m(k):

y _M(k)＝Sy _cor(k)

Wherein, S is the state-transition matrix on N × N rank,

C. calculate pressure in oil gas water horizontal three-phase separator at M continuous controlling increment Δ u (k) ..., the prediction output valve y under Δ u (k+M-1) _pM, concrete grammar is:

y _PM(k)＝y _P0(k)+AΔu _M(k)

Wherein,

$y_{\mathrm{PM}}\left(k\right)=\left[\begin{array}{c}{y}_M(k+1|k)\\ y_M(k+2|k)\\ .\\ .\\ .\\ y_M(k+P|k)\end{array}\right],y_{P0}\left(k\right)=\left[\begin{array}{c}{y}_0(k+1|k)\\ y_0(k+2|k)\\ .\\ .\\ .\\ y_0(k+P|k)\end{array}\right],{\mathrm{\&Delta;u}}_M\left(k\right)=\left[\begin{array}{c}\mathrm{\&Delta;u}\left(k\right)\\ \mathrm{\&Delta;u}(k+1)\\ .\\ .\\ .\\ \mathrm{\&Delta;u}(k+M-1)\end{array}\right]$

Y _p0(k) be y _m(k) front P item, y _m(k+1|k), y _m(k+2|k) ..., y _m(k+P|k) be pressure in oil gas water horizontal three-phase separator in the k moment to k+1, k+2 ..., the model prediction output valve in k+P moment.

D. time domain M=1 is controlled in order, and chooses the objective function J (k) of oil gas water horizontal three-phase separator internal pressure, and J (k) form is as follows:

$\min J\left(k\right)={{\left|\right|(\mathrm{ref}\left(k\right)-y_{\mathrm{PM}}\left(k\right))\left|\right|}_Q}^2+{\left|\right|\mathrm{\&Delta;u}\left(k\right)\left|\right|}_r^2=Q{(\mathrm{ref}\left(k\right)-y_{P0}\left(k\right)-\mathrm{A\&Delta;u}\left(k\right))}^2+{\mathrm{r\&Delta;u}}^2\left(k\right)$

ref(k)＝[ref ₁(k),ref ₂(k),…,ref _P(k)] ^Τ

Q＝diag(q ₁,q ₂,…q _P)

r＝diag(r ₁,r ₂,…r _M)

ref _i(k)＝β ⁱy(k)+(1-β ⁱ)c(k)

Wherein, Q is error weighting matrix, q ₁, q ₂..., q _pfor the weighting coefficient of error weighting matrix; β is softening coefficient, and c (k) is the setting value of oil gas water horizontal three-phase separator internal pressure; R=diag (r ₁, r ₂... r _m) for controlling weighting matrix, r ₁, r ₂... r _mfor controlling the weighting coefficient of weighting matrix; Ref (k) is the reference locus of k moment oil gas water horizontal three-phase separator internal pressure, ref _i(k) be the value of i reference point in reference locus.

E. the controlled quentity controlled variable u (k) of vent valve valve opening in oil gas water horizontal three-phase separator is converted:

e(k)＝c(k)-y(k)

u(k)＝u(k-1)+K _p(k)(e(k)-e(k-1))+K _i(k)e(k)-K _f(k)(y(k)-y(k-1)-Kd(y(k)-2y(k-1)+y(k-2))＝u(k-1)+K _p(k)(e(k)-e(k-1))+K _i(k)e(k)-K _f(k)(y(k)-y(k-1)-Kd(y(k)-y(k-1))+Kd(y(k-1)-y(k-2))

U (k) is further processed, can obtain

u(k)＝u(k-1)+w(k) ^ΤE(k)

Wherein,

w(:,k)＝[K _p(k)+K _i(k),-K _p(k),-K _f(k)-K _d(k),K _d(k)]

E(k)＝(e(k),e(k-1),y(k)-y(k-1),y(k-1)-y(k-2)) ^Τ

Kp (k), K _i(k), K _f(k), K _d(k) be respectively the ratio of the ratio of PI-PD controller outer shroud, the integration of outer shroud, interior ring, the differential parameter of interior ring, e (k) is the error between k moment reference locus value and real output value, Τ is transpose of a matrix symbol, w (:, k) be four lines k column matrix.

F. u (k) is updated in the objective function in steps d, solves the parameter in PI-PD controller, can obtain:

$w(:,k)=\frac{(\mathrm{ref}\left(k\right)-y_{p0}\left(k\right))\mathrm{QAE}}{(A^T\mathrm{QA}+r)E^{T}E}$

Further can obtain:

K _p(k)＝w(1,k)+w(2,k)

K _i(k)＝-w(2,k)

K _f(k)＝-w(3,k)-w(4,k)

K _d(k)＝w(4,k)

G. obtain the parameter K of PI-PD controller _p(k), K _i(k), K _f(k), K _d(k) after, form controlled quentity controlled variable u (k) and act on oil gas water horizontal three-phase separator

u(k)＝u(k-1)+K _p(k)(e(k)-e(k-1))+K _i(k)e(k)-K _f(k)(y(k)-y(k-1)-Kd(y(k)-2y(k-1)+y(k-2))

H. in next moment, according to b, to the step in g, continue to solve the parameter k that PI-PD controller is new _p(k+1), k _i(k+1), k _f(k+1), k _d(k+1) value, acts on controlled device, and circulation successively.

Claims (1)

1. the oil gas water horizontal three-phase separator compress control method that dynamic matrix control is optimized, is characterized in that the concrete steps of the method are:

Step (1). by the real-time step response data of process object, set up the model of controlled device, concrete grammar is:

1-a., to step input signal of controlled device, records the step response curve of controlled device;

The step response curve that 1-b. obtains step 1-a carries out filtering processing, then fits to a smooth curve, records step response data corresponding to each sampling instant on smooth curve, and first sampling instant is T _s, adjacent two sampling instant interludes are T _s, sampling instant order is T _s, 2T _s, 3T _sthe step response of controlled device will be at some moment t _nafter=NT, tend to be steady, work as a _iwith a _nerror and measuring error while having the identical order of magnitude, can think a _nbe approximately equal to the steady-state value of step response, i > N; Set up the model vector a of object:

a＝[a ₁,a ₂,…a _N] ^Τ

Wherein Τ is transpose of a matrix symbol, a _ibe the data of process object step response, N is modeling time domain;

Step (2). the PIPD controller of design controlled device, concrete grammar is:

2-a sets up the dynamic matrix of controlled device

Utilize the model vector a of step 1-b acquisition, set up the dynamic gating matrix of controlled device, its form is as follows:

$A=\left(\begin{array}{cccc}{a}_1& 0& ...& 0\\ a_2& a_1& ...& 0\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ a_P& a_{P-1}& ...& a_{P-M+1}\end{array}\right)$

Wherein, A is P × M rank dynamic matrix of controlled device, the optimization time domain that P is Dynamic array control algorithm, the control time domain that M is Dynamic array control algorithm, M < P < N;

2-b. calculates the model prediction initial response value y in current k moment of controlled device _m(k)

1.. calculate the model predication value y after k-1 moment access control increment Delta u (k-1) _p(k-1):

y _P(k-1)＝y _M(k-1)+A ₀Δu(k-1)

Wherein,

$y_P(k-1)=\left[\begin{array}{c}{y}_1\left(k\right|k-1)\\ y_1(k+1|k-1)\\ .\\ .\\ .\\ y_1(k+N-1|k-1)\end{array}\right],y_M\left(k\right)=\left[\begin{array}{c}{y}_0\left(k\right|k-1)\\ y_0(k+1|k-1)\\ .\\ .\\ .\\ y_0(k+N-1|k-1)\end{array}\right],A_0=\left[\begin{array}{c}{a}_1\\ a_2\\ .\\ .\\ .\\ a_N\end{array}\right]$

Y ₁(k|k-1), y ₁(k+1|k-1) ..., y ₁(k+N-1|k-1) represent respectively controlled device in the k-1 moment to k, k+1 ..., the model predication value after k+N-1 moment access control increment Delta u (k-1), y ₀(k|k-1), y ₀(k|k-1) ... y ₀(k+N-1|k-1) represent that the k-1 moment is to k, k+1 ..., the initial predicted value in k+N-1 moment, A ₀for the matrix that step response data is set up, Δ u (k-1) is the input control increment in k-1 moment;

2.. calculate the model predictive error value e (k) of k moment controlled device:

ess(k)＝y(k)-y ₁(k|k-1)

Wherein, the real output value of the controlled device that y (k) the expression k moment records, y ₁(k|k-1) represent to have added after controlling increment Δ u (k-1), controlled device is the model predication value to the k moment in the k-1 moment;

3.. calculate the modified value y of k moment model output _cor(k):

y _cor(k)＝y _M(k-1)+h*ess(k)

Wherein,

$y_{\mathrm{cor}}\left(k\right)=\left[\begin{array}{c}{y}_{\mathrm{cor}}\left(k\right|k)\\ y_{\mathrm{cor}}(k+1|k)\\ .\\ .\\ .\\ y_{\mathrm{cor}}(k+N-1|k)\end{array}\right],h=\left[\begin{array}{c}1\\ \mathrm{\&alpha;}\\ .\\ .\\ .\\ \mathrm{\&alpha;}\end{array}\right]$

Y _cor(k|k), y _cor(k+1|k) ... y _cor(k+N-1|k) represent respectively the modified value of controlled device at k moment forecast model, the weight matrix that h is error compensation, α is error correction coefficient;

4.. calculate the initial response value y of the model prediction in k moment _m(k):

y _M(k)＝Sy _cor(k)

Wherein, S is the state-transition matrix on N × N rank,

2-c. calculates controlled device at M continuous controlling increment Δ u (k) ..., the prediction output valve y under Δ u (k+M-1) _pM, concrete grammar is:

y _PM(k)＝y _p0(k)+AΔu _M(k)

Wherein,

$y_{\mathrm{PM}}\left(k\right)=\left[\begin{array}{c}{y}_M(k+1|k)\\ y_M(k+2|k)\\ .\\ .\\ .\\ y_M(k+P|k)\end{array}\right],y_{P0}\left(k\right)=\left[\begin{array}{c}{y}_0(k+1|k)\\ y_0(k+2|k)\\ .\\ .\\ .\\ y_0(k+P|k)\end{array}\right],{\mathrm{\&Delta;u}}_M\left(k\right)=\left[\begin{array}{c}\mathrm{\&Delta;u}\left(k\right)\\ \mathrm{\&Delta;u}(k+1)\\ .\\ .\\ .\\ \mathrm{\&Delta;u}(k+M-1)\end{array}\right]$

Y _m(k+1|k), y _m(k+2|k) ..., y _m(k+P|k) be the k moment to k+1, k+2 ..., the model prediction output valve in k+P moment, y ₀(k+1|k), y ₀(k+2|k-1) ... y ₀(k+N|k) represent that the k moment is to k+1, k+2 ... the initial predicted value in k+P moment;

2-d. makes the control time domain M=1 of controlled device, chooses the objective function J (k) of controlled device, and form is as follows:

$\min J\left(k\right)={{\left|\right|(\mathrm{ref}\left(k\right)-y_{\mathrm{PM}}\left(k\right))\left|\right|}_Q}^2+{\left|\right|\mathrm{\&Delta;u}\left(k\right)\left|\right|}_r^2=Q{(\mathrm{ref}\left(k\right)-y_{P0}\left(k\right)-\mathrm{A\&Delta;u}\left(k\right))}^2+{\mathrm{r\&Delta;u}}^2\left(k\right)$

ref(k)＝[ref ₁(k),ref ₂(k),…,ref _P(k)] ^Τ

Q＝diag(q ₁,q ₂,…q _P)

r＝diag(r ₁,r ₂,…r _M)

ref _i(k)＝β ⁱy(k)+(1-β ⁱ)c(k)

Wherein, Q is error weighting matrix, q ₁, q ₂..., q _pfor the weighting coefficient of weighting matrix; β is softening coefficient, the setting value that c (k) is process object; R is for controlling weighting matrix, r ₁, r ₂... r _mfor controlling the weighting coefficient of weighting matrix, the reference locus that ref (k) is system, ref _i(k) be the value of i reference point in reference locus;

2-e. converts controlled quentity controlled variable u (k):

e(k)＝c(k)-y(k)

u(k)＝u(k-1)+K _p(k)(e(k)-e(k-1))+K _i(k)e(k)-K _f(k)(y(k)-y(k-1)-Kd(y(k)-2y(k-1)+y(k-2))＝u(k-1)+K _p(k)(e(k)-e(k-1))+K _i(k)e(k)-K _f(k)(y(k)-y(k-1)-Kd(y(k)-y(k-1))+Kd(y(k-1)-y(k-2))

U (k) is further processed, can obtain

u(k)＝u(k-1)+w(k) ^ΤE(k)

Wherein,

w(:,k)＝[K _p(k)+K _i(k),-K _p(k),-K _f(k)-K _d(k),K _d(k)] ^Τ

E(k)＝(e(k),e(k-1),y(k)-y(k-1),y(k-1)-y(k-2)) ^Τ

Kp (k), K _i(k), K _f(k), K _d(k) be respectively the ratio of k moment PI-PD controller outer shroud, the integration of outer shroud, the ratio of interior ring, the differential parameter of interior ring, e (k) is the error between k moment reference locus value and real output value, Τ is transpose of a matrix symbol, w (:, k) be four lines k column matrix;

The objective function that 2-f. is updated to u (k) in step 2-d solves the parameter in PI-PD controller, can obtain:

$w(:,k)=\frac{(\mathrm{ref}\left(k\right)-y_{p0}\left(k\right))\mathrm{QAE}}{(A^T\mathrm{QA}+r)E^{T}E}$

Further obtain:

K _p(k)＝w(1,k)+w(2,k)

K _i(k)＝-w(2,k)

K _f(k)＝-w(3,k)-w(4,k)

K _d(k)＝w(4,k)

2-g. obtains the parameter K of PI-PD controller _p(k), K _i(k), K _f(k), K _d(k) after, form controlled quentity controlled variable u (k) and act on controlled device

u(k)＝u(k-1)+K _p(k)(e(k)-e(k-1))+K _i(k)e(k)-K _f(k)(y(k)-y(k-1)-K _d(y(k)-2y(k-1)+y(k-2))；

2-h., in next moment, continues to solve to 2-g the parameter k that PI-PD controller is new according to step 2-b _p(k+1), k _i(k+1), k _f(k+1), k _d(k+1) value, successively circulation.