**1. Introduction**

In order to improve labor productivity and reduce investment costs, large-capacity prebaked reduction cells are currently used in various enterprises. Due to their high efficiency and low energy consumption, pre-baked cells of 400–600 kA have gradually become the mainstream cell type in China's aluminum electrolytic industry [1]. The capacity of the reduction cell is continuously increasing, while the auxiliary facilities and the intelligent control technology of aluminum electrolytic are relatively backward. Thus, the problems of the local anode effect and local precipitation in the large reduction cell have increasingly become the main factors of instability in the production process [2], which has caused serious economic losses and some casualties. The main cause of these problems is the uneven distribution of alumina concentration in the anode bottom surface of the large aluminum reduction cell [3]. Through some studies and experiments, it is known that the concentration of alumina is generally controlled between 1.5% and 3.5%. Currently, the change of groove resistance and the concentration of alumina is basically linear and easy to identify, and the current efficiency is also the highest [4]. If the concentration of alumina is excessively high, there will be problems such as increased energy consumption, fluctuation of aluminum liquid layer, etc., and even induced precipitation at the bottom of

**Citation:** Cui, J.; Wang, P.; Li, X.; Huang, R.; Li, Q.; Cao, B.; Lu, H. Multipoint Feeding Strategy of Aluminum Reduction Cell Based on Distributed Subspace Predictive Control. *Machines* **2022**, *10*, 220. https://doi.org/10.3390/ machines10030220

Academic Editor: Xiang Li

Received: 23 February 2022 Accepted: 17 March 2022 Published: 21 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the cell and crusting on the cell side which reduce the service life of the reduction cell [5]. An excessively low concentration of alumina will lead to frequent anode effects in the reduction cell. Once the anode effect occurs, the cell voltage rises sharply. Meanwhile, Moxnes et al. [6] found that when the alumina concentration distribution is more uniform, the current efficiency of the reduction cell is higher, and the probability of abnormal cell conditions such as the anode effect is lower. Therefore, adjusting the feeding interval of each feeder of the large-scale aluminum reduction cell and accurately controlling the amount of alumina feeding to make the alumina concentration evenly distributed has become a general concern and urgen<sup>t</sup> issue that has extremely important practical significance for the further development of large and super large aluminum electrolytic technology [5].

The uniform distribution of alumina concentration plays a vital role in the stable operation and efficient production of large aluminum electrolysis. In order to achieve this goal, many scholars have carried out in-depth research on alumina. In [7], a method combining fuzzy control theory and expert experience was proposed to control the alumina concentration by changing the feeding interval. In [4], by improving the crust breaking and feeding device and control system of the reduction cell, the feeding interval was set separately for each feeding point, and the single-point precise feeding control is achieved. However, for large aluminum reduction cells with multiple feeders, single-point feeding control method cannot effectively control the uniform distribution of the alumina concentration. With the development of soft measurement technology [8,9] and distributed data measurement technology, more and more scholars have begun to try to apply soft-sensing the aluminum electrolysis industry. The least squares support vector machine method for the alumina concentration soft measurement model is established in [10]. In [11], in order to obtain more accurate results, a soft-sensor model of alumina concentration was proposed that introduces time series to optimize the input parameters of a deep belief network (DBN). In [12], a KPI was developed through a probabilistic soft sensing based on maximizing the coefficient of determination to estimate the alumina concentration. An improved Kalman filter for the soft sensing of alumina concentration is presented in [13]. Therefore, intelligent control methods based on the soft measurement model of the alumina concentration emerge endlessly. The generalized regression neural network (GRNN) was adopted to identify the alumina concentration model and a fuzzy cerebellar model neural network (FCMAC) controller was proposed for the feeding equipment to control the alumina concentration in [14]. A data-driven intelligent control system based on the least squares support vector machine alumina concentration soft sensing model was proposed in [15]. An extended Kalman filter (EKF) was used to estimate the local alumina concentration to design a multivariable blanking control strategy in [16]. However, these control methods ignore the influence of each feeding port on the alumina concentration near other feeding ports due to the flow of electrolytes during the aluminum reduction process, and only consider the overall parameters such as cell voltage and cell resistance, and do not fully utilize important distribution parameters. In order to fully understand the distribution of an alumina concentration in aluminum reduction cells, the dissolution and diffusion of alumina were studied in [17–20]. In [21], a multi-coupling distributed alumina simulation model was constructed by ANSYS, and the uniform distribution of alumina concentration was achieved through the simulation model. However, it is still very difficult to establish an accurate model through mechanism analysis. Some scholars have studied the production process using a large amount of data. At present, the most advanced technology was that the relationship between the distributed current and the feeding rate was established using random forest, and the stability of the aluminum reduction cell was maintained by controlling the distributed current to be consistent in [22]. However, the final simulation results did not verify the uniform distribution of alumina concentration.

In order to solve the problem of the uneven distribution of alumina concentration and adapt to the aluminum electrolysis process with difficulty in mechanism modeling, a distributed subspace predictive control data-driven method is applied to large aluminum reduction cells, and the subsystem model is established by fully using distributed data. The distributed control algorithm is used to realize the distributed feeding of multiple feeders in large electrolytic cell considering the influence of each feeder, and the simulation is carried out in MATLAB. The simulation results show that the application of this method in a large aluminum reduction cell is feasible, and it has important guiding significance for realizing the uniform distribution of alumina concentration in a large aluminum reduction cell in industrial processes.

Compared with the existing research on alumina concentration control strategy, the contributions of this paper are as follows:


The rest of this paper is organized as follows. In Section 2, a distributed feeding control strategy for aluminum electrolysis is proposed. Section 3 introduces the distributed subspace predictive control algorithm and discusses the implementation of the proposed algorithm in aluminum electrolysis in detail. In Section 4, the feasibility of the algorithm is verified by MATLAB simulation. In Section 5, relevant conclusions are given.

### **2. Design of Distributed Feeding Control Scheme for Aluminum Electrolysis**

The top view of a 400 kA large aluminum reduction cell in an aluminum plant is shown in Figure 1. FD1, FD2, FD3, FD4, FD5, and FD6 are the six feed ports of the cell. Twenty-four anode guides are on the B side. Due to the flow of electrolytes caused by carbon dioxide and carbon monoxide gas produced by anode, electromagnetic field, temperature and concentration difference, when feeding at any feeding port, the alumina concentration in other regions will be affected to some extent. Taking the six feeding ports as centers, the aluminum reduction cell is divided into six subsystems. For aluminum electrolysis, which is a complex large system composed of multiple mutually influencing subsystems, distributed control cannot only better consider the impact of feeding between subsystems but the computation complexity is greatly reduced compared to centralized control. As predictive control has been widely used in engineering applications and has high control accuracy, distributed model predictive control has received more attention from scholars.

**Figure 1.** Top view of 400 kA large aluminum reduction cell.

The main idea of the distributed predictive control algorithm is to transform a largescale online optimization problem into a small-scale distributed optimization of each subsystem, and at the same time, each subsystem communicates and shares information

through the network, thereby improving the control performance of the system. For the general distributed model predictive control method, the design of the controller must be based on accurate modeling. Due to the complexity of the aluminum reduction process, there are many difficulties in the research on the multi-point distributed feeding control of the aluminum electrolysis process. Firstly, the aluminum reduction process is a dynamic system with complex physical and chemical reactions, multi-coupling and large delay [24]. Secondly, it is difficult to establish a precise distributed multi-point feeding mechanism model due to the extremely complex environment such as high temperature and strong corrosion in industrial aluminum reduction cells. Finally, it is difficult to obtain the coupling relationship between each subsystem. The subspace identification method is not limited to the prior structure information and mechanism model of the system but directly uses the historical input and output data to solve the prediction model [25,26], and can obtain the prediction model of each subsystem through parameter decomposition. This method is more suitable for complex large systems composed of multiple subsystems and it is difficult to establish an accurate mechanism model [27,28]. Therefore, the data-driven subspace identification method and predictive control are designed in a control system design framework, which is applied to the aluminum electrolysis system with model uncertainty and has better control performance [29–31].

For large aluminum reduction cells, since the distributed alumina concentration data cannot be obtained in real-time, the development of distributed current measurement technology provides the basis for its soft measurement. According to the relevant mechanism of an aluminum reduction cell, there is a close relationship between the distributed alumina concentration and distributed current. This research group has also conducted in-depth research on the soft measurement of distributed alumina concentration [32,33]. Therefore, the distributed alumina concentration data required in this paper can be obtained by using the soft sensing model, to carry out the follow-up work.

For the 400 kA aluminum reduction cell, the structure of the predictive control principle of one of the subsystems is shown in Figure 2. It is mainly composed of a subspace prediction model and distributed controller.

**Figure 2.** Predictive control diagram of aluminum electrolysis subsystem.

The subspace prediction model is a model of the entire system identified by the input and output data. After parameter decomposition, the prediction model of the subsystem can be obtained. The prediction model of each subsystem includes the influence of other subsystems on itself. The aluminum electrolysis system shown in Figure 1 can be described as due to the flow of electrolyte, whilst other subsystems cause changes in alumina concentration to a subsystem. Each subsystem has a separate controller to control the feeder responsible for supplying the alumina powder, and a distributed control algorithm is designed under the condition that the subsystems can communicate with each

other. At each moment, each subsystem solves the optimal control signal of its system when the optimal control signals of other subsystems are known, and the global optimality is also guaranteed in the case of achieving local optimality. For the aluminum electrolysis system, the advantage of this method is that the six feeding ports can change the original group feeding or timing feeding strategy so that the six feeding ports can consider the influence of other feeding ports. The purpose of distributed feeding is to make the alumina concentration distribution more uniform, reduce the occurrence of local precipitation and local anode effect, ensure the stable and efficient operation of the entire cell, and improve the production efficiency of the aluminum plant.

### **3. Distributed Subspace Predictive Control**

The basic idea of the data-driven distributed subspace predictive controller design is to first obtain the input and output data of length n, then use these data to solve the distributed prediction model, and finally use the obtained distributed prediction model to design the controller [24]. According to the actual data collection situation on-site, the input variable *u*1, *u*2, ··· , *um*(*m* = 6) is determined as the alumina feeding amount of the six subsystems, and the output variable *y*1, *y*2, ··· , *ym*(*m* = 6) is determined as the alumina concentration for the six subsystems. Among them, the alumina feeding amount data are obtained by combining the dissolution and consumption mechanism of alumina and the feeding interval. According to the relevant mechanism of an aluminum reduction cell, there is a close relationship between the distributed alumina concentration and distributed current [22]. The soft sensor model of the distributed current and distributed alumina concentration is established by using the current data of a single anode guide rod and the alumina concentration data of different areas collected by field test. The alumina concentration data and alumina discharge data of six subsystems with n = 1000 can be obtained.

### *3.1. Data-Driven Distributed Prediction Model*

For an aluminum reduction cell with six subsystems, the output prediction model can be described as: at time *k*, the relationship between the output prediction vector *y*<sup>ˆ</sup>*t*(*k*)=[*y*ˆf−<sup>1</sup>(*k*) T, ··· , *y*ˆf−*<sup>m</sup>*(*k*) T] T composed of the alumina concentration prediction values of each subsystem and the input and output data is [30]:

$$
\hat{y}\_{\mathbf{f}}(k) = L\_w \cdot w\_{\mathbf{P}}(k) + L\_u \cdot u\_{\mathbf{f}}(k) \tag{1}
$$

where *Lw* ∈ *RmN*×2*mN* and *Lu* ∈ *RmN*×*mN* are the unknown parameter matrix, which is obtained by subspace identification, *N* is the length of the prediction window, *<sup>u</sup>*f(*k*) is the input vector composed of the future alumina feeding amount of each subsystem, and *<sup>w</sup>*p(*k*) is the vector composed of the past input and output data, respectively, defined as follows:

$$w\_{\mathbb{P}}(k) \triangleq \begin{bmatrix} w\_{\mathbb{P}^{-1}}(k) \\ \vdots \\ w\_{\mathbb{P}^{-m}}(k) \end{bmatrix}; \boldsymbol{u}\_{\mathbb{P}}(k) \triangleq \begin{bmatrix} u\_{\mathbb{f}\_{-1}}(k) \\ \vdots \\ u\_{\mathbb{f}\_{-m}}(k) \end{bmatrix}; \boldsymbol{w}\_{\mathbb{p}^{j}}(k) \triangleq \begin{bmatrix} u\_{\mathbb{f}\_{-1}}(k) \\ \vdots \\ u\_{\mathbb{f}\_{-m}}(k) \end{bmatrix}; \boldsymbol{z}\_{\mathbb{p}^{j}}(k) \triangleq \begin{bmatrix} u\_{\mathbb{p}^{j}}(k)^{\mathsf{T}}y\_{\mathbb{p}^{j}}(k)^{\mathsf{T}} \end{bmatrix}^{\mathsf{T}}; \boldsymbol{z}\_{\mathbb{p}^{j}}(k) \triangleq \begin{bmatrix} \boldsymbol{w}\_{\mathbb{p}^{j}}(k) \\ \boldsymbol{w}\_{\mathbb{p}^{j}}(k) \end{bmatrix}; \boldsymbol{z}\_{\mathbb{p}^{j}}(k) \triangleq \begin{bmatrix} \boldsymbol{g}\_{i}(k) \\ \vdots \\ \boldsymbol{u}\_{i}(k+N-1) \end{bmatrix}; \boldsymbol{z}\_{\mathbb{p}^{j}}(k) \triangleq \begin{bmatrix} \boldsymbol{g}\_{i}(k) \\ \vdots \\ \boldsymbol{g}\_{i}(k+N-1) \end{bmatrix}.$$

where *i* = 1, 2, ... , 6, the subscripts "*p*" and "*f*" represent the past and the future, respectively. In order to realize distributed feeding control, a distributed prediction model needs to be established. Equation (1) is decomposed into the following form:

$$
\begin{pmatrix} \mathfrak{g}\_1(k) \\ \mathfrak{g}\_2(k) \\ \vdots \\ \mathfrak{g}\_m(k) \end{pmatrix} = \begin{bmatrix} L\_{\mathfrak{u}}(1) \\ L\_{\mathfrak{u}}(2) \\ \vdots \\ L\_{\mathfrak{u}}(m) \end{bmatrix} \cdot w\_{\mathfrak{p}}(k) + \begin{bmatrix} L\_{\mathfrak{u}(1,1)} & L\_{\mathfrak{u}(1,2)} & \cdots & L\_{\mathfrak{u}(1,m)} \\ L\_{\mathfrak{u}(2,1)} & L\_{\mathfrak{u}(2,2)} & \cdots & L\_{\mathfrak{u}(2,m)} \\ \vdots & \vdots & \vdots & \vdots \\ L\_{\mathfrak{u}(m,1)} & L\_{\mathfrak{u}(m,2)} & \cdots & L\_{\mathfrak{u}(m,m)} \end{bmatrix} \cdot \begin{pmatrix} u\_1(k) \\ u\_2(k) \\ \vdots \\ u\_m(k) \end{pmatrix} \tag{2}
$$

Based on this decomposition, the alumina concentration prediction model of each subsystem can be obtained:

$$\mathcal{Y}\_{i}(k) = L\_{\mathfrak{u}(i,i)}\mathfrak{u}\_{i}(k) + \sum\_{\substack{j=1, j\neq i \atop \text{Educts of } j \text{ other subsystems,} \, i=1,2,...,6}}^{\text{m}} \frac{L\_{\mathfrak{u}(i,j)}\mathfrak{u}\_{j}(k)}{\text{m}^{\mathfrak{u}}(k)} + L\_{\mathfrak{w}(i)} \cdot \mathfrak{w}\_{\mathbb{P}}(k) \tag{3}$$

Each prediction model includes the influence of the feeding number of other subsystems on itself.

### *3.2. Design of Distributed Predictive Controller for Aluminum Reduction Cell System*

The control target of the distributed aluminum reduction cell system is to control the alumina concentration to track the reference value. After certain research and experiments, the reference alumina concentration value of each subsystem can be obtained:

$$R\_{ref} = \begin{bmatrix} \ r\_1 & r\_2 & \cdots & r\_m \end{bmatrix}^T \tag{4}$$

The control strategy of the distributed aluminum reduction cell system is shown in Figure 2. According to Equation (3), predicting the output of any subsystem requires knowing the input and output data of all subsystems at the previous time [28]. Therefore, the controller of each subsystem must have a communication function to transmit its information to other subsystems. The difference between distributed control and centralized control is that the global performance index can be expressed as the sum of the performance indexes of all subsystems [34]:

> *J*

$$l = \sum\_{i=1}^{m} J\_i \tag{5}$$

then the performance index of the *i*th subsystem can be expressed as

$$\min\_{\mu\_{\underline{t}\_{i\cdot i}}(k)} J\_{\bar{t}} = \left[ y\_{\underline{t}\dots i}^{\gamma}(k) - r\_{\bar{i}}(k) \right]^{\mathrm{T}} Q\_{\bar{t}} \left[ y\_{\underline{t}\dots i}^{\gamma}(k) - r\_{\bar{i}}(k) \right] + u\_{\underline{t}\dots i}(k)^{\mathrm{T}} R\_{\bar{i}} u\_{\underline{t}\dots i}(k) \tag{6}$$

where *ri*(*k*) is the reference input signal vector of the subsystem; and *Qi* and *Ri* are the positive definite weighting matrixes. Substituting Equation (3) into Equation (6), we can obtain: *m*

$$\begin{aligned} J\_i &= \left[ L\_{w(i)} w\_\mathcal{P} - r\_i(k) + \sum\_{j=1}^n L\_{u(i,j)} \cdot u\_j(k) \right]^\mathcal{I} Q\_i \\ \left[ L\_{w(i)} w\_\mathcal{P} - r\_i(k) + \sum\_{j=1}^m L\_{u(i,j)} \cdot u\_j(k) \right] &+ u\_{t-i}(k)^\mathcal{T} \mathcal{R}\_i u\_{t-i}(k) \end{aligned} \tag{7}$$

Differentiate the objective function to find the extremum:

$$\frac{\partial f\_i}{\partial u\_{\ell\_-i}(k)} = 0 \tag{8}$$

The controller for each subsystem can be obtained

$$\boldsymbol{\mu}\_{\boldsymbol{t}\_{-}(i)}(\boldsymbol{k}) = -\left[\boldsymbol{R}\_{i} + \boldsymbol{L}\_{\boldsymbol{u}(i,i)}^{\mathrm{T}} \boldsymbol{Q}\_{i} \boldsymbol{L}\_{\boldsymbol{u}(i,i)}\right]^{-1} \boldsymbol{L}\_{\boldsymbol{u}(i,i)}^{\mathrm{T}} \boldsymbol{Q}\_{i} \left[\boldsymbol{L}\_{\boldsymbol{w}(i)} \boldsymbol{w}\_{\mathbb{P}}(\boldsymbol{k}) - \boldsymbol{r}\_{i}(\boldsymbol{k}) + \sum\_{j=1, j\neq i}^{m} \boldsymbol{L}\_{\boldsymbol{u}(i,j)} \cdot \boldsymbol{u}\_{\boldsymbol{t}\_{-}(i)}(\boldsymbol{k})\right] \tag{9}$$

:

In actual control, we only put *<sup>u</sup>*f−*<sup>i</sup>* the first component of the controlled input into the future input data matrix and pass this matrix to other subsystems. When performing predictive control iteration and rolling optimization in the above steps, the optimal control quantity is always calculated. Therefore, the actual output alumina concentration of each subsystem of aluminum electrolysis can be synchronized with the reference alumina concentration to achieve the control target.

### *3.3. Determination of Parameters of Aluminum Reduction Cell Prediction Model and Design of Data-Driven Distributed Predictive Control Algorithm*

In order to obtain the prediction model of Equation (3), the above input and output data of n = 1000 are used to solve the parameter matrix *Lw* and *Lu* in Equation (1). The prediction step *N* is set to 5. Input data at time (0, 1, ··· , *N* − 1) and output data at time (0, 1, ··· , 2 *N* − 1) are used to predict the output at time (*<sup>N</sup>*, *N* + 1, ··· , 2 *N* − 1) according to Equation (1), as shown in Figure 3. After that, we move the time window and the input data at time (1, 2, ··· , 2 *N*) and output data at time (1, 2, ··· , *N*) are used to predict the output at time (*N* + 1, *N* + 2, ··· , 2 *N*) according to Equation (1), as shown in Figure 4.

**Figure 3.** Output predictions ( *N* = 5).

**Figure 4.** Step 2: Output prediction.

In order to solve the parameter matrix *Lw* and *Lu*, Equation (1) is rewritten into the Hankel matrix: 

$$\mathbf{Y\_{f}} = L\_{\text{d}\mathbb{P}} \cdot \mathbf{W\_{P}} + L\_{\text{u}} \cdot \mathbf{l} \\ \mathbf{I\_{f}} = \begin{bmatrix} \ \ L\_{\text{u}} & \ L\_{\text{u}} & \end{bmatrix} \begin{bmatrix} \ \mathbf{W\_{P}} \\ \ \mathbf{l} \mathbf{I\_{f}} \end{bmatrix} \tag{10}$$

where *Y*f, *<sup>W</sup>*p and *U*f are the Hankel matrix composed of input and output data, defined as follows:

$$\begin{aligned} \mathrm{^i\hat{Y}\_{\mathrm{f}}} = \left[ \begin{array}{c} \hat{Y}\_{\mathrm{f}\_{-\mathrm{f}}} \\ \vdots \\ \hat{Y}\_{\mathrm{f}\_{-\mathrm{m}}} \end{array} \right]; \mathrm{U}\_{\mathrm{f}} = \left[ \begin{array}{c} \boldsymbol{\mathcal{U}}\_{\mathrm{f}\_{-\mathrm{i}}} \\ \vdots \\ \boldsymbol{\mathcal{U}}\_{\mathrm{f}\_{-\mathrm{m}}} \end{array} \right]; \mathrm{Y}\_{\mathrm{P}} = \left[ \begin{array}{c} \boldsymbol{\mathcal{Y}}\_{\mathrm{P}\_{-\mathrm{i}}} \\ \vdots \\ \boldsymbol{\mathcal{Y}}\_{\mathrm{p}\_{-\mathrm{m}}} \end{array} \right]; \mathrm{U}\_{\mathrm{P}} = \left[ \begin{array}{c} \boldsymbol{\mathcal{U}}\_{\mathrm{P}\_{-\mathrm{i}}} \\ \vdots \\ \boldsymbol{\mathcal{U}}\_{\mathrm{p}\_{-\mathrm{m}}} \end{array} \right]; \mathrm{W}\_{\mathrm{P}} = \left[ \begin{array}{c} \boldsymbol{\mathcal{W}}\_{\mathrm{P}\_{-\mathrm{i}}} \\ \vdots \\ \boldsymbol{\mathcal{W}}\_{\mathrm{p}\_{-\mathrm{m}}} \end{array} \right] \end{aligned}$$

*<sup>W</sup>*p−*<sup>i</sup>* Δ= *U*Tp−*<sup>i</sup> <sup>Y</sup>*Tp−*<sup>i</sup>*<sup>T</sup> ∈ *R*2*N*×*<sup>M</sup> <sup>U</sup>*p−*<sup>i</sup>* Δ= ⎡ ⎢⎢⎢⎣ *ui*(0) *ui*(1) ··· *ui*(*<sup>M</sup>* − 1) *ui*(1) *ui*(2) ··· *ui*(*M*) . . . . . . ··· . . . *ui*(*<sup>N</sup>* − 1) *ui*(*N*) ··· *ui*(*<sup>N</sup>* + *M* − 2) ⎤ ⎥⎥⎥⎦; *U*f−*<sup>i</sup>* Δ= ⎡ ⎢⎢⎢⎣ *ui*(*N*) *ui*(*<sup>N</sup>* + 1) ··· *ui*(*<sup>N</sup>* + *M* − 1) *ui*(*<sup>N</sup>* + 1) *ui*(*<sup>N</sup>* + 2) ··· *ui*(*<sup>N</sup>* + *M*) . . . . . . ··· . . . *ui*(2*<sup>N</sup>* − 1) *ui*(2*<sup>N</sup>*) ··· *ui*(2*<sup>N</sup>* + *M* − 2) ⎤ ⎥⎥⎥⎦; *<sup>Y</sup>*p−*<sup>i</sup>* Δ= ⎡ ⎢⎢⎢⎣ *yi*(0) *yi*(1) ··· *yi*(*<sup>M</sup>* − 1) *yi*(1) *yi*(2) ··· *yi*(*M*) . . . . . . ··· . . . *yi*(*<sup>N</sup>* − 1) *yi*(*N*) ··· *yi*(*<sup>N</sup>* + *M* − 2) ⎤ ⎥⎥⎥⎦;*Y*ˆf−*i* Δ= ⎡ ⎢⎢⎢⎣ *y* ˆ *i*(*N*) *y*<sup>ˆ</sup>*i*(*<sup>N</sup>* + 1) ··· *y*<sup>ˆ</sup>*i*(*<sup>N</sup>* + *M* − 1) *y* ˆ *i*(*<sup>N</sup>* + 1) *y*<sup>ˆ</sup>*i*(*<sup>N</sup>* + 2) ··· *y*<sup>ˆ</sup>*i*(*<sup>N</sup>* + *M*) . . . . . . ··· . . . *y* ˆ *i*(2*<sup>N</sup>* − 1) *y*<sup>ˆ</sup>*i*(2*<sup>N</sup>*) ··· *y*<sup>ˆ</sup>*i*(2*<sup>N</sup>* + *M* − 2) ⎤ ⎥⎥⎥⎦

The problem is solved by least squares method:

$$\min\_{L\_{\text{uv}}, L\_{\text{uv}}} \left\lVert \mathcal{Y}\_{\text{f}} - \left[ \begin{array}{cc} L\_{\text{uv}} & L\_{\text{u}} \end{array} \right] \right\rVert \left\lVert \begin{array}{cc} \mathcal{W}\_{\text{P}} \\ \mathcal{U}\_{\text{f}} \end{array} \right\rVert \Big|\_{\text{F}}^{2} \tag{11}$$

This problem can be solved by orthogonal projection. According to the subspace projection theorem, *Y*f projects to the column space of *<sup>W</sup>*p and *U*f:

$$\hat{Y}\_{\rm f} = Y\_{\rm f} / \left[\begin{array}{c} \mathcal{W}\_{\rm P} \\ \mathcal{U}\_{\rm f} \end{array} \right] = Y\_{\rm f} \left[\begin{array}{c} \mathcal{W}\_{\rm P} \\ \mathcal{U}\_{\rm f} \end{array} \right]^{\rm T} \left( \left[\begin{array}{c} \mathcal{W}\_{\rm P} \\ \mathcal{U}\_{\rm f} \end{array} \right] \left[\begin{array}{c} \mathcal{W}\_{\rm P} \\ \mathcal{U}\_{\rm f} \end{array} \right]^{\rm T} \right)^{+} \left[\begin{array}{c} \mathcal{W}\_{\rm P} \\ \mathcal{U}\_{\rm f} \end{array} \right] \tag{12}$$

where the symbol "+" stands for pseudo-inverse; "/" stands for data space projection, then:

$$
\begin{bmatrix} L\_w & L\_w \end{bmatrix} = \mathbf{Y}\_\mathbf{f} \begin{bmatrix} \mathbf{W}\_\mathbf{P} \\ \mathbf{U}\_\mathbf{f} \end{bmatrix}^\mathrm{T} \left( \begin{bmatrix} \mathbf{W}\_\mathbf{P} \\ \mathbf{U}\_\mathbf{f} \end{bmatrix} \begin{bmatrix} \mathbf{W}\_\mathbf{P} \\ \mathbf{U}\_\mathbf{f} \end{bmatrix}^\mathrm{T} \right)^\dagger \tag{13}
$$

substituting Equation (13) into Equation (2), we obtain:

$$
\hat{Y}\_{\rm f} = \chi\_{\rm f} \begin{bmatrix} \ \mathcal{W}\_{\rm P} \\ \ \mathcal{U}\_{\rm f} \end{bmatrix}^{\dagger} \begin{bmatrix} \ \mathcal{W}\_{\rm P} \\ \ \mathcal{U}\_{\rm f} \end{bmatrix} \tag{14}
$$

QR decomposition of matrix ⎡⎣ *<sup>W</sup>*p *U*f *Y*f ⎤⎦ is shown as

$$
\begin{bmatrix} W\_{\text{P}} \\ \mathcal{U}\_{\text{f}} \\ \mathcal{Y}\_{\text{f}} \end{bmatrix} = \boldsymbol{R}^{\text{T}} \boldsymbol{Q}^{\text{T}} = \begin{bmatrix} \mathcal{R}\_{11} & \mathcal{0} & \mathcal{0} \\ \mathcal{R}\_{21} & \mathcal{R}\_{22} & \mathcal{0} \\ \mathcal{R}\_{31} & \mathcal{R}\_{32} & \mathcal{R}\_{33} \end{bmatrix} \boldsymbol{Q}^{\text{T}} \tag{15}
$$

then, Equation (12) can be rewritten as

$$\begin{aligned} \hat{Y}\_{\text{f}} &= \begin{bmatrix} R\_{31} & R\_{32} & R\_{33} \end{bmatrix} Q^{\text{T}} \left[ \begin{bmatrix} R\_{11} & 0 & 0 \\ R\_{21} & R\_{22} & 0 \end{bmatrix} Q^{\text{T}} \right]^{\dagger} \cdot \begin{bmatrix} W\_{\text{p}} \\ \mathcal{U}\_{\text{f}} \end{bmatrix} = \\ \begin{bmatrix} R\_{31} & R\_{32} & R\_{33} \end{bmatrix} Q^{\text{T}} Q^{\text{T}^{+}} \begin{bmatrix} R\_{11} & 0 & 0 \\ R\_{21} & R\_{22} & 0 \end{bmatrix}^{\dagger} \cdot \begin{bmatrix} W\_{\text{p}} \\ \mathcal{U}\_{\text{f}} \end{bmatrix} = \begin{bmatrix} R\_{31} & R\_{32} \end{bmatrix} \begin{bmatrix} R\_{11} & 0 \\ R\_{21} & R\_{22} \end{bmatrix}^{+} \begin{bmatrix} W\_{\text{p}} \\ \mathcal{U}\_{\text{f}} \end{bmatrix} \end{aligned} \tag{16}$$

comparing Equations (2) and (16), we obtain the solution of *Lw* and *Lu*:

$$
\begin{bmatrix} L\_{\varpi} & L\_{\Psi} \end{bmatrix} = \begin{bmatrix} R\_{31} & R\_{32} \end{bmatrix} \begin{bmatrix} R\_{11} & 0 \\ R\_{21} & R\_{22} \end{bmatrix}^{+} \tag{17}
$$

furthermore, we obtain:

$$L\_w = \begin{bmatrix} L\_w(1) \\ L\_w(2) \\ \vdots \\ L\_w(m) \end{bmatrix}; L\_u = \begin{bmatrix} L\_{u(1,1)} & L\_{u(1,2)} & \cdots & L\_{u(1,m)} \\ L\_{u(z,1)} & L\_{u(2,z)} & \cdots & L\_{u(z,m)} \\ \vdots & \vdots & \ddots & \vdots \\ L\_{u(m,1)} & L\_{u(m,2)} & \cdots & L\_{u(m,m)} \end{bmatrix} \tag{18}$$

where *Lw* is a 30 × 60 matrix, *Lu* is a 30 × 30 matrix, then after decomposition, we obtain *Lw*(*i*) and *Lu*(*<sup>i</sup>*,*j*):

$$L\_{w(i)} = L\_w(i:m:mN-m+i,:)$$

$$L\_{w(ij)} = L\_w(i:m:mN-m+i,j:m:mN-m+j)$$

where *i* = 1, 2, ··· , *m*, *j* = 1, 2, ··· , *m*, *m* = 6, *N* = 5. The prediction model of six subsystems of an aluminum reduction cell is further obtained.

The Nash optimal method is used to design the distributed control algorithm. The definition of Nash optimal is as follows.

For a complex large system with m subsystems, if there is a vector solution that *u<sup>N</sup>* = *<sup>u</sup>N*1 , ··· , *uNi* , ··· , *uNm* satisfies the following inequalities for all subsystems *ui*(*<sup>i</sup>* = 1, 2, ··· , *m*) [35]:

$$\begin{array}{c} J\_i\left(u\_1^N, \dots, u\_i^N, \cdot, \cdot, u\_m^N\right) \lessapprox \\ J\_i\left(u\_i^N, \cdot, \cdot, u\_{i-1}^N, u\_i, u\_{i+1}^N, \cdot, \cdot, u\_m^N\right) \end{array} \tag{19}$$

then, the vector is *u<sup>N</sup>* = *<sup>u</sup>N*1 , ··· , *uNi* , ··· , *uNm* called the optimal Nash solution of the system. This solution optimizes the control performance of the entire large system, and all subsystems will not change this control decision.

As mentioned in the algorithm introduction above, in distributed predictive control, each subsystem must know the optimal control signal of other subsystems before solving its own optimal control signal, but each controller solves the optimal control signal at every moment simultaneously. In order to make all subsystems optimal at the same time, the iterative method is usually adopted. At each sampling moment, an iterative calculation is carried out to obtain the optimal control input signal of each subsystem at the sampling moment and determine whether it satisfies the Nash optimal solution. At the same time, the control signal is transmitted to other subsystems. When the iterative values of all subsystems meet the conditions, the iteration ends, so the global optimization of a complex large system is realized. The detailed steps of the iterative algorithm are as follows:


Flowchart is shown in Figure 5:

**Figure 5.** Flow chart of distributed predictive control algorithm.
