Research on Unsteady Inverse Heat Conduction Based on Dynamic Matrix Control

: For the unsteady multi-boundary inverse heat conduction problem, a real-time solution method for boundary heat ﬂux based on dynamic matrix control is proposed in the paper. The method solves the heat ﬂux at the boundary in real-time by measuring the temperature information at the measurement points of the heat transfer system. A two-dimensional direct heat conduction model of the heat transfer system is established in the paper, and is solved by the ﬁnite difference method to obtain the temperature information of the measurement points under any heat ﬂux boundary. Then, the correspondence between the heat ﬂux of boundary and the temperature information is presented by means of a step-response model. The regularization parameters are introduced into the method to improve the stability of the inversion process, and the effect of real-time inversion on the heat ﬂux of the boundary is achieved through rolling optimization. The experimental results show that the proposed method can achieve real-time inversion of the heat ﬂuxes of the two-dimensional boundary with good accuracy.


Introduction
Heat transfer is important physical information in key parts of many special industrial equipment and can determine the production process or equipment life. Due to the complexity of the equipment operating environment, precise heat transfer information and boundary conditions are often difficult to measure or are unknown. For example, the temperature of the external thermal insulation tiles of a spacecraft cannot be directly measured when they re-enter the atmosphere and the heat flux at the solid-liquid interface cannot be measured during the growth of silicon single crystals. This difficult-to-obtain information on heat transfer can have a significant impact on the quality of the product and the safety of the equipment. In this case, the temperature information of the thermal system is commonly used to provide an effective estimate of the unknown boundary thermal property parameters, which are also known as inverse heat conduction problem (IHCP). As methods for solving inverse heat conduction problems continue to improve and their application in industry becomes more and more in demand, the study of efficient, accurate, and robust methods for solving IHCP is a valuable research task.
In recent decades, the inverse heat conduction problem has attracted a great deal of scholarly attention. Zhou et al. studied the identification of aircraft surface heat flux and inversed the heat flux of an aircraft surface from the temperature information of some measuring points [1]. Han et al. investigated the local stresses in pipelines caused by stratification of heat flux in the pipeline system of a nuclear power plant, and inverted the thermal stresses of the pipeline based on known temperature information to ensure the structural integrity of the pipeline during service [2]. Research on inverse heat conduction problem is widespread in the fields of solar energy, manufacturing, power engineering, aerospace, etc. [3][4][5].
At present, research on inverse heat transfer algorithms is divided into two main categories. One is for global domain inversion algorithms and the other for sequential inversion algorithms. The global domain inversion algorithm is based on the full temperature information to invert the unknown boundary over the entire time domain. For example, Cheng et al. proposed a spatio-temporal radial polynomial configuration method for solving a two-dimensional heat conduction inverse problem [6]. Pacheco et al. proposed to optimize the regularization parameters based on the generalized cross-validation method and solve linear and non-linear heat conduction inverse problems by the Tikhonov regularization method [7]. Yu et al. used the conjugate gradient method (CGM) to estimate the heat flux absorbed by a brake disc in a braking system [8]. The sequential inversion algorithm estimates the boundary information for the next moment based on the results of the previous moment's inversion, rolling forward until the unknown boundary information for the entire time domain is obtained. For example, Dong et al. combined the Broyden combination method with SFSM to reconstruct the temperature field during the heating of billets [9]. The literature [10,11] applied the Kalman technique to the study of inverse heat conduction problems. In response to the low accuracy of the first-order linearized approximation of the Kalman technique, an unscented Kalman filter (UKF) was proposed and combined with unscented Rauch-Tung-Striebel smoother (URTSS) to be applied to non-linear inverse heat conduction problems. In addition, there are also scholars using a hybrid method of a particle swarm optimization algorithm, normal distribution method, and finite element method to solve the inverse heat conduction problem and estimate thermal conductivity [12]. Zhang et al. used a neural network approach to solve the heat transfer coefficient at the interface between the casting and the metal mould [13]. Deng et al. combined Hopfield neural networks with BP neural networks to solve the heat conduction inverse problem, identifying unknown boundary conditions [14]. Therefore, sequential inversion algorithms and full-domain inversion algorithms are the two main types of methods for solving heat conduction inverse problems, where sequential inversion algorithms invert the boundary information based on the measurement information prior to the current moment and keep forwarding the inverse time domain, and are, therefore, suitable for online monitoring scenarios. However, full-domain inversion requires measurement information over the entire time domain and is suitable for offline calculation scenarios.
In recent years, some scholars gradually introduced the idea of feedback control into the solution of heat conduction inverse problems. For example, Wan et al. used the principle of proportional-integral differential (PID) control in feedback control theory to establish an inverse method for unsteady heat conduction, transforming the inverse problem of unsteady heat conduction into a PID control problem [15]. After that, Wan et al. combined a single artificial neuron with a PID inverse algorithm to adjust PID parameter weights adaptively and established a single neuron adaptive PID (SNA-PID) inverse method to estimate thermal boundary conditions for poor adaptive ability [16]. The literature [17] combines decentralized fuzzy inference with sequential quadratic programming and proposes a hybrid optimization technique to solve the temperature reconstruction problem of the thermal field. The PID control algorithm has the advantages of simple structure and strong anti-interference ability, and can meet the requirements of real time [18,19]. However, this method is not effective when there are multiple thermal boundary conditions in the heat transfer system. In recent years, other methods have been proposed. For example, Ahn et al. used an efficient iterative hybrid parameter selection algorithm to obtain stable inverse solutions and address the ill-posed nature of the inverse heat conduction problem [20]. Mostajeran et al. proposed a deep neural network approach for solving the heat transfer problem [21]. Wang et al. proposed an algorithm to identify non-smooth heat fluxes in a one-dimensional heat transfer system based on minimizing the data-match term and penalty term alternatively [22].
In this paper, the inverse problem of multi-boundary heat conduction is solved based on the dynamic matrix control method. Firstly, a two-dimensional linear unsteady direct heat transfer model is established, and the model is solved by the finite difference method to obtain the temperature field information under the action of boundary heat flux. Then, based on the dynamic step-response of the heat transfer system, the prediction model is established and the temperature information in the future time is predicted. Thirdly, the rolling optimization model is established by the sensitivity coefficient of the measuring point. By combining it with the rolling optimization theory, the predicted heat flux value at the current moment is optimized based on the deviation between the predicted temperature and the measured temperature. This approach aims to simultaneously invert the heat flux at multiple boundaries in the time domain. Finally, the inversion results of the SFSM method are compared with those of the proposed method. The influence of the location, number of measuring points, and future time step on the inversion results is analyzed. Through this analysis, the effectiveness of the proposed method is demonstrated. In this paper, the concept of control is introduced to solve the heat conduction inverse problem using dynamic matrix control. Based on the mathematical model of predictive control and rolling optimization, a predictive model for the measured point temperature and a rolling optimization model for the boundary heat flux are established. The proposed approach provides an effective theory and method for solving the two-dimensional heat conduction inverse problem.

Heat Conduction Model
Consider a heat transfer system with multiple heat flux boundaries, as illustrated in Figure 1. In this paper, the inverse problem of multi-boundary heat conduction is solved based on the dynamic matrix control method. Firstly, a two-dimensional linear unsteady direct heat transfer model is established, and the model is solved by the finite difference method to obtain the temperature field information under the action of boundary heat flux. Then, based on the dynamic step-response of the heat transfer system, the prediction model is established and the temperature information in the future time is predicted. Thirdly, the rolling optimization model is established by the sensitivity coefficient of the measuring point. By combining it with the rolling optimization theory, the predicted heat flux value at the current moment is optimized based on the deviation between the predicted temperature and the measured temperature. This approach aims to simultaneously invert the heat flux at multiple boundaries in the time domain. Finally, the inversion results of the SFSM method are compared with those of the proposed method. The influence of the location, number of measuring points, and future time step on the inversion results is analyzed. Through this analysis, the effectiveness of the proposed method is demonstrated. In this paper, the concept of control is introduced to solve the heat conduction inverse problem using dynamic matrix control. Based on the mathematical model of predictive control and rolling optimization, a predictive model for the measured point temperature and a rolling optimization model for the boundary heat flux are established. The proposed approach provides an effective theory and method for solving the two-dimensional heat conduction inverse problem.

Heat Conduction Model
Consider a heat transfer system with multiple heat flux boundaries, as illustrated in Figure 1. In Figure 1, 1 ( ) and 2 ( ) represent the time varying heat flux in two directions, respectively. All other boundaries are adiabatic. The two-dimensional heat conduction equation is shown as follows: The boundary conditions are shown in Equations (2)- (5). In Figure 1, q 1 (t) and q 2 (t) represent the time varying heat flux in two directions, respectively. All other boundaries are adiabatic. The two-dimensional heat conduction equation is shown as follows: The boundary conditions are shown in Equations (2)-(5).
−λ T(r, h, t) The initial condition is where c is the specific heat capacity, ρ is the density, λ is the thermal conductivity, T(r, h, t) is the temperature at a point in time t, h is the height of the heat conduction systems, r is the radius of the heat conduction systems, and q(t) is the boundary heat flux at time t.

Finite Difference Method
In this paper, the finite difference method is used to solve the heat conduction model. This method replaces the continuous solution domain with a finite number of grid nodes, replacing the differential in the heat conduction equation with the difference quotient of the function values at the grid nodes, thus, discretizing the continuous equation. The time term is discretized in forward difference form and the spatial term is discretized in central difference form for the computational domain, as shown in Equations (7)- (10).
∂T(r, h, t) ∂r Disregarding the truncation error, the finite difference form of the two-dimensional heat conduction equation can be derived from the above equation as shown.
where ∆t denotes the time step, then the number of nodes on the time axis is n t = t ∆t , ∆h denotes the spatial step on the H axis, then the number of grid nodes on the H axis is n h = H/∆h, ∆r denotes the spatial step on the R axis, then the number of grid nodes on the R axis is n r = R ∆r , T(i, j, k) = T r i , h j , t k , T r i , h j , t k denotes the temperature at time t k , the R axis position is r i , and the H axis position is h j The temperature information at the current time can be derived from the temperature information at the previous time, allowing us to obtain the temperature distribution of the thermal conduction system under the influence of boundary heat flux. Since the two-dimensional heat conduction equation is obtained by spatial discretization in an explicit finite difference format, the time and space steps must satisfy the following conditions in order to obtain a stable and convergent solution.

Predictive Model of Temperature
A predictive model can predict future output information using historical input information and future input information of a system. In dynamic matrix control algorithms, the predictive model plays a crucial role. The adjustment of boundary heat flux estimates is achieved through rolling optimization. The predicted values of the measured point for the next p time steps, starting from the moment k, are obtained by the prediction model. Therefore, the accuracy of the prediction model has a large impact on the accuracy of the inversion results. The prediction model is obtained from the unit step-response of the heat transfer system. Assuming the heat conduction system with M temperature measurement points, N boundary heat flux, the predicted future time step is P and the current moment is k, the prediction model is shown in Equation (13).
where φ p m,n (m = 1, 2 · · · , M, n = 1, 2, · · · , N, p = 1, 2 · · · , P) represents the dynamic stepresponse coefficient of the measurement point, also called the sensitivity coefficient, which indicates the degree of influence of the boundary heat flux on the temperature of the measurement point. ∆q k n denotes the increment of the boundary heat flux at moment k. (T k m ) pre indicates the predicted value of the temperature at the measurement point when pre can be estimated from the predicted value of the temperature under the impact of the boundary heat flux increment ∆q k−1 at the previous moment, according to the prediction model, as shown in Equation (14).
Translate Equation (13) into the form of a control increment as shown in Equation (15).
Translate Equation (15) into the form of a vector equation.
where T pre is the predicted estimated temperature,T pre is the initial temperature value, A is the dynamic step response matrix, and q is the heat flux. The expressions are as follows: where 0 = [0] M×N , ∆φ k represents the incremental sensitivity factor and is expressed as shown in Equation (22).
The expression for q is shown in Equation (23).
The expression for q k is shown in Equation (24).
Based on the above analysis, the prediction model shown in Equation (16) was developed using the dynamic step-response of the system. The model predicts future temperature estimates based on the optimal boundary heat flux estimate. The predicted temperature approximates the measured temperature to obtain heat flux for the full time domain.

Establishing a Rolling Optimization Objective Function
Dynamic matrix control determines the control increments of the system through rolling optimization to ensure that the system achieves optimal control targets. Optimal heat flux is obtained by rolling optimization. At time k, the boundary heat fluxes q k , q k+1 , . . . . . . , q k+p−1 are obtained based on a rolling optimization objective function using the measured temperature and the predicted temperature. However, this is from taking only q k as the boundary heat flux at moment k. At moment k + 1, the p boundary heat fluxes q k+1 , q k+2 , . . . . . . ,q k+p in the optimized time domain t k+1 t k+p are obtained in the same way, repeated to obtain the boundary heat flux in the time domain. The rolling optimization objective function is as shown in Equation (25). minJ q k , q k+1 , · · · , q k+p−1 = where α is the regularization parameter. The regularization parameter ensures the stability of the inversion results [23]. The rolling optimization objective function shown in Equation (25) is converted to a vector Equation (26).
is the temperature vector of measurement point m in the optimized time domain.
The above analysis was combined to establish a rolling optimization objective function for the heat conduction inverse problem. The objective function takes into account the square of the deviation between the predicted temperature and the measured temperature, as well as the stability of the inverse heat flux. The optimal estimate of the boundary heat flux is obtained by solving the extreme value of the objective function, dJ(Q) dQ = 0. The results are shown as: where E is a unit matrix of order p. At time k, only the best estimate at the current moment is taken. Then the optimal estimate at moment k is shown as: where Φ = [1, 0, · · · , 0] is the calculation taking the first component of the timing vector.

Solving for the Measurement Point Sensitivity Coefficient
The sensitivity coefficient is an important parameter in the solution algorithm proposed in this paper. Sensitivity coefficients are required for both rolling optimization and prediction models. The mathematical model of the sensitivity coefficient is shown as follows: where H n (r, h, t) represents the degree of influence of the boundary heat flux on the temperature at the measurement point, also known as the sensitivity coefficient. The sensitivity coefficient satisfies the dynamic step-response equation as shown in Equation (30) within time t k t k+p−1 .
The constraints are as follows: The dynamic step-response equations shown in Equations (30)-(34) are solved using the finite difference method in t k t k+p−1 . Forward differencing is used for the time term, and central differencing is used for the spatial term. The dynamic step-response coefficients of the temperature measurement points at different locations for each boundary heat flux in t k t k+p−1 are obtained by solving the problem, as shown in Equation (35).

Solving for Optimal Regularization Parameters by the Residual Principle
The regularization parameter is in the rolling optimization objective function. By introducing the regularization parameter, it ensures that the heat flux estimates are as smooth as possible and improves the stability of the inverse heat flux. In a real production environment, temperature measurement equipment is inevitably affected by noise. As a result, the measured temperature value is the result of the superposition of the true temperature value and the measurement error, as shown in Equation (37).
where ω is a set of random numbers that follow a standard normal distribution over the interval [−2.576, 2.576]. σ is the standard deviation of the temperature measurement error.
The residuals of the boundary heat fluxes in the time domain are shown as follows: where (q k ) exa and (q k ) pre are the exact and predicted values of the boundary heat flux, respectively. The inverse boundary heat flux is applied to the system to obtain the predicted temperature at the measurement point, and the temperature residuals of the inverse results are calculated by Equation (40).
Ideally, the predicted temperatures and exact temperatures should be equal, T pre m = T exa m . The optimal regularization parameter is: Therefore, different regularization parameters can be chosen to obtain inverse residuals in practical applications. The optimal regularization parameter is selected according to Equation (41).

Steps of Inverse Heat Transfer Algorithm Based on Dynamic Matrix Control
In Figure 2, b represents the number of iterations, q represents estimated heat flux, η represents the threshold, k represents the current moment, and K represents the full inversion moment. Based on the above analysis, the dynamic matrix control inversion of the boundary heat fluxes proceeds as follows: Energies 2023, 16 (14) The optimal estimate is obtained by rolling optimization according to Equation (28) Solve for the sensitivity factor according to Equation ( Step 1: Set the thermal property parameters for the heat conduction equation; Step 2: Solve the dynamic step response equation to obtain the sensitivity coefficient matrix; Step 3: The optimal estimate of the prediction boundary heat flux is obtained by rolling optimization; Step 4: According to the prediction model, obtain the predicted value of the temperature of the measurement point under the predicted heat flux; Step 5: Determine whether convergence is achieved. If convergence is achieved, the predicted value at the current time is taken as the optimal estimate, and the predicted value of the boundary heat flux at other times in the future time step is taken as the initial heat flux value at the next time. If the convergence criterion is not reached, return to step 3 until the convergence criterion is reached.

Verification of Algorithm Validity
To verify the effectiveness of inverse heat transfer algorithms based on dynamic matrix control, we consider the inverse heat transfer problem with multiple heat flux boundaries, as shown in Figure 1  Step 1: Set the thermal property parameters for the heat conduction equation; Step 2: Solve the dynamic step response equation to obtain the sensitivity coefficient matrix; Step 3: The optimal estimate of the prediction boundary heat flux is obtained by rolling optimization; Step 4: According to the prediction model, obtain the predicted value of the temperature of the measurement point under the predicted heat flux; Step 5: Determine whether convergence is achieved. If convergence is achieved, the predicted value at the current time is taken as the optimal estimate, and the predicted value of the boundary heat flux at other times in the future time step is taken as the initial heat flux value at the next time. If the convergence criterion is not reached, return to step 3 until the convergence criterion is reached.

Verification of Algorithm Validity
To verify the effectiveness of inverse heat transfer algorithms based on dynamic matrix control, we consider the inverse heat transfer problem with multiple heat flux boundaries, as shown in Figure 1. Where L x = 0.4 m, L y = 1 m. The number of grids is 100 × 40. The thermal conductivity λ = 40/m · k, Density ρ = 2700 kg/m 3 . Three different forms of boundary heat flux are applied into the heat transfer system. The temperature variation of the measuring point is obtained by solving the heat transfer equations, which serve as the input. The inversion results of boundary heat flux are compared with the actual boundary heat flux to verify their effectiveness.
The first boundary heat flux is the exponential form: The second boundary heat flux is the sinusoidal form: The third boundary heat flux is the square form: The three different forms of boundary heat flux mentioned above are applied to the heat conduction system to obtain the temperature changes at two measurement points, denoted as Tc1, Tc2. r Tc1 = 0.2 m, h Tc1 = 0.1 m, r Tc2 = 0.3 m, h Tc2 = 0.4 m. The temperature changes at measurement points Tc1, Tc2 are shown in Figure 3.
The temperature values are used as input to invert the boundary heat flux q 1 (t) and q 2 (t). The results are compared with those obtained by the sequential function specification method (SFSM) to verify the effectiveness of the algorithm [24].
In practical engineering, the production environment often introduces unavoidable disturbances that can lead to errors in the measuring equipment. These errors are superimposed on the measured object's temperature, resulting in deviations between the true and measured temperatures. Therefore, to enhance the reliability of the experiment, it is important to simulate the measurement errors that occur in real engineering scenarios.
The mathematical model is shown as where σ is the standard deviation of the measurement error and ω is a standard normal distributed random number with a confidence coefficient of 99% within the confidence interval [−2.576, 2.576].
In order to express the accuracy and effectiveness of the proposed inverse heat transfer algorithm, the RMS error of heat flux and the RMS error of temperature are defined in Equations (49) and (50), respectively: where, S q represents the RMS error of the estimated boundary heat flux and S T represents the RMS error of the measured point temperature.   As shown in Figure 4, the dynamic matrix-controlled inverse solution method of heat transfer proposed in the paper can accurately inverse the three heat flux q1 and q2. Compared with the SFSM method, the algorithm proposed in this paper exhibits smaller thermal boundary inversion errors. In the case of square, the inversion results obtained using SFSM are oscillatory and fail to accurately capture the boundary heat flux when there are significant changes in the actual heat flux. The algorithm proposed in this paper maintains the accuracy of the inversion results at times of large changes in the heat flux. In sinusoidal form, the heat flux starts to decrease at t = 90 s and the SFSM algorithm fails to track the real heat flux in a timely manner, resulting in a significant increase in error. The DMC algorithm is able to follow the boundary heat flux accurately as the heat flux decreases. A comparison of the inversion errors of the boundary heat flux using the two different algorithms is shown in Table 1. The table shows that the inverse error of the boundary heat flux obtained using dynamic matrix control is lower than that of the SFSM method for the same heat flux boundary. This indicates the effectiveness of the algorithm proposed in this paper.

Measurement Point Location Effects on Inversion Results
Two different temperature measurement points are selected in the discretized grid to invert the boundary heat flux, and the effect of the measurement point locations on the inversion results is investigated. In the first group, measurement point one is located at (97,20) in the discretized grid, and measurement point two is located at (87, 37). In the second group, measurement point one is located at (92, 20), and measurement point two is located at (87, 33). The time step r = 3, the maximum measurement error ∆T max = 3.
As can be seen from Figure 5, the location of the measurement points can influence the inversion results. The inversion errors from the first set of measurement points are smaller than the second set. The reason for this is that the first group of measurement points are placed closer to the thermal boundary than the second group. Consequently, the measurement points in the first group are more sensitive to the changes in the boundary heat flux and can better track the variations. Therefore, careful consideration should be given to the selection of measurement point locations in practical applications.

Influence of the Number of Measurement Points on the Inversion Results
In this experiment, the influence of the number of measurement points on the inversion results is investigated by inverting the boundary heat flux using three different numbers of measurement points. The number of temperature measurement points taken is M = 2, M = 4, and M = 6. The time step r = 3, the maximum measurement error ∆T max = 5. The results of the inversion are shown in Figure 6. As can be seen from Figure 6, when the number of measurement points M = 6, the inversion error is significantly smaller than when the number of measurement points M = 2. When the number of measurement points M = 4, the error in the inversion of the three different forms of boundary heat flux ranges between M = 2 and M = 6. Therefore, under the same conditions, the higher the number of temperature measurement points, the smaller the inversion error of the boundary heat flux. With the increase in the number of measuring points, more information about temperature distribution is included in the system, which can more accurately indicate the deviation between the predicted temperature and the measured temperature. At the same time, increasing the number of measuring points can reduce the error and improve the accuracy of the inversion. However, it is found in the experiment that as the number of measurement points increases, the algorithm execution time increases accordingly, which reduces the real-time capability of the algorithm. Therefore, in practical applications, increasing the number of temperature measurement points without affecting the real-time performance of the algorithm can effectively improve the inversion accuracy of the boundary heat flux.

Influence of Future Time Steps on Inversion Results
In this experiment, to research the influence of future time steps on the inversion results, we set r = 4, r = 5 and r = 6. We set the number of temperature measurement points M = 2 and the maximum measurement error ∆T max . The inversion results are shown in Figure 7. As can be seen from Figure 7, the inversion error is higher when taking time step r = 4 than when taking time step r = 5 and r = 6. For time step r = 6, the inversion error is the lowest and the inverse boundary heat flux is the closest to the actual value. It can be seen that the longer the time step, the more temperature information the prediction model has about the future time step, and the higher the inversion accuracy. Therefore, an appropriate increase in the time step can improve the inversion accuracy. In the experiments, it is found that when the boundary heat flux changes more drastically, the increase in time step lead to oscillations in the iterative process, resulting in unstable inversion results. Therefore, in the process of inversing the multi-boundary heat flux using dynamic matrix control, if the future time step is too short, the inversion accuracy is reduced. On the other hand, if the future time step is too large, although the inversion accuracy can be improved, it reduces the robustness of the inversion process and increases the computational time of the algorithm.

Conclusions
In this paper, a real-time inverse algorithm based on dynamic matrix control is proposed for the inverse heat conduction problem of multi-boundary heat flux. Firstly, a two-dimensional direct heat conduction model is established, and solved by finite difference method. Combined with the dynamic matrix control principle, a prediction model is established. Rolling optimization is used to invert the boundary heat flux in real-time, and a regularization parameter is introduced to improve the stability of the inversion process. The influence of measurement point location, the number of measurement points, and future time step on the inversion results is analyzed. The experimental results show that the higher of measurement points and the longer the future time step, the higher the accuracy of the boundary heat flux inversion. However, having too many measurement points and a too long time step can affect the real-time capability and stability of the algorithm. When the number of measurement points and the time step length are not chosen properly, both can lead to algorithm divergence. In the future, we will study the heat transfer inverse problem based on a three-dimensional model, and consider various heat transfer methods. Additionally, we will explore the application of a fuzzy control approach to identify unknown boundary heat fluxes.

Conflicts of Interest:
The authors declare no conflict of interest.