CN106200379A - A kind of distributed dynamic matrix majorization method of Nonself-regulating plant - Google Patents

Abstract

The invention discloses a kind of distributed dynamic matrix majorization method of Nonself-regulating plant.The present invention first passes through the matrix model vector of the multivariable process gathering step response data foundation containing Nonself-regulating plant, then the on-line optimization implementation issue of multivariable process changes into the optimal enforcement problem of each small-scale subsystem.Then suitable performance indications are chosen, the Nash optimization solution of each intelligent body is obtained by continuous iteration, and then obtain the parameter of each intelligent body dynamic matrix controller, each intelligent body is implemented the instant control law in this moment again, and time domain is rolled to subsequent time, repeat above-mentioned optimization process, thus complete the optimization task of whole system.The present invention is on the premise of ensureing relatively high control precision and stability, it is possible to effectively compensate for tradition DDMC method deficiency in the multivariable process containing Nonself-regulating plant controls, and meets the demand of actual industrial process.

Description

A kind of distributed dynamic matrix majorization method of Nonself-regulating plant

Technical field

The invention belongs to technical field of automation, relate to the distributed dynamic matrix majorization (DDMC) of a kind of Nonself-regulating plant Method.

Background technology

Real process is widely present the large scale system of large amount of complex higher-dimension, uses centralized integrated solution to meter The performance of calculation machine and processing speed etc. often require that the highest, disagree with the economy that must take in actual industrial system.Point Cloth dynamic matrix control (DDMC) is led to as a Main Branches of Distributed Predictive Control (DMPC), comprehensive utilization computer Letter technology and control theory, be distributed to the line solver problem of a complicated large scale system in subsystems go distribution real Existing, effectively reduce scale and the complexity of problem, can well control to exist multivariate, close coupling, uncertain controlled right As, improve system control performance.But in actual industrial process, there is many multivariable processes containing Nonself-regulating plant, Such as a part of storage tank, boiler drum level, rectifying column liquid level etc..It is typical owing to the transmission function of Nonself-regulating plant containing Integral element, and then cause controlled device response under definite value step to tend to infinite, this allows for tradition DDMC algorithm cannot Directly application.If able to tradition DDMC method is improved in real process, just can effectively make up tradition DDMC method Deficiency in the multivariable process containing Nonself-regulating plant controls so that DMPC is extended the most further and sends out Exhibition.

Summary of the invention

The present invention seeks to for tradition DDMC method containing Nonself-regulating plant multivariable process control in deficiency it Place, it is proposed that a kind of DDMC method of Nonself-regulating plant.

The method first passes through the matrix model of the multivariable process gathering step response data foundation containing Nonself-regulating plant Vector, excavates basic plant characteristic, then the on-line optimization implementation issue of multivariable process changes into each son on a small scale The optimal enforcement problem of system, and combine the theory in multiple agent and thought, each subsystem under network environment is regarded as It is an intelligent body, between each intelligent body, carries out material, energy and information communication by network.Then by by one for non- Method and a kind of new error calibration method of self regulating plant improvement transfer matrix combine, and choose suitable performance indications, Obtained the Nash optimization solution of each intelligent body by continuous iteration based on receiving the thought of assorted optimization, and then obtain the dynamic square of each intelligent body The parameter of battle array controller, more each intelligent body is implemented the instant control law in this moment, and time domain is rolled to subsequent time, weight Multiple above-mentioned optimization process, thus complete the optimization task of whole system.

The technical scheme is that and set up by data acquisition, model, predict the means such as mechanism, optimization, establish one Plant the distributed dynamic matrix majorization method of Nonself-regulating plant, utilize the method before ensureing relatively high control precision and stability Put, it is possible to effectively compensate for tradition DDMC method deficiency in the multivariable process containing Nonself-regulating plant controls, and meet The demand of actual industrial process.

The step of the inventive method includes:

Step 1. sets up corresponding dynamic matrix model vector by the real-time step response data of Nonself-regulating plant, specifically Method is:

1.1 according to Distributed Predictive Control thought, by the large scale system dispersion of a N input N output Nonself-regulating plant For N number of intelligent body subsystem；

1.2 under steady state operating conditions, for input, i-th intelligent body output is carried out step with jth intelligent body controlled quentity controlled variable Response experiment, records jth (1≤j≤N) the individual input step response curve to i-th (1≤i≤N) individual output respectively；

1.3 step response curves step 1.2 obtained are filtered processing, and then fit to a smooth curve, note The step response data that on record smooth curve, each sampling instant is corresponding, first sampling instant is T_s, during adjacent two samplings The interval time carved is T_s, sampling instant order is T_s、2T_s、3T_s……；The step response data of controlled device will be at some Moment t_L=I_ijT_sStart to present and determine slope rising, with the data in this momentFor starting point, data before are denoted as respectivelySet up jth input to the step response model vector a between i-th output_ij:

$a_{ij}={\[a_1^{ij},a_2^{ij},\mathrm{...},a_{I_{ij}-1}^{ij},a_{I_{ij}}^{ij}+\mathrm{\δ},a_{I_{ij}}^{ij}+2\mathrm{\δ},...,a_{L_{ij}}^{ij}\]}^T$

$a_{L_{ij}}^{ij}=a_{I_{ij}}^{ij}+(L_{ij}-I_{ij})\mathrm{\δ}$

Wherein T is the transposition symbol of matrix, δ be step response data be constant-slope rise after adjacent two data it Between constant difference, L_ijModel length i-th exported for the jth input set, L_ij≥I_ij+1。

Step 2. designs the dynamic matrix controller of i-th intelligent body, and concrete grammar is:

2.1 utilize the model vector a that step 1 obtains_ijSetting up the dynamic matrix of controlled device, its form is as follows:

Wherein A_ijFor jth intelligent body input P × M rank dynamic matrix to the output of i-th intelligent body, a_ijK () is jth The step response data that i-th is exported by individual input, P is the optimization time domain of Dynamic array control algorithm, and M is dynamic matrix control The control time domain of algorithm, L_ij=L (1≤i≤3,1≤j≤3), M < P < L, N are input and output number；

2.2 model prediction initial communication values y obtaining the i-th intelligent body current k moment_i,0(k)

First, controlling increment △ u is added in the k-1 moment₁(k-1),△u₂(k-1),…,△u_n(k-1) i-th intelligence, is obtained The model predication value y of energy body_i,P(k-1):

$y_{i,P}(k-1)=y_{i,0}(k-1)+A_{ii,0}{\mathrm{\&Delta;u}}_i(k-1)+\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij,0}{\mathrm{\&Delta;u}}_j(k-1)$

Wherein,

y_i,P(k-1)=[y_i,1(k|k-1),y_i,1(k+1|k-1),…,y_i,1(k+L-1|k-1)]^T

y_i,0(k-1)=[y_i,0(k|k-1),y_i,0(k+1|k-1),…,y_i,0(k+L-1|k-1)]^T,

A_ii,0=[a_ii(1),a_ii(2),…,a_ii(L)]^T,A_ij,0=[a_ij(1),a_ij(2),…,a_ij(L)]^T

y_i,1(k|k-1),y_i,1(k+1|k-1),…,y_i,1(k+L-1 | k-1) represent that i-th intelligent body is in the k-1 moment respectively To k, k+1 ..., the model predication value in k+L-1 moment, y_i,0(k|k-1),y_i,0(k+1|k-1),…,y_i,0(k+L-1 | k-1) table Show the k-1 moment to k, k+1 ..., the initial prediction in k+L-1 moment, A_ii,0,A_ij,0It is respectively i-th intelligent body and jth intelligence The matrix that the step response data to the output of i-th intelligent body is set up, △ u can be inputted by body₁(k-1),△u₂(k-1),…,△u_n (k-1) it is the input controlled quentity controlled variable of k-1 moment each intelligent body；

It is then possible to obtain model predictive error value e of k moment i-th intelligent body_i(k):

e_i(k)=y_i(k)-y_i,1(k|k-1)

Wherein y_iThe real output value of k i-th intelligent body that () expression k moment records；

Obtain k moment revised model output valve y further_i,cor(k):

y_i,cor(k)=y_i,0(k-1)+h₁*e_i(k)+h₂*e_i(k)

Wherein,

y_i,cor(k)=[y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k+L-1|k)]^T,

h₁=[1, α ..., α]^T,h₂=[0,1 ..., L-1]^T

y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k+L-1 | k) represent that i-th intelligent body is at k moment model respectively Correction value, h₁And h₂For the weight matrix of error compensation, α is error correction coefficient, 0 < α≤1；

Finally obtain initial communication value y of the model prediction in i-th intelligent body k moment_i,0(k):

y_i,0(k)=Sy_i,cor(k)

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

2.3 obtain i-th intelligent body at M continuous print controlling increment △ u according to step 2.1_i(k),△u_i(k+1),…, △u_i(k+M-1) prediction output valve y under_i,PM, concrete grammar is:

$y_{i,PM}\left(k\right)=y_{i,P0}\left(k\right)+A_{ii}{\mathrm{\&Delta;u}}_{i,M}\left(k\right)+\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij}{\mathrm{\&Delta;u}}_{j,M}\left(k\right)$

Wherein,

y_i,PM(k)=[y_i,M(k+1|k),y_i,M(k+2|k),…,y_i,M(k+P|k)]^T

y_i,P0(k)=[y_i,0(k+1|k),y_i,0(k+2|k),…,y_i,0(k+P|k)]^T

△u_i,M(k)=[△ u_i(k),△u_i(k+1),…,△u_i(k+M-1)]^T

△u_j,M(k)=[△ u_j(k),△u_j(k+1),…,△u_j(k+M-1)]^T

y_i,P0K () is y_i,0The front P item of (k), y_i,0(k+1|k),y_i,0(k+2|k),…,y_i,0(k+P | k) it is that the k moment is to k+ 1, k+2 ..., the model prediction output valve in k+P moment；

2.4 performance indications J setting up Nonself-regulating plant i-th intelligent body dynamic matrix controller_i(k) and reference locus ω_i (k), form is as follows:

minJ_i(k)=(ω_i(k)-y_i,PM(k))^TQ_i(ω_i(k)-y_i,PM(k))+△u_i,M(k)^TR_i△u_i,M(k)

ω_i(k)=[ω_i(k+1),ω_i(k+2),…,ω_i(k+P)]^T

ω_i(k+ ε)=β^εy(k)+(1-β^ε) c (k) (ε=1,2 ..., P)

WhereinFor error weighting matrix,For controlling weighting square Battle array,WithIt is respectively Q_i,R_iIn weight coefficient, ω_iK () is the reference locus of i-th intelligent body, β Softening coefficient for reference locus；

The thought of 2.5 foundation Nash optimization, is obtained receiving of i-th intelligent body current k moment by performance indications in step 2.4 Assorted optimal solution:

${\mathrm{\&Delta;u}}_{i,M}^{*}\left(k\right)=D_{ii}\left({\mathrm{\&omega;}}_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij}{\mathrm{\&Delta;u}}_{j,M}^{*}(k\left)\right)$

Wherein:

2.6 can be obtained by step 2.2 to 2.5 in the new round iteration optimal solution of k moment intelligent body i be:

${\mathrm{\&Delta;u}}_{i,M}^{l+1}\left(k\right)=D_{ii}\left({\mathrm{\&omega;}}_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij}{\mathrm{\&Delta;u}}_{j,M}^l(k\left)\right)$

Obtain the whole system optimal control law in the k moment further:

${\mathrm{\&Delta;u}}_M^{l+1}\left(k\right)=D_1\left(\mathrm{\&omega;}\right(k)-y_{P0}(k\left)\right)-D_0{\mathrm{\&Delta;u}}_M^l\left(k\right)$

Wherein:

${\mathrm{\&Delta;u}}_M^{l+1}\left(k\right)={\&lsqb;{\mathrm{\&Delta;u}}_{1,M}^{l+1}\left(k\right),{\mathrm{\&Delta;u}}_{2,M}^{l+1}\left(k\right),...,{\mathrm{\&Delta;u}}_{n,M}^{l+1}\left(k\right)\&rsqb;}^T$

${\mathrm{\&Delta;u}}_M^l\left(k\right)={\&lsqb;{\mathrm{\&Delta;u}}_{1,M}^l\left(k\right),{\mathrm{\&Delta;u}}_{2,M}^l\left(k\right),...,{\mathrm{\&Delta;u}}_{n,M}^l\left(k\right)\&rsqb;}^T$

ω (k)=[ω₁(k),ω₂(k),…,ω_n(k)]^T, y_P0(k)=[y_1,P0(k),y_2,P0(k),…,y_n,P0(k)]^T

2.7 using the Nash optimization solution first term in i-th intelligent body k moment as instant control law △ u_iK (), obtains intelligent body Actual controlled quentity controlled variable u of i_i(k)=u_i(k-1)+△u_iK () acts on i-th intelligent body；

2.8 at subsequent time, repeats step 2.2 to 2.7 and continues to solve the instant control law △ u of i-th intelligent body_i(k+ 1), and then obtain the optimal solution △ u (k+1) of whole system, and circulate successively.

The present invention proposes a kind of DDMC method of Nonself-regulating plant.The method, will on the basis of tradition DDMC method A kind of method for Nonself-regulating plant improvement transfer matrix and a kind of new error calibration method combine, and are ensureing higher control On the premise of precision and stability, effectively compensate for tradition DDMC method in the multivariable process containing Nonself-regulating plant controls Deficiency, and meet the demand of actual industrial process.

Detailed description of the invention

As a example by general predictive control:

Drum Water Level Control System for Boiler is the multivariate Nonself-regulating plant of a typical band integral element, regulating measure Use and control water-supply valve valve opening.

Step 1. by the real-time step response data of boiler drum level object set up corresponding dynamic matrix model to Amount, concrete grammar is:

1.1 according to Distributed Predictive Control thought, by the extensive system of one 3 input 3 output boiler drum level object System is separated into 3 subsystems；

1.2 under steady state operating conditions, enters i-th boiler drum level with jth boiler feed valve valve opening for input Row step response is tested, and records jth (1≤j≤3) the individual input step response curve to i-th (1≤i≤3) individual output respectively；

1.3 step response curves step 1.2 obtained are filtered processing, and then fit to a smooth curve, note The step response data that on record smooth curve, each sampling instant is corresponding, first sampling instant is T_s, during adjacent two samplings The interval time carved is T_s, sampling instant order is T_s、2T_s、3T_s……；The step response data of boiler drum level will be at certain One moment t_L=I_ijT_sStart to present and determine slope rising, with the data in this momentFor starting point, data before are remembered respectively DoSet up the input of jth boiler to the step response model vector a between the output of i-th boiler_ij:

$a_{ij}={\&lsqb;a_1^{ij},a_2^{ij},\mathrm{...},a_{I_{ij}-1}^{ij},a_{I_{ij}}^{ij}+\mathrm{\&delta;},a_{I_{ij}}^{ij}+2\mathrm{\&delta;},...,a_{L_{ij}}^{ij}\&rsqb;}^T$

$a_{L_{ij}}^{ij}=a_{I_{ij}}^{ij}+(L_{ij}-I_{ij})\mathrm{\&delta;}$

Wherein T is the transposition symbol of matrix, δ be step response data be constant-slope rise after adjacent two data it Between constant difference, L_ijModel length i-th exported for the jth input set, L_ij≥I_ij+1。

Step 2. designs the dynamic matrix controller of i-th boiler, specifically:

2.1 utilize the model vector a that step 1 obtains_ijSetting up the dynamic matrix of boiler drum level, its form is as follows:

Wherein A_ijFor jth boiler input P × M rank dynamic matrix to the output of i-th boiler, a_ijK () is jth pot The stove input step response data to the output of i-th boiler, P is the optimization time domain of Dynamic array control algorithm, and M is dynamic matrix The control time domain of control algolithm, L_ij=L (1≤i≤3,1≤j≤3), M < P < L, N=3 are input and output number；

2.2 model prediction initial communication values y obtaining the i-th boiler current k moment_i,0(k)

First, controlling increment △ u is added in the k-1 moment₁(k-1),△u₂(k-1),…,△u_n(k-1) (n=3), obtains The model predication value y of i-th boiler_i,P(k-1):

$y_{i,P}(k-1)=y_{i,0}(k-1)+A_{ii,0}{\mathrm{\&Delta;u}}_i(k-1)+\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij,0}{\mathrm{\&Delta;u}}_j(k-1)$

Wherein,

y_i,P(k-1)=[y_i,1(k|k-1),y_i,1(k+1|k-1),…,y_i,1(k+L-1|k-1)]^T

y_i,0(k-1)=[y_i,0(k|k-1),y_i,0(k+1|k-1),…,y_i,0(k+L-1|k-1)]^T,

A_ii,0=[a_ii(1),a_ii(2),…,a_ii(L)]^T,A_ij,0=[a_ij(1),a_ij(2),…,a_ij(L)]^T

y_i,1(k|k-1),y_i,1(k+1|k-1),…,y_i,1(k+L-1 | k-1) represent that i-th boiler is in the k-1 moment pair respectively K, k+1 ..., the model predication value in k+L-1 moment, y_i,0(k|k-1),y_i,0(k+1|k-1),…,y_i,0(k+L-1 | k-1) represent The k-1 moment to k, k+1 ..., the initial prediction in k+L-1 moment, A_ii,0,A_ij,0It is respectively i-th boiler and jth boiler is defeated Enter the matrix that the step response data to the output of i-th boiler is set up, △ u₁(k-1),△u₂(k-1),…,△u_n(k-1) it is k- 1 moment each boiler input to water valve valve opening increment；

It is then possible to obtain model predictive error value e of k moment i-th boiler_i(k):

e_i(k)=y_i(k)-y_i,1(k|k-1)

Wherein y_iThe real output value of k i-th boiler that () expression k moment records；

Obtain k moment revised model output valve y further_i,cor(k):

y_i,cor(k)=y_i,0(k-1)+h₁*e_i(k)+h₂*e_i(k)

Wherein,

y_i,cor(k)=[y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k+L-1|k)]^T,

h₁=[1, α ..., α]^T,h₂=[0,1 ..., L-1]^T

y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k+L-1 | k) represent that i-th boiler is at k moment model respectively Correction value, h₁And h₂For the weight matrix of error compensation, α is error correction coefficient, 0 < α≤1；

Finally obtain initial communication value y of the model prediction in i-th boiler k moment_i,0(k):

y_i,0(k)=Sy_i,cor(k)

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

2.3 obtain i-th boiler at M continuous print controlling increment △ u according to step 2.1_i(k),△u_i(k+1),…,△ u_i(k+M-1) prediction output valve y under_i,PM, concrete grammar is:

$y_{i,PM}\left(k\right)=y_{i,P0}\left(k\right)+A_{ii}{\mathrm{\&Delta;u}}_{i,M}\left(k\right)+\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij}{\mathrm{\&Delta;u}}_{j,M}\left(k\right)$

Wherein,

y_i,PM(k)=[y_i,M(k+1|k),y_i,M(k+2|k),…,y_i,M(k+P|k)]^T

y_i,P0(k)=[y_i,0(k+1|k),y_i,0(k+2|k),…,y_i,0(k+P|k)]^T

△u_i,M(k)=[△ u_i(k),△u_i(k+1),…,△u_i(k+M-1)]^T

△u_j,M(k)=[△ u_j(k),△u_j(k+1),…,△u_j(k+M-1)]^T

y_i,P0K () is y_i,0The front P item of (k), y_i,0(k+1|k),y_i,0(k+2|k),…,y_i,0(k+P | k) it is that the k moment is to k+ 1, k+2 ..., the model prediction output valve in k+P moment；

2.4 performance indications J setting up boiler drum level object i-th boiler dynamic matrix controller_i(k) and reference rail Mark ω_i(k), form is as follows:

minJ_i(k)=(ω_i(k)-y_i,PM(k))^TQ_i(ω_i(k)-y_i,PM(k))+△u_i,M(k)^TR_i△u_i,M(k)

ω_i(k)=[ω_i(k+1),ω_i(k+2),…,ω_i(k+P)]^T

ω_i(k+ ε)=β^εy(k)+(1-β^ε) c (k) (ε=1,2 ..., P)

WhereinFor error weighting matrix,For controlling weighting square Battle array,WithIt is respectively Q_i,R_iIn weight coefficient, ω_iK () is the reference locus of i-th boiler, β is The softening coefficient of reference locus；

The thought of 2.5 foundation Nash optimization, is obtained receiving of i-th boiler current k moment by performance indications in step 2.4 assorted Optimal solution:

${\mathrm{\&Delta;u}}_{i,M}^{*}\left(k\right)=D_{ii}\left({\mathrm{\&omega;}}_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij}{\mathrm{\&Delta;u}}_{j,M}^{*}(k\left)\right)$

Wherein:

2.6 can be obtained by step 2.2 to 2.5 in the new round iteration optimal solution of k moment boiler i be:

${\mathrm{\&Delta;u}}_{i,M}^{l+1}\left(k\right)=D_{ii}\left({\mathrm{\&omega;}}_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij}{\mathrm{\&Delta;u}}_{j,M}^l(k\left)\right)$

Obtain the whole system optimal control law in the k moment further:

${\mathrm{\&Delta;u}}_M^{l+1}\left(k\right)=D_1\left(\mathrm{\&omega;}\right(k)-y_{P0}(k\left)\right)-D_0{\mathrm{\&Delta;u}}_M^l\left(k\right)$

Wherein:

${\mathrm{\&Delta;u}}_M^{l+1}\left(k\right)={\&lsqb;{\mathrm{\&Delta;u}}_{1,M}^{l+1}\left(k\right),{\mathrm{\&Delta;u}}_{2,M}^{l+1}\left(k\right),...,{\mathrm{\&Delta;u}}_{n,M}^{l+1}\left(k\right)\&rsqb;}^T$

${\mathrm{\&Delta;u}}_M^l\left(k\right)={\&lsqb;{\mathrm{\&Delta;u}}_{1,M}^l\left(k\right),{\mathrm{\&Delta;u}}_{2,M}^l\left(k\right),...,{\mathrm{\&Delta;u}}_{n,M}^l\left(k\right)\&rsqb;}^T$

ω (k)=[ω₁(k),ω₂(k),…,ω_n(k)]^T, y_P0(k)=[y_1,P0(k),y_2,P0(k),…,y_n,P0(k)]^T

2.7 using the Nash optimization solution first term in i-th boiler k moment as instant control law △ u_iK (), obtains boiler i's Actual water-supply valve valve opening u_i(k)=u_i(k-1)+△u_iK () acts on i-th boiler；

2.8 at subsequent time, repeats step 2.2 to 2.7 and continues to solve the instant control law △ u of i-th boiler_i(k+1), And then obtain the optimal control law △ u (k+1) of whole system, and circulate successively.

Claims (1)

1. the distributed dynamic matrix majorization method of a Nonself-regulating plant, it is characterised in that the method comprises the following steps:

Step 1. sets up corresponding dynamic matrix model vector by the real-time step response data of Nonself-regulating plant, specifically:

The large scale system of one N input N output Nonself-regulating plant, according to Distributed Predictive Control thought, is separated into N number of by 1.1 Intelligent body subsystem；

1.2 under steady state operating conditions, for input, i-th intelligent body output is carried out step response with jth intelligent body controlled quentity controlled variable Experiment, records jth (1≤j≤N) the individual input step response curve to i-th (1≤i≤N) individual output respectively；

1.3 step response curves step 1.2 obtained are filtered processing, and then fit to a smooth curve, recording light The step response data that on sliding curve, each sampling instant is corresponding, first sampling instant is T_s, adjacent two sampling instants Interval time is T_s, sampling instant order is T_s、2T_s、3T_s……；The step response data of controlled device will be in some moment t_L=I_ijT_sStart to present and determine slope rising, with the data in this momentFor starting point, data before are denoted as respectivelySet up jth input to the step response model vector a between i-th output_ij:

$a_{ij}={\&lsqb;a_1^{ij},a_2^{ij},...,a_{I_{ij}-1}^{ij},a_{I_{ij}}^{ij}+\mathrm{\&delta;},a_{I_{ij}}^{ij}+2\mathrm{\&delta;},...,a_{L_{ij}}^{ij}\&rsqb;}^T$

$a_{L_{ij}}^{ij}=a_{I_{ij}}^{ij}+(L_{ij}-I_{ij})\mathrm{\&delta;}$

Wherein T is the transposition symbol of matrix, δ be step response data be after constant-slope rises between adjacent two data Constant difference, L_ijModel length i-th exported for the jth input set, L_ij≥I_ij+1；

Step 2. designs the dynamic matrix controller of i-th intelligent body, specifically:

2.1 utilize the model vector a that step 1 obtains_ijSetting up the dynamic matrix of controlled device, its form is as follows:

Wherein A_ijFor jth intelligent body input P × M rank dynamic matrix to the output of i-th intelligent body, a_ijK () is jth input Step response data to i-th output, P is the optimization time domain of Dynamic array control algorithm, and M is Dynamic array control algorithm Control time domain, L_ij=L (1≤i≤3,1≤j≤3), M < P < L, N are input and output number；

2.2 model prediction initial communication values y obtaining the i-th intelligent body current k moment_i,0(k)

First, controlling increment △ u is added in the k-1 moment₁(k-1),△u₂(k-1),…,△u_n(k-1) i-th intelligent body, is obtained Model predication value y_i,P(k-1):

$y_{i,P}(k-1)=y_{i,0}(k-1)+A_{ii,0}{\mathrm{\&Delta;u}}_i(k-1)+\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij,0}{\mathrm{\&Delta;u}}_j(k-1)$

Wherein,

y_i,P(k-1)=[y_i,1(k|k-1),y_i,1(k+1|k-1),…,y_i,1(k+L-1|k-1)]^T

y_i,0(k-1)=[y_i,0(k|k-1),y_i,0(k+1|k-1),…,y_i,0(k+L-1|k-1)]^T,

A_ii,0=[a_ii(1),a_ii(2) ..., a_ii(L)]^T,A_ij,0=[a_ij(1),a_ij(2),…,a_ij(L)]^T

y_i,1(k|k-1),y_i,1(k+1|k-1),…,y_i,1(k+L-1 | k-1) represent respectively i-th intelligent body in the k-1 moment to k, K+1 ..., the model predication value in k+L-1 moment, y_i,0(k|k-1),y_i,0(k+1|k-1),…,y_i,0(k+L-1 | k-1) represent k-1 Moment to k, k+1 ..., the initial prediction in k+L-1 moment, A_ii,0,A_ij,0It is respectively i-th intelligent body and jth intelligent body is defeated Enter the matrix that the step response data to the output of i-th intelligent body is set up, △ u₁(k-1),△u₂(k-1),…,△u_n(k-1) it is The input controlled quentity controlled variable of k-1 moment each intelligent body；

Then, model predictive error value e of k moment i-th intelligent body is obtained_i(k):

e_i(k)=y_i(k)-y_i,1(k|k-1)

Wherein y_iThe real output value of k i-th intelligent body that () expression k moment records；

Obtain k moment revised model output valve y further_i,cor(k):

y_i,cor(k)=y_i,0(k-1)+h₁*e_i(k)+h₂*e_i(k)

Wherein,

y_i,cor(k)=[y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k+L-1|k)]^T,

h₁=[1, α ..., α]^T,h₂=[0,1 ..., L-1]^T

y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k+L-1 | k) represent i-th intelligent body repairing at k moment model respectively On the occasion of, h₁And h₂For the weight matrix of error compensation, α is error correction coefficient, 0 < α≤1；

Finally obtain initial communication value y of the model prediction in i-th intelligent body k moment_i,0(k):

y_i,0(k)=Sy_i,cor(k)

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

2.3 obtain i-th intelligent body at M continuous print controlling increment △ u according to step 2.1_i(k),△u_i(k+1),…,△u_i (k+M-1) prediction output valve y under_i,PM:

$y_{i,PM}\left(k\right)=y_{i,P0}\left(k\right)+A_{ii}{\mathrm{\&Delta;u}}_{i,M}\left(k\right)+\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij}{\mathrm{\&Delta;u}}_{j,M}\left(k\right)$

Wherein,

y_i,PM(k)=[y_i,M(k+1|k),y_i,M(k+2|k),…,y_i,M(k+P|k)]^T

y_i,P0(k)=[y_i,0(k+1|k),y_i,0(k+2|k),…,y_i,0(k+P|k)]^T

△u_i,M(k)=[△ u_i(k),△u_i(k+1),…,△u_i(k+M-1)]^T

△u_j,M(k)=[△ u_j(k),△u_j(k+1),…,△u_j(k+M-1)]^T

y_i,P0K () is y_i,0The front P item of (k), y_i,0(k+1|k),y_i,0(k+2|k),…,y_i,0(k+P | k) it is that the k moment is to k+1, k+ 2 ..., the model prediction output valve in k+P moment；

2.4 performance indications J setting up Nonself-regulating plant i-th intelligent body dynamic matrix controller_i(k) and reference locus ω_i(k), Form is as follows:

minJ_i(k)=(ω_i(k)-y_i,PM(k))^TQ_i(ω_i(k)-y_i,PM(k))+△u_i,M(k)^TR_i△u_i,M(k)

ω_i(k)=[ω_i(k+1),ω_i(k+2),…,ω_i(k+P)]^T

ω_i(k+ ε)=β^εy(k)+(1-β^ε) c (k) (ε=1,2 ..., P)

WhereinFor error weighting matrix,For controlling weighting matrix,WithIt is respectively Q_i,R_iIn weight coefficient, ω_iK () is the reference locus of i-th intelligent body, β is The softening coefficient of reference locus；

2.5 according to the thought of Nash optimization, by performance indications in step 2.4 obtain the i-th intelligent body current k moment receive assorted Excellent solution:

${\mathrm{\&Delta;u}}_{i,M}^{*}\left(k\right)=D_{ii}\left({\mathrm{\&omega;}}_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij}{\mathrm{\&Delta;u}}_{j,M}^{*}\left(k\right))$

Wherein:

2.6 obtained by step 2.2 to 2.5 in the new round iteration optimal solution of k moment intelligent body i be:

${\mathrm{\&Delta;u}}_{i,M}^{l+1}\left(k\right)=D_{ii}\left({\mathrm{\&omega;}}_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij}{\mathrm{\&Delta;u}}_{j,M}^l\left(k\right))$

Obtain the whole system optimal control law in the k moment further:

${\mathrm{\&Delta;u}}_M^{l+1}\left(k\right)=D_1\left(\mathrm{\&omega;}\right(k)-y_{P0}(k\left)\right)-D_0{\mathrm{\&Delta;u}}_M^l\left(k\right)$

Wherein:

${\mathrm{\&Delta;u}}_M^{l+1}\left(k\right)={\&lsqb;{\mathrm{\&Delta;u}}_{1,M}^{l+1}\left(k\right),{\mathrm{\&Delta;u}}_{2,M}^{l+1}\left(k\right),...,{\mathrm{\&Delta;u}}_{n,M}^{l+1}\left(k\right)\&rsqb;}^T$

${\mathrm{\&Delta;u}}_M^l\left(k\right)={\&lsqb;{\mathrm{\&Delta;u}}_{1,M}^l\left(k\right),{\mathrm{\&Delta;u}}_{2,M}^l\left(k\right),...,{\mathrm{\&Delta;u}}_{n,M}^l\left(k\right)\&rsqb;}^T$

ω (k)=[ω₁(k),ω₂(k),…,ω_n(k)]^T, y_P0(k)=[y_1,P0(k),y_2,P0(k),…,y_n,P0(k)]^T

2.7 using the Nash optimization solution first term in i-th intelligent body k moment as instant control law △ u_iK (), obtains the reality of intelligent body i Border controlled quentity controlled variable u_i(k)=u_i(k-1)+△u_iK () acts on i-th intelligent body；

2.8 at subsequent time, repeats step 2.2 to 2.7 and continues to solve the instant control law △ u of i-th intelligent body_i(k+1), enter And obtain the optimal solution △ u (k+1) of whole system, and circulate successively.