CN105955014A

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

Info

Publication number: CN105955014A
Application number: CN201610308520.8A
Authority: CN
Inventors: 张日东; 汪大卫; 张俊锋
Original assignee: Hangzhou Dianzi University
Current assignee: Hangzhou Dianzi University
Priority date: 2016-05-11
Filing date: 2016-05-11
Publication date: 2016-09-21

Abstract

本发明公开了一种分布式动态矩阵控制优化的焦炭炉炉膛压力控制方法。本发明首先通过采集阶跃响应数据建立多变量过程的输入输出模型，再将多变量过程的在线优化实施问题转化成各个小规模子系统的优化实施问题，把网络环境下的每个子系统看作为一个智能体，同时各智能体之间通过网络进行物质、能量及信息通信，以提高整个系统的控制性能。然后再选取合适的性能指标通过不断迭代计算求取各智能体的纳什最优解，进而得到PID控制器参数，再对每个智能体实施该时刻的即时控制律，并将时域滚动到下一时刻，重复上述优化过程，从而完成优化任务。本发明能很好的处理多变量耦合系统，在一定程度上弥补了模型的不确定性，改善了系统的控制性能。The invention discloses a method for controlling pressure in a coke oven chamber optimized by distributed dynamic matrix control. The present invention first establishes the input-output model of the multivariable process by collecting step response data, and then converts the online optimization implementation problem of the multivariable process into the optimization implementation problem of each small-scale subsystem, and regards each subsystem in the network environment as a An intelligent body, and at the same time, all intelligent bodies communicate with each other through the network in terms of matter, energy and information to improve the control performance of the entire system. Then select the appropriate performance index to obtain the Nash optimal solution of each agent through continuous iterative calculation, and then obtain the PID controller parameters, and then implement the real-time control law for each agent at this moment, and scroll the time domain to the next For a moment, the above optimization process is repeated to complete the optimization task. The invention can well deal with the multi-variable coupling system, compensates the uncertainty of the model to a certain extent, and improves the control performance of the system.

Description

Coke oven furnace pressure control method optimized by distributed dynamic matrix control

技术领域technical field

本发明属于自动化技术领域，涉及一种基于分布式动态矩阵控制(DDMC)优化的焦炭炉炉膛压力比例积分微分(PID)控制方法。The invention belongs to the technical field of automation, and relates to a proportional-integral-derivative (PID) control method for coke oven furnace pressure optimization based on distributed dynamic matrix control (DDMC).

背景技术Background technique

在实际工业过程中，由于PID控制结构简单，操作方便，已经广泛地应用于实际过程控制系统中。但是，随着科学技术的发展，现代过程系统呈现出规模庞大，耦合强烈，以及约束众多的特点，给传统PID控制的应用带来了诸多挑战。DDMC作为预测控制自然的延伸和发展，综合利用计算机通信技术和控制理论，能很好的控制存在多变量、强耦合、不确定的被控对象，提高了系统的控制性能。如果能够在实际过程中将DDMC和PID控制技术相结合，系统控制性能将得到进一步改善，同时又能保证控制结构的形式简单。In the actual industrial process, due to the simple structure and convenient operation of PID control, it has been widely used in the actual process control system. However, with the development of science and technology, the modern process system presents the characteristics of large scale, strong coupling and numerous constraints, which brings many challenges to the application of traditional PID control. As a natural extension and development of predictive control, DDMC comprehensively utilizes computer communication technology and control theory, which can well control the controlled objects with multi-variable, strong coupling and uncertainty, and improve the control performance of the system. If the DDMC and PID control technology can be combined in the actual process, the system control performance will be further improved, and at the same time, the form of the control structure can be kept simple.

发明内容Contents of the invention

本发明目的是针对传统PID控制在多变量过程控制的不足之处，提出了一种基于DDMC优化的焦炭炉炉膛压力PID控制方法。该方法结合DDMC和传统PID控制算法，弥补了传统PID控制的不足，又保证了良好的控制性能。本发明方法首先通过采集阶跃响应数据建立多变量过程的输入输出模型，再将多变量过程的在线优化实施问题转化成各个小规模子系统的优化实施问题，把网络环境下的每个子系统看作为一个智能体，同时各智能体之间通过网络进行物质、能量及信息通信，以提高整个系统的控制性能。然后再选取合适的性能指标通过不断迭代计算求取各智能体的纳什最优解，进而得到PID控制器参数，再对每个智能体实施该时刻的即时控制律，并将时域滚动到下一时刻，重复上述优化过程，从而完成整个系统的优化任务。The purpose of the present invention is to propose a PID control method for coke oven furnace pressure based on DDMC optimization aiming at the shortcomings of traditional PID control in multivariable process control. This method combines DDMC and traditional PID control algorithm to make up for the shortcomings of traditional PID control and ensure good control performance. The method of the present invention first establishes the input-output model of the multivariable process by collecting step response data, and then converts the online optimization implementation problem of the multivariable process into the optimization implementation problem of each small-scale subsystem, and sees each subsystem under the network environment As an intelligent body, at the same time, each intelligent body communicates matter, energy and information through the network to improve the control performance of the entire system. Then select the appropriate performance index to obtain the Nash optimal solution of each agent through continuous iterative calculation, and then obtain the PID controller parameters, and then implement the real-time control law for each agent at this moment, and scroll the time domain to the next For a moment, the above optimization process is repeated to complete the optimization task of the entire system.

本发明的技术方案是通过数据采集、模型建立、预测机理、优化等手段，确立了一种基于DDMC优化的焦炭炉炉膛压力PID控制方法，利用该方法能很好的处理多变量耦合系统，在一定程度上弥补了模型的不确定性，改善了系统的控制性能。The technical solution of the present invention is to establish a coke oven furnace pressure PID control method based on DDMC optimization by means of data collection, model establishment, prediction mechanism, optimization, etc., using this method to handle multi-variable coupling systems well. To a certain extent, it makes up for the uncertainty of the model and improves the control performance of the system.

本发明方法的步骤包括：The steps of the inventive method comprise:

步骤1.通过焦炭炉炉膛压力对象的实时阶跃响应数据建立被控对象的模型，具体方法是：Step 1. Establish the model of the controlled object through the real-time step response data of the coke oven furnace pressure object. The specific method is:

1.1在稳态工况下，以u_j为输入对输出y_i进行阶跃响应实验，分别记录第j(1≤j≤3)个输入对第i(1≤i≤3)个输出的阶跃响应曲线；1.1 Under steady-state conditions, take u _j as the input to conduct step response experiments on the output y _i , and record the order of the j-th (1≤j≤3) input to the i-th (1≤i≤3) output jump response curve;

1.2将步骤1.1得到的阶跃响应曲线进行滤波处理，然后拟合成一条光滑曲线，记录光滑曲线上每个采样时刻对应的阶跃响应数据，第一个采样时刻为T_s，采样时刻顺序为T_s、2T_s、3T_s……；被控对象的阶跃响应将在某一个时刻t_L＝L_ijT_s后趋于平稳，当与的误差和测量误差有相同的数量级时，即可认为近似等于阶跃响应的稳态值。建立第j个输入对第i个输出之间的阶跃响应模型向量a_ij：1.2 Filter the step response curve obtained in step 1.1, and then fit it into a smooth curve, record the step response data corresponding to each sampling time on the smooth curve, the first sampling time is T _s , and the order of sampling time is T _s , 2T _s , 3T _s ...; the step response of the controlled object will tend to be stable after a certain moment t _L =L _ij T _s , when and When the error and measurement error have the same order of magnitude, it can be considered It is approximately equal to the steady-state value of the step response. Establish the step response model vector a _ij between the jth input and the ith output:

${a a}_{i i j j} = = {[[{a a}_{11}^{i i j j},, {a a}_{22}^{i i j j},, ... ...,, {a a}_{{L L}_{i i j j}}^{i i j j}]]}^{T T}$

其中T为矩阵的转置符号,L_ij为第j个输入对第i个输出的建模时域。Where T is the transpose symbol of the matrix, L _ij is the modeling time domain of the jth input to the ith output.

步骤2.设计第i个智能体的PID控制器，具体方法是：Step 2. Design the PID controller of the i-th agent, the specific method is:

2.1利用步骤1获得的模型向量a_ij建立被控对象的动态矩阵，其形式如下：2.1 Use the model vector a _ij obtained in step 1 to establish the dynamic matrix of the controlled object, the form of which is as follows:

其中A_ij为第j个智能体输入对第i个智能体的P×M阶动态矩阵，a_ij(k)为第j个输入对第i个输出阶跃响应的数据，P为动态矩阵控制算法的优化时域，M为动态矩阵控制算法的控制时域，N＝3为输入输出个数，为方便起见记L_ij＝L(1≤i≤N,1≤j≤N),M<P<L；Among them, A _ij is the P×M order dynamic matrix of the j-th agent input to the i-th agent, a _ij (k) is the data of the step response of the j-th input to the i-th output, and P is the dynamic matrix control The optimization time domain of the algorithm, M is the control time domain of the dynamic matrix control algorithm, N=3 is the number of input and output, for convenience, record L _ij =L(1≤i≤N,1≤j≤N),M<P<L;

2.2获取第i个智能体当前k时刻的模型预测初始响应值y_i,0(k)2.2 Obtain the model prediction initial response value y _i,0 (k) of the i-th agent at the current moment k

首先，在k-1时刻加入控制增量△u₁(k-1),△u₂(k-1),…,△u_n(k-1)，得到第i个智能体的模型预测值y_i,P(k-1)：First, add control increments △u ₁ (k-1), △u ₂ (k-1),...,△u _n (k-1) at time k-1 to get the model prediction value of the i-th agent y _i,P (k-1):

${y the y}_{i i,, P P} ((k k - - 11)) = = {y the y}_{i i,, 00} ((k k - - 11)) + + {A A}_{i i i i,, 00} {Δu Δ u}_{i i} ((k k - - 11)) + + {Σ Σ}_{j j = = 11,, j j &NotEqual; &NotEqual; i i}^{n no} {A A}_{i i j j,, 00} {Δu Δu}_{j j} ((k k - - 11))$

其中，in,

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 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)分别表示第i个智能体在k-1时刻对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)表示k-1时刻对k,k+1,…,k+L-1时刻的初始预测值，A_ii,0,A_ij,0分别为第i个智能体和第j个智能体对第i个智能体的阶跃响应数据建立的矩阵，△u₁(k-1),△u₂(k-1),…,△u_n(k-1)为k-1时刻各智能体的输入控制量；y _i,1 (k|k-1), y _i,1 (k+1|k-1),…,y _i,1 (k+L-1|k-1) represent the i-th agent respectively At time k-1, the model prediction value at 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 forecast value at k-1 time to k,k+1,...,k+L-1 time, A _ii,0 ,A _{ij , 0} is the matrix established by the step response data of the i-th agent and the j-th agent to the i-th agent, △u ₁ (k-1), △u ₂ (k-1),…, △u _n (k-1) is the input control amount of each agent at time k-1;

然后，可以得到k时刻第i个智能体的模型预测误差值e_i(k)：Then, the model prediction error value e _i (k) of the i-th agent at time k can be obtained:

e_i(k)＝y_i(k)-y_i,1(k|k-1)e _i (k)=y _i (k)-y _i,1 (k|k-1)

其中y_i(k)表示k时刻测得的第i个智能体的实际输出值；Where y _i (k) represents the actual output value of the i-th agent measured at time k;

进一步得到k时刻修正后的模型输出值y_i,cor(k)：Further obtain the corrected model output value y _i,cor (k) at time k:

y_i,cor(k)＝y_i,0(k-1)+h*e_i(k)y _i,cor (k)=y _i,0 (k-1)+h*e _i (k)

其中，in,

y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k+L-1|k)分别表示第i个智能体在k时刻模型的修正值，h为误差补偿的权矩阵，α为误差校正系数；y _{i, cor} (k|k), y _{i, cor} (k+1|k),..., y _{i, cor} (k+L-1|k) respectively represent the correction of the model of the i-th agent at time k value, h is the weight matrix of error compensation, and α is the error correction coefficient;

最后得到k时刻的模型预测的初始响应值y_i,0(k)：Finally, the initial response value y _i,0 (k) predicted by the model at time k is obtained:

y_i,0(k)＝Sy_i,cor(k)y _i,0 (k)=Sy _i,cor (k)

其中，S为L×L阶的状态转移矩阵，Among them, S is the state transition matrix of L×L order,

2.3计算第i个炉膛在M个连续的控制增量△u_i(k),△u_i(k+1),…,△u_i(k+M-1)下的预测输出值y_i,PM，具体方法是：2.3 Calculate the predicted output value y _{i of} the i-th furnace under M continuous control increments △u _i (k), △u _i (k+1),..., △u _i (k+M-1), _PM , the specific method is:

${y the y}_{i i,, P P M m} ((k k)) = = {y the y}_{i i,, P P 00} ((k k)) + + {A A}_{i i i i} {Δu Δ u}_{i i,, M m} ((k k)) + + {Σ Σ}_{j j = = 11,, j j &NotEqual; &NotEqual; i i}^{n no} {A A}_{i i j j} {Δu Δu}_{j j,, M m} ((k k))$

其中，in,

△u_i,M(k)＝[△u_i(k),△u_i(k+1),…,△u_i(k+M-1)]^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 △u _j,M (k)＝[△u _j (k),△u _j (k+1),…,△u _j (k+M-1)] ^T

2.4选取第i个智能体的性能指标J_i(k)，形式如下：2.4 Select the performance index J _i (k) of the i-th agent, 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)minJ _i (k)＝(ω _i (k)-y _i,PM (k)) ^T Q _i (ω _i (k)-y _i,PM (k))+△u _i,M (k) ^T R _i △ u _i,M (k)

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

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

其中为误差加权矩阵，为控制加权矩阵，和分别为Q_i,R_i中的权重系数，ω_i(k)为第i个智能体的参考轨迹，β为参考轨迹的柔化系数；in is the error weighting matrix, For the control weight matrix, and are the weight coefficients in Q _i and R _i respectively, ω _i (k) is the reference trajectory of the i-th agent, and β is the softening coefficient of the reference trajectory;

2.5将第i个智能体的控制量进行变换：2.5 Transform the control amount of the i-th agent:

$\begin{matrix} {u u}_{i i,, M m} ((k k)) = = {u u}_{i i,, M m} ((k k - - 11)) + + {K K}_{p p}^{i i,, M m} ((k k)) (({e e}_{i i,, M m} ((k k)) - - {e e}_{i i,, M m} ((k k - - 11)))) + + {K K}_{I I}^{i i,, M m} ((k k)) {e e}_{i i,, M m} ((k k)) \\ + + {K K}_{d d}^{i i,, M m} ((k k)) (({e e}_{i i,, M m} ((k k)) - - 22 {e e}_{i i,, M m} ((k k - - 11)) + + {e e}_{i i,, M m} ((k k - - 22)))) \end{matrix}$

其中，in,

u_i,M(k)＝[u_i(k),u_i(k+1),…,u_i(k+M-1)]^T 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 e _i,M (k)=[e _i (k),e _i (k+1),…,e _i (k+M-1)] ^T

${K K}_{p p}^{i i,, M m} ((k k)) = = B B l l o o c c k k d d i i a a g g (({K K}_{p p}^{i i} ((k k)),, {K K}_{p p}^{i i} ((k k + + 11)),, ... ...,, {K K}_{p p}^{i i} ((k k + + M m - - 11))))$

${K K}_{I I}^{i i,, M m} ((k k)) = = B B l l o o c c k k d d i i a a g g (({K K}_{I I}^{i i} ((k k)),, {K K}_{I I}^{i i} ((k k + + 11)),, ... ...,, {K K}_{I I}^{i i} ((k k + + M m - - 11))))$

${K K}_{d d}^{i i,, M m} ((k k)) = = B B l l o o c c k k d d i i a a g g (({K K}_{d d}^{i i} ((k k)),, {K K}_{d d}^{i i} ((k k + + 11)),, ... ...,, {K K}_{d d}^{i i} ((k k + + M m - - 11))))$

Block diag为对角矩阵块，分别为k时刻第i个智能体的PID控制器的比例系数、积分系数和微分系数，e_i,M(k)为k时刻第i个智能体的参考轨迹和实际输出的误差；Block diag is a diagonal matrix block, are the proportional coefficient, integral coefficient and differential coefficient of the PID controller of the i-th agent at time k, respectively, and e _i,M (k) is the error between the reference trajectory and the actual output of the i-th agent at time k;

2.6进一步可将步骤2.5中控制量转换为如下形式：2.6 Further, the control quantity in step 2.5 can be converted into the following form:

u_i,M(k)＝u_i,M(k-1)+Φ_i,M(k)E_i,M(k)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 E _i,M (k)=[e _i,M (k),e _i,M (k-1),e _i,M (k-2)] ^T

其中，in,

2.7考虑第i个智能体的控制增量形式如下：2.7 Consider the control increment form of the i-th agent as follows:

△u_i,M(k)＝Φ_i,M(k)E_i,M(k)△u _i,M (k)=Φ _i,M (k)E _i,M (k)

将其代入步骤2.4，对性能指标J_i(k)求导，可得到在k时刻智能体i的新一轮迭代最优解为Substituting it into step 2.4, and deriving the performance index J _i (k), the optimal solution for a new round of iteration of agent i at time k can be obtained as

${Φ Φ}_{i i,, M m}^{l l + + 11} ((k k)) = = \frac{{(({ω ω}_{i i} ((k k)) - - {y the y}_{i i,, P P 00} ((k k)) - - {Σ Σ}_{j j = = 11,, j j &NotEqual; &NotEqual; i i}^{n no} {A A}_{i i j j} {Φ Φ}_{j j,, M m}^{l l} ((k k)) {E E.}_{j j,, M m} ((k k))))}^{T T} {Q Q}_{i i} {A A}_{i i i i} {E E.}_{i i,, M m} ((k k))}{(({A A}_{i i i i}^{T T} {Q Q}_{i i} {A A}_{i i i i} + + {R R}_{i i})) {E E.}_{i i,, M m} {((k k))}^{T T} {E E.}_{i i,, M m} ((k k))}$

依据纳什最优的思想，通过不断迭代计算可以得到智能体i的纳什最优解：According to the idea of Nash optimality, the Nash optimal solution of agent i can be obtained through continuous iterative calculation:

${Φ Φ}_{i i,, M m}^{* *} ((k k)) = = \frac{{(({ω ω}_{i i} ((k k)) - - {y the y}_{i i,, P P 00} ((k k)) - - {Σ Σ}_{j j = = 11,, j j &NotEqual; &NotEqual; i i}^{n no} {A A}_{i i j j} {Φ Φ}_{j j,, M m}^{* *} ((k k)) {E E.}_{j j,, M m} ((k k))))}^{T T} {Q Q}_{i i} {A A}_{i i i i} {E E.}_{i i,, M m} ((k k))}{(({A A}_{i i i i}^{T T} {Q Q}_{i i} {A A}_{i i i i} + + {R R}_{i i})) {E E.}_{i i,, M m} {((k k))}^{T T} {E E.}_{i i,, M m} ((k k))}$

进一步可以得到：Further can get:

最后将得到的第i个智能体的PID控制器的参数构成控制量u_i,M(k)，并作用于第i个智能体；Finally, the parameters of the PID controller of the i-th agent will be obtained Constitute the control quantity u _i,M (k), and act on the i-th agent;

2.8在下一时刻，重复步骤2.2到2.7继续求解第i个智能体PID控制器新的参数并依次循环。2.8 At the next moment, repeat steps 2.2 to 2.7 to continue to solve the new parameters of the i-th agent PID controller And cycle in turn.

本发明提出了一种基于DDMC优化的焦炭炉炉膛压力PID控制方法。该方法通过采集阶跃响应数据建立多变量过程的输入输出模型，设计了一种新型的PID控制器，在保证控制结构形式简单的同时，弥补了传统PID控制的不足，有效地降低了多变量过程控制问题的规模和复杂性，改善了控制效果。The invention proposes a coke oven furnace pressure PID control method based on DDMC optimization. This method establishes the input-output model of the multivariable process by collecting step response data, and designs a new type of PID controller. The scale and complexity of the process control problem improves control effectiveness.

具体实施方式detailed description

以焦炭炉炉膛压力控制为例：Take the coke oven furnace pressure control as an example:

焦炭炉炉膛压力控制系统是一个典型的多变量耦合过程，调节手段采用控制烟道挡板的阀门开度。The coke oven furnace pressure control system is a typical multi-variable coupling process, and the adjustment means is to control the valve opening of the flue baffle.

步骤1.通过焦炭炉炉膛压力对象的实时阶跃响应数据建立被控对象的模型，具体方法是：Step 1. Establish the model of the controlled object through the real-time step response data of the coke oven furnace pressure object, the specific method is:

1.1在稳态工况下，以第j个炉膛输入为输入对第i个炉膛输出进行阶跃响应实验，分别记录第j(1≤j≤3)个输入对第i(1≤i≤3)个输出的阶跃响应曲线；1.1 Under steady-state conditions, take the j-th furnace input as the input to conduct a step response experiment on the i-th furnace output, and record the j-th (1≤j≤3) input to the i-th (1≤i≤3) input ) output step response curves;

1.2将步骤1.1得到的阶跃响应曲线进行滤波处理，然后拟合成一条光滑曲线，记录光滑曲线上每个采样时刻对应的阶跃响应数据，第一个采样时刻为T_s，采样时刻顺序为T_s、2T_s、3T_s ……；烟道挡板的阀门开度的阶跃响应将在某一个时刻t_L＝L_ijT_s后趋于平稳，当与的误差和测量误差有相同的数量级时，即可认为近似等于阶跃响应的稳态值。建立第j个输入对第i个输出之间的阶跃响应模型向量a_ij：1.2 Filter the step response curve obtained in step 1.1, and then fit it into a smooth curve, record the step response data corresponding to each sampling time on the smooth curve, the first sampling time is T _s , and the order of sampling time is T _s , 2T _s , 3T _s ... ; the step response of the valve opening of the flue damper will tend to be stable after a certain moment t _L =L _ij T _s , when and When the error and measurement error have the same order of magnitude, it can be considered It is approximately equal to the steady-state value of the step response. Establish the step response model vector a _ij between the jth input and the ith output:

其中T为矩阵的转置符号,L_ij为第j个输入对第i个输出的建模时域；Where T is the transpose symbol of the matrix, L _ij is the modeling time domain of the jth input to the ith output;

步骤2.设计第i个炉膛的PID控制器，具体方法是：Step 2. Design the PID controller of the i-th furnace, the specific method is:

2.1利用步骤1获得的模型向量a_ij建立焦炭炉炉膛压力对象的动态矩阵，其形式如下：2.1 Use the model vector a _ij obtained in step 1 to establish the dynamic matrix of the coke oven furnace pressure object, and its form is as follows:

其中A_ij为第j个炉膛输入对第i个炉膛的P×M阶动态矩阵，a_ij(k)为第j个输入对第i个输出阶跃响应的数据，P为动态矩阵控制算法的优化时域，M为动态矩阵控制算法的控制时域，N＝3为输入输出个数，为方便起见记L_ij＝L(1≤i≤N,1≤j≤N),M<P<L；Among them, A _ij is the P×M order dynamic matrix of the input of the j-th furnace to the i-th furnace, a _ij (k) is the data of the step response of the j-th input to the i-th output, and P is the dynamic matrix control algorithm Optimizing the time domain, M is the control time domain of the dynamic matrix control algorithm, N=3 is the number of input and output, for convenience, record L _ij =L(1≤i≤N,1≤j≤N), M<P<L;

2.2获取第i个炉膛当前k时刻的模型预测初始响应值y_i,0(k)2.2 Obtain the model predicted initial response value y _i,0 (k) of the i-th furnace at the current moment k

先得到k-1时刻加入控制增量△u_i(k-1)(1≤i≤3),后第i个炉膛的模型预测值y_i,P(k-1)：First get the control increment △u _i (k-1) (1≤i≤3) at time k-1, and then the model prediction value y _i,P (k-1) of the i-th furnace:

${y the y}_{i i,, P P} ((k k - - 11)) = = {y the y}_{i i,, 00} ((k k - - 11)) + + {A A}_{i i i i,, 00} {Δu Δu}_{i i} ((k k - - 11)) + + {Σ Σ}_{j j = = 11,, j j &NotEqual; &NotEqual; i i}^{n no} {A A}_{i i j j,, 00} {Δu Δu}_{j j} ((k k - - 11))$

其中，in,

y_i,1(k|k-1),y_i,1(k+1|k-1),…,y_i,1(k+L-1|k-1)分别表示第i个智能体在k-1时刻对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)表示k-1时刻对k,k+1,…,k+L-1时刻的初始预测值，A_ii,0,A_ij,0分别为第i个智能体和第j个智能体对第i个智能体的阶跃响应数据建立的矩阵，△u_i(k-1)(1≤i≤3)为k-1时刻各智能体的输入控制量；y _i,1 (k|k-1), y _i,1 (k+1|k-1),…,y _i,1 (k+L-1|k-1) represent the i-th agent respectively At time k-1, the model prediction value at 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 forecast value at k-1 time to k,k+1,...,k+L-1 time, A _ii,0 ,A _{ij , 0} is the matrix established by the step response data of the i-th agent and the j-th agent to the i-th agent, △u _i (k-1)(1≤i≤3) is the k-1 moment The input control amount of each agent;

接着得到k时刻第i个炉膛的模型预测误差值e_i(k)：Then get the model prediction error value e _i (k) of the i-th furnace at time k:

e_i(k)＝y_i(k)-y_i,1(k|k-1)e _i (k)=y _i (k)-y _i,1 (k|k-1)

其中y_i(k)表示k时刻测得的第i个炉膛的实际输出值；Among them, y _i (k) represents the actual output value of the i-th furnace measured at time k;

进一步得到k时刻修正的模型输出值y_i,cor(k)：Further obtain the corrected model output value y _i,cor (k) at time k:

y_i,cor(k)＝y_i,0(k-1)+h*e_i(k)y _i,cor (k)=y _i,0 (k-1)+h*e _i (k)

其中，in,

y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k+L-1|k)分别表示第i个炉膛在k时刻模型的修正值，h为误差补偿的权矩阵，α为误差校正系数；y _{i, cor} (k|k), y _{i, cor} (k+1|k),..., y _{i, cor} (k+L-1|k) respectively represent the correction value of the i-th furnace at time k , h is the weight matrix of error compensation, α is the error correction coefficient;

y_i,0(k)＝Sy_i,cor(k)y _i,0 (k)=Sy _i,cor (k)

${y the y}_{i i,, P P M m} ((k k)) = = {y the y}_{i i,, P P 00} ((k k)) + + {A A}_{i i i i} {Δu Δu}_{i i,, M m} ((k k)) + + {Σ Σ}_{j j = = 11,, j j &NotEqual; &NotEqual; i i}^{n no} {A A}_{i i j j} {Δu Δu}_{j j,, M m} ((k k))$

其中，in,

其中y_i,P0(k)是y_i,0(k)的前P项，y_i,0(k+1|k),y_i,0(k+2|k),…,y_i,0(k+P|k)为k时刻对k+1,k+2,…,k+P时刻的模型预测输出值；Where y _i,P0 (k) is the first P item of y i, ₀ (k), y _i,0 (k+1|k),y _i,0 (k+2|k),…,y _{i, 0} (k+P|k) is the model prediction output value at time k for k+1, k+2,...,k+P;

2.4选取第i个炉膛的性能指标J_i(k)，形式如下：2.4 Select the performance index J _i (k) of the i-th furnace, the form is as follows:

其中ω_i(k)为第i个炉膛的参考轨迹，为误差加权矩阵，为控制加权矩阵，和分别为Q_i,R_i中的权重系数，β为参考轨迹的柔化系数；where ω _i (k) is the reference trajectory of the i-th furnace, is the error weighting matrix, For the control weight matrix, and are the weight coefficients in Q _i and R _i respectively, and β is the softening coefficient of the reference track;

2.5将第i个炉膛的阀门开度进行变换：2.5 Transform the valve opening of the i-th furnace:

其中，in,

分别为k时刻第i个炉膛的PID控制器的比例系数、积分系数和微分系数，e_i,M(k)为k时刻第i个炉膛的参考轨迹和实际输出的误差； are the proportional coefficient, integral coefficient and differential coefficient of the PID controller of the i-th furnace at time k, respectively, e _i,M (k) is the error between the reference trajectory and the actual output of the i-th furnace at k time;

2.6进一步可将步骤2.5中阀门开度转换为如下形式：2.6 Further, the valve opening in step 2.5 can be converted into the following form:

其中，in,

2.7考虑第i个炉膛的控制增量形式如下：2.7 Consider the control increment form of the i-th furnace as follows:

△u_i,M(k)＝Φ_i,M(k)E_i,M(k)△u _i,M (k)=Φ _i,M (k)E _i,M (k)

将其代入步骤2.4，对性能指标J_i(k)求导，可得到在k时刻炉膛i的新一轮迭代最优解为Substituting it into step 2.4, and deriving the performance index J _i (k), the optimal solution for a new round of iteration of furnace i at time k can be obtained as

依据纳什最优的思想，通过不断迭代计算可以得到炉膛i的纳什最优解：According to the idea of Nash optimality, the Nash optimal solution of furnace i can be obtained through continuous iterative calculation:

进一步可以得到：Further can get:

最后得到第i个炉膛PID控制器的参数以后构成阀门开度u_i,M(k)作用于第i个炉膛；Finally get the parameters of the i-th furnace PID controller Afterwards, the valve opening u _i,M (k) acts on the i-th furnace;

2.8在下一时刻，重复步骤2.2到2.7继续求解第i个炉膛PID控制器新的参数并依次循环。2.8 At the next moment, repeat steps 2.2 to 2.7 to continue to solve the new parameters of the i-th furnace PID controller And cycle in turn.

Claims

1. The coke oven hearth pressure control method of distributed dynamic matrix control optimization, it is characterized in that the method comprises the following steps;

Step 1. Establish the model of the controlled object through the real-time step response data of the coke oven furnace pressure object, specifically:

1.1 Under steady-state conditions, take u _j as the input to conduct step response experiments on the output y _i , and record the order of the j-th (1≤j≤3) input to the i-th (1≤i≤3) output jump response curve;

1.2 Filter the step response curve obtained in step 1.1, and then fit it into a smooth curve, record the step response data corresponding to each sampling time on the smooth curve, the first sampling time is T _s , and the order of sampling time is T _s , 2T _s , 3T _s ...; the step response of the controlled object will tend to be stable after a certain moment t _L =L _ij T _s , when and When the error and measurement error have the same order of magnitude, it can be considered Approximately equal to the steady-state value of the step response; establish the step response model vector a _ij between the jth input and the ith output:

{a a}_{i i j j} = = {[[{a a}_{11}^{i i j j},, {a a}_{22}^{i i j j},, ... ...,, {a a}_{{L L}_{i i j j}}^{i i j j}]]}^{T T}

Where T is the transpose symbol of the matrix, L _ij is the modeling time domain of the jth input to the ith output;

Step 2. Design the PID controller of the i-th agent, specifically:

2.1 Use the model vector a _ij obtained in step 1 to establish the dynamic matrix of the controlled object, the form of which is as follows:

Among them, A _ij is the P×M order dynamic matrix of the j-th agent input to the i-th agent, a _ij (k) is the data of the step response of the j-th input to the i-th output, and P is the dynamic matrix control The optimization time domain of the 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=3 is the number of input and output;

2.2 Obtain the model prediction initial response value y _i,0 (k) of the i-th agent at the current moment k

First, add control increments Δu ₁ (k-1), Δu ₂ (k-1),...,Δu _n (k-1) at time k-1 to obtain the model prediction value y _{i of the i-th agent, P} (k-1):

{y the y}_{i i,, P P} ((k k - - 11)) = = {y the y}_{i i,, 00} ((k k - - 11)) + + {A A}_{i i i i,, 00} {Δu Δu}_{i i} ((k k - - 11)) + + {Σ Σ}_{j j = = 11,, j j &NotEqual; &NotEqual; i i}^{n no} {A A}_{i i j j,, 00} {Δu Δ u}_{j j} ((k k - - 11))

in,

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 the i-th agent respectively At time k-1, the model prediction value at 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 forecast value at k-1 time to k,k+1,...,k+L-1 time, A _ii,0 ,A _{ij , 0} is the matrix established by the step response data of the i-th agent and the j-th agent to the i-th agent, Δu ₁ (k-1), Δu ₂ (k-1),…,Δu _n (k-1) is the input control amount of each agent at k-1 time;

Then, get the model prediction error value e _i (k) of the i-th agent at time k:

e _i (k)=y _i (k)-y _i,1 (k|k-1)

Where y _i (k) represents the actual output value of the i-th agent measured at time k;

Further obtain the corrected model output value y _i,cor (k) at time k:

y _i,cor (k)=y _i,0 (k-1)+h*e _i (k)

in,

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 represent the correction of the model of the i-th agent at time k value, h is the weight matrix of error compensation, and α is the error correction coefficient;

Finally, the initial response value y _i,0 (k) predicted by the model at time k is obtained:

y _i,0 (k)=Sy _i,cor (k)

Among them, S is the state transition matrix of L×L order,

2.3 Calculate the predicted output value y _{i, PM} of the i-th furnace under M continuous control increments Δu _i (k), Δu _i (k+1),..., Δu _i (k+M-1), specifically yes:

{y the y}_{i i,, P P M m} ((k k)) = = {y the y}_{i i,, P P 00} ((k k)) + + {A A}_{i i i i} {Δu Δu}_{i i,, M m} ((k k)) + + {Σ Σ}_{j j = = 11,, j j &NotEqual; &NotEqual; i i}^{n no} {A A}_{i i j j} {Δu Δ u}_{j j,, M m} ((k k))

in,

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 the first P item of y i, ₀ (k), y _i,0 (k+1|k),y _i,0 (k+2|k),…,y _i,0 (k+P|k) is the model prediction output value at time k for k+1, k+2,...,k+P;

2.4 Select the performance index J _i (k) of the i-th agent, the form is as follows:

minJ _i (k)＝(ω _i (k)-y _i,PM (k)) ^T Q _i (ω _i (k)-y _i,PM (k))+Δu _i,M (k) ^T R _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)

in is the error weighting matrix, For the control weight matrix, and are the weight coefficients in Q _i and R _i respectively, ω _i (k) is the reference trajectory of the i-th agent, and β is the softening coefficient of the reference trajectory;

2.5 Transform the control amount of the i-th agent:

\begin{matrix} {u u}_{i i,, M m} ((k k)) = = {u u}_{i i,, M m} ((k k - - 11)) + + {K K}_{p p}^{i i,, M m} ((k k)) (({e e}_{i i,, M m} ((k k)) - - {e e}_{i i,, M m} ((k k - - 11)))) + + {K K}_{I I}^{i i,, M m} ((k k)) {e e}_{i i,, M m} ((k k)) \\ + + {K K}_{d d}^{i i,, M m} ((k k)) (({e e}_{i i,, M m} ((k k)) - - 22 {e e}_{i i,, M m} ((k k - - 11)) + + {e e}_{i i,, M m} ((k k - - 22)))) \end{matrix}

in,

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 K}_{p p}^{i i,, M m} ((k k)) = = B B l l o o c c k k d d i i a a g g (({K K}_{p p}^{i i} ((k k)),, {K K}_{p p}^{i i} ((k k + + 11)),, ... ...,, {K K}_{p p}^{i i} ((k k + + M m - - 11))))

{K K}_{I I}^{i i,, M m} ((k k)) = = B B l l o o c c k k d d i i a a g g (({K K}_{I I}^{i i} ((k k)),, {K K}_{I I}^{i i} ((k k + + 11)),, ... ...,, {K K}_{I I}^{i i} ((k k + + M m - - 11))))

{K K}_{d d}^{i i,, M m} ((k k)) = = B B l l o o c c k k d d i i a a g g (({K K}_{d d}^{i i} ((k k)),, {K K}_{d d}^{i i} ((k k + + 11)),, ... ...,, {K K}_{d d}^{i i} ((k k + + M m - - 11))))

Block diag is a diagonal matrix block, are the proportional coefficient, integral coefficient and differential coefficient of the PID controller of the i-th agent at time k, respectively, and e _i,M (k) is the error between the reference trajectory and the actual output of the i-th agent at time k;

2.6 Further convert 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

in,

e _i,M (k)=[e _i (k),e _i (k+1),…,e _i (k+M-1)] ^T

2.7 Consider the control increment form of the i-th agent as follows:

Δu _i,M (k)=Φ _i,M (k)E _i,M (k)

Substituting it into step 2.4, deriving the performance index J _i (k), and obtaining the optimal solution of the new iteration of agent i at time k is

{Φ Φ}_{i i,, M m}^{l l + + 11} ((k k)) = = \frac{{(({ω ω}_{i i} ((k k)) - - {y the y}_{i i,, P P 00} ((k k)) - - {Σ Σ}_{j j = = 11,, j j &NotEqual; &NotEqual; i i}^{n no} {A A}_{i i j j} {Φ Φ}_{j j,, M m}^{l l} ((k k)) {E E.}_{j j,, M m} ((k k))))}^{T T} {Q Q}_{i i} {A A}_{i i i i} {E E.}_{i i,, M m} ((k k))}{(({A A}_{i i i i}^{T T} {Q Q}_{i i} {A A}_{i i i i} + + {R R}_{i i})) {E E.}_{i i,, M m} {((k k))}^{T T} {E E.}_{i i,, M m} ((k k))}

According to the idea of Nash optimality, the Nash optimal solution of agent i can be obtained through continuous iterative calculation:

{Φ Φ}_{i i,, M m}^{* *} ((k k)) = = \frac{{(({ω ω}_{i i} ((k k)) - - {y the y}_{i i,, P P 00} ((k k)) - - {Σ Σ}_{j j = = 11,, j j &NotEqual; &NotEqual; i i}^{n no} {A A}_{i i j j} {Φ Φ}_{j j,, M m}^{* *} ((k k)) {E E.}_{j j,, M m} ((k k))))}^{T T} {Q Q}_{i i} {A A}_{i i i i} {E E.}_{i i,, M m} ((k k))}{(({A A}_{i i i i}^{T T} {Q Q}_{i i} {A A}_{i i i i} + + {R R}_{i i})) {E E.}_{i i,, M m} {((k k))}^{T T} {E E.}_{i i,, M m} ((k k))}

Further get:

Finally, the parameters of the PID controller of the i-th agent will be obtained Constitute the control quantity u _i,M (k), and act on the i-th agent;

2.8 At the next moment, repeat steps 2.2 to 2.7 to continue to solve the new parameters of the i-th agent PID controller And cycle in turn.