CN105955014A - Method for controlling coke furnace chamber pressure based on distributed dynamic matrix control optimization - Google Patents

Abstract

The invention discloses a method for controlling the coke furnace chamber pressure based on distributed dynamic matrix control optimization. The method comprises the steps of firstly building an input-output model of a multivariable process through acquiring step response data, transforming an online optimization implementation problem of the multivariable process into an optimization implementation problem of each small-scale subsystem, regarding each subsystem in a network environment as an intelligent body, and carrying out material, energy and information communication among the intelligent bodies through a network so as to improve the control performance of the whole system; then selecting appropriate performance indexes, solving a Nash optimization solution of each intelligent body through continuous iterative computation, thus acquiring parameters of a PID controller, then implementing a real-time control law of the moment to each intelligent body, rolling the time domain to the next moment, and repeating the above optimization process so as to complete the optimization task. The method disclosed by the invention can well process a multivariable coupling system, makes up the uncertainty of the model to a certain extent, and improves the control performance of the system.

Description

Distributed dynamic matrix control optimized coke furnace pressure control method

Technical Field

The invention belongs to the technical field of automation, and relates to a coke oven hearth pressure proportional-integral-derivative (PID) control method based on Distributed Dynamic Matrix Control (DDMC) optimization.

Background

In an actual industrial process, the PID control system is simple in structure and convenient to operate, and is widely applied to an actual process control system. However, with the development of scientific technology, modern process systems exhibit large scale, strong coupling, and numerous constraints, which present many challenges to the application of conventional PID control. DDMC is used as the natural extension and development of predictive control, and comprehensively utilizes the computer communication technology and the control theory, so that multivariable, strong-coupling and uncertain controlled objects can be well controlled, and the control performance of the system is improved. If the DDMC and PID control technology can be combined in the actual process, the system control performance can be further improved, and meanwhile, the simple form of the control structure can be ensured.

Disclosure of Invention

The invention aims to provide a DDMC (direct data multi-carrier) optimization-based coke oven hearth pressure PID control method aiming at the defect of multivariable process control of the traditional PID control. The method combines DDMC and a traditional PID control algorithm, makes up the defects of the traditional PID control, and ensures good control performance. The method firstly establishes an input/output model of the multivariable process by collecting step response data, then converts the online optimization implementation problem of the multivariable process into the optimization implementation problem of each small-scale subsystem, considers each subsystem under the network environment as an intelligent agent, and simultaneously carries out material, energy and information communication between the intelligent agents through the network so as to improve the control performance of the whole system. Then selecting proper performance indexes, calculating Nash optimal solutions of the agents through continuous iterative calculation, further obtaining PID controller parameters, implementing the instant control law of the time for each agent, rolling the time domain to the next time, and repeating the optimization process, thereby completing the optimization task of the whole system.

The technical scheme of the invention is that a coke oven hearth pressure PID control method based on DDMC optimization is established by means of data acquisition, model establishment, prediction mechanism, optimization and the like, a multivariable coupling system can be well processed by the method, the uncertainty of the model is made up to a certain extent, and the control performance of the system is improved.

The method comprises the following steps:

step 1, establishing a model of a controlled object through real-time step response data of a coke oven hearth pressure object, wherein the specific method comprises the following steps:

1.1 in steady-state operating conditions, in u_jFor input to output y_iPerforming a step response experiment, and respectively recording a step response curve of j (j is more than or equal to 1 and less than or equal to 3) th input to i (i is more than or equal to 1 and less than or equal to 3) th output;

1.2, filtering the step response curve obtained in the step 1.1, fitting the step response curve into a smooth curve, and recording step response data corresponding to each sampling moment on the smooth curve, wherein the first sampling moment is T_sThe sampling time sequence is T_s、2T_s、3T_s… …, respectively; the step response of the controlled object will be at a certain time t_L＝L_ijT_sThen tends to be smooth whenAndwhen the error of (A) and the measurement error are of the same order of magnitude, it can be considered thatApproximately equal to the steady state value of the step response. Establishing a step response model vector a between the jth input and the ith output_ij：

$a_{ij}={\[a_1^{ij},a_2^{ij},...,a_{L_{ij}}^{ij}\]}^T$

Where T is the transposed symbol of the matrix, L_ijThe time domain for the jth input versus the ith output is modeled.

Step 2, designing a PID controller of the ith intelligent agent, wherein the specific method comprises the following steps:

2.1 Using the model vector a obtained in step 1_ijEstablishing a dynamic matrix of the controlled object, wherein the form of the dynamic matrix is as follows:

wherein A is_ijInputting the P × M order dynamic matrix for the jth agent to the ith agent, a_ij(k) For j input data corresponding to i output step response, P is optimized time domain of dynamic matrix control algorithm, M is control time domain of dynamic matrix control algorithm, N is 3 number of input and output, and L is recorded for convenience_ij＝L(1≤i≤N,1≤j≤N),M<P<L；

2.2 obtaining model prediction initial response value y of ith intelligent agent at current k moment_i,0(k)

First, control increment △ u is added at time k-1₁(k-1),△u₂(k-1),…,△u_n(k-1) obtaining a model predicted value y of the ith agent_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) represents the model prediction value of the ith agent at the time k-1 to the time k, k +1, …, k + L-1, y_i,0(k|k-1),y_i,0(k+1|k-1),…,y_i,0(k + L-1| k-1) represents the initial predicted value at time k-1 versus time k, k +1, …, k + L-1, A_ii,0,A_ij,0Matrices built for the ith agent and the jth agent's step response data for the ith agent, △ u, respectively₁(k-1),△u₂(k-1),…,△u_n(k-1) is the input control quantity of each agent at the time of k-1;

then, the model prediction error value e of the ith agent at time k can be obtained_i(k)：

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

Wherein y is_i(k) Representing the actual output value of the ith intelligent agent measured at the k moment;

further obtaining a model output value y after the k moment is corrected_i,cor(k)：

y_i,cor(k)＝y_i,0(k-1)+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

y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k + L-1| k) respectively represents the corrected value of the model of the ith intelligent agent at the moment k, h is a weight matrix of error compensation, and α is an error correction coefficient;

finally, obtaining an initial response value y of model prediction at the moment k_i,0(k)：

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

Wherein S is a state transition matrix of L x L order,

2.3 calculate M consecutive control increments △ u for the ith furnace_i(k),△u_i(k+1),…,△u_iPredicted output value y at (k + M-1)_i,PMThe specific method comprises the following steps:

$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,P0(k) is y_i,0(k) The first P term, y_i,0(k+1|k),y_i,0(k+2|k),…,y_i,0(k + P | k) is a model prediction output value of the k moment to the k +1, k +2, … and k + P moment;

2.4 selecting the performance index J of the ith agent_i(k) The 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)

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

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

ω_i(k+)＝β y(k)+(1-β )c(k)(＝1,2,…,P)

whereinIn order to be a matrix of error weights,in order to control the weighting matrix,andare respectively Q_i,R_iWeight coefficient of middle, ω_i(k) The reference track of the ith agent is β, and the softening coefficient of the reference track is β;

2.5, the control quantity of the ith agent is transformed:

$\begin{array}{c}{u}_{i,M}\left(k\right)=u_{i,M}(k-1)+K_p^{i,M}\left(k\right)\left(e_{i,M}\right(k)-e_{i,M}(k-1\left)\right)+K_I^{i,M}\left(k\right)e_{i,M}\left(k\right)\\ +K_d^{i,M}\left(k\right)\left(e_{i,M}\right(k)-2e_{i,M}(k-1)+e_{i,M}(k-2\left)\right)\end{array}$

wherein,

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

e_i,M(k)＝[e_i(k),e_i(k+1),…,e_i(k+M-1)]^T

$K_p^{i,M}\left(k\right)=Blockdiag\left(K_p^i\right(k),K_p^i(k+1),...,K_p^i(k+M-1\left)\right)$

$K_I^{i,M}\left(k\right)=Blockdiag\left(K_I^i\right(k),K_I^i(k+1),...,K_I^i(k+M-1\left)\right)$

$K_d^{i,M}\left(k\right)=Blockdiag\left(K_d^i\right(k),K_d^i(k+1),...,K_d^i(k+M-1\left)\right)$

the Block diag is a diagonal matrix Block,proportional coefficient, integral coefficient and differential coefficient of PID controller of ith intelligent agent at k moment_i,M(k) The error of the reference track and the actual output of the ith intelligent agent at the moment k;

2.6 the control quantity in step 2.5 can be further converted into the following form:

u_i,M(k)＝u_i,M(k-1)+Φ_i,M(k)E_i,M(k)

E_i,M(k)＝[e_i,M(k),e_i,M(k-1),e_i,M(k-2)]^T

wherein,

e_i,M(k)＝[e_i(k),e_i(k+1),…,e_i(k+M-1)]^T

2.7 consider the form of control increments for the ith agent as follows:

△u_i,M(k)＝Φ_i,M(k)E_i,M(k)

substituting it into step 2.4 to obtain performance index J_i(k) Derivation is carried out to obtain a new round of iterative optimal solution of the agent i at the moment k

${\mathrm{\&Phi;}}_{i,M}^{l+1}\left(k\right)=\frac{{\left({\mathrm{\&omega;}}_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij}{\mathrm{\&Phi;}}_{j,M}^l(k\left)E_{j,M}\right(k\left)\right)}^T{Q}_i{A}_{ii}{E}_{i,M}\left(k\right)}{(A_{ii}^T{Q}_i{A}_{ii}+R_i)E_{i,M}{\left(k\right)}^T{E}_{i,M}\left(k\right)}$

According to the Nash optimal thought, the Nash optimal solution of the agent i can be obtained through continuous iterative computation:

${\mathrm{\&Phi;}}_{i,M}^{*}\left(k\right)=\frac{{\left({\mathrm{\&omega;}}_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij}{\mathrm{\&Phi;}}_{j,M}^{*}(k\left)E_{j,M}\right(k\left)\right)}^T{Q}_i{A}_{ii}{E}_{i,M}\left(k\right)}{(A_{ii}^T{Q}_i{A}_{ii}+R_i)E_{i,M}{\left(k\right)}^T{E}_{i,M}\left(k\right)}$

further, it is possible to obtain:

finally, the obtained parameters of the PID controller of the ith intelligent agentConstitutes a control quantity u_i,M(k) And acts on the ith agent;

2.8 at the next moment, repeating the steps 2.2 to 2.7 to continuously solve the new parameter of the ith intelligent agent PID controllerAnd are cycled in sequence.

The invention provides a DDMC optimization-based coke furnace hearth pressure PID control method. The method establishes an input/output model of the multivariable process by collecting step response data, designs a novel PID controller, ensures simple control structure form, makes up the defects of the traditional PID control, effectively reduces the scale and complexity of the multivariable process control problem, and improves the control effect.

Detailed Description

Taking the coke oven hearth pressure control as an example:

the coke oven hearth pressure control system is a typical multivariable coupling process, and the regulating means adopts the valve opening degree of a control flue damper.

Step 1, establishing a model of a controlled object through real-time step response data of a coke oven hearth pressure object, wherein the specific method comprises the following steps:

1.1 under a steady-state working condition, taking the jth hearth input as an input to carry out a step response experiment on the ith hearth output, and respectively recording a step response curve of the jth input (j is more than or equal to 1 and less than or equal to 3) to the ith output (i is more than or equal to 1 and less than or equal to 3);

1.2, filtering the step response curve obtained in the step 1.1, fitting the step response curve into a smooth curve, and recording step response data corresponding to each sampling moment on the smooth curve, wherein the first sampling moment is T_sThe sampling time sequence is T_s、2T_s、3T_s… …, respectively; the step response of the valve opening of the flue damper will be at a certain time t_L＝L_ijT_sThen tends to be smooth whenAndwhen the error of (A) and the measurement error are of the same order of magnitude, it can be considered thatApproximately equal to the steady state value of the step response. Establishing a step response model vector a between the jth input and the ith output_ij：

$a_{ij}={\&lsqb;a_1^{ij},a_2^{ij},\mathrm{...},a_{L_{ij}}^{ij}\&rsqb;}^T$

Where T is the transposed symbol of the matrix, L_ijModeling the time domain for the jth input versus the ith output;

step 2, designing a PID controller of the ith hearth, wherein the specific method comprises the following steps:

2.1 Using the model vector a obtained in step 1_ijEstablishing a dynamic matrix of coke oven hearth pressure objects, which is in the form of:

wherein A is_ijP × M order dynamic matrix, a, for jth furnace input to ith furnace_ij(k) For j input data corresponding to i output step response, P is optimized time domain of dynamic matrix control algorithm, M is control time domain of dynamic matrix control algorithm, N is 3 number of input and output, and L is recorded for convenience_ij＝L(1≤i≤N,1≤j≤N),M<P<L；

2.2 obtaining model prediction initial response value y of ith hearth at current k moment_i,0(k)

First obtain the k-1 time to add the control increment △ u_i(k-1) (i is more than or equal to 1 and less than or equal to 3), and the model predicted value y of the ith furnace_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. mu. k-1) each representsThe model predicted value y of the ith agent at the moment k-1 to the moment k, k +1, … and k + L-1_i,0(k|k-1),y_i,0(k+1|k-1),…,y_i,0(k + L-1| k-1) represents the initial predicted value at time k-1 versus time k, k +1, …, k + L-1, A_ii,0,A_ij,0Matrices built for the ith agent and the jth agent's step response data for the ith agent, △ u, respectively_i(k-1) (i is more than or equal to 1 and less than or equal to 3) is the input control quantity of each agent at the moment of k-1;

then obtaining a model prediction error value e of the ith hearth at the moment k_i(k)：

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

Wherein y is_i(k) Representing the actual output value of the ith hearth measured at the moment k;

further obtaining a model output value y corrected at the moment k_i,cor(k)：

y_i,cor(k)＝y_i,0(k-1)+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

y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k + L-1| k) respectively represents the corrected value of the model of the ith hearth at the moment k, h is a weight matrix of error compensation, and α is an error correction coefficient;

finally, obtaining an initial response value y of model prediction at the moment k_i,0(k)：

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

Wherein S is a state transition matrix of L x L order,

2.3 calculate M consecutive control increments △ u for the ith furnace_i(k),△u_i(k+1),…,△u_iPredicted output value y at (k + M-1)_i,PMThe specific method comprises the following steps:

$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

wherein y is_i,P0(k) Is y_i,0(k) The first P term, y_i,0(k+1|k),y_i,0(k+2|k),…,y_i,0(k + P | k) is a model prediction output value of the k moment to the k +1, k +2, … and k + P moment;

2.4 selecting the performance index J of the ith hearth_i(k) The 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)

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

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

ω_i(k+)＝β y(k)+(1-β )c(k)(＝1,2,…,P)

wherein ω is_i(k) Is a reference track of the ith furnace chamber,in order to be a matrix of error weights,in order to control the weighting matrix,andare respectively Q_i,R_iβ is the softening coefficient of the reference track;

2.5, the opening degree of the valve of the ith hearth is changed:

$\begin{array}{c}{u}_{i,M}\left(k\right)=u_{i,M}(k-1)+K_p^{i,M}\left(k\right)\left(e_{i,M}\right(k)-e_{i,M}(k-1\left)\right)+K_I^{i,M}\left(k\right)e_{i,M}\left(k\right)\\ +K_d^{i,M}\left(k\right)\left(e_{i,M}\right(k)-2e_{i,M}(k-1)+e_{i,M}(k-2\left)\right)\end{array}$

wherein,

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

e_i,M(k)＝[e_i(k),e_i(k+1),…,e_i(k+M-1)]^T

$K_p^{i,M}\left(k\right)=Blockdiag\left(K_p^i\right(k),K_p^i(k+1),...,K_p^i(k+M-1\left)\right)$

$K_I^{i,M}\left(k\right)=Blockdiag\left(K_I^i\right(k),K_I^i(k+1),...,K_I^i(k+M-1\left)\right)$

$K_d^{i,M}\left(k\right)=Blockdiag\left(K_d^i\right(k),K_d^i(k+1),...,K_d^i(k+M-1\left)\right)$

proportional coefficient, integral coefficient and differential coefficient of PID controller of ith furnace at k moment respectively, e_i,M(k) The error of the reference track and the actual output of the ith hearth at the moment k;

2.6 the valve opening in step 2.5 can be further converted to the following form:

u_i,M(k)＝u_i,M(k-1)+Φ_i,M(k)E_i,M(k)

E_i,M(k)＝[e_i,M(k),e_i,M(k-1),e_i,M(k-2)]^T

wherein,

e_i,M(k)＝[e_i(k),e_i(k+1),…,e_i(k+M-1)]^T

2.7 consider the form of control increments for the ith furnace as follows:

△u_i,M(k)＝Φ_i,M(k)E_i,M(k)

substituting it into step 2.4 to obtain performance index J_i(k) Derivation is carried out, and a new round of iterative optimal solution of the hearth i at the moment k can be obtained

${\mathrm{\&Phi;}}_{i,M}^{l+1}\left(k\right)=\frac{{\left({\mathrm{\&omega;}}_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij}{\mathrm{\&Phi;}}_{j,M}^l(k\left)E_{j,M}\right(k\left)\right)}^T{Q}_i{A}_{ii}{E}_{i,M}\left(k\right)}{(A_{ii}^T{Q}_i{A}_{ii}+R_i)E_{i,M}{\left(k\right)}^T{E}_{i,M}\left(k\right)}$

According to the Nash optimal idea, the Nash optimal solution of the hearth i can be obtained through continuous iterative computation:

${\mathrm{\&Phi;}}_{i,M}^{*}\left(k\right)=\frac{{\left({\mathrm{\&omega;}}_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{A}_{ij}{\mathrm{\&Phi;}}_{j,M}^{*}(k\left)E_{j,M}\right(k\left)\right)}^T{Q}_i{A}_{ii}{E}_{i,M}\left(k\right)}{(A_{ii}^T{Q}_i{A}_{ii}+R_i)E_{i,M}{\left(k\right)}^T{E}_{i,M}\left(k\right)}$

further, it is possible to obtain:

finally obtaining the parameters of the ith furnace PID controllerThe valve opening u is formed later_i,M(k) Acting on the ith hearth;

2.8 at the next moment, repeating the steps 2.2 to 2.7 to continuously solve the new parameters of the ith furnace PID controllerAnd are cycled in sequence.

Claims (1)

1. The distributed dynamic matrix control optimized coke oven hearth pressure control method is characterized by comprising the following steps;

step 1, establishing a model of a controlled object through real-time step response data of a coke oven hearth pressure object, specifically:

1.1 in steady-state operating conditions, in u_jFor input to output y_iPerforming a step response experiment, and respectively recording a step response curve of j (j is more than or equal to 1 and less than or equal to 3) th input to i (i is more than or equal to 1 and less than or equal to 3) th output;

1.2 Steps from step 1.1Filtering the response curve, fitting the response curve into a smooth curve, and recording step response data corresponding to each sampling time on the smooth curve, wherein the first sampling time is T_sThe sampling time sequence is T_s、2T_s、3T_s… …, respectively; the step response of the controlled object will be at a certain time t_L＝L_ijT_sThen tends to be smooth whenAndwhen the error of (A) and the measurement error are of the same order of magnitude, it can be considered thatApproximately equal to the steady state value of the step response; establishing a step response model vector a between the jth input and the ith output_ij：

$a_{ij}={\&lsqb;a_1^{ij},a_2^{ij},...,a_{L_{ij}}^{ij}\&rsqb;}^T$

Where T is the transposed symbol of the matrix, L_ijModeling the time domain for the jth input versus the ith output;

step 2, designing a PID controller of the ith intelligent agent, specifically:

2.1 Using the model vector a obtained in step 1_ijEstablishing a dynamic matrix of the controlled object, wherein the form of the dynamic matrix is as follows:

wherein A is_ijInputting the P × M order dynamic matrix for the jth agent to the ith agent, a_ij(k) For the jth input to ith output step response data, P is the optimized time domain of the dynamic matrix control algorithm, M is the control time domain of the dynamic matrix control algorithm, L_ij＝L(1≤i≤3,1≤j≤3),M<P<L, N is 3 and is the number of input and output;

2.2 obtaining model prediction initial response value y of ith intelligent agent at current k moment_i,0(k)

First, a control increment Δ u is added at time k-1₁(k-1),Δu₂(k-1),…,Δu_n(k-1) obtaining a model predicted value y of the ith agent_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) represents the model prediction value of the ith agent at the time k-1 to the time k, k +1, …, k + L-1, y_i,0(k|k-1),y_i,0(k+1|k-1),…,y_i,0(k + L-1| k-1) represents the initial predicted value at time k-1 versus time k, k +1, …, k + L-1, A_ii,0,A_ij,0Matrices, Deltau u, respectively established for the ith agent and the jth agent for the ith agent's step response data₁(k-1),Δu₂(k-1),…,Δu_n(k-1) is the input control quantity of each agent at the time of k-1;

then, obtaining a model prediction error value e of the ith intelligent agent at the moment k_i(k)：

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

Wherein y is_i(k) Representing the actual output value of the ith intelligent agent measured at the k moment;

further obtaining a model output value y after the k moment is corrected_i,cor(k)：

y_i,cor(k)＝y_i,0(k-1)+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

y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k + L-1| k) respectively represents the corrected value of the model of the ith intelligent agent at the moment k, h is a weight matrix of error compensation, and α is an error correction coefficient;

finally, obtaining an initial response value y of model prediction at the moment k_i,0(k)：

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

Wherein S is a state transition matrix of L x L order,

2.3 calculate M consecutive control increments Deltau for the ith furnace_i(k),Δu_i(k+1),…,Δu_iPredicted output value y at (k + M-1)_i,PMThe method specifically comprises the following steps:

$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,P0(k) is y_i,0(k) The first P term, y_i,0(k+1|k),y_i,0(k+2|k),…,y_i,0(k + P | k) is a model prediction output value of the k moment to the k +1, k +2, … and k + P moment;

2.4 selecting the performance index J of the ith agent_i(k) The 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)

Δu_i,M(k)＝[Δu_i(k),Δu_i(k+1),…,Δu_i(k+M-1)]^T

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

ω_i(k+)＝β y(k)+(1-β )c(k)(＝1,2,…,P)

whereinIn order to be a matrix of error weights,in order to control the weighting matrix,andare respectively Q_i,R_iWeight coefficient of middle, ω_i(k) The reference track of the ith agent is β, and the softening coefficient of the reference track is β;

2.5, the control quantity of the ith agent is transformed:

$\begin{array}{c}{u}_{i,M}\left(k\right)=u_{i,M}\left(k-1\right)+K_p^{i,M}\left(k\right)\left(e_{i,M}\left(k\right)-e_{i,M}\left(k-1\right)\right)+K_I^{i,M}\left(k\right)e_{i,M}\left(k\right)\\ +K_d^{i,M}\left(k\right)\left(e_{i,M}\left(k\right)-2e_{i,M}\left(k-1\right)+e_{i,M}\left(k-2\right)\right)\end{array}$

wherein,

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

e_i,M(k)＝[e_i(k),e_i(k+1),…,e_i(k+M-1)]^T

$K_p^{i,M}\left(k\right)=Blockdiag\left(K_p^i\right(k),K_p^i(k+1),...,K_p^i(k+M-1\left)\right)$

$K_I^{i,M}\left(k\right)=Blockdiag\left(K_I^i\right(k),K_I^i(k+1),...,K_I^i(k+M-1\left)\right)$

$K_d^{i,M}\left(k\right)=Blockdiag\left(K_d^i\right(k),K_d^i(k+1),...,K_d^i(k+M-1\left)\right)$

the Block diag is a diagonal matrix Block,proportional coefficient, integral coefficient and differential coefficient of PID controller of ith intelligent agent at k moment_i,M(k) The error of the reference track and the actual output of the ith intelligent agent at the moment k;

2.6 further converts the control quantity in step 2.5 into the following form:

u_i,M(k)＝u_i,M(k-1)+Φ_i,M(k)E_i,M(k)

E_i,M(k)＝[e_i,M(k),e_i,M(k-1),e_i,M(k-2)]^T

wherein,

e_i,M(k)＝[e_i(k),e_i(k+1),…,e_i(k+M-1)]^T

2.7 consider the form of control increments for the ith agent as follows:

Δu_i,M(k)＝Φ_i,M(k)E_i,M(k)

substituting it into step 2.4 to obtain performance index J_i(k) Derivation is carried out to obtain a new round of iterative optimal solution of the agent i at the moment k

${\mathrm{\&Phi;}}_{i,M}^{l+1}\left(k\right)=\frac{{\left({\mathrm{\&omega;}}_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\&NotEqual;i}{\overset{n}{\mathrm{\&Sigma;}}}{A}_{ij}{\mathrm{\&Phi;}}_{j,M}^l(k\left)E_{j,M}\right(k\left)\right)}^T{Q}_i{A}_{ii}{E}_{i,M}\left(k\right)}{(A_{ii}^T{Q}_i{A}_{ii}+R_i)E_{i,M}{\left(k\right)}^T{E}_{i,M}\left(k\right)}$

According to the Nash optimal thought, the Nash optimal solution of the agent i can be obtained through continuous iterative computation:

${\mathrm{\&Phi;}}_{i,M}^{*}\left(k\right)=\frac{{\left({\mathrm{\&omega;}}_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\&NotEqual;i}{\overset{n}{\mathrm{\&Sigma;}}}{A}_{ij}{\mathrm{\&Phi;}}_{j,M}^{*}(k\left)E_{j,M}\right(k\left)\right)}^T{Q}_i{A}_{ii}{E}_{i,M}\left(k\right)}{(A_{ii}^T{Q}_i{A}_{ii}+R_i)E_{i,M}{\left(k\right)}^T{E}_{i,M}\left(k\right)}$

further obtaining:

finally, the obtained parameters of the PID controller of the ith intelligent agentConstitutes a control quantity u_i,M(k) And acts on the ith agent;

2.8 at the next moment, repeating the steps 2.2 to 2.7 to continuously solve the new parameter of the ith intelligent agent PID controllerAnd are cycled in sequence.