Nonlinear multivariable model-free adaptive predictive control method for refrigeration station

Abstract

The invention relates to the technical field of energy-saving control. The objective of the invention is to provide a nonlinear multivariable model-free adaptive prediction control method for a refrigeration station. Comprising the following steps: S1, analyzing the energy efficiency of each loop of a refrigeration station, and determining a to-be-optimized control variable; s2, constructing a control system; and S3, executing control. According to the method, a refrigeration station system tight format dynamic linearization model is constructed based on a data driving method, and state online estimation is carried out through system input and output data by using a pseudo-Jacobian matrix in the model; then, the dynamic relation between each control variable of the refrigeration station and the energy efficiency of the refrigeration unit is deeply analyzed, an energy-saving target control item is introduced into an optimization target, so that the system runs in a high-energy-efficiency state, and a time-varying proportion-integral control item is introduced to compensate information loss in the linearization process; meanwhile, a multi-step prediction rolling optimization mechanism is designed to overcome the influence of large inertia of the system; in addition, a grey wolf optimization algorithm is utilized, a fitness function is designed and improved, and controller parameters are set; finally, the convergence and stability of the method are proved theoretically.

Description

Technical Field

[0001] The present invention relates to the technical field of energy-saving control, and in particular to a nonlinear multivariable model-free adaptive predictive control method for a refrigeration station. Background Art

[0002] Refrigeration station is a key system that provides cold source for industry, commerce or large buildings. Its structure is as follows Figure 1 As shown, it mainly consists of three parts, namely the refrigeration unit, the chilled water circuit and the cooling water circuit. The refrigeration unit is considered the core of the refrigeration station. The liquid refrigerant enters the evaporator through the expansion valve, absorbs the heat of the chilled water and evaporates into vapor. It then passes through the compressor and becomes high-temperature and high-pressure steam before entering the condenser, transferring its own heat to the cooling water and liquefying it, thus completing the heat transfer cycle of the refrigeration station. The chilled water circuit mainly consists of a chilled water pump, chilled water pipes and an air handling unit. The chilled water flows through the air handling unit to remove the heat from the indoor air and exchanges the heat with the refrigerant in the evaporator part of the refrigeration unit, thus completing the chilled water absorption and release cycle. The cooling water circuit mainly consists of a cooling tower, a cooling water pump and cooling water pipes. When the cooling water flows through the cooling tower, it exchanges heat with the air under the action of the high-speed rotation of the fan at the top of the cooling tower, dissipating the heat into the atmosphere and removing the heat from the refrigerant vapor in the condenser of the refrigeration unit, thus completing the cooling water absorption and release cycle.

[0003] Currently, PID (proportional, integral, and breeze) control is widely used in refrigeration plant control due to its simple structure, fast response, and ease of implementation. Improved PID control methods, such as nonlinear PID controllers and hybrid PID controllers, have also been applied in refrigeration plant systems. However, for complex nonlinear systems, single-loop local PID control is not only difficult to tune but also lacks the ability to improve the overall performance of the entire system. Furthermore, the controller's robustness is limited in the presence of external environmental disturbances. To improve system dynamic performance and enhance energy efficiency, a small number of researchers have proposed model-based refrigeration plant control methods, including model predictive control, optimal control, linear quadratic Gaussian control, and decoupled multivariable control. However, these control methods are highly model-dependent, making it difficult to accurately establish mathematical models for complex nonlinear systems such as refrigeration plants. Machine learning-based model construction methods rely heavily on the quality of training data. On the one hand, they require data that covers the complete dynamic characteristics of the system, making implementation complex. On the other hand, the control performance implicit in the data itself often fails to achieve the desired energy-saving effect, limiting the effectiveness of the control algorithm. Furthermore, control performance cannot be guaranteed when system equipment performance degrades or when external environmental parameters change significantly compared to historical data.

[0004] Model-free adaptive control (MFAC) is a data-driven adaptive control method proposed by Hou Zhongsheng et al. It does not rely on a system model and adaptively changes the control strategy solely through input and output data, enabling real-time adaptation to system dynamics. Because this method does not require prior knowledge of the system's precise internal structure or dynamic characteristics, it has low computational complexity and simplifies the controller design process. MFAC and its offshoots have been successfully applied in a variety of fields, such as density balance control of highway and auxiliary road systems, nitrate and dissolved oxygen concentration treatment in wastewater, low-speed and tracking control of motor systems, dynamic event-triggered safety tracking control of subway systems, and trajectory tracking control of continuum manipulators, demonstrating the effectiveness of MFAC in complex system control. In recent years, some scholars have conducted related research in the field of central air-conditioning control. Chen et al. applied MFAC to the control of a vapor-compression refrigeration system. They designed an optimization objective function with the goal of minimizing tracking error and preventing excessive input changes, and used a BP neural network to self-tune the MFAC parameters. Jia et al. introduced a control rate step size factor into MFAC to improve the heating rate and energy efficiency of air conditioning heat pump systems, optimizing algorithm stability and convergence speed. Jiang et al. combined a first-principles model with a full-format dynamic linearized data-driven model of MFAC for a multi-zone HVAC system, constructing a system hybrid model to describe the multi-zone climate dynamics and designing an MFAC scheme for multi-zone climate control. Feng et al. addressed the time-delay uncertainty of air conditioning systems caused by equipment aging, nonlinearity, and external factors. They designed an MFAC control method to adapt to the changing characteristics of HVAC systems and studied MFAC parameter optimization methods. Their results further clarified the scope of MFAC's applicability. However, these studies primarily focused on the output error and input variation of air conditioning systems, failing to fully consider external interference factors and the need to improve system energy efficiency. Consequently, some researchers have improved traditional MFAC schemes. To address the complex nonlinearities and strong coupling characteristics of absorption refrigeration systems, Dong et al. incorporated an output error rate into the MFAC objective function and introduced time-varying proportional and integral control terms to improve the system's tracking speed and denoising capabilities. However, the sum of the output error coefficient and the coefficient of the output error differential was fixed at 1, resulting in insufficient algorithm adaptability. Wang et al. proposed a partial MFAC method to address the room temperature control problem in agricultural greenhouses. By introducing an energy penalty term and optimizing sensitive parameters, they achieved an energy saving rate of 12.35%. However, the room temperature control system only involves single-loop control, while the refrigeration station system involves the internal circuit of the chiller, the chilled water circuit, and the cooling water circuit. The circuits are highly coupled, and the controlled variables and control quantities of each circuit have a significant impact on the system's energy efficiency. Therefore, the design of an optimization objective function oriented towards energy efficiency improvement is relatively difficult. Summary of the Invention

[0005] The purpose of the present invention is to provide a nonlinear multivariable model-free adaptive predictive control method for a refrigeration station, the control target of which is to improve the overall energy efficiency of the system while meeting the building load requirements.

[0006] To achieve the above-mentioned purpose, the technical solution adopted by the present invention is: a nonlinear multivariable model-free adaptive predictive control method for a refrigeration station; comprising the following steps:

[0007] S1. Analyze the energy efficiency of each circuit in the refrigeration station and determine the control variables to be optimized as chilled water supply temperature, cooling water flow rate, and chilled water flow rate;

[0008] S2. Build a control system: The control system adopts a structure that combines an upper optimization layer and a field control layer. The upper optimization layer includes a model-free adaptive predictive controller, PJM prediction, PJM estimation, optimization objectives, time-varying proportional control terms, and time-varying integral control terms. The upper optimization layer can calculate the optimal set values ​​of the control variables of each loop. The field control layer can use the PID method to adjust the operation of the field actuator equipment so that the chilled water supply temperature, cooling water flow rate, and chilled water flow rate follow the optimal set values ​​of the control variables calculated by the upper optimization layer.

[0009] S3. Execution control.

[0010] Preferably, in step S2, when the upper optimization layer calculates the optimal setting value of each loop control variable, it includes:

[0011] S2.1: Design the optimization objective function, the optimization objective function is as shown in formula (3);

[0012]

[0013] Where λ>0 is the weight factor, ζ>0 is the energy saving factor; y * (k+i)∈R m×1 is the expected value of the controlled variable of the system at time k+i, y(k+i)∈R m×1 is the actual value of the controlled variable at time k+i, i=1,2,…N, N is the system output prediction time domain; u * (k+r) is the expected value of each control variable of the system at time k+r, u(k+r) is the actual value of each control variable of the system at time k+r, r=1,2,…N u , N u is the system control time domain; where Δu(k+j)=u(k+j)-u(k+j-1);

[0014] S2.2: Use the MFAPC method with energy-saving control items to calculate the optimal set values ​​of the control variables of each loop:

[0015] Consider the following discrete-time MIMO system:

[0016] y(k+1)=f(y(k),…,y(kn y ),u(k),…,u(kn u )) (4)

[0017] Where u(k),y(k)∈R m , represents the control quantity and controlled variable of the system at time k; n u and n y is an unknown positive integer; is an unknown nonlinear function;

[0018] The following nonlinear dynamics assumptions are given:

[0019] Assumption 1: f i The partial derivative of (…)(i=1,…,m) with respect to the controlled variable u(k) is continuous; this is the basic condition for controller design;

[0020] Assumption 2: System equation (4) satisfies the generalized Lipschitz condition, indicating that the relationship between the controlled variables and the control variables of the refrigeration station system is constrained and the rate of change is limited; as shown in equation (5):

[0021] ‖y(k1+1)-y(k2+1)‖≤q‖u(k1)-u(k2)‖ (5)

[0022] Among them, for any time k1≠k2, k1, k2≥0, and u(k1)≠u(k2); q>0 is a constant;

[0023] Lemma 1: For the refrigeration station system (4) that satisfies Assumptions 1 and 2, there exists a pseudo-Jacobian matrix Φ(k) that can be transformed into the following compact form dynamic linearized data model:

[0024] y(k+1)=y(k)+Φ(k)Δu(k) (6)

[0025]

[0026] And for any time k, Φ(k) is bounded;

[0027] In order to facilitate the algorithm stability analysis later, the following assumptions are made:

[0028] Assumption 3: PJM Φ(k) satisfies the diagonal advantage condition, i.e. |φ ij (k)| <s1,s2≤|φ ii (k)|≤αs2, i=1,2,…,m,j=1,2,…,m,i≠j, and the signs of the elements in Φ(k) remain unchanged; s1 and s2 are two positive numbers, and s2=βs1, where β is a positive number;

[0029] The N-step forward prediction value of the system controlled variable can be obtained from formula (6):

[0030]

[0031] Considering j>N u When Δu(k+j-1)=0, the control variable increment coefficient of the controlled variable prediction equation is obtained as:

[0032]

[0033] The N-step prediction equation of the system controlled variable can be written as:

[0034]

[0035] Among them, Y N (k+1)=[y T (k+1),y T (k+2),…,y T (k+N)] T ∈R Nm×1 ;E y (k)=[I m×m ,…,I m×m ] T ∈R Nm×m ;

[0036] because

[0037] u(k)=u(k-1)+Δu(k) (11)

[0038] according to And formula (11) can be used to get the future N u The optimal control variable for the step:

[0039]

[0040] Then the incremental coefficient of the control variable prediction equation is:

[0041]

[0042] The N of the system control variables can be obtained u Step prediction equation:

[0043]

[0044] in

[0045] make According to formulas (10) and (14), the optimization objective function (3) can be rewritten as:

[0046]

[0047] Substitute formula (10) and formula (14) into formula (15) and use the optimization condition Under the conditions of satisfying Assumptions 1, 2 and Lemma 1, the following control law can be obtained:

[0048]

[0049] Observe formula (16), where the item that plays a role in tracking the expected value of the controlled variable is Define it as:

[0050]

[0051] We can get:

[0052]

[0053] S2.3: S4, use the Grey Wolf optimization algorithm to perform parameter tuning;

[0054] Preferably, in step S2.2, the incremental PID method is combined with the MFAPC method, and additional tracking error information is introduced to compensate for the reduction in system information accuracy caused by over-linearization, thereby obtaining the ICF-MFAPC method, and the control law is designed as follows:

[0055]

[0056] Among them, kp≥0, ki>0 are weighting coefficients;

[0057] Formula (19) is simplified to:

[0058]

[0059] Then the control quantity can be obtained:

[0060]

[0061] in

[0062] Preferably, in step 2.3, the fitness function of the gray wolf optimization algorithm is as follows:

[0063]

[0064] Where fitness>0; Represents the expected EER during the optimization process, which is generally a constant. y1(k+1) represents the EER at time k+1.

[0065] Preferably, the control strategy in step 3 execution control includes:

[0066] S3.1 Pseudo-Jacobian matrix estimation;

[0067] S3.2 Pseudo-Jacobian matrix prediction;

[0068] S3.3 Model-free adaptive control law update.

[0069] Preferably, the Jacobian matrix estimation in S3.1 includes the following algorithm:

[0070]

[0071] in, is the estimated value of PJM; η∈(0,2] is the step size factor of PJM update, which can enhance the flexibility of the algorithm, μ>0; and:

[0072]

[0073] To improve the tracking ability of the algorithm and satisfy Assumption 3, PJM estimates that the following provisions need to be made:

[0074] like or or

[0075] like or

[0076] Preferably, the Jacobian matrix prediction in S3.2 includes the following algorithm:

[0077]

[0078] Where θ i (k)∈R m×m ,i=1,…,N p is the order of the multi-layer hierarchical forecasting method; j = 1,…,N u -1;

[0079] θ(k) update:

[0080]

[0081] Where δ∈(0,1]; and must satisfy: θ(k)=θ(0), if ‖θ(k)‖≥M θ , M θ is a positive number;

[0082] Preferably, the model-free adaptive control law update in S3.3 includes the following algorithm:

[0083]

[0084]

[0085] According to (23), (26) and (27), the parameters that need to be optimized for GWO are obtained. is also unknown, so (22) can be written as follows:

[0086]

[0087] The beneficial effects of the present invention are concentrated in the following aspects: the method adopts a tight dynamic linearization method to construct an input / output data-driven model of a refrigeration station; the method deeply analyzes the dynamic relationship between the flow rate and temperature of chilled water and cooling water and the energy efficiency of the refrigeration unit, and determines the energy efficiency improvement optimization target item; at the same time, the output error rate is added to the control rate update function, and time-varying proportional control items and time-varying integral control items are introduced into the optimization objective function to obtain additional tracking error information, compensate for the data-driven model error caused by linearization, and improve the dynamic performance and robustness of the system under different load conditions; at the same time, the prediction rolling optimization mechanism in the model predictive control method is introduced to overcome the influence of large lag and large inertia characteristics of the refrigeration station system, and further improve the dynamic performance of the system; on this basis, the gray wolf optimization algorithm is introduced, and the fitness function is improved according to the control objective design, and the controller parameters are adjusted. BRIEF DESCRIPTION OF THE DRAWINGS

[0088] Figure 1 This is the schematic diagram of the refrigeration station system;

[0089] Figure 2 This is a graph showing the effect of chilled water supply temperature on the COP of the refrigeration unit;

[0090] Figure 3 This is a graph showing the effect of chilled water flow on the COP of the refrigeration unit;

[0091] Figure 4 This is a graph showing the effect of cooling water flow on the COP of the refrigeration unit;

[0092] Figure 5 This is the structure diagram of the model-free adaptive predictive control for the refrigeration station;

[0093] Figure 6 This is the model-free adaptive predictive control flow chart for the refrigeration station system;

[0094] Figure 7 This is a schematic diagram of the neural network model structure of the refrigeration station;

[0095] Figure 8 Statistical chart of EER prediction for neural network modeling;

[0096] Figure 9 Develop a statistical chart for predicting cooling capacity using neural network modeling;

[0097] Figure 10 This is a statistical chart showing the improvement of EER by the ICF-MFAPC algorithm;

[0098] Figure 11 It is the ICF-MFAPC algorithm continuous follow-up cooling capacity statistics chart;

[0099] Figure 12 This is the EER improvement diagram;

[0100] Figure 13 It is the cooling capacity following situation diagram;

[0101] Figure 14 Chilled water supply temperature result graph;

[0102] Figure 15 This is the chilled water flow result diagram;

[0103] Figure 16 Cooling water flow results graph. DETAILED DESCRIPTION

[0104] The present invention will be described in detail below in conjunction with the structural form and specific applications of the refrigeration station system.

[0105] 1. Energy efficiency analysis of refrigeration station system

[0106] Refrigeration station is a key system that provides cold source for industry, commerce or large buildings. Its structure is as follows: Figure 1As shown, it mainly consists of three parts, namely the refrigeration unit, the chilled water circuit and the cooling water circuit. The refrigeration unit is considered the core of the refrigeration station. The liquid refrigerant enters the evaporator through the expansion valve, absorbs the heat of the chilled water and evaporates into vapor. It then passes through the compressor and becomes high-temperature and high-pressure steam before entering the condenser, transferring its own heat to the cooling water and liquefying it, thus completing the heat transfer cycle of the refrigeration station. The chilled water circuit mainly consists of a chilled water pump, chilled water pipes and an air handling unit. The chilled water flows through the air handling unit to remove the heat from the indoor air and exchanges the heat with the refrigerant in the evaporator part of the refrigeration unit, thus completing the chilled water absorption and release cycle. The cooling water circuit mainly consists of a cooling tower, a cooling water pump and cooling water pipes. When the cooling water flows through the cooling tower, it exchanges heat with the air under the action of the high-speed rotation of the fan at the top of the cooling tower, dissipating the heat into the atmosphere and removing the heat from the refrigerant vapor in the condenser of the refrigeration unit, thus completing the cooling water absorption and release cycle.

[0107] Energy efficiency analysis of each circuit of the refrigeration station:

[0108] Since the energy consumption of the refrigeration unit accounts for more than half of the energy consumption of the refrigeration station system, improving the coefficient of performance (COP) of the refrigeration unit can effectively improve the overall energy efficiency ratio (EER) of the refrigeration station system. COP is the ratio of the cooling capacity of the refrigeration unit to the input power per unit time, and EER is the ratio of the cooling capacity of the refrigeration system to the input power per unit time:

[0109]

[0110]

[0111] Among them, Q ch is the cooling capacity; P chill P is the operating power of the refrigeration unit; pc and P pc are the operating power of the chilled water pump and the cooling water pump respectively; P tower The operating power of the cooling tower.

[0112] The main variables of a refrigeration unit include the chilled water supply temperature, chilled water flow rate, cooling water flow rate, and cooling water supply temperature (i.e., the temperature of the cooling water flowing back to the refrigeration unit from the cooling tower). The lower the cooling water supply temperature, the higher the heat exchange efficiency of the refrigerant in the condenser, and therefore the higher the COP of the refrigeration unit. The cooling water supply temperature is primarily affected by the outdoor wet-bulb temperature, the cooling tower fan speed, and the cooling tower structure. Its minimum is the outdoor wet-bulb temperature + the cooling tower approximation. The time constant between the cooling water supply temperature and the cooling tower fan speed is relatively short. Therefore, the present invention sets the cooling water supply temperature setpoint to the outdoor wet-bulb temperature + 3°C to 5°C approximation, and sets a separate PID control loop for the cooling tower fan for regulation, without incorporating it into the ICF-MFAPC control. Therefore, the optimization variables of the present invention are determined to be the chilled water supply temperature, cooling water flow rate, and chilled water flow rate. To improve the COP of the refrigeration unit and thus the EER of the refrigeration station, it is necessary to analyze the impact of these three variables on the COP of the refrigeration unit.

[0113] Figures 2 to 4 This is the relationship curve between each variable and the COP of the refrigeration unit obtained through experimental means. Figure 2 It can be seen that when the cooling water supply temperature, chilled water flow rate and cooling water flow rate are constant, the chilled water supply temperature increases, and the COP of the refrigeration unit increases accordingly. Moreover, when the chilled water supply temperature increases by 1°C, the COP of the refrigeration unit will increase by 1.5% to 3%. This is because when the evaporation temperature increases, the compression ratio of the compressor decreases, and the power consumption of the compressor decreases; and increasing the evaporation temperature will reduce the specific volume of the refrigerant vapor and increase the cooling capacity per unit volume. Therefore, when the condensing temperature is constant, increasing the chilled water supply temperature can improve the COP of the refrigeration unit. Figure 3 As shown in the figure, when the chilled water supply temperature, cooling water supply temperature and cooling water flow rate are constant, the three chilled water variable flow rate curves are very close, causing the COP to drop within 2% under low load conditions; under different chilled water supply and return water temperature differences, the constant flow rate COP drops within 5%, which shows that the COP change of the refrigeration unit caused by the chilled water variable flow rate is small. The reason is that on the one hand, when the chilled water flow rate decreases, the heat exchange efficiency on the evaporator side decreases; on the other hand, when the chilled water supply temperature is constant, the change in the chilled water return temperature will affect the evaporation temperature; under partial load, the average water temperature of the evaporator of the variable flow system is higher than that of the constant flow system, and the evaporation temperature is relatively high, thereby improving the COP of the unit. These two aspects offset each other, resulting in a small impact of the chilled water variable flow rate on the COP of the refrigeration unit. As shown in the figure, when the chilled water flow rate decreases, the heat exchange efficiency on the evaporator side decreases, and on the other hand, when the chilled water supply temperature is constant, the change in the chilled water return temperature will affect the evaporation temperature; under partial load, the average water temperature of the evaporator of the variable flow system is higher than that of the constant flow system, and the evaporation temperature is relatively high, thereby improving the COP of the unit. Figure 4As shown in the figure, when the chilled water supply temperature, chilled water flow rate and cooling water supply temperature are constant, the cooling water flow rate curve is similar to the chilled water flow rate curve. However, the cooling water flow rate has a greater impact on the COP. The COP at a constant flow rate drops by 3% to 8%. Properly increasing the cooling water flow rate can reduce the condensing temperature and improve the COP, but the water pump energy consumption must be comprehensively considered; too low a flow rate will lead to an increase in the condensing temperature and reduce the efficiency of the refrigeration unit.

[0114] Therefore, in the operation of the refrigeration unit, simply increasing or decreasing a certain control variable cannot guarantee an increase in COP. Although the chilled water supply temperature can be increased as much as possible under the premise of meeting the building load demand, the cooling water flow rate and the chilled water flow rate are affected by the supply and return water temperature difference and the external environment, making it difficult to ensure that the refrigeration unit load rate is maintained at Figure 3 and Figure 4 Therefore, based on the above analysis results, the present invention comprehensively considers the impact of various variables on the COP of the refrigeration unit and adds an energy-saving target control item to the MFAPC optimization objective function to achieve system EER improvement.

[0115] 2. Control scheme design

[0116] 1. Control system structure

[0117] The proposed nonlinear multivariable model-free adaptive predictive control structure for the refrigeration station is as follows: Figure 5 As shown in the figure, the upper optimization layer uses the ICF-MFAPC algorithm to calculate the optimal set values ​​for each loop control variable, including the chilled water supply temperature, chilled water flow rate, and cooling water flow rate. The field control layer uses the PID method to adjust the cooling water pump speed and the chilled water pump speed so that the cooling water flow rate and the chilled water flow rate keep up with the set values ​​calculated by the optimization layer. Due to the complex structure and high cost of the refrigeration unit, the manufacturer does not allow direct control of it. After the upper optimization layer transmits the chilled water supply temperature set value to the field control layer, the refrigeration unit's internal controller can realize the regulation of the compressor and expansion valve to ensure that the chilled water supply temperature keeps up with the optimal set value.

[0118] Figure 5 The upper and middle optimization layers consist of a model-free adaptive predictive controller, pseudo-Jacobian matrix (PJM) prediction, PJM estimation, optimization objectives, time-varying proportional control terms, and time-varying integral control terms. In the figure, y(k+1) is the controlled variable, i.e., the cooling capacity generated by the refrigeration station; u(k) is the optimal value of the control variables, including cooling water flow, chilled water flow, and chilled water supply temperature. The cooling capacity required by the refrigeration station, that is, the set value of the system cooling capacity; The desired values ​​of the control variables cooling water flow, chilled water flow, and chilled water supply temperature are controlled to ensure that each control variable reaches the ideal range while meeting the cooling demand, thereby making the system EER as high as possible; is the PJM estimate, is the PJM forecast value; E y (k) and E u (k) is the transformation matrix.

[0119] 2. Optimization objective function design

[0120] The control strategy proposed in this paper aims to meet building load requirements while simultaneously improving the overall energy efficiency (EER) of the cooling station system. Complex energy transfer and coupling relationships exist between various devices in a cooling station. For example, reducing cooling water flow can save cooling water pump energy, but this may also reduce the COP of the cooling unit. Since cooling units account for over 50% of energy consumption, achieving the control objective requires maximizing the COP of the cooling unit. To achieve this, the system optimization objective function is designed:

[0121]

[0122] Where λ>0 is the weight factor, ζ>0 is the energy saving factor; y * (k+i)∈R m×1 is the expected value of the system controlled variable (cooling capacity) at time k+i, y(k+i)∈R m×1 is the actual value of the controlled variable at time k+i, i=1,2,…N, N is the system output prediction time domain; u * (k+r) is the expected value of each control variable of the system at time k+r, u(k+r) is the actual value of each control variable of the system at time k+r, r=1,2,…N u , N u is the system control time domain; Δu(k+j)=u(k+j)-u(k+j-1).

[0123] The first term in formula (3) is to ensure the tracking accuracy of the system controlled variable; the second term can limit the rate of change of the controlled variable, thereby ensuring the stability of the controlled variable input of the refrigeration station; the third term is to ensure that the controlled variable is as close to the ideal target value as possible, thereby improving the system EER.

[0124] 3. Algorithm construction

[0125] Consider the following discrete-time MIMO system:

[0126] y(k+1)=f(y(j),…,y(kn y ),u(k),…,u(kn u )) (4)

[0127] Where u(k),y(k)∈R m , represents the control quantity and controlled variable of the system at time k; n u and n y is an unknown positive integer; is an unknown nonlinear function.

[0128] The following nonlinear dynamics assumptions are given:

[0129] Assumption 1: f i The partial derivatives of (…)(i=1,…,m) with respect to the control variable u(k) are continuous.

[0130] Assumption 2: System (4) satisfies the generalized Lipschitz condition:

[0131] ‖y(k1+1)-y(k2+1)‖≤q‖u(k1)-u(k2)‖ (5)

[0132] Among them, for any time k1≠k2, k1, k2≥0, and u(k1)≠u(k2); q>0 is a constant.

[0133] Note 1: Assumption 1 is the basic condition for controller design; Assumption 2 indicates that there are constraints on the relationship between the controlled variables and the control variables of the refrigeration station system, and the rate of change is limited, so the characteristics of the refrigeration station system meet the above assumptions.

[0134] Lemma 1: For the refrigeration station system (4) that satisfies Assumptions 1 and 2, there exists a pseudo-Jacobian matrix Φ(k) that can be transformed into the following compact form dynamic linearized data model:

[0135] y(k+1)=y(k)+Φ(k)Δu(k) (6)

[0136]

[0137] And for any time k, Φ(k) is bounded, see the literature for proof.

[0138] In order to facilitate the algorithm stability analysis later, the following assumptions are made:

[0139] Assumption 3: PJMΦ(k) satisfies the diagonal advantage condition, that is, |φ ij (k)| <s1,s2≤|φ ii (k)|≤αs2, i=1,2,…,m,j=1,2,…,m,i≠j, and the signs of the elements in Φ(k) remain unchanged. Where s1 and s2 are two positive numbers, and s2=βs1, where β is a positive number.

[0140] Note 2: Because the relationship between the control variables and the controlled variables of the refrigeration station system is bounded, its characteristics meet Assumption 3.

[0141] The N-step forward prediction value of the system controlled variable can be obtained from formula (6):

[0142]

[0143] Considering j>N u When Δu(k+j-1)=0, the control variable increment coefficient of the controlled variable prediction equation is obtained as:

[0144]

[0145] The N-step prediction equation of the system controlled variable can be written as:

[0146]

[0147] Among them, Y N (k+1)=[y T (k+1),y T (k+2),…,y T (k+N)] T ∈R Nm×1 ;E y (k)=[I m×m ,…,I m×m ] T ∈R Nm×m ;

[0148] because:

[0149] u(k)=u(k-1)+Δu(k) (11)

[0150] according to And formula (11) can be used to get the future N u The optimal control variable for the step:

[0151]

[0152] Then the incremental coefficient of the control variable prediction equation is:

[0153]

[0154] The N of the system control variables can be obtained u Step prediction equation:

[0155]

[0156] in

[0157] make According to formulas (10) and (14), the optimization objective function (3) can be rewritten as:

[0158]

[0159] Substitute (10) and (14) into (15) and apply the optimization condition Under the conditions of satisfying Assumptions 1, 2 and Lemma 1, the following control law can be obtained:

[0160]

[0161] Observe (16), one of the items that plays a role in tracking the expected value of the controlled variable is Define it as:

[0162]

[0163] We can get:

[0164]

[0165] As can be seen from (18), the control law of the MFAPC method with the energy-saving control term includes two items, one is the time-varying error integral control term, and the other is the error control term between the desired control variable and the actual control variable. The former plays the role of making the controlled variable track the set value. The research results show that in the absence of external interference, the traditional MFAC scheme can converge the system controlled variable error to zero, but for complex nonlinear systems such as refrigeration stations, the excessive linearization of MFAC will cause the system information accuracy to decrease, thereby generating overshoot, and the noise and measurement interference in the system can easily cause the controlled variable error to converge to a non-zero constant. Therefore, the disturbance in the system operation and the existence of measurement interference may significantly deteriorate the control performance of MFAC. To solve this problem, the present invention combines incremental PID with MFAPC, introduces additional tracking error information to compensate for the reduction in system information accuracy caused by excessive linearization, and designs the control law as follows:

[0166]

[0167] Among them, kp≥0, ki>0 are weighting coefficients.

[0168] In formula (19), the present invention uses Replace (18) Since (18) is derived from the objective function (15), simply replacing it will cause the algorithm to lose its original calculation logic. By observing (18) and the objective function (16), it can be seen that The coefficient of is related to the coefficient of the term where Err(k) exists. When kp=0, The previous coefficient θ becomes ki. Similarly, when kp≠0, θ>0, so that the control variable increment calculated by the result of (19) stabilizes the system. Here, θ is an adjustable parameter whose purpose is to make the control algorithm more general. When the dimension of the system I / O data is large, the matrix inversion calculation is very time-consuming. Therefore, Equation (19) can be simplified using a heuristic simplification method to:

[0169]

[0170] Then the control quantity can be obtained:

[0171]

[0172] in

[0173] It can be seen that when kp=0, ki=1, θ=1, ζ=0, the control law (21) degenerates into the MFAPC method control law, that is, the original MFAPC method is a special case of the improved method.

[0174] 4. Gray Wolf Optimization Algorithm Parameter Tuning

[0175] The Grey Wolf Optimizer (GWO) is a swarm intelligence optimization algorithm that simulates the social behavior and hunting strategies of gray wolves. Proposed by Seyedali Mirjalili et al. in 2014, the GWO algorithm seeks the optimal solution to an optimization problem by simulating the social structure and hunting behavior of gray wolves in nature. Its core concept is derived from the group cooperation and leadership division of gray wolves. In gray wolf hunting, the pack employs a strategy of encircling prey, gradually narrowing the encirclement until the prey is captured and shared. This collaborative hunting behavior provides inspiration for the algorithm's search for the optimal solution.

[0176] The standard GWO algorithm uses a problem-specific fitness function, which typically measures the quality of the solution based on the objective function value. To improve the algorithm's performance in specific application scenarios, this paper improves on the traditional fitness function. Considering the energy-saving optimization needs of refrigeration stations and the cooling capacity tracking effect, the ICF-MFAPC parameter tuning Grey Wolf Optimization Algorithm fitness function is designed:

[0177]

[0178] Where fitness>0; Represents the expected EER during the optimization process, which is generally a constant. y1(k+1) represents the EER at time k+1. represents the cooling demand at time k+1, and y2(k+1) represents the cooling capacity at time k+1; eer >0 is the penalty coefficient related to EER.

[0179] In order to obtain the minimum fitness, improve EER and meet the cooling load requirements, and k eer A reasonable value needs to be set. It needs to be set according to the actual operating conditions of the refrigeration station system and the maximum EER that can be achieved during operation; eer If it is too large, the algorithm will only pursue energy efficiency improvement while ignoring the cooling load demand. If it is too small, the energy efficiency improvement effect will be poor. At the same time, the order of the first item in formula (22) is irreplaceable with the actual EER. Its meaning is that when When , EER does not reach the expected value, and fitness increases; When , EER reaches the expected value and fitness decreases.

[0180] 5. Refrigeration station control strategy

[0181] According to the above control law algorithms (20) and (21), as well as the existing PJM estimation algorithm and PJM prediction algorithm, the present invention designs a model-free predictive control strategy for the refrigeration station as follows:

[0182] 1) Pseudo-Jacobian matrix estimation

[0183]

[0184] in, is the estimated value of PJM; η∈(0,2] is the step size factor of PJM update, which can enhance the flexibility of the algorithm, μ>0; and:

[0185]

[0186] To improve the tracking ability of the algorithm and satisfy Assumption 3, PJM estimates that the following provisions need to be made:

[0187] like or or

[0188] like or

[0189] 2) Pseudo-Jacobian matrix prediction

[0190]

[0191] Where θ i (k)∈R m×m ,i=1,…,N p is the order of the multi-layer hierarchical forecasting method; j = 1,…,N u -1;

[0192] θ(k) update:

[0193]

[0194] Where δ∈(0,1]; and must satisfy: θ(k)=θ(0), if ‖θ(k)‖≥M θ , M θ Is a positive number.

[0195] 3) Model-free adaptive control law update

[0196]

[0197] According to (23), (26) and (27), the parameters that need to be optimized for GWO are obtained. is also unknown, so (22) can be written as follows:

[0198]

[0199] The algorithm flow is as follows Figure 6 shown.

[0200] 6. Algorithm stability analysis

[0201] 6.1 Prerequisites

[0202] In order to analyze the stability and convergence of the ICF-MFAPC algorithm, the following assumptions and theorems are given:

[0203] Assumption 4 For a given bounded expected signal y * (k+1), there is always a bounded control variable signal u * (k), so that the controlled variable signal of the system y(k+1)=y * (k+1).

[0204] Theorem 1 For discrete-time nonlinear system (4) that satisfies assumptions 1 to 4, when y * (k+1)=y r =const,u * (k)=u r =const, for assumption 3, s2 = βs1, when β satisfies certain conditions, using the ICF-MFAPC method, there is always λ min>0, so that λ>λ min When:

[0205] (1) The tracking error of the controlled variable of the system is convergent, that is, Where γ is a positive constant.

[0206] (2) The sequences of controlled variables {y(k)} and control variables {u(k)} are bounded.

[0207] The proof consists of two parts. The first part proves that the tracking error of the controlled variable of the system is bounded. The second part proves that the sequence of controlled variables {y(k)} and control variables {u(k)} is bounded. For the convenience of proof, the norms used in the following text are all matrix infinite norms.

[0208] 6.2 The tracking error of the system controlled variable is bounded

[0209] The system tracking error is specified as e(k+1)=y r -y(k+1), when y * (k+1)=y r =const,u * (k)=u r =const, according to (27) and (28):

[0210]

[0211] because:

[0212]

[0213] have:

[0214]

[0215] Substitute equation (6) into e(k+1)=y r -y(k+1) gives:

[0216] e(k+1)=y r -y(k)-Φ(k)Δu(k) (35)For the convenience of proof, define:

[0217]

[0218] From (34), (35) and (36), we can get:

[0219] e(k+1)=a1(k)e(k)+a2(k)e(k-1)+a3(k)[u(k-1)-u r ] (37)

[0220] From (28), (34) and (36), we can get:

[0221] u(k)-u r =b1(k)[u(k-1)-u r ]+b2(k)e(k)-b3(k)e(k-1) (38)

[0222] Take the infinite norm on both sides of formula (37):

[0223]

[0224] Where χ is any positive number. To prove that e(k+1) converges, we need to know the infinite norm of the matrix in equation (36), and convergence.

[0225] (1) Infinity norm of the correlation matrix

[0226] Assume z 1,i 、z 2,i and z 3,i (i=1,2,…,m) is the sum of the absolute values ​​of the elements in the i-th row of matrices a1(k), a2(k) and a3(k); and assuming that z 4,i 、z 5,i and z 6,i are the sum of the absolute values ​​of the elements in row i of matrices b1(k), b2(k) and b3(k). We can get:

[0227]

[0228] Where (40) and (43) are min When exists:

[0229]

[0230]

[0231] According to Assumption 3: |φ ij (k)| <s1,s2≤|φ ii (k)|≤αs2, i=1,2,…,m,j=1,2,…,m,i≠j,s2=βs1, we can get:

[0232] ①z 1,i

[0233] The first term is the absolute value of the diagonal elements of a1(k), according to (46):

[0234]

[0235] The second term can be expressed as:

[0236]

[0237] but:

[0238]

[0239] ②z 2,i

[0240] The first term is also the absolute value of the diagonal elements of a2(k):

[0241]

[0242] Because z 2,i The second term of z 1,i Similarly, we can get:

[0243]

[0244] ③z 3,i

[0245]

[0246] ④z 4,i

[0247] From (47) we can get:

[0248]

[0249] ⑤z 5,i

[0250]

[0251] ⑥z 6,i

[0252]

[0253] Consider the relationship between the infinity norm of a matrix and the sum of the absolute values ​​of the elements in the i-th row of the matrix:

[0254]

[0255] (2) Bounded

[0256] Taking the infinite norm on both sides of formula (38) yields:

[0257]

[0258] Since ‖b2(k)‖>‖b3(k)‖, according to (57), formula (58) can be rewritten as:

[0259]

[0260] make From (57) and (59), we can obtain:

[0261]

[0262] Assume:

[0263]

[0264] ①

[0265] From (53), there exists λ min >0, when λ > λ min then

[0266] ②d4 < 1

[0267] From (54), there exists λ min >0, when λ > λ min then d4 < 1.

[0268] ③

[0269] That is, 2d3d5 < χ(1 - d4). According to (53) - (55), we can obtain:

[0270]

[0271] Since χ is an arbitrary positive number, for any ζ, there exists λ min >0, when λ > λ min then

[0272] ④‖a1(k)‖ + ‖a2(k)‖ + χ < d1 + d2 + χ < d7 < 1

[0273] According to equations (50) and (52), we can obtain:

[0274]

[0275] When That is then (63) is a parabola opening upwards. At this time, there always exists β such that 0 < d1 + d2 < 1 holds. Then there exists λ min >0, when λ > λ min then there always exists a positive number χ such that d1 + d2 + χ = d7 < 1.

[0276] Then (39) is satisfied:

[0277] ‖e(k+1)‖≤d7max{‖e(k)‖,‖e(k-1)‖,h(k-1)} <max{‖e(k)‖,‖e(k-1)‖,h(k-1)} (64)

[0278] ‖e(k)‖ <max{‖e(k-1)‖,‖e(k-2)‖,h(k-2)} (65)

[0279] Substituting (65) into (64) yields:

[0280] ‖e(k)‖ <max{h(k-1),h(k-2),…,h(0),‖e(1)‖,‖e(0)‖} (66)

[0281] Substituting (66) into (60), we can obtain the following equation according to d4<1 in equation (61):

[0282]

[0283] Assume the following function:

[0284]

[0285] There are two situations:

[0286] ①f(k-1)≤max{f(k-2),f(k-3),…,f(0),‖e(1)‖,‖e(0)‖}:

[0287]

[0288] ②f(k-1)>max{f(k-2),f(k-3),…,f(0),‖e(1)‖,‖e(0)‖}:

[0289]

[0290] have:

[0291]

[0292] From (61) and (72), we can see that f(k) is bounded, that is, h(k) is bounded, and

[0293] (3) e(k+1) converges

[0294] When any ‖e(k)‖>h(k) holds, equation (64) exists:

[0295] ‖e(k)‖ <d7max{‖e(k-1)‖,‖e(k-2)‖} (73)

[0296] ‖e(k-1)‖ <d7max{‖e(k-2)‖,‖e(k-3)‖} (74)

[0297] Merge (73) and (74):

[0298] max{‖e(k)‖,‖e(k-1)‖}≤d7max{‖e(k-2)‖,‖e(k-3)‖} (75)

[0299] We can get:

[0300]

[0301] It can be seen that ‖e(k+1)‖ converges. but

[0302] 6.3 The sequences of controlled variables {y(k)} and control variables {u(k)} are bounded

[0303] Since e(k)=y r -y(k) and e(k) and y r is bounded, so the controlled variable sequence {y(k)} is bounded.

[0304] Since e(k) and u r Bounded. From (36) and (37), we can see that the control variable sequence {u(k)} is bounded.

[0305] 3. Experimental Research

[0306] 1. Construction of neural network model of refrigeration station system

[0307] This study investigated the effectiveness of the proposed method using a refrigeration system in a public building in Jiangsu Province. The building is a commercial complex with an area of ​​approximately 13,000 square meters. Detailed information on the refrigeration system components is provided in Table 1. This refrigeration system is a typical public building configuration, and the experimental results offer general guidance.

[0308] Table 1 Refrigeration station system equipment information

[0309] Tab.1 Equipment Information of Refrigeration Station System

[0310]

[0311] At present, the refrigeration station system adopts local loop PID to control the chilled water supply temperature, chilled water flow and cooling water flow to meet the cooling load demand of the building, but does not consider the improvement of the overall energy efficiency of the system. Since it is not suitable to conduct online debugging of unverified control algorithms during the operation of the refrigeration station system, the present invention first establishes a refrigeration station system prediction model to simulate the control system to verify the effectiveness of the proposed method. Due to the complex structure of the refrigeration station system, changes in outdoor temperature, humidity and indoor heat source will affect the cooling effect and EER. At the same time, there is strong coupling between the components within the system, and the actuator has nonlinearity, which makes mechanism modeling more difficult. Kolmogorov et al. have proved that a three-layer feedforward neural network model can represent any nonlinear relationship. Therefore, the present invention adopts a three-layer feedforward neural network to construct a dynamic model of the refrigeration station system.

[0312] In order to ensure the accuracy of the model, the model input needs to include the control variable, namely the chilled water flow rate f pch , Chilled water supply temperature T sucw and cooling water flow rate f pc . In addition, the outdoor dry-bulb temperature T and relative humidity RH as environmental conditions will affect the cooling effect of chilled water and the workload of the system, and have a great impact on the operating efficiency of the refrigeration station; EER and cooling capacity Q reflect the operating status of the refrigeration station system. Therefore, the input of the neural network model of the refrigeration station includes chilled water temperature and flow, cooling water flow, dry-bulb temperature and relative humidity, as well as EER and cooling capacity at the current moment. Considering the system control objectives, the model output includes two variables, namely EER and cooling capacity. The model structure is as follows Figure 7 shown.

[0313] The range of the number of neurons in the hidden layer of the neural network is determined by the following formula:

[0314]

[0315] where m hid is the number of neurons in the hidden layer; m in is the number of neurons in the input layer; m out is the number of neurons in the output layer; m ran is a positive integer, and 1 <m ran <10. The number of neurons in the hidden layer is calculated to be 12.

[0316] The present invention collects the refrigeration station system operation data and environmental data through sensors with a collection interval of 5 minutes. Some of the data are shown in Table 2.

[0317] Table 2 Partial operating data of refrigeration station system

[0318] Tab.2 Partial operating data of refrigeration station system

[0319]

[0320] The neural network training algorithm uses the Bayesian regularization method, which combines model uncertainty, adaptive regularization and prior knowledge through a probabilistic framework to avoid overfitting and improve the generalization ability of the model. The final neural network modeling effect is shown in Figure 8 and Figure 9 , where the horizontal axis is the sampling time point. The relative errors of EER and cooling capacity prediction are both within ±5%, which can meet the needs of engineering and control system simulation experiments.

[0321] 2. Control system simulation experiment

[0322] Based on the fitness function (30), the present invention uses GWO to optimize the parameters of the ICF-MFAPC algorithm (the parameter setting results are shown in Table 3), and simulates the data of the refrigeration station running continuously for one month. * (k+1) and GWO algorithm The cooling demand value at time k+1 of the original data is taken from the GWO algorithm. Take the EER value of the original PID control at time k+1, and the simulation results of EER and cooling capacity are as follows: Figure 10 and Figure 11 As shown. After calculation, we know that Figure 10 The average EER of ICF-MFAPC is 3.9998, and the average EER of the original refrigeration station local loop PID control is 3.4387. The EER of the proposed control strategy is improved by 16.3%, and similar results are obtained in multiple experiments. Figure 11 It can be seen that the ICF-MFAPC algorithm can ensure the stable operation of the system. It can well meet the building cooling load demand during the system startup phase when the cooling capacity is low and the system operation phase when the cooling capacity is stable. The following error in the stable phase does not exceed 5%.

[0323] In order to compare the control effect of the proposed method with that of the MFAPC method, the present invention uses the data of one day of operation of the refrigeration station and adopts GWO combined with the fitness function (30) to optimize the parameters of the MFAPC and ICF-MFAPC methods respectively. The adjusted parameters are shown in Table 3, where u is used to ensure the consistency of the order of magnitude of the input and output information of the algorithm. * is the normalized value. The simulation results of EER and cooling capacity are shown in Figure 12 and Figure 13As shown in Figure 2, when the GWO fitness is minimum, MFAPC has two situations: good cooling load tracking but energy consumption, and lagging cooling load tracking but energy saving. That is, it is difficult to take into account both cooling load tracking and energy saving. Therefore, MFAPC-1 and MFAPC-2 are used here to represent the above two situations respectively.

[0324] Table 3 Controller parameters

[0325] Tab.3 Controller parameters

[0326]

[0327] From the experimental results, it can be concluded that the cooling capacity following of MFAPC-2 has oscillations, while MFAPC-1 has obvious cooling capacity following lag. Due to the introduction of time-varying proportional-integral control terms in the optimization objective function, the ICF-MFAPC method can make the cooling capacity provided by the refrigeration station better follow the cooling capacity required by the building. The maximum relative error is 4.04%, which is lower than 6.19% and 8.05% of MFAPC-1 and MFAPC-2, and the tracking accuracy is improved by 30.7% compared with MFAPC-2, which has a good tracking effect on demand cooling capacity. After the system stabilizes, the EER of the ICF-MFAPC method can reach up to about 4.206, which is higher than the 3.766 of the original PID control and 4.054 of MFAPC-1, and the EER is improved by 12.98% and 6.03% on average. At the same time, the traditional MFAPC method has obvious cooling capacity following lag in the case of MFAPC-1, and the cooling capacity following has oscillations in the case of MFAPC-2, indicating that the basic MFAPC method is difficult to take into account both cooling capacity demand and EER improvement. Figures 14-16 It can be seen that the control variables under different strategies can meet the building cooling demand, but the EER results are different. Under the ICF-MFAPC framework, the GWO algorithm obtains the optimal setting value u of the three control variables of chilled water supply temperature, chilled water flow and cooling water flow by optimization. * , then improve the MFAPC algorithm to keep the control variable at u * , thereby achieving the effect of improving EER. The energy efficiency of refrigeration plant equipment affects each other, requiring comprehensive consideration of control variables to optimize coordinated operation, balance load distribution, and achieve optimal overall energy efficiency. However, traditional PID control and MFAPC methods lack the ability to adjust control variables, making it difficult to achieve the goal of meeting cooling demand while improving EER.

[0328] Because refrigeration plants are complex, energy-intensive, and nonlinear systems, model-based control methods face difficulties in modeling the plant's system mechanisms. Data-driven modeling is heavily dependent on data quality, and basic MFAC methods struggle to balance dynamic performance and energy efficiency optimization. To address this issue, this paper proposes an ICF-MFAPC method. By analyzing the impact of chilled and cooling water flow and temperature on the energy efficiency of the refrigeration unit, an energy-saving target control term is constructed. The system output error is introduced into the control rate update function, and a time-varying proportional-integral control term is designed to effectively compensate for information loss during the linearization process, enhancing the system's dynamic performance and robustness under varying load conditions. Furthermore, the rolling optimization mechanism of model predictive control is combined to overcome the effects of large lags and inertia in the refrigeration plant system. To address the ICF-MFAPC parameter tuning issue, the GWO algorithm is introduced, and an improved GWO fitness function is designed based on the control objective. Finally, this paper provides a detailed mathematical proof of the closed-loop stability and convergence of the ICF-MFAPC algorithm. A neural network prediction model is constructed for a commercial building refrigeration plant system in Jiangsu Province, and simulation experiments are conducted to demonstrate the proposed method. The results show that compared with the traditional PID and MFAPC methods, the proposed ICF-MFAPC method improves the dynamic and steady-state performance of the system. Compared with the original PID control of the cooling station, the system EER is improved by 16.3% while meeting the building load demand.

Claims

1. A nonlinear multivariable model-free adaptive predictive control method for a refrigeration station; characterized in that: The following steps are involved: S1. Analyze the energy efficiency of each circuit in the refrigeration station and determine the control variables to be optimized as chilled water supply temperature, cooling water flow rate, and chilled water flow rate; S2. Build a control system: The control system adopts a structure that combines an upper optimization layer and a field control layer. The upper optimization layer includes a model-free adaptive predictive controller, PJM prediction, PJM estimation, optimization objectives, time-varying proportional control terms, and time-varying integral control terms. The upper optimization layer can calculate the optimal set values ​​of the control variables of each loop. The field control layer can use the PID method to adjust the operation of the field actuator equipment so that the chilled water supply temperature, cooling water flow rate, and chilled water flow rate follow the optimal set values ​​of the control variables calculated by the upper optimization layer. S3. Execution control.

2. The nonlinear multivariable model-free adaptive predictive control method for a refrigeration station according to claim 1, characterized in that: In step S2, the upper optimization layer calculates the optimal setting value of each loop control variable, including: S2.1: Design the optimization objective function, the optimization objective function is as shown in formula (3); Where λ>0 is the weight factor, ζ>0 is the energy saving factor; y * (k+i)∈R m×1 is the expected value of the controlled variable of the system at time k+i, y(k+i)∈R m×1 is the actual value of the controlled variable at time k+i, i=1,2,…N, N is the system output prediction time domain; u * (k+r) is the expected value of each control variable of the system at time k+r, u(k+r) is the actual value of each control variable of the system at time k+r, r=1,2,…N u , N u is the system control time domain; where Δu(k+j)=u(k+j)-u(k+j-1); S2.2: Use the MFAPC method with energy-saving control items to calculate the optimal set values ​​of the control variables of each loop: Consider the following discrete-time MIMO system: y(k+1)=f(y(k),…,y(kn y ),u(k),…,u(kn u )) (4)where u(k),y(k)∈R m , represents the control quantity and controlled variable of the system at time k; n u and n y is an unknown positive integer; is an unknown nonlinear function; The following nonlinear dynamics assumptions are given: Assumption 1: f i The partial derivative of (…)(i=1,…,m) with respect to the controlled variable u(k) is continuous; this is the basic condition for controller design; Assumption 2: System equation (4) satisfies the generalized Lipschitz condition, indicating that the relationship between the controlled variables and the control variables of the refrigeration station system is constrained and the rate of change is limited; as shown in equation (5): ‖y(k1+1)-y(k2+1)‖≤q‖u(k1)-u(k2)‖ (5)where, for any time k1≠k2, k1, k2≥0, and u(k1)≠u(k2); q>0 is a constant; Lemma 1: For the refrigeration station system (4) that satisfies Assumptions 1 and 2, there exists a pseudo-Jacobian matrix Φ(k) that can be transformed into the following compact form dynamic linearized data model: y(k+1)=y(k)+Φ(k)Δu(k) (6) And for any time k, Φ(k) is bounded; In order to facilitate the algorithm stability analysis later, the following assumptions are made: Assumption 3: PJMΦ(k) satisfies the diagonal advantage condition, that is, |φ ij (k)| <s1,s2≤|φ ii (k)|≤αs2, i=1,2,…,m,j=1,2,…,m,i≠j, and the signs of the elements in Φ(k) remain unchanged; s1 and s2 are two positive numbers, and s2=βs1, where β is a positive number; The N-step forward prediction value of the system controlled variable can be obtained from formula (6): Considering j>N u When Δu(k+j-1)=0, the control variable increment coefficient of the controlled variable prediction equation is obtained as: The N-step prediction equation of the system controlled variable can be written as: Among them, Y N (k+1)=[y T (k+1),y T (k+2),…,y T (k+N)] T ∈R Nm×1 ;E y (k)=[I m×m ,…,I m×m ] T ∈R Nm ×m ; because u(k)=u(k-1)+Δu(k) (11) according to And formula (11) can be used to get the future N u The optimal control variable for the step: Then the incremental coefficient of the control variable prediction equation is: The N of the system control variables can be obtained u Step prediction equation: in make According to formulas (10) and (14), the optimization objective function (3) can be rewritten as: Substitute formula (10) and formula (14) into formula (15) and use the optimization condition Under the conditions of satisfying Assumptions 1, 2 and Lemma 1, the following control law can be obtained: Observe formula (16), where the item that plays a role in tracking the expected value of the controlled variable is Define it as: We can get: S2.3: S4, use the Grey Wolf optimization algorithm to adjust the parameters.

3. The nonlinear multivariable model-free adaptive predictive control method for a refrigeration station according to claim 2, characterized in that: In step S2.2, the incremental PID method is combined with the MFAPC method to introduce additional tracking error information to compensate for the reduction in system information accuracy caused by over-linearization, thereby obtaining the ICF-MFAPC method. The control law is designed as follows: Among them, kp≥0, ki>0 are weighting coefficients; Formula (19) is simplified to: Then the control quantity can be obtained: in 4. The nonlinear multivariable model-free adaptive predictive control method for a refrigeration station according to claim 3, characterized in that: In step 2.3, the fitness function of the gray wolf optimization algorithm is as follows: Where fitness>0; Represents the expected EER during the optimization process, which is generally a constant. y1(k+1) represents the EER at time k+1.

5. The nonlinear multivariable model-free adaptive predictive control method for a refrigeration station according to claim 4, characterized in that: The control strategy in step 3 execution control includes: S3.1 Pseudo-Jacobian matrix estimation; S3.2 Pseudo-Jacobian matrix prediction; S3.3 Model-free adaptive control law update.

6. The nonlinear multivariable model-free adaptive predictive control method for a refrigeration station according to claim 5, characterized in that: The Jacobian matrix estimation in S3.1 includes the following algorithms: in, is the estimated value of PJM; η∈(0,2] is the step size factor of PJM update, which can enhance the flexibility of the algorithm, μ>0; and: To improve the tracking ability of the algorithm and satisfy Assumption 3, PJM estimates that the following provisions need to be made: like or or i=1,2,…,m like or i,j=1,2,…,m,i≠j Preferably, the Jacobian matrix prediction in S3.2 includes the following algorithm: Where θ i (k)∈R m×m ,i=1,…,N p is the order of the multi-layer hierarchical forecasting method; j = 1,…,N u -1; θ(k) update: Where δ∈(0,1]; and must satisfy: θ(k)=θ(0), if ‖θ(k)‖≥M θ , M θ Is a positive number.

7. The nonlinear multivariable model-free adaptive predictive control method for a refrigeration station according to claim 6, characterized in that: The model-free adaptive control law update in S3.3 includes the following algorithms: According to (23), (26) and (27), the parameters that need to be optimized for GWO are obtained. is also unknown, so (22) can be written as follows:

Images