Multidimensional Interpolation Decoupling Strategy for CD Basis Weight of Papermaking Process

: With a focus on the multivariable coupling characteristics in a cross-directional (CD) basis weight control system, we study the coupling characteristics of a CD control system and decoupling control, and we propose a novel multivariable interpolation decoupling control strategy and a real-time decomposition algorithm in this paper. Based on a model of the CD basis weight proﬁle, a system non-square interaction matrix of high-dimensional data is analyzed by experimental studies and numerical simulation. Along the diagonal of the interaction matrix, the matrix block method is adopted to reduce the system dimension and convert it into a square system. A multivariable control system with high dimensionality is divided into several subsystems. For the high-dimensional Toeplitz symmetric subsystem with small-scale coupling characteristics, an interpolation decoupling algorithm is proposed. Then, a decoupling compensator with the structure of a symmetric Toeplitz matrix was obtained. Compared with the conventional diagonal decoupling matrix, the branch number of the new decoupling network is reduced from 2408 to 186, which realizes the fast decoupling of multivariable systems. The results were even better when we used a double size interaction matrix obtained by interpolation between actual values. By designing the diagonalized controller for the new decoupled system, a decouped CD control system for corrugated paper with a basis weight of 133 g / m 2 was implemented in an actual project in a paper mill. 64-dimensional a compensator with the structure of a symmetric Toeplitz was Based on the multi-actuator parallel a decoupled multivariable system is calculation this is on the data collected from the actual An interpolation decoupling controller is designed to transform an MIMO system into multiple systems. Experimental results show that the method has an obvious decoupling control e ﬀ ect, and this control strategy can provide help for the actual project. In addition, the multidimensional interpolation decoupling control strategy and a real-time decomposition algorithm can be applied to various processes whose input and output numbers are di ﬀ erent.


Introduction
Paper-making is by far the most complex and complicated of all the pulp and paper processes. An important quality characteristic in the production of paper sheets is the range of functions in the basis weight, which should be kept as small as possible. The quality of the manufactured paper is assessed in two dimensions: one is the maintenance of the average sheet property profile along the paper sheet as it moves while it is being manufactured, which is referred to as machine directional (MD) control, and the other is the maintenance of the sheet profile across the width of the paper machine, which is referred to as cross-directional (CD) control. In the process of modern high-speed paper machines, it is necessary to monitor and control the performance of paper continuously to ensure that the quality of paper meets the requirements of MD and CD [1]. The CD basis weight profile is a very important quality index of paper in the papermaking process. Control of the CD basis weight on paper machines has long been known as an especially difficult problem. This is due to several difficulties including the high dimensionality of the cross-directional system, the high lateral space interaction between actuators, the uncertainty of the model, and the limited control authority of the actuator [2].
The primary reason for the difficulty of CD control is the complexity of the profile response to a slice adjustment [3]. The modern hydraulic headbox is equipped with dilution water supply equipment for basis weight profile control. The adjustment of the CD profile by dilution water is superior to the traditional slice adjustment [4], and the pursuit of better paper quality has placed new demands on CD control systems. In a CD profile control system, the number of sheet horizontal measuring points ranges from 200 to 2000 [5], and the number of actuators is up to 300. A large number of horizontal measuring points and actuators will inevitably cause coupling problems between adjacent dilution water valves. In addition, because of the impact characteristics of the fluid flow, it will affect several adjacent areas on both sides. The final result is that when a dilution valve is activated, it will affect the CD profile value of the adjacent measurement points [6]. This leads to the control system becoming a high-dimensional, multivariable, and strongly coupled hysteretic system that is difficult to control.
For the coupling problem, the classical decoupling control includes inverse matrix decoupling and a dynamic complementary decoupling method [7]. The inverse matrix decoupling method depends on the accuracy of the mathematical model of the controlled object. The dynamic complementary decoupling may be used in the design of non-singular as well as a singular process control system. In addition, adaptive decoupling control provides a means of pattern recognition converting a linguistic control strategy based on expert knowledge into an automatic control strategy, which realizes the decoupling control of a linear time-varying MIMO (Multiple-Input Multiple-Output) system. However, it requires a large amount of calculation.
In view of the characteristics of strong multivariable coupling in a CD basis weight control system, it may be difficult for a conventional decoupling control strategy to meet the actual needs of a CD basis weight control system. Thus, it is necessary to study the effective decoupling or compensation control strategy.
To address these issues, this paper employs both matrix analysis and multiprocessor decomposition design techniques to address the actuator coupling. Based on the matrix block method and the interpolation decoupling strategy, a general method of transforming the control problem of a high-dimensional large-scale system into a single loop group control problem is proposed.
For accurate simulation, the coupling calculation is based on the data collected from the actual paper mill. The edge effect is fully considered in the coupling matrix. Based on the block matrix method, a multivariable system with high dimensionality is divided into several subsystems. For the high-dimensional Toeplitz symmetric subsystem with small-scale coupling characteristics, an interpolation decoupling controller with the structure of a symmetric Toeplitz matrix is designed to transform an MIMO system into multiple SISO (Single-Input Single-Output) systems. In additional, the one-to-one correspondence between the measurement points and the actuators is solved by a new constraint on control variables. Finally, a multiprocessor system is designed and applied to a paper mill. The multidimensional interpolation decoupling control strategy and a real-time decomposition algorithm can be applied to various processes whose numbers of input and output are different.
The rest of this paper is organized as follows. The materials of the technological process for basis weight control in the paper machine are described in Section 2. The CD control system modeling and the numerical simulation of coupling characteristics are described in Section 3. In Section 4, an interpolation decoupling strategy is designed for the Toeplitz symmetric interaction matrix, and the theoretical results show and prove each step of the algorithm. The simulation and application are presented in Section 5, and finally, we summarize the results.

Materials
The control process is carried out in a hydraulic headbox with diluted water. The operation of the hydraulic headbox is shown in Figure 1. The box body is equipped with glass holes, and the dilution water was colored with red ink so that the flow of fluid from the point of addition to the lip Symmetry 2020, 12, 149 3 of 22 of the headbox can be observed. The dilution water regulating device is shown in Figure 2, which includes a dilution water cone pipe and multiple dilution water valves with automatic control units. One interface of the dilution water valve is connected with the dilution water cone pipe, and the other interface is provided with a water conveying hose that is connected to the pulping branch pipe. The amount of white water injected into the pulp is controlled by adjusting the opening of the dilution valve, and the consistency of the pulp can thereby be adjusted, thus realizing the local fine-tuning of the cross-directional basis weight of the paper.

Materials
The control process is carried out in a hydraulic headbox with diluted water. The operation of the hydraulic headbox is shown in Figure 1. The box body is equipped with glass holes, and the dilution water was colored with red ink so that the flow of fluid from the point of addition to the lip of the headbox can be observed. The dilution water regulating device is shown in Figure 2, which includes a dilution water cone pipe and multiple dilution water valves with automatic control units. One interface of the dilution water valve is connected with the dilution water cone pipe, and the other interface is provided with a water conveying hose that is connected to the pulping branch pipe. The amount of white water injected into the pulp is controlled by adjusting the opening of the dilution valve, and the consistency of the pulp can thereby be adjusted, thus realizing the local fine-tuning of the cross-directional basis weight of the paper.

Measuring Method of the Process Data
The process of CD control for a Fourdrinier machine equipped with a hydraulic headbox is shown in Figure 3. A scanning frame, located between the calender and the winder of the paper machine, supports the basis weight on-line detecting device that scans the paper sheet that moves through the opening of the frame. The system continuously obtains the instantaneous sample value of the basis weight along the machine direction and its CD profiles. These measurements can only be made on paper that is dry enough at the end of the line. This causes a long time delay between the actuators on the hydraulic headbox and measurements from the scanner at the end of the paper machine.

Materials
The control process is carried out in a hydraulic headbox with diluted water. The operation of the hydraulic headbox is shown in Figure 1. The box body is equipped with glass holes, and the dilution water was colored with red ink so that the flow of fluid from the point of addition to the lip of the headbox can be observed. The dilution water regulating device is shown in Figure 2, which includes a dilution water cone pipe and multiple dilution water valves with automatic control units. One interface of the dilution water valve is connected with the dilution water cone pipe, and the other interface is provided with a water conveying hose that is connected to the pulping branch pipe. The amount of white water injected into the pulp is controlled by adjusting the opening of the dilution valve, and the consistency of the pulp can thereby be adjusted, thus realizing the local fine-tuning of the cross-directional basis weight of the paper.

Measuring Method of the Process Data
The process of CD control for a Fourdrinier machine equipped with a hydraulic headbox is shown in Figure 3. A scanning frame, located between the calender and the winder of the paper machine, supports the basis weight on-line detecting device that scans the paper sheet that moves through the opening of the frame. The system continuously obtains the instantaneous sample value of the basis weight along the machine direction and its CD profiles. These measurements can only be made on paper that is dry enough at the end of the line. This causes a long time delay between the actuators on the hydraulic headbox and measurements from the scanner at the end of the paper machine.

Measuring Method of the Process Data
The process of CD control for a Fourdrinier machine equipped with a hydraulic headbox is shown in Figure 3. A scanning frame, located between the calender and the winder of the paper machine, supports the basis weight on-line detecting device that scans the paper sheet that moves through the opening of the frame. The system continuously obtains the instantaneous sample value of the basis weight along the machine direction and its CD profiles. These measurements can only be made on paper that is dry enough at the end of the line. This causes a long time delay between the actuators on the hydraulic headbox and measurements from the scanner at the end of the paper machine.
The dilution actuator set-points are determined by the weight profile obtained from the scanning weight sensor. Control actions are based on the profile estimates, which are calculated from the measurements provided by the scanning sensor. The scanner and sensor interface collects analog sensor and scanner data from the basis weight measurement system and builds the basis weight profile to suit the paper machine and actuator requirements. The dilution actuator set-points are determined by the weight profile obtained from the scanning weight sensor. Control actions are based on the profile estimates, which are calculated from the measurements provided by the scanning sensor. The scanner and sensor interface collects analog sensor and scanner data from the basis weight measurement system and builds the basis weight profile to suit the paper machine and actuator requirements.
The sheet typically moves at a velocity on the order of 10 m/s, while the sensors travels back and forth across the paper web at speed from 0 to 1 m/s. This change is usually measured by installing a gauge on a frame so that the gauge can be scanned back and forth throughout the process. When the sheet passes under the gauge, the measured track follows the zig-zag path, as shown in Figure 4, where the angle 'α ' may be only approximately 1°. The horizontal parts in the path correspond to the scanner remaining at the paper edge for a few seconds before changing the direction. Moreover, the gauge stops at one side of the sheet (for unpredictable time intervals) for recalibrating itself. From this sparse data, the whole profile can be reconstructed. The experimental process was carried out in a paper mill in Wugong city, Shaanxi Province of China in September 2019. The CD basis weight on-line measurement is applied in the paper machine 12 (PM 12), the width of the PM 12 is 5400 mm, and the design speed is 650 m/min. PM 12 produces corrugated paper. In this paper mill, the data of paper basis weight is obtained by a QCS (Quality Control System) through scanning sensors. The data exchange between the OPC (OLE (Object Linking and Embedding) for Process Control) toolbox and the configuration software of WinCC is carried out by MATLAB. The control algorithm is based on the data obtained by MATLAB. The results of the calculations are fed back to WinCC. A S7-300 PLC is adopted for the controller implementation. The interface is driven by an industry standard PLC and provides the programming flexibility to adapt to most measurement systems. The transmission of profiles to the operator station is done through industrial Ethernet communication. The PLC (Programmable Logic Controller) and supporting electronics form a compact unit for ease of installation and maintenance.  The sheet typically moves at a velocity on the order of 10 m/s, while the sensors travels back and forth across the paper web at speed from 0 to 1 m/s. This change is usually measured by installing a gauge on a frame so that the gauge can be scanned back and forth throughout the process. When the sheet passes under the gauge, the measured track follows the zig-zag path, as shown in Figure 4, where the angle 'α' may be only approximately 1 • . The horizontal parts in the path correspond to the scanner remaining at the paper edge for a few seconds before changing the direction. Moreover, the gauge stops at one side of the sheet (for unpredictable time intervals) for recalibrating itself. From this sparse data, the whole profile can be reconstructed. The dilution actuator set-points are determined by the weight profile obtained from the scanning weight sensor. Control actions are based on the profile estimates, which are calculated from the measurements provided by the scanning sensor. The scanner and sensor interface collects analog sensor and scanner data from the basis weight measurement system and builds the basis weight profile to suit the paper machine and actuator requirements.

Modeling and Block Decomposition of CD Control System
The sheet typically moves at a velocity on the order of 10 m/s, while the sensors travels back and forth across the paper web at speed from 0 to 1 m/s. This change is usually measured by installing a gauge on a frame so that the gauge can be scanned back and forth throughout the process. When the sheet passes under the gauge, the measured track follows the zig-zag path, as shown in Figure 4, where the angle 'α ' may be only approximately 1°. The horizontal parts in the path correspond to the scanner remaining at the paper edge for a few seconds before changing the direction. Moreover, the gauge stops at one side of the sheet (for unpredictable time intervals) for recalibrating itself. From this sparse data, the whole profile can be reconstructed. The experimental process was carried out in a paper mill in Wugong city, Shaanxi Province of China in September 2019. The CD basis weight on-line measurement is applied in the paper machine 12 (PM 12), the width of the PM 12 is 5400 mm, and the design speed is 650 m/min. PM 12 produces corrugated paper. In this paper mill, the data of paper basis weight is obtained by a QCS (Quality Control System) through scanning sensors. The data exchange between the OPC (OLE (Object Linking and Embedding) for Process Control) toolbox and the configuration software of WinCC is carried out by MATLAB. The control algorithm is based on the data obtained by MATLAB. The results of the calculations are fed back to WinCC. A S7-300 PLC is adopted for the controller implementation. The interface is driven by an industry standard PLC and provides the programming flexibility to adapt to most measurement systems. The transmission of profiles to the operator station is done through industrial Ethernet communication. The PLC (Programmable Logic Controller) and supporting electronics form a compact unit for ease of installation and maintenance.  The experimental process was carried out in a paper mill in Wugong city, Shaanxi Province of China in September 2019. The CD basis weight on-line measurement is applied in the paper machine 12 (PM 12), the width of the PM 12 is 5400 mm, and the design speed is 650 m/min. PM 12 produces corrugated paper. In this paper mill, the data of paper basis weight is obtained by a QCS (Quality Control System) through scanning sensors. The data exchange between the OPC (OLE (Object Linking and Embedding) for Process Control) toolbox and the configuration software of WinCC is carried out by MATLAB. The control algorithm is based on the data obtained by MATLAB. The results of the calculations are fed back to WinCC. A S7-300 PLC is adopted for the controller implementation. The interface is driven by an industry standard PLC and provides the programming flexibility to adapt to most measurement systems. The transmission of profiles to the operator station is done through industrial Ethernet communication. The PLC (Programmable Logic Controller) and supporting electronics form a compact unit for ease of installation and maintenance.

Modeling of CD Control System
The mathematical model of the CD control system is typically formed as follows: the CD process, i.e., the actuator input-output relationship, is described with a fixed, approximate non-square interaction matrix model, usually assuming identical and equally spaced actuators. It is further assumed that the response is separable, and therefore, the dynamic response of the CD actuator set can be described with a scalar first-order plus delay model. For the different actuator sets, the width of the actuators differ and their spatial responses are very different. The dynamic responses of the actuator sets vary significantly and are dependent on the machine speed and the distance between the actuating and measurement locations.
Like many other sheet and film processes, the simulation model assumes the first-order plus time delay dynamics between the actuator and the sensor. For the input array to the output array, a Toeplitz symmetric input-output interaction matrix is used.
The standard paper machine CD process model is given by where Y(s) ∈ C m and U(s) ∈ C n are the measurement profile and the actuator set-point profile, respectively, and G 0 ∈ R m×n is the spatial interaction coupling matrix that describes the response of the process. It also describes the mapping from the actuators to measurements; each component of the actuator vector U corresponds to a column of G 0 . m and n are the numbers of sheet properties and controlled actuator arrays, respectively. G 0 is a block diagonal matrix whose number of columns is equal to the number of actuators, and the number of rows is equal to the number of measurements, where g ij refers to the influence coefficient of the j th actuator on the i th measurement point. For a typical paper machine, if all the CD drive arrays and all the control characteristics are included in a control system, there may be up to 600 individual drives and 6000 profile measurements. In other words, the dimensions of most group process models G(s) in Equation (1) could be 6000 × 600. h(s) is the first-order-plus-time-delay transfer function of the process, and τ is a delay time of the process. The expressions of Y(s), U(s), G 0 ∈ R m×n , and h(s) are as follows.
The model is based on three assumptions related to in-compressible fluid, irrotational flow, and small amplitude waves. As shown in Figure 5, the static CD profile of the sheet is a two-dimensional wave. Amplitude is a function of position x and time t. It shows the pulp with surface waves on the wire. The x-axis shows the CD, and the y-axis shows the thickness of the pulp. In order to establish the spatial CD responses model, edge effects are studied based on the actuator response model of lip adjustment. The process of establishing the response model of the lip adjustment process is as follows.
For abscissa value x (measuring point position) and ordinate value (time) y, assume that ) , ( t x η is the function of the "S" type wave of the CD deviation profile. The origin of the coordinate system is placed at the interface of the fluid and the bed(wire). The wire represents the porous bed. The amplitude of the wave is shown as g, and the mean value is shown as h.
In order to establish the spatial CD responses model, edge effects are studied based on the actuator response model of lip adjustment. The process of establishing the response model of the lip adjustment process is as follows.
For abscissa value x (measuring point position) and ordinate value (time) y, assume that η(x, t) is the function of the "S" type wave of the CD deviation profile.
Regarding the assumptions outlined above, the slurry flow is analyzed according to the incompressible fluid. The conservation of mass on a surface "S" that surrounds the incompressible moving fluid implies the continuity equation.
It implies that the sum of the velocity gradients in any direction should be zero. Assuming an irrotational flow, then the velocity potential φ exits, the gradient of which is the velocity of the fluid, and the continuity Equation (3) reduces to a Laplace equation: Newton's second law of motion is then applied to fluid inside "S" The boundary condition at the surface is The dynamic boundary condition is After linearizing the dynamic boundary condition, Equation (7) can be simplified as Finally, the function of the "S" type wave of the CD deviation profile can be solved by a multi-objective model. The objective functions are The boundary conditions are The boundary condition is used to separate variables. Consider the following form of surface disturbance. η The real part of the Equation (11) is only considered When the CD profile does not vary with time, The response modeling method and identification method of the slice lip actuator (continuum device) and dilution actuator (discrete device) are similar.
So, the response model of the dilution water valve of g(x, t) has a similar form of expression with η(x, t). As the CD basis weight does not vary with time, g(x, t) can be simplified to g(x). For a continuous sheet-forming process, the relationship between the two-dimensional sheet change g(x) and CD control action to the j th CD actuator can be expressed as Equation (14).
It is shown that the spatial response of a single controller can be represented by the extension of a set of orthogonal functions, where the spatial coordinate x is a scalar real number, γ is a gain parameter, and ζ is a width parameter. The parameters γ and ζ define the linear transformation of the response by stretching it vertically and horizontally, respectively. The attenuation parameter α changes the size of the negative lobes of the response. The divergence parameter β defines the presence of two maxima in the response and the distance between these two maxima.
The shape of the spatial response represented by Equation (14) is illustrated in Figure 6. For the line board shape, . For the normal shape,

Identification of Coupling Matrix
In the CD control system, the measured dimension m is usually much higher than that of the actuator setpoints dimension n (generally, ). In the experimental simulation, taking a dilution water hydraulic headbox with 64 dilution water valves as an example, the measurements and the actuator form a large-scale interaction coupling matrix with dimensions of 64 320× .
Experiments have shown that when a single dilution water valve is operated, the affected area includes 26 measured data points on both sides of the valve [8].
The number of measurements is an integer multiple of the number of actuators in this example, such as n m × = 5 , which means that one actuator maps to five measured data points. Each actuator affects 26 measured data points on each side. The actuator has the strongest influence on the measured data point at the downstream position, and the coupling coefficient is the largest. As the distance increases, the impact will gradually decrease. If there are more than 26 lateral intervals, the effect is weak [9]. The mapping relationship between the actuators and the measured data is shown

Identification of Coupling Matrix
In the CD control system, the measured dimension m is usually much higher than that of the actuator setpoints dimension n (generally, 3n < m < 10n). In the experimental simulation, taking a dilution water hydraulic headbox with 64 dilution water valves as an example, the measurements and the actuator form a large-scale interaction coupling matrix with dimensions of 320 × 64. Experiments have shown that when a single dilution water valve is operated, the affected area includes 26 measured data points on both sides of the valve [8].
The number of measurements is an integer multiple of the number of actuators in this example, such as m = 5 × n, which means that one actuator maps to five measured data points. Each actuator affects 26 measured data points on each side. The actuator has the strongest influence on the measured data point at the downstream position, and the coupling coefficient is the largest. As the distance increases, the impact will gradually decrease. If there are more than 26 lateral intervals, the effect is weak [9]. The mapping relationship between the actuators and the measured data is shown in Figure 7, where the solid line indicates the greatest impact, and the dotted line indicates the gradual weakening.  where 26 is the half-width of the response of one actuator and m represents the scanning space between adjacent actuators.  Based on such an example, a large-scale interaction coupling matrix G 0 can be described as shown in Equation (15). The parameters g 1 , · · · , g 26 describe the spatial response of the process, where Symmetry 2020, 12, 149 9 of 22 26 is the half-width of the response of one actuator and m represents the scanning space between adjacent actuators.
As the number of n actuators in the system, the process model can be given by where one of the columns g j ∈ R m represents the sample space response of the j th actuator, which is determined by c j is the measured data coordinate of the spatial response centre for the j th actuator. Based on the Gaussian shape, the coupling matrix G 0 can be calculated. G 0 is a sparse strip matrix with a large number of zero elements (see Figure 8), where the non-zero data reflect the coupling coefficient of each CD position and the strip width reflects the strength of the coupling. Figures 9 and 10 show the non-zero data distribution. It can be seen that the elements on the main diagonal are larger, and the remaining elements are small. This indicates that G 0 is a non-square, high-dimensional, symmetric interaction matrix with small-range coupling characteristics.
As the symmetric inherent in the physical architecture of the web forming processes [10], the structure of a Toeplitz symmetric matrix is usually illustrated in the interaction matrix G 0 . The same constant element is repeated along each diagonal of the matrix in a Toeplitz symmetric matrix. The constant square band-diagonal interaction matrix G 0 with m independent spatial coupling partials in Equation (15) can be expressed as a Toeplitz symmetric matrix: G 0 = Toeplitz(g 1 , · · · , g m , 0, · · · , 0). Symmetry 2020, 12, x FOR PEER REVIEW 10 of 23 diagonal are larger, and the remaining elements are small. This indicates that 0 G is a non-square, high-dimensional, symmetric interaction matrix with small-range coupling characteristics.

Block and Decomposition Control of Coupling Matrix
For non-square high-dimensional coupling systems, it is necessary to adopt a suitable method to reduce the system dimension and convert it into a square system for decoupling. According to the strip coupling width of the above interaction matrix, G 0 is standardized as shown in Figure 11, where p is the coupling width, which indicates the number of affected actuators on the adjacent side when adjusting an actuator (p = 2 in this paper), and 2p + 1 is the strip width.

Block and Decomposition Control of Coupling Matrix
For non-square high-dimensional coupling systems, it is necessary to adopt a suitable method to reduce the system dimension and convert it into a square system for decoupling. According to the strip coupling width of the above interaction matrix, 0 G is standardized as shown in Figure 11, where p is the coupling width, which indicates the number of affected actuators on the adjacent side when adjusting an actuator ( 2 = p in this paper), and 1 2 + p is the strip width. Figure 11. The standardized incidence matrix.

Coupling Matrix Blocking
For simplifying the analysis, a blocking-based triangulation strategy is adopted for the coupling matrix 0 G . Along the diagonal, the coupling matrix can be divided into q q × block matrices. Thus, a new block gain matrix of 0 G can be partitioned as [H] represents the smallest integer, which is no less than H [11]. 1 G is an upper triangular matrix. The structures of 1 G and 2 G are given by Figure 11. The standardized incidence matrix.

Coupling Matrix Blocking
For simplifying the analysis, a blocking-based triangulation strategy is adopted for the coupling matrix G 0 . Along the diagonal, the coupling matrix can be divided into q × q block matrices. Thus, a new block gain matrix ofĜ 0 can be partitioned aŝ where q = m 2p+1 = 320 2×2+1 = 64 and [ ] is the rounding up function; for example; for example, [H] represents the smallest integer, which is no less than H [11]. G 1 is an upper triangular (2p + 1) × (2p + 1) matrix, and G 2 is a symmetric (2p + 1) × (2p + 1) matrix. The structures of G 1 and G 2 are given by For the new block gain matrix ofĜ 0 , the output variables are different from the original output of Y(s). The new CD process model is shown in Equation (20).
whereŶ(s) ∈ C n and U(s) ∈ C n ,Ŷ(s) indicates the low-dimensional measurement profile obtained by transformation matrix of F.Ŷ Based on such a division, the original CD control system with the dimensions of m × n is divided into q groups, and each group contains 2p + 1 subsystems, where the model of the subsystem can be expressed asŶ For the j th subsystem,Ŷ j (t) is the output of the measurement profile at time t. These measurements can only be made on paper that is dry enough at the end of the line. This causes a long time delay between the actuators on the hydraulic headbox and measurements from the scanner at the end of the paper machine.û j (t − τ) is the input of the j th CD actuator set-point profile at t − τ time. u j−1 (t − τ) and u j+1 (t − τ) represent the ( j − 1) th and ( j + 1) th CD actuator set-point profile at t − τ time. E(t) is the noise sequence. The symmetric parts of 2 and 3 in Equation (22) correspond to the paper CD basis weight boundaries. The block matrix G 1 , G 1 T is located at the boundary of the interaction matrix, and G 2 is located on the main diagonal of the interaction matrix. When the effect controlled by G 2 is non-zero and the other parts of the control effect are zero, G 0 can be converted to a main diagonal matrix. Then, the system can be simplified as a main diagonal system, and the coupling can be reduced. Based on the input variable adjustment, the decomposition algorithm is adopted to eliminate the control effect of sub-blocks 2 and 3 in this paper.

Decomposition Algorithm Control
By a blocking-based triangulation strategy, q multivariable subsystems that have dimensions of 2p + 1 are obtained. To eliminate the edges of each subsystem, a decomposition algorithm is adopted in this part. The decomposition algorithm steps are as follows.
Divide the control variables into groups according to the horizontal position, each group containing five control points.
Construct a new constraint applied to the control input U.
The construction principle of the new constraint is that the control action cannot be applied to all the actuators in the same group at the same time [12], such as that the u j actuator cannot be set for simultaneous actions in subsystem j, although the actuators in different groups can be set simultaneously.
Cyclically input the control variables in chronological order. Actuators belonging to different groups can act simultaneously [13]. For example, from the time of t to t + 2p + 1, the following constraints are imposed on each subsystem control input, as shown in Table 1. Table 1. The form of the subsystem control input.
Based on the real mapping matrix method introduced by Duncan [14], the transformation matrix of F can be simplified as described in Equation (24). The Equation (21) can be described in Equation (25).

of 22
For the new low-dimensional square system, the interaction matrix with a number of dimensions equal to the number of actuators and a one-to-one correspondence between the CD measuring point and the actuator is solved. The new interaction matrix G 0 can be given as follows Moreover, the new interaction matrix G 0 is a Toeplitz symmetric matrix. G 0 = Toeplitz( f ) [15], where f = [1, 0.7, 0.3, 0 · · · 0].

Interpolation Decoupling Strategy for CD Control System
For the high-dimensional Toeplitz symmetric interation matrix with small-scale coupling characteristics, a decoupling network is added between the controller and the controlled object. By the decoupling network, the controlled object is transformed into a principal diagonally dominant system. The control structure with a decoupler is shown in Figure 12. the decoupling network, the controlled object is transformed into a principal diagonally dominant system. The control structure with a decoupler is shown in Figure 12. The coupling matrix is a static matrix, which does not change with time. Then, the decoupler is also a static matrix [16]. For the new associated coupling matrix obtained above, the generalized inverse of a 64-dimensional Toeplitz matrix can be calculated. The product of a matrix and its adjoint matrix is a diagonal matrix in mathematics; that is, From the point of view of simplified design, the decoupler can be designed as (29) The detailed structure is shown in Equation (30).
It can be seen that the decoupling matrix of K is a 64-dimensional symmetric matrix, and each element is not zero. Then, the number of decoupling network branches that need to be calculated is The coupling matrix is a static matrix, which does not change with time. Then, the decoupler is also a static matrix [16]. For the new associated coupling matrix obtained above, the generalized inverse of a 64-dimensional Toeplitz matrix can be calculated. The product of a matrix and its adjoint matrix is a diagonal matrix in mathematics; that is, From the point of view of simplified design, the decoupler can be designed as Symmetry 2020, 12, 149

of 22
The detailed structure is shown in Equation (30).
It can be seen that the decoupling matrix of K is a 64-dimensional symmetric matrix, and each element is not zero. Then, the number of decoupling network branches that need to be calculated is 64 × 64/2 = 2048, and the decoupling speed is slower with the increase of the system dimensions. In order to achieve fast decoupling of the controlled system, it is necessary to further analyze the system and simplify the decoupling network.
The number of basis weight measurements is typically about 360, and the number of actuators is about 64 (i.e., five data per actuator). In order to reduce the number of calculated measurements, the corresponding deviation of each actuator is obtained by averaging the corresponding five data.
However, the average for each actuator may not catch a peak if there is one between the actuators (see Figure 13, the dilution valve i exists between i 0 and i 1 ). In this case, the dilution valve is interpolated, and the more dilution valves that are interpolated, the smaller the control deviation after the decoupling matrix. Based on the coupling matrix mentioned above, the interpolation matrix method is proposed to decouple the coupling system quickly. It is assumed that a point of the dilution water valve is inserted in each equidistant area to obtain twice the number of dilution water valves. Strictly, (twice the number of actuators-1). Compared with the control method of the decoupling matrix in Figure 14, the average section width of the control deviation is halved, as shown in Figure 15. Then, the peak value of the basis weight profile between the actual actuators has a great impact on the average control deviation [17]. In addition, the size of the interference factor matrix is also doubled. This control deviation technique by the interpolated decoupling matrix is named the 'Interpolated Decoupling Matrix Method' [18].  It is assumed that a point of the dilution water valve is inserted in each equidistant area to obtain twice the number of dilution water valves. Strictly, (twice the number of actuators-1). Compared with the control method of the decoupling matrix in Figure 14, the average section width of the control deviation is halved, as shown in Figure 15. Then, the peak value of the basis weight profile between the actual actuators has a great impact on the average control deviation [17]. In addition, the size of the interference factor matrix is also doubled. This control deviation technique by the interpolated decoupling matrix is named the 'Interpolated Decoupling Matrix Method' [18]. It can be seen that the decoupling matrix of K is a 64-dimensional symmetric matrix, and each 348 element is not zero. Then the number of decoupling network branches need to be calculated is 64 × 349 64/2=2048, and the decoupling speed is slower with the increase of system dimension. In order to 350 achieve fast decoupling of the controlled system, it is necessary to further analyze the system and 351 simplify the decoupling network.

352
The number of basis weight measurements is typically about 360, and the number of actuators is 353 about 64 (i.e., 5 data per actuator). In order to reduce the number of calculated measurements, the 354 corresponding deviation of each actuator is obtained by averaging the corresponding 5 data.   For the output of the dilution valve actuator U(s) = [u 1 (s), u 2 (s), · · · , u n (s)] T , assuming the output of the i 0 th dilution valve is u i0 , when the i 0 th dilution valve operates, it is necessary to adjust the opening of the left and right dilution valves in a certain proportion, while for the i 0 th region, there is a combined control regulation. The output u i0 of the i 0 th dilution valve can be expressed as where a j is the interference factor obtained for each existing actuator. P j is the control variable obtained from control deviation by the PI (Proportion-Integral) controller for each existing actuator.
In the interpolated decoupling matrix method, the coupling factor matrix is twice the original. The output u i0 of the i 0 th dilution valve is expressed as follows.
where b j is the interference factor obtained for the actuator, assuming that the number of actuators is equal to the number of actual actuators ×2. P j is the manipulated variable obtained from control deviation, assuming that the number of actuators is equal to twice the number of actual actuators. After the interpolation input of the controller, when the opening of the dilution valve changes, the influence on the left and right sides of the slurry flow will be reduced [19]. Therefore, the correlation matrix can be further transformed into a new weak coupling correlation matrix G 0 , which is expressed as For the new coupling matrix, the new 64-dimensional Toeplitz matrix is singularity judged by det( G 0 ) 0, and the diagonal decoupler can be designed as A new decoupling matrix is obtained by truncating the new decoupling matrix K , the value less than 10 −3 is approximated to 0, and the banded symmetric matrix K as shown in Equation (35) is obtained.
The non-zero data distribution of K is shown in Figures 16 and 17. It can be shown that the main diagonal elements are larger and the other elements are smaller.
The non-zero data distribution of ' K is shown in Figures 16 and 17. It can be shown that the main diagonal elements are larger and the other elements are smaller.  It can be seen that the decoupling network is also a Toeplitz matrix. There are only 436 non-zero branches in the strip of ' K , and the decoupling matrix ' K is a symmetric matrix. It means that only The non-zero data distribution of ' K is shown in Figures 16 and 17. It can be shown that the main diagonal elements are larger and the other elements are smaller.  It can be seen that the decoupling network is also a Toeplitz matrix. There are only 436 non-zero branches in the strip of ' K , and the decoupling matrix ' K is a symmetric matrix. It means that only It can be seen that the decoupling network is also a Toeplitz matrix. There are only 436 non-zero branches in the strip of K , and the decoupling matrix K is a symmetric matrix. It means that only 186 branches need to be calculated in the process of decoupling calculation. The decoupling branch is greatly reduced, which improves the decoupling speed of the system. This result indicates that the original complex decoupling network is transformed into a simple and easy-to-implement decoupling network based on the interpolation method.

Results of Simulation and Application
For the new system with 64 dimensions, there are 64 actuators mapping 64 regions, which means that the coupling between each region is eliminated. Thus, the control of a high-dimensional large-scale system can be transformed into multiple single-loop control. The block diagram of subsystem parallel control is shown in Figure 18. 186 branches need to be calculated in the process of decoupling calculation. The decoupling branch is greatly reduced, which improves the decoupling speed of the system. This result indicates that the original complex decoupling network is transformed into a simple and easy-to-implement decoupling network based on the interpolation method.

Results of Simulation and Application
For the new system with 64 dimensions, there are 64 actuators mapping 64 regions, which means that the coupling between each region is eliminated. Thus, the control of a high-dimensional large-scale system can be transformed into multiple single-loop control. The block diagram of subsystem parallel control is shown in Figure 18. A decoupled CD control system for corrugated paper with a basis weight of 133 g/m 2 was implemented in an actual project of a paper mill.
The basis weight of corrugated paper is 140 g/m 2 , and the allowable basis weight fluctuation range is 7 ± g/m 2 . In order to reduce the pulp consumption in the paper production process, the actual project stabilizes the paper basis weight at 133 g/m 2 for automatic control. Based on the criterion of equal basis weight in each region, the following model is established.  A decoupled CD control system for corrugated paper with a basis weight of 133 g/m 2 was implemented in an actual project of a paper mill.
The basis weight of corrugated paper is 140 g/m 2 , and the allowable basis weight fluctuation range is ±7 g/m 2 . In order to reduce the pulp consumption in the paper production process, the actual project stabilizes the paper basis weight at 133 g/m 2 for automatic control. Based on the criterion of equal basis weight in each region, the following model is established.
g 11 u 1 + g 12 u 2 + g 13 u 3 = b g 21 u 1 + g 22 u 2 + g 23 u 3 + g 24 u 4 = b g 31 u 1 + g 32 u 2 + g 33 u 3 + g 34 u 4 + g 35 u 5 = b . . . where U = u 1 u 2 · · · u 63 u 64 is the output of the dilution valve actuator and b is the normalized setting value of the CD basis weight.
According to the calculated regulation of each dilution valve, the opening of the dilution valve is changed as the initial value of the input variable of the control system changes.
In the decoupled CD system, each single loop controlled by 64 actuators has similar physical characteristics. It is the first-order plus time delay (FOPTD) dynamic response of the process. Each loop controller design is relatively simple. Take a two-loop adjustment for example. For each single-loop design controller, the response curve under the action of the controller is shown in Figures 19 and 20. According to the figures, the first-step input produces a good output y 1 (t) when applied to the system, and y 2 (t) is almost 0. The effect is found to be similar when the second input is used alone. The loops do not affect each other. It is proven that the coupling relationship has been eliminated very well. Moreover, the control quality is improved.
According to the calculated regulation of each dilution valve, the opening of the dilution valve is changed as the initial value of the input variable of the control system changes.
In the decoupled CD system, each single loop controlled by 64 actuators has similar physical characteristics. It is the first-order plus time delay (FOPTD) dynamic response of the process. Each loop controller design is relatively simple. Take a two-loop adjustment for example. For each single-loop design controller, the response curve under the action of the controller is shown in Figures 19 and 20. According to the figures, the first-step input produces a good output ) ( 1 t y when applied to the system, and ) ( 2 t y is almost 0. The effect is found to be similar when the second input is used alone. The loops do not affect each other. It is proven that the coupling relationship has been eliminated very well. Moreover, the control quality is improved.
According to the calculated regulation of each dilution valve, the opening of the dilution valve is changed as the initial value of the input variable of the control system changes.
In the decoupled CD system, each single loop controlled by 64 actuators has similar physical characteristics. It is the first-order plus time delay (FOPTD) dynamic response of the process. Each loop controller design is relatively simple. Take a two-loop adjustment for example. For each single-loop design controller, the response curve under the action of the controller is shown in Figures 19 and 20. According to the figures, the first-step input produces a good output ) ( 1 t y when applied to the system, and ) ( 2 t y is almost 0. The effect is found to be similar when the second input is used alone. The loops do not affect each other. It is proven that the coupling relationship has been eliminated very well. Moreover, the control quality is improved.  Setting the paper basis weight at 133 g/m 2 , the dilution valves are controlled by the S7-300 PLC. When the working condition is stable, a monitoring picture of a certain period is shown in Figure 21, and 50 scanning points are measured randomly. The basis weight deviation curve is shown in Figure 22. It can be seen that the fluctuation deviation is small, so that the algorithm is proven to be very practical. Setting the paper basis weight at 133 g/m 2 , the dilution valves are controlled by the S7-300 PLC. When the working condition is stable, a monitoring picture of a certain period is shown in Figure 21, and 50 scanning points are measured randomly. The basis weight deviation curve is shown in Figure  22. It can be seen that the fluctuation deviation is small, so that the algorithm is proven to be very practical.  To further verify the effectiveness of the algorithm, a low-dimensional CD system was set up in the experiment, with the CD measured data points 24  m and the number of dilution valves being the same as m ( m n  ). The coupling width of the interaction matrix was 2, and the basis weight fluctuations after the decoupling control were compared. As shown in Figures 23 and 24, it can be seen that the basis weight fluctuation range is significantly reduced, and the product quality is improved based on the decoupling control strategy.  Setting the paper basis weight at 133 g/m 2 , the dilution valves are controlled by the S7-300 PLC. When the working condition is stable, a monitoring picture of a certain period is shown in Figure 21, and 50 scanning points are measured randomly. The basis weight deviation curve is shown in Figure  22. It can be seen that the fluctuation deviation is small, so that the algorithm is proven to be very practical.  To further verify the effectiveness of the algorithm, a low-dimensional CD system was set up in the experiment, with the CD measured data points 24 = m and the number of dilution valves being the same as m ( m n = ). The coupling width of the interaction matrix was 2, and the basis weight fluctuations after the decoupling control were compared. As shown in Figures 23 and 24, it can be seen that the basis weight fluctuation range is significantly reduced, and the product quality is improved based on the decoupling control strategy.  To further verify the effectiveness of the algorithm, a low-dimensional CD system was set up in the experiment, with the CD measured data points m = 24 and the number of dilution valves being the same as m (n = m). The coupling width of the interaction matrix was 2, and the basis weight fluctuations after the decoupling control were compared. As shown in Figures 23 and 24, it can be seen that the basis weight fluctuation range is significantly reduced, and the product quality is improved based on the decoupling control strategy. Setting the paper basis weight at 133 g/m 2 , the dilution valves are controlled by the S7-300 PLC. When the working condition is stable, a monitoring picture of a certain period is shown in Figure 21, and 50 scanning points are measured randomly. The basis weight deviation curve is shown in Figure  22. It can be seen that the fluctuation deviation is small, so that the algorithm is proven to be very practical.  To further verify the effectiveness of the algorithm, a low-dimensional CD system was set up in the experiment, with the CD measured data points 24 = m and the number of dilution valves being the same as m ( m n = ). The coupling width of the interaction matrix was 2, and the basis weight fluctuations after the decoupling control were compared. As shown in Figures 23 and 24, it can be seen that the basis weight fluctuation range is significantly reduced, and the product quality is improved based on the decoupling control strategy.

Conclusions
In view of the characteristics of strong multivariable coupling in a cross-directional (CD) basis weight control system, a conventional decoupling control strategy may be difficult to meet the actual needs of a CD basis weight control system. Thus, it is necessary to study the effective decoupling or compensation control strategy. In this paper, a general method of transforming the control problem of a high-dimensional large-scale system into a single-loop group control problem is proposed based on the matrix block method and the interpolation decoupling strategy. Based on a model of the CD basis weight profile, the coupling properties of a CD control system are analyzed by numerical simulation. The matrix block method and model reduction were applied in the interaction coupling matrix of dimensions of 64 320× , and the original 64 320× non-square CD control system can be converted to a new low-dimensional CD system with dimensions of only 64 64× . For the 64-dimensional Toeplitz symmetric subsystem with small-scale coupling characteristics, an interpolation decoupling algorithm is adopted, and then a decoupling compensator with the structure of a symmetric Toeplitz matrix was obtained. Based on the multi-actuator parallel control strategy, a decoupled multivariable system is achieved. The coupling calculation in this paper is based on the data collected from the actual paper mill. An interpolation decoupling controller is designed to transform an MIMO system into multiple SISO systems. Experimental results show that the method has an obvious decoupling control effect, and this control strategy can provide help for the actual project. In addition, the multidimensional interpolation decoupling control strategy and a real-time decomposition algorithm can be applied to various processes whose input and output numbers are different. We sincerely thank them for the funding of the project.