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

Distributed dynamic matrix control method for non-self-balancing object

Technical Field

The invention belongs to the technical field of automation, and relates to a Distributed Dynamic Matrix Control (DDMC) method of a non-self-balancing object.

Background

In the actual process, a large number of complex and high-dimensional large-scale systems widely exist, the performance, the processing speed and the like of a computer are often required to be high by adopting centralized integral solution, and the requirements are contrary to the economical efficiency which must be considered in an actual industrial system. The Distributed Dynamic Matrix Control (DDMC) is used as a main branch of the distributed predictive control (DMPC), the computer communication technology and the control theory are comprehensively utilized, the online solving problem of a complex large-scale system is dispersed into each subsystem to be distributed and realized, the scale and the complexity of the problem are effectively reduced, multivariable, strong coupling and uncertain controlled objects can be well controlled, and the control performance of the system is improved. In actual industrial processes, however, there are many multivariable processes that contain non-self-balancing objects, such as a portion of the storage tank, boiler drum level, rectifier level, and so forth. Due to the fact that a transfer function of a non-self-balanced object contains a typical integral link, the response of a controlled object under a constant value step tends to be infinite, and therefore the traditional DDMC algorithm cannot be directly applied. If the traditional DDMC method can be improved in an actual process, the defects of the traditional DDMC method in multivariable process control containing non-self-balancing objects can be effectively overcome, and the DMPC is further extended and developed in the actual application.

Disclosure of Invention

The invention aims to provide a DDMC method of a non-self-balancing object aiming at the defects of the traditional DDMC method in multivariable process control containing the non-self-balancing object.

The method comprises the steps of firstly establishing a matrix model vector of a multivariable process containing non-self-balancing objects by collecting step response data, excavating basic object characteristics, then converting an online optimization implementation problem of the multivariable process into an optimization implementation problem of each small-scale subsystem, considering each subsystem under a network environment as an intelligent agent by combining theories and ideas in a plurality of intelligent agents, and carrying out material, energy and information communication among the intelligent agents through a network. Then, a method for improving the transfer matrix aiming at the non-self-balance object is combined with a new error correction method, a proper performance index is selected, Nash optimal solutions of all the agents are obtained through continuous iteration based on the Nash optimization idea, parameters of the dynamic matrix controller of each agent are further obtained, an instant control law of the moment is implemented on each agent, the time domain is rolled to the next moment, the optimization process is repeated, and therefore the optimization task of the whole system is completed.

The technical scheme of the invention is that a distributed dynamic matrix control method of a non-self-balancing object is established by means of data acquisition, model establishment, prediction mechanism, optimization and the like, and the method can effectively make up the defects of the traditional DDMC method in multivariable process control containing the non-self-balancing object on the premise of ensuring higher control precision and stability, and meet the requirements of the actual industrial process.

The method comprises the following steps:

step 1, establishing a corresponding dynamic matrix model vector through real-time step response data of a non-self-balancing object, wherein the specific method comprises the following steps:

1.1, dispersing a large-scale system with N input and N output non-self-balancing objects into N intelligent agent subsystems according to a distributed predictive control idea;

1.2 under a steady-state working condition, taking the jth intelligent agent control quantity as an input to carry out a step response experiment on the ith intelligent agent output quantity, and respectively recording step response curves of the jth input (j is more than or equal to 1 and is less than or equal to N) to the ith output (i is more than or equal to 1 and is less than or equal to N);

1.3, filtering the step response curve obtained in the step 1.2, 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 interval time between two adjacent sampling time is T_sThe sampling time sequence is T_s、2T_s、3T_s… …, respectively; the step response data of the controlled object will be at a certain time t_L＝I_ijT_sStarts to show a constant slope rise, and data at the momentAs a starting point, the previous data are respectively notedEstablishing a step response model vector a between the jth input and the ith 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 transposed sign of the matrix, L is the constant difference between two adjacent data after the step response data rises with a constant slope_ijFor a set model length of jth input to ith output, L_ij≥I_ij+1。

Step 2, designing a dynamic matrix 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_ijP × M order dynamic matrix, a, for jth agent input versus ith agent output_ij(k) Step response data of the jth input to the ith output, P is an optimized time domain of the dynamic matrix control algorithm, M is a control time domain of the dynamic matrix control algorithm, and L_ij＝L(1≤i≤3,1≤j≤3),M<P<L and N are input and output numbers;

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) denotes the initial time of k-1 versus time of k, k +1, …, k + L-1Predicted value, A_ii,0,A_ij,0Matrices built for the ith agent and jth agent input versus the ith agent output step response data, △ 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)+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) represents the correction value of the model of the ith agent at time k, h₁And h₂For error compensation, weight matrix α is error correction coefficient, 0<α≤1；

Finally, obtaining the initial response value y of model prediction of the ith intelligent agent at the moment k_i,0(k)：

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

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

2.3 obtaining M consecutive control increments △ u for the ith agent according to step 2.1_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 establishing Performance index J of No self-balance object ith agent dynamic matrix controller_i(k) And a reference track omega_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)

ω_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 according to the Nash optimal idea, obtaining the Nash optimal solution of the ith intelligent agent at the current k moment according to the performance indexes in the step 2.4:

${\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 from steps 2.2 to 2.5, a new iteration of the optimal solution for agent i at time k can be obtained as follows:

${\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)$

further obtaining the optimal control law of the whole system at the moment k:

${\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 taking the Nash optimal solution initial term of the ith agent at the moment k as the instant control law △ u_i(k) Obtaining the actual control quantity u of the intelligent agent i_i(k)＝u_i(k-1)+△u_i(k) Acting on the ith agent;

2.8 at the next moment, repeating steps 2.2 to 2.7 to continuously solve the instant control law △ u of the ith agent_i(k +1), and then the optimal solution △ u (k +1) of the whole system is obtained, and the steps are circulated in sequence.

The invention provides a DDMC method of a non-self-balancing object. The method combines a method for improving the transfer matrix aiming at the non-self-balancing object with a new error correction method on the basis of the traditional DDMC method, effectively makes up the defects of the traditional DDMC method in the multivariable process control containing the non-self-balancing object on the premise of ensuring higher control precision and stability, and meets the requirements of the actual industrial process.

Detailed Description

Taking boiler drum water level control as an example:

the boiler drum water level control system is a typical multivariable non-self-balancing object with an integral link, and the regulating means adopts the control of the opening degree of a water supply valve.

Step 1, establishing a corresponding dynamic matrix model vector through real-time step response data of a boiler drum water level object, wherein the specific method comprises the following steps:

1.1, dispersing a large-scale system of a 3-input and 3-output boiler drum water level object into 3 subsystems according to a distributed predictive control idea;

1.2 under a steady-state working condition, taking the opening degree of a jth boiler feed water valve as an input to carry out a step response experiment on the water level of an ith boiler drum, and respectively recording step response curves of the jth input (j is more than or equal to 1 and less than or equal to 3) to the ith (i is more than or equal to 1 and less than or equal to 3) output;

1.3, filtering the step response curve obtained in the step 1.2, 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 interval time between two adjacent sampling time is T_sThe sampling time sequence is T_s、2T_s、3T_s… …, respectively; the step response data of the boiler drum water level will be at a certain time t_L＝I_ijT_sStarts to show a constant slope rise, and data at the momentAs a starting point, the previous data are respectively notedEstablishing a step response model vector a between the jth boiler input and the ith boiler output_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 transposed sign of the matrix, L is the constant difference between two adjacent data after the step response data rises with a constant slope_ijFor a set model length of jth input to ith output, L_ij≥I_ij+1。

Step 2, designing a dynamic matrix controller of the ith boiler, which specifically comprises the following steps:

2.1 Using the model vector a obtained in step 1_ijEstablishing a dynamic matrix of boiler drum water levels, wherein the form of the dynamic matrix is as follows:

wherein A is_ijP × M order dynamic matrix for jth boiler input to ith boiler output, a_ij(k) Inputting step response data output by the ith boiler for the jth boiler, P being the optimized time domain of the dynamic matrix control algorithm, M being 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 boiler 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) (n is 3), and a model predicted value y of the ith boiler is obtained_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 boiler 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 i-th boiler and the j-th boiler input to the i-th boiler output step response data, △ u, respectively₁(k-1),△u₂(k-1),…,△u_n(k-1) the valve opening increment of the input feed water valve of each boiler at the time of k-1;

then, a model prediction error value e of the ith boiler at the moment 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 boiler measured at the moment k;

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)+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) represents the correction value of the model at the moment k for the ith boiler, h₁And h₂Weight matrix for error compensationα is the error correction coefficient, 0<α≤1；

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

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

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

2.3 obtaining M consecutive control increments △ u for the ith boiler according to step 2.1_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 establishing a performance index J of the ith boiler dynamic matrix controller of a boiler drum water level object_i(k) And a reference track omega_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)

ω_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 boiler is shown, and β is the softening coefficient of the reference track;

2.5 according to Nash optimal thought, obtaining Nash optimal solution of the ith boiler at the current k moment according to the performance indexes in the step 2.4:

${\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 from steps 2.2 to 2.5, a new iteration of the optimal solution for boiler i at time k can be obtained as:

${\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)$

further obtaining the optimal control law of the whole system at the moment k:

${\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 taking the Nash optimal solution first term of the ith boiler at the time k as an instant control law △ u_i(k) Obtaining the actual opening u of the feed valve of the boiler i_i(k)＝u_i(k-1)+△u_i(k) Acting on the ith boiler;

2.8 at the next moment, repeating steps 2.2 to 2.7 to continuously solve the immediate control law △ u of the ith boiler_i(k +1), and further obtaining an optimal control law △ u (k +1) of the whole system, and sequentially circulating.

Claims (1)

1. A method for distributed dynamic matrix control of a non-self-balancing object, the method comprising the steps of:

step 1, establishing a corresponding dynamic matrix model vector through real-time step response data of a non-self-balancing object, specifically:

1.1, dispersing a large-scale system with N input and N output non-self-balancing objects into N intelligent agent subsystems according to a distributed predictive control idea;

1.2 under a steady-state working condition, taking the jth intelligent agent control quantity as an input to carry out a step response experiment on the ith intelligent agent output quantity, and respectively recording step response curves of the jth input (j is more than or equal to 1 and is less than or equal to N) to the ith output (i is more than or equal to 1 and is less than or equal to N);

1.3, filtering the step response curve obtained in the step 1.2, 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 interval time between two adjacent sampling time is T_sThe sampling time sequence is T_s、2T_s、3T_s… …, respectively; the step response data of the controlled object will be at a certain time t_L＝I_ijT_sStarts to show a constant slope rise, and data at the momentAs a starting point, the previous data are respectively notedEstablishing 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_{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 transposed sign of the matrix, L is the constant difference between two adjacent data after the step response data rises with a constant slope_ijFor a set model length of jth input to ith output, L_ij≥I_ij+1；

Step 2, designing a dynamic matrix controller of the ith intelligent agent, which specifically 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_ijP × M order dynamic matrix, a, for jth agent input versus ith agent output_ij(k) Step response data of the jth input to the ith output, P is an optimized time domain of the dynamic matrix control algorithm, M is a control time domain of the dynamic matrix control algorithm, and L_ij＝L(1≤i≤3,1≤j≤3),M<P<L and N are input and output numbers;

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 jth agent input versus the ith agent output step response data, △ u, respectively₁(k-1),△u₂(k-1),…,△u_n(k-1) is the time k-1Input control quantity of each agent;

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)+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) represents the correction value of the model of the ith agent at time k, h₁And h₂For error compensation, weight matrix α is error correction coefficient, 0<α≤1；

Finally, obtaining the initial response value y of model prediction of the ith intelligent agent at the moment k_i,0(k)：

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

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

2.3 obtaining M consecutive control increments △ u for the ith agent according to step 2.1_i(k),△u_i(k+1),…,△u_iPredicted output value y at (k + M-1)_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,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 time k vs. time k +1, k +2, …, k + PA model prediction output value;

2.4 establishing Performance index J of No self-balance object ith agent dynamic matrix controller_i(k) And a reference track omega_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)

ω_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 according to the Nash optimal idea, obtaining the Nash optimal solution of the ith intelligent agent at the current k moment according to the performance indexes in the step 2.4:

${\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 from steps 2.2 to 2.5, a new iteration of the optimal solution for agent i at time k is:

${\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))$

further obtaining the optimal control law of the whole system at the moment k:

${\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 taking the Nash optimal solution initial term of the ith agent at the moment k as the instant control law △ u_i(k) Obtaining the actual control quantity u of the intelligent agent i_i(k)＝u_i(k-1)+△u_i(k) Acting on the ith agent;

2.8 at the next moment, repeating steps 2.2 to 2.7 to continuously solve the instant control law △ u of the ith agent_i(k +1), and then the optimal solution △ u (k +1) of the whole system is obtained, and the steps are circulated in sequence.